From: martin@math.rwth-aachen.de (  Martin Schoenert)
Subject: Re: Probabilistic Prime Numbers
Organization: Rechnerbetrieb Informatik  /  RWTH Aachen
Date: 12 Nov 92 08:34:49 GMT

paj@uk.co.gec-mrc (Paul Johnson) writes:

    I am trying to write a small prime number generator using Algorithm P
    in Knuth.

    Note for those without a copy of Knuth: this algorithm uses a randomly
    generated number to perform a primality test.  It can return one of
    two results: "Not prime" or "Probably prime".  If the latter, then the
    chance of error is < 0.25.  Repeated tests with different random
    numbers can reduce this probability very quickly.

    I cannot understand the discussion of this test in Knuth, nor can I
    understand what his variables are.  Is `n' the number under test, or
    is it `q'?  If the former, why the statement "Let n = 1 + (2^k)q" just
    before the statement of Algorithm P.

n is the number to test.  The statement "Let  n = 1  + (2^k)q" means that
you  should  find  an integer k and an odd  integer q such that the above
equation holds.  For example, if n = 273, then k = 4 and q = 17.

He coninues:

    What I would really like is some C or Pascal code which implements
    Algorithm P.  Failing that, a restatement with a more careful
    declaration of variables would be gratefully received.  Sedgwick and
    Numerical Recipies do not deal with primes.  Can anyone help?

I can't  provide  C or Pascal code,  but  it  shouldn't be  difficult  to
translate the  following into C or  Pascal. (Note  that you will  need  a
large integer  arithmetic.  If  you don't  have that you  are limited  to
numbers less  than 2^32 and are better of just trial dividing  by all odd
numbers up to 2^16).

#############################################################################
##
#V  Primes[]  . . . . . . . . . . . . . . . . . . . . . . . .  list of primes
##
##  'Primes' is a set, i.e., sorted list, of the 168 primes less than 1000.
##
##  This is used in 'IsPrimeInt' and 'FactorsInt' to cast  out  small  primes
##  quickly.
##
Primes := [   2,  3,  5,  7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53,
     59, 61, 67, 71, 73, 79, 83, 89, 97,101,103,107,109,113,127,131,137,139,
    149,151,157,163,167,173,179,181,191,193,197,199,211,223,227,229,233,239,
    241,251,257,263,269,271,277,281,283,293,307,311,313,317,331,337,347,349,
    353,359,367,373,379,383,389,397,401,409,419,421,431,433,439,443,449,457,
    461,463,467,479,487,491,499,503,509,521,523,541,547,557,563,569,571,577,
    587,593,599,601,607,613,617,619,631,641,643,647,653,659,661,673,677,683,
    691,701,709,719,727,733,739,743,751,757,761,769,773,787,797,809,811,821,
    823,827,829,839,853,857,859,863,877,881,883,887,907,911,919,929,937,941,
    947,953,967,971,977,983,991,997 ];


#############################################################################
##
#F  PowerModInt(<r>,<e>,<m>)  . . . . . . power of one integer modulo another
##
PowerModInt := function ( r, e, m )
    local   pow, f;

    # reduce r initially
    r := r mod m;

    # handle special case
    if e = 0  then
        return 1;
    fi;

    # now use the repeated squaring method (right-to-left)
    pow := 1;
    f := 2 ^ (LogInt( e, 2 ) + 1);
    while 1 < f  do
        pow := (pow * pow) mod m;
        f := QuoInt( f, 2 );
        if f <= e  then
            pow := (pow * r) mod m;
            e := e - f;
        fi;
    od;

    # return the power
    return pow;
end;


#############################################################################
##
#F  IsPrimeInt( <n> ) . . . . . . . . . . . . . . . . . . .  test for a prime
##
##  'IsPrimeInt' returns 'true' if the  integer <n>  is  a prime  and 'false'
##  otherwise.  By convention 'IsPrimeInt( 0 ) = IsPrimeInt( 1 ) = false' and
##  we define 'IsPrimeInt( -<n> ) = IsPrimeInt( <n> )'.
##
IsPrimeInt := function ( n )
    local  p, e, o, x, i, d;

    # make $n$ positive and handle trivial cases
    if n < 0         then n := -n;       fi;
    if n in Primes   then return true;   fi;
    if n <= 1000     then return false;  fi;

    # do trial divisions by the primes less than 1000
    # faster than anything fancier because $n$ mod <small int> is very fast
    for p  in Primes  do
        if n mod p = 0  then return false;  fi;
        if n < (p+1)^2  then return true;   fi;
    od;

    # find $e$ and $o$ odd such that $n-1 = 2^e * o$
    e := 0;  o := n-1;   while o mod 2 = 0  do e := e+1;  o := o/2;  od;

    # try bases [2,3,5,7]
    for d  in [2,3,5,7]  do

        # look at the seq $2^o, 2^{2 o}, 2^{4 o}, .., 2^{2^e o}=2^{n-1}$
        x := PowerModInt( d, o, n );
        i := 0;
        while i < e  and x <> 1  and x <> n-1  do
            x := x * x mod n;
            i := i + 1;
        od;

        # if it is not of the form $.., -1, 1, 1, ..$ then $n$ is composite
        if not (x = n-1 or (i = 0 and x = 1))  then
            return false;
        fi;

        # there are no strong pseudo-primes to base 2 smaller than 2047
        if d = 2  and n < 2047  then
            return true;
        fi;

    od;

    # after those tests $n$ is probably a prime
    return "n maybe a prime";
end;

Martin.

-- .- .-. - .. -.  .-.. --- ...- . ...  .- -. -. .. -.- .-
Martin Sch"onert,   Martin.Schoenert@Math.RWTH-Aachen.DE,  +49 241 804551
Lehrstuhl D f"ur Mathematik, Templergraben 64, RWTH, D 51 Aachen, Germany