\\% Pocklington-Lehmer test.
\\ Author: Joerg Arndt
\\ License: GPL version 3 or later
\\ online at http://www.jjj.de/pari/
\\ version: 2011-January-19 (12:49)

pocklington_lehmer(F, u, c=10000)=
{ /* Pocklington-Lehmer test for the primality of p=f*u+1.
   * Return last successful base, else zero.
   * F must be the factorization of f.
   * Test bases a=2...c
   * Must have u<f.
   */
    local(n, f, C, p, t, ct);
    n = matsize(F)[1];
    f = prod(j=1, n, F[j,1]^F[j,2]);
    if ( f<=u, return(0) );
    p = f*u + 1;
    C = vector(n, j, (p-1)/F[j,1]);
    ct = n;  \\ number remaining prime divisors of f
    for (a=2, c,
\\    forprime (a=2, c,  \\ seems to be faster
        if ( 1==Mod(a,p)^(p-1),
            for (j=1, n,
                if ( C[j]!=0,  \\ skip entries already removed
\\                    t = component( Mod(a,p)^C[j], 2); 
                    t = lift( Mod(a,p)^C[j] ); 
                    if( 1==gcd(t-1, p),
                        C[j] = 0;  \\ remove entry
                        ct -= 1;   \\ number of remaining entries
                    );
                );
            );
            if ( ct==0, return(a) );
        );
    );
    return( 0 );
} /* ----- */

\\ ==== end of file ====