feat: add chudnovsky algorithm

src/main/java/me/hatter/math/ChudnovskyAlgorithm.java

@@ -0,0 +1,35 @@
+package me.hatter.math;
+
+import me.hatter.math.util.MathUtil;
+
+import java.math.BigDecimal;
+import java.math.BigInteger;
+import java.math.RoundingMode;
+
+// https://zh.wikipedia.org/zh-cn/楚德诺夫斯基算法
+public class ChudnovskyAlgorithm {
+
+    public static void main(String[] args) {
+        // 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679
+        System.out.println(computePi());
+    }
+
+    public static BigDecimal computePi() {
+        BigDecimal c = BigDecimal.valueOf(426880).multiply(MathUtil.sqrt(BigDecimal.valueOf(10005), 100));
+        BigInteger m = BigInteger.valueOf(1);
+        BigInteger l = BigInteger.valueOf(13591409);
+        BigInteger x = BigInteger.valueOf(1);
+        BigInteger k = BigInteger.valueOf(6);
+        BigDecimal s = BigDecimal.valueOf(l.longValue());
+
+        for (int i = 0; i < 10; i++) {
+            m = (k.pow(3).subtract(k.multiply(BigInteger.valueOf(16)))).multiply(m).divide(BigInteger.valueOf(i + 1).pow(3));
+            l = l.add(BigInteger.valueOf(545140134));
+            x = x.multiply(BigInteger.valueOf(-262537412640768000L));
+            s = s.add(new BigDecimal(m.multiply(l)).divide(new BigDecimal(x), 100, RoundingMode.FLOOR));
+            k = k.add(BigInteger.valueOf(12));
+        }
+
+        return c.divide(s, 100, RoundingMode.FLOOR);
+    }
+}

src/main/java/me/hatter/math/util/MathUtil.java

@@ -41,4 +41,28 @@ public class MathUtil {
         }
         return v.add(BigDecimal.valueOf(nums[0]));
     }
+
+    // Newton-Raphson法によるBigDecimalの平方根計算
+    public static BigDecimal sqrt(BigDecimal n, int scale) {
+        if (n.compareTo(BigDecimal.ZERO) < 0) {
+            throw new IllegalArgumentException("Cannot calculate square root of a negative number");
+        }
+
+        if (n.compareTo(BigDecimal.ZERO) == 0) {
+            return BigDecimal.ZERO;
+        }
+
+        BigDecimal x = n.divide(new BigDecimal("2"), scale, RoundingMode.HALF_UP);
+        BigDecimal precision = new BigDecimal("1").movePointLeft(scale); // 誤差の許容範囲
+        BigDecimal lastX;
+
+        do {
+            lastX = x;
+            x = n.divide(x, scale, RoundingMode.HALF_UP)
+                    .add(x)
+                    .divide(new BigDecimal("2"), scale, RoundingMode.HALF_UP);
+        } while (x.subtract(lastX).abs().compareTo(precision) > 0);
+
+        return x;
+    }
 }