是否有一种算法可以在亚线性时间内计算第n个斐波那契数?
是否有一种算法可以在亚线性时间内计算第n个斐波那契数?
使用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
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))
我们首先应该注意到,斐波那契数列(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
乘以2并可能加1来使得 x=n
。这可以通过迭代 n
的位,并检查它是0还是1来完成。
思路是我们可以保持 F(x)
与 x
同步。在每次这样的迭代中,当我们将 x
倍增并可能将1添加到 x
时,我们还可以使用上述方程式使用先前的 F(x)
和 F(x+1)
计算出 F(x)
的新值。
由于迭代次数将是 n
的对数,因此总的(大数)操作也是 n
的对数。
我已经了解了一些计算斐波那契数列的高效时间复杂度方法,以下是其中一些:
方法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;
}
方法二-(使用矩阵{{1,1},{1,0}}的幂)
这是一个O(n)算法,它依赖于以下事实:如果我们将矩阵M = {{1,1},{1,0}}乘以自身n次(换句话说,计算power(M,n)),那么我们得到的结果矩阵中行和列(0,0)的元素就是第(n+1)个斐波那契数。这个解决方案的时间复杂度为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);
}
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);
}
方法三(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) 这个公式是如何工作的呢? 该公式可以从上述矩阵方程中推导出来。 fibonaccimatrix
两边取行列式,我们得到 (-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];
}
方法四 - 使用公式 在这种方法中,我们直接实现斐波那契数列中第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
你可以使用奇怪的平方根方程式得到一个精确的答案。原因是$\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
c++14示例实现:
uint64_t fib(unsigned n) {
static struct FIB_NUMBERS {
uint64_t numbers[100] {0, 1};
unsigned max_n = 2;
constexpr FIB_NUMBERS() {
for (; max_n<100; ++max_n) {
numbers[max_n] = numbers[max_n-1] + numbers[max_n-2];
if (numbers[max_n] < numbers[max_n-1]) {
// overflow
numbers[max_n--] = 0;
break;
}
}
}
} fib;
if (n <= fib.max_n)
return fib.numbers[n];
return 0;
}
编译器将生成所有斐波那契数,直到n = 93 运行时只是一个查找。 为了检测来自调用者的溢出:
auto res = fib(n);
if (res==0 && n>0) overflow_msg();