Implementazione della formula Bailey-Borwein-Plouffe in C++?

Question

EDIT: Il requisito era vago e invece di calcolare l'n-esima cifra di pi desideravano solo l'ennesima cifra che non andasse oltre la limitazione dei float, quindi il metodo della forza bruta ha funzionato per i requisiti.

Ho bisogno di calcolare PI la n-esima cifra e volevo provare ad usare lo BBP formula ma ho delle difficoltà. L'equazione che ho digitato non sembra darmi PI in modo corretto.

(1/pow(16,n))((4/(8 * n + 1)) - (2/(8 * n + 4)) - (1/(8 * n + 5)) - (1/(8 * n + 6))) 

ho avuto successo con il modo in cui la forza bruta di trovare PI ma che è solo in modo accurato e trovare il numero n-esimo è difficile.

(4 - (4/3) + (4/5) - (4/7)...) 

volevo sapere se qualcuno ha avuto una migliore idea di come fare questo o forse aiutare con la mia equazione BBP su quello che ho incasinato?

Grazie,
LF4

Funzionale ma non molto preciso fino a quando un paio di iterazioni in e poi si deve disreguard ultimi.

#include <iostream> 

using namespace std; 

int main() 
{ 
    int loop_num = 0; 
    cout << "How many digits of pi do you want?: "; 
    cin >> loop_num; 

    double my_pi = 4.0; 
    bool add_check = false; 
    int den = 3; 
    for (int i = 0; i < loop_num; i++) 
    { 
     if (add_check) 
     { 
      my_pi += (4.0/den); 
      add_check = false; 
      den += 2; 
     } 
     else 
     { 
      my_pi -= (4.0/den); 
      add_check = true; 
      den += 2; 
     } 
    } 
    cout << "Calculated PI is: " << my_pi << endl; 
    system("pause"); 

    return 0; 
} 

Quello che spero sarebbe un programma migliore.

#include <iostream> 
#include <cmath> 

using namespace std; 

const double PI_BASE = 16.0; 

int main() 
{ 
    int loop_num = 0; 
    cout << "How many digits of pi do you want?: "; 
    cin >> loop_num; 

    double my_pi = 0.0; 
    for (int i = 0; i <= loop_num; i++) 
    { 
     my_pi += (1.0/pow(PI_BASE,i))((4.0/(8.0 * i + 1.0)) - 
              (2.0/(8.0 * i + 4.0)) - 
              (1.0/(8.0 * i + 5.0)) - 
              (1.0/(8.0 * i + 6.0))); 
    } 
    cout << "Calculated PI is: " << my_pi << endl; 
    system("pause"); 

    return 0; 
} 

Answer 1

Indipendentemente da ciò che la formula si usa, avrete bisogno arbitraria precisione aritmetica di ottenere più di 16 cifre. (Poiché "double" ha solo 16 cifre di precisione).

La formula Chudnovsky è la formula più rapida per calcolare Pi e converge a 14 cifre per trimestre. Tuttavia, è estremamente difficile da implementare in modo efficiente.

A causa della complessità di questa formula, non è necessario utilizzare Pi per meno di qualche migliaio di cifre. Quindi non usarlo a meno che tu non sia pronto a fare il tutto esaurito con aritmetica di precisione arbitraria.

una buona implementazione open-source della Formula Chudnovsky utilizzando la libreria GMP è qui: http://gmplib.org/pi-with-gmp.html

Answer 2

Sembra che si sta cercando di calcolare le cifre decimali del π quando la formula BBP viene utilizzato principalmente per il calcolo esadecimale arbitrario cifre di π. In sostanza, la formula BBP può essere usata per calcolare la n^esimo cifra esadecimale di π senza calcolare le cifre precedenti, cifre esadecimali 0, 1, ..., n - 1.

David H. Bailey (Bailey di Bailey-Borwein-Plouffe) ha scritto C and Fortran code per calcolare il n^th cifra esadecimale di π utilizzando la formula BBP. Su una macchina con IEEE 754 doppio aritmetica, la precisione fino a n   ≈   1.18   × contare da 0; cioè π = (3.243F6A8 ...) così l'output del programma quando n = 3 inizia con “F”:

 
position = 3 
fraction = 0.963509103793105 
hex digits = F6A8885A30 

