algoritmo pigreco

[gurutech]

ciao a tutti,
sto cercando di calcolare pigreco, e per ora mi è venuto in mente solo quello che vedete qui sotto.
Probabilmente è una delle procedure computazionalmente più imprecise, perchè non riesco ad andare oltre la 20ma iterazione, poi diventa incoerente.
Conoscete altri metodi migliori?

#!/usr/bin/perl

$l=1;

for ($i=0;$i<20;$i++) {
$m=1 - sqrt (1 - ($l**2)/4);
$l=sqrt($m*2);
$pigreco=6*(2**$i)*$l;
print "$i * $pigreco\n";
}

[gurutech]

dimenticavo, i valori che calcolo sono
0 * 3.10582854123025
1 * 3.13262861328124
2 * 3.13935020304687
3 * 3.14103195089053
4 * 3.14145247228534
5 * 3.14155760791162
6 * 3.14158389214894
7 * 3.14159046323676
8 * 3.14159210604305
9 * 3.14159251658815
10 * 3.14159261864079
11 * 3.14159264532122
12 * 3.14159264532122
13 * 3.14159264532122
14 * 3.14159264532122
15 * 3.14159230381174
16 * 3.14159230381174
17 * 3.14158683965504
18 * 3.14158683965504
19 * 3.14167426502176

[gurutech]

da un punto di vista matematico, ho trovato questo
http://www.cecm.sfu.ca/pi/formulas.ps

e secondo questo
http://www.cecm.sfu.ca/projects/ISC/.../C/pi10000.txt

l'algoritmo che ho fatto ci azzecca alla 7ma cifra fino al 14mo ciclo poi comincia a sbagliare

[gurutech]

oops
secondo questa fonte
http://numbers.computation.free.fr/Constants/Pi/pi.html

sono rimasto tra i babilonesi e archimede...

[kingv]

Originariamente inviato da gurutech
oops
secondo questa fonte
http://numbers.computation.free.fr/Constants/Pi/pi.html

sono rimasto tra i babilonesi e archimede...

potresti in alternativa misurare la circonferenza di un cerchio e dividerlo per il diametro, solo che devi essere mooooolto preciso .

[gurutech]

Originariamente inviato da kingv
potresti in alternativa misurare la circonferenza di un cerchio e dividerlo per il diametro, solo che devi essere mooooolto preciso .

è un'idea. si rischia di finire tra due estremi:

-A- Stato dell'Indiana : nel 1897 hanno provato a legiferare pi=3.2
http://faqs.jmas.co.jp/FAQs/sci-math-faq/indianabill
-B- Ramanujan

[kingv]

Originariamente inviato da gurutech


-A- Stato dell'Indiana : nel 1897 hanno provato a legiferare pi=3.2
http://faqs.jmas.co.jp/FAQs/sci-math-faq/indianabill

io spero che sia una bufala perche' in caso contrario tanti miei dubbi sulla sanità mentale della classe politica americana sarebbero giustificati.



posso chiederti perche' stai cercandoti un algoritmo "implementabile"? Vuoi farci su qualche tipo di benchmark?

[gurutech]

Originariamente inviato da kingv
io spero che sia una bufala perche' in caso contrario tanti miei dubbi sulla sanità mentale della classe politica americana sarebbero giustificati.

no, no giustificati pure

posso chiederti perche' stai cercandoti un algoritmo "implementabile"? Vuoi farci su qualche tipo di benchmark?

spippolamento? insanità mentale?

The amazing number pi Peter Borwein
[...]
Why compute the digits of PI? Sometimes it is necessary to do so, though hardly ever more than 6 or so digits that Archimedes computed several thousand years ago are needed for physical applications. Even far fetched computations like the volume of a spherical universe only require a few dozen digits. There is also the "Everest Hypothesis" ("because its there").
[...]

[gurutech]

ora, nonostante sia l'una passata, non sono riuscito a resistere e ho costruito un'altro piccolo algoritmo per trovare pi per eccesso (quello di prima era per difetto):

#!/usr/bin/perl

