Baillie_PSW素性判定的mathematica代码

mathematica · 发表于 2019-2-27 14:37:40

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

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

×

mathematica本身有了大整数加减乘除取模，所以在mathematica之上要简单很多，很久以前就有实现这个算法的心愿，今年终于完成了

Clear["Global`*"];(*Clear all variables*)
Lucas[n0_]:=Module[
(*指定局部变量*)
{n=n0,m,s,P,Q,d,JS,mi2,UU,VV,uutemp,vvtemp,Utemp,Vtemp,kk},
If[Mod[n,2]==0&&n>2,Return[False]];(*排除偶数*)
If[IntegerQ[Sqrt[n]],Return[False]];(*如果是完全平方数,返回False*)
(*写成n+1=2^s*m的形式*)
m=n+1;
(*s=0;While[Mod[m,2]==0,m=m/2;s=s+1];*)
(*根据P=3 4 5 6 7 Q=1,以及雅克比符号等于-1来找到P,如果d与n有约数,则返回false,且d不应该是n的倍数*)
(*此处如果n很小,而d较大且是n的倍数,这时可能存在n是质数,而雅克比符号等于零的情况,这种情况需要特殊处理一下*)
P=3;Q=1;d=P^2-4*Q;JS=JacobiSymbol[d,n];If[JS==0,Return[False]];
While[JS==1,P=P+1;d=P^2-4*Q;JS=JacobiSymbol[d,n];If[JS==0,Return[False]]];
mi2=IntegerDigits[m,2];(*把m写成二进制的方式*)
UU=1;VV=P;(*分别是lucas序列的U(1)与V(1)的值*)
Do[
Utemp=Mod[UU*VV,n];
Vtemp=Mod[VV^2-2,n];(*此处Q=1*)
If[mi2[[kk]]==1,
uutemp=(P*Utemp+Vtemp);
vvtemp=(d*Utemp+P*Vtemp);
(*uutemp,vvtemp都可能是奇数,如果是奇数,则加上一个n,这样就是偶数了,
下面才能除以2得到整数,至于为什么要这么做,我也不是太清楚为什么,
可能是由于n是奇数,模的时候,减去偶数个n奇偶性不变,而减去奇数个n奇偶性变了*)
If[OddQ[uutemp],uutemp=uutemp+n];
If[OddQ[vvtemp],vvtemp=vvtemp+n];
UU=Mod[uutemp/2,n];
VV=Mod[vvtemp/2,n],
UU=Utemp;
VV=Vtemp
],
{kk,2,Length@mi2}];(*此处必须从2开始,程序没有错误!*)
If[UU==0&&Mod[VV,n]==2,Return[True],Return[False]]
]
(*miller rabin测试,n0是被测试的整数,a0是选择的基*)
MillerRabin[n0_,a0_]:=Module[{n=n0,a=a0,s,m,t1,k},
s=0;m=n-1;While[Mod[m,2]==0,m=m/2;s=s+1];
t1=PowerMod[a,m,n];
If[t1==1,Return[True]];
k=0;While[k<s-1&&t1!=n-1,k=k+1;t1=Mod[t1^2,n]];
If[t1==n-1,Return[True],Return[False]]
]

复制代码

一些细节问题都写在代码注释里了

https://bbs.emath.ac.cn/forum.ph ... 757&fromuid=865
https://bbs.emath.ac.cn/forum.ph ... 734&fromuid=865

共享给感兴趣的人，也当做玩mathematica软件的一个练习吧！

mathematica · 发表于 2020-3-2 10:47:12

Clear["Global`*"];(*Clear all variables*)
Lucas[n0_]:=Module[
(*指定局部变量*)
{n=n0,m,s,P,Q,d,JS,mi2,UU,VV,UUtemp,VVtemp,kk},
If[Mod[n,2]==0&&n>2,Return[False]];(*排除偶数*)
If[IntegerQ[Sqrt[n]],Return[False]];(*如果是完全平方数,返回False*)
(*写成n+1=2^s*m的形式*)
m=n+1;
(*s=0;While[Mod[m,2]==0,m=m/2;s=s+1];*)
(*根据P=3 4 5 6 7 Q=1,以及雅克比符号等于-1来找到P,如果d与n有约数,则返回false,且d不应该是n的倍数*)
(*此处如果n很小,而d较大且是n的倍数,这时可能存在n是质数,而雅克比符号等于零的情况,这种情况需要特殊处理一下*)
P=3;Q=1;d=P^2-4*Q;JS=JacobiSymbol[d,n];If[JS==0,Return[False]];
While[JS==1,P=P+1;d=P^2-4*Q;JS=JacobiSymbol[d,n];If[JS==0,Return[False]]];
mi2=IntegerDigits[m,2];(*把m写成二进制的方式*)
UU=1;VV=P;(*分别是lucas序列的U(1)与V(1)的值*)
Do[
UU=Mod[UU*VV,n];
VV=Mod[VV^2-2,n];(*此处Q=1*)
If[mi2[[kk]]==1,
UUtemp=(P*UU+VV);
VVtemp=(d*UU+P*VV);
(*UUtemp,VVtemp都可能是奇数,如果是奇数,则加上一个n,这样就是偶数了,
下面才能除以2得到整数,至于为什么要这么做,我也不是太清楚为什么,
可能是由于n是奇数,模的时候,减去偶数个n奇偶性不变,而减去奇数个n奇偶性变了*)
If[OddQ[UUtemp],UUtemp=UUtemp+n];
If[OddQ[VVtemp],VVtemp=VVtemp+n];
UU=Mod[UUtemp/2,n];
VV=Mod[VVtemp/2,n]
],
{kk,2,Length@mi2}];(*此处必须从2开始,程序没有错误!*)
If[UU==0&&Mod[VV,n]==2,Return[True],Return[False]]
]
(*miller rabin测试,n0是被测试的整数,a0是选择的基*)
MillerRabin[n0_,a0_]:=Module[{n=n0,a=a0,s,m,t1,k},
s=0;m=n-1;While[Mod[m,2]==0,m=m/2;s=s+1];
t1=PowerMod[a,m,n];
If[t1==1,Return[True]];
k=0;While[k<s-1&&t1!=n-1,k=k+1;t1=Mod[t1^2,n]];
If[t1==n-1,Return[True],Return[False]]
]

复制代码

进一步把代码简化了

mathematica · 发表于 2020-3-2 12:49:32

本帖最后由 mathematica 于 2020-3-2 12:57 编辑

这个文章的说明了为什么要把lucas伪素数与强伪素数相结合！

账号		自动登录	找回密码
密码			欢迎注册

[原创] Baillie_PSW素性判定的mathematica代码

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

点评