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

Question

有没有什么算法可以在亚线性时间内计算出第 n 个斐波那契数？

Answer 1

根据 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和Λ使得Λ是对角的并且 p>
 M = U Λ U-1Mn = (U Λ U-1)  n
    = U Λ U-1U Λ U-1U Λ U-1 ... n 次
    = U Λ Λ Λ ... U-1
    = U Λ nU-1
将对角矩阵 Λ 提高到 n 次方只是将 Λ 中的每个元素提高到 n i>th，因此这给出了将 M 提高到 n 次方的 O(1) 方法。但是，Λ中的值不可能是整数，所以会出现一些错误。

为我们的 2x2 矩阵定义 Λ 为


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

要找到每个λ，我们解决

 |M - λI| = 0

给了

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

λ² - λ - 1 = 0

使用二次公式

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

如果你已经阅读了 Jason 的回答，你可以看到它的发展方向。

求解特征向量X₁和X₂：

如果 X₁ = [ X_1,1, X_{1 ,2} ] M.X_{1 1} = λ₁X1 sub>X_1,1 + X_1,2 = λ₁X_1,1X_1,1 = λ₁X_1,2 => X₁ = [ Φ, 1 ] X₂ = [ 1-Φ, 1 ]

这些向量给出U：

U = [ X_1,1, X_2,2 ] [ X_1,1, X_2,2 ] = [Φ, 1-Φ] [ 1, 1 ]

使用

反转U A = [ a b ] [ c d ] => A^-1 = ( 1 / |A| ) [ d -b ] [-c一]

所以 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

所以健全性检查成立。

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

Mⁿ = UΛⁿU^-1 = ( √5 )^-1 [ Φ Ψ ] 。 [ Φⁿ 0 ] 。 [ 1 -Ψ ] [ 1 1 ] [ 0 Ψⁿ ] [ -1 Φ ] = ( √5 )^-1 [ Φ Ψ ] 。 [ Φⁿ -ΨΦⁿ ] [ 1 1 ] [ -Ψⁿ ΨⁿΦ ] = ( √5 )^-1 [ Φ Ψ ] 。 [ Φⁿ Φ^n-1 ] [ 1 1 ] [ -Ψⁿ -Ψ^n-1 ] 作为 ΨΦ = -1 = ( √5 )^-1 [ Φⁿ⁺¹-Ψn+1 sup> Φⁿ-Ψⁿ ] [ Φⁿ-Ψⁿ Φn-1 sup>-Ψ^n-1 ]

所以

fib(n) = Mⁿ1,2子> = ( Φⁿ - (1-Φ)ⁿ ) / √5

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

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

Answer 2

nth 斐波那契数由下式给出

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

在哪里

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

假设原始数学运算（+、-、* 和 /）是 O(1)，您可以使用此结果在 O(log n) 时间计算 nth 斐波那契数（ 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);
}

Answer 3

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

不可能！

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

fib n = 楼层 (phiⁿ/√5 + ¹/₂)

fib n ~= phiⁿ/√5

fib n 是多少位数？

numDigits (fib n) = log (fib n) = log (phiⁿ/√5) = log phiⁿ - log √5 = n * log phi - 记录√5

numDigits (fib n) = n * const + const

是 O(n)

由于请求的结果是O(n)，所以不能小于O(n ) 时间。

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

Answer 4

其中一个exercises in SICP 是关于这个的，其中有here. 描述的答案

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

函数 Fib(count) 一个 ← 1 b ← 0 p ← 0 q ← 1 当 计数 > 0 做 如果 偶数(count) 那么 p ← p² + q² q ← 2pq + q² 计数 ← 计数 ÷ 2 其他 a ← bq + aq + ap b ← bp + aq 计数 ← 计数 - 1 如果结束 结束时 返回 b 结束函数

Answer 5

您也可以通过对整数矩阵求幂来做到这一点。如果你有矩阵

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

如果[] 是矩阵下标并且^ 是矩阵求幂，那么(M^n)[1, 2] 将等于nth 斐波那契数。对于固定大小的矩阵，正整数幂的取幂可以在 O(log n) 时间内以与实数相同的方式完成。

编辑：当然，根据您想要的答案类型，您也许可以使用恒定时间算法。像其他公式所示，nth 斐波那契数随着n 呈指数增长。即使使用 64 位无符号整数，您也只需要一个 94 条目的查找表即可覆盖整个范围。

第二次编辑：首先使用特征分解进行矩阵指数与下面的 JDunkerly 解决方案完全相同。这个矩阵的特征值是(1 + sqrt(5))/2和(1 - sqrt(5))/2。

