如何高效计算斐波那契数

happysxyf

自己写的斐波那契数计算,只能计算到第50万项,以及前50万项的和,cmd下运行:

1. fib 500000

复制代码

一共用时17秒,感觉太慢了,有啥高效的计算斐波那契数的方法吗?
我觉得这个肯定能压到1秒内出结果,但是不知到代码该怎么改进?

fib.c

1. #include    <math.h>
   
2. #include   <stdio.h>
   
3. #include  <stdlib.h>
   
4. 
   
5. /*斐波那契函数*/
   
6. void fib(int n, int k)
   
7. {
   
8.         int  Fa[110000];
   
9.         int  Fb[110000];
   
10.         Fa[1]=Fb[1]=1;
    
11.         int  i,j,tmp=0,add;
    
12.         if(n==0 &&k==0){printf("fib(0) = 0");}
    
13.         else if(n==0 &&k==1){printf("∑fib(0) = 0");}
    
14.         else if(n<=2 &&k==0){printf("fib(%d) = 1",n);}
    
15.         else if(n<=2 &&k==1){printf("∑fib(%d) = %d",n,n);}
    
16.         else{
    
17.                 if(k==1){n+=2;}
    
18.                 int N=ceil(n*log10((1+sqrt(5))/2)-log10(5)/2);
    
19.                 for(i=3;i<=n;i++){
    
20.                         add=0;
    
21.                         for(j=1;j<=N/9+1;j++){
    
22.                                 tmp=Fa[j]+Fb[j]+add;
    
23.                                 Fa[j]=Fb[j];
    
24.                                 add=tmp/1000000000;
    
25.                                 Fb[j]=tmp%1000000000;
    
26.                         }
    
27.                 }
    
28.                 if(n<45){
    
29.                         if(k==1){printf("∑fib(%d) = ",n-2);}else{printf("fib(%d) = ",n);}
    
30.                         if(k==1){printf("%d",Fb[1]-1);}else{printf("%d",Fb[1]);}
    
31.                 }
    
32.                 else if(n>=45){
    
33.                         if(k==1){printf("∑fib(%d) = ",n-2);}else{printf("fib(%d) = ",n);}
    
34.                         printf("%d",Fb[N/9+1]);
    
35.                         for(j=N/9;j>1;j--){printf("%09d",Fb[j]);}
    
36.                         if(k==1){printf("%09d",Fb[1]-1);}else{printf("%09d",Fb[1]);}
    
37.                 }
    
38.         }
    
39. }
    
40. 
    
41. /*命令行下运行*/
    
42. int main(int argc, char *argv[])
    
43. {
    
44.         fib(atoi(argv[1]),0); //求斐波那契数
    
45.         printf("\n");
    
46.         fib(atoi(argv[1]),1); //求斐波那契数前n项和
    
47.         return 0;
    
48. }

复制代码

zeroieme

$a_{n}=\frac{\sqrt{5}}{5} \cdot [(\frac{1 + \sqrt{5}}{2})^{n} - (\frac{1 - \sqrt{5}}{2})^{n}]$
求前n项和也就是两等比数列的和

kastin

曾经用Matlab编写过计算斐波那契数，方法一般有三种：
1）递归法（最耗时间）
2）迭代法（不方便“并行”）
3）矩阵法（方便编程和“并行”）
4）通项公式法（有浮点误差，但比较小）
显然是后面三种方法较快，只是它们之间各有优缺点。

矩阵法这里简单介绍一下。
若用 `F_n` 表示第 `n` 项斐波那契数，定义式可改写为矩阵形式$$\begin{pmatrix}F_{n+1}\\F_n\end{pmatrix}=\begin{pmatrix}1 & 1 \\ 1 & 0\end{pmatrix} \begin{pmatrix}F_n\\F_{n-1}\end{pmatrix}$$可记 `A=\left(\begin{smallmatrix}1 & 1\\1 & 0\end{smallmatrix}\right)`，计算 `B=A^{n-1}`，那么 `F_n=B_{11}`.

当然，为了提速，和“并行”（即计算多个斐波那契数），矩阵的n次方，根据矩阵乘法满足结合律，可以用 `O(\log n)` 的算法求：二分法或者幂的二进制法（二者本质是一样的）。

二分法指的是计算矩阵 `A` 的 `n` 幂 `A^n`：
若指数满足 `2|n`，那么可以得到 `A^n=A^{\frac n2}A^{\frac n2}`;否则 `A^n=A^{\lfloor\frac n2\rfloor}A^{\lfloor\frac n2\rfloor}A`.
反复运用上述算法即可快速得到结果，不必一个个地连乘。

