Error in π Approximation Causing Divergence on 29th Digit
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日
我认为你复制公式时有错误。我猜测你的代码实现了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可以轻松处理,我相信还有其他软件包也可以处理。