I am trying to re-implement SuperPI myself in Rust, but the results I get are not very accurate. SuperPI computes pi using the Gauss-Legendre algorithm.
The Gauss–Legendre algorithm is quite simple, but the problem is how to store values with full precision. I store all numbers as arbitrary sized rationals, but I am still getting weird results.
Pseudocode:
a0 = 1/1;
b0 = 470832/665857; // an approximation to 1/√2
t0 = 1/4;
p0 = 1;
for i in 0..n {
a_i+1 = (a + b)/2;
squared_b_i+1 = a * b;
b_i+1 = sqrt(squared_b_i+1); // calculated as sqrt(numerator)/sqrt(denominator)
t_i+1 = t - p * (a_i - a_i+1)^2;
p_i+1 = 2 * p;
}
let result = (a_n+b_n)^2 / (4t_n);
(actual code may be found here)
My implementation has at least two problems:
- It doesn't matter how good
1/√2
approximation I take (470832/665857 in this example), the algorithm doesn't converge to the correct value. - Calculating √(a/b) as √a/√b is not very accurate, because it does ⌊√a⌋/⌊√b⌋ which can introduce a large error.
As a result, my implementation calculates 11 true digits per 20 iterations, while SuperPI does 1M digits.
How do I implement the Gauss–Legendre algorithm efficiently, given the above challenges?