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

use arpfloat::{Float, FP256};

///! Calculate the value of PI using the Chudnovsky_algorithm.
///!  cargo run --example calc_pi --release

fn main() {
    // https://en.wikipedia.org/wiki/Chudnovsky_algorithm
    let iterations = 5;

    // Constants:
    let c1 = Float::from_u64(FP256, 10005).sqrt();
    let c2 = Float::from_u64(FP256, 545140134);
    let c3 = Float::from_i64(FP256, -262537412640768000);
    let c16 = Float::from_u64(FP256, 16);
    let c12 = Float::from_u64(FP256, 12);

    // Initial state.
    let mut kc = Float::from_u64(FP256, 6);
    let mut m = Float::from_u64(FP256, 1);
    let mut l = Float::from_u64(FP256, 13591409);
    let mut x = Float::from_u64(FP256, 1);
    let mut s = Float::from_u64(FP256, 13591409);

    for q in 1..iterations + 1 {
        let q3 = Float::from_u64(FP256, q * q * q);
        let k3 = &kc * &(&kc * &kc);
        m = (k3 - (&kc * &c16)) * m / q3;
        l += &c2;
        x *= &c3;
        s += &(&m * &l) / &x;
        kc += &c12;
    }
    let pi = Float::from_u64(FP256, 426880) * (c1 / s);
    println!("pi = {}", pi);
    assert_eq!(pi.as_f64(), std::f64::consts::PI);
}