是否有任何算法可以在亚线性时间内计算第 n 个斐波那契数?
16 回答
根据 Pillsy 对矩阵求幂的引用,对于矩阵
M = [1 1] [1 0]
然后
fib ( n ) = M n 1,2
使用重复乘法来提高矩阵的幂不是很有效。
矩阵求幂的两种方法是分治法,在O ( ln n ) 步中产生M n,或者是常数时间的特征值分解,但由于有限的浮点精度可能会引入错误。
如果您想要一个大于浮点实现精度的精确值,则必须使用基于此关系的 O ( ln n ) 方法:
M n = ( M n /2 ) 2如果n偶数 = M · M n -1如果n是奇数
M上的特征值分解找到两个矩阵U和Λ,使得Λ是对角线并且
M = U Λ U -1 M n = ( U Λ U -1 ) n = U Λ U -1 U Λ U -1 U Λ U -1 ... n 次 = U Λ Λ Λ ... U -1 = U Λ n U -1将对角矩阵Λ提高到n次方是将 Λ中的每个元素提高到n次的简单问题,因此这给出了将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 , λ 2 } = { Φ, 1-Φ } 其中 Φ = ( 1 + √5 ) / 2
如果您已阅读 Jason 的答案,您可以看到它的发展方向。
求解特征向量X 1和X 2:
如果X 1 = [ X 1,1 , X 1,2 ] 米。X 1 1 = λ 1 X 1 X 1,1 + X 1,2 = λ 1 X 1,1 X 1,1 = λ 1 X 1,2 => X 1 = [ Φ, 1 ] X 2 = [ 1-Φ, 1 ]
这些向量给出U:
U = [ X 1,1 , X 2,2 ] [ X 1,1 , X 2,2 ] = [Φ, 1-Φ] [ 1, 1 ]
使用反转U
A = [ ab ] [光盘] => A -1 = ( 1 / | A | ) [ d -b ] [-ca]
所以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 n 1,2所需的一切:
M n = U Λ n U -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 ) = M n 1,2 = ( Φ n - (1-Φ) n ) / √5
这与其他地方给出的公式一致。
您可以从递归关系中推导出它,但在工程计算和模拟中,计算大型矩阵的特征值和特征向量是一项重要的活动,因为它提供了方程组的稳定性和谐波,并允许有效地将矩阵提升到高次幂。
n
斐波那契数由下式给出
f(n) = Floor(phi^n / sqrt(5) + 1/2)
在哪里
phi = (1 + sqrt(5)) / 2
假设原始数学运算(、 和+
)-
是您可以使用此结果及时计算第 斐波那契数(因为公式中的幂运算)。*
/
O(1)
n
O(log n)
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);
}
如果您想要确切的数字(这是“bignum”,而不是 int/float),那么恐怕
不可能!
如上所述,斐波那契数的公式是:
fib n = 楼层 (phi n /√5 + 1 / 2 )
fib n ~= phi n /√5
是多少位数fib n
?
numDigits (fib n) = log (fib n) = log (phi n /√5) = log phi n - log √5 = n * log phi - log √5
numDigits (fib n) = n * const + const
它是O ( n )
由于请求的结果为O ( n ),因此无法在小于O ( n ) 的时间内计算出来。
如果您只想要答案的低位数字,则可以使用矩阵求幂方法在亚线性时间内计算。
您也可以通过对整数矩阵求幂来做到这一点。如果你有矩阵
/ 1 1 \
M = | |
\ 1 0 /
then(M^n)[1, 2]
将等于n
第 Fibonacci 数,如果[]
是矩阵下标并且^
是矩阵求幂。对于固定大小的矩阵,正整数幂的取幂可以在 O(log n) 时间内以与实数相同的方式完成。
编辑:当然,根据您想要的答案类型,您可能可以使用恒定时间算法。像其他公式所示,n
第 斐波那契数随 呈指数增长n
。即使使用 64 位无符号整数,您也只需要一个 94 条目的查找表即可覆盖整个范围。
第二次编辑:首先使用特征分解进行矩阵指数与下面的 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 库。
从一组斐波那契数列开始。使用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 个斐波那契数等进行了测试。
这是我递归 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 是奇数,您可以fib(n),fib(n-1)
使用fib(n-1),fib(n-2)
计算,如果 n 是偶数,您可以fib(n),fib(n-1)
使用计算fib(n/2),fib(n/2-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)]
简化右手边导致偶数情况。
使用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
定点算术是不准确的。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
.
这是一个单线计算 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) 的整数是合理的,因为这与答案的大小相当。
为了理解这一点,设 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) modulo (X^2-X-1)
。然后该pow
函数可以替换为标准 Pythonpow
函数(它方便地包括第三个参数z
,用于计算结果模数z
。X
选择的是2<<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))
我们从平凡的字典(斐波那契数列的前两个值)开始,并不断地将斐波那契值添加到字典中。
前 100000 个斐波那契值大约需要 0.7 秒(Intel Xeon CPU E5-2680 @ 2.70 GHz,16 GB RAM,Windows 10-64 位操作系统)
在此处查看分而治之算法
该链接具有该问题的其他一些答案中提到的矩阵求幂的伪代码。
您可以使用奇怪的平方根方程来获得准确的答案。原因是 $\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
我遇到了一些以有效时间复杂度计算斐波那契的方法,其中一些是 -
方法 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 次(换句话说,计算幂(M, n )),那么我们得到第 (n+1) 个斐波那契数作为结果矩阵中行和列 (0, 0) 的元素。该解决方案将有 O(n) 时间。
矩阵表示给出了斐波那契数的以下封闭表达式:fibonaccimatrix
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) 时间复杂度下工作。我们可以在前面的方法中进行递归乘法得到幂(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(来源:维基)
为了得到要证明的公式,我们只需要做以下事情 如果 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
我们首先应该注意,斐波那契数(F(n))
增长得非常快,n
并且不能用64 位表示n
大于 93。因此计算它们的程序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 = n
通过迭代加倍并可能将 1 添加到x
. 这可以通过遍历 的位n
并检查它是 0 还是 1 来完成。
这个想法是,我们可以保持F(x)
同步x
。在每次这样的迭代中,当我们将 加倍x
并可能加 1 时x
,我们还可以F(x)
使用 和 的早期值计算新值,F(x)
并F(x+1)
使用上述等式。
Since the number of iterations will be logarithmic in n
, the total (large-number) operations are also logarithmic in n
.