0

I'm trying to get as many correct digits of pi as possible without sacrificing too much speed , so i found a converging series by David Chudnovsky and Gregory Chudnovsky based on Ramanujan's work, when i run it i get up to 3.14159265358973 after 3 the digits are not correct no matter the value of k I also noticed that the counter = 1 , meaning that the loop doesn't work? I would really appreciate the help. My code is:

#include <iostream>
#include <cmath>
#include <iomanip>
using namespace std;


long double  factorial(unsigned long long int x, unsigned long long int n);


int main()
{
    long double big_sum = 0.0;
    unsigned long long int k = 0;
    unsigned long long int n = 0;
    long double numerator;
    long double denominator;
    long double turm;


    turm = (long double)12.0;

    cout << "TYPE NUMBER OF DIGITS : ";
    cin >> n;
    cout << endl;

    /* unsigned long long int counter = 0; */

    for (k = 0; k < pow(10, 10000); k++)
    {
        /* counter++; */
    
        cout << setprecision(n) << fixed;

        numerator = pow((-1.0), k) * factorial(6 * k, n) * (545140134 * k + 13591409);

        denominator = factorial(3 * k, n) * pow(factorial(k, n), 3) * pow(640320, 3 * k + (long 
        double)(3.0 / 2));

        big_sum += turm * (numerator / denominator);

        cout << "PI IS = " << pow(big_sum, -1);
    
        /* cout << "  counter = " << counter; */
    }

    return 0;
}

long double factorial(unsigned long long int x, unsigned long long int n)
{
    if (x == 0) return 1;

    long long int factorial_of_x = 1;

    for (unsigned long long int i = 1; i <= x; x++)
    {
        cout << setprecision(n) << fixed;

        factorial_of_x *= i;
    }

    return factorial_of_x;
}
463035818_is_not_an_ai
  • 109,796
  • 11
  • 89
  • 185
