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

mathematica

把我的代码的最后一行改进一下:

1. (*问题:fibonacci(10^14) % 1234567891011是多少?*)
2. (*问题网址:http://bbs.emath.ac.cn/thread-4054-1-1.html*)
3. (*参考资料:http://en.wikipedia.org/wiki/Modular_exponentiation*)
4. Clear["Global`*"];(*Clear all variables*)
5. result={{1,0},{0,1}};(*单位矩阵*)(*也是最终的求解的模幂的结果矩阵*)
6. base={{1,1},{1,0}};(*模幂的底,此处是矩阵*)
7. modulus=1234567891011;(*模*)
8. n=10^14;
9. exponent=n-1;(*模幂的指数*)
10. (*循环计算,此处是矩阵的乘法运算符号 . ,而不是 * 这个运算符号*)
11. While[
12. exponent>0,
13. If[ Mod[exponent,2]==1,
14. (*如果是奇数按照下列方式处理,以及处理指数*)
15. result=Mod[result.base,modulus];exponent=exponent-1,
16. (*如果是偶数,按照下列方式处理,以及处理指数*)
17. base=Mod[base.base,modulus];exponent=exponent/2
18. ]
19. ];
20. out=result.{{1},{1}};(*模幂的最后的矩阵再乘以列向量*)
24. Mod[out[[2,1]],modulus](*第二行第一列是所要求得到的结果*)

复制代码

mathematica

最最正确的代码!!!如下所示

1. (*问题:fibonacci(10^14) % 1234567891011是多少?*)
2. (*问题网址:http://bbs.emath.ac.cn/thread-4054-1-1.html*)
3. (*参考资料:http://en.wikipedia.org/wiki/Modular_exponentiation*)
4. Clear["Global`*"];(*Clear all variables*)
5. result={{1,0},{0,1}};(*单位矩阵*)(*也是最终的求解的模幂的结果矩阵*)
6. base={{1,1},{1,0}};(*模幂的底,此处是矩阵*)
7. modulus=1234567891011;(*模*)
8. n=10^14;
9. exponent=n-1;(*模幂的指数*)
10. (*循环计算,此处是矩阵的乘法运算符号 . ,而不是 * 这个运算符号*)
11. While[
12. exponent>0,
13. If[ Mod[exponent,2]==1,
14. (*如果是奇数按照下列方式处理,以及处理指数*)
15. result=Mod[result.base,modulus];exponent=exponent-1,
16. (*如果是偶数,按照下列方式处理,以及处理指数*)
17. base=Mod[base.base,modulus];exponent=exponent/2
18. ]
19. ];
20. out=result.{{1},{1}};(*模幂的最后的矩阵再乘以列向量*)
21. Mod[out[[2,1]],modulus](*第二行第一列是所要求得到的结果*)

复制代码

mathematica

为什么只有精华帖,却没有精华回复呢? 我觉得我的回复还算是比较精华的吧????????????

mathematica

调用http://bbs.emath.ac.cn/viewthread.php?tid=4448&extra= 处的子函数的求解

1. base={{1,1},{1,0}};(*模幂的底,此处是矩阵*)
2. modulus=1234567891011;(*模*)
3. n=10^14;
4. exponent=n-1;(*模幂的指数*)
5. out=myPowerMod[base,exponent,modulus].{{1},{1}};(*模幂的最后的矩阵再乘以列向量*)
6. Mod[out[[2,1]],modulus](*第二行第一列是所要求得到的结果*)

复制代码

结果是921144120792

nyy

我直接用lucas数列的递推关系来求解这个问题。

1. Clear["Global`*"];(*Clear all variables*)
   
2. n=1234567891011;(*模*)
   
3. mi2=IntegerDigits[10^14,2];(*把m写成二进制的方式*)
   
4. UU=1;VV=1;(*分别是lucas序列的U(1)与V(1)=P的值*)
   
5. P=1;
   
6. QQ=Q=-1;(*这个其实是Q的值*)
   
7. d=(-1)^2-4*(-1);(*判别式的计算*)
   
8. Do[
   
9.     UU=Mod[UU*VV,n];
   
10.     VV=Mod[VV^2-2*QQ,n];
    
11.     QQ=Mod[QQ^2,n];(*这行必须加上*)
    
12.     If[mi2[[kk]]==1,
    
13.         UUtemp=(P*UU+VV);
    
14.         VVtemp=(d*UU+P*VV);
    
15.         (*UUtemp,VVtemp都可能是奇数,如果是奇数,则加上一个n,这样就是偶数了,
    
16.         下面才能除以2得到整数,至于为什么要这么做,我也不是太清楚为什么,
    
17.         可能是由于n是奇数,模的时候,减去偶数个n奇偶性不变,而减去奇数个n奇偶性变了.
    
18.         20220824补充:除以2,其实是乘以2模n的逆,而2模n的逆=(n+1)/2,这里n肯定是奇数
    
19.         当UUtemp=2k(偶数的时候),2k*(n+1)/2=k*(n+1)=k(mod n),相当于直接除以2,
    
20.         当UUtemp=2k+1(奇数的时候),(2k+1)*(n+1)/2=2k*(n+1)/2+(n+1)/2=k+(n+1)/2=(2k+1+n)/2(mod n),相当于加上n后再除以2,
    
21.         VVtemp同理
    
22.         *)
    
23.         If[OddQ[UUtemp],UUtemp=UUtemp+n];
    
24.         If[OddQ[VVtemp],VVtemp=VVtemp+n];
    
25.         UU=Mod[UUtemp/2,n];
    
26.         VV=Mod[VVtemp/2,n];
    
27.         QQ=Mod[QQ*Q,n](*这个必须加上*)
    
28.     ],
    
29.     {kk,2,Length@mi2}];(*此处必须从2开始,程序没有错误!*)
    
30. Print[{UU,VV}]
    

复制代码

求解结果
{921144120792,712346208233}第一个是求解结果

相关数学链接见
https://handwiki.org/wiki/Lucas_sequence

nyy

mathe 发表于 2012-2-20 16:26
这里用这个公式用处不大，要用gxqcn的方法
${(F_{n+1}=F_n+F_{n-1}),(F_n=F_n):}$
所以得到

用https://bbs.emath.ac.cn/forum.ph ... =4054&pid=41742的办法

\[[(F_{n}),(F_{n-1})]=[(1,1),(1,0)]^{n-2}[(F_2),(F_1)]=[(1,1),(1,0)]^{n-1}[(1),(1)]\]

1. Clear["Global`*"];(*Clear all variables*)
   
2. mat={{1,1},{1,0}};(*递推矩阵*)
   
3. n=10^14-2;(*矩阵的n次方*)
   
4. m=1234567891011;(*模*)
   
5. aa=Algebra`MatrixPowerMod[mat,n,m](*计算矩阵的n次方模m的结果*)
   
6. bb=aa.{{1},{1}}(*矩阵乘以向量*)
   
7. cc=Mod[bb[[1,1]],m](*计算模m后的结果*)
   

复制代码

计算结果

{{512884989226, 408259131566}, {408259131566, 104625857660}}

{{921144120792}, {512884989226}}

921144120792 这个就是要的结果

参考资料
https://mathematica.stackexchang ... xpower-with-modulus

nyy

斐波那契数列模n是有周期性的，而m=10^14明显比n大，
因此可以减小m，再考虑质因数分解，再考虑孙子定理，
因此这个问题肯定可以被大幅度的化简，
暂时手机不方便回答问题，下周一看看。
参考资料。
https://blog.csdn.net/weixin_43817941/article/details/128770480

nyy

nyy 发表于 2024-1-7 09:51
斐波那契数列模n是有周期性的，而m=10^14明显比n大，
因此可以减小m，再考虑质因数分解，再考虑孙子定理，
...

直接从通项公式干起!

\[\displaystyle{ F_n = \cfrac{1}{\sqrt{5}}\left(\cfrac{1+\sqrt{5}}{2}\right)^n-\cfrac{1}{\sqrt{5}}\left(\cfrac{1-\sqrt{5}}{2}\right)^n
=\frac{\sqrt{5}((1+\sqrt{5})^n-(1+\sqrt{5})^n)}{5\times2^n}}\]

先算出分母模m的逆，
再算出分子模m的剩余（先算出\((1+\sqrt{5})^n\)模m的剩余，再算出\((1-\sqrt{5})^n\)模m的剩余，前者减去后者，再乘以根号5），
再乘以上面的逆，对这个乘积模m，就能得到答案！

全部代码

1. Clear["Global`*"];(*Clear all variables*)
   
2. (*子函数,根号d的二次域模m的乘法.xx={x1,y1},表示x1+y1*sqrt(d),yy同,d,m为整数*)
   
3. dmult[xx_,yy_,d_,m_]:=Module[{x1=xx[[1]],y1=xx[[2]],x2=yy[[1]],y2=yy[[2]]},Mod[{x1*x2+d*y1*y2,x2*y1+x1*y2},m]]
   
4. (*子函数,根号d的二次域的模m的模幂算法.base={x1,y1},表示x1+y1*sqrt(d),其余参数为整数*)
   
5. dpowermod[base_,d_,n_,m_]:=Module[{result,n2},
   
6.     n2=IntegerDigits[n,2];(*n次方,写成二进制形式*)
   
7.     result=base;
   
8.     Do[result=dmult[result,result,d,m];(*平方一下*)
   
9.        If[n2[[k]]==1,result=dmult[result,base,d,m]](*如果是1,那么就再乘以base一下*)
   
10.     ,{k,2,Length[n2]}];
    
11.     Return[result](*返回最终结果*)
    
12. ]
    
13. (*初始赋值*)
    
14. n=10^14;(*幂*)
    
15. m=1234567891011;(*模*)
    
16. aa=PowerMod[5,-1,m]*PowerMod[2,-n,m](*5*2^n模m的逆,5*2^n位于通项公式的分母上*)
    
17. bb=dpowermod[{1, 1},5,n,m](*通项公式的第1部分的(1+根5)^n mod m*)
    
18. cc=dpowermod[{1,-1},5,n,m](*通项公式的第2部分的(1-根5)^n mod m*)
    
19. dd=dmult[{0,1},bb-cc,5,m](*第1部分减去第2部分,再乘以根号5*)
    
20. ee=Mod[aa*dd[[1]],m](*最后结果*)
    

复制代码

求解结果
673708543183564435712198

{229821218902, 838772581899}

{229821218902, 395795309112}

{980318472924, 0}

921144120792这个就是最终需要求的结果

nyy

nyy 发表于 2024-1-5 16:08
用https://bbs.emath.ac.cn/forum.ph ... =4054&pid=41742的办法

\[[(F_{n ...

还可以减少几行代码

1. Clear["Global`*"];(*Clear all variables*)
   
2. mat={{1,1},{1,0}};(*递推矩阵*)
   
3. n=10^14;(*矩阵的n次方*)
   
4. m=1234567891011;(*模*)
   
5. aa=Algebra`MatrixPowerMod[mat,n,m](*计算矩阵的n次方模m的结果*)
   

复制代码

求解结果如下
\[\left(
\begin{array}{cc}
199461219007 & 921144120792 \\
921144120792 & 512884989226 \\
\end{array}
\right)\]

从上面可以知道，副对角线上的元素就是要求解的结果！

利用了结果
\[\displaystyle{ \begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix}^n = \begin{pmatrix} F_{n+1} & F_n \\ F_n & F_{n-1} \end{pmatrix}. }\]

参考资料
https://handwiki.org/wiki/Fibonacci%20number

nyy

mathe 发表于 2012-2-21 08:48
其实写成$1/{sqrt{5}}(\omega^n+(-\omega)^{-n})$也可以用模幂运算
只是这时需要符号运算比较好。在模不是5 ...

用你的办法！

pari gp的代码

1. m=1234567891011 /*模*/
   
2. r1=Mod(Mod(1,m)*x,x^2-x-1)^(10^14) /*10^14次方*/
   
3. print(lift(r1))
   

复制代码

输出结果

1.    [logfile is "pari.log"]
   
2. (12:08) gp > m=1234567891011
   
3. %1 = 1234567891011
   
4. (12:09) gp > r1=Mod(Mod(1,m)*x,x^2-x-1)^(10^14)
   
5. %2 = Mod(Mod(921144120792, 1234567891011)*x + Mod(512884989226, 1234567891011), x^2 - x - 1)
   
6. (12:09) gp > print(lift(r1))
   
7. Mod(921144120792, 1234567891011)*x + Mod(512884989226, 1234567891011)
   

复制代码

根据这个结果，得到
于是x^(10^14)=921144120792*x+512884989226 (mod x^2-x-1) 多项式系数已经对m=1234567891011取模

再把x分别用(1 +Sqrt[5])/2、(1 - Sqrt[5])/2替换一下，得到

1. Clear["Global`*"];(*Clear all variables*)
   
2. f=921144120792*x+512884989226
   
3. aa=(f/.x->((1+Sqrt[5])/2))-(f/.x->((1-Sqrt[5])/2))//Simplify
   
4. 
   

复制代码

得到
\[921144120792 \sqrt{5}\]
再除以根号5，
得到921144120792（严格地说，这个数需要模一下才是最终结果！）