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

 

Debian Testing vs Ubuntu 12.04 in HP DC 5100 SFF

Certo che considerare hardware obsoleto una macchina con un processsore P4 3.2 GHz e 512 Mb di Ram puo' essere fuori luogo ma si deve ricordare che la macchina era in commercio circa 7 anni fa

Il calcolatore HP DC 5100 SFF (dove SFF sta per Small Form Factor) e' decisamente uno dei computer da tavolo che preferisco principalmente per la buona ingegnerizzazione interna e per la presenza di ben 8 porte USB oltre ad una scheda di rete da 1Gb (la scheda video non e' eccezionale)

La prima prova e' stata fatta con una Ubuntu 12.04 dato che pensavo il computer in grado di reggere questa distribuzione e poi volevo farmi un'idea personale di Unity di cui ho sentito sempre parlare un gran male

Ubuntu oramai e' diventata ingestibile su hardware anche un minimo datato ed e' stata subito rimossa senza particolari rimpianti. Una curiosita': se si vogliono installare le applicazioni server dal Software Center (per esempio un server FTP o SSH) e' sostanzialmente impossibile perche' non compaiono..sono dovuto passare da apt-get

Dato che la macchina e' molto recente (almeno rispetto ai miei standard) ho provato a mettere su Debian Testing iniziando dal cd di installazione con Gnome....con sommo rammarico l'installazione non e' riuscita a far partire correttamente il server X (vedere Debian che non gestisce un chipset grafico di almeno 7 anni prima fa pensare su quale [non] sara' il futuro della distribuzione almeno per il desktop)

Usando il cd di installazione di Debian Testing con XFCE+LXDE l'installazione e' arrivata a buon fine quindi mediante apt ho installato GNOME 3. E' un po' spiazzante rispetto al passato ma devo ammettere che non e' poi cosi' male se ci si prende la mano
C'e' una comoda funzione tipo Expose' posizionando il cursore del mouse nell'angolo in alto a sinistra

Expose' su GNOME

Di fatto non esiste piu' un menu' a tendina sostituito da una enorme finestra che sembra essere stata presa a prestito dall'interfaccia dei tablet

Menu' di Gnome

Sul processore P4 640 non funzionano purtroppo le estensione di QEmu-Kvm (a questo punto mi sto chiedendo quali siano i processori che le supportano)

Il processore invece supporta completamente i 64 Bit. In basso l'esecuzione del Live CD 64 Bit di Scientific Linux

Live CD Scientific Linux 6.3 64 Bit

Sqrt e compilazione di math.h

Una cosa curiosa...dopo circa 20 anni mi sono ritrovato a scrivere un programma per la ricerca di numeri primi piu' che altro per valutare l'IDE di Geany e Dev-cpp

Il programma che riporto in calce include la libreria math.h e si compila senza troppi problemi in Dev-cpp mentre portato su Linux e compilato  con

gcc -Wall  -o num_primi num_primi.c

riporta che la funzione sqrt non ha riferimenti
se invece si compila come

g++ -Wall  -o num_primi num_primi.c

curioso no???

la soluzione e' che C++ non ha necessita' di includere la funzione radice quadrata come libreria esterna mentre C si' per cui la giusta stringa di compilazione e'

gcc -Wall -lm  -o num_primi num_primi.c

Geany

Configurazione di compilazione di Geany

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

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define MAX 100000

int indice;
int t,s;
int matrice[MAX];

int calcola(int numero)
{
int i,c;

      c = 0;
      for (i=1;i<=indice;i++)
      {
          if (((numero%matrice[i]) == 0) || (sqrt(matrice[i]) >
numero)) //se il modulo e' zero e' un divisore allora esce subito
             {
             c = 1;
             break;
             }
      }
if (c==0)
         {
             printf("%i\n\r",numero);
             return 1;
             }
             else
             return 0;
}

int main(int argc, char *argv[])
{
  matrice[0] = 1;
  matrice[1] = 2;
  indice = 1;
  for (t=3;t<=MAX;t++)
      {
      if (t%2!=0) // calcola solo se e' un numero dispari
         {
          s = calcola(t); //calcola se e' primo
          if (s==1) //se e' primo inseriscilo nella matrice
                   {
                    indice++;
                    matrice[indice] = t;
                    }
          }
      }
  return 0;
}

Switch di compilazione in Dev-cpp

Per includere le librerie durante la compilazione in Dev-C++, oltre ad avere gia' incluso il DevPack per esempio di Sqlite, si deve indicare anche in Progetto/Opzioni Progetto/Parametri lo switch appropriato

Nell'esempio della figura e' incluso lo switch -lsqlite3 per Sqlite

Esempio Mandelbrot con PngWriter/C

