\section{Probabalistic Primality Testing}
%$Log: Prime.lhs,v $
%Revision 1.1  2004/08/05 11:13:39  malcolm
%Add a regression testsuite for the nhc98 compiler.  It isn't very good,
%but it is better than nothing.  I've been using it for about four years
%on nightly builds, so it's about time it entered the repository!  It
%includes a slightly altered version of the nofib suite.
%Instructions are in the README.
%
%Revision 1.2  1996/07/25 21:32:57  partain
%Bulk of final changes for 2.01
%
%Revision 1.1  1996/01/08 20:04:21  partain
%Initial revision
%
%Revision 1.2  92/06/30  18:34:22  dlester
%Bug fix; now works for 1.
%
%Revision 1.1  92/06/30  15:55:20  dlester
%Initial revision
%

> module Prime (multiTest) where
> import IntLib
> import Random

The function @multiTest@ tests the integer @n@ for primality @k@ times
using the random numbers supplied in @rs@; it returns a @Bool@ and a
depleted sequence of random numbers. When the @Bool@ is @True@ the
number is prime with probability $1 - 4^{-@k@}$; the number is
definitely composite when the @Bool@ is @False@.

> multiTest :: Int -> [Int] -> Integer -> (Bool, [Int])
> multiTest k rs n
>  = if n <= 1 || even n then (n==2, rs) else mTest k rs
>    where mTest 0 rs = (True, rs)
>          mTest k rs = if t then mTest (k-1) rs' else (False, rs')
>                           where (t, rs') = singleTest n (findKQ n) rs

The function @findKQ@ takes an odd integer $n$ and returns the tuple
$(k,q)$ such that $n = q2^k+1$.

> findKQ :: Integer -> (Integer, Integer)
> findKQ n = f (0, (n-1))
>            where f (k,q) = if r == 0 then f (k+1, d) else (k, q)
>                            where (d, r) = q `divMod` 2

The function @singleTest@ takes an odd, positive, @Integer@ @n@ and a
pair of @Integer@'s derived from @n@ by the @findKQ@ function

> singleTest :: Integer -> (Integer, Integer) -> [Int] -> (Bool, [Int])
> singleTest n kq rs
>  = (singleTestX n kq (2+x), rs')
>    where (x, rs')       = random (n-2) rs

Given a value @x@ we can test whether we have a witness to @n@'s
compositeness using @singleTestX@.

> singleTestX n (k, q) x
>  = t == 1 || t == n-1 || witness ts
>    where (t:ts)         = take (fromInteger k) (iterate square (powerMod x q n))
>          witness []     = False
>          witness (t:ts) = if t == n-1 then True       else
>                           if t == 1   then False      else
>                                            witness ts
>          square x       = (x*x) `mod` n

The @random@ function takes a number @n@ and a list of pseudo-random
numbers @rs@ and returns a tuple consisting of an @Integer@ $x$ in the
range $0 \leq x < @n@$, and the remainder of the pseudo-random
numbers.

> random :: Integer -> [Int] -> (Integer, [Int])
> random n rs = (makeNumber 65536 (uniform ns rs1), rs2)
>               where ns        = chop 65536 n
>                     (rs1,rs2) = splitAt (length ns) rs

The @uniform@ function generates a sequence of @Integer@'s such that,
when considered as a sequence of digits, we generate a number uniform
in the range @0..ns@ from the random numbers @rs@.

> uniform :: [Integer] -> [Int] -> [Integer]
> uniform [n]    [r]    = [toInteger r `mod` n]
> uniform (n:ns) (r:rs) = if t == n then t: uniform ns rs
>                                   else t: map ((`mod` 65536). toInteger) rs
>                         where t  = toInteger r `mod` (n+1)

This concludes the Prime testing algorithm.