是否有任何算法来计算子线性时间内的第n个斐波纳契数?
n
th斐波纳契数由下式给出
f(n) = Floor(phi^n / sqrt(5) + 1/2)
哪里
phi = (1 + sqrt(5)) / 2
假设原始数学运算(+
,-
,*
和/
)是O(1)
,你可以使用这个结果计算n
时间内的O(log n)
th斐波纳契数(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);
}
除了通过数学方法进行微调之外,最好的解决方案之一(我相信)是使用字典以避免重复计算。
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))
我们从平凡的字典(Fibonacci序列的前两个值)开始,并不断将Fibonacci值添加到字典中。
首个100000斐波纳契值(Intel Xeon CPU E5-2680 @ 2.70 GHz,16 GB RAM,Windows 10-64位操作系统)花了大约0.7秒
看看分而治之的算法here
该链接具有用于该问题的一些其他答案中提到的矩阵求幂的伪代码。
定点算术不准确。 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
这是在O(log n)算术运算中使用大小为O(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。在这个计算中唯一棘手的步骤是用(1 + phi)代替phi ^ 2,因为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个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
表示phi,可以执行(a+b*phi)(c+d*phi)
作为整数运算(a+bX)(c+dX) modulo (X^2-X-1)
。然后pow
函数可以被标准Python pow
函数替换(它方便地包括第三个参数z
,它计算结果模z
。选择的X
是2<<i
。
根据Pillsy对矩阵求幂的引用,对于矩阵
M = [1 1] [1 0]
然后
fib(n) = Mn1,2
使用重复乘法将矩阵提升到幂并不是非常有效。
矩阵求幂的两种方法是分而治之,其以O(ln n)步骤产生Mn,或者是恒定时间的特征值分解,但是由于有限的浮点精度而可能引入误差。
如果希望精确值大于浮点实现的精度,则必须使用基于此关系的O(ln n)方法:
Mn = (Mn/2)2 if n even = M·Mn-1 if n is odd
M上的特征值分解找到两个矩阵U和Λ,使得Λ是对角线和
M = U Λ U-1 Mn = ( U Λ U-1) n = U Λ U-1 U Λ U-1 U Λ U-1 ... n times = U Λ Λ Λ ... U-1 = U Λ n U-1Raising a the diagonal matrix Λ to the nth power is a simple matter of raising each element in Λ to the nth, so this gives an O(1) method of raising M to the nth power. However, the values in Λ are not likely to be integers, so some error will occur.
为我们的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
如果您已经阅读了杰森的答案,您可以看到这将会发生什么。
求解特征向量X1和X2:
if 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 Φ ] 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 ] = 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 ] as ΨΦ = -1 = ( √5 )-1 [ Φn+1-Ψn+1 Φn-Ψn ] [ Φn-Ψn Φn-1-Ψn-1 ]
所以
fib(n) = Mn1,2 = ( Φn - (1-Φ)n ) / √5
这与其他地方给出的公式一致。
您可以从重复关系推导出它,但在工程计算和模拟中计算大矩阵的特征值和特征向量是一项重要的活动,因为它给出了方程组的稳定性和谐波,并允许有效地将矩阵提高到高功率。
如果你想要确切的数字(这是一个“bignum”,而不是一个int / float),那我恐怕
不可能!
如上所述,斐波纳契数的公式为:
fib n = floor(phin /√5+ 1/2)
fib n~ = filter /√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)的时间内计算。
如果您只想要答案的低位数,则可以使用矩阵求幂方法在子线性时间内计算。
其中一个exercises in SICP是关于这个,其答案描述here.
在命令式的风格,程序看起来像
Function Fib(count) a ← 1 b ← 0 p ← 0 q ← 1 While count > 0 Do If Even(count) Then p ← p² + q² q ← 2pq + q² count ← count ÷ 2 Else a ← bq + aq + ap b ← bp + aq count ← count - 1 End If End While Return b End Function
您也可以通过对整数矩阵进行取幂来实现。如果你有矩阵
/ 1 1 \
M = | |
\ 1 0 /
那么(M^n)[1, 2]
将等于n
th斐波那契数,如果[]
是矩阵下标而^
是矩阵求幂。对于固定大小的矩阵,可以在O(log n)时间内以与实数相同的方式对正整数幂进行取幂。
编辑:当然,根据您想要的答案类型,您可以使用恒定时间算法。与其他公式显示的一样,n
th斐波纳契数与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个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)
计算fib(n-1),fib(n-2)
,如果n是奇数,如果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