In questo esempio l'insieme di Mandelbrot e' stato generato con la libreria PngWriter che permette di creare dei file png.

La libreria e' gia' pacchettizzata in Debian ed e' di facile utilizzo perche' si deve solo indicare le dimensioni del file png da creare e le coordinate del punto da disegnare con il colore in formato RGB

Per compilare il file di esempio deve essere utilizzata la seguente stringa

g++ mand_png.c -o mand_png `freetype-config --cflags` -I/usr/local/include  -L/usr/local/lib -lpng -lpngwriter -lz -lfreetype

I file del progetto sono reperibili a questo Link


Per visualizzare il file png da linea di comando puo' essere utile installare il pacchetto ImageMagick e digitare

display mandelbrot.png

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

#include <pngwriter.h>

#define SCREEN_WIDTH 1024
#define SCREEN_HEIGHT 768

float re_min = -2.0;
float im_min = -1.2;

float re_max = 1.0;
float im_max = 1.2;

double iterazioni = 255;


double r;

float a,b;
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,"mandelbrot.png");

float re_factor = (re_max-re_min);
float im_factor = (im_max-im_min);

for (i=0;i<SCREEN_HEIGHT;i++)
{
for (j=0;j<SCREEN_WIDTH;j++)
{
a = re_min+(j*re_factor/SCREEN_WIDTH);
b = im_min+(i*im_factor/SCREEN_HEIGHT);

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

}

Esempio GUI con RAD C++ in Windows

Un ulteriore esempio di Gui utilizzando stavolta delle librerie Rad-C++ che sono disponibili solo per Windows

Il programma e' stato compilato usando Dev-C++ (in versione portable per altro)
 

Dev-C++ Schermata principale

Per utilizzare la libreria deve essere aggiunto l'apposito DevPak  dal menu Strumenti/Package Manager

Finestra di Package Manager

-----------------------------------------------------------
#include <radc++.h>

Form form1("Progress Bar",100,100,215,150);

ProgressBar pgb1(-1,10,10,190, 20,form1,RCP_HORIZONTAL);
Label label1("0",-1,90,40 ,30,20,form1);
Track track(AUTO_ID,10,70,190,30,form1);



FormProcedure proc(FormProcArgs) {
    ON_CLOSE() Application.close();
   
    ON_TRACK_CHANGE(track) {      
        int pos = track.position;
        label1.caption = str(pos);
        pgb1.percent=pos;
    }
    return 0;
}

rad_main()
    form1.procedure = proc;
    track.minRange = 0;
    track.maxRange = 100;
rad_end()

FreeDos in Qemu

Per installare FreeDos in QEmu la procedura e' abbastanza semplice

prima si crea un file della dimensione voluta del futuro C:\ (in questo caso 400 Mb) e si scarica il file fdfullcd.iso dal sito di FreeDos

dd if=/dev/zero of=freedosfile.img bs=1024 count=400000

successivamente si lancia l'emulatore indicando il supporto di installazione nel cdrom e specificando di fare il boot da D

qemu -cdrom fdfullcd.iso -hda freedosfile.img -boot d

si seguono quindi le istruzioni a schermo

Per avviare la macchina una volta terminato il processo e' sufficiente

qemu -hda freedosfile.img -boot c

Uno dei problemi con QEmu/FreeDos e' quello di non poter indicare in modo semplice una cartella condivisa tra il sistema Guest e quello Host in modo da inserire i file di installazione dei futuri programmi

Leggendo su Internet molti settano QEmu per gestire una rete virtuale ed effettuano un collegamento di rete tra Host e Guest

per me e' stata piu' semplice inserire i file dei programmi su una chiavetta USB e poi creare una immagine img della chiavetta mediante il seguente comando

dd if=/dev/sdb of=/home/luca/qemu/usb.img bs=1024

a questo punto si puo' lanciare QEmu indicando nel disco D: l'immagine della chiavetta

qemu -hda freedos.img -hdb usb.img -boot c

Prompt di FreeDos

Nell'immagine di FreeDos completa sono inclusi molti software compreso OpenGem. Prima di poterlo utilizzare deve essere pero' installato. Per farlo si deve andare nella cartella FDos/OpenGem e far partire l'installazione

Open Gem

Ovviamente poter usare nell'emulatore un programma che e' stato tuo compagno per almeno 5 anni non ha prezzo e quindi via con l'installazione di Turbo Pascal 5

 Per il funzionamento del Turbo Pascal 5 e' importante selezionare l'opzione 1 al boot di FreeDos

 Dopo l'installazione dei dischetti il rassicurante schermo blu e bianco della IDE di Turbo Pascal

Fa abbastanza impressione vedere funzionare dopo tanto tempo l'esempio BGIDemo

BGIDemo Turbo Pascal

Quasi quasi provo ad installare Fractint.........................