George122V
  • 13
  • 6
  • 1
    what counter = 1? The one in comments? Then please turn the comments into code. Comments are not executed. Strictly speaking there is no `counter` in your code – 463035818_is_not_an_ai Nov 10 '20 at 10:05
  • 1
    i didnt check the details but whenever the factorial appears in a formula you must not keep a term for the factorial alone, this will overflow. Keep a term for `x = (numerator / denominator)` and update that in each iteraton, ie not seperate. The individual terms will surely overflow (no matter what type you use) but the whole `x` will stay somewhere close to `1` – 463035818_is_not_an_ai Nov 10 '20 at 10:08
  • 1
    ...and please do not tag unrelated languages. Even if you do not care about C and C++ being different languages, I do not want to search for a C++ question and then find that the answers are in a different language – 463035818_is_not_an_ai Nov 10 '20 at 10:11
  • Yes the counter in the comments , i know it's not executed i ckecked earlier and its equal to 1 , then i commented to remind me to delete it later , – George122V Nov 10 '20 at 10:12
  • `double` has limited precision. Are you sure you haven't exceeded its bounds? – tadman Nov 10 '20 at 10:12
  • 1
    You are iterating for `pow(10, 10000)` itereations... those are probably too many, and just loop forever. – rodrigo Nov 10 '20 at 10:12
  • idclev 463035818 , thank you i'll try it rn – George122V Nov 10 '20 at 10:12
  • @rodrigo It's not so bad if you have a 10e1000GHz CPU. – tadman Nov 10 '20 at 10:14
  • It's worth noting there's elegant methods for [calculating digits of pi individually](https://cs.uwaterloo.ca/~alopez-o/math-faq/mathtext/node12.html) which side-steps the whole IEEE floating point mess entirely. – tadman Nov 10 '20 at 10:17
  • @rodrigo it doesn't , it looks like it only loops once(counter == 1) since i get the result almost instantly, that's what i'm trying to figure out – George122V Nov 10 '20 at 10:17
  • 1
    Ah, but that is because you are printing from inside the loop. And your `factorial` never ends with `x > 0` because you have `x++` instead of `i++` in the loop increment. – rodrigo Nov 10 '20 at 10:21
  • 1
    `factorial(6 * k, n) ` overflows at k=4. (The factorial grows *fast*.) – molbdnilo Nov 10 '20 at 10:23
  • @rodrigo oh that was it , didn't even notice it , thank you it gets 2 more digits now and then pi is = -nan(ind) , i guess thats the overflow right? – George122V Nov 10 '20 at 10:31
  • Yes... actually with this algorithm it is unclear what will happen first: overflow your long- longs or reaching the precision limit of your long-doubles. – rodrigo Nov 10 '20 at 10:47
  • see [Baking-Pi Challenge - Understanding & Improving](https://stackoverflow.com/a/22295383/2521214) the 2 lines of code at the bottom will give you 800 digits in few [ms] ... but the first [BPP approach](https://stackoverflow.com/a/56035284/2521214) is the way to go... Also as mentioned [factorials grow fast](https://stackoverflow.com/a/18333853/2521214) – Spektre Nov 10 '20 at 11:34

3 Answers3

2

At least these problems

Increment i:

// for (unsigned long long int i = 1; i <= x; x++)
for (unsigned long long int i = 1; i <= x; i++)

Reduce huge iteration limit as factorial(6 * k, n) quickly exceeds unsigned long long range.

//for (k = 0; k < pow(10, 10000); k++) {
for (k = 0; k < 5; k++) {

I reached a stable result after 2 iterations. (Spaces added for clarity)

PI IS = 3.1415926535897 342 091255279861173
PI IS = 3.1415926535897 932 400306920008859
PI IS = 3.1415926535897 932 400306920008859

π Ref = 3.1415926535897 932 384626433832795...

Sample fixed factorial

long double factorial(unsigned x) {
    if (x == 0) return 1;
    long double factorial_of_x = 1.0;
    for (unsigned i = 1; i <= x; i++) {
        factorial_of_x *= i;
    }
    return factorial_of_x;
}
chux - Reinstate Monica
  • 143,097
  • 13
  • 135
  • 256
0

You have those terms in your loop:

numerator = pow((-1.0), k) * factorial(6 * k, n) * (545140134 * k + 13591409);

denominator = factorial(3 * k, n) * pow(factorial(k, n), 3) * pow(640320, 3 * k + (long 
    double)(3.0 / 2));

big_sum += turm * (numerator / denominator);

The problem is that numerator and denominator both grow extremely fast. No matter what type you use, the will overflow at some point.

On the other hand, the term that you add in each iteration turm * (numerator / denominator); will converge towards 1 (because, assuming you made no mistake the whole iteration is set up such that big_sum converges to 1/pi).

Instead of calculating the complete terms in each iteration, only consider the change in the next iteration. I don't know the formula you are using, so only for the sake of illustration I will use extremely simplified expressions (fact(k) is factorial of k):

num = pow((-1.0), k) * fact(k) * (k);
den = fact(k);
big_sum += turm * (numerator / denominator);

Now consider only what changes in each iteration:

num_i+1 = num_i * (-1.0) * (i+1) * (i+1)/(i);
den_i+1 = den_i * (i+1);

This alone won't help of course. Its just a different way to write down the same, the values would be the same. The trick is to realize that while num and den grow fast, dividing them yields a relatively small number. So instead of keeping the seperate terms...

dn_i+1 = num_i+1 / den_i+1
       = num_i/den_i * (-1.0) * (i+1)/(i);

I hope it is possible to grasp the idea. As I said before, I used (wrong) simplified expressions. Though with the real expressions the strategy is the same: Instead of calculating big numbers, accumulate only the factor you apply to the final result. This factor is guaranteed to converge to 1 (because otherwise the iteration would not converge to some finite result).

463035818_is_not_an_ai
  • 109,796
  • 11
  • 89
  • 185
  • Thank you, i think i understand the idea , i'll try it out – George122V Nov 10 '20 at 12:45
  • @George122V fwiw, I don't like the answer myself. Problem was that I realized only while writing that the expressions are more involved than I initially thought. Also I confused adding / multiplying when making the argument that the factor converges to `1`. The point is just that if the iteration converges then the change you apply to `big_sum` will get smaller and smaller with each iteration, while the intermediate results will get bigger and bigger. The whole idea is just to avoid those growing intermediate terms. Hope you can get something from the answer... – 463035818_is_not_an_ai Nov 10 '20 at 12:50
-1

i saw your post you can just use the equation here it is π/4 = 3 tan^-1(1/4) + tan^-1(1/20) + tan^-1(1/1985)

or you can do what Harvard does : a physics sim

here is a 3blue1brown video https://www.youtube.com/watch?v=HEfHFsfGXjs

Pislify
  • 1
  • 1