logo标签列表

子线性时间内求解第n个斐波那契数。

performancealgorithmmathtime-complexityfibonacci
77

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

4
可以说这与算法有关，因为原帖中含糊地提到了算法复杂度...... 不过我仍然很想知道是哪个算法。 - Matthew Scharley
2
下面的两个答案给出了正确的公式。至于这个问题是否与编程有关：它是计算机科学的一部分。用于推导公式的装置称为“生成函数”，在算法分析中扮演着重要的角色。 - azheglov
1
@azheglov：虽然生成函数很有用，但是不需要它们来推导出斐波那契数列的闭合形式表达式。 - jason
7
有时为了某种原因，您希望高效地解决问题。有时所需的洞察力可能是一种新的实现，有时是算法，有时是数学。每当出现后者时，没有必要将情况贬低为“与编程无关”。 - ShreevatsaR
7
结果的大小与n成线性关系。因此，不存在这样的算法。当然，这并不否定下面计算斐波那契数列使用O(log n)算术运算的好答案。 - Accipitridae
显示剩余3条评论
17个回答
101

参考 Pillsy 对矩阵指数的讨论，假设给定矩阵

M = [1 1] 
    [1 0]

那么 

fib(n) = Mn1,2

将矩阵连乘以产生幂并不高效。

有两种方法可以进行矩阵指数：分而治之法可在 O(ln n) 步内得到 Mn，或特征值分解法可以在常数时间内得到结果，但由于浮点精度有限可能会出现错误。

如果您希望获得大于浮点实现精度的确切值，则必须使用基于以下关系的 O ( ln n ) 方法：

Mn = (Mn/2)2 如果 n 是偶数
   = M·Mn-1 如果 n 是奇数

对 M 进行特征值分解可找到两个矩阵 U 和 Λ，使得 Λ 是对角矩阵，并且有 

 M  = U Λ U-1 
 Mn = ( U Λ U-1) n
    = U Λ U-1 U Λ U-1 U Λ U-1 ... n 次
    = U Λ Λ Λ ... U-1 
    = U Λ n U-1
将对角矩阵 Λ 的幂次方为 n，只需要将 Λ 中的每个元素取 n 次方即可，因此这提供了一种 O(1) 方法来计算 M 的 n 次方。但是，Λ 中的值不太可能是整数，因此会出现一些误差。

定义我们2x2矩阵的Λ为

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

要找到每个λ，我们解决

 |M - λI| = 0

得到

 |M - λI| = -λ ( 1 - λ ) - 1
λ² - λ - 1 = 0

使用二次公式 

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

如果你已经读了Jason的答案，你就会知道发生了什么。

解出特征向量X1和X2:

如果 X1 = [ X1,1, X1,2 ]
 M.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   Φ ]
令 Ψ = 1-Φ，另一个特征值
因为 Φ 是方程 λ²-λ-1=0 的根
所以  -ΨΦ = Φ²-Φ = 1
并且 Ψ+Φ = 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 ]
= M

因此，检验结果是正确的。

现在我们有了计算Mn1,2所需的一切：

Mn = 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 ] 因为ΨΦ=-1
= ( √5 )-1 [ Φn+1-Ψn+1      Φn-Ψn ]
              [ Φn-Ψn      Φn-1-Ψn-1 ]

所以

 fib(n) = Mn1,2
        = ( Φn - (1-Φ)n ) / √5

这与其他地方给出的公式相符。

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

+1 - 像往常一样，真是太棒了。你用什么排版工具来排版的呢？LaTeX吗？ - duffymo
这是从Gilbert Strang的《代数学》或其他优秀的线性代数书籍中复制粘贴的内容。 - alinsoar
1
@alinsoar 这不是“复制粘贴”，而是为了检查我是否还记得我的Lin A，参考了Open University课程笔记和维基百科。 - Pete Kirkham
我跟着 Gilbert Strang 学习了线性代数课程，那里面也有这个问题。确实，通过矩阵分解表达递归的问题是经典的，可以在任何好的教材/课程中找到。 - alinsoar
66

