c++ - 生成平方根的连续分数

我写了这段代码来生成平方根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个!

有人可以帮我吗？

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

最佳答案

根本问题是您不能将非平方的平方根精确地表示为浮点数。

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

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

因此，我们看到，由于0 < ξ - a < 1，所以逼近和实际值之间的差的绝对值会增加。每当出现较大的部分商(ξ - a接近0)时，差异就会增加很大的比例。一旦差的(绝对值)为1或更大，就可以保证下一个计算出的部分商是错误的，但是第一个错误的部分商很可能会更早出现。

Charles mentioned使用n正确数字的原始近似值，您可以计算出连续分数的n部分商。这是一个很好的经验法则，但是正如我们看到的那样，任何大的部分商的成本都更高，从而减少了可获得的部分商的数量，并且有时您会更早地得到错误的部分商。
√139的情况是一个周期相对较长且带有几个较大的部分商的情况，因此第一个错误计算的部分商在周期结束之前出现就不足为奇了(我很惊讶它没有更早出现) 。

使用浮点算法，无法避免这种情况。

但是对于二次波动，我们可以仅使用整数算术来避免该问题。假设您要计算的连续分数扩展

ξ = (√D + P) / Q

其中Q将D - P²除以D > 1并不是一个完美的正方形(如果不满足除数条件，您可以将D替换为D*Q²，将P替换为P*Q，将Q替换为Q²;您的情况为P = 0, Q = 1，这是平凡的满足条件)。将完整的商写为

ξ_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。

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

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

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

因此上述算法的代码变为

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的连续分数。

关于c++ - 生成平方根的连续分数，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/12182701/