用mathematica写的Lucas伪素数判定源代码

mathematica

1. Clear["Global`*"];(*Clear all variables*)
   
2. Lucas[n0_]:=Module[{n=n0,m,s},
   
3.     If[Mod[n,2]==0&&n>2,Return[False]];(*排除偶数*)
   
4.     If[IntegerQ[Sqrt[n]],Return[False]];(*如果是完全平方数,返回False*)
   
5.     (*写成n+1=2^s*m的形式*)
   
6.     m=n+1;
   
7.     (*s=0;While[Mod[m,2]==0,m=m/2;s=s+1];*)
   
8.     (*根据P=3 4 5 6 7 Q=1,以及雅克比符号等于-1来找到P,如果d与n有约数,则返回false,且d不应该是n的倍数*)
   
9.     P=3;Q=1;d=P^2-4*Q;JS=JacobiSymbol[d,n];If[JS==0,Return[False]];
   
10.     While[JS==1,P=P+1;d=P^2-4*Q;JS=JacobiSymbol[d,n];If[JS==0,Return[False]]];
    
11.     mi2=IntegerDigits[m,2];(*把m写成二进制的方式*)
    
12.     UU=1;VV=P;(*分别是lucas序列的U(1)与V(1)的值*)
    
13.     Do[
    
14.         Utemp=Mod[UU*VV,n];
    
15.         Vtemp=Mod[VV^2-2,n];(*此处Q=1*)
    
16.         If[mi2[[kk]]==1,
    
17.             uutemp=(P*Utemp+Vtemp);
    
18.             vvtemp=(d*Utemp+P*Vtemp);
    
19.             (*uutemp,vvtemp都可能是奇数,如果是奇数,则加上一个n,这样就是偶数了,
    
20.             下面才能除以2得到整数,至于为什么要这么做,我也不是太清楚为什么,
    
21.             可能是由于n是奇数,模的时候,减去偶数个n奇偶性不变,而减去奇数个n奇偶性变了*)
    
22.             If[OddQ[uutemp],uutemp=uutemp+n];
    
23.             If[OddQ[vvtemp],vvtemp=vvtemp+n];
    
24.             UU=Mod[uutemp/2,n];
    
25.             VV=Mod[vvtemp/2,n],
    
26.             UU=Utemp;
    
27.             VV=Vtemp
    
28.         ],
    
29.         {kk,2,Length@mi2}];(*此处必须从2开始,程序没有错误!*)
    
30.     If[UU==0&&Mod[VV,n]==2,Return[True],Return[False]]
    
31. ]
    
32. Lucas[97]
    
33. Lucas[57]
    
34. Lucas[37]
    
35. (*下面这个强伪素数被判定成合数*)
    
36. Lucas[803837457453639491257079614341942108138837688287558145837488917522\
    
37. 2974273765333652186502336163960045457915042023603208766569966760987284\
    
38. 0439654082329287387918508691668573282677617710293896977394701670823042\
    
39. 8687109997439976544144845341155872450633409279022275296229414984230688\
    
40. 1685404326457534018329786111298960644845216191652872597534901]
    

复制代码