Pensavo da tempo ad una nuova formula per calcolare in modo più efficiente le cifre decimali del numero Pi Greco. Prendendo spunto dalla pagina Wikipedia Calcolo Pi Greco e da altre fonti, sembra che la formula più valida sia la Formula di Bailey–Borwein–Plouffe, scritta in questa forma:
Così ho pensato a qualcosa di nuovo: introdurre un nuovo parametro e calibrarlo, per minimizzare ulteriormente l'errore (o meglio, ricalibrare tutti i parametri). Scrivendo un programma in C ho applicato il metodo di Newton-Raphson per sistemi di equazioni, non essendo semplice come il caso scalare (ovvero equazione singola) faccio prima una precisazione:
- si definiscono K derivate parziali quanti sono i K parametri (dal punto di vista computazionale, limite del rapporto incrementale con incremento molto piccolo, es. 10-15)
- un ciclo esegue N volte in successione (non assieme!!) le seguenti operazioni (prima ovviamente il solito controllo denominatore diverso da zero):
- Ak+1=Ak-(f(A,B,...)-PI)/(df/dA)
- Bk+1=Bk-(f(A,B,...)-PI)/(df/dB)
- e così via fino all'ultimo parametro, in questo caso Ik+1 quindi 9 parametri da calibrare
- infine si ottengono tutti i parametri calibrati e si stampa il risultato
La formula BBP è un algoritmo iterativo che converge in modo spaventosamente veloce, vista l'enormità di cifre decimali l'ho testata solamente con un numero di iterazioni N=5, davvero molto piccolo (si vede dunque l'enorme efficacia!) e per confronto ho dovuto recuperare 50 cifre decimali di Pi Greco (vedi Prime 100000 cifre di Pi greco) dato che la libreria <math.h> contiene la costante M_PI con "sole" 20 cifre decimali, troppo limitante per i nostri ambiziosi scopi 🙂
Se la formula BBP aveva un'enorme efficacia, grazie alla mia nuova calibrazione con un parametro aggiuntivo, l'efficacia è aumentata a dismisura. Giusto per intenderci, ecco i risultati (nello screenshot con N=10, mentre nell'elenco puntato con N=5):
- PiGreco 50cifre = 3.14159265358979311599796346854418516159057617187500
- BBP_originale = 3.14159264546033645260081357264425605535507202148438
- BBP_Giulio_M = 3.14159265356988193218512606108561158180236816406250
Li ho incolonnati in modo da vedere ad occhio l'accuratezza. Ripeto che ho eseguito SOLO 5 iterazioni altrimenti per assurdo con N=10 la mia formula era già accurata oltre la 50° cifra decimale di Pi Greco!!
In grassetto l'ultima cifra corretta, come si vede con N=5 la formula BBP originale arriva alla settima cifra, la mia variante arriva alla decima cifra. Con N=10 BBP originale arriva alla quattordicesima cifra e come detto la mia formula invece oltre la cinquantesima cifra corretta, per un'accuratezza maggiore dovrei scrivere un altro programma ad-hoc per il calcolo.
Dato che sono bravo, vi dico quali sono tutti i nuovi parametri (come li ho definiti io) e in seguito tutto il codice del programma, scritto in un centinaio di righe.
Dalla formula scritta all'inizio, ho chiamato A,B,C quelli che hanno rispettivamente coefficiente pari a 120, 151, 47; D,E,F,G,H invece erano rispettivamente 512, 1024, 712, 194, 15; nel mio caso ho introdotto un nuovo coefficiente I che moltiplica k5
Seguendo questa convenzione, i miei coefficienti sono (per semplificare, ne ha tarati diversi come uguali a 1, calibrando il valore degli altri):
- A=1
- B=1
- C= 2,90120737
- D=1
- E=1
- F=1
- G=1
- H=0,93906863
- I=1
Ecco infine tutto il codice in C 🙂
//migliorare BBP con nuovo parametro, partendo da formula in base10
// https://it.wikipedia.org/wiki/Formula_di_Bailey%E2%80%93Borwein%E2%80%93Plouffe
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define N 5
#define EPS 0.000000000000001
#define PI 3.14159265358979323846264338327950288419716939937510 //50cifre
#define DEN (D*i*i*i*i+E*i*i*i+F*i*i+G*i+H+I*i*i*i*i*i)
double bbp(unsigned char n, double A, double B, double C, double D, double E, double F, double G, double H, double I);
double dfA(unsigned char i, double A, double B, double C, double D, double E, double F, double G, double H, double I);
double dfB(unsigned char i, double A, double B, double C, double D, double E, double F, double G, double H, double I);
double dfC(unsigned char i, double A, double B, double C, double D, double E, double F, double G, double H, double I);
double dfD(unsigned char i, double A, double B, double C, double D, double E, double F, double G, double H, double I);
double dfE(unsigned char i, double A, double B, double C, double D, double E, double F, double G, double H, double I);
double dfF(unsigned char i, double A, double B, double C, double D, double E, double F, double G, double H, double I);
double dfG(unsigned char i, double A, double B, double C, double D, double E, double F, double G, double H, double I);
double dfH(unsigned char i, double A, double B, double C, double D, double E, double F, double G, double H, double I);
double dfI(unsigned char i, double A, double B, double C, double D, double E, double F, double G, double H, double I);
int main() {
register double sum=0;
register unsigned char i;
register double A,B,C,D,E,F,G,H,I;
//INIT
A=1.0;B=1.0;C=1.0;D=1.0;E=1.0;F=1.0;G=1.0;H=1.0;I=1.0;
//sistema lineare per determinare incognite: f()=bbp()-PI=0
for(i=0;i<N;i++){
if(DEN!=0){
if(dfA(N,A,B,C,D,E,F,G,H,I)!=0){
A=A-(bbp(N,A,B,C,D,E,F,G,H,I)-PI)/dfA(N,A,B,C,D,E,F,G,H,I);
}
if(dfB(N,A,B,C,D,E,F,G,H,I)!=0){
B=B-(bbp(N,A,B,C,D,E,F,G,H,I)-PI)/dfB(N,A,B,C,D,E,F,G,H,I);
}
if(dfC(N,A,B,C,D,E,F,G,H,I)!=0){
C=C-(bbp(N,A,B,C,D,E,F,G,H,I)-PI)/dfC(N,A,B,C,D,E,F,G,H,I);
}
if(dfD(N,A,B,C,D,E,F,G,H,I)!=0){
D=D-(bbp(N,A,B,C,D,E,F,G,H,I)-PI)/dfD(N,A,B,C,D,E,F,G,H,I);
}
if(dfE(N,A,B,C,D,E,F,G,H,I)!=0){
E=E-(bbp(N,A,B,C,D,E,F,G,H,I)-PI)/dfE(N,A,B,C,D,E,F,G,H,I);
}
if(dfF(N,A,B,C,D,E,F,G,H,I)!=0){
F=F-(bbp(N,A,B,C,D,E,F,G,H,I)-PI)/dfF(N,A,B,C,D,E,F,G,H,I);
}
if(dfG(N,A,B,C,D,E,F,G,H,I)!=0){
G=G-(bbp(N,A,B,C,D,E,F,G,H,I)-PI)/dfG(N,A,B,C,D,E,F,G,H,I);
}
if(dfH(N,A,B,C,D,E,F,G,H,I)!=0){
H=H-(bbp(N,A,B,C,D,E,F,G,H,I)-PI)/dfH(N,A,B,C,D,E,F,G,H,I);
}
if(dfI(N,A,B,C,D,E,F,G,H,I)!=0){
I=I-(bbp(N,A,B,C,D,E,F,G,H,I)-PI)/dfI(N,A,B,C,D,E,F,G,H,I);
}
}
}
printf("A=%.8f|B=%.8f|C=%.8f|D=%.8f|E=%.8f|F=%.8f|G=%.8f|H=%.8f|I=%.8f\n",A,B,C,D,E,F,G,H,I);
printf("bbpGM=%.50f\n",bbp(N,A,B,C,D,E,F,G,H,I));
printf("bbp_o=%.50f\n",bbp(N,120,151,47,512,1024,712,194,15,0));
printf("___PI=%.50f\n",PI);
return 0;
}
double bbp(unsigned char n, double A, double B, double C, double D, double E, double F, double G, double H, double I){
/*return 1.0/pow(16,i)*(A*i*i+B*i+C)/(D*i*i*i*i+E*i*i*i+F*i*i+G*i+H+I*i*i*i*i*i);*/
double ris=0;
unsigned char i;
for(i=0;i<n;i++){
ris+=1.0/pow(16,i)*(A*i*i+B*i+C)/(D*i*i*i*i+E*i*i*i+F*i*i+G*i+H+I*i*i*i*i*i);
}
return ris;
}
double dfA(unsigned char i, double A, double B, double C, double D, double E, double F, double G, double H, double I){
return (bbp(i,A+EPS,B,C,D,E,F,G,H,I)-bbp(i,A,B,C,D,E,F,G,H,I))/EPS;
}
double dfB(unsigned char i, double A, double B, double C, double D, double E, double F, double G, double H, double I){
return (bbp(i,A,B+EPS,C,D,E,F,G,H,I)-bbp(i,A,B,C,D,E,F,G,H,I))/EPS;
}
double dfC(unsigned char i, double A, double B, double C, double D, double E, double F, double G, double H, double I){
return (bbp(i,A,B,C+EPS,D,E,F,G,H,I)-bbp(i,A,B,C,D,E,F,G,H,I))/EPS;
}
double dfD(unsigned char i, double A, double B, double C, double D, double E, double F, double G, double H, double I){
return (bbp(i,A,B,C,D+EPS,E,F,G,H,I)-bbp(i,A,B,C,D,E,F,G,H,I))/EPS;
}
double dfE(unsigned char i, double A, double B, double C, double D, double E, double F, double G, double H, double I){
return (bbp(i,A,B,C,D,E+EPS,F,G,H,I)-bbp(i,A,B,C,D,E,F,G,H,I))/EPS;
}
double dfF(unsigned char i, double A, double B, double C, double D, double E, double F, double G, double H, double I){
return (bbp(i,A,B,C,D,E,F+EPS,G,H,I)-bbp(i,A,B,C,D,E,F,G,H,I))/EPS;
}
double dfG(unsigned char i, double A, double B, double C, double D, double E, double F, double G, double H, double I){
return (bbp(i,A,B,C,D,E,F,G+EPS,H,I)-bbp(i,A,B,C,D,E,F,G,H,I))/EPS;
}
double dfH(unsigned char i, double A, double B, double C, double D, double E, double F, double G, double H, double I){
return (bbp(i,A,B,C,D,E,F,G,H+EPS,I)-bbp(i,A,B,C,D,E,F,G,H,I))/EPS;
}
double dfI(unsigned char i, double A, double B, double C, double D, double E, double F, double G, double H, double I){
return (bbp(i,A,B,C,D,E,F,G,H,I+EPS)-bbp(i,A,B,C,D,E,F,G,H,I))/EPS;
}
Pensandoci alcune funzioni le potevo anche compattare, ma va bene così 🙂 per il caso sistema di equazioni, un'alternativa al metodo di Newton-Raphson, tendenzialmente più semplice da gestire, sarebbe il caso della bisezione, di cui avevo parlato qui con un esempio (che tra l'altro avevo inizialmente difficoltà ad implementare e la soluzione è arrivata grazie ad uno spunto di ChatGPT 😁 ): Metodo della bisezione per sistemi di equazioni
L'efficienza di questi algoritmi è importantissima quando si tratta di voler calcolare molte cifre, attualmente il record mondiale, appartenente a Google Cloud (e non a me) ha raggiunto 100.000 miliardi di cifre, nel 2022 (approfondimento); una curiosità interessante: <<la lettura di tutte le 100 trilioni di cifre ad alta voce alla velocità di una al secondo richiederebbe più di 3,1 milioni di anni>>.
Cosa ne pensate? Vi piace l'idea? Ovviamente il concetto è estendibile altrove, per migliorare ancora un'ipotetica nuova formula più efficiente! 🙂
