From beac62667923c4afd5d3d885a7f46aeed8ea7071 Mon Sep 17 00:00:00 2001
From: Hatter Jiang <jht5945@gmail.com>
Date: Tue, 22 Jul 2025 00:06:40 +0800
Subject: [PATCH] feat: add chudnovsky algorithm

---
 README.md                                     |  1 +
 .../me/hatter/math/ChudnovskyAlgorithm.java   | 35 +++++++++++++++++++
 .../java/me/hatter/math/util/MathUtil.java    | 24 +++++++++++++
 3 files changed, 60 insertions(+)
 create mode 100644 src/main/java/me/hatter/math/ChudnovskyAlgorithm.java

diff --git a/README.md b/README.md
index 83b4bb6..9a67cc9 100644
--- a/README.md
+++ b/README.md
@@ -10,6 +10,7 @@
 | LeibnizFormula          | $\pi$ 的莱布尼茨公式      | [leibniz-formula](https://hatter.ink/static/resource/muboard/?id=leibniz-formula&version=1)                                 |
 | EulerNumber             | 欧拉数                | [eulers-number](https://hatter.ink/static/resource/muboard/?id=eulers-number&version=latest&full=1)                         |
 | BirthdayParadox         | 生日悖论               | [birthday-paradox](https://hatter.ink/static/resource/muboard/?id=birthday-paradox&version=latest&full=1)                   |
+| ChudnovskyAlgorithm     | 楚德诺夫斯基算法           | [chudnovsky-algorithm](https://hatter.ink/static/resource/muboard/?id=chudnovsky-algorithm&version=latest&full=1)           |
 
 
 <br>
diff --git a/src/main/java/me/hatter/math/ChudnovskyAlgorithm.java b/src/main/java/me/hatter/math/ChudnovskyAlgorithm.java
new file mode 100644
index 0000000..e38d861
--- /dev/null
+++ b/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);
+    }
+}
diff --git a/src/main/java/me/hatter/math/util/MathUtil.java b/src/main/java/me/hatter/math/util/MathUtil.java
index 62ec09e..33f61ea 100644
--- a/src/main/java/me/hatter/math/util/MathUtil.java
+++ b/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;
+    }
 }