venerdì 21 settembre 2012

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
 

Nessun commento:

Posta un commento

Perche' investire su Unix

 Un libro trovato nel libero scambio alla Coop su cio' che poteva essere e non e' stato...interessante la storia su Unix del primo c...