亚线性时间内的第 n 个斐波那契数

StackOverflow https://stackoverflow.com/questions/1525521

题

是否有任何算法可以在亚线性时间内计算第 n 个斐波那契数？

有帮助吗？

解决方案

在nth斐波那契数由

给定 
f(n) = Floor(phi^n / sqrt(5) + 1/2)

，其中

phi = (1 + sqrt(5)) / 2

假设原始数学运算（+，-，*和/）的O(1)可以使用该结果来计算在n时间O(log n)th Fibonacci数（因为式中的求幂的O(log n)）。

在C＃：

static double inverseSqrt5 = 1 / Math.Sqrt(5);
static double phi = (1 + Math.Sqrt(5)) / 2;
/* should use 
   const double inverseSqrt5 = 0.44721359549995793928183473374626
   const double phi = 1.6180339887498948482045868343656
*/

static int Fibonacci(int n) {
    return (int)Math.Floor(Math.Pow(phi, n) * inverseSqrt5 + 0.5);
}

其他提示

遵循 Pillsy 对矩阵求幂的参考，这样对于矩阵

中号 = [1 1] 
    [1 0]

然后 

谎言(n) = 中号n1,2

使用重复乘法求矩阵幂的效率不是很高。

矩阵求幂的两种方法是分而治之，得出 中号n 在 氧(ln）步骤，或特征值分解，其时间恒定，但由于浮点精度有限，可能会引入误差。

如果您想要一个大于浮点实现精度的精确值，则必须使用基于以下关系的 O ( ln n ) 方法： 

中号n = (中号n/2)2 if n even
   = 中号·中号n-1 if n is odd

特征值分解 中号 找到两个矩阵 U 和 Λ 这样 Λ 是对角线并且 

 中号  = U Λ U-1 
 中号n = ( U Λ U-1) n
    = U Λ U-1 U Λ U-1 U Λ U-1 ... n times
    = U Λ Λ Λ ... U-1 
    = U Λ n U-1
求对角矩阵 Λ 到 n次幂是一个简单的问题，即提高每个元素 Λ 到 nth，所以这给出了一个 O(1) 的提升方法 中号 到 n次幂。然而，中的值 Λ 不太可能是整数，因此会出现一些错误。

定义 Λ 对于我们的 2x2 矩阵

Λ = [ λ1 0 ]
  = [ 0 λ2 ]

去寻找每一个 λ, ，我们解决 

 |中号 - λ我| = 0

这使

 |中号 - λ我| = -λ ( 1 - λ ) - 1

λ² - λ - 1 = 0

使用二次公式 

λ    = ( -b ± √ ( b² - 4ac ) ) / 2a
     = ( 1 ± √5 ) / 2
 { λ1, λ2 } = { Φ, 1-Φ } where Φ = ( 1 + √5 ) / 2

如果您读过 Jason 的回答，您就可以知道接下来会发生什么。

求解特征向量 X1 和 X2:

if X1 = [ X1,1, X1,2 ]

 中号.X1 1 = λ1X1

 X1,1 + X1,2 = λ1 X1,1
 X1,1      = λ1 X1,2

=>
 X1 = [ Φ,   1 ]
 X2 = [ 1-Φ, 1 ]

这些向量给出 U:

U = [ X1,1, X2,2 ]
    [ X1,1, X2,2 ]

  = [ Φ,   1-Φ ]
    [ 1,   1   ]

反相 U 使用

A   = [  a   b ]
      [  c   d ]
=>
A-1 = ( 1 / |A| )  [  d  -b ]
                   [ -c   a ]

所以 U-1 是（谁）给的

U-1 = ( 1 / ( Φ - ( 1 - Φ ) )  [  1  Φ-1 ]
                               [ -1   Φ  ]
U-1 = ( √5 )-1  [  1  Φ-1 ]
               [ -1   Φ  ]

完整性检查：

UΛU-1 = ( √5 )-1 [ Φ   1-Φ ] . [ Φ   0 ] . [ 1  Φ-1 ] 
                     [ 1   1  ]   [ 0  1-Φ ]   [ -1   Φ ]

let Ψ = 1-Φ, the other eigenvalue

as Φ is a root of λ²-λ-1=0 
so  -ΨΦ = Φ²-Φ = 1
and Ψ+Φ = 1

UΛU-1 = ( √5 )-1 [ Φ   Ψ ] . [ Φ   0 ] . [  1  -Ψ ] 
                 [ 1   1 ]   [ 0   Ψ ]   [ -1   Φ ]

       = ( √5 )-1 [ Φ   Ψ ] . [ Φ   -ΨΦ ] 
                 [ 1   1 ]   [ -Ψ  ΨΦ ]

       = ( √5 )-1 [ Φ   Ψ ] . [ Φ    1 ] 
                 [ 1   1 ]   [ -Ψ  -1 ]

       = ( √5 )-1 [ Φ²-Ψ²  Φ-Ψ ] 
                  [ Φ-Ψ      0 ]

       = [ Φ+Ψ   1 ]    
         [ 1     0 ]

       = [ 1     1 ] 
         [ 1     0 ]

       = 中号

所以健全性检查成立。

现在我们已经拥有了计算所需的一切 中号n1,2:

中号n = UΛnU-1
   = ( √5 )-1 [ Φ   Ψ ] . [ Φn  0 ] . [  1  -Ψ ] 
              [ 1   1 ]   [ 0   Ψn ]   [ -1   Φ ]

   = ( √5 )-1 [ Φ   Ψ ] . [  Φn  -ΨΦn ] 
              [ 1   1 ]   [ -Ψn   ΨnΦ ]

   = ( √5 )-1 [ Φ   Ψ ] . [  Φn   Φn-1 ] 
              [ 1   1 ]   [ -Ψn  -Ψn-1 ] as ΨΦ = -1

   = ( √5 )-1 [ Φn+1-Ψn+1      Φn-Ψn ]
              [ Φn-Ψn      Φn-1-Ψn-1 ]

所以

 谎言(n) = 中号n1,2
        = ( Φn - (1-Φ)n ) / √5

这与其他地方给出的公式一致。

您可以从递归关系中推导出它，但在工程计算和仿真中，计算大型矩阵的特征值和特征向量是一项重要的活动，因为它给出了方程组的稳定性和谐波，并允许有效地将矩阵提升为高幂。

如果你想要确切的数字（这是一个“bignum”，而不是一个 int/float），那么恐怕

不可能！

如上所述，斐波那契数列的公式为：

fib n = 下限 (phin/√5 + 1/2)

fib n ~= phin/√5

是多少位数字 fib n?

numDigits (fib n) = log (fib n) = log (phin/√5) = 对数 phin - log √5 = n * log phi - log √5

numDigits (fib n) = n * const + const

它是 氧(n)

由于请求的结果是 氧(n），不能用小于来计算 氧(n） 时间。

如果您只想要答案的低位数字，则可以使用矩阵求幂方法在亚线性时间内进行计算。

中的一个 SICP练习 是关于这个的，里面有答案描述 这里。

在命令式风格中，程序看起来像这样

功能 斐波那契(数数)
    A ← 1
    乙 ← 0
    p ← 0
    q ← 1

    尽管 数数 > 0 做
        如果 Even(数数) 然后
             p ← p² + q²
             q ← 2pq + q²
             数数 ← 数数 ÷ 2
        别的
             A ← 巴克 + 水 + 美联社
             乙 ← BP + 水
             数数 ← 数数 - 1
        万一
    结束时

    返回 乙
结束功能

可以通过幂整数的矩阵，以及这样做。如果你有矩阵

    / 1  1 \
M = |      |
    \ 1  0 /

然后(M^n)[1, 2]将是等于nth斐波那契数，如果[]是一个矩阵标和^是矩阵求幂。对于一个固定大小的矩阵，幂一个正整数功效可以在O（log n）的时间以同样的方式与实数来完成。

编辑：当然，这取决于你想要的答案的类型，你可以逃脱一个常数时间算法。像其他公式显示，该nth Fibonacci数呈指数增长n。即使使用64位无符号整数，你只需要一个94输入的查找表，以便覆盖整个范围。

<强> SECOND编辑：否则矩阵指数与第一特征分解是完全等同于下面JDunkerly的溶液。这个矩阵的特征值是(1 + sqrt(5))/2和(1 - sqrt(5))/2。

维基百科具有闭合形式解 http://en.wikipedia.org/wiki/Fibonacci_number

或者在C＃：

    public static int Fibonacci(int N)
    {
        double sqrt5 = Math.Sqrt(5);
        double phi = (1 + sqrt5) / 2.0;
        double fn = (Math.Pow(phi, N) - Math.Pow(1 - phi, N)) / sqrt5;
        return (int)fn;
    }

有关非常大的，这个递归函数的工作。它使用以下等式：

F(2n-1) = F(n-1)^2 + F(n)^2
F(2n) = (2*F(n-1) + F(n)) * F(n)

您需要一个库，让你大整数的工作。我使用 https://mattmccutchen.net/bigint/ 中的BigInteger库。

开始与斐波那契数的阵列。使用的FIB [0] = 0，FIB中[1] = 1，的FIB [2] = 1，FIB中[3] = 2，的FIB [4] = 3，等等。在本例中，我使用的第一个501的阵列（计数为0）。你可以在这里找到的第一个500非零斐波那契数： http://home.hiwaay.net /~jalison/Fib500.html 。这需要一点点的编辑把它以正确的格式，但不是太硬了。

后，可以看到在使用该功能的任何Fibonacci数（在C）：

BigUnsigned GetFib(int numfib)
{
int n;
BigUnsigned x, y, fib;  

if (numfib < 501) // Just get the Fibonacci number from the fibs array
    {
       fib=(stringToBigUnsigned(fibs[numfib]));
    }
else if (numfib%2) // numfib is odd
    {
       n=(numfib+1)/2;
       x=GetFib(n-1);
       y=GetFib(n);
       fib=((x*x)+(y*y));
    }
else // numfib is even
    {
       n=numfib/2;
       x=GetFib(n-1);
       y=GetFib(n);
       fib=(((big2*x)+y)*y);
   }
return(fib);
}

我为第25,000 Fibonacci数等的测试此。

下面是我认为递归的log（n）次递归版本。我认为这是最简单的调用形式为：

def my_fib(x):
  if x < 2:
    return x
  else:
    return my_fib_helper(x)[0]

def my_fib_helper(x):
  if x == 1:
    return (1, 0)
  if x % 2 == 1:
    (p,q) = my_fib_helper(x-1)
    return (p+q,p)
  else:
    (p,q) = my_fib_helper(x/2)
    return (p*p+2*p*q,p*p+q*q)

它的工作原理，因为可以使用fib(n),fib(n-1)如果n是奇数，当n为偶数，则可以使用fib(n-1),fib(n-2)计算fib(n),fib(n-1)计算fib(n/2),fib(n/2-1)。

基础案例和奇数的情况下是简单的。为了导出偶数的情况下，开始一个，B，C作为连续斐波纳契值（例如，8,5,3），并在基质中写，其中a = B + C。注意：

[1 1] * [a b]  =  [a+b a]
[1 0]   [b c]     [a   b]

从这一点，我们看到，前三个Fibonacci数的矩阵，倍任何三个连续Fibonacci数的矩阵，等于下。因此，我们知道：

      n
[1 1]   =  [fib(n+1) fib(n)  ]
[1 0]      [fib(n)   fib(n-1)]

所以：

      2n                        2
[1 1]    =  [fib(n+1) fib(n)  ]
[1 0]       [fib(n)   fib(n-1)]

简化右手侧导致甚至情况。

使用ř 

l1 <- (1+sqrt(5))/2
l2 <- (1-sqrt(5))/2

P <- matrix(c(0,1,1,0),nrow=2) #permutation matrix
S <- matrix(c(l1,1,l2,1),nrow=2)
L <- matrix(c(l1,0,0,l2),nrow=2)
C <- c(-1/(l2-l1),1/(l2-l1))

k<-20 ; (S %*% L^k %*% C)[2]
[1] 6765

除了由数学方法，最好的最优解（I相信），以避免重复计算使用的字典。的一个微调

import time

_dict = {1:1, 2:1}

def F(n, _dict):
    if n in _dict.keys():
        return _dict[n]
    else:
        result = F(n-1, _dict) + F(n-2, _dict)
        _dict.update({n:result})
        return result

start = time.time()

for n in range(1,100000):
    result = F(n, _dict) 

finish = time.time()

print(str(finish - start))

我们先从琐碎字典（斐波纳契数列的头两个值）和不断增加斐波纳契值到词典中。

花了大约0.7秒的第一100000个斐波纳契值（的Intel Xeon CPU E5-2680 @ 2.70千兆赫，16 GB RAM，视窗10-64位OS）

定点算术是不准确的。 Jason的C＃代码给出对于n不正确的答案= 71（308061521170130代替308061521170129）及以后。

有关正确答案，可以使用一个计算代数系统。 Sympy就是Python这样的库。有在 http://live.sympy.org/ 交互式控制台。复制并粘贴该功能

phi = (1 + sqrt(5)) / 2
def f(n):
    return floor(phi**n / sqrt(5) + 1/2)

然后计算

>>> f(10)
55

>>> f(71)
308061521170129

您可能想尝试检查phi。

您可以用怪异的树根平方公式得到一个确切的答案。其原因是，$ \的sqrt（5）$末下降了，你只需要跟踪系数用自己的乘法格式。

def rootiply(a1,b1,a2,b2,c):
    ''' multipy a1+b1*sqrt(c) and a2+b2*sqrt(c)... return a,b'''
    return a1*a2 + b1*b2*c, a1*b2 + a2*b1

def rootipower(a,b,c,n):
    ''' raise a + b * sqrt(c) to the nth power... returns the new a,b and c of the result in the same format'''
    ar,br = 1,0
    while n != 0:
        if n%2:
            ar,br = rootiply(ar,br,a,b,c)
        a,b = rootiply(a,b,a,b,c)
        n /= 2
    return ar,br

def fib(k):
    ''' the kth fibonacci number'''
    a1,b1 = rootipower(1,1,5,k)
    a2,b2 = rootipower(1,-1,5,k)
    a = a1-a2
    b = b1-b2
    a,b = rootiply(0,1,a,b,5)
    # b should be 0!
    assert b == 0
    return a/2**k/5

if __name__ == "__main__":
    assert rootipower(1,2,3,3) == (37,30) # 1+2sqrt(3) **3 => 13 + 4sqrt(3) => 39 + 30sqrt(3)
    assert fib(10)==55

下面是一个单行，计算F（n）时，使用尺寸为O（n）的整数，在O（log n）的算术运算：

for i in range(1, 50):
    print(i, pow(2<<i, i, (4<<2*i)-(2<<i)-1)//(2<<i))

使用大小为O（n）是合理的整数，因为这是相媲美的答案的大小。

要了解这一点，让披是黄金比例（最大的溶液为x ^ 2 = X + 1）和F（n）是第n个斐波那契数，其中F（0）= 0，F（1 ）= F（2）= 1

现在，披^ N = F（N-1）+ F（n）的披

  

通过感应证明：披^ 1 = 0 + 1 *披= F（0）+ F（1）披。如果披^ N =   F（N-1）+ F（n）的PHI，然后披^（N + 1）= F（N-1）披+ F（n）的披^ 2 = F（N-1）+披   F（n）的（PHI + 1）= F（N）+（F（N）+ F（N-1））φ= F（N）+ F（N + 1）披。在这个计算中唯一棘手的步骤是由（1个+ PHI），它遵循因为Phi是黄金比例代替披^ 2之一。

另外数字形式（A + B * PHI），其中a，b是整数，是根据乘法关闭的。

  

证明：（P0 + P1 * PHI）（Q0 + Q1 * PHI）= p0q0 +（p0q1 + q1p0）披+ p1q1 *披^ 2 =   p0q0 +（p0q1 + q1p0）披+ p1q1 *（PHI + 1）=（p0q0 + p1q1）+   （p0q1 + q1p0 + p1q1）*披

使用这种表示，可以计算披^ n的在O（log n）的整数运算使用指数通过平方。其结果将是F（N-1）+ F（n）的PHI，从中可以读出第n个Fibonacci数。

def mul(p, q):
    return p[0]*q[0]+p[1]*q[1], p[0]*q[1]+p[1]*q[0]+p[1]*q[1]

def pow(p, n):
    r=1,0
    while n:
        if n&1: r=mul(r, p)
        p=mul(p, p)
        n=n>>1
    return r

for i in range(1, 50):
    print(i, pow((0, 1), i)[1])

请注意，大多数这些代码是一个标准的乘幂逐平方函数。

要得到一衬里开始该答案，可以注意到，由一个足够大的整数X表示披，可以作为所述整数操作(a+b*phi)(c+d*phi)执行(a+bX)(c+dX) modulo (X^2-X-1)。然后pow功能可以通过标准Python pow函数来代替（其方便地包括一第三个参数z其计算结果的模z。选择的X是2<<i。

许可以下： CC-BY-SA 和 归因
不隶属于 StackOverflow
scroll top