$R=2/(sqrt(3));
for ($i=1;$i<20;$i++) {
$l=sqrt(($R-1)/($R+1));
$R=sqrt((2*$R)/($R+1));
$pigreco=$l*6*(2**$i);
print "$pigreco\n";
}

pigreco sta matematicamente in mezzo tra i due algoritmi...

3.21539030917347
3.15965994209751
3.14608621513146
3.14271459964532
3.14187304997995
3.14166274705598
3.14161017660122
3.14159703431985
3.14159374879501
3.14159292765212
3.14159272140683
3.14159267301376
3.14159268159282
3.14159264371695
3.14159230341067
3.14159230371147
3.14158683962997
3.14158683964877
3.14167426502019

[Frank1962]

ma che linguaggio stai usando? ...php ?

[cionci]

Perl

[gurutech]

ecco, ho messo tutto qui
http://www.gurutech.it/polimi/pigreco.pdf
compreso come ho fatto a tirar fuori le formule.
Adesso è meglio che vada a studiare sul serio, vi lascio in omaggio questo link dove trovate la stessa cosa che ho fatto io in java

http://www.math.utah.edu/~alfeld/Arc...rchimedes.html

[Mixmar]

Se vi interessano solo le cifre, consiglio:

http://mathworld.wolfram.com/BBPFormula.html

Questo sì che è avere tempo da perd... ehm, dedicare con profitto a ricerche superiori!

[gurutech]

Originariamente inviato da Mixmar
Questo sì che è avere tempo da perd... ehm, dedicare con profitto a ricerche superiori!

ok ok, è solo una curiosità!

però leggiucchiando qua è là ho scoperto che sul calcolo di pi ci si può anche lucrare.
mi spiego: per calcolare 1 miliardo di cifre di pigreco in un tempo umano e senza errori con un computer, bisogna inventarsi un metodi di fare conti molto in fretta che poi può essere applicato a un'infinità di altri campi.
In questo senso il modo che abbiamo tutti imparato alle elementari per fare le moltiplicazioni, applicato ad un computer è una schifezza.

[/\/\@®¢Ø]

Una serie molto efficace e' quella di Borwein e Borwein (1987), che aggiunge ad ogni iterazione 25 (!) cifre esatte.
La formula non e' proprio elegantissima, ma ci si puo' accontentare...

