ghc/lib/std/Random.lhs

   1
   2 This module implements a (good) random number generator.
   3
   4 The June 1988 (v31 #6) issue of the Communications of the ACM has an
   5 article by Pierre L'Ecuyer called, "Efficient and Portable Combined
   6 Random Number Generators".  Here is the Portable Combined Generator of
   7 L'Ecuyer for 32-bit computers.  It has a period of roughly 2.30584e18.
   8
   9 Transliterator: Lennart Augustsson
  10
  11 \begin{code}
  12 module Random
  13         (
  14          random,
  15          randomIO
  16         ) where
  17
  18 import CPUTime (getCPUTime)
  19 import Time (getClockTime, ClockTime(..))
  20
  21 randomIO :: (Integer, Integer) -> IO [Integer]
  22 randomIO lh = do
  23     ct          <- getCPUTime
  24     (TOD sec _) <- getClockTime
  25     return (random lh (toInteger sec * 12345 + ct))
  26
  27 random :: (Integer, Integer) -> Integer -> [Integer]
  28 random (l, h) s =
  29     if l > h then error "Random.random: Empty interval" else
  30     if s < 0 then random (l, h) (-s) else
  31         let (q, s1) = s `divMod` 2147483562
  32             s2 = q `mod` 2147483398
  33             k = h-l + 1
  34             b = 2147483561
  35             n = iLogBase b k
  36             f is = let (xs, is') = splitAt n is
  37                    in  foldr (\ i r -> fromInt i + r * b) 0 xs `mod` k + l : f is'
  38         in  f (randomInts (toInt (s1+1)) (toInt (s2+1)))
  39
  40 iLogBase b i = if i < b then 1 else 1 + iLogBase b (i `div` b)
  41
  42 -- Use seeds s1 in 1..2147483562 and s2 in 1..2147483398 to generate
  43 -- an infinite list of random Ints.
  44 randomInts :: Int -> Int -> [Int]
  45 randomInts s1 s2 =
  46     if 1 <= s1 && s1 <= 2147483562 then
  47         if 1 <= s2 && s2 <= 2147483398 then
  48             rands s1 s2
  49         else
  50             error "randomInts: Bad second seed."
  51     else
  52         error "randomInts: Bad first seed."
  53
  54 rands :: Int -> Int -> [Int]
  55 rands s1 s2 = z' : rands s1'' s2''
  56         where   z'   = if z < 1 then z + 2147483562 else z
  57                 z    = s1'' - s2''
  58
  59                 k    = s1 `quot` 53668
  60                 s1'  = 40014 * (s1 - k * 53668) - k * 12211
  61                 s1'' = if s1' < 0 then s1' + 2147483563 else s1'
  62
  63                 k'   = s2 `quot` 52774
  64                 s2'  = 40692 * (s2 - k' * 52774) - k' * 3791
  65                 s2'' = if s2' < 0 then s2' + 2147483399 else s2'
  66
  67 \end{code}