rustedbytes_pi/
lib.rs

1use dashu::integer::IBig;
2use dashu_float::round::mode::Zero;
3use dashu_float::{Context, FBig};
4use dashu_macros::{dbig, ibig};
5
6/*
7https://en.wikipedia.org/wiki/Chudnovsky_algorithm
8
9#Note: For extreme calculations, other code can be used to run on a GPU, which is much faster than this.
10import decimal
11
12def binary_split(a, b):
13    if b == a + 1:
14        Pab = -(6*a - 5)*(2*a - 1)*(6*a - 1)
15        Qab = 10939058860032000 * a**3
16        Rab = Pab * (545140134*a + 13591409)
17    else:
18        m = (a + b) // 2
19        Pam, Qam, Ram = binary_split(a, m)
20        Pmb, Qmb, Rmb = binary_split(m, b)
21
22        Pab = Pam * Pmb
23        Qab = Qam * Qmb
24        Rab = Qmb * Ram + Pam * Rmb
25    return Pab, Qab, Rab
26
27
28def chudnovsky(n):
29    """Chudnovsky algorithm."""
30    P1n, Q1n, R1n = binary_split(1, n)
31    return (426880 * decimal.Decimal(10005).sqrt() * Q1n) / (13591409*Q1n + R1n)
32
33
34print(f"1 = {chudnovsky(2)}")  # 3.141592653589793238462643384
35
36decimal.getcontext().prec = 100 # number of digits of decimal precision
37for n in range(2,10):
38    print(f"{n} = {chudnovsky(n)}")  # 3.14159265358979323846264338...
39 */
40
41/// Compute the binary split for the Chudnovsky algorithm.
42/// This function recursively computes the values of Pab, Qab, and Rab
43/// for the given range [a, b].
44/// The binary split method is used to optimize the computation.
45pub fn binary_split(a: usize, b: usize) -> (IBig, IBig, IBig) {
46    // Base case: if the range contains only one element
47    if b - a == 1 {
48        // Compute Pab, Qab, and Rab for the single element
49        let pab =
50            -(ibig!(6) * a - ibig!(5)) * (ibig!(2) * a - ibig!(1)) * (ibig!(6) * a - ibig!(1));
51        let qab = ibig!(10939058860032000) * a.pow(3);
52        let rab = &pab * (ibig!(545140134) * a + ibig!(13591409));
53
54        // Return the computed values
55        (pab, qab, rab)
56    } else {
57        // Recursive case: split the range into two halves
58        let m = (a + b) / 2;
59
60        // Use rayon for parallel computation to improve performance
61        let ((pam, qam, ram), (pmb, qmb, rmb)) = rayon::join(
62            || binary_split(a, m), // Compute for the left half
63            || binary_split(m, b), // Compute for the right half
64        );
65
66        // Combine the results from the two halves
67        let pab = &pam * &pmb;
68        let qab = &qam * &qmb;
69        let rab = &qmb * &ram + &pam * &rmb;
70
71        // Return the combined results
72        (pab, qab, rab)
73    }
74}
75
76/// Compute the Chudnovsky algorithm for Pi with the specified number of iterations and digits.
77pub fn chudnovsky(iterations: usize, digits: usize) -> FBig<Zero, 10> {
78    // Ensure at least 2 iterations to get a valid result
79    let iterations = iterations.max(2);
80
81    // More precision is needed in order to avoid error propagation
82    let precision = digits + 2;
83
84    let context = Context::<Zero>::new(precision);
85
86    // Use rayon for parallel computation to improve performance
87    let (q, (_p1n, q1n, r1n)) = rayon::join(
88        || context.sqrt(dbig!(10005).repr()).value(),
89        || binary_split(1, iterations),
90    );
91
92    // Perform division in parallel to further optimize
93    let (n, d) = rayon::join(
94        || ibig!(426880) * &q * &q1n,
95        || ibig!(13591409) * &q1n + &r1n,
96    );
97
98    n / d
99}
100
101/// Calculate the number of iterations needed for the Chudnovsky algorithm
102/// to achieve the desired number of digits.
103pub fn chudnovsky_iterations(digits: usize) -> usize {
104    // Each iteration gives approximately 14.181647462 decimal digits
105    const DIGITS_PER_TERM: f64 = 14.181647462;
106    ((digits as f64) / DIGITS_PER_TERM).ceil() as usize
107}
108
109/// Compute the digits of Pi using the Chudnovsky algorithm.
110/// The result is a string representation of Pi with the specified number of digits.
111pub fn compute_pi(digits: usize) -> String {
112    // Compute the number of iterations needed for the Chudnovsky algorithm
113    // to achieve the desired number of digits
114    let iterations = chudnovsky_iterations(digits);
115
116    // Compute pi using the Chudnovsky algorithm
117    let pi = chudnovsky(iterations, digits);
118
119    // Use a more efficient string slicing approach
120    let pi_str = pi.to_string();
121    pi_str.chars().take(digits + 2).collect()
122}
rustedbytes_pi/lib.rs

rustedbytes_pi/
lib.rs
