fibonacci(10^14) % 1234567891011是多少?

mathe

(22:31) gp > m=[1,1;1,0] %3 = [1 1] [1 0] (22:35) gp > Mod(m,1234567891011) %4 = [Mod(1, 1234567891011) Mod(1, 1234567891011)] [Mod(1, 1234567891011) Mod(0, 1234567891011)] (22:35) gp > %4^(10^14-1) %5 = [Mod(921144120792, 1234567891011) Mod(512884989226, 1234567891011)] [Mod(512884989226, 1234567891011) Mod(408259131566, 1234567891011)] (22:35) gp > %5.[1,1]~ *** unknown member function: %5.[1,1]~ ^------ (22:36) gp > %5*[1,1]~ %6 = [Mod(199461219007, 1234567891011), Mod(921144120792, 1234567891011)]~ (22:36) gp > 所以结果是921144120792

wayne

11# mathe 这么快！！ 我还以为要运行10^14次！ 不亲自用Gp运行 我还不相信呢！

wayne

大家都给出算法思想了，我的反应的确有点慢了，呵呵！ 特别谢谢mathe 3次提醒，

1. int M[2][2] = {{1,0}{0,1}}
3. int fib(int n)
4. {
5. matpow(n-1);
6. return M[0][0];
7. }
9. void matpow(int n)
10. {
11. if (n > 1) {
12. matpow(n/2);
13. M = M*M;
14. }
15. if (n is odd) M = M*{{1,1}{1,0}}
16. }

复制代码

wayne

引入了矩阵 之所以快 主要是因为二分式的模幂运算 避免了大量无谓的计算。

mathe

其实写成$1/{sqrt{5}}(\omega^n+(-\omega)^{-n})$也可以用模幂运算 只是这时需要符号运算比较好。在模不是5的倍数下，我们知道$\omega$是方程$x^2-x-1=0$的根 也就是说，我们同样可以计算$\omega^n(mod m)$,只是这时，我们认为$\omega^2 -=\omega+1(mod m)$ 于是任何时候，$\omega^k -=u_k+v_k*\omega (mod m)$ 于是得出$\omega^{2k} -= (u_k+v_k*\omega)^2(mod m) -=u_k^2+2u_kv_k*\omega+v_k^2(\omega+1) -=(u_k^2+v_k^2)+(2u_kv_k+v_k^2)\omega(mod m)$

gxqcn

方法很妙，但最后还有个$1/sqrt5$如何处理？

gxqcn

在32位系统下，可将1234567891011分解成(3*7*13*67*107)*630803 = 630803 * 1957137， 分别用 630803、1957137 为模计算，然后再用孙子定理，效率应该更高。

wayne

15# mathe 我也正在想能不能针对fibonacci的特性，对模逆算法做一些定制， mathe刚好正在做这一步，。。。 太精妙了！ 可以用建一个数据结构，避免符号运算，另外

在模不是5的倍数下，我们知道ω是方程x2-x-1=0的根

这个在任何情况下都成立的吧

mathematica

(22:31) gp > m=[1,1;1,0] %3 = [1 1] [1 0] (22:35) gp > Mod(m,1234567891011) %4 = [Mod(1, 1234567891011) Mod(1, 1234567891011)] [Mod(1, 1234567891011) Mod(0, 1234567891011)] (22:35) gp ... mathe 发表于 2012-2-20 22:37

完全能看懂，第一次知道pari/gp可以计算矩阵的模幂， 我以前只是用它计算整数的模幂，跟mathe学到了以前 不知道的知识

mathematica

这个其实不就是lucas序列吗？其实郭的素性判定算法中有，我在第一个回复的时候就说了，只是郭愿不愿意显露他的身手而已！