mi piace modificare leggermente la versione C in modo che n (denominato id nel codice) può essere sovrascritta da un argomento della riga di comando:

 
--- piqpr8.c.orig 2011-10-08 14:54:46.840423000 -0400 
+++ piqpr8.c 2011-10-08 15:04:41.524437000 -0400 
@@ -14,14 +14,18 @@ 
/* David H. Bailey  2006-09-08 */ 

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

-main() 
+int main(int argc, char *argv[]) 
{ 
    double pid, s1, s2, s3, s4; 
    double series (int m, int n); 
    void ihex (double x, int m, char c[]); 
    int id = 1000000; 
+ if (argc == 2) { 
+ id = atoi(argv[1]); 
+ } 
#define NHX 16 
    char chx[NHX]; 

@@ -36,6 +40,8 @@ 
    ihex (pid, NHX, chx); 
    printf (" position = %i\n fraction = %.15f \n hex digits = %10.10s\n", 
    id, pid, chx); 
+ 
+ return EXIT_SUCCESS; 
} 

void ihex (double x, int nhx, char chx[]) 

Answer 3

La formula BBP non è adatta per una facile individuazione di n-esima cifra decimale come si ritorna facilmente esadecimale e solo cifre esadecimali. Quindi per ricalcolare in decimali dovrai raccogliere tutte le cifre esadecimali.

È molto meglio usare Newton formula:

Pi/2 = 1 + 1/3 + 1 * 2/3 * 5 + 1 * 2 * 3/3 * 5 * 7 + ... . n!/(2n + 1) !! + ....

sprofonda per schema di Horner:

Pi/2 = 1 + 1/3 * (1 + 2/5 * (1 + 3/7 * (1 + .... .. n/(2n + 1) * (1) .....)))

Così hai scritto Pi come una serie posizionale in cui su ogni posizione frazionaria hai una base diversa usata (n/(2n + 1)), e tutte le cifre sono uguali a 2. Ovviamente converge, poiché quella base è inferiore a 1/2, quindi per calcolare Pi fino a n decimali significativi non è necessario nient'altro che log_2 (10) * n termini (N = 10 * n/3 + 1 è roba perfetta).

Si inizia con la matrice di elementi interi N, tutti uguali 2, e ripetutamente, n volte, effettuare le seguenti operazioni:

1.) moltiplicare tutti gli elementi per 10.

2.) ricalcola ciascun elemento [k] (da N fino a 1) per avere una "cifra" inferiore al denominatore (2 * k + 1),
ma allo stesso tempo è necessario spostare un qoutient alla posizione sinistra, quindi:
q = elemento [k]/(2 * k + 1);
elemento [k]% = (2 * k + 1);
elemento [k-1] + = q * k; // k è il contatore, quindi non aspettare di moltiplicare.

3.) prendere l'elemento [0]. È uguale a 10 * prima cifra, quindi è necessario emettere l'elemento [0]/10 e memorizzare l'elemento
[0]% = 10;

MA c'è un indizio: la somma massima per le cifre massime possibili (2 * n) della formula di Newton è 2. Quindi è possibile ottenere fino a 19/10 dall'elemento [1]. Quando si aggiunge all'elemento [0] (moltiplicato per 10 nel passo 1.) è possibile ottenere 90 + 19 = 109. Così a volte succede cifra digitata sarà [10]. In tal caso, si sa che la cifra corretta è 0 e che 1 deve essere aggiunto alla cifra precedentemente emessa.

Ci sono due modi per risolvere questo problema:

1.) Non uscita viene calcolata l'ultima cifra fino a quando il prossimo. Inoltre, memorizza il numero di nove consecutive e li emette come nove o come 1 seguito da zero in base al primo non 9 cifre.

2.) Inserire le cifre emesse nell'array dei risultati, in modo da poter aggiungere facilmente 1 se [10] si verifica.

Sul mio PC posso calcolare (in Java) 10.000 cifre decimali in 10 s. La complessità è O (n^2).

I valori dell'elemento [k] non superano mai i 12 * k, quindi utilizzando il tipo lungo 64-bit sulla macchina veloce è possibile calcolare più di 10^15 cifre (molto robusto circa).