# Implementing Gauss–Legendre algorithm using arbitrary-length rationals

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:

1. 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.
2. 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?

• @D.W. ok, I always was thinking that programming languages are proper subset of pseudocode languages, but I removed all language-specific details. The main idea resides: for some reason rationals (even with very big numerators and denominators) gives very bad results. – Alex Zhukovskiy Jul 11 '18 at 16:34
• Thanks for all the edits and improvements! Check my edit to see whether it correctly reflects your intent. – D.W. Jul 11 '18 at 22:16

It sounds like you are using arbitrary-precision arithmetic, i.e., you represent each rational number as a ratio $a/b$ where $a,b$ are stored as arbitrary size integers (bignums). If that's the case, that's probably the source of your problem.
Don't do that. That's highly inefficient, because you'll end up with very large integers $a,b$, that may be far larger than needed for the desired precision level; and it'll be imprecise, because of the ⌊√a⌋/⌊√b⌋ problem you mention.