java - π近似误差导致第29位发散

Question

我有一个我一直在研究的 π 近似程序。它使用 Rananujen-Chudnovsky 级数 ( https://stackoverflow.com/a/531/2058221 ) 来近似 π。我测试了通过将 2 个 π 变体从我测试的不同来源静态加载到内存中而生成的数字。由于某种原因，我计算的经验版本始终无法使超过 29 位的数字正确。我已经将 MathContext 用于几乎所有我可能会用到的地方。注意:broadcastSystemMessage() 本质上是一个 System.out.println() 但它使格式“很好”，本质上将第一个参数视为标题 + : 和其次是正文。

算法实现:

/**
     * The much faster Ramanujan-Chudnovsky algorithm.
     */
    RAMANUJAN_CHUDNOVSKY {

                private final long k1 = 545140134, k2 = 13591409, k3 = 640320,
                k4 = 100100025, k5 = 327843840, k6 = 53360;

                @Override
                public String getAlgorithmName() {
                    return "Ramanujan-Chudnovsky";
                }

                @Override
                public String getAlgorithmFormula() {
                    return "S=Σ from 0 to ∞(-1^n *(6n)! * (k2 + n * k1) / ((n!) ^ 3 * (3n)! * 8 * k4 * k5) ^ n)"
                    + "π = k6 * sqrt(k3) / S";
                }

                @Override
                public BigDecimal initCalculation(long iterations, long significantDigits) {
                    //God, if you're real, please forgive me for this; Java didn't give me a choice.
                    MathContext context = new MathContext((int) significantDigits);
                    BigDecimal s = new BigDecimal(0, context);
                    for (int n = 0; n < iterations; n++) {
                        s = s.add(new BigDecimal(Math.pow(-1, n), context).multiply(new BigDecimal(factorial(6 * n), context), context)
                                .multiply(new BigDecimal(BigInteger.valueOf(k2).add(BigInteger.valueOf(n).multiply(BigInteger.valueOf(k1))),
                                                context), context).divide(new BigDecimal(factorial(n).pow(3), context).multiply(
                                                new BigDecimal(factorial(new BigDecimal(3 * n, context).toBigInteger()), context), context)
                                        .multiply(new BigDecimal(BigInteger.valueOf(8).multiply(BigInteger.valueOf(k4))
                                                        .multiply(BigInteger.valueOf(k5)), context), context).pow(n, context), context), context);
                    }
                    Main.brodcastSystemMessage("Check", k6 + " || + " + k3 + " || " + sqrt(new BigDecimal(k3, context), context));
                    return new BigDecimal(k6, context).multiply(sqrt(new BigDecimal(k3, context), context),
                            context).divide(s, context);
                    //Square Root of k3 approximation: 800.19997500624804755833750301086
                }
            }

π:

public static final BigDecimal π = new BigDecimal("3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196442881097566593344612847564823378678316527120190914564856692346034861045432664821339360726024914127372458700660631558817488152092096282925409171536436789259036001133053054882046652138414695194151160943305727036575959195309218611738193261179310511854807446237996274956735188575272489122793818301194912983367336244065664308602139494639522473719070217986094370277053921717629317675238467481846766940513200056812714526356082778577134275778960917363717872146844090122495343014654958537105079227968925892354201995611212902196086403441815981362977477130996051870721134999999837297804995105973173281609631859502445945534690830264252230825334468503526193118817101000313783875288658753320838142061717766914730359825349042875546873115956286388235378759375195778185778053217122680661300192787661119590921642019893809525720106548586327886593615338182796823030195203530185296899577362259941389124972177528347913151557485724245415069595082953311686172785588907509838175463746493931925506040092770167113900984882401285836160356370766010471018194295559619894676783744944825537977472684710404753464620804668425906949129331367702898915210475216205696602405803815019351125338243003558764024749647326391419927260426992279678235478163600934172164121992458631503028618297455570674983850549458858692699569092721079750930295532116534498720275596023648066549911988183479775356636980742654252786255181841757467289097777279380008164706001614524919217321721477235014144197356854816136115735255213347574184946843852332390739414333454776241686251898356948556209921922218427255025425688767179049460165346680498862723279178608578438382796797668145410095388378636095068006422512520511739298489608412848862694560424196528502221066118630674427862203919494504712371378696095636437191728746776465757396241389086583264599581339");

调用方式:

BigDecimal π = PIAlgorithms.RAMANUJAN_CHUDNOVSKY.initCalculation(10,
            1000);
    brodcastSystemMessage("π Calculation Result", π + "");
    brodcastSystemMessage("", StringUtils.getCommonPrefix(π.toPlainString(), PIAlgorithms.π.toPlainString()).length() + "");

杂项。算法中使用的函数:

/**
 * Custom factorial function.
 *
 * @param integer The integer to use
 * @return The factorial of the number
 */
private static BigInteger factorial(int integer) {
    return factorial(BigInteger.valueOf(integer));
}

/**
 * Custom factorial function.
 *
 * @param integer The integer to use
 * @return The factorial of the number
 */
private static BigInteger factorial(@NotNull BigInteger integer) {
    if (integer.equals(BigInteger.ZERO)) {
        return BigInteger.ONE;
    } else if (integer.compareTo(BigInteger.ZERO) < 0) {
        throw new IllegalArgumentException("Can't take the factorial of a number less than zero");
    }
    BigInteger i = integer.equals(BigInteger.ONE) ? integer : integer.multiply(factorial(integer.subtract(BigInteger.ONE)));
    System.out.println(integer + "! --> " + i);
    return i;
} 
private static BigDecimal sqrt(BigDecimal number, MathContext context) {
    BigDecimal first = new BigDecimal("0", context);
    BigDecimal second = new BigDecimal(Math.sqrt(number.doubleValue()), context);
    while (!first.equals(second)) {
        first = second;
        second = number.divide(first, context);
        second = second.add(first, context);
        second = second.divide(new BigDecimal("2"), context);

    }
    return second;
}

编辑:前 5 项:3.141592653589793238462643383587297242678前 10 项:3.141592653589793238462643383587297242678前 15 项:3.141592653589793238462643383587297242678比较<π>:3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628

Answer 1

我认为您复制的公式不正确。我假设您的代码实现了 getAlgorithmFormula 中显示的公式。该公式与引用文献中显示的公式不同——被加数的分母括号不正确。

在 Maxima 的输出中，g 是您的公式，g2 是正确的公式:

(%i118) g(N);
          N - 1
          ====                                n
          \     (545140134 n + 13591409) (- 1)  (6 n)!
           >    --------------------------------------
          /                         n   3 n       n
          ====    262537412640768000  n!    (3 n)!
          n = 0
(%o118)   --------------------------------------------
                       426880 sqrt(10005)
(%i119) g2(N);
          N - 1
          ====                                n
          \     (545140134 n + 13591409) (- 1)  (6 n)!
           >    --------------------------------------
          /                           n   3
          ====      262537412640768000  n!  (3 n)!
          n = 0
(%o119)   --------------------------------------------
                       426880 sqrt(10005)

您的公式的前两项正确(n = 0 和 n = 1)。这 2 个足以让你正确地给出 28 位数字。从那时起 (n >= 2) 项不正确，因此您只能使用 28 位数字。

有必要用Java吗？处理任意精度数字非常笨拙。我知道Maxima可以很容易地处理，我相信还有其他包。