PrimeQ神奇,是否有重大bug?

无心人

388342位十进制的整数，只能用概率算法做的，一切确定性方法，都会完蛋~

dlpg070

在寻找新的最大素数过程中,
对于 N = k * 2^n +/- 1 (对于奇数 k 和 n > 1)
似乎多数用确定型素性测试
例如:
https://primes.utm.edu/primes/page.php?id=121330


N= 2618163402417*2^1290000 - 1(388342位十进制)
Proof-code(s): (*):        L927 : Brown1, TwinGen, PrimeGrid, LLR
(A proof code list the persons, programs and projects involved in a proof. )
验证用: Proth prime test & Lucas Lehmer Riesel prime test
确为素数!!!
mathematica 确实无能为力

mathematica

dlpg070 发表于 2019-9-18 14:34
在寻找新的最大素数过程中,
对于 N = k * 2^n +/- 1 (对于奇数 k 和 n > 1)
似乎多数用确定型素性测 ...

baillie psw，这个算法最好了，再加上几十次miller rabin，就可以了

mathematica

dlpg070 发表于 2019-9-18 14:34
在寻找新的最大素数过程中,
对于 N = k * 2^n +/- 1 (对于奇数 k 和 n > 1)
似乎多数用确定型素性测 ...

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.             If[OddQ[UUtemp],UUtemp=UUtemp+n];
    
26.             If[OddQ[VVtemp],VVtemp=VVtemp+n];
    
27.             UU=Mod[UUtemp/2,n];
    
28.             VV=Mod[VVtemp/2,n]
    
29.         ],
    
30.         {kk,2,Length@mi2}];(*此处必须从2开始,程序没有错误!*)
    
31.     If[UU==0&&Mod[VV,n]==2,Return[True],Return[False]]
    
32. ]
    
33. 
    
34. (*miller rabin测试,n0是被测试的整数,a0是选择的基*)
    
35. MillerRabin[n0_,a0_]:=Module[{n=n0,a=a0,s,m,t1,k},
    
36.     s=0;m=n-1;While[Mod[m,2]==0,m=m/2;s=s+1];
    
37.     t1=PowerMod[a,m,n];
    
38.     If[t1==1,Return[True]];
    
39.     k=0;While[k<s-1&&t1!=n-1,k=k+1;t1=Mod[t1^2,n]];
    
40.     If[t1==n-1,Return[True],Return[False]]
    
41. ]
    
42. 
    

复制代码

我的baillie psw的mathematica版本算法，只要你实现了大整数的加减乘除，这个就比较简单了！