Loe*_*rio 9 c++ math square-root
我写了这个代码用于生成平方根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);
Dan*_*her 18
根本问题是您无法将非正方形的平方根精确表示为浮点数.
如果?是精确值和x近似值(假设仍然非常好,特别是floor(?) = a = floor(x)仍然保持),则继续分数算法的下一步之后的差异是
?' - x' = 1/(? - a) - 1/(x - a) = (x - ?) / ((? - a)*(x - a)) ? (x - ?) / (? - a)^2
因此,我们看到在每一步中,近似值与实际值之间的差值的绝对值增加,因为0 < ? - a < 1.每当出现大的部分商(? - a接近0)时,差异就会增大.一旦(绝对值)差值为1或更大,下一个计算出的部分商保证是错误的,但很可能先发生了第一个错误的部分商.
查尔斯 提到了近似值,通过使用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为a 时完成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 的连续分数.
罪魁祸首不是floor。罪魁祸首是计算A= 1.0 / (A - B);深入挖掘，罪魁祸首是计算机用来表示实数的 IEEE 浮点机制。减法和加法会失去精度。当你的算法重复执行时，重复减法会失去精度。
当您计算出连分数项 {11,1,3,1,3,7,1,1,2,11,2} 时，A 的 IEEE 浮点值仅适用于 6 位，而不是人们预计十五或十六。当您到达 {11,1,3,1,3,7,1,1,2,11,2,1,1,7,3,1,3,1} 时，您的 A 值纯粹是垃圾。它已经失去了所有的精确性。