My program for calculating pi using Chudnovsky in C++ precision problem

Question

My code:

#include <iostream>
#include <iomanip>
#include <cmath>

long double fac(long double num) {
    long double result = 1.0;
    for (long double i=2.0; i<num; i++)
       result *= i;
    return result;
}

int main() {
    using namespace std;
    long double pi=0.0;
    for (long double k = 0.0; k < 10.0; k++) {
        pi += (pow(-1.0,k) * fac(6.0 * k) * (13591409.0 + (545140134.0 * k))) 
            / (fac(3.0 * k) * pow(fac(k), 3.0) * pow(640320.0, 3.0 * k + 3.0/2.0));
    }
    pi *= 12.0;
    cout << setprecision(100) << 1.0 / pi << endl;
    return 0;
}

My output:

3.1415926535897637228433865175247774459421634674072265625

The problem with this output is that it outputed 56 digits instead of 100; How do I fix that?

Answer 1

First of all your factorial is wrong the loop should be for (long double i=2.0; i<=num; i++) instead of i<num !!!

As mentioned in the comments double can hold only up to ~16 digits so your 100 digits is not doable by this method. To remedy this there are 2 ways:

1. use high precision datatype

   there are libs for this, or you can implement it on your own you need just few basic operations. Note that to represent 100 digits you need at least

   ceil(100 digits/log10(2)) = 333 bits
   

   of mantisa or fixed point integer while double has only 53

   53*log10(2) = 15.954589770191003346328161420398 digits
   

2. use different method of computation of PI

   For arbitrary precision I recommend to use BPP However if you want just 100 digits you can use simple taylor seriesbased like this on strings (no need for any high precision datatype nor FPU):

    //The following 160 character C program, written by Dik T. Winter at CWI, computes pi  to 800 decimal digits. 
    int a=10000,b=0,c=2800,d=0,e=0,f[2801],g=0;main(){for(;b-c;)f[b++]=a/5;
    for(;d=0,g=c*2;c-=14,printf("%.4d",e+d/a),e=d%a)for(b=c;d+=f[b]*a,f[b]=d%--g,d/=g--,--b;d*=b);}
   

Aside the obvious precision limits Your implementation is really bad from both performance and precision aspects that is why you lost precision way sooner as you hitting double precision limits in very low iterations of k. If you rewrite the iterations so the subresults are as small as can be (in terms of bits of mantisa) and not use too much unnecessary computations here few hints:

1. why are you computing the same factorials again and again

   You have k! in loop where k is incrementing why not just multiply the k to some variable holding actual factorial instead? for example:

   //for (    k=0;k<10;k++){              ... fac(k) ... }
     for (f=1,k=0;k<10;k++){ if (k) f*=k; ... f      ... }
   

2. why are you divide by factorials again and again

   if you think a bit about it then if (a>b) you can compute this instead:

   a! / b! = (1*2*3*4*...*b*...*a) / (1*2*3*4*...*b)
   a! / b! = (b+1)*(b+2)*...*(a)
   

3. I would not use pow at all for this

   pow is "very complex" function causing further precision and performance losses for example pow(-1.0,k) can be done like this:

   //for (     k=0;k<10;k++){       ... pow(-1.0,k) ... }
     for (s=+1,k=0;k<10;k++){ s=-s; ... s           ... }
   

   Also pow(640320.0, 3.0 * k + 3.0/2.0)) can be computed in the same way as factorial, pow(fac(k), 3.0) you can 3 times multipply the variable holding fac(k) instead ...

4. the therm pow(640320.0, 3.0 * k + 3.0/2.0) outgrows even (6k)!

   so you can divide it by it to keep subresults smaller...

These few simple tweaks will enhance the precision a lot as you will overflow the double precision much much latter as the subresults will be much smaller then the naive ones as factorials tend to grow really fast

Putting all together leads to this:

double pi_Chudnovsky() // no pow,fac lower subresult
    {                   // https://en.wikipedia.org/wiki/Chudnovsky_algorithm
    double pi,s,f,f3,k,k3,k6,p,dp,q,r;
    for (pi=0.0,s=1.0,f=f3=1,k=k3=k6=0.0,p=640320.0,dp=p*p*p,p*=sqrt(p),r=13591409.0;k<27.0;k++,s=-s)
        {
        if (k)  // f=k!,  f3=(3k)!, p=pow(640320.0,3k+1.5)*(3k)!/(6k)!, r=13591409.0+(545140134.0*k)
            {
            p*=dp; r+=545140134.0;
            f*=k; k3++; f3*=k3; k6++; p/=k6; p*=k3;
                  k3++; f3*=k3; k6++; p/=k6; p*=k3;
                  k3++; f3*=k3; k6++; p/=k6; p*=k3;
                                k6++; p/=k6;
                                k6++; p/=k6;
                                k6++; p/=k6;
            }
        q=s*r; q/=f; q/=f; q/=f; q/=p; pi+=q;
        }
    return 1.0/(pi*12.0);
    }

as you can see k goes up to 27, while your naive method can go only up to 18 on 64 bit doubles before overflow. However the result is the same as the double mantissa is saturated after 2 iterations ...

Answer 2

I am feeling happy due to following code :)

/*
    I have compiled using cygwin
    change "iostream...using namespace std" OR iostream.h based on your compiler at related OS.
*/
#include <iostream>
#include <iomanip>
#include <cmath>
using namespace std;
long double fac(long double num)
{
    long double result = 1.0;
    for (long double i=2.0; num > i; ++i)
    {
        result *= i;
    }
    return result;
}
int main()
{
    long double pi=0.0;
    for (long double k = 0.0; 10.0 > k; ++k)
    {
        pi += (pow(-1.0,k) * fac(6.0 * k) * (13591409.0 + (545140134.0 * k)))
        / (fac(3.0 * k) * pow(fac(k), 3.0) * pow(640320.0, 3.0 * k + 3.0/2.0));
    }
    pi *= 12.0;
    cout << "BEFORE USING setprecision VALUE OF DEFAULT PRECISION " << cout.precision() << "\n";
    cout << setprecision(100) << 1.0 / pi << endl;
    cout << "AFTER  USING setprecision VALUE OF CURRENT PRECISION   WITHOUT USING fixed " << cout.precision() << "\n";
    cout << fixed;
    cout << "AFTER  USING setprecision VALUE OF CURRENT PRECISION   USING fixed " << cout.precision() << "\n";
    cout << "USING fixed PREVENT THE EARTH'S ROUNDING OFF INSIDE OUR UNIVERSE :)\n";
    cout << setprecision(100) << 1.0 / pi << endl;
    return 0;
}
/*
$ # Sample output:
$ g++  73256565.cpp -o ./a.out;./a.out
$ ./a.out
BEFORE USING setprecision VALUE OF DEFAULT PRECISION 6
3.14159265358976372457810999350158454035408794879913330078125
AFTER   USING setprecision VALUE OF CURRENT PRECISION   WITHOUT USING fixed 100
AFTER   USING setprecision VALUE OF CURRENT PRECISION   USING fixed 100
USING fixed PREVENT THE EARTH'S ROUNDING OFF INSIDE OUR UNIVERSE :)
3.1415926535897637245781099935015845403540879487991333007812500000000000000000000000000000000000000000
*/