1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
use rug::{Float, ops::Pow};
use std::convert::TryFrom;

/// Calculates the value of pi to a specified number of decimal places using
/// the Gauss-Legendre algorithm.
///
/// # Arguments
///
/// * `digits` - The number of decimal places of pi to calculate.
///
/// # Returns
///
/// A `Float` representing the calculated value of pi to the specified
/// number of decimal places.
///
pub fn compute_pi(digits: usize) -> Float {
    let precision = ((digits as f64) * 3.3219280948874).ceil() as u32 + 10;
    let threshold = Float::with_val(precision, 10).pow(-i32::try_from(digits).unwrap());
    let mut a = Float::with_val(precision, 1);
    let two = Float::with_val(precision, 2);
    let mut b = Float::with_val(precision, 1.0 / two.sqrt());
    let mut t = Float::with_val(precision, 0.25);
    let mut p = Float::with_val(precision, 1);
    let mut pi_old = Float::with_val(precision, 0);

    let pi_result = loop {
        let sum = a.clone() + b.clone();
        let a_next = Float::with_val(precision, &sum / 2.0);
        let product = Float::with_val(precision, &a * &b);
        b = product.sqrt();
        let difference = Float::with_val(precision, &a - &a_next);
        let difference_squared = difference.pow(2);
        t -= &p * difference_squared;
        a = a_next;
        p *= 2;
        let denominator = Float::with_val(precision, &t * 4.0);
        let numerator = Float::with_val(precision, (&sum).pow(2));
        let pi = numerator / denominator;
        let pi_diff = Float::with_val(pi.prec(), &pi - &pi_old).abs();
        if pi_diff < threshold {
            break pi;
        }
        pi_old = pi;
    };
    pi_result
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn test_compute_pi_10_digits() {
        let pi_computed = compute_pi(10);
        let pi_expected = Float::with_val(256, 3.1415926535);
        assert!((pi_computed - &pi_expected).abs() < Float::with_val(256, 1e-10));
    }
}