找回密码
 欢迎注册
查看: 1634|回复: 5

[原创] Baillie_PSW素性判定的mathematica代码(有详细解释)

[复制链接]
发表于 2022-8-24 15:32:32 | 显示全部楼层 |阅读模式

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

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

×
  1. Clear["Global`*"];(*Clear all variables*)
  2. Lucas[n0_]:=Module[
  3.     (*指定局部变量*)
  4.     {n=n0,m,s,P,Q,d,JS,mi2,UU,VV,UUtemp,VVtemp,kk},
  5.     If[Mod[n,2]==0&&n>2,Return[False]];(*排除偶数*)
  6.     If[IntegerQ[Sqrt[n]],Return[False]];(*如果是完全平方数,返回False*)
  7.     (*写成n+1=2^s*m的形式*)
  8.     m=n+1;
  9.     (*s=0;While[Mod[m,2]==0,m=m/2;s=s+1];*)
  10.     (*根据P=3 4 5 6 7 Q=1,以及雅克比符号等于-1来找到P,如果d与n有约数,则返回false,且d不应该是n的倍数*)
  11.     (*此处如果n很小,而d较大且是n的倍数,这时可能存在n是质数,而雅克比符号等于零的情况,这种情况需要特殊处理一下*)
  12.     P=3;Q=1;d=P^2-4*Q;JS=JacobiSymbol[d,n];If[JS==0,Return[False]];
  13.     While[JS==1,P=P+1;d=P^2-4*Q;JS=JacobiSymbol[d,n];If[JS==0,Return[False]]];
  14.     mi2=IntegerDigits[m,2];(*把m写成二进制的方式*)
  15.     UU=1;VV=P;(*分别是lucas序列的U(1)与V(1)的值*)
  16.     Do[
  17.         UU=Mod[UU*VV,n];
  18.         VV=Mod[VV^2-2,n];(*此处Q=1*)
  19.         If[mi2[[kk]]==1,
  20.             UUtemp=(P*UU+VV);
  21.             VVtemp=(d*UU+P*VV);
  22.             (*UUtemp,VVtemp都可能是奇数,如果是奇数,则加上一个n,这样就是偶数了,
  23.             下面才能除以2得到整数,至于为什么要这么做,我也不是太清楚为什么,
  24.             可能是由于n是奇数,模的时候,减去偶数个n奇偶性不变,而减去奇数个n奇偶性变了.
  25.             20220824补充:除以2,其实是乘以2模n的逆,而2模n的逆=(n+1)/2,这里n肯定是奇数
  26.             当UUtemp=2k(偶数的时候),2k*(n+1)/2=k*(n+1)=k(mod n),相当于直接除以2,
  27.             当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,
  28.             VVtemp同理
  29.             *)
  30.             If[OddQ[UUtemp],UUtemp=UUtemp+n];
  31.             If[OddQ[VVtemp],VVtemp=VVtemp+n];
  32.             UU=Mod[UUtemp/2,n];
  33.             VV=Mod[VVtemp/2,n]
  34.         ],
  35.         {kk,2,Length@mi2}];(*此处必须从2开始,程序没有错误!*)
  36.     If[UU==0&&Mod[VV,n]==2,Return[True],Return[False]]
  37. ]

  38. (*miller rabin测试,n0是被测试的整数,a0是选择的基*)
  39. MillerRabin[n0_,a0_]:=Module[{n=n0,a=a0,s,m,t1,k},
  40.     s=0;m=n-1;While[Mod[m,2]==0,m=m/2;s=s+1];
  41.     t1=PowerMod[a,m,n];
  42.     If[t1==1,Return[True]];
  43.     k=0;While[k<s-1&&t1!=n-1,k=k+1;t1=Mod[t1^2,n]];
  44.     If[t1==n-1,Return[True],Return[False]]
  45. ]

复制代码
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2022-8-24 15:33:23 | 显示全部楼层
  1.             (*UUtemp,VVtemp都可能是奇数,如果是奇数,则加上一个n,这样就是偶数了,
  2.             下面才能除以2得到整数,至于为什么要这么做,我也不是太清楚为什么,
  3.             可能是由于n是奇数,模的时候,减去偶数个n奇偶性不变,而减去奇数个n奇偶性变了.
  4.             20220824补充:除以2,其实是乘以2模n的逆,而2模n的逆=(n+1)/2,这里n肯定是奇数
  5.             当UUtemp=2k(偶数的时候),2k*(n+1)/2=k*(n+1)=k(mod n),相当于直接除以2,
  6.             当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,
  7.             VVtemp同理
  8.             *)
复制代码


