\\pseudoprime base a
{
ps(n,a)=local(out);out=0;
if(Mod(a,n)^(n-1)==Mod(1,n),
  out=1
  );
out
}
\\strong pseudoprime base a
{
sps(n,a)=local(out,r,b,c);out=0;r=valuation(n-1,2);b=(n-1)/2^r;c=Mod(a,n)^b;
if(c==Mod(1,n),
    out=1,
    for(j=0,r-1,
       if(c==Mod(-1,n),
           out=1;break,
	   c=c^2
	   )
       )
  );
out
}
\\Euler pseudo prime base a (Solovay-Strassen test)
{
eps(n,a)=if(Mod(a,n)^((n-1)/2)==Mod(kronecker(a,n),n),1,0)
}
\\smallest witness for n, warning n must be odd and composite
{
witness(n)=local(a);a=1;while(sps(n,a),a++);a
}
{
witness2(n)=local(a);a=2;while(ps(n,a),a++);a
}

\\record=0;for(n=2,10000000,if(n%2==0||isprime(n),,k=witness(n);
\\if(k>record,record=k;print(n," ",record))))

\\The set S of Miller-Rabin bases
{
set_s(n)=local(v);
v=[];for(a=1,n-1,if(gcd(a,n)==1 && sps(n,a)==1,v=concat(v,[a])));v
}
\\The set G of Solovay-Strassen bases
{
set_g(n)=local(v);
v=[];for(a=1,n-1,if(gcd(a,n)==1 && eps(n,a)==1,v=concat(v,[a])));v
}
\\A formula for the size of G
{
size_g(n)=local(m,d,r,del,rmin);m=factor(n);r=valuation(n-1,2);
rmin=valuation(m[1,1]-1,2);
for(j=2,matsize(m)[1],rj=valuation(m[j,1]-1,2);if(rj<rmin,rmin=rj));
del=1;if(rmin==r,del=2,
for(j=1,matsize(m)[1],rj=valuation(m[j,1]-1,2);
if(rj<r && m[j,2]%2,del=1/2;break)));
d=del*prod(j=1,matsize(m)[1],gcd((n-1)/2,m[j,1]-1));
d
}
\\A formula for the size of S.
{
size_s(n)=local(m,d,r,del,rmin,b);m=factor(n);r=valuation(n-1,2);b=(n-1)/2^r;
rmin=valuation(m[1,1]-1,2);
for(j=2,matsize(m)[1],rj=valuation(m[j,1]-1,2);if(rj<rmin,rmin=rj));
d=(1+sum(j=0,rmin-1,2^(matsize(m)[1]*j)))*prod(j=1,matsize(m)[1],gcd(b,m[j,1]-1));
d
}