\\lucas2.gp - Fast computation of Lucas sequences.
\\You can get a copy of this file  via anonymous ftp at 
\\math.math.ucl.ac.be in the directory /pub/joye/pari.

\\lucas2(k,P,Q) returns the k-th terms of the Lucas sequences
\\{U_k(P,Q)} and {V_k(P,Q)} with parameters P and Q. 

\\Usage:
\\   lucas2(k,P,Q)=k-th term of U_k(P,Q) and V_k(P,Q)

{
        lucas2(k,P,Q,     n,s,j,Ql,Qh,Uh,Vl,Vh) =
        n = floor(log(k)/log(2)) + 1;
        s = 0;
        while( bittest(k,s) == 0, s = s + 1 );
        Uh = 1; Vl = 2; Vh = P;
        Ql = 1; Qh = 1;
        forstep( j = n-1, s+1, -1,
                Ql = Ql * Qh;
                if( bittest(k,j),
                        Qh = Ql * Q;
                        Uh = Uh * Vh;
                        Vl = Vh * Vl - P * Ql;
                        Vh = Vh * Vh - 2 * Qh,

                        Uh = Uh * Vl - Ql;
                        Vh = Vh * Vl - P * Ql;
                        Vl = Vl * Vl - 2 * Ql;
                        Qh = Ql;
                );
        );
        Ql = Ql * Qh; Qh = Ql * Q;
        Uh = Uh * Vl - Ql;
        Vl = Vh * Vl - P * Ql;
        Ql = Ql * Qh;
        for( j = 1, s,
                Uh = Uh * Vl;
                Vl = Vl * Vl - 2 * Ql;
                Ql = Ql * Ql;
        );
        [Uh,Vl]
}


\\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
}
\\Baillie-Pomerance-Selfridge-Wagstaff test

{
psw(n)=local(out,d);out=1;if(n%2==0 && n != 2, out =0);
if(n%2, 
if(issquare(n),out=0,out=sps(n,2));
if(out==1,
  forstep(j=5,n+2,2,d=(-1)^((j-1)/2)*j;if(kronecker(d,n)==-1,break));
  if(lucas2(n+1,Mod(1,n),Mod((1-d)/4,n))[1]!=0,out=0)
  ));
out
}