生成平方根的连续分数

我写了这段代码来生成平方根 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*Q²代替D，用P*Q代替P，用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;
}

它很容易计算周期长度为 182 的 √7981的连续分数(例如(。

罪魁祸首不是floor。罪魁祸首是计算A= 1.0 / (A - B); 深入挖掘，罪魁祸首是您的计算机用来表示实数的 IEEE 浮点机制。减法和加法会失去精度。在算法反复执行时重复减法会失去精度。

当你计算出连续分数项 {11,1,3,1,3,7,1,1,2,11,2} 时，你的 IEEE 浮点值 A 只好到六位，而不是人们期望的十五或十六位。当你到达{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 函数，算法略有不同。

但是，问题还没有解决。是的，通常，√n的 CF 系数采用重复回文的形式，以第一个非零系数的双倍结束。比如√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]则得到一个糟糕的近似值(如√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]

在这种情况下，一旦获得必要的精度，就停止迭代是合理的。正如我在开头提到的，有两种方法。

到有理数算法从任何十进制数中获取简单的连续分数。这将使 CF 精确解析为该小数，而没有任何浮点错误。有关此算法的详细说明，您可以查看我以前的答案。
碰巧有一种更直接的√n算法，它基本上与 1 相同，但针对平方根进行了调整。在这种情况下，您不提供√n而是n。

思路如下。我们必须为输入定义一个通用形式，它涉及分子处的平方根值。然后，我们尝试在连续分数部分中达到相同的表达式，以便能够迭代。

让我们的输入是以下形式

q + √n
______
   p

对于简单的平方根运算，我们可以假设q是0的，p是1的。如果我们能在下一阶段建立这种形式，那么我们就可以很容易地迭代。

从初始阶段开始，q = 0、p = 1、m是√n的整数部分，1/x是浮点部分，我们的目标是将x变成(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计算x和下一个m。算法的广义形式变为;

    √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

此时p可被n - (q - pm)^2整除。这现在是稳定的，我们可以根据需要扩展它。让我们为q和p做新的任务;

q' = pm-q;
p' = (n - q'^2)/p;
     √n + q'
x' = ______
        p'
m' = Math.floor(x')

注意，当p'变成1(n - q'^2 = p(时，我们就在回文的末尾。但是，为了决定在哪里停止，我使用的是与上面的备选方案 1 中链接的 toRational 算法中描述的相同机制。一旦达到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};
}

这两种算法略有不同。

~60% 的时间它们会产生相同的分数。
~27% 的时间toRational在 JS 浮点分辨率内给出较小的分子和分母。但是，如果我们在 JS 中多一个数字，那肯定不是一个更好的近似值。
~13% 的时间rationalSqrt(这个(在 JS 浮点分辨率内给出较小的分子和分母。
rationalSqrt将产生精确的系数，正如人们期望的平方根一样，但一旦分辨率足够，就会被截断。
toRational给出了预期的系数，但最后一个系数可能与您对平方根级数的期望完全无关。

一个这样的例子是;

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形式？这可能是一项有趣的任务。
使用 JS 十进制分辨率实际限制 CF 将在 CF 算术中产生类似的精度，您将浪费时间。为了获得更高的平方根分辨率，我建议while (Math.abs(nr-n2/d2) > Number.EPSILON)条件替换为while (cs.length < 30)如果 30 CF 分辨率系数足以满足您的应用。这个数字取决于你。

我使用了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,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;
}

相关内容

最新更新

热门标签：