Answer 6

维基百科有一个封闭形式的解决方案 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;
    }

Answer 7

对于真正大的，这个递归函数有效。它使用以下等式：

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);
}

我已经针对第 25,000 个斐波那契数等进行了测试。

Answer 8

这是我递归 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)

它之所以有效，是因为如果 n 是奇数，如果 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 作为连续的斐波那契值（例如，8,5,3）开始，并将它们写入矩阵，其中 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)]

简化右手边导致偶数情况。

Answer 9

使用R

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

Answer 10

定点运算不准确。 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。

Answer 11

这是一个单行程序，它使用大小为 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])

请注意，此代码的大部分是标准的平方乘幂函数。

要找到开始这个答案的单行，可以注意到用足够大的整数 X 表示 phi，可以将 (a+b*phi)(c+d*phi) 执行为整数运算 (a+bX)(c+dX) modulo (X^2-X-1)。然后pow 函数可以替换为标准 Python pow 函数（它方便地包括第三个参数z，它计算结果以z 为模。选择的X 是2<<i。

Answer 12

除了通过数学方法进行微调之外，最好的最佳解决方案之一（我相信）是使用字典以避免重复计算。

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

我们从琐碎的字典（斐波那契数列的前两个值）开始，并不断将斐波那契值添加到字典中。

前 100000 个斐波那契值大约需要 0.7 秒（Intel Xeon CPU E5-2680 @ 2.70 GHz，16 GB RAM，Windows 10-64 位操作系统）

Answer 13

见分治算法here

该链接包含此问题的其他一些答案中提到的矩阵求幂的伪代码。

Answer 14

您可以使用奇怪的平方根方程来获得准确的答案。原因是 $\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

Answer 15

我遇到了一些计算斐波那契的有效时间复杂度的方法，以下是其中的一些 -

方法 1 - 动态规划 现在这里的子结构是众所周知的，因此我将直接跳到解决方案 -

static int fib(int n) 
{ 
int f[] = new int[n+2]; // 1 extra to handle case, n = 0 
int i; 

f[0] = 0; 
f[1] = 1; 

for (i = 2; i <= n; i++) 
{ 
    f[i] = f[i-1] + f[i-2]; 
} 

return f[n]; 
}

上面的空间优化版本可以如下完成 -

static int fib(int n) 
 { 
    int a = 0, b = 1, c; 
    if (n == 0) 
        return a; 
    for (int i = 2; i <= n; i++) 
    { 
        c = a + b; 
        a = b; 
        b = c; 
    } 
    return b; 
} 

方法2-（利用矩阵{{1,1},{1,0}}的幂）

这是一个 O(n) 依赖于这样一个事实：如果我们将矩阵 M = {{1,1},{1,0}} 乘以 n 次（换句话说，计算 power(M, n ) )，然后我们得到第 (n+1) 个斐波那契数作为结果矩阵中行和列 (0, 0) 处的元素。这个解决方案需要 O(n) 时间。

矩阵表示给出了斐波那契数的以下封闭表达式： 斐波那契矩阵

static int fib(int n) 
{ 
int F[][] = new int[][]{{1,1},{1,0}}; 
if (n == 0) 
    return 0; 
power(F, n-1); 

return F[0][0]; 
} 

/*multiplies 2 matrices F and M of size 2*2, and 
puts the multiplication result back to F[][] */
static void multiply(int F[][], int M[][]) 
{ 
int x = F[0][0]*M[0][0] + F[0][1]*M[1][0]; 
int y = F[0][0]*M[0][1] + F[0][1]*M[1][1]; 
int z = F[1][0]*M[0][0] + F[1][1]*M[1][0]; 
int w = F[1][0]*M[0][1] + F[1][1]*M[1][1]; 

F[0][0] = x; 
F[0][1] = y; 
F[1][0] = z; 
F[1][1] = w; 
} 

/*function that calculates F[][] raise to the power n and puts the 
result in F[][]*/
static void power(int F[][], int n) 
{ 
int i; 
int M[][] = new int[][]{{1,1},{1,0}}; 

// n - 1 times multiply the matrix to {{1,0},{0,1}} 
for (i = 2; i <= n; i++) 
    multiply(F, M); 
} 

这可以优化为在 O(Logn) 时间复杂度下工作。我们可以在前面的方法中进行递归乘法得到power(M, n)。

static int fib(int n) 
{ 
int F[][] = new int[][]{{1,1},{1,0}}; 
if (n == 0) 
    return 0; 
power(F, n-1); 

return F[0][0]; 
} 

