\\ 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<numprimes, 
		j=j+1; 	
		lastp=(lastp+2*halfprimediff[j]) 
	);
	lastp;
}

\\ Create the extra modulo 30 successive differences starting at the last prime
\\ extra sieve is modulo 30=2*3*5

numextras=(2-1)*(3-1)*(5-1); \\ Units mod 30
extradiffs=vector(numextras);

{
	current=lastp%30;j=current+1;

	i=0;
	while(i<numextras, 
		i=i+1; 
		while( (j%2)*(j%3)*(j%5) == 0,
			j=j+1
		); \\ The mod 30 primality!
		extradiffs[i]=(j-current)%30;
		current=j;
		j=j+1 
	);
}

\\ Trialdivision assuming that n<=(largest prime we have)^2
\\ Only upto a bound B=(1.1)*(last prime)
trialbound=lastp+lastp\10;
trialtop=trialbound^2;

trialdiv(n)=
{
	local(i,divisor,result);

	i=0;divisor=0;result=[-1,0]~;
	until( (result[1]<divisor) || (divisor > 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<k,
				if( p+1 == n, break); \\ "No info"
				p=(p*p)%n;
				e=e+1;
				if( p == 1, return(a));
			);
			if( p+1 != n, return(a)); \\ a^(n-1)/2 != -1 means composite
		);
		c=c-1;
	);
	return(0);
}

\\ Next prime with probability > 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];
	);
}