# mandelbrot fractal

folkert folkert at vanheusden.com
Mon Sep 7 11:58:03 UTC 2015

```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);

}

}

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