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

}

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

Pi Greco con OpenMP e GMP

Per testare il calcolo parallelo e l'Hyperthreading avevo la necessita' di un algoritmo parallelizzabile e con una convergenza piuttosto lenta per poter apprezzare le differenze di tempo di calcolo nella versione parallela e non parallela.

Scartato il calcolo di Pi Greco mediante la formula di Gauss (gia' utilizzata in un precedente post) perche' ricorsiva e quindi non parallelizzabile e scartato l'insieme di Mandelbrot perche' comunque un minimo complesso , ho provato a determinare il valore di Pi Greco mediante lo sviluppo in serie di Taylor mediante la formula

La sommatoria converge con estrema lentezza al valore di Pi Greco

Standard (tempo di esecuzione 0.049s)
-----------------------------------------------------

#include <cmath>
#include <iostream>

using namespace std;

int main()
{

float pi_s = 0.0;

for (int n=0; n<100000;++n)
{
pi_s = pi_s + (pow(-1,n)/(2*n+1));
}
cout << "Normale " << pi_s*4 << endl <<endl;

return 0;
}
-----------------------------------------------------


OMP (tempo di esecuzione 0.05s)
-----------------------------------------------------

#include <cmath>
#include <iostream>
#include <omp.h>

using namespace std;

int main()
{

float pi_s = 0.0;
#rpragma omp parallel
{
#pragma omp for reduction(+:pi_s) nowait
for (int n=0; n<100000;++n)
{
pi_s = pi_s + (pow(-1,n)/(2*n+1));
}
}
cout << "Normale " << pi_s*4 << endl <<endl;

return 0;
}

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


GMP (tempo di esecuzione 0.1s)
-----------------------------------------------------

#include <cmath>
#include <iostream>
#include <gmp.h>

using namespace std;

int main()
{
mpf_t pi;
mpf_t transi;
mpf_t print;


int sopra;
int sotto;
//float pi_s = 0.0;

mpf_init_set_ui(pi,0);
mpf_init_set_ui(transi,1);
mpf_init_set_ui(print,0);

for (int n=0; n<100000000;++n)
{
sopra = pow(-1,n);
sotto = (2*n)+1;
mpf_set_si(transi,sopra);
mpf_div_ui(transi,transi,sotto);
mpf_add(pi,pi,transi);
mpf_mul_ui(print,pi,4);

//Calcolo eseguito senza GMP
//pi_s = pi_s + (pow(-1,n)/(2*n+1)); 
//cout << "Normale " << pi_s*4 << endl <<endl;
}
gmp_printf("Pi %.Ff\n",print);

return 0;
}



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

L'integrazione di OMP e GMP non e' banale perche' se si uniscono i due esempi sopra riportati il compilatore riporta che la variabile pi has invalid type for ' reduction'

Calcolo in precisione arbitraria di Pi Greco con metodo Gauss-Legendre usando GMP

Per calcolare il valore di Pi Greco e' stato impiegato il metodo di Gauss-Legendre che viene descritto in Wikipedia e a questo Link

per semplicita' il codice e' stato calcolato in Python
------------------------------------------------
import math

a=1
b=1/math.sqrt(2)
t=0.25
p = 1

for s in range(0,3):
    x = (a+b)/2
    y = math.sqrt(a*b)
    t = t - p*(a-x)*(a-x)
    a = x
    b = y
    p =2*p

    pi = ((a+b)*(a+b))/(4*t)

    print(pi)
 ------------------------------------------------

3.1405792505221686
3.141592646213543
3.141592653589794
3.141592653589794

come si vede gia' alla quarta iterazione l'algoritmo esce dal campo di precisione di Python

Per questo motivo e' stata usata la libreria GMP ed il C++ per il calcolo

per compilare l'esempio sotto riportato la stringa e' la seguente

g++  -Wall -O3  -o pi_greco main.cpp -lgmp

Attenzione: nell'uso di GMP si puo' trovarsi di fronte a segmentation fault abbastanza bizzarri che non sono risolvibili usando Gdb. In generale si tratta di errori derivanti da variabili non inizializzate (mpf_init)
-------------------------------------------

#include <cstdlib>
#include <iostream>
#include <iomanip>
#include <gmp.h>

using namespace std;

int main()
{
mpf_t a;
mpf_t b;   
mpf_t t;
mpf_t p;
mpf_t x;
mpf_t y;
mpf_t transi;
mpf_t pi_greco;

mpf_set_default_prec(100000);              //precisione

// a0 = 1
mpf_init_set_ui (a,1);

// b0=1/sqrt(2)
mpf_init (b);
mpf_sqrt_ui (b,2);
mpf_ui_div(b,1,b);

// t0 =1/4
mpf_init_set_ui(t,4);
mpf_ui_div(t,1,t);

// p0 = 1
mpf_init(p);
mpf_set_ui(p,1);

mpf_init(pi_greco);
mpf_init(x);
mpf_init(y);
mpf_init(transi);       
        
for (int f=0;f<2;f++)                   
    {
        cout << f << endl;

        //x(n+1)= [a(n)+b(n)]/2
        mpf_add(x,a,b);
mpf_div_ui(x,x,2);

        //y(n+1) = sqrt[a(n)*b(n)]
        mpf_mul(y,a,b);
        mpf_sqrt (y,y);

        //t(n+1)=t(n)-p(n)*[a(n)-a(n+1)]*[a(n)-x(n+1)]
        mpf_sub(transi,a,x);
        mpf_pow_ui (transi,transi,2);
        mpf_mul(transi,p,transi);
        mpf_sub(t,t,transi);

        //a(n+1) = x(n+1)
        mpf_set(a,x);
               
        //b(n+1) = y(n+1)
        mpf_set(b,y);

        //p(n+1) = 2*p(n)
        mpf_mul_ui(p,p,2);
        
        // pi = {[a(n)+b(n)]*[a(n)+b(n)]}/4*t(n)
        mpf_add(pi_greco,a,b);
        mpf_pow_ui(pi_greco,pi_greco,2);
        mpf_mul_ui(transi,t,4);
        mpf_div(pi_greco,pi_greco,transi);
        
        // stampa il risultato con 20 cifre decimali
        gmp_printf ("%.*Ff \n", 40,pi_greco);
        cout <<endl;
     }

mpf_clear(a);
mpf_clear(b);   
mpf_clear(t);
mpf_clear(p);
mpf_clear(x);
mpf_clear(y);
mpf_clear(transi);
mpf_clear(pi_greco);
return 0;
}

-------------------------------------------
Si riporta di seguito i valori di Pi Greco ad ogni ciclo per un massimo di 1000 valori. Come si osserva gia' dopo 10 cicli viene superata la precisione di 1000 cifre

Ciclo 1
3.14057925052216824831133126897582331177344023751294833564348669334558275803490290782728762155276690054600542214681392392
6603771131652739450259928202137372667423710233461082924341315072787268802546564470528512459515402613483972677323877167912
2732426840482031390025375091177446774967798597085294165910336986203880016912836454106752850429503718532198839982181412863
912905298516564960825702694249040704828539096482646778049496792800995214649820385219718428430430420345783725624986800708
2320226036193321178238482766848584469890521983908330650712909235160675519465414462413539002454813618809982307509964977286
0945401707367366970165459000732289584509236172938129093091466761897251126475917737086756783991546941700851243087017229326
2787616014502438293686095830323226109995848436301136757362913889190508632794887150785016586705141804572465034411751747905
4769631536442522283224388018096944998415551524594221964891234773771474776768528798266179701167763677386772878322823597320
1580750379529048814636645504142030 

Ciclo 2
3.14159264621354228214934443198269577431443722334560279455953948482143476722079526469464344891799130587916462170553518844
2692995943470362111923739681179958736576363907084342931450942394899921183673857971703487633920451060508862477138697445011
975595819591178663703119106652209062422531024241796665754348697348004518821485689267478610281209486160824640075197621208
806439587870791716258988628640907629010918139149013250231997157310590323054416533109635187266992895461083024753433971913
590216260360220405827412847758511214263777802150234914081972370755553638111308997125116036733933640929504739848476056542
804617168825299904305178140800993980593101220914985990549079984548368668619772815061607537106543380426568377560862604710
028374349561731696112951135747426484672197156765386651866428687106800590814836408817914514460647065520971296744821284921
720344206810148741218937154284703784790977421658554307341204200703086089758217520703729558275445098275689876636417983778
2007316946310718306439624939183091860788 

Ciclo 3
3.1415926535897932382795127748018639743812255048354469357873307020263821378389273990314169420434690584298425938468896459
54864005006320361387620067064201557799722915516970065949879812904955325404160658122432921584567412480132658234102425242
52426935034784627342871330628976094908535984940539663050301158950668753595626492715190641080426968106271396859272506557
77383647103787084578459179460310441980410882901610316263215368106231882953875707276225717049376715701039242696961178724
68297792388897328074464470104523028664638809256347199787994658372921417210334354534760789618530384500740608769448526311
60497525882112144066278960327177038166937417479222925842128868808625708523725592929485345308935574361079933702310062997
34979411025447337010283942857551033393532404002081853728714425478248662613970237530603743281096758930221096865760730328
07230351586017812155184446963418523456826552224176908013157755987027979431887869515157263235786247135869497704335546561
5309826339353100069069940892554438783063770913930 

Ciclo 4
3.141592653589793238462643383279502884197114678283648921556617106976026764500643061711006577726598068436361664148276914
16485454070719194016483154461773916689351138620438227940063974546931698156751035142469532101098611678905891188447050048
40090776596878992645170200657786944482215874884883621955402241744349806434167484323444943567596290489661641779152307648
7057640235437621389894094577707569183081654726000204122101060172752158793184361160223821024384600827202314993101818301
6458931744683792684563455883344392627096230917639257443852962863640654199563674088888287525746364238807767091890744306
1596693050772814190050355302838023182354046894783704083764725765267606545185060017656214287346299504077787125420507849
8816395323458883044734431089512463732038595977319790434755545771852082698521409782762852655532500542032553534992266610
96970488727123112218301554150531703404329404312275370200264167678391290552856574222870116270034273962952041947837770191
244768782797615392945520337668843695160923519839477112 

Ciclo 5
3.141592653589793238462643383279502884197169399375105820974944592307816406286208998625628703211672035909744439104150382
53620268912062548831722602181688300582667339844722503666496802561322840114141427517200670102131796896649838756135898733
63134215664300849320442846223244617797140809888885869453795196266152568107315141375849579292331286408679818068363480572
67705700469057209966524902310739676149712910773110507495364911803154947737522943642558604421490091433775926896945313449
827825968763428595957916678812624666318103855903736966421693991108625831384244812413387260931014465915680197497137423850
67249756935679099308129897318592665552743505791045891022240492376337141259853397905526768204019376716993960932751372667
541533303467320056807628017515230977454574006708476024107441318242391183003472414561170642354846396860045026634375831113
99618590688385206236665907589311452199164182618422376237670021444031990005842921760493388970950982068663048012270214220
506475390670669130082388422407589713638444992693 

Ciclo 6
3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230
664709384460955058223172535940812848111745028410270193621252484471023268213360964757747198064531869104652513001491340
5588839210803172323333291749755803274563100575361033879380180768113969946384820292073800055880157131655655454116091736
7526885083086450945682848450384209928567222128844931376229119068712911801907137991678216927948027219089415361922552985
9801855506287714898820522711638405658662397885516317778343608838451175986983344754061347501445992910446021579718854502
0471957544902329417938579109369879229630170688796199507333475267101573269259532026529682068250364326951112679643978548
7269202286585409655597371510536159149755134869808797084458146151256003086993120824899628883650952024633891283348086235
5457839971678511309996024984107193774373370039677361749416101038482033732440027837551488999780875713601471119840859050
76515887310454305599434081033650724454834932852976277999364 

Ciclo 7
3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306
64709384460955058223172535940812848111745028410270193852110555964462294895493038196442881097566593344612847564823378678
31652712019091456485669234603486104543266482133936072602491412737245870066063155881748815209209628292540917143057413534
50440667382513751358565292034707378014393798889366622982253340138274403690736382942540509418318531315363994956407718315
44363343252145872768081739126150762386054529342043033646781684473982919765195990119462357320544935385089619528373019010
375317231932464284885948015032624536588936698124039903107582956066528891986708861709457326630471519235091610969082746708
96684798666105739203571338294861172168169123548084509639332249592932980298335568020884918653678899221443061377827509174
90772979140210763743333962419814109961106955155594716965286044244839339469753269421606409455499598959642733696815862352
0683440365069765290751214023128834675375542410936 

Ciclo 8
3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066
47093844609550582231725359408128481117450284102701938521105559644622948954930381964428810975665933446128475648233786783
165271201909145648566923460348610454326648213393607260249141273724587006606315588174881520920962829254091715364367892590
360011330530548820466521384146951941511609433057270365759591953092186117381932611793105118548074462379962749567351885752
7248912279381830119491298336733624406566430860213949463952247371907021798609437027705392171762931767523846748184676694051
3200056812714526356082778577134275778960917363717872146844090122495343014654958537105079227968914782804866418612955899291
9690811348745019844323375317385085646339351403203780272140712314635426544617219666936810924934065012116028927817479289827
4544904121263955420283220833303191816870251963489867129139945633356295410588944254610943288050001312506092665981545168221
277057553914981371823626270813911105710 

Ciclo 9
3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647
0938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652
71201909145648566923460348610454326648213393607260249141273724587006606315588174881520920962829254091715364367892590360011
33053054882046652138414695194151160943305727036575959195309218611738193261179310511854807446237996274956735188575272489122
79381830119491298336733624406566430860213949463952247371907021798609437027705392171762931767523846748184676694051320005681
27145263560827785771342757789609173637178721468440901224953430146549585371050792279689258923542019956112129021960864034418
15981362977477130996051870721134999999837297804995105973173281609631859502445945534690830264252230825334468503526193118817
10100031378387528865875332083814206171776691473035982534904287554687311595628638823537875937519577818577805321712268066130
019278766111959092164201989 

Ciclo 10
3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647
0938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652
7120190914564856692346034861045432664821339360726024914127372458700660631558817488152092096282925409171536436789259036001
1330530548820466521384146951941511609433057270365759591953092186117381932611793105118548074462379962749567351885752724891
2279381830119491298336733624406566430860213949463952247371907021798609437027705392171762931767523846748184676694051320005
6812714526356082778577134275778960917363717872146844090122495343014654958537105079227968925892354201995611212902196086403
4418159813629774771309960518707211349999998372978049951059731732816096318595024459455346908302642522308253344685035261931
1881710100031378387528865875332083814206171776691473035982534904287554687311595628638823537875937519577818577805321712268
066130019278766111959092164201989 

di seguito il calcolo di Wolfram Alpha