7

I have been trying to write a simple program to calculate pi using the Chudnovsky algorithm however I keep getting the wrong value output. The latest code i have written is below and outputs:

9.642715619298075837448823278218780086541162343253084414940204168864066834806498471622628399332216456e11

Can anyone tell me where I went wrong.

As Peter de Rivaz pointed out I was discarding the value of b with that fixed the output is now: -1.76779979383639157654764981441635890608880847407921749358841620214761790018058‌​3600120191582474909093e-2 

    Apfloat sum = new Apfloat(0);

    for(int k = 0; k < n; k++) {
        int thrk= 3*k;

        Apfloat a = ApintMath.factorial(6*k); //(6k)! * (-1)^k
        a = a.multiply(ApintMath.pow(new Apint(-1),k));

        Apfloat b = new Apfloat(545140134);
        b = b.multiply(new Apfloat(k));
        b = b.add(new Apfloat(13591409)); // 13591409 + 545140134k

        Apfloat c = ApintMath.factorial(thrk); // (3k!)

        Apfloat d = ApintMath.factorial(k);
        d = ApfloatMath.pow(d, 3); // (k!)^3
        Apfloat e = new Apfloat(640320); 
        e = ApfloatMath.pow(e,(thrk)); // (640320)^(3k)

        a = a.multiply(b); // a is know the numerator
        c = c.multiply(d).multiply(e); // c in know the denominator 

        Apfloat div = a.divide(c.precision(digits));
        sum = sum.add(div);
    }

    Apfloat f = new Apfloat(10005, digits);// Works out the constant sqrt part
    f = ApfloatMath.sqrt(f);
    f = f.divide(new Apfloat(42709344*100));
    Apfloat pi = ApfloatMath.pow(sum.multiply(f), -1);

    System.out.println(pi);
Rudi Kershaw
  • 12,332
  • 7
  • 52
  • 77
Kingkaz
  • 73
  • 5

2 Answers2

5

PROBLEM 1

One problem is in the line:

b.add(new Apfloat(13591409));

this will add 13591409 to b, and discard the result.

Try:

b = b.add(new Apfloat(13591409));

PROBLEM 2

There is a problem in the line:

f = f.divide(new Apfloat(42709344*100));

The problem is that numbers in Java are 32bit integers by default, so 42709344*100 overflows.

Try:

f = f.divide(new Apfloat(42709344*100L));
Peter de Rivaz
  • 33,126
  • 4
  • 46
  • 75
  • Thanks, can believe I missed that know its out putting the same value a previous version was -1.767799793836391576547649814416358906088808474079217493588416202147617900180583600120191582474909093e-2 – Kingkaz Feb 04 '14 at 20:27
2

The denominator for the Chudnovsky algorithm involves 640320^(3k + 3/2) - you only use 640320^(3k).

dcsohl
  • 7,186
  • 1
  • 26
  • 44
  • I have pulled the 3/2 out of the sequence hence the square root part – Kingkaz Feb 04 '14 at 21:25
  • 1
    Ok, I can see that - since 640320^(3k + 3/2) = 640320^(3k) * 640320^(3/2) you can pull out the 640320^(3/2)... but I don't get the 42709344*100 part. 640320^(3/2) = 640320*8*sqrt(10005) = 5122560*sqrt(10005). Shouldn't that be what you pull out? Combine that with the 12 and you have a constant factor of 1/(426880*sqrt(10005))... did I do that right? – dcsohl Feb 04 '14 at 21:33
  • 1
    And now i see that 1/(426880*sqrt(10005)) == sqrt(10005)/(4270934400). So that's ok then. What isn't ok is that 4270934400 is very nearly 2^32. It's way way over the bounds that an `int` can hold, and would in fact probably come out as a negative number. – dcsohl Feb 04 '14 at 21:45