1/pi = 12*somma(n da 0 a inf.)
( (-1)^n *(6n!)[ 212175710912 sqrt(61) + 1657145277365+n(13773980892672 sqrt(61) + 107578229802750)] / (n)!^3 (3n)! [5280 (236674 + 30303 sqrt(61)]^(3n+3/2)

[/\/\@®¢Ø]

Originariamente inviato da gurutech
ora, nonostante sia l'una passata, non sono riuscito a resistere e ho costruito un'altro piccolo algoritmo per trovare pi per eccesso (quello di prima era per difetto):

#!/usr/bin/perl

$R=2/(sqrt(3));
for ($i=1;$i<20;$i++) {
$l=sqrt(($R-1)/($R+1));
$R=sqrt((2*$R)/($R+1));
$pigreco=$l*6*(2**$i);
print "$pigreco\n";
}

pigreco sta matematicamente in mezzo tra i due algoritmi...

3.21539030917347
3.15965994209751
3.14608621513146
3.14271459964532
3.14187304997995
3.14166274705598
3.14161017660122
3.14159703431985
3.14159374879501
3.14159292765212
3.14159272140683
3.14159267301376
3.14159268159282
3.14159264371695
3.14159230341067
3.14159230371147
3.14158683962997
3.14158683964877
3.14167426502019

Devi fare attenzione ad una cosa importante, ovvero alla precisione dei tipi che usi per i conti. In questi casi e' praticamente necessario utilizzare numeri di precisione arbitraria. A meno di ricorrere a librerie pero', la maggior parte dei linguaggi supporta precisione arbitraria _al piu'_ sugli interi. Il fatto di usare funzioni come sqrt nel tuo caso probabilmente ammazza la precisione. Meglio piuttosto usare un algoritmo che operi sui razionali, anche se magari piu' lento, come il seguente


Il risultato che ho ottenuto io (in haskell) e' il seguente:

2.0
2.6666666666666665
2.933333333333333
3.0476190476190474
3.0984126984126985
3.1215007215007216
3.1321567321567323
3.137129537129537
3.139469680646151
3.140578169680337
3.1411060216013778
3.1413584725201362
3.1414796489611403
3.141537993173476
3.141566159344948
3.1415797881375958
3.1415863960370616
3.1415896055882304
3.141591166991502
3.141591927675147
3.1415922987403397
3.1415924799582244
3.141592568553635
3.1415926119088358
3.141592633144036
3.141592643553448
3.141592648659952
3.141592651166781
3.141592652398206
3.1415926530034826
3.1415926533011596
3.1415926534476357
3.141592653519747
3.1415926535552643
3.141592653572766
3.141592653581393
3.141592653585648
3.1415926535877468
3.141592653588783
3.141592653589294
3.1415926535895466
3.1415926535896714
3.141592653589733
3.1415926535897634
3.1415926535897785
3.141592653589786
3.1415926535897896
3.1415926535897913
3.1415926535897922
3.1415926535897927
3.141592653589793

come vedi converge in modo migliore, a differenza del tuo codice, che dopo un po' "sbarella" (ci sono comunque alcuni errori dovuti alla conversione da razionale a virgola mobile che provero' ad evitare)

Qui comunque trovi diverse formule, qualcuna puo' fare al caso tuo.

[gurutech]

Originariamente inviato da /\/\@®¢Ø
A meno di ricorrere a librerie pero', la maggior parte dei linguaggi supporta precisione arbitraria _al piu'_ sugli interi.

potrei provare a passare il tutto a bc, dici che funzionerà?

Il fatto di usare funzioni come sqrt nel tuo caso probabilmente ammazza la precisione. Meglio piuttosto usare un algoritmo che operi sui razionali, anche se magari piu' lento,

be' ho provato anche con leibnitz
pi/4 = 1 - 1/3 + 1/5 - 1/7 ...
ma per tirar fuori 7 cifre ci mette un milione di iterazioni :eheh:

[/\/\@®¢Ø]

Originariamente inviato da gurutech
potrei provare a passare il tutto a bc, dici che funzionerà?

boh non l'ho mai usato.




be' ho provato anche con leibnitz
pi/4 = 1 - 1/3 + 1/5 - 1/7 ...
ma per tirar fuori 7 cifre ci mette un milione di iterazioni :eheh: [/quote]
Beh, non troppo semplice
La formula che ho usato io ad esempio va gia' abbastanza bene...
forse ancora meglio una delle seguenti, tutte dovute a Ramanujan e che trovi nel link che ho riportato:

[/\/\@®¢Ø]

come non detto quelle tre forse non sono cosi' semplici (uhm... termine a pedice... )
la' puoi comuinque trovarne di opportune.

[gurutech]

ecco con bc(1) è molto meglio

bash-2.05a$ time bc pigreco.bc 
bc 1.06
Copyright 1991-1994, 1997, 1998, 2000 Free Software Foundation, Inc.
This is free software with ABSOLUTELY NO WARRANTY.
For details type `warranty'. 
3.1415926535897932384626433832795028841971693993751058209749445*02977\
93979003832755187039605537355163056713811295792041015968426903519036\
19418561079235059815597867847277845963233428011049306481123254672547\
06311502208770628628413205845869436023049732563569828695132669875998\
199952912012796343769788055552

real    0m12.504s
user    0m12.198s
sys     0m0.003s

sono giuste le prime 61 cifre (fino all'*). il codice è questo:

scale=300
l=1
for (i=0; i<100; i++) {
m = 1 - sqrt(1-(l^2)/4)
l=sqrt(2*m)
pigreco=6*(2^i)*l
}
print pigreco
print "\n"

quit

(1)

man bc
...
NAME
bc - An arbitrary precision calculator language
...
DESCRIPTION
bc is a language that supports arbitrary precision numbers with interactive execution of
statements. There are some similarities in the syntax to the C programming language. A
standard math library is available by command line option.
...