Mandelbrot con GMP

Avevo gia' utilizzato la libreria GMP in questo post per il calcolo di PiGreco
Questa volta la utilizzero' per effettuare zoom estremi dell'insieme di Mandelbrot grazie alla sua caratteristica di poter utilizzare una precisione arbitraria ovvero di poter settare il numero di bit che rappresentano una quantita' in virgola mobile. Il tipo float cambia a seconda da compilatore a compilatore e da processore a processore ma diciamo che in generale puo' essere considerato a 32 bit, il che va bene quando il numero di cifre significative non e' molto esteso

Per mostrare la differenza nell'uso di GMP rispetto al solo tipo float ho scritto due programmi ed ho effettuato due zoom dell'insieme di Mandelbrot

Zoom 1:

In questo caso il risultato e' identico

GMP

Float

Zoom 2
Riducendo la finestra di visualizzazione di centesimo rispetto allo zoom precedente GMP continua a lavorare in modo corretto mentre il programma scritto in C++ puro "esplode"

GMP

Float

Naturalmente tutto cio' ha un prezzo in termini di memoria consumata (la precisione e' vero che e' arbitraria ma deve comunque fare i conti con la ram disponibile) ed il tempo di calcolo (ovviamente piu' lento con GMP).Curiosamente provando a lanciare questo programma in macchina virtuale VMWare, il programma mi diceva che non aveva abbastanza memoria disponibile e sono passato ad un calcolatore reale

esempio con l'utilizzo di GMP
----------------------------------------------------------

#include <stdio.h>
#include <math.h>
#include <gmp.h>
#include <iostream>

using namespace std;

#include <pngwriter.h>

#define SCREEN_WIDTH 1280
#define SCREEN_HEIGHT 960


int main(void){

mpf_t A,B, Im, DRe, DIm,test;
mpf_t X,Y,XN,YN,XX,YY,Xs,Ys;
float test2;

int j,k,i;

long int precision = 256;
double r;




mpf_set_default_prec(precision);

mpf_init2(A,precision);
mpf_init2(B,precision);
mpf_init2(Im,precision);
mpf_init2(DRe,precision);
mpf_init2(DIm,precision);
mpf_init2(test,precision);
mpf_init2(X,precision);
mpf_init2(Y,precision);
mpf_init2(XN,precision);
mpf_init2(YN,precision);
mpf_init2(XX,precision);
mpf_init2(YY,precision);
mpf_init2(Xs,precision);
mpf_init2(Ys,precision);

// punto vertice
// Re = -0.74364085
// Im = 0.13182733
// Diametro = 0.00012068
// Zoom = 25497

pngwriter png(SCREEN_WIDTH,SCREEN_HEIGHT,0,"mandelbrot2.png");
double iterazioni = 1024;

mpf_set_d(A,-0.74364085);
//mpf_set_d(A,-0.74370119);
mpf_set_d(B,0.13182733);
mpf_set_d(Im,0.13182733);
mpf_set_d(DRe,0.000000009428125);
mpf_set_d(DIm,0.00000001257083);



mpf_set_d(X,0.0);
mpf_set_d(Y,0.0);
mpf_set_d(Xs,0.0);
mpf_set_d(Ys,0.0);
mpf_set_d(XX,0.0);
mpf_set_d(YY,0.0);
mpf_set_d(XN,0.0);
mpf_set_d(YN,0.0);

cout << "Inizio" << endl;

for (i=0;i<SCREEN_HEIGHT;i++)
    {
    mpf_add(A,A,DRe);
    mpf_set(B,Im);
    for (j=0;j<SCREEN_WIDTH;j++)
    {
    mpf_add(B,B,DIm);
    mpf_set_d(X,0.0);
    mpf_set_d(Y,0.0);

    for (k=0;k<=iterazioni;k++)
        {
        mpf_pow_ui(Xs,X,2);
        mpf_pow_ui(Ys,Y,2);
        mpf_sub(XN,Xs,Ys);
        mpf_add(XN,XN,A); //XN = (x*x)-(y*y) +a
        mpf_mul(YN,X,Y);
        mpf_mul_ui(YN,YN,2);
        mpf_add(YN,YN,B); //YN = 2*X*Y

        mpf_pow_ui(XX,XN,2);
        mpf_pow_ui(YY,YN,2);
        mpf_add(test,XX,YY);// test = (XN*XN)+(YN*YN)
        test2 = mpf_get_d(test);
        if (test2 > 4.0)
            {
            r = k%2;
            png.plot(j,i,r,r,r);
            break;
            }
        mpf_set(X,XN);
        mpf_set(Y,YN);
        }
    }
    }

png.close();

mpf_clear(A);
mpf_clear(B);
mpf_clear(A);
mpf_clear(Im);
mpf_clear(DRe);
mpf_clear(DIm);
mpf_clear(X);
mpf_clear(Y);
mpf_clear(Xs);
mpf_clear(Ys);
mpf_clear(XX);
mpf_clear(YY);
mpf_clear(XN);
mpf_clear(YN);
mpf_clear(test);

cout << "Fine" << endl;

return 0;
}


----------------------------------------------------------
esempio in C++ puro
----------------------------------------------------------
#include <pngwriter.h>

#define SCREEN_WIDTH 1280
#define SCREEN_HEIGHT 960



// punto vertice
// Re = -0.74364085
// Im = 0.13182733
// Diametro = 0.00012068
// Zoom = 25497

double iterazioni = 1024;

float a = -0.74364085;
float b = 0.13182733;
float im = 0.13182733;

float dre = 0.000000009428125;
float dim = 0.00000001257083;


double r;

float x,y,x_new,y_new;

int test;

int k,j,i;

int keypress = 0;


int main() {

pngwriter png(SCREEN_WIDTH, SCREEN_HEIGHT,0,"mand_norm2.png");


for (i=0;i<SCREEN_HEIGHT;i++)
 {
 a=a+dre;
 b=im;
 for (j=0;j<SCREEN_WIDTH;j++)
  {
  b = b+dim;

  x = 0;
  y = 0;
  test = 0;

  for (k=0;k<=iterazioni;k++)
   {
   x_new = (x*x)-(y*y)+a;
   y_new = (2*x*y)+b;
   if (((x_new*x_new)+(y_new*y_new))>4)
    {
    r = k%2;
    png.plot(j,i, r, r, r);
    break;
    }
   x = x_new;
   y = y_new;
   }

  }

 }
png.close();
return(0);

}

----------------------------------------------------------