斐波那契数列的第n项为

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

其中

phi = (1 + sqrt(5)) / 2
假设基本的数学运算符（+、-、* 和 /）的时间复杂度为 O(1)，你可以使用这个结果以 O(log n) 的时间复杂度计算第 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);
}
7
@Json 我没有给你的回答点踩，但其他人可能会这样做，因为你的回答暗示第N个斐波那契数可以在O(log n)时间内计算，这是错误的。你的代码正在计算一个近似值。在任意精度下，你的代码至少需要O(n)的时间，因为答案的长度是O(n)。 - PeterAllenWebb
10
提供的公式不是近似值。第n个斐波那契数等于phi^n / sqrt(5) + 1/2的下取整，其中phi = (1 + sqrt(5)) / 2。这是一个事实。其次，我理解其他人关于答案长度为O(n)的观点，但我已经在我的答案中添加了一条注释，假设原始的数学运算花费恒定时间（我知道除非限制输入，否则它们不是）。我的观点是我们可以用O(log n)个算术运算找到第n个斐波那契数。 - jason
4
@Jason：假设幂运算的时间复杂度也是O(1)，这将使整个算法的时间复杂度变为O(1)。虽然这很好，但实际上幂运算和其他基本数学操作的时间复杂度都不是O(1)。因此简言之，这个公式虽好，但不能在次线性时间内计算结果。 - yairchu
12
@Jason：这个公式不是一个近似值，但是代码是一个近似值（除了在想象中的C#实现中，Math.Pow(…)具有无限精度的情况下，此时代码是O(n)）。 - ShreevatsaR
15
@Jason: 不行。在n=1000的情况下运行您的代码（斐波那契数[43466...849228875] （http://www.gutenberg.org/cache/epub/302/pg302.txt）仅有209位数字），然后告诉我您是否能完全正确地得到所有数字。要使Math.Floor获得正确的整数部分，必须通过Math.Pow准确计算那么多位数字。事实上，在我的C ++实现中，即使可以用long long精确表示整数130496954492865，也会计算出16位F_{74} = 130496954492865不正确，如果C＃得到的数字超过这个数量，我会感到惊讶。 - ShreevatsaR
显示剩余11条评论
57

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

是不可能的！

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

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

fib n ~= phin/√5

fib n 有多少位数字？

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

numDigits (fib n) = n * const + const

它是O（n）

由于所请求的结果是O（n）, 因此它不能在小于O（n）的时间内计算。

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

2
@yairchu：让我重新表述一下，如果我理解正确的话。理论上，计算fib_n需要计算n位数字，因此对于任意的n，它将花费O(n)的时间。然而，如果fib_n < sizeof(long long)，那么我们可以在O(log n)的时间内计算fib_n，因为机器架构提供了一个并行设置位的机制。（例如，int i = -1;需要设置32位，但在32位机器上，所有32位可以在恒定时间内设置。 - Sumit
8
如果你只想支持适合32位的结果，那么你也可以为这个序列的前48个结果建立查找表。显然这是O(1)，但是：对于有界N进行大O分析是愚蠢的，因为你总是可以将任何内容合并到常量因子中。因此，我的答案是针对无界输入的。 - yairchu
1
@yairchu：你能演示一下你的逻辑吗？以一个众所周知的例子为例，比如对于一个由n个数字组成的序列进行基于比较的排序，其中每个数字有O(log n)位，时间复杂度为O(n*log n)。 - jfs
1
这取决于你对“时间”的理解是什么。对于排序（或哈希表查找），“时间”意味着比较次数。在问题中，它可能意味着算术运算。在这个答案中，它被理解为类似逐位操作的东西。 - Paul Hankin
4
在以sqrt(2)为底数的进位制下，整数确实会有有限的表示，但在奇数位上它只是零，相当于二进制。如果以sqrt(2)为底数的任何奇数位上不为零，那么你就得到了一个无理数。一个可能需要用到以黄金比例为底数的情况是在ADC将连续信号转换为模拟信号时。据我所知，这是黄金比例“工业”应用的一种，用于在舍入信号时减少粗略化。个人而言，我使用以黄金比例和斐波那契编码作为一种符号方便的方式来处理辫群的斐波那契任意子表示。 - saolof
显示剩余9条评论
35

