mandelbrot fractal

[folkert]

Hi,

I'm trying to calculate the mandelbrot fractal using GMP.
I've got it working and it is amazingly fast, even for data-types of
256 bits in size!
Now what I'm curious about: can this piece of code be improved in some
way? From a GMP-lib usage point of view I mean.

// written by folkert at vanheusden.com
// released under agpl v3.0
#include <stdio.h>
#include <gmpxx.h>

void my_init_mpf(mpf_t *const what, const int n_bits)
{
	mpf_init(*what);
	mpf_set_prec(*what, n_bits);
}

int main(int argc, char *argv[])
{
	const int n_bits = 256;

	mpf_t x1, x2, y1, y2;
	mpf_init_set_d(x1, -2.0);
	mpf_init_set_d(y1, -2.0);
	mpf_init_set_d(x2, 2.0);
	mpf_init_set_d(y2, 2.0);

	const int max_it = 256;
	const int nx = 640;
	const int ny = 480;

	printf("P6 %d %d %d\n", nx, ny, max_it - 1);

	mpf_t two;
	mpf_init_set_d(two, 2.0);

	mpf_t nxg, nyg;
	mpf_init_set_ui(nxg, nx); // dummyx = xres
	mpf_init_set_ui(nyg, ny); // dummyy = yres

	mpf_t tx, ty;
	mpf_inits(tx, ty, NULL);
	mpf_sub(tx, x2, x1); // d2x = x2 - x1
	mpf_sub(ty, y2, y1); // d2y = y2 - y1

	mpf_t dx, dy;
	my_init_mpf(&dx, n_bits);
	my_init_mpf(&dy, n_bits);
	mpf_div(dx, tx, nxg); // dx = d2x / dummyx
	mpf_div(dy, ty, nyg); // dy = d2y / dummyy

	mpf_t cx, cy;
	my_init_mpf(&cx, n_bits);
	my_init_mpf(&cy, n_bits);

	mpf_t t, wx, wy, sqx, sqy;
	my_init_mpf(&t, n_bits);
	my_init_mpf(&wx, n_bits);
	my_init_mpf(&wy, n_bits);
	my_init_mpf(&sqx, n_bits);
	my_init_mpf(&sqy, n_bits);

	mpf_set(cy, y1);
	for(int y=0; y<ny; y++)
	{
		mpf_set(cx, x1);

		for(int x=0; x<nx; x++)
		{
			mpf_set_d(wx, 0.0);
			mpf_set_d(wy, 0.0);

			mpf_set_d(sqx, 0.0);
			mpf_set_d(sqy, 0.0);

			int it = 0;

			do {
				mpf_mul(t, wy, two);	// y = 2 * y
				mpf_mul(t, t, wx);	// * x
				mpf_add(wy, t, cy);	// + cy

				//

				mpf_sub(wx, sqx, sqy);	// x = x^2 - y^2
				mpf_add(wx, wx, cx);	// x = z + cx

				//

				mpf_mul(sqx, wx, wx);	// sqx = x^2
				mpf_mul(sqy, wy, wy);	// sqy = y^2

				//

				mpf_add(t, sqx, sqy);	// t = x^2 + y^2
			}
			while(mpf_cmp_d(t, 16.0) < 0 && ++it < max_it);

			double dwx = mpf_get_d(wx);
			double dwy = mpf_get_d(wy);

			if (it < max_it)
				printf("%c%c%c", int(it - dwx * 2.0), max_it - it, int(it + dwy * 2.0));
			else
				printf("%c%c%c", 0, 0, 0);

			mpf_add(cx, cx, dx);
		}

		mpf_add(cy, cy, dy);
	}

	return 0;
}


Please note that it emits a .pnm-file to stdout. So if you want to run
it yourself, invoke it like this:

./a.out > test.pnm

This test.pnm-file can be viewed with any Linux graphics viewer.


regards,

Folkert van Heusden

-- 
You've probably gotten really fed up with never winning in the Mega-
Millions lottery. Well, weep no longer: www.smartwinning.info tells
you everything that might help you deciding what numbers to choose.
With nice graphs and pretty animations!