dashu_float/math/
consts.rs

1use crate::{
2    error::assert_limited_precision,
3    fbig::FBig,
4    repr::{Context, Word},
5    round::{Round, Rounded},
6};
7use dashu_base::{BitTest, Sign, UnsignedAbs};
8use dashu_int::{IBig, UBig};
9
10impl<R: Round> Context<R> {
11    /// Calculate π using the Chudnovsky algorithm with binary splitting.
12    ///
13    /// The Chudnovsky algorithm is one of the most efficient methods for
14    /// high-precision π calculation, providing ~14.18 decimal digits per term.
15    ///
16    /// # Methodology
17    /// We use Binary Splitting to evaluate the series. This technique transforms
18    /// the linear-time summation into a recursive tree evaluation. By combining
19    /// terms into large products, it allows the library to leverage fast
20    /// multiplication algorithms (like Toom-3 or FFT) as the numbers grow,
21    /// leading to significant performance gains over simple iterative summation.
22    ///
23    /// // TODO: consider adding a static cache for π at common precisions.
24    #[must_use]
25    pub fn pi<const B: Word>(&self) -> Rounded<FBig<R, B>> {
26        assert_limited_precision(self.precision);
27
28        // Calculate required bits based on target precision in base B.
29        // bits = ceil(precision * log2(B))
30        let bits = if B.is_power_of_two() {
31            self.precision.saturating_mul(B.ilog2() as usize)
32        } else {
33            self.precision.saturating_mul(B.ilog2() as usize + 1)
34        };
35
36        let num_terms = (bits * 100 / 4708) + 1;
37        let guard_bits = num_terms.bit_len() + 32;
38        let work_bits = bits + guard_bits;
39
40        // Evaluate the series components using binary splitting
41        let (_p, q, t) = chudnovsky_bs(0, num_terms);
42
43        // Final formula: pi = (426880 * sqrt(10005) * Q) / T
44
45        // Convert work bits back to base B precision.
46        // precision_B = ceil(work_bits / log2(B))
47        let work_precision = if B == 2 {
48            work_bits
49        } else {
50            work_bits / B.ilog2() as usize + 1
51        };
52        let work_context = Self::new(work_precision);
53
54        let q_f = work_context.convert_int::<B>(q.into()).value();
55        let t_f = work_context.convert_int::<B>(t).value();
56
57        let sqrt_10005 = work_context
58            .sqrt(&work_context.convert_int::<B>(10005.into()).value().repr)
59            .value();
60        let constant = work_context.convert_int::<B>(426_880.into()).value();
61
62        let pi = (constant * sqrt_10005 * q_f) / t_f;
63        pi.with_precision(self.precision)
64    }
65}
66
67/// Binary splitting implementation for the Chudnovsky series.
68/// Returns (P, Q, T) for the range [a, b).
69fn chudnovsky_bs(a: usize, b: usize) -> (UBig, UBig, IBig) {
70    if b - a == 1 {
71        // Base case: calculate single term
72        if a == 0 {
73            return (UBig::ONE, UBig::ONE, IBig::from_parts_const(Sign::Positive, 13_591_409));
74        }
75
76        let k = a as u64;
77        let p = UBig::from(6 * k - 5) * (2 * k - 1) * (6 * k - 1);
78        let q = UBig::from(k).pow(3) * UBig::from_u64(10_939_058_860_032_000);
79        let t_val = IBig::from_parts_const(Sign::Positive, 13_591_409)
80            + IBig::from_parts_const(Sign::Positive, 545_140_134) * k;
81        let t_abs = &p * t_val.unsigned_abs();
82        let t = IBig::from(t_abs) * Sign::from(a % 2 == 1);
83        return (p, q, t);
84    }
85
86    // Recursive step
87    let mid = (a + b) / 2;
88    let (p_l, q_l, t_l) = chudnovsky_bs(a, mid);
89    let (p_r, q_r, t_r) = chudnovsky_bs(mid, b);
90
91    let p = &p_l * &p_r;
92    let q = &q_l * &q_r;
93    // T = T_L * Q_R + T_R * P_L
94    let t = IBig::from(q_r) * t_l + IBig::from(p_l) * t_r;
95    (p, q, t)
96}
97
98impl<R: Round, const B: Word> FBig<R, B> {
99    /// Calculate π with the given precision and the default rounding mode.
100    #[inline]
101    #[must_use]
102    pub fn pi(precision: usize) -> Self {
103        Context::<R>::new(precision).pi().value()
104    }
105}
dashu_float/math/consts.rs

dashu_float/math/
consts.rs
