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

sppq(n,a)=
{ /* Return whether n is a strong pseudoprime to base a (Rabin Miller) */
    local(q, t, b, e);
    q = n-1;  t = 0;  while ( 0==bitand(q,1), q\=2; t+=1 );
    /* here  n==2^t*q+1 */
    b = Mod(a, n)^q;
\\    print1("  b=", lift(b) );
    if ( 1==b, return(1) );
    e = 1;
    while ( e<t,
\\        if ( e>1, print1("  b^",2^(e-1),"=", lift(b)); );
        if( (b==1) || (b==n-1), break(); );
        b *= b;
        e++;
    );
    return( if ( b!=(n-1), 0, 1 ) );
} /* -------- */

rabin_miller(n, na=20)=
{ /* Rabin Miller test */
    local(a);
    for (u=1, na,
        a = prime(u);
        if ( a>n-2,  break() );
        if ( 0==sppq(n, a), return(0) );  /* proven composite */
    );
    return(1); /* composite with probability less than 0.25^na */
} /* -------- */

rabin_miller_v(n, va)=
{
    local(a);
    for (u=1, #va,
        a = va[u];
        if ( a>n-2,  next() );
        if ( 0==sppq(n, a), return(0) );  /* proven composite */
    );
    return(1);
} /* -------- */


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