Bailey–Borwein–Plouffe formula implementation in C++?

Question

EDIT: The requirement was vague and instead of calculating the n-th digit of pi they just wanted pi to the n-th digit not going beyond floats limitation so the brute force way worked for the requirements.

I need to calculate PI the the n-th digit and I wanted to try using the BBP formula but am having difficulties. The equation I typed up doesn't seem to be giving me PI correctly.

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

I was successful with the brute force way of finding PI but that is only so accurate and finding the nth number is difficult.

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

I wanted to find out if anyone had a better idea of how to do this or maybe help with my BBP equation on what I messed up?

Thank you,
LF4

Functional but not very accurate until a few iterations in and then you have to disreguard the last few.

#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;
}

What I'm hoping would be a better program.

#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

The BBP formula is not suitable for easy finding of n-th decimal digit as it easily returns hex and only hex digits. So to recalculate into decimals you will need to collect all hex digits.

It is much better to use Newton formula:

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

It collapses to Horner's scheme:

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

So you have Pi written as a positional series where on each fractional position you have different base used (n/(2n+1)), and all digits are equal to 2. It obviously converges, as that base is less than 1/2, so to calculate Pi up to n significant decimal dgits you need no more than log_2(10)*n terms (N = 10*n/3+1 is perfect stuff).

You start with the array of N integer elements, all equals 2, and repeatedly, n times, do the following:

1.) Multiply all elements by 10.

2.) Recalculate each element[k] (from N down to 1) to have a "digit" less then denominator (2*k+1),
but at the same time you need to move a qoutient to the left position, so:
q = element[k] / (2*k+1);
element[k] %= (2*k+1);
element[k-1] += q * k; //k is the counter, so don't forgrt to multiply.

3.) take element[0]. It equals 10 * first digit, so you need to output element[0] / 10 and store
element[0] %= 10;

BUT there is a clue: the maximal sum for maximal possible digits (2*n) of Newton formula is 2. So you can obtain as much as 19/10 from element[1]. When adding to element[0] (multiplied by 10 in step 1.) you can obtain 90+19=109. So it sometimes happens outputted digit will be [10]. In such a case you know, that the correct digit is 0, and the 1 must be added to the previously outputted digit.

There are two ways of solving this issue:

1.) Not to output the last digit until the next one is calculated. Moreover, store the number of consecutive nines and output them as nines or as 1 followed by zeros depending on first non 9 digit.

2.) Put outputted digits into result array, so you can easily add 1 if [10] occurs.

On my PC I can calculate (in Java) 10,000 decimal digits in 10 s. The complexity is O(n^2).

The values of element[k] never exceeds 12*k, so using 64-bit long type on fast machine you can calculate more than 10^15 digits (very robust approx.).

Answer 2

Regardless of what formula you use, you will need arbitrary precision arithmetic to get more than 16 digits. (Since "double" only has 16 digits of precision).

The Chudnovsky Formula is the fastest known formula for computing Pi and converges at 14 digits per term. However, it is extremely difficult to implement efficiently.

Due to the complexity of this formula, there's no point in using to compute Pi to less than a few thousand digits. So don't use it unless you're ready to go all-out with arbitrary precision arithmetic.

A good open-sourced implementation of the Chudnovsky Formula using the GMP library is here: http://gmplib.org/pi-with-gmp.html

Answer 3

It looks like you are trying to calculate the decimal digits of π when the BBP formula is mainly used to calculate arbitrary hexadecimal digits of π. Basically, the BBP formula can be used to calculate the n^th hexadecimal digit of π without computing the previous digits, hex digits 0, 1, ..., n - 1.

David H. Bailey (the Bailey of Bailey–Borwein–Plouffe) has written C and Fortran code to calculate the n^th hexadecimal digit of π using the BBP formula. On a machine with IEEE 754 double arithmetic, it is accurate up to n ≈ 1.18 × 10⁷ counting from 0; i.e. π = (3.243F6A8...)₁₆ so the output of the program when n = 3 begins with “F”:

 position = 3
 fraction = 0.963509103793105
 hex digits =  F6A8885A30

I like to modify the C version slightly so that n (named id in the code) can be overridden by a command line argument:

--- 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[])