I’m now using the polynomial $f = X^3 – 3X -1$. Given a number to test $p$, things are good if $f$ is irreducible over $F_{p}$. I’m not too sure why this is so, but it’s also the case for the standard quadratic Frobenius test, that one works in proper extensions of $F_{p}$. Empirically, the primes which do not pass the cubic Frobenius test are those in OEIS sequence A129805, and this makes up 1/3 of the primes (based on the remainder modulo 9): https://oeis.org/A129805 . I suspect that for those primes in sequence A129805, the polynomial $f$ isn’t irreducible over $F_{p}$. The basic Galois-theoretic argument is set up with $\alpha$, $\beta$ and $\gamma$ being three roots where the Frobenius autmorphism permutes the three roots cyclically. This means the roots are not in $F_p$. So one should not worry that these primes do not pass this test. In any case, some Galois theory gives three congruences, one for each of $V_{p}$, $V_{p+1}$, $V_{p+2}$ modulo $p$ where $V_k$ is the sum of the $kth$ powers of the three roots of $f$. Earlier today, I asked a question at Mathematics Stack Exchange, and user Aphelli appears to have solved my problem, showing that $V_{p+2} == 3 (mod p) $ or $V_{p+2}== -6 (mod p)$ and here is a link: https://math.stackexchange.com/questions/4993773/how-does-one-compute-the-value-of-alphap2-betap2-gammap2-in/4993807#4993807.
The cubic Frobenius tests run about 6 times slower than PARI/gp’s ispseudoprime function, but more than 300 times faster than the deterministic, foolproof, test from the isprime function. (on numbers of a few hundred digits). Below,I post the latest code for the cubic Frobenius primality test.
It would be nice to have a variant that works (passes) primes that the current test doesn’t pass. It would also be nice to have proofs relating irreducibility of $f$ over $F_{p}$ to the residue class of $p$ modulo 9. Finally, it would be nice to generalize the cubic Frobenius test to encompass polynomials $f$ of the more general form $f = X^3 – rX -s$. [Note: Computations support the proposition that $f$ is irreducible over $F_{p}$ iff $p=3$ or $p==1 (mod 9)$ or $p==8 (mod 9)$, for an odd prime $p$.]
? isfrob3spec
%761 = (p)->init3(p,3,1);Mpp=power(M,p);output=Mpp*v3;test=(output[2,1]==Mod(-3,p));test=test&&(output[3,1]==Mod(0,p));test=test&&((output[1,1]==Mod(3,p))||(output[1,1]==Mod(-6,p)));return(test)
? init3
%762 = (p,r,s)->M[1,1]=Mod(0,p);M[1,2]=Mod(r,p);M[1,3]=Mod(s,p);M[2,1]=Mod(1,p);M[2,2]=M[1,1];M[2,3]=M[1,1];M[3,1]=M[1,1];M[3,2]=Mod(1,p);M[3,3]=M[1,1];v3[1,1]=Mod(2*r,p);v3[2,1]=Mod(0,p);v3[3,1]=Mod(3,p)
? power
%763 = (a,n)->v=digits(n,2);l=length(v);curpow=a;for(Y=2,l,ap2=curpow*curpow;if(v[Y],ap2=a*ap2);curpow=ap2;);return(curpow)
===================================
functions as of Nov 5:
? \u
disc165649 =
(X)->val=!simple37(X);return(val)
disc3721 =
(X)->val=!simple61(X);return(val)
disc4225 =
(X)->val=!simple13(X);return(val)
disc49 =
(X)->val=!simple7(X);return(val)
disc81 =
(X)->val=!simple9(X);return(val)
init3 =
(p,r,s)->M[1,1]=Mod(0,p);M[1,2]=Mod(r,p);M[1,3]=Mod(s,p);M[2,1]=Mod(1,p);M[2,2]=M[1,1];M[2,3]=M[1,1];M[3,1]=M[1,1];M[3,2]=Mod(1,p);M[3,3]=M[1,1];v3[1,1]=Mod(2*r,p);v3[2,1]=Mod(0,p);v3[3,1]=Mod(3,p)
isfrob3gen =
(p,r,s)->init3(p,r,s);Mpp=power(M,p);output=Mpp*v3;test=(output[2,1]==Mod(-r,p));test=test&&(output[3,1]==Mod(0,p));return(test)
isfrob3select =
(p)->found=0;if(disc49(p),r=7;s=7;found=1);if((!found)&&disc81(p),r=3;s=1;found=1);if((!found)&&disc165649(p),r=37;s=37;found=1);if((!found)&&disc4225(p),r=13;s=13;found=1);if(disc3721(p),r=61;s=183;found=1);if(!found,return(-1));if(found,init3(p,r,s);Mpp=power(M,p);output=Mpp*v3;test=(output[2,1]==Mod(-r,p));test=test&&(output[3,1]==Mod(0,p));return(test))
isfrob3spec =
(p)->init3(p,3,1);Mpp=power(M,p);output=Mpp*v3;test=(output[2,1]==Mod(-3,p));test=test&&(output[3,1]==Mod(0,p));test=test&&((output[1,1]==Mod(3,p))||(output[1,1]==Mod(-6,p)));return(test)
isfrob3spec13 =
(p)->init3(p,13,13);Mpp=power(M,p);output=Mpp*v3;flag=simple13(p);if(!flag,test=(output[2,1]==Mod(-13,p)));if(flag,test=(output[2,1]==Mod(26,p)));;test=test&&(output[3,1]==Mod(0,p));return(test)
p37irr =
(N)->polisirreducible(Mod(1,N)*X^3-Mod(37,N)*X-Mod(37,N))
p7irr =
(N)->polisirreducible(Mod(1,N)*X^3-Mod(7,N)*X-Mod(7,N))
p9irr =
(N)->polisirreducible(Mod(1,N)*X^3-Mod(3,N)*X-Mod(1,N))
power =
(a,n)->v=digits(n,2);l=length(v);curpow=a;for(Y=2,l,ap2=curpow*curpow;if(v[Y],ap2=a*ap2);curpow=ap2;);return(curpow)
simple13 =
(X)->rem=X%13;val=(rem==1)||(rem==5)||(rem==8)||(rem==12);return(val)
simple37 =
(X)->rem=X%37;val=(rem==1)||(rem==6)||(rem==8)||(rem==11)||(rem==14)||(rem==23)||(rem==26)||(rem==29)||(rem==31)||(rem==36)||(rem==10)||(rem==27);return(val)
simple61 =
(X)->rem=X%61;val=(rem==1)||(rem==3)||(rem==8)||(rem==9)||(rem==11)||(rem==23)||(rem==27)||(rem==28)||(rem==33)||(rem==34)||(rem==38)||(rem==50)||(rem==52)||(rem==53)||(rem==58)||(rem==60)||(rem==9)||(rem==52)||(rem==24)||(rem==37)||(rem==20)||(rem==41)||(rem==0);return(val)
simple7 =
(X)->rem=X%7;val=((rem==1)||(rem==6));return(val)
simple9 =
(X)->rem=X%9;val=((rem==1)||(rem==8));return(val)
?
======================================================
functions Nov 5 updated 8:30 pm
? \u
disc10220809 =
(X)->val=!simple139(X);return(val)
disc112225 =
(X)->val=!simple67(X);return(val)
disc165649 =
(X)->val=!simple37(X);return(val)
disc16605625 =
(X)->val=!simple163(X);return(val)
disc17689 =
(X)->val=!simple19(X);return(val)
disc1803649 =
(X)->val=!simple79(X);return(val)
disc3396649 =
(X)->val=!simple97(X);return(val)
disc3721 =
(X)->val=!simple61(X);return(val)
disc4225 =
(X)->val=!simple13(X);return(val)
disc49 =
(X)->val=!simple7(X);return(val)
disc81 =
(X)->val=!simple9(X);return(val)
init3 =
(p,r,s)->M[1,1]=Mod(0,p);M[1,2]=Mod(r,p);M[1,3]=Mod(s,p);M[2,1]=Mod(1,p);M[2,2]=M[1,1];M[2,3]=M[1,1];M[3,1]=M[1,1];M[3,2]=Mod(1,p);M[3,3]=M[1,1];v3[1,1]=Mod(2*r,p);v3[2,1]=Mod(0,p);v3[3,1]=Mod(3,p)
isfrob3gen =
(p,r,s)->init3(p,r,s);Mpp=power(M,p);output=Mpp*v3;test=(output[2,1]==Mod(-r,p));test=test&&(output[3,1]==Mod(0,p));return(test)
isfrob3select =
(p)->found=0;if(disc49(p),r=7;s=7;found=1);if((!found)&&disc81(p),r=3;s=1;found=1);if((!found)&&disc165649(p),r=37;s=37;found=1);if((!found)&&disc4225(p),r=13;s=13;found=1);if(disc3721(p),r=61;s=183;found=1);if(disc3396649(p),r=97;s=97;found=1);if(disc1803649(p),r=79;s=79;found=1);if(disc17689(p),r=19;s=19;found=1);if(disc10220809(p),r=139;s=139;found=1);if(disc16605625(p),r=163;s=163;found=1);if(disc112225(p),r=67;s=201;found=1);if(!found,return(-1));if(found,init3(p,r,s);Mpp=power(M,p);output=Mpp*v3;test=(output[2,1]==Mod(-r,p));test=test&&(output[3,1]==Mod(0,p));return(test))
isfrob3spec =
(p)->init3(p,3,1);Mpp=power(M,p);output=Mpp*v3;test=(output[2,1]==Mod(-3,p));test=test&&(output[3,1]==Mod(0,p));test=test&&((output[1,1]==Mod(3,p))||(output[1,1]==Mod(-6,p)));return(test)
isfrob3spec13 =
(p)->init3(p,13,13);Mpp=power(M,p);output=Mpp*v3;flag=simple13(p);if(!flag,test=(output[2,1]==Mod(-13,p)));if(flag,test=(output[2,1]==Mod(26,p)));;test=test&&(output[3,1]==Mod(0,p));return(test)
p37irr =
(N)->polisirreducible(Mod(1,N)*X^3-Mod(37,N)*X-Mod(37,N))
p7irr =
(N)->polisirreducible(Mod(1,N)*X^3-Mod(7,N)*X-Mod(7,N))
p9irr =
(N)->polisirreducible(Mod(1,N)*X^3-Mod(3,N)*X-Mod(1,N))
power =
(a,n)->v=digits(n,2);l=length(v);curpow=a;for(Y=2,l,ap2=curpow*curpow;if(v[Y],ap2=a*ap2);curpow=ap2;);return(curpow)
simple13 =
(X)->rem=X%13;val=(rem==1)||(rem==5)||(rem==8)||(rem==12);return(val)
simple139 =
(X)->rem=X%139;if(rem==0,val=1);if(rem>0,val=tab139[rem]);return(val)
simple163 =
(X)->rem=X%163;if(rem==0,val=0);if(rem>0,val=tab163[rem]);return(val)
simple19 =
(X)->rem=X%19;if(rem==0,val=0);if(rem>0,val=(rem==1)||(rem==7)||(rem==8)||(rem==11)||(rem==12)||(rem==18));return(val)
simple37 =
(X)->rem=X%37;val=(rem==1)||(rem==6)||(rem==8)||(rem==11)||(rem==14)||(rem==23)||(rem==26)||(rem==29)||(rem==31)||(rem==36)||(rem==10)||(rem==27);return(val)
simple61 =
(X)->rem=X%61;val=(rem==1)||(rem==3)||(rem==8)||(rem==9)||(rem==11)||(rem==23)||(rem==27)||(rem==28)||(rem==33)||(rem==34)||(rem==38)||(rem==50)||(rem==52)||(rem==53)||(rem==58)||(rem==60)||(rem==9)||(rem==52)||(rem==24)||(rem==37)||(rem==20)||(rem==41)||(rem==0);return(val)
simple67 =
(X)->rem=X%67;if(rem==0,val=1);if(rem>0,val=tab67[rem]);return(val)
simple7 =
(X)->rem=X%7;val=((rem==1)||(rem==6));return(val)
simple79 =
(X)->rem=X%79;if(rem==0,val=1);if(rem>0,val=tab79[rem]);return(val)
simple9 =
(X)->rem=X%9;val=((rem==1)||(rem==8));return(val)
simple97 =
(X)->rem=X%97;if(rem==0,val=0);if(rem>0,val=tab97[rem]);return(val)
?
