Baillie_PSW素性判定的mathematica代码（有详细解释）

nyy 发表于 2022-8-24 15:32:32

Clear["Global`*"];(*Clear all variables*)
Lucas:=Module[
(*指定局部变量*)
{n=n0,m,s,P,Q,d,JS,mi2,UU,VV,UUtemp,VVtemp,kk},
If==0&&n>2,Return];(*排除偶数*)
If],Return];(*如果是完全平方数,返回False*)
(*写成n+1=2^s*m的形式*)
m=n+1;
(*s=0;While==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;If];
While;If]];
mi2=IntegerDigits;(*把m写成二进制的方式*)
UU=1;VV=P;(*分别是lucas序列的U(1)与V(1)的值*)
Do[
   UU=Mod;
   VV=Mod;(*此处Q=1*)
   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
   ],
   {kk,2,Length@mi2}];(*此处必须从2开始,程序没有错误!*)
If==2,Return,Return]
]

(*miller rabin测试,n0是被测试的整数,a0是选择的基*)
MillerRabin:=Module[{n=n0,a=a0,s,m,t1,k},
s=0;m=n-1;While==0,m=m/2;s=s+1];
t1=PowerMod;
If];
k=0;While];
If,Return]
]

nyy 发表于 2022-8-24 15:33:23

         (*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同理
         *)

感谢@无心人让我彻底搞明白了我曾经搞不明白的！

nyy 发表于 2022-8-26 15:41:03

今年彻底整明白lucas算法，很开心!

nyy 发表于 2022-8-26 16:28:49

ExtraLucas:=Module[
(*指定局部变量*)
{n=n0,m,s,P,Q,d,JS,mi2,UU,VV,UUtemp,VVtemp,kk,k},
If==0&&n>2,Return];(*排除偶数*)
If],Return];(*如果是完全平方数,返回False*)
(*写成n+1=2^s*m的形式,其中m是奇数*)
m=n+1;
s=0;While==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;If];
While;If]];
mi2=IntegerDigits;(*把m写成二进制的方式*)
UU=1;VV=P;(*分别是lucas序列的U(1)与V(1)的值*)
Do[
   UU=Mod;
   VV=Mod;(*此处Q=1*)
   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
   ],
   {kk,2,Length@mi2}];(*此处必须从2开始,程序没有错误!*)
(*UU=0的同时VV必须是正负2,才返回True*)
If==0,Mod==0]],Return];
k=0;While];(*计算VV的值*)
If,Return](*判定VV的值是否等于零*)
]

nyy 发表于 2022-8-26 16:40:11

本帖最后由 nyy 于 2022-8-26 16:41 编辑

论文链接
https://www.ams.org/journals/mcom/2001-70-234/S0025-5718-00-01197-2/S0025-5718-00-01197-2.pdf

郭先强，你用的是
k=0;While];(*计算VV的值*)
（此处运行后k有可能达到s-1）
还是
k=0;While];(*计算VV的值*)

你用的是哪个？
@gxqcn

nyy 发表于 2022-8-26 16:48:42

(*全部是Strong Lucas pseudoprimes*)
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};
bbb=Lucas[#]&/@aaa
ccc=ExtraLucas[#]&/@aaa

(*全部是Lucas pseudoprimes*)
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}
ccc=ExtraLucas[#]&/@aaa
Tally

(*全部是Extra strong Lucas pseudoprimes*)
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}
ccc=ExtraLucas[#]&/@aaa
Tally

页: [1]

数学研发论坛's Archiver

Baillie_PSW素性判定的mathematica代码（有详细解释）