用Chudnovsky算法计算pi的java错误结果

Question

我试图用Chudnovsky algorithm来计算π。我做了一些研究，使用了大小数。但当我运行代码时，它给了我错误的值。我检查了我的代码很多次，但没有发现任何改进

public class PI {
    private static void Chudnovsky(int numSigDigits) {
        long time1 = System.nanoTime();
        int numDigits = numSigDigits + 10;
        MathContext mc = new MathContext(numDigits,RoundingMode.HALF_EVEN), mcx=new MathContext(numDigits*20, RoundingMode.HALF_EVEN);

        BigDecimal ratC = new BigDecimal(426880L, mc),
                irratC = new BigDecimal(10005L, mcx),
                incrL = new BigDecimal(545140134L, mc),
                l = new BigDecimal(13591409L, mcx),
                mulX = new BigDecimal(-262537412640768000L, mc),
                x = new BigDecimal(1L, mcx),
                incrK = new BigDecimal(12L, mc),
                k = new BigDecimal(6L, mc),
                m = new BigDecimal(1L, mc),
                b16 = new BigDecimal(16L, mc), j;

        BigDecimal C = ratC.add(irratC.sqrt(mcx), mcx);
        BigDecimal sum = new BigDecimal(0, mcx);

        for (int i = 0; i < numDigits;i++) {
            j = new BigDecimal(i + 1, mc);

            sum = sum.add(m.multiply(l, mc).divide(x, mc), mc);

            l = l.add(incrL, mcx);
            x = x.multiply(mulX, mcx);
            m = m.multiply(k.pow(3, mc).subtract(b16.multiply(k, mc), mc).divide(j.pow(3, mc), mc), mc);
            k = k.add(incrK, mc);
        }
        
        long time2 = System.nanoTime();
        System.out.println((double) (time2 - time1) / 1000000000);

        System.out.println(C.divide(sum, mcx).setScale(numSigDigits, RoundingMode.HALF_EVEN).toEngineeringString());
    }

    public static void main(String[] args) {
        Chudnovsky(15);
        System.out.println(Math.PI);
        //System.out.println(str);
    }
}

输出：

0.270327757  //time
0.031415434926348  //my value of pi
3.141592653589793  //Math.PI

我从https://en.wikipedia.org/wiki/Chudnovsky_algorithm中获取了这个算法

Answer 1

你应该改变

BigDecimal C = ratC.add(irratC.sqrt(mcx), mcx);

到

BigDecimal C = ratC.multiply(irratC.sqrt(mcx), mcx);