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);
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.