mandelbrot fractal

[Vincent Diepeveen]

hi Folkert - you probably want to consider building a small piece of gpgpu 
code for it in CUDA. The gpu's have a very small instruction set, yet they 
do have quite some conditional instructions. So should be able to 
calculate for a streamcore the next pixel even if other cores still are 
busy with their current one, without 'if then else' type constructions.

AMD gpu's only have 24 bits available to you and the too limited 
instruction set of OpenCL from conditional viewpoint seen with 32 bits is 
basically 4x slower at AMD gpu's (and the top bits multiplication with 24 
bits would also involve 4 units as there is no OpenCL equivalent for it 
other than the one casted at 32x32 bits), whereas very cheap bit older 
Nvidia gpu's have fast 32 bits integer calculations and you can pickup 
even an older generation Tesla with 512 streamcores for just 192 dollar or 
so on ebay.

If you are just interested in the 32 bits integer math then a gamers chip 
would also be interesting f you underclock it a little.

That's teraflops of calculation power for under 200 euro and just a bit 
of puzzling with conditional instructions to move to the next pixel if 
current one has been calculated.

This is not too many lines of codes - one weekend of coding.

No CPU library can rival that.

Regards,
Vincent

BTW very interested in the program you got now to see how with more bits 
precision mandelbrot looks like - where is it for download?

On Mon, 7 Sep 2015, folkert wrote:

> 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!
> _______________________________________________
> gmp-discuss mailing list
> gmp-discuss at gmplib.org
> https://gmplib.org/mailman/listinfo/gmp-discuss
>