使用Chudnovsky算法计算pi时出现错误 - Java

Question

使用Chudnovsky算法计算pi时出现错误 - Java

7

我一直在尝试使用Chudnovsky算法编写一个简单的程序来计算pi，但是我始终得到错误的输出值。我最近编写的代码如下，并输出：

9.642715619298075837448823278218780086541162343253084414940204168864066834806498471622628399332216456e11

有人能告诉我哪里出了问题吗？

正如Peter de Rivaz指出的那样，我丢弃了b的值，修复后输出为：-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);

- Kingkaz

2个回答

2

Chudnovsky算法中的分母涉及到640320^(3k + 3/2)，但实际上只需要使用640320^(3k)。

- dcsohl

我已经从序列中提取了3/2，因此平方根部分。 - Kingkaz

1

好的，我可以看到 - 因为 640320^(3k + 3/2) = 640320^(3k) * 640320^(3/2)，所以你可以提取出 640320^(3/2)... 但是我不明白 42709344100 部分。640320^(3/2) = 6403208sqrt(10005) = 5122560sqrt(10005)。难道这不是你要提取的吗？将其与12相结合，你就有了一个常数因子 1/(426880*sqrt(10005))... 我做对了吗？ - dcsohl

1

现在我知道 1/(426880*sqrt(10005)) == sqrt(10005)/(4270934400)。所以那没问题。但不好的是，4270934400 几乎等于 2^32。它远远超出了 int 可以容纳的范围，并且实际上可能会变成负数。 - dcsohl

没错！我没有意识到，谢谢你，现在它可以工作了 :) - Kingkaz

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Peter de Rivaz · Accepted Answer

问题1

一个问题出现在以下这行代码中：

b.add(new Apfloat(13591409));

这将把 13591409 添加到 b 中，并丢弃结果。

尝试：

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

问题 2

这条线路存在问题：

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

问题在于Java中默认使用32位整数，因此42709344*100会溢出。

尝试如下：

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