ghc/lib/hbc/Random.hs

   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
  12 module Random(randomInts, randomDoubles, normalRandomDoubles) where
  13 -- Use seeds s1 in 1..2147483562 and s2 in 1..2147483398 to generate
  14 -- an infinite list of random Ints.
  15 randomInts :: Int -> Int -> [Int]
  16 randomInts s1 s2 =
  17     if 1 <= s1 && s1 <= 2147483562 then
  18         if 1 <= s2 && s2 <= 2147483398 then
  19             rands s1 s2
  20         else
  21             error "randomInts: Bad second seed."
  22     else
  23         error "randomInts: Bad first seed."
  24
  25 rands :: Int -> Int -> [Int]
  26 rands s1 s2 = z' : rands s1'' s2''
  27         where   z'   = if z < 1 then z + 2147483562 else z
  28                 z    = s1'' - s2''
  29
  30                 k    = s1 `quot` 53668
  31                 s1'  = 40014 * (s1 - k * 53668) - k * 12211
  32                 s1'' = if s1' < 0 then s1' + 2147483563 else s1'
  33
  34                 k'   = s2 `quot` 52774
  35                 s2'  = 40692 * (s2 - k' * 52774) - k' * 3791
  36                 s2'' = if s2' < 0 then s2' + 2147483399 else s2'
  37
  38 -- Same values for s1 and s2 as above, generates an infinite
  39 -- list of Doubles uniformly distibuted in (0,1).
  40 randomDoubles :: Int -> Int -> [Double]
  41 randomDoubles s1 s2 = map (\x -> fromIntegral x * 4.6566130638969828e-10) (randomInts s1 s2)
  42
  43 -- The normal distribution stuff is stolen from Tim Lambert's
  44 -- M*****a version
  45
  46 -- normalRandomDoubles is given two seeds and returns an infinite list of random
  47 -- normal variates with mean 0 and variance 1.  (Box Muller method see
  48 -- "Art of Computer Programming Vol 2")
  49 normalRandomDoubles :: Int -> Int -> [Double]
  50 normalRandomDoubles s1 s2 = boxMuller (map (\x->2*x-1) (randomDoubles s1 s2))
  51
  52 -- boxMuller takes a stream of uniform random numbers on [-1,1] and
  53 -- returns a stream of normally distributed random numbers.
  54 boxMuller :: [Double] -> [Double]
  55 boxMuller (x1:x2:xs) | r <= 1    = x1*m : x2*m : rest
  56                      | otherwise = rest
  57                                    where r = x1*x1 + x2*x2
  58                                          m = sqrt(-2*log r/r)
  59                                          rest = boxMuller xs