(22:31) gp > m=
%3 =
(22:35) gp > Mod(m,1234567891011)
%4 =
(22:35) gp > %4^(10^14-1)
%5 =
(22:35) gp > %5.~
*** unknown member function: %5.~
^------
(22:36) gp > %5*~
%6 = ~
(22:36) gp >
所以结果是921144120792
11# mathe
:b:
这么快！！
我还以为要运行10^14次！
不亲自用Gp运行 我还不相信呢！
:L
大家都给出算法思想了，我的反应的确有点慢了，呵呵！
特别谢谢mathe 3次提醒，:loveliness:int M = {{1,0}{0,1}}
int fib(int n)
{
matpow(n-1);
return M;
}
void matpow(int n)
{
if (n > 1) {
matpow(n/2);
M = M*M;
}
if (n is odd) M = M*{{1,1}{1,0}}
}
引入了矩阵 之所以快主要是因为二分式的模幂运算 避免了大量无谓的计算。
其实写成$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)$
方法很妙，但最后还有个1/sqrt5如何处理？
在32位系统下，可将1234567891011分解成(3*7*13*67*107)*630803 = 630803 * 1957137，
分别用 630803、1957137 为模计算，然后再用孙子定理，效率应该更高。
15# mathe
我也正在想能不能针对fibonacci的特性，对模逆算法做一些定制，
mathe刚好正在做这一步，。。。
太精妙了！
可以用建一个数据结构，避免符号运算，另外
在模不是5的倍数下，我们知道ω是方程x2-x-1=0的根
这个在任何情况下都成立的吧
(22:31) gp > m=
%3 =
(22:35) gp > Mod(m,1234567891011)
%4 =
(22:35) gp...
mathe 发表于 2012-2-20 22:37 https://bbs.emath.ac.cn/static/image/common/back.gif
完全能看懂，第一次知道pari/gp可以计算矩阵的模幂，
我以前只是用它计算整数的模幂，跟mathe学到了以前
不知道的知识
这个其实不就是lucas序列吗？其实郭的素性判定算法中有，我在第一个回复的时候就说了，只是郭愿不愿意显露他的身手而已！