Dopo il caso interessante Calcolo lista numeri primi fino ad N, vediamo una lista di soluzioni numeriche per trovare le cifre decimali di pi greco. Vediamo alcuni degli algoritmi, con l'implementazione in linguaggio C (migliore efficienza), poi commentiamo i risultati.
- metodi geometrici: l'idea è quella di creare una griglia NxN, aggiungere 1 se la cella i,j si trova all'interno del cerchio (equazione matematica), tramite la condizione
if(N*N-i*i-j*j>0){sum+=1;}
. L'area del cerchio è A=pi*D4, da questa formula posso ricavare pi (pi greco) quindi tramite il valore approssimato dell'area (cerchio inscritto nella griglia NxN, diametro N). Non riporto tutto il codice, si può leggermente ottimizzare calcolando pi/4 di un quadrante (anziché tutta la griglia, vista la simmetria!) e infine moltiplicare per 4. Comunque sia è un metodo estremamente inefficiente
- integrale di Gauss: dal caso Benchmark performance C, Python, MATLAB - GNU Octave, JavaScript, PHP, in cui (per ben altri scopi!) avevamo preso come riferimento di esempio la soluzione numerica di una Gaussiana, il risultato è comunque funzione di pi greco. Anche questo risulta comunque estremamente inefficiente rispetto ad altre soluzioni che vediamo dopo
- sommatorie di termini: qui ci sono soluzioni davvero interessanti, ad esempio:
- funzione continua generalizzata
- serie di Gregory-Leibniz
- serie di Nilakantha
- formula di Bailey–Borwein–Plouffe (BBP)
Vediamo uno ad uno questi casi, per la spiegazione e approfondimenti la maggior parte di essi si trovano nella pagina Wikipedia Calcolo pi greco.
Frazione continua
Il codice è il seguente:
#include <stdio.h>
#include <stdlib.h>
#define N 20
int main() {
register double sum=3;
register short i;
for(i=N;i>0;i--){
sum=6+(2*i-1)*(2*i-1)/sum;
}
sum-=3;
printf("somma=%.20f\n",sum);
return 0;
}
Serie di Gregory-Leibniz
Qui si poteva optare per diverse scelte, applicare una struttura if-else con una formula diversa per il caso pari-dispari, oppure usare una funzione (-1)N oppure creare un'espressione unica, facendo poi un salto di 2 del contatore anziché di 1; è stata scelta la prima soluzione, struttura if-else.
Il codice è il seguente:
#include <stdio.h>
#include <stdlib.h>
#define N 1000
int main() {
register double sum=0;
register short i;
for(i=0;i<N;i++){
if(i%2==0){
sum+=1.0/(2*i+1);
}else{
sum+=-1.0/(2*i+1);
}
}
sum*=4;
printf("somma=%.20f\n",sum);
return 0;
}
Serie di Nilakantha
Il codice è il seguente:
#include <stdio.h>
#include <stdlib.h>
#define N 1000
int main() {
register double sum=0;
register short i;
for(i=1;i<=N;i+=2){
sum+=1.0/(i*(i+1.0)*(2.0*i+1.0))-1.0/((i+1.0)*(i+2.0)*(2.0*i+3.0));
}
sum+=3;
printf("somma=%.20f\n",sum);
return 0;
}
Serie di Bailey–Borwein–Plouffe (Formula BBP)
Il codice è il seguente:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define N 20
int main() {
register double sum=0;
register unsigned short i;
for(i=0;i<N;i++){
sum+=1.0/pow(16,i)*(4.0/(8*i+1)-2.0/(8*i+4)-1.0/(8*i+5)-1.0/(8*i+6));
}
printf("somma=%.20f\n",sum);
return 0;
}
Risultati, comparazione algoritmi
Vediamo ora l'aspetto più interessante, un confronto e commento di questi risultati:
indichiamo con "N" le iterazioni totali e "n" l'n-esima cifra corretta nell'approssimazione di pi greco (es. N=20 e n=4 significa che con 20 iterazioni abbiamo un'accuratezza di 4 cifre decimali di pi greco).
- frazione continua: efficienza complessivamente buona, terzo gradino del podio 🥉
N=20, n=4
N=50, n=5
N=100, n=6
N=200, n=7
- serie di Gregory-Leibniz: bocciata, anche su Wikipedia possiamo leggere <<la sua velocità di convergenza è troppo lenta perché sia di interesse pratico>>
N=20, n=1
N=50, n=1
N=100, n=1
N=200, n=1
N=1000, n=2
- serie di Nilakantha: risultato molto valido, medaglia d'argento 🥈
N=20, n=4
N=50, n=5
N=100, n=6
N=200, n=7
N=1000, n=9
- Formula di Bailey–Borwein–Plouffe: abbiamo il vincitore! BBP vince nettamente la medaglia d'oro! 🥇
N=20, n=15
Questi sono i risultati del benchmark fra alcuni dei principali algoritmi considerati. Per il confronto ha senso considerare l'accuratezza nelle cifre decimali piuttosto che il tempo di esecuzione del codice, dato che almeno come ordine di grandezza questo è fondamentalmente simile a parità di numero di iterazioni del ciclo. La formula BBP mostra 15 cifre decimali corrette con appena N=20 ovvero 20 iterazioni totali (mi sono fermato a N=20 per non dover fare una verifica a mano, visto il numero di cifre 😅); per confronto, il secondo posto spetta a Nilakantha che arriva a 9 cifre decimali corrette con ben 1000 iterazioni!