在 SICP 中，有一个与此相关的练习，答案在这里。

用命令式风格编写的程序将类似于：

函数 Fib(count)
    a ← 1
    b ← 0
    p ← 0
    q ← 1
 当 count > 0 时做
        如果 偶数(count) 则
             p ← p² + q²
             q ← 2pq + q²
             count ← count ÷ 2
        否则
             a ← bq + aq + ap
             b ← bp + aq
             count ← count - 1
        结束如果
    结束循环

    返回 b
结束函数
这是一个Python的实现（可与“twisted”框架一起使用），链接如下：https://github.com/zed/txfib/blob/41ea022cc8cffc1d4b63996d313e644b494be7dd/fibonacci.py#L139 - jfs
"If Even(count) Then" 应该翻译为 "如果奇数(count)则"。 - Md. Monirul Islam
@MonirulIslamMilon if even(count) 是正确的。该序列从零开始（第零个斐波那契数是零）：0,1,1,2,3,5,8,13,... - jfs
书籍链接现在为：https://mitpress.mit.edu/sites/default/files/sicp/full-text/book/book-Z-H-11.html#%_thm_1.19 - Lee D
晚了点评论，但变量p和a在用于计算q和b之前被覆盖。为避免此问题，预先计算项并更改p和q分配的顺序: | qq = q·q | q = 2·p·q + qq | p = p·p + qq | ... | aq = a·q | a = b·q + aq + a·p | b = b·p + aq |. - rcgldr
24

您也可以通过对一个整数矩阵进行指数运算来达到同样的效果。如果您有以下矩阵：

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

如果[]是矩阵下标符号，^是矩阵幂运算符，则 (M^n)[1,2] 将等于第 n 个斐波那契数。对于固定大小的矩阵，幂运算可以像实数一样在O(log n)时间内完成。

编辑： 当然，根据您想要的回答类型，您可能能够使用常量时间算法轻松解决问题。正如其他公式所示，第 n 个斐波那契数呈指数增长。即使使用64位无符号整数，您只需要一个94个条目的查找表来涵盖整个范围。

第二次编辑：首先进行特征分解来做矩阵幂等价于JDunkerly在下面的解法。该矩阵的特征值为(1 + sqrt(5))/2和(1 - sqrt(5))/2。

3
利用M的特征值分解来高效地计算M的n次幂。 - Pete Kirkham
1
提议的方法适用于整数计算（可能需要长算术）。使用特征分解的方法并不有趣：如果您不需要整数计算，则可以使用Jason答案中的公式。 - Konstantin Tenzin
1
@Konstantin，Jason答案中的公式是通过特征分解得出的结果，所以你自相矛盾了。 - Pete Kirkham
@Pete Kirkham 这个公式可以通过几种方法得到：特征方程，特征分解，归纳证明。我不确定特征分解是最简单的方法。无论如何，它是众所周知的，而且立即使用它更容易。 - Konstantin Tenzin
5

维基百科提供了斐波那契数列的封闭式解法。 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;
    }
2
当n为非负整数时，您可以通过使用“|1-phi|^n / sqrt(5) < 1/2”的事实来避免计算两个指数。 - jason
不知道调整总是使用另一种形式，但这是一个不错的优化。 - JDunkerley
1
正确解决方案的近似结果涉及矩阵乘法。 - cerkiewny
4
对于非常大的数字，可以使用以下递归函数来计算。它使用以下等式：
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库。
首先，用斐波那契数列的数组开始。将fibs[0]=0，fibs[1]=1，fibs[2]=1，fibs[3]=2，fibs[4]=3等等。在这个例子中，我使用了前501个（从0开始计数）的数组。您可以在此处找到前500个非零斐波那契数： http://home.hiwaay.net/~jalison/Fib500.html。需要一些编辑才能将其放入正确的格式中，但这并不太难。
然后，您可以使用以下函数（在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);
}

