我编写了这段代码来生成平方根 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 或更大，则保证下一个计算的部分商是错误的，但很可能第一个错误的部分商出现得更早。
查尔斯 提到了一个近似值，即使用具有n正确数字的原始近似值，您可以计算出n连分数的部分商。这是一个很好的经验法则，但是正如我们所看到的，任何大的部分商都会花费更多的精度，从而减少可获得的部分商的数量，并且有时你会更早地得到错误的部分商。
的情况√139是一个相对较长的周期和几个大的部分商，所以在周期完成之前出现第一个错误计算的部分商也就不足为奇了（我很惊讶它没有更早发生）。
使用浮点运算，没有办法阻止这种情况。
但是对于二次曲线的情况，我们可以通过仅使用整数算术来避免这个问题。假设您要计算的连分数展开式
ξ = (√D + P) / Q
其中Q除D - P²且D > 1不是完全平方（如果不满足可除性条件，您可以D用D*Q²、PwithP*Q和Qwith替换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 的 (eg) 的连分数。
罪魁祸首不是floor。罪魁祸首是计算A= 1.0 / (A - B);深入挖掘，罪魁祸首是你的计算机用来表示实数的 IEEE 浮点机制。减法和加法失去精度。当你的算法反复做减法时，反复减法会失去精度。
当您计算出连分数项 {11,1,3,1,3,7,1,1,2,11,2} 时，您的 A 的 IEEE 浮点值仅适用于六位而不是十五或十六一个人会期望。当您到达 {11,1,3,1,3,7,1,1,2,11,2,1,1,7,3,1,3,1} 时，您的 A 值是纯垃圾. 它已经失去了所有的精确度。
数学中的 sqrt 函数并不精确。您可以使用任意高精度的 sympy 代替。这是一个非常简单的代码，可以计算 sympy 中包含的任何平方根或数字的连分数：
from __future__ import division #only needed when working in Python 2.x
import sympy as sp
p=sp.N(sp.sqrt(139), 5000)
n=2000
x=range(n+1)
a=range(n)
x[0]=p
for i in xrange(n):
a[i] = int(x[i])
x[i+1]=1/(x[i]-a[i])
print a[i],
我已将您的数字的精度设置为 5000，然后在此示例代码中计算了 2000 个连分数系数。
如果有人试图用一种没有整数的语言来解决这个问题，这里是来自接受的答案的代码，适用于JavaScript.
注意~~添加了两个（楼层操作员）。
export const squareRootContinuedFraction = D =>{
let R = ~~Math.sqrt(D);
let f = [];
f.push(R);
if (R*R === D) {
return f;
}
let a = R, P = 0, Q = 1;
do {
P = a*Q - P;
Q = ~~((D - P *P)/Q);
a = ~~((R + P)/Q);
f.push(a);
} while (Q != 1);
return f;
};
我在电子表格中使用了你的算法，我也得到了 12，我认为你的算法一定犯了错误，我尝试了 253 个值，但 B 没有达到它的最终值。
你能试着解释一下算法应该做什么以及它是如何工作的吗？
我想我得到了你的算法，而你在你的问题中犯了一个错误，应该是 12。为了将来参考，算法可以在这个页面http://en.wikipedia.org/wiki/Continued_fraction找到，它很容易如果反数值非常接近您要四舍五入的整数，则会出现十进制/数值计算问题。
在 Excel 下做原型时，我无法重现 3.245 的 wiki 页面示例，因为在某些时候 Floor() 将数字设为 3 而不是 4，因此需要进行一些边界检查以检查准确性......
在这种情况下，您可能想要添加最大迭代次数，检查退出条件的容差（退出条件应该是 A 等于 B btw）
您的代码没有计算n. 它尝试计算已计算的连分数√n。我的意思是没关系，但是，如果它是正确的，你的方法更适合一般十进制到理性的转换。然而，对于常规（简单）连分数（所有分子都是 1）的 sqrt 函数，算法略有不同。
然而问题还没有结束。是的，作为一项规则，CF 系数√n是重复回文的形式，它以第一个非零系数的双倍结尾。比如√31 =[5;1,1,3,5,3,1,1,10,1,1,3,5,3,1,1,10..]。现在没有简单的方法来查询每个给定的回文长度n。有一些已知的模式，但它们远未定义所有的通用模式n。所以在第一个回文结束时停止迭代是一种非常不确定的方法。想象
__
√226 =[15;30]
尽管
____________________________________________________
√244 =[15;1,1,1,1,1,2,1,5,1,1,9,1,6,1,9,1,1,5,1,2,1,1,1,1,1,30]
如果您决定在2*f[0]大多数情况下停止迭代，您会得到一个像 in 这样的错误近似值，或者像 in case√226这样一个过度计算的近似值。√244此外，一旦n成长，追逐回文的尽头变得更加没有意义，因为你永远不需要这样的精确度。
___________________________________________________________________________________________________________________________________________________________________________
√7114 = [84;2,1,9,3,1,10,2,23,1,1,1,1,1,2,1,27,2,1,1,3,1,2,1,1,1,16,4,3,1,3,2,1,6,18,1,1,2,6,11,11,6,2,1,1,18,6,1,2,3,1,3,4,16,1,1,1,2,1,3,1,1,2,27,1,2,1,1,1,1,1,23,2,10,1,3,9,1,2,168]
在这种情况下，一旦获得必要的精度就停止迭代是合理的。正如我在开头提到的，有两种方法。
√n它与 1 基本相同，但针对平方根进行了调整。在这种情况下，您不提供√nbut n。思路如下。我们必须为在分子处包含平方根值的输入定义一个通用形式。然后我们尝试在连分数部分达到相同的表达式以便能够迭代。
让我们的输入是形式
q + √n
______
p
对于简单的平方根运算，我们可以假设qis0和pis 1。如果我们可以在下一阶段建立这种形式，那么我们可以轻松地进行迭代。
从初始阶段开始，其中,q = 0是的整数部分和是浮点部分，我们的目标是带入形式；p = 1m√n1/xx(q + √n) / p
1 1 1 (√n + m) √n + m
√n = m + ___ ⇒ x = _______ ⇒ x = ________ . ________ ⇒ x = ________
x √n - m (√n - m) (√n + m) n - m^2
现在√n在分子处，我们有以下形式；
√n + q
x = ______
p
哪里q = m和p = n - m^2。就在这一点上，您可以计算x，然后m'按地板计算x。它成为了;
√n + q 1 p p(√n - (q - pm')) p(√n + (pm' - q))
x = ______ = m' + ___ ⇒ x' = ______________ = _________________ = _________________
p x' √n + (q - pm') n - (q - pm')^2 n - (q - pm')^2
此时n - (q - pm')^2分母的 at 可以被分子的 at 整除p。现在这是稳定的，我们可以根据需要扩展它。让我们为qand分配新的任务p；
q' = pm'-q;
p' = (n - q'^2)/p;
√n + q'
x' = ______
p'
请注意，当p'变为 1 ( n - q'^2 = p) 时，我们处于回文的末尾。然而，为了决定在哪里停止，我使用了与我的 toRational 算法中描述的机制相同的机制，正如上面备选方案 1 中链接的那样。一旦达到 JS 浮点分辨率，它基本上就会停止。JavaScript 代码如下；
function rationalSqrt(n){
var nr = Math.sqrt(n),
m = Math.floor(nr),
p = n-m**2,
q = m,
cs = [m],
n0 = 1,
d0 = 0,
n1 = m,
d1 = 1,
n2 = 0,
d2 = 1;
if (nr === m) return {n:m,d:1,cs};
while (Math.abs(nr-n2/d2) > Number.EPSILON){
m = Math.floor((nr+q)/p);
q = m*p-q;
p = (n-q**2)/p;
cs.push(m);
n2 = m*n1+n0;
d2 = m*d1+d0;
n0 = n1;
d0 = d1;
n1 = n2;
d1 = d2;
}
return {n:n2,d:d2,cs};
}
这两种算法略有不同。
toRational在 JS 浮点分辨率内给出更小的分子和分母。rationalSqrt（这一次）在 JS 浮点分辨率内给出更小的分子和分母。rationalSqrt将产生精确的系数，正如人们期望的平方根一样，但一旦分辨率足够，就会被截断。toRational给出了预期的 couse 系数，但最后一个可能与您对平方根系列的期望完全无关。一个这样的例子是；
rationalSqrt(511); //returns
{ n : 882184734
, d : 39025555
, cs: [22,1,1,1,1,6,1,14,4,1,21,1,4,14,1,6,1]
}
尽管
toRational(Math.sqrt(511));
{ n : 1215746799
, d : 53781472
, cs: [22,1,1,1,1,6,1,14,4,1,21,1,4,14,1,10]
}
进一步的想法：考虑给我们RCF系数。我们可以反转rationalSqrt算法来获得它的(q + √n) / p形式吗？这可能是一项有趣的任务。
我使用 Surd Storage 类型来获得 n 的平方根的无限精度。
(b * \sqrt(n) + d)/c
=
(b * c * sqrt(n) - c * d + a_i * c^2) / (b^2 * n - d^2 - (a_i * c)^2 + 2* a_i * c * d)
sqrt(n) 的底值只使用一次。之后剩余的迭代存储为 surd 类型。这避免了其他算法中出现的舍入误差，并且可以实现无限（内存受限）分辨率。
a_0 = sqrt (n) 的底值
a_i = (b_i * a_0 + d_i) / c_i
b_i+1 = b_i * c
c_i+1 = (b_i)^2 * n - (d_i)^2 - (a_i * c_i)^2 + 2 * a_i * c_i * d_i
d_i+1 = a_i * (c_i)^2 - c_i * d_i
g = gcd(b_i+1 , c_i+1 , d_i+1)
b_i+1 = b_i+1 / g
c_i+1 = c_i+1 / g
d_i+1 = d_i+1 / g
a_i+1 = (b_i+1 * x + d_i+1) / c_i+1
然后对于 i=0 到 i=Maximum_terms 产生一个以 [a_0;a_1,a_2 ... ,2*a_0] 开头的连分数
当 a_i 项等于 a_0 的 2 倍时，我终止分数。这是序列重复的点。
数学是由 Electro World 完成的，关于数学的非常好的视频可以在这里找到 https://youtu.be/GFJsU9QsytM
下面提供了用 Java 编写的 BigInteger 源代码。希望你喜欢。
如果找到重复序列，则返回布尔值 true；如果未找到所需精度的重复序列，则返回布尔值。
可以根据Maximum_terms轻松修改精度以适应。
平方根 139 [11;1,3,1,3,7,1,1,2,11,2,1,1,7,3,1,3,1,22] 重复长度 18
15 [3;1,6] 的平方根重复长度 2
2501 [50;100] 的平方根重复长度 1
10807 的平方根 [103;1,22,9,2,2,5,4,1,1,1,6,15,1,5,2,1,3,6,34,2,34, 6,3,1,2,5,1,15,6,1,1,1,4,5,2,2,9,22,1,206] 重复长度 40
一个可能的两倍加速将是查看系列的回文性质。在本例中为 34、2、34。只有一半的序列需要确定。
public static Boolean SquareRootConFrac(BigInteger N) {
BigInteger A,B=BigInteger.ONE,C=B,D=BigInteger.ZERO;
BigInteger A0=N.sqrt(),Bi=B,Ci=C,Di=D,G;
BigInteger TwoA0 = BigInteger.TWO.multiply(A0);
int Frac_Length=0, Maximum_terms=10000; //Precision 10000 terms
String str="";
Boolean Repeat=false, Success=false, Initial_BCD=true;
while(!Repeat) {
Frac_Length++; Success=!(Frac_Length==Maximum_terms);
A=((B.multiply(A0)).add(D)).divide(C); Repeat=A.equals(TwoA0)||!Success;
Bi=B.multiply(C);
Ci=(B.multiply(B).multiply(N)).subtract(D.multiply(D)).subtract(A.multiply(A).multiply(C).multiply(C)).add(BigInteger.TWO.multiply(A).multiply(C).multiply(D));
Di=(A.multiply(C).multiply(C)).subtract(C.multiply(D));
G=Bi.gcd(Ci).gcd(Di);
B=Bi.divide(G);C=Ci.divide(G);D=Di.divide(G);
if(Initial_BCD) {str="["+A+";";System.out.print(str);Initial_BCD=false;}
else {str=""+A;System.out.print(str);if(!Repeat){str=",";System.out.print(str);}}
}
str="]";System.out.println(str);
str="repeat length ";System.out.print(str);
if(Success) {str=""+(Frac_Length-1);System.out.println(str);}
else {str="not found";System.out.println(str);}
return Success;
}