diff --git a/src/main.rs b/src/main.rs index 6aea6bd..87a606a 100644 --- a/src/main.rs +++ b/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, } -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 { +/// Calculate Pi to `n` decimal digits using the Chudnovsky algorithm (binary splitting). +pub fn calculate_pi_chudnovsky(n: u32) -> Result { 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 { 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 { 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"); } }