我已经测试了25000个斐波那契数及其类似数字。

这段代码不太高效。假设fibs[]数组只有10个元素，你调用Fib(101)。Fib(101)会调用Fib(51)和Fib(50)。Fib(51)会调用Fib(26)和Fib(25)。Fib(50)会调用Fib(25)和Fib(24)。所以Fib(25)被调用了两次，造成了浪费。即使使用了高达500的fibs，Fib(100000)也会出现同样的问题。 - Eyal
3

这是我的递归版本，递归了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)

这段话是关于 IT 技术的内容。它讲述了如何计算斐波那契数列，具体方法是：如果 n 是奇数，则可以使用 fib(n-1) 和 fib(n-2) 来计算 fib(n) 和 fib(n-1)；如果 n 是偶数，则可以使用 fib(n/2) 和 fib(n/2-1) 来计算 fib(n) 和 fib(n-1)。这个方法的原理是利用了斐波那契数列的规律。要推导出偶数情况，需要将连续的斐波那契数值写成一个矩阵，其中 a = b+c。基本情况和奇数情况比较简单。
[1 1] * [a b]  =  [a+b a]
[1 0]   [b c]     [a   b]

从这个公式中，我们可以看出前三个斐波那契数列的矩阵乘以任意三个连续的斐波那契数列的矩阵等于下一个斐波那契数。因此我们知道：

      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)]

简化右侧的内容会导致偶数情况。
我想在这里强调，您希望根据F(n)和F(n-1)的值计算出F(2n)和F(2n+1)的值。您没有说明您想要做什么。 - alinsoar
1

以下是一行代码，使用O(n)大小的整数，在O(log n)算术运算中计算F(n)：

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

使用大小为 O(n) 的整数是合理的，因为它与答案的大小相当。
要理解这一点，让 phi 成为黄金比率（x^2=x+1 的最大解），F(n) 是第 n 个斐波那契数，其中 F(0)=0，F(1)=F(2)=1。
现在，phi^n = F(n-1) + F(n)phi。
归纳证明：phi^1 = 0 + 1*phi = F(0) + F(1)phi。如果 phi^n = F(n-1) + F(n)phi，则 phi^(n+1) = F(n-1)phi + F(n)phi^2 = F(n-1)phi + F(n)(phi+1) = F(n) + (F(n)+F(n-1))phi = F(n) + F(n+1)phi。这个计算中唯一棘手的步骤是将 phi^2 替换为 (1+phi)，这是因为 phi 是黄金比率。
另外，形如 (a+b*phi) 的数字，其中 a、b 是整数，在乘法下是封闭的。
证明：(p0+p1*phi)(q0+q1*phi) = p0q0 + (p0q1+q1p0)phi + p1q1*phi^2 = p0q0 + (p0q1+q1p0)phi + p1q1*(phi+1) = (p0q0+p1q1) + (p0q1+q1p0+p1q1)*phi。
利用这个表示方法，可以使用平方取幂法在O(log n)个整数运算中计算出phi^n。结果将是F(n-1)+F(n)phi，从中可以读出第n个斐波那契数。
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])

请注意，这段代码的大部分是标准的平方指数函数。
要得到本答案开始的一行代码，可以注意到phi可以用足够大的整数X来表示，可以将(a+b*phi)(c+d*phi)作为整数运算(a+bX)(c+dX)模(X^2-X-1)。然后pow函数可以被标准Python pow函数所替换（它方便地包括第三个参数z，该参数计算结果模z）。所选择的X是2 << i。
1
固定点算术不准确。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。
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接