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

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 2

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