\\ Testing of primes in pari \\ Read in the list of primes stored \\ as half the difference with the previous prime. \r savedprimes.gp; \\ Keep track of the last prime we have lastp= { lastp=0;j=0; while(j trialbound) , if( i < numprimes, \\ Have we exhausted all our primes? i=i+1; divisor=divisor+2*halfprimediff[i], \\ cycle over extradiffs i=(i%numextras)+1; divisor=divisor+extradiffs[i] ); result=divrem(n,divisor); if( result[2] == 0, return(divisor); ); ); return(0); } \\ Left-right exponentiation lrexp(base,expon,modulus)= { local(e,u); e=length(expon)*32; \\ This is the maximum number of bits in expon u=1; while(e>=0, u=(u*u)%modulus; if(bittest(expon,e) == 1, u=(u*base)%modulus ); e=e-1 ); return(u); } \\ Miller-Rabin test run on n c times millrabtest(n,c)= { local(q,k,a,e,p); k=0;q=(n-1)\2; while((q%2) == 0, k=k+1; q=q\2); \\ Now n=1+q2^(k+1) \\ a=0; until(c <= 0, \\ a=a+2*halfprimediff[i]; Why not take prime bases? Don't know! a=random(); p=lrexp(a,q,n); if( p != 1, e=0;while( e 1 - (1/4)^c findnext(n,c)= { \\ local(pos,resid); if( (n%2) == 0, n=n+1); \\ Make it odd while( (n%3)*(n%5) == 0, n=n+2); \\ make it prime to 30 resid=7;pos=0; while( 1, if( (n%30)==resid, break, pos=pos+1;resid=(resid+extradiffs[pos])%30 ) ); \\ Found the position in extradiffs while( 1, if( trialdiv(n) == 0, if ( n < trialtop, return(n) ); if ( millrabtest(n,c) == 0, return(n) ); ); pos=(pos%numextras)+1; n=n+extradiffs[pos]; ); }