试图在 C++ 中逼近 pi 数字

Question

我试图在不牺牲太多速度的情况下获得尽可能多的正确数字，所以我发现 David Chudnovsky 和 ​​Gregory Chudnovsky 基于 Ramanujan 的工作的收敛系列，当我运行它时，无论 k 的值如何，3 位数字都不正确后，我会达到 3.14159265358973 我还注意到 counter = 1 ，这意味着循环不起作用？我非常感谢您的帮助。 我的代码是：

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

Answer 1

至少这些问题

递增i:

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

减少巨大的迭代限制，因为factorial(6 * k, n) 迅速超过unsigned long long 范围。

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

经过 2 次迭代，我达到了稳定的结果。 （为清楚起见添加了空格）

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

π Ref = 3.1415926535897 932 384626433832795...

---

样本固定阶乘

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

Answer 2

你的循环中有这些术语：

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

问题在于numerator 和denominator 都增长得非常快。无论你使用什么类型，都会在某个时候溢出。

另一方面，您在每次迭代中添加的项 turm * (numerator / denominator); 将收敛到 1（因为，假设您没有弄错，整个迭代的设置使得 big_sum 收敛到 1/pi） .

每次迭代不计算完整项，只考虑下一次迭代的变化。我不知道你使用的公式，所以为了说明，我将使用极其简化的表达式（fact(k) 是k 的阶乘）：

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

现在只考虑每次迭代中的变化：

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

当然，仅此一项是无济于事的。它只是写下相同内容的不同方式，值将相同。诀窍是要意识到虽然num 和den 增长很快，但将它们分开会产生相对较小的数字。所以不要保留单独的条款......

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

我希望有可能掌握这个想法。正如我之前所说，我使用了（错误的）简化表达式。尽管对于实数表达式，策略是相同的：不要计算大数，只累积您应用于最终结果的因子。这个因子保证收敛到1（因为否则迭代不会收敛到某个有限的结果）。

Answer 3

我看到了你的帖子，你可以使用方程式 这里是 π/4 = 3 tan^-1(1/4) + tan^-1(1/20) + tan^-1(1/1985)

或者你可以做哈佛所做的事情：物理模拟

这是一个 3blue1brown 视频 https://www.youtube.com/watch?v=HEfHFsfGXjs