Generating continued fractions for square roots

Question

I wrote this code for generating Continued Fraction of a square root N.
But it fails when N = 139.
The output should be {11,1,3,1,3,7,1,1,2,11,2,1,1,7,3,1,3,1,22}
Whilst my code gives me a sequence of 394 terms... of which the first few terms are correct but when it reaches 22 it gives 12!

Can somebody help me with this?

vector <int> f;
int B;double A;
A = sqrt(N*1.0);
B = floor(A);
f.push_back(B);                 
while (B != 2 * f[0])) {
    A = 1.0 / (A - B);
    B =floor(A);                            
    f.push_back(B);     
}
f.push_back(B);

Answer 1

The root problem is that you cannot exactly represent the square root of a non-square as a floating-point number.

If ξ is the exact value and x the approximation (which is supposed to be still quite good, so that in particular floor(ξ) = a = floor(x) still holds), then the difference after the next step of the continued fraction algorithm is

ξ' - x' = 1/(ξ - a) - 1/(x - a) = (x - ξ) / ((ξ - a)*(x - a)) ≈ (x - ξ) / (ξ - a)^2

Thus we see that in each step the absolute value of the difference between the approximation and the real value increases, since 0 < ξ - a < 1. Every time a large partial quotient occurs (ξ - a is close to 0), the difference increases by a large factor. Once (the absolute value of) the difference is 1 or greater, the next computed partial quotient is guaranteed to be wrong, but very probably the first wrong partial quotient occurs earlier.

Charles mentioned the approximation that with an original approximation with n correct digits, you can compute about n partial quotients of the continued fraction. That is a good rule of thumb, but as we saw, any large partial quotients cost more precision and thus reduce the number of obtainable partial quotients, and occasionally you get wrong partial quotients much earlier.

The case of √139 is one with a relatively long period with a couple of large partial quotients, so it's not surprising that the first wrongly computed partial quotient appears before the period is completed (I'm rather surprised that it doesn't occur earlier).

Using floating-point arithmetic, there's no way to prevent that.

But for the case of quadratic surds, we can avoid that problem by using integer arithmetic only. Say you want to compute the continued fraction expansion of

ξ = (√D + P) / Q

where Q divides D - P² and D > 1 is not a perfect square (if the divisibility condition is not satisfied, you can replace D with D*Q², P with P*Q and Q with Q²; your case is P = 0, Q = 1, where it is trivially satisfied). Write the complete quotients as

ξ_k = (√D + P_k) / Q_k (with ξ_0 = ξ, P_0 = P, Q_0 = Q)

and denote the partial quotients a_k. Then

ξ_k - a_k = (√D - (a_k*Q_k - P_k)) / Q_k

and, with P_{k+1} = a_k*Q_k - P_k,

ξ_{k+1} = 1/(ξ_k - a_k) = Q_k / (√D - P_{k+1}) = (√D + P_{k+1}) / [(D - P_{k+1}^2) / Q_k],

so Q_{k+1} = (D - P_{k+1}^2) / Q_k — since P_{k+1}^2 - P_k^2 is a multiple of Q_k, by induction Q_{k+1} is an integer and Q_{k+1} divides D - P_{k+1}^2.

The continued fraction expansion of a real number ξ is periodic if and only if ξ is a quadratic surd, and the period is completed when in the above algorithm, the first pair (P_k, Q_k) repeats. The case of pure square roots is particularly simple, the period is completed when first Q_k = 1 for a k > 0, and P_k, Q_k are always nonnegative.

With R = floor(√D), the partial quotients can be calculated as

a_k = floor((R + P_k) / Q_k)

so the code for the above algorithm becomes

std::vector<unsigned long> sqrtCF(unsigned long D) {
    // sqrt(D) may be slightly off for large D.
    // If large D are expected, a correction for R is needed.
    unsigned long R = floor(sqrt(D));
    std::vector<unsigned long> f;
    f.push_back(R);
    if (R*R == D) {
        // Oops, a square
        return f;
    }
    unsigned long a = R, P = 0, Q = 1;
    do {
        P = a*Q - P;
        Q = (D - P*P)/Q;
        a = (R + P)/Q;
        f.push_back(a);
    }while(Q != 1);
    return f;
}

which easily calculates the continued fraction of (e.g.) √7981 with a period length of 182.

Answer 2

The culprit isn't floor. The culprit is the calculation A= 1.0 / (A - B); Digging deeper, the culprit is the IEEE floating point mechanism your computer uses to represent real numbers. Subtraction and addition lose precision. Repeatedly subtracting as your algorithm is doing repeatedly loses precision.

By the time you have calculated the continued fraction terms {11,1,3,1,3,7,1,1,2,11,2}, your IEEE floating point value of A is good to only six places rather than the fifteen or sixteen one would expect. By the time you get to {11,1,3,1,3,7,1,1,2,11,2,1,1,7,3,1,3,1} your value of A is pure garbage. It has lost all precision.

Answer 3

The sqrt function in math is not precise. You can use sympy instead with arbitrarily high precision. Here is a very easy code to calculate continued fractions for any square root or number included in sympy:

from __future__ import division #only needed when working in Python 2.x
import sympy as sp

p=sp.N(sp.sqrt(139), 5000)

n=2000
x=range(n+1)
a=range(n)
x[0]=p

for i in xrange(n):
    a[i] = int(x[i])
    x[i+1]=1/(x[i]-a[i])
    print a[i],

I have set the precision of your number to 5000 and then calculated 2000 continued fraction coefficients in this example code.

Answer 4

In case anyone is trying to solve this is in a language without integers, here is the code from the accepted answer adapted for JavaScript.

Note two ~~ (floor operators) have been added.

export const squareRootContinuedFraction = D =>{
    let R = ~~Math.sqrt(D);
    let f = [];
    f.push(R);
    if (R*R === D) {
        return f;
    }
    let a = R, P = 0, Q = 1;
    do {
        P = a*Q - P;
        Q = ~~((D - P *P)/Q);
        a = ~~((R + P)/Q);
        f.push(a);
    } while (Q != 1);
    return f;
};

Answer 5

I used you algo in a spreadsheet and I get 12 as well, I think you must have made a mistake in your algo, I tried for 253 values, and B did not reach the final value of it.

Can you try to explain a bit more what the algo should do and how it would work ?

I think I got your algo and you made a mistake in your question, it should be 12. For future reference, the algo can be found at this page http://en.wikipedia.org/wiki/Continued_fraction and it is very prone to issues with decimal/numerical computational issues if the inverse value is very close to the integer you are trying to round at.

When doing the prototype under Excel, I could not reproduce the example of the wiki page for 3.245, because at some point Floor() floored the number to 3 instead of 4, so some boundary checking to check for accuracy is required ...

In this case you probably want to add a maximum number of iteration, a tolerance for checking the exit condition (the exit condition should be that A is equal to B btw)