Latest PARI/gp code for the cubic Frobenius primality test

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)

disc16785409 =
  (X)->val=!simple241(X);return(val)

disc17689 =
  (X)->val=!simple19(X);return(val)

disc1792921 =
  (X)->val=!simple103(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)

isfrob3seldisc =
  (p)->found=0;if((p==2)||(p==3),return(1));if(((p%2)==1)&&((p%3)>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(disc1792921(p),r=103;s=309;found=1);if(disc16785409(p),r=241;s=1205;found=1);if(!found,return(-1));if(found,return(1));return(0));return(0)

isfrob3select =
  (p)->found=0;if((p==2)||(p==3),return(1));if(((p%2)==1)&&((p%3)>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(disc1792921(p),r=103;s=309;found=1);if(disc16785409(p),r=241;s=1205;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));return(0));return(0)

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)

simple103 =
  (X)->rem=X%103;if(rem==0,val=1);if(rem>0,val=tab103[rem]);return(val)

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)

simple241 =
  (X)->rem=X%241;if(rem==0,val=1);if(rem>0,val=tab241[rem]);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)||(rem==0));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=0;if(rem==0,val=1);if(rem>0,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)