二进制法：
将幂指数 `n` 分解为二进制数 ，比如 `50=(110010)_2=2^5+2^4+2^1`，那么那么 `A^{50}= A^{32}A^{16}A^{2}` .
因此先计算出 `A^2`，然后再用它自乘得到 `A^4`，`A^4`自乘得到 `A^8`，...，得到 `A^16`，自乘得到 `A^32`，于是有了 `A^{32},A^{16},A^2`，结果就出来了。其中一共用到5+2=7次矩阵乘法而已。

用Mahtematica来检验

1. Fibonacci[500000]; // Timing

复制代码

{0.015600,Null}
计算才花费0.0156 s

它是怎么实现的呢？在《关于内部实现的一些注释》中，提到了Fibonacci[n] 使用的是基于 n 的二进制序列的迭代方法，也就是上面提到的矩阵法中的二进制法.

happysxyf

zeroieme 发表于 2016-5-15 12:46
$a_{n}=\frac{\sqrt{5}}{5} \cdot [(\frac{1 + \sqrt{5}}{2})^{n} - (\frac{1 - \sqrt{5}}{2})^{n}]$
求前 ...

多谢,但是求和不是关键,你没发现
a(n+2)=a(n+1)+a(n)
a(n+2)=[a(n)+a(n-1)]+a(n)=a(n)+[a(n-1)+a(n)]
a(n+2)=[a(n-1)+a(n-2)]+[a(n-1)+a(n)]=a(n-1)+[a(n-2)+a(n-1)+a(n)]
... ...
最后
a(n+2)=1+[a(1)+a(2)+...+a(n)]
所以前n项和=a(n+2)-1。所以归根结底还是求a(n)
我的代码就是先展示一下a(n)的值，再展示a(n+2)-1就是前n项和。
单独算a(50万)或前50万的和都是8秒左右。我希望在1秒内算出a(50万)。
a(50万)足有11万位左右。单纯通项公式后边的几百位基本就是无用的值。还得实现n次幂的大数算法，不现实。

happysxyf

kastin 发表于 2016-5-15 14:12
曾经用Matlab编写过计算斐波那契数，方法一般有三种：
1）递归法（最耗时间）
2）迭代法（不方便“并行” ...

快速幂，指数二进制化迭代，这个我也考虑过，但是第50万项太大了有11万多位，里边涉及到的幂乘都是大数，还得处理大数的乘方、乘法运算，不借助GMP肯定慢。{0.015600,Null}，Mahtematica为何会显示null，她显示出来11万位的结果了吗？用时0.0156秒绝对是其乘法乘方库的优化做到的。

happysxyf

math_humanbeing 发表于 2016-5-15 15:27
http://bbs.emath.ac.cn/blog-4014-50.html
也是个办法

细看我的代码就是那个方法，我很早就推出了和=a(n+2)-1,还推出了一个项数定律，归根结底还是求an慢。

happysxyf

math_humanbeing 发表于 2016-5-15 17:07

pascal,我用freepascal编译了,输了三个数,但是pascal子过程很慢.溢出为负.

happysxyf

kastin 发表于 2016-5-15 14:12
曾经用Matlab编写过计算斐波那契数，方法一般有三种：
1）递归法（最耗时间）
2）迭代法（不方便“并行” ...

Mahtematica可能预存了某些基础a(n)值,当Fibonacci[N]很大时,她用基础a(n)直接快速幂迭代了Fibonacci[N],软件可能采用了偷懒的方式完成了这项工作,

chyanog

https://www.felix021.com/blog/read.php?2111

1. def fibonacci(n):
   
2.     def iter(a, b, p, q, n):
   
3.         if n == 0:
   
4.             return b
   
5.         elif n % 2 == 0:
   
6.             return iter(a, b, p * p + q * q, 2 * p * q + q * q, n / 2)
   
7.         else:
   
8.             return iter(a * (p + q) + b * q, a * q + b * p, p, q, n - 1)
   
9.     return iter(1, 0, 0, 1, n)

复制代码

1. {1, 0, 0, 1, 10^7} //. {a_, b_, p_, q_, n_} :>
   
2.     Switch[n, 0,
   
3.      b, _?EvenQ, {a, b, p^2 + q^2, 2 p q + q^2,
   
4.       n/2}, _, {a*(p + q) + b*q, a q + b p, p, q, n - 1}] //
   
5.   IntegerLength // AbsoluteTiming

复制代码

happysxyf

chyanog 发表于 2016-5-28 23:18
https://www.felix021.com/blog/read.php?2111

这个幂迭代思路不错，但是C语言实现起来还是要构造大数运算。
最后发现GMP大数运算库自带“斐波那契”函数用的也该方法，不过gmp库太大了。