Error in π Approximation Causing Divergence on 29th Digit

20 浏览
0 Comments

Error in π Approximation Causing Divergence on 29th Digit

我编写了一个近似π的程序。它使用拉南納森-楚德諾夫斯基級數(https://stackoverflow.com/a/531/2058221)来近似π。我通过在内存中静态加载来自不同来源的两个π的变量进行测试。我对生成的数字进行测试。我使用了 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

admin 更改状态以发布 2023年5月22日
0
0 Comments

我认为你复制公式时有错误。我猜测你的代码实现了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)。这两项足以给你28位正确的结果。从那之后(n>=2),项就不正确了,所以你被限制在了28位。

需要使用Java吗?处理任意精度数字非常笨拙。我知道Maxima可以轻松处理,我相信还有其他软件包也可以处理。

0