生成平方根的连续分数 [英] Generating continued fractions for square roots

问题描述

我写了这段代码来生成平方根N的连续分数.
但是当N = 139时失败.
输出应为{11,1,3,1,3,7,1,1,2,11,2,1,1,7,3,1,3,1,22}
虽然我的代码给了我一个394个词的序列……其中前几个词是正确的，但是当它达到22个时，它给出了12个！

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!

有人可以帮我吗?

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.

如果ξ是精确值，而x是近似值(假定仍然相当好，因此特别是floor(ξ) = a = floor(x)仍然成立)，则继续分数算法的下一步骤后的差是

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

因此，我们看到，由于0 < ξ - a < 1，逼近值与实际值之差的绝对值增加.每次出现较大的商数(ξ - a接近0)时，差值都会大幅度增加.一旦差的(绝对值)等于1或更大，就可以保证下一个计算出的部分商是错误的，但是第一个错误的部分商很可能会更早出现.

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.

√139的情况是一个周期相对较长的情况，其中包含两个较大的部分商，因此不足以在周期完成之前出现第一个错误计算的部分商也就不足为奇了(我很惊讶它没有不会更早发生.

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

其中Q除以D - P²且D > 1不是一个完美的正方形(如果不满足除数条件，则可以将D替换为D*Q²，将P替换为P*Q和与Q²一起使用；您的情况是P = 0, Q = 1，这是微不足道的).将完整的商写为

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)

并表示部分商a_k.然后

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

，然后使用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],

所以Q_{k+1} = (D - P_{k+1}^2) / Q_k—因为P_{k+1}^2 - P_k^2是Q_k的倍数，所以归纳法Q_{k+1}是整数，而Q_{k+1}除以D - P_{k+1}^2.

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.

当且仅当ξ是二次surd，并且在上述算法中第一对(P_k, Q_k)重复时，周期才结束，实数ξ的连续分数扩展是周期性的.纯平方根的情况特别简单，当k > 0的第一个Q_k = 1和P_k, Q_k始终为非负数时，该周期就结束了.

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.

使用R = floor(√D)，部分商可以计算为

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;
}

可以轻松计算周期长度为182的(例如)√7981的连续分数.

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

这篇关于生成平方根的连续分数的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！