感谢@无心人 让我彻底搞明白了我曾经搞不明白的!
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2022-8-26 15:41:03 | 显示全部楼层
今年彻底整明白lucas算法,很开心!
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2022-8-26 16:28:49 | 显示全部楼层
  1. ExtraLucas[n0_]:=Module[
  2.     (*指定局部变量*)
  3.     {n=n0,m,s,P,Q,d,JS,mi2,UU,VV,UUtemp,VVtemp,kk,k},
  4.     If[Mod[n,2]==0&&n>2,Return[False]];(*排除偶数*)
  5.     If[IntegerQ[Sqrt[n]],Return[False]];(*如果是完全平方数,返回False*)
  6.     (*写成n+1=2^s*m的形式,其中m是奇数*)
  7.     m=n+1;
  8.     s=0;While[Mod[m,2]==0,m=m/2;s=s+1];
  9.     (*根据P=3 4 5 6 7 Q=1,以及雅克比符号等于-1来找到P,如果d与n有约数,则返回false,且d不应该是n的倍数*)
  10.     (*此处如果n很小,而d较大且是n的倍数,这时可能存在n是质数,而雅克比符号等于零的情况,这种情况需要特殊处理一下*)
  11.     P=3;Q=1;d=P^2-4*Q;JS=JacobiSymbol[d,n];If[JS==0,Return[False]];
  12.     While[JS==1,P=P+1;d=P^2-4*Q;JS=JacobiSymbol[d,n];If[JS==0,Return[False]]];
  13.     mi2=IntegerDigits[m,2];(*把m写成二进制的方式*)
  14.     UU=1;VV=P;(*分别是lucas序列的U(1)与V(1)的值*)
  15.     Do[
  16.         UU=Mod[UU*VV,n];
  17.         VV=Mod[VV^2-2,n];(*此处Q=1*)
  18.         If[mi2[[kk]]==1,
  19.             UUtemp=(P*UU+VV);
  20.             VVtemp=(d*UU+P*VV);
  21.             (*UUtemp,VVtemp都可能是奇数,如果是奇数,则加上一个n,这样就是偶数了,
  22.             下面才能除以2得到整数,至于为什么要这么做,我也不是太清楚为什么,
  23.             可能是由于n是奇数,模的时候,减去偶数个n奇偶性不变,而减去奇数个n奇偶性变了.
  24.             20220824补充:除以2,其实是乘以2模n的逆,而2模n的逆=(n+1)/2,这里n肯定是奇数
  25.             当UUtemp=2k(偶数的时候),2k*(n+1)/2=k*(n+1)=k(mod n),相当于直接除以2,
  26.             当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,
  27.             VVtemp同理
  28.             *)
  29.             If[OddQ[UUtemp],UUtemp=UUtemp+n];
  30.             If[OddQ[VVtemp],VVtemp=VVtemp+n];
  31.             UU=Mod[UUtemp/2,n];
  32.             VV=Mod[VVtemp/2,n]
  33.         ],
  34.         {kk,2,Length@mi2}];(*此处必须从2开始,程序没有错误!*)
  35.     (*UU[m]=0的同时VV[m]必须是正负2,才返回True*)
  36.     If[And[UU==0,Or[Mod[VV-2,n]==0,Mod[VV+2,n]==0]],Return[True]];
  37.     k=0;While[k<s-1&&VV!=0,k=k+1;VV=Mod[VV^2-2,n]];(*计算VV[m*2^k]的值*)
  38.     If[VV==0,Return[True],Return[False]](*判定VV[m*2^k]的值是否等于零*)
  39. ]
复制代码
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2022-8-26 16:40:11 | 显示全部楼层
本帖最后由 nyy 于 2022-8-26 16:41 编辑

论文链接
https://www.ams.org/journals/mco ... 5718-00-01197-2.pdf

郭先强,你用的是
  1.     k=0;While[k<s-1&&VV!=0,k=k+1;VV=Mod[VV^2-2,n]];(*计算VV[m*2^k]的值*)
复制代码
(此处运行后k有可能达到s-1)
还是
  1.     k=0;While[k<s-2&&VV!=0,k=k+1;VV=Mod[VV^2-2,n]];(*计算VV[m*2^k]的值*)
复制代码


你用的是哪个?
@gxqcn
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2022-8-26 16:48:42 | 显示全部楼层
  1. (*全部是Strong Lucas pseudoprimes*)
  2. aaa={5459,5777,10877,16109,18971,22499,24569,25199,40309,58519,75077,97439,100127,113573,115639,130139,155819,158399,161027,162133,176399,176471,189419,192509,197801,224369,230691,231703,243629,253259,268349,288919,313499,324899};
  3. bbb=Lucas[#]&/@aaa
  4. ccc=ExtraLucas[#]&/@aaa

  5. (*全部是Lucas pseudoprimes*)
  6. aaa={323,377,1159,1829,3827,5459,5777,9071,9179,10877,11419,11663,13919,14839,16109,16211,18407,18971,19043,22499,23407,24569,25199,25877,26069,27323,32759,34943,35207,39059,39203,39689,40309,44099,46979,47879}
  7. ccc=ExtraLucas[#]&/@aaa
  8. Tally[ccc]

  9. (*全部是Extra strong Lucas pseudoprimes*)
  10. aaa={989,3239,5777,10877,27971,29681,30739,31631,39059,72389,73919,75077,100127,113573,125249,137549,137801,153931,155819,161027,162133,189419,218321,231703,249331,370229,429479,430127,459191,473891,480689,600059,621781,632249,635627}
  11. ccc=ExtraLucas[#]&/@aaa
  12. Tally[ccc]
复制代码
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

小黑屋|手机版|数学研发网 ( 苏ICP备07505100号 )

GMT+8, 2024-3-29 17:31 , Processed in 0.044577 second(s), 16 queries .

Powered by Discuz! X3.5

© 2001-2024 Discuz! Team.

快速回复 返回顶部 返回列表