static void multiply(int F[][], int M[][]) 
{ 
int x =  F[0][0]*M[0][0] + F[0][1]*M[1][0]; 
int y =  F[0][0]*M[0][1] + F[0][1]*M[1][1]; 
int z =  F[1][0]*M[0][0] + F[1][1]*M[1][0]; 
int w =  F[1][0]*M[0][1] + F[1][1]*M[1][1]; 

F[0][0] = x; 
F[0][1] = y; 
F[1][0] = z; 
F[1][1] = w; 
} 

static void power(int F[][], int n) 
{ 
if( n == 0 || n == 1) 
  return; 
int M[][] = new int[][]{{1,1},{1,0}}; 

power(F, n/2); 
multiply(F, F); 

if (n%2 != 0) 
   multiply(F, M); 
} 

方法3（O(log n)时间） 下面是一个更有趣的递归公式，可用于在 O(log n) 时间内找到第 n 个斐波那契数。

如果 n 为偶数，则 k = n/2： F(n) = [2*F(k-1) + F(k)]*F(k)

如果 n 是奇数，则 k = (n + 1)/2 F(n) = F(k)*F(k) + F(k-1)*F(k-1) 这个公式是如何工作的？ 该公式可以从上述矩阵方程推导出来。 斐波那契矩阵

取两边的行列式，我们得到 (-1)n = Fn+1Fn-1 – Fn2 此外，由于任意方阵 A 的 AnAm = An+m，因此可以推导出以下恒等式（它们是从矩阵乘积的两个不同系数中获得的）

FmFn + Fm-1Fn-1 = Fm+n-1

通过设置 n = n+1，

FmFn+1 + Fm-1Fn = Fm+n

设 m = n

F2n-1 = Fn2 + Fn-12

F2n = (Fn-1 + Fn+1)Fn = (2Fn-1 + Fn)Fn（来源：Wiki）

为了得到要证明的公式，我们只需要执行以下操作 如果 n 是偶数，我们可以设 k = n/2 如果 n 是奇数，我们可以设 k = (n+1)/2

public static int fib(int n) 
{ 

    if (n == 0) 
        return 0; 

    if (n == 1 || n == 2) 
        return (f[n] = 1); 

    // If fib(n) is already computed 
    if (f[n] != 0) 
        return f[n]; 

    int k = (n & 1) == 1? (n + 1) / 2 
                        : n / 2; 

    // Applyting above formula [See value 
    // n&1 is 1 if n is odd, else 0. 
    f[n] = (n & 1) == 1? (fib(k) * fib(k) +  
                    fib(k - 1) * fib(k - 1)) 
                   : (2 * fib(k - 1) + fib(k))  
                   * fib(k); 

    return f[n]; 
} 

方法 4 - 使用公式 在这种方法中，我们直接实现了斐波那契数列中第 n 项的公式。时间 O(1) 空间 O(1) Fn = {[(√5 + 1)/2] ^ n} / √5

static int fib(int n) { 
double phi = (1 + Math.sqrt(5)) / 2; 
return (int) Math.round(Math.pow(phi, n)  
                    / Math.sqrt(5)); 
} 

参考：http://www.maths.surrey.ac.uk/hosted-sites/R.Knott/Fibonacci/fibFormula.html

Answer 16

我们首先应该注意，斐波那契数 (F(n)) 与 n 增长得非常快，并且对于大于 93 的 n 不能以 64 位 表示。因此，计算它们的程序这样的n 需要使用额外的机制来对这些大数进行操作。现在，仅考虑（大量）操作的计数，顺序计算它们的算法将需要线性操作数。

我们可以从以下关于斐波那契数的恒等式中受益：

F(2m) = 2*F(m)*F(m+1) − (F(m))^2

F(2m+1) = (F(m))^2 + (F(m+1))^2

（像 A^2 这样的符号表示 A 的平方）。

所以，如果我们知道F(m)和F(m+1)，我们可以直接计算F(2m)和F(2m+1)。

考虑n 的二进制表示。请注意，从x = 1 开始，我们可以通过迭代加倍并可能在x 上加1 来生成x = n。这可以通过遍历n 的位来完成，并检查它是 0 还是 1。

这个想法是，我们可以保持F(x) 与x 同步。在每次这样的迭代中，当我们将x 加倍并可能在x 上加1 时，我们还可以使用F(x) 和F(x+1) 的早期值，通过上述等式计算F(x) 的新值。

由于在n 中迭代次数将是对数，所以在n 中总（大数）操作也是对数。