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.


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?

  • $\begingroup$ @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. $\endgroup$ – Alex Zhukovskiy Jul 11 '18 at 16:34
  • $\begingroup$ Thanks for all the edits and improvements! Check my edit to see whether it correctly reflects your intent. $\endgroup$ – 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.

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.

  • $\begingroup$ Unfortunately, my language doesn't have FP arithmetics other than integers, so I have to implement it myself. AFAIK all operators should be implemented on top of arbitrary sized integers, but I need to fix all of them to make it work correctly with point that differs from 0. Is there any literature that could help me with design? I would appreciate it. $\endgroup$ – Alex Zhukovskiy Jul 14 '18 at 9:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.