[project @ 1996-07-25 21:02:03 by partain]
[nofib.git] / spectral / primetest / Prime.lhs
1 \section{Probabalistic Primality Testing}
2 %$Log: Prime.lhs,v $
3 %Revision 1.2  1996/07/25 21:32:57  partain
4 %Bulk of final changes for 2.01
5 %
6 %Revision 1.1  1996/01/08 20:04:21  partain
7 %Initial revision
8 %
9 %Revision 1.2  92/06/30  18:34:22  dlester
10 %Bug fix; now works for 1.
11 %
12 %Revision 1.1  92/06/30  15:55:20  dlester
13 %Initial revision
14 %
15
16 > module Prime (multiTest) where
17 > import IntLib
18 > import Random
19
20 The function @multiTest@ tests the integer @n@ for primality @k@ times
21 using the random numbers supplied in @rs@; it returns a @Bool@ and a
22 depleted sequence of random numbers. When the @Bool@ is @True@ the
23 number is prime with probability $1 - 4^{-@k@}$; the number is
24 definitely composite when the @Bool@ is @False@.
25
26 > multiTest :: Int -> [Int] -> Integer -> (Bool, [Int])
27 > multiTest k rs n
28 >  = if n <= 1 || even n then (n==2, rs) else mTest k rs
29 >    where mTest 0 rs = (True, rs)
30 >          mTest k rs = if t then mTest (k-1) rs' else (False, rs')
31 >                           where (t, rs') = singleTest n (findKQ n) rs
32
33 The function @findKQ@ takes an odd integer $n$ and returns the tuple
34 $(k,q)$ such that $n = q2^k+1$.
35
36 > findKQ :: Integer -> (Integer, Integer)
37 > findKQ n = f (0, (n-1))
38 >            where f (k,q) = if r == 0 then f (k+1, d) else (k, q)
39 >                            where (d, r) = q `divMod` 2
40
41 The function @singleTest@ takes an odd, positive, @Integer@ @n@ and a
42 pair of @Integer@'s derived from @n@ by the @findKQ@ function
43
44 > singleTest :: Integer -> (Integer, Integer) -> [Int] -> (Bool, [Int])
45 > singleTest n kq rs
46 >  = (singleTestX n kq (2+x), rs')
47 >    where (x, rs')       = random (n-2) rs
48
49 Given a value @x@ we can test whether we have a witness to @n@'s
50 compositeness using @singleTestX@.
51
52 > singleTestX n (k, q) x
53 >  = t == 1 || t == n-1 || witness ts
54 >    where (t:ts)         = take (fromInteger k) (iterate square (powerMod x q n))
55 >          witness []     = False
56 >          witness (t:ts) = if t == n-1 then True       else
57 >                           if t == 1   then False      else
58 >                                            witness ts
59 >          square x       = (x*x) `mod` n
60
61 The @random@ function takes a number @n@ and a list of pseudo-random
62 numbers @rs@ and returns a tuple consisting of an @Integer@ $x$ in the
63 range $0 \leq x < @n@$, and the remainder of the pseudo-random
64 numbers.
65
66 > random :: Integer -> [Int] -> (Integer, [Int])
67 > random n rs = (makeNumber 65536 (uniform ns rs1), rs2)
68 >               where ns        = chop 65536 n
69 >                     (rs1,rs2) = splitAt (length ns) rs
70
71 The @uniform@ function generates a sequence of @Integer@'s such that,
72 when considered as a sequence of digits, we generate a number uniform
73 in the range @0..ns@ from the random numbers @rs@.
74
75 > uniform :: [Integer] -> [Int] -> [Integer]
76 > uniform [n]    [r]    = [toInteger r `mod` n]
77 > uniform (n:ns) (r:rs) = if t == n then t: uniform ns rs
78 >                                   else t: map ((`mod` 65536). toInteger) rs
79 >                         where t  = toInteger r `mod` (n+1)
80
81 This concludes the Prime testing algorithm.