1 
    2 
    3 
    4 
    5 
    6 
    7 
    8 
    9 
   10 
   11 
   12 
   13 
   14 
   15 
   16 
   17 
   18 
   19 
   20   module Prime (multiTest) where
   21   import IntLib
   22   import MyRandom
   23 
   24 
   25 
   26 
   27 
   28 
   29 
   30   multiTest :: Int -> [Int] -> Integer -> (Bool, [Int])
   31   multiTest k rs n
   32    = if n <= 1 || even n then (n==2, rs) else mTest k rs
   33      where mTest 0 rs = (True, rs)
   34            mTest k rs = if t then mTest (k-1) rs' else (False, rs')
   35                             where (t, rs') = singleTest n (findKQ n) rs
   36 
   37 
   38 
   39 
   40   findKQ :: Integer -> (Integer, Integer)
   41   findKQ n = f (0, (n-1))
   42              where f (k,q) = if r == 0 then f (k+1, d) else (k, q)
   43                              where (d, r) = q `divMod` 2
   44 
   45 
   46 
   47 
   48   singleTest :: Integer -> (Integer, Integer) -> [Int] -> (Bool, [Int])
   49   singleTest n kq rs
   50    = (singleTestX n kq (2+x), rs')
   51      where (x, rs')       = random (n-2) rs
   52 
   53 
   54 
   55 
   56   singleTestX n (k, q) x
   57    = t == 1 || t == n-1 || witness ts
   58      where (t:ts)         = take (fromInteger k) (iterate square (powerMod x q n))
   59            witness []     = False
   60            witness (t:ts) = if t == n-1 then True       else
   61                             if t == 1   then False      else
   62                                              witness ts
   63            square x       = (x*x) `mod` n
   64 
   65 
   66 
   67 
   68 
   69 
   70   random :: Integer -> [Int] -> (Integer, [Int])
   71   random n rs = (makeNumber 65536 (uniform ns rs1), rs2)
   72                 where ns        = chop 65536 n
   73                       (rs1,rs2) = splitAt (length ns) rs
   74 
   75 
   76 
   77 
   78 
   79   uniform :: [Integer] -> [Int] -> [Integer]
   80   uniform [n]    [r]    = [toInteger r `mod` n]
   81   uniform (n:ns) (r:rs) = if t == n then t: uniform ns rs
   82                                     else t: map ((`mod` 65536). toInteger) rs
   83                           where t  = toInteger r `mod` (n+1)
   84 
   85