谁能找到通过下面lucas素性测试的Lucas伪素数？

nyy

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

复制代码

nyy

1. (*全部是Lucas pseudoprimes*)
   
2. 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}
   
3. bbb=StrongLucas[#]&/@aaa
   
4. Tally[bbb](*统计判定结果*)
   

复制代码

用lucas伪素数测试，全部给出false

nyy

1. (*生成第2000个到第2200个素数*)
   
2. ccc=Prime@Range[2000,2200];
   
3. ddd=StrongLucas[#]&/@ccc(*对素数进行判定*)
   
4. Tally[ddd](*统计判定结果*)
   

复制代码

用201个素数给出测试，结果证明很可能算法是正确的！
全部给出true

nyy

1. (*全部是以2位底的强伪素数*)
   
2. aaa={2047,3277,4033,4681,8321,15841,29341,42799,49141,52633,65281,74665,80581,85489,88357,90751,104653,130561,196093,220729,233017,252601,253241,256999,271951,280601,314821,357761,390937,458989,476971,486737}
   
3. bbb=StrongLucas[#]&/@aaa(*对素数进行判定*)
   
4. Tally[bbb](*统计判定结果*)
   

复制代码

用以2位底的强伪素数
结果全部false

nyy

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=StrongLucas[#]&/@aaa
   
4. Tally[bbb]
   

复制代码

用Strong Lucas pseudoprimes结果全部给出false

nyy

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

复制代码

用Extra strong Lucas pseudoprimes
结果全部给出false

nyy

@无心人 来找个试试看！

nyy

这些Lucas强伪素数有什么规律?
https://bbs.emath.ac.cn/forum.ph ... 8&fromuid=14149
(出处: 数学研发论坛)
http://ntheory.org/data/slpsps.txt
这儿有178118个lucas强伪素数，居然没有一个通过测试！

代码如下

1. (*读取强伪素数,一行一个,所以需要Flatten一下,第一行数据仅仅只是为了测试import函数*)
   
2. data=Flatten@Import["C:/Users/Administrator/Desktop/0008/123.txt","Data"]
   
3. data=Flatten@Import["C:/Users/Administrator/Desktop/0008/slpsps.txt","Data"];
   
4. aaa=Select[data,StrongLucas](*选择通过测试的强伪素数*)
   

复制代码

nyy

http://ntheory.org/data/slpsps-baillie.txt

这儿的474971个lucas强伪素数，同样全部被干掉了！

1. (*读取强伪素数,一行一个,所以需要Flatten一下*)
   
2. data=Flatten@Import["C:/Users/Administrator/Desktop/0008/123.txt","Data"]
   
3. data=Flatten@Import["C:/Users/Administrator/Desktop/0008/slpsps.txt","Data"];
   
4. data=Flatten@Import["C:/Users/Administrator/Desktop/0008/slpsps-baillie.txt","Data"];
   
5. aaa=Select[data,StrongLucas[#]&](*选择通过测试的强伪素数*)
   
6. aaa=Select[data[[1;;20]],StrongLucas[#]&](*选择通过测试的强伪素数*)
   

复制代码

nyy

http://ntheory.org/data/spsps.txt
这儿是以2为底的强伪素数。
总共有419489个，其中我还添加了两个素数在里面，就为了测试我写的算法能否把素数捡出来，
测试结果告诉我：所有的强伪素数都被干掉了，唯独留下了两个素数！

1. data=Flatten@Import["C:/Users/Administrator/Desktop/0008/spsps.txt","Data"];
   
2. aaa=Select[data,StrongLucas[#]&](*选择通过测试的强伪素数*)
   

复制代码