Square root calculation using continued fractions to n bits of precision

Question

This is an unsolved problem from my past arbitrary-precision rational numbers C++ assignment.

For calculation, I used this expression from Wikipedia (a being the initial guess, r being its remainder):

I ended up, just by guessing from experiments, with this approach:

1. Use an integer square root function on the numerator/denominator, use that as the guess
2. Iterate the continued fraction until the binary length of the denominator was at least the target precision

This worked well enough to get me through the official tests, however, from my testing, the precision was too high (sometimes almost double) – i.e. the code was inefficient – and I had no proof it worked on any input (and hence no confidence in the code).

A simplified excerpt from the code (natural/rational store arbitrary length numbers, assume all operations return fractions in their simplest form):

rational sqrt(rational input, int precision) {
  rational guess(isqrt(input.numerator), isqrt(input.denominator));  // a
  rational remainder = input - power(guess, 2);                      // r
  rational result = guess;

  rational expansion;
  while (result.denominator.size() <= precision) {
    expansion = remainder / (2 * guess + expansion);
    result = guess + expansion;

    // Handle rational results
    if (power(root, 2) == input) {
      break;
    }
  }
  return result;
}

Can it be done better? If so, how?

Answer 1

Square roots can easily and very accurately be calculated by the General Continued Fractions (GCF). Being general means it can have any positive number as the numerator in contrast to the Regular or Simple Continued Fractions (RCF) where the numerators are all 1s. In order to comprehend the answer as a whole, it is best to start from the beginning.

The method used to solve the square root of any positive number n by a GFC (a + x) whereas a being the integral and x being the continued fractional part, is;

                                                              n − a^2
√n = a + x ⇒ n = a^2 + 2ax + x^2 ⇒ n − a^2 = x(2a + x) ⇒ x = _______
                                                              2a + x

Right at this moment you have a GCF since x nicely gets placed at the denominator and once you replace x with it's definition you get an indefinitely extending definition of x. Regarding a, you are free to choose it among integers which are less than the √n. So if you want to find √11 then a can be chosen among 1, 2 or 3. However it's always better to chose the biggest one in order to be able to simplify the GCF into an RCF at the next stage.

Remember that x = (n − a^2) / (2a + x) and n = 11 and a = 3. Now if we write the first two terms then we may simplify the GCF to RCF with all numerators as 1.

      2         2        divide both            1
x = _____ ⇒ _________ ⇒ numerator and    ⇒ _________ = x
    6 + x   6 +   2      denominator by 2   3 +   1
                _____                           _____
                6 + x                           6 + x

Accordingly our RCF for √11 is;

                        1            ___
√11 = 3 + x ⇒ 3 + _____________ = [3;3,6]
                          1
                  3 + _________
                            1
                      6 + _____
                              1
                          3 + _....
                              6

Notice the coefficient notation [3; 3, 6, 3, 6, ...] which in this particular case resembles an infinite array. This is how RCF's are expressed in coefficient notation, the first item being the a and the tail after ; are the RCF coefficients of x. These two are sufficient since we already know that in RCF all numerators are fixed to 1.

Coming back to your precision question. You now have √11 = 3 + x where x is your RCF as [3;3,6,3,6,3,6...]. Normally you can try by picking a depth and reducing from right like [3,3,6,3,6,3,6...].reduceRight((p,c) => c + 1/p) as it would be done in JS. Not a precise enough result.? Then try it again from another depth. This is in fact how it is descriped in the linked wikipedia topic as bottom up. However it would be much efficient to go from left to right (top to bottom) by calculating the intermediate convergents one after the other, at a single pass. Every next intermediate convergent yields a better precision for you to test and decide weather to stop or continue. When you reach to a coefficient sufficient enough just stop there. Having said that, once you reach to the desired coefficient you may still do some fine tuning by increasing or decreasing that coefficient. Decreasing the coefficients at even indices or increasing the ones at odd indices would decrease the convergent and vice versa.

So in order to be able to do a left to right (top to bottom) analysis there is a special rule as

n2/d2 = (xn * n1 + n0)/(xn * d1 + d0)

1. We need to know last two interim convergents (n0/d0 and n1/d1) along with the current coefficient xn in order to be able calculate the next convergent (n2/d2).
2. We will start with two initial convergents as Infinity (n0/d0 = 1/0) and the a that we've chosen above (Remember √n = a + x) which is 3 so (n1/d1 = 3/1). Knowing that the 3 before the semicolon is in fact a, our first xn is the 3 right after the semicolon in our coefficients array [3;»» 3 ««,6,3,6,3,6...].
3. After we calculate n2/d2 and do our test, if need be, for the next step we will shift our convergents to the left so that we have the last two ready to calculate the next convergent. n0/d0 <- n1/d1 <- n2/d2

Here i present the table for the n2/d2 = (xn * n1 + n0)/(xn * d1 + d0) rule.

n0/d0  n1/d1   xn  index   n2/d2    decimal val.
_____  ______  __  _____  ________  ____________
 1/0     3/1    3  1 odd    10/3    3.33333333..
 3/1    10/3    6  2 evn    63/19   3.31578947..
10/3    63/19   3  3 odd   199/60   3.31666666..
63/19  199/60   6  4 evn  1257/379  3.31662269..
  .       .     .    .        .           .
  .       .     .    .        .           .

So as you may notice we are very quickly approaching to √11 which is 3.31662479... Note that the odd indices overshoot and evens undershoot due to cascading reciprocals. Since √11 is an irrational this will continue convergining indefinitely up until we say enough.

Remember, as mentioned earlier, once you reach to the desired coefficient you may still do some fine tuning by increasing or decreasing that coefficient (xn). Decreasing the coefficients at even indices or increasing the ones at odd indices would decrease the convergent and vice versa.

The problem here is, not all √n can simply be turned into RCF by a simple division as shown above. For a more generalized way to generate RCF from any √n you may check a more recent answer of mine.