用mathematica写的Lucas伪素数判定源代码
Clear["Global`*"];(*Clear all variables*)Lucas:=Module[{n=n0,m,s},
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的倍数*)
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[
Utemp=Mod;
Vtemp=Mod;(*此处Q=1*)
If]==1,
uutemp=(P*Utemp+Vtemp);
vvtemp=(d*Utemp+P*Vtemp);
(*uutemp,vvtemp都可能是奇数,如果是奇数,则加上一个n,这样就是偶数了,
下面才能除以2得到整数,至于为什么要这么做,我也不是太清楚为什么,
可能是由于n是奇数,模的时候,减去偶数个n奇偶性不变,而减去奇数个n奇偶性变了*)
If,uutemp=uutemp+n];
If,vvtemp=vvtemp+n];
UU=Mod;
VV=Mod,
UU=Utemp;
VV=Vtemp
],
{kk,2,Length@mi2}];(*此处必须从2开始,程序没有错误!*)
If==2,Return,Return]
]
Lucas
Lucas
Lucas
(*下面这个强伪素数被判定成合数*)
Lucas[803837457453639491257079614341942108138837688287558145837488917522\
2974273765333652186502336163960045457915042023603208766569966760987284\
0439654082329287387918508691668573282677617710293896977394701670823042\
8687109997439976544144845341155872450633409279022275296229414984230688\
1685404326457534018329786111298960644845216191652872597534901]
页:
[1]