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

Question

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

Answer 1

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

如果ξ是精确值，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 .

实数的连分数展开ξ是周期性的当且仅当ξ是二次surd，当上述算法中第一对(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;
}

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