Implement Chudnovsky binary-splitting Pi calculator (replace BBP)

src/main.rs

@@ -1,7 +1,5 @@
 use clap::Parser;
-use rug::Float;
-use rug::ops::Pow;
-use rayon::prelude::*;
+use rug::{Float, Integer, ops::Pow};
 use std::fs::File;
 use std::io::Write;
 use std::path::PathBuf;
@@ -12,7 +10,7 @@ struct Args {
     /// Number of digits of Pi to calculate (digits after the decimal point)
     n: u32,
 
-    /// Number of threads to use
+    /// Number of threads to use (kept for compatibility; Chudnovsky is CPU bound)
     #[arg(short, long, default_value_t = 4)]
     threads: usize,
 
@@ -21,58 +19,61 @@ struct Args {
     output: Option<PathBuf>,
 }
 
-fn bbp_term(k: u32, prec: u32) -> Float {
-    // Compute one BBP term at precision `prec`.
-    let mut term = Float::with_val(prec, 4);
-    term /= Float::with_val(prec, 8 * k + 1);
-
-    let mut term2 = Float::with_val(prec, 2);
-    term2 /= Float::with_val(prec, 8 * k + 4);
-    term -= term2;
-
-    let mut term3 = Float::with_val(prec, 1);
-    term3 /= Float::with_val(prec, 8 * k + 5);
-    term -= term3;
-
-    let mut term4 = Float::with_val(prec, 1);
-    term4 /= Float::with_val(prec, 8 * k + 6);
-    term -= term4;
-
-    let sixteen = Float::with_val(prec, 16);
-    term /= sixteen.pow(k as i32);
-
-    term
+// Binary splitting for the Chudnovsky algorithm.
+// Returns (P, Q, T) as big integers for the interval [a, b)
+fn bs(a: u64, b: u64) -> (Integer, Integer, Integer) {
+    if b - a == 1 {
+        if a == 0 {
+            // P = 1, Q = 1, T = 13591409
+            return (Integer::from(1), Integer::from(1), Integer::from(13591409));
+        }
+        let a_i = Integer::from(a as i128);
+        let p: Integer = (Integer::from(6 * a as i128 - 5)
+            * Integer::from(2 * a as i128 - 1)
+            * Integer::from(6 * a as i128 - 1))
+            .into();
+        let q: Integer = (Integer::from(a as i128).pow(3) * Integer::from(640320i128).pow(3)).into();
+        let mut t: Integer = (p.clone() * Integer::from(13591409i128 + 545140134i128 * a_i)).into();
+        if a % 2 == 1 {
+            t = -t;
+        }
+        return (p, q, t);
+    }
+    let m = (a + b) / 2;
+    let (p1, q1, t1) = bs(a, m);
+    let (p2, q2, t2) = bs(m, b);
+    let p = (&p1 * &p2).into();
+    let q = (&q1 * &q2).into();
+    let t1q2: Integer = (&t1 * &q2).into();
+    let p1t2: Integer = (&p1 * &t2).into();
+    let t = t1q2 + p1t2;
+    (p, q, t)
 }
 
-/// Calculate Pi to `n` decimal digits using a parallelized BBP summation.
-/// Returns a decimal string containing Pi truncated to `n` digits after the decimal point.
-pub fn calculate_pi(n: u32, num_threads: usize) -> Result<String, String> {
+/// Calculate Pi to `n` decimal digits using the Chudnovsky algorithm (binary splitting).
+pub fn calculate_pi_chudnovsky(n: u32) -> Result<String, String> {
     if n == 0 {
         return Err("n must be > 0".into());
     }
-    if num_threads == 0 {
-        return Err("threads must be > 0".into());
-    }
 
-    // Bits of precision: log2(10) ~= 3.321928. Add some guard bits.
+    // Each term of Chudnovsky yields ~14.181647462725477 decimal digits
+    let digits_per_term = 14.181647462725477;
+    let terms = ((n as f64) / digits_per_term).ceil() as u64 + 1;
+
+    // Bits of precision: log2(10) ~= 3.321928. Add guard bits.
     let prec = (n as f64 * 3.3219280948873626).ceil() as u32 + 20;
 
-    // BBP converges in base-16; use a modest overestimate for term count.
-    let num_terms = (n as usize / 1) + 20; // conservative
+    let (_p, q, t) = bs(0, terms);
 
-    // Use rayon thread pool to control threads for parallel work.
-    let pool = rayon::ThreadPoolBuilder::new()
-        .num_threads(num_threads)
-        .build()
-        .map_err(|e| format!("Failed to build thread pool: {}", e))?;
+    // Convert big integers to high-precision floats
+    let prec_u = prec as u32;
+    let qf = Float::with_val(prec_u, q);
+    let tf = Float::with_val(prec_u, t);
 
-    let pi = pool.install(|| {
-        // Parallel iterator over term indices.
-        (0..num_terms as u32)
-            .into_par_iter()
-            .map(|k| bbp_term(k, prec))
-            .reduce(|| Float::with_val(prec, 0), |a, b| a + b)
-    });
+    // C = 426880 * sqrt(10005)
+    let c = Float::with_val(prec_u, 426880) * Float::with_val(prec_u, 10005).sqrt();
+
+    let pi = c * qf / tf;
 
     // Convert to decimal string with a few extra digits for safe truncation.
     let extra = 10usize;
@@ -85,7 +86,6 @@ pub fn calculate_pi(n: u32, num_threads: usize) -> Result<String, String> {
     let out = if pi_string.len() >= end_pos {
         pi_string[..end_pos].to_string()
     } else {
-        // If not enough digits were produced, pad with zeros.
         let mut s = pi_string;
         if !s.contains('.') {
             s.push('.');
@@ -102,7 +102,7 @@ pub fn calculate_pi(n: u32, num_threads: usize) -> Result<String, String> {
 fn main() {
     let args = Args::parse();
 
-    match calculate_pi(args.n, args.threads) {
+    match calculate_pi_chudnovsky(args.n) {
         Ok(pi_str) => {
             if let Some(path) = args.output {
                 match File::create(&path) {
@@ -123,11 +123,11 @@ fn main() {
 
 #[cfg(test)]
 mod tests {
-    use super::calculate_pi;
+    use super::calculate_pi_chudnovsky;
 
     #[test]
     fn pi_10_digits() {
-        let pi = calculate_pi(10, 2).expect("calculation failed");
+        let pi = calculate_pi_chudnovsky(10).expect("calculation failed");
         assert_eq!(pi, "3.1415926535");
     }
 }