找回密码
 欢迎注册
楼主: happysxyf

[求助] 如何高效计算斐波那契数

[复制链接]
发表于 2016-5-15 11:06:31 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?欢迎注册

×
自己写的斐波那契数计算,只能计算到第50万项,以及前50万项的和,cmd下运行:
  1. fib 500000
复制代码

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

fib.c
  1. #include    <math.h>
  2. #include   <stdio.h>
  3. #include  <stdlib.h>

  4. /*斐波那契函数*/
  5. void fib(int n, int k)
  6. {
  7.         int  Fa[110000];
  8.         int  Fb[110000];
  9.         Fa[1]=Fb[1]=1;
  10.         int  i,j,tmp=0,add;
  11.         if(n==0 &&k==0){printf("fib(0) = 0");}
  12.         else if(n==0 &&k==1){printf("∑fib(0) = 0");}
  13.         else if(n<=2 &&k==0){printf("fib(%d) = 1",n);}
  14.         else if(n<=2 &&k==1){printf("∑fib(%d) = %d",n,n);}
  15.         else{
  16.                 if(k==1){n+=2;}
  17.                 int N=ceil(n*log10((1+sqrt(5))/2)-log10(5)/2);
  18.                 for(i=3;i<=n;i++){
  19.                         add=0;
  20.                         for(j=1;j<=N/9+1;j++){
  21.                                 tmp=Fa[j]+Fb[j]+add;
  22.                                 Fa[j]=Fb[j];
  23.                                 add=tmp/1000000000;
  24.                                 Fb[j]=tmp%1000000000;
  25.                         }
  26.                 }
  27.                 if(n<45){
  28.                         if(k==1){printf("∑fib(%d) = ",n-2);}else{printf("fib(%d) = ",n);}
  29.                         if(k==1){printf("%d",Fb[1]-1);}else{printf("%d",Fb[1]);}
  30.                 }
  31.                 else if(n>=45){
  32.                         if(k==1){printf("∑fib(%d) = ",n-2);}else{printf("fib(%d) = ",n);}
  33.                         printf("%d",Fb[N/9+1]);
  34.                         for(j=N/9;j>1;j--){printf("%09d",Fb[j]);}
  35.                         if(k==1){printf("%09d",Fb[1]-1);}else{printf("%09d",Fb[1]);}
  36.                 }
  37.         }
  38. }

  39. /*命令行下运行*/
  40. int main(int argc, char *argv[])
  41. {
  42.         fib(atoi(argv[1]),0); //求斐波那契数
  43.         printf("\n");
  44.         fib(atoi(argv[1]),1); //求斐波那契数前n项和
  45.         return 0;
  46. }
复制代码

评分

参与人数 1金币 +20 收起 理由
gxqcn + 20 首贴奖励,欢迎常来。

查看全部评分

毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2016-5-15 12:46:03 | 显示全部楼层
本帖最后由 zeroieme 于 2016-5-15 12:47 编辑

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

点评

$sqrt(5)$符号化与二进制法求幂就可以整数计算。  发表于 2016-5-15 22:52
计算机精度不是无限的,这个公式会有浮点误差。  发表于 2016-5-15 13:21
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2016-5-15 14:12:44 | 显示全部楼层
本帖最后由 kastin 于 2016-5-15 14:15 编辑

曾经用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 的二进制序列的迭代方法,也就是上面提到的矩阵法中的二进制法.

点评

只是我不知道怎么获取它dll的接口函数 用的是mathlink  发表于 2016-5-28 20:43
gmp不是一个公共库吗?如果mathematica使用的是原始的,估计接口就是gmp提供的。如果Mathematica做了修改和改进,那么这个就很难知道了,得问他们公司的开发人员。  发表于 2016-5-16 13:36
我把mathematic软件肢解了,提取了它的gmp.dll,原来是靠gmp库在撑腰,只是我不知道怎么获取它dll的接口函数  发表于 2016-5-16 11:28
我刚下载Mahtematica.试了一下,确实无所不能,计算50万项她都不用时间,很逆天,就是计算阶乘不及HugeCalc,HugeCalc也太精于算法,5000万都能啃下,primes 40亿内基本不足1秒,虽然素范围小了点,各有千秋!  发表于 2016-5-15 20:50
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2016-5-15 15:55:51 | 显示全部楼层
本帖最后由 happysxyf 于 2016-5-15 16:22 编辑
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次幂的大数算法,不现实。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2016-5-15 16:07:42 | 显示全部楼层
kastin 发表于 2016-5-15 14:12
曾经用Matlab编写过计算斐波那契数,方法一般有三种:
1)递归法(最耗时间)
2)迭代法(不方便“并行” ...

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

点评

肯定在大数乘法里面做了优化啊,这是基础运算操作,通过内核固化了的,所以它不提。  发表于 2016-5-15 19:32
我没让它显示,显示很长,所以就屏蔽了。  发表于 2016-5-15 19:30
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2016-5-15 16:10:07 | 显示全部楼层
本帖最后由 happysxyf 于 2016-5-15 16:27 编辑


细看我的代码就是那个方法,我很早就推出了和=a(n+2)-1,还推出了一个项数定律,归根结底还是求an慢。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2016-5-15 20:37:57 | 显示全部楼层


pascal,我用freepascal编译了,输了三个数,但是pascal子过程很慢.溢出为负.
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2016-5-15 21:07:35 | 显示全部楼层
kastin 发表于 2016-5-15 14:12
曾经用Matlab编写过计算斐波那契数,方法一般有三种:
1)递归法(最耗时间)
2)迭代法(不方便“并行” ...

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

点评

有可能。  发表于 2016-5-16 10:17
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2016-5-28 23:18:04 | 显示全部楼层
本帖最后由 chyanog 于 2016-5-28 23:20 编辑

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
复制代码

2016-05-28_231718.png
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2016-5-29 10:10:10 | 显示全部楼层


这个幂迭代思路不错,但是C语言实现起来还是要构造大数运算。
最后发现GMP大数运算库自带“斐波那契”函数用的也该方法,不过gmp库太大了。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

小黑屋|手机版|数学研发网 ( 苏ICP备07505100号 )

GMT+8, 2024-3-29 06:00 , Processed in 0.052559 second(s), 20 queries .

Powered by Discuz! X3.5

© 2001-2024 Discuz! Team.

快速回复 返回顶部 返回列表