把我的代码的最后一行改进一下:(*问题:fibonacci(10^14) % 1234567891011是多少?*)
(*问题网址:http://bbs.emath.ac.cn/thread-4054-1-1.html*)
(*参考资料:http://en.wikipedia.org/wiki/Modular_exponentiation*)
Clear["Global`*"];(*Clear all variables*)
result={{1,0},{0,1}};(*单位矩阵*)(*也是最终的求解的模幂的结果矩阵*)
base={{1,1},{1,0}};(*模幂的底,此处是矩阵*)
modulus=1234567891011;(*模*)
n=10^14;
exponent=n-1;(*模幂的指数*)
(*循环计算,此处是矩阵的乘法运算符号 . ,而不是 * 这个运算符号*)
While[
exponent>0,
If[ Mod==1,
(*如果是奇数按照下列方式处理,以及处理指数*)
result=Mod;exponent=exponent-1,
(*如果是偶数,按照下列方式处理,以及处理指数*)
base=Mod;exponent=exponent/2
]
];
out=result.{{1},{1}};(*模幂的最后的矩阵再乘以列向量*)
Mod],modulus](*第二行第一列是所要求得到的结果*)
最最正确的代码!!!如下所示(*问题:fibonacci(10^14) % 1234567891011是多少?*)
(*问题网址:http://bbs.emath.ac.cn/thread-4054-1-1.html*)
(*参考资料:http://en.wikipedia.org/wiki/Modular_exponentiation*)
Clear["Global`*"];(*Clear all variables*)
result={{1,0},{0,1}};(*单位矩阵*)(*也是最终的求解的模幂的结果矩阵*)
base={{1,1},{1,0}};(*模幂的底,此处是矩阵*)
modulus=1234567891011;(*模*)
n=10^14;
exponent=n-1;(*模幂的指数*)
(*循环计算,此处是矩阵的乘法运算符号 . ,而不是 * 这个运算符号*)
While[
exponent>0,
If[ Mod==1,
(*如果是奇数按照下列方式处理,以及处理指数*)
result=Mod;exponent=exponent-1,
(*如果是偶数,按照下列方式处理,以及处理指数*)
base=Mod;exponent=exponent/2
]
];
out=result.{{1},{1}};(*模幂的最后的矩阵再乘以列向量*)
Mod],modulus](*第二行第一列是所要求得到的结果*)
为什么只有精华帖,却没有精华回复呢?
我觉得我的回复还算是比较精华的吧????????????
调用http://bbs.emath.ac.cn/viewthread.php?tid=4448&extra=
处的子函数的求解base={{1,1},{1,0}};(*模幂的底,此处是矩阵*)
modulus=1234567891011;(*模*)
n=10^14;
exponent=n-1;(*模幂的指数*)
out=myPowerMod.{{1},{1}};(*模幂的最后的矩阵再乘以列向量*)
Mod],modulus](*第二行第一列是所要求得到的结果*)
结果是921144120792
我直接用lucas数列的递推关系来求解这个问题。
Clear["Global`*"];(*Clear all variables*)
n=1234567891011;(*模*)
mi2=IntegerDigits;(*把m写成二进制的方式*)
UU=1;VV=1;(*分别是lucas序列的U(1)与V(1)=P的值*)
P=1;
QQ=Q=-1;(*这个其实是Q的值*)
d=(-1)^2-4*(-1);(*判别式的计算*)
Do[
UU=Mod;
VV=Mod;
QQ=Mod;(*这行必须加上*)
If]==1,
UUtemp=(P*UU+VV);
VVtemp=(d*UU+P*VV);
(*UUtemp,VVtemp都可能是奇数,如果是奇数,则加上一个n,这样就是偶数了,
下面才能除以2得到整数,至于为什么要这么做,我也不是太清楚为什么,
可能是由于n是奇数,模的时候,减去偶数个n奇偶性不变,而减去奇数个n奇偶性变了.
20220824补充:除以2,其实是乘以2模n的逆,而2模n的逆=(n+1)/2,这里n肯定是奇数
当UUtemp=2k(偶数的时候),2k*(n+1)/2=k*(n+1)=k(mod n),相当于直接除以2,
当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,
VVtemp同理
*)
If,UUtemp=UUtemp+n];
If,VVtemp=VVtemp+n];
UU=Mod;
VV=Mod;
QQ=Mod(*这个必须加上*)
],
{kk,2,Length@mi2}];(*此处必须从2开始,程序没有错误!*)
Print[{UU,VV}]
求解结果
{921144120792,712346208233}第一个是求解结果
相关数学链接见
https://handwiki.org/wiki/Lucas_sequence
本帖最后由 nyy 于 2024-1-5 16:10 编辑
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.php?mod=redirect&goto=findpost&ptid=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)]\]
Clear["Global`*"];(*Clear all variables*)
mat={{1,1},{1,0}};(*递推矩阵*)
n=10^14-2;(*矩阵的n次方*)
m=1234567891011;(*模*)
aa=Algebra`MatrixPowerMod(*计算矩阵的n次方模m的结果*)
bb=aa.{{1},{1}}(*矩阵乘以向量*)
cc=Mod],m](*计算模m后的结果*)
计算结果
{{512884989226, 408259131566}, {408259131566, 104625857660}}
{{921144120792}, {512884989226}}
921144120792 这个就是要的结果
参考资料
https://mathematica.stackexchange.com/questions/123236/matrixpower-with-modulus
斐波那契数列模n是有周期性的，而m=10^14明显比n大，
因此可以减小m，再考虑质因数分解，再考虑孙子定理，
因此这个问题肯定可以被大幅度的化简，
暂时手机不方便回答问题，下周一看看。
参考资料。
https://blog.csdn.net/weixin_43817941/article/details/128770480
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，就能得到答案！
全部代码
Clear["Global`*"];(*Clear all variables*)
(*子函数,根号d的二次域模m的乘法.xx={x1,y1},表示x1+y1*sqrt(d),yy同,d,m为整数*)
dmult:=Module[{x1=xx[],y1=xx[],x2=yy[],y2=yy[]},Mod[{x1*x2+d*y1*y2,x2*y1+x1*y2},m]]
(*子函数,根号d的二次域的模m的模幂算法.base={x1,y1},表示x1+y1*sqrt(d),其余参数为整数*)
dpowermod:=Module[{result,n2},
n2=IntegerDigits;(*n次方,写成二进制形式*)
result=base;
Do;(*平方一下*)
If]==1,result=dmult](*如果是1,那么就再乘以base一下*)
,{k,2,Length}];
Return(*返回最终结果*)
]
(*初始赋值*)
n=10^14;(*幂*)
m=1234567891011;(*模*)
aa=PowerMod*PowerMod(*5*2^n模m的逆,5*2^n位于通项公式的分母上*)
bb=dpowermod[{1, 1},5,n,m](*通项公式的第1部分的(1+根5)^n mod m*)
cc=dpowermod[{1,-1},5,n,m](*通项公式的第2部分的(1-根5)^n mod m*)
dd=dmult[{0,1},bb-cc,5,m](*第1部分减去第2部分,再乘以根号5*)
ee=Mod],m](*最后结果*)
求解结果
673708543183564435712198
{229821218902, 838772581899}
{229821218902, 395795309112}
{980318472924, 0}
921144120792这个就是最终需要求的结果
本帖最后由 nyy 于 2024-1-8 09:23 编辑
nyy 发表于 2024-1-5 16:08
用https://bbs.emath.ac.cn/forum.php?mod=redirect&goto=findpost&ptid=4054&pid=41742的办法
\[[(F_{n ...
还可以减少几行代码
Clear["Global`*"];(*Clear all variables*)
mat={{1,1},{1,0}};(*递推矩阵*)
n=10^14;(*矩阵的n次方*)
m=1234567891011;(*模*)
aa=Algebra`MatrixPowerMod(*计算矩阵的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 于 2024-1-8 12:17 编辑
mathe 发表于 2012-2-21 08:48
其实写成$1/{sqrt{5}}(\omega^n+(-\omega)^{-n})$也可以用模幂运算
只是这时需要符号运算比较好。在模不是5 ...
用你的办法！
pari gp的代码
m=1234567891011 /*模*/
r1=Mod(Mod(1,m)*x,x^2-x-1)^(10^14) /*10^14次方*/
print(lift(r1))
输出结果
(12:08) gp > m=1234567891011
%1 = 1234567891011
(12:09) gp > r1=Mod(Mod(1,m)*x,x^2-x-1)^(10^14)
%2 = Mod(Mod(921144120792, 1234567891011)*x + Mod(512884989226, 1234567891011), x^2 - x - 1)
(12:09) gp > print(lift(r1))
Mod(921144120792, 1234567891011)*x + Mod(512884989226, 1234567891011)
根据这个结果，得到
于是x^(10^14)=921144120792*x+512884989226 (mod x^2-x-1) 多项式系数已经对m=1234567891011取模
再把x分别用(1 +Sqrt)/2、(1 - Sqrt)/2替换一下，得到
Clear["Global`*"];(*Clear all variables*)
f=921144120792*x+512884989226
aa=(f/.x->((1+Sqrt)/2))-(f/.x->((1-Sqrt)/2))//Simplify
得到
\
再除以根号5，
得到921144120792（严格地说，这个数需要模一下才是最终结果！）