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.