CODEWITHSUNDEEP
c++
stack122adam
stack122adam

Reputation: 13

Trying to approximate pi digits in C++

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

Upvotes: 0

Views: 332

Answers (3)

463035818_is_not_an_ai
463035818_is_not_an_ai

Reputation: 123450

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).

Upvotes: 0

chux
chux

Reputation: 154582

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

Upvotes: 2

Pislify
Pislify

Reputation: 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

Upvotes: -1

Related Questions