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:


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

Instead, use fixed-precision arithmetic: each number is represented directly as a number (not a ratio) using a fixed-precision representation.  That will be much more efficient and more accurate.  It won't have the problem of ⌊√a⌋/⌊√b⌋, since you're not limited to a ratio of integers.  You might need to find a library that implements fixed-precision arithmetic, or to implement that part yourself, and you'll probably need to estimate the appropriate precision level based on the number of iterations and the amount of precision achievable with that given precision level.