ghc/lib/std/Random.lhs

   1 %
   2 % (c) The GRASP/AQUA Project, Glasgow University, 1995-99
   3 %
   4
   5
   6 The June 1988 (v31 #6) issue of the Communications of the ACM has an
   7 article by Pierre L'Ecuyer called, "Efficient and Portable Combined
   8 Random Number Generators".  Here is the Portable Combined Generator of
   9 L'Ecuyer for 32-bit computers.  It has a period of roughly 2.30584e18.
  10
  11 Transliterator: Lennart Augustsson
  12
  13 sof 1/99 - code brought (kicking and screaming) into the new Random
  14 world..
  15
  16 \begin{code}
  17 module Random
  18         (
  19           RandomGen(next, split)
  20         , StdGen
  21         , mkStdGen
  22         , Random ( random,   randomR,
  23                    randoms,  randomRs,
  24                    randomIO, randomRIO )
  25         , getStdRandom
  26         , getStdGen
  27         , setStdGen
  28         , newStdGen
  29         ) where
  30
  31 import CPUTime (getCPUTime)
  32 import PrelST
  33 import PrelRead
  34 import PrelIOBase
  35 import PrelNumExtra ( float2Double, double2Float )
  36 import PrelBase
  37 import PrelArr
  38 import Char ( isSpace, chr, ord )
  39 import Time (getClockTime, ClockTime(..))
  40
  41 \end{code}
  42
  43 \begin{code}
  44 class RandomGen g where
  45    next  :: g -> (Int, g)
  46    split :: g -> (g, g)
  47
  48 \end{code}
  49
  50 \begin{code}
  51 data StdGen
  52  = StdGen Int Int
  53
  54 instance RandomGen StdGen where
  55   next  = rand1
  56   split = splitStdGen
  57
  58 instance Show StdGen where
  59   showsPrec p (StdGen s1 s2) =
  60      showSignedInt p s1 .
  61      showSpace          .
  62      showSignedInt p s2
  63
  64 instance Read StdGen where
  65   readsPrec p = \ r ->
  66      case try_read r of
  67        r@[_] -> r
  68        _   -> [(unsafePerformIO mkStdRNG,r)] -- because it shouldn't ever fail.
  69     where
  70       try_read r = do
  71          (s1, r1) <- readDec (dropWhile isSpace r)
  72          (s2, r2) <- readDec (dropWhile isSpace r1)
  73          return (StdGen s1 s2, r2)
  74
  75 \end{code}
  76
  77 \begin{code}
  78 mkStdGen :: Int -> StdGen -- why not Integer ?
  79 mkStdGen s
  80  | s < 0     = mkStdGen (-s)
  81  | otherwise = StdGen (s1+1) (s2+1)
  82       where
  83         (q, s1) = s `divMod` 2147483562
  84         s2      = q `mod` 2147483398
  85
  86 createStdGen :: Integer -> StdGen
  87 createStdGen s
  88  | s < 0     = createStdGen (-s)
  89  | otherwise = StdGen (toInt (s1+1)) (toInt (s2+1))
  90       where
  91         (q, s1) = s `divMod` 2147483562
  92         s2      = q `mod` 2147483398
  93
  94 \end{code}
  95
  96 \begin{code}
  97
  98 -- Q: do all of these merit class membership?
  99 class Random a where
 100   randomR :: RandomGen g => (a,a) -> g -> (a,g)
 101   random  :: RandomGen g => g -> (a, g)
 102
 103   randomRs :: RandomGen g => (a,a) -> g -> [a]
 104   randoms  :: RandomGen g => g -> [a]
 105
 106   randomRIO :: (a,a) -> IO a
 107   randomIO  :: IO a
 108
 109   randoms  g      = x : randoms g' where (x,g') = random g
 110   randomRs ival g = x : randomRs ival g' where (x,g') = randomR ival g
 111
 112   randomIO         = getStdRandom random
 113   randomRIO range  = getStdRandom (randomR range)
 114
 115 \end{code}
 116
 117 \begin{code}
 118 instance Random Int where
 119   randomR (a,b) g = randomIvalInteger (toInteger a, toInteger b) g
 120   random g        = randomR (minBound,maxBound) g
 121
 122 instance Random Char where
 123   randomR (a,b) g =
 124       case (randomIvalInteger (toInteger (ord a), toInteger (ord b)) g) of
 125         (x,g) -> (chr x, g)
 126   random g        = randomR (minBound,maxBound) g
 127
 128 instance Random Bool where
 129   randomR (a,b) g =
 130       case (randomIvalInteger (toInteger (bool2Int a), toInteger (bool2Int b)) g) of
 131         (x, g) -> (int2Bool x, g)
 132        where
 133          bool2Int False = 0
 134          bool2Int True  = 1
 135
 136          int2Bool 0     = False
 137          int2Bool _     = True
 138
 139   random g        = randomR (minBound,maxBound) g
 140
 141 instance Random Integer where
 142   randomR ival g = randomIvalInteger ival g
 143   random g       = randomR (toInteger (minBound::Int), toInteger (maxBound::Int)) g
 144
 145 instance Random Double where
 146   randomR ival g = randomIvalDouble ival id g
 147   random g       = randomR (0::Double,1) g
 148
 149 -- hah, so you thought you were saving cycles by using Float?
 150 instance Random Float where
 151   randomR (a,b) g = randomIvalDouble (float2Double a, float2Double b) double2Float g
 152   random g        = randomIvalDouble (0::Double,1) double2Float g
 153
 154 \end{code}
 155
 156
 157 \begin{code}
 158 mkStdRNG :: IO StdGen
 159 mkStdRNG = do
 160     ct          <- getCPUTime
 161     (TOD sec _) <- getClockTime
 162     return (createStdGen (sec * 12345 + ct))
 163
 164 randomIvalInteger :: (RandomGen g, Num a) => (Integer, Integer) -> g -> (a, g)
 165 randomIvalInteger (l,h) rng
 166  | l > h     = randomIvalInteger (h,l) rng
 167  | otherwise = case (f n 1 rng) of (v, rng') -> (fromInteger (l + v `mod` (k+1)), rng')
 168      where
 169        k = h - l + 1
 170        b = 2147483561
 171        n = iLogBase b k
 172
 173        f 0 acc g = (acc, g)
 174        f n acc g =
 175           let
 176            (x,g')   = next g
 177           in
 178           f (n-1) (fromInt x + acc * b) g'
 179
 180 randomIvalDouble :: (RandomGen g, Fractional a) => (Double, Double) -> (Double -> a) -> g -> (a, g)
 181 randomIvalDouble (l,h) fromDouble rng
 182   | l > h     = randomIvalDouble (h,l) fromDouble rng
 183   | otherwise =
 184        case (randomIvalInteger (toInteger (minBound::Int), toInteger (maxBound::Int)) rng) of
 185          (x, rng') ->
 186             let
 187              scaled_x =
 188                 fromDouble l +
 189                 fromDouble (h-l) *
 190                  (fromIntegral (x::Int) * 4.6566130638969828e-10)
 191                   -- magic number stolen from old HBC code (Random.randomDoubles.)
 192             in
 193             (scaled_x, rng')
 194
 195 iLogBase :: Integer -> Integer -> Integer
 196 iLogBase b i = if i < b then 1 else 1 + iLogBase b (i `div` b)
 197
 198 rand1 :: StdGen -> (Int, StdGen)
 199 rand1 (StdGen s1 s2) = (z', StdGen s1'' s2'')
 200         where   z'   = if z < 1 then z + 2147483562 else z
 201                 z    = s1'' - s2''
 202
 203                 k    = s1 `quot` 53668
 204                 s1'  = 40014 * (s1 - k * 53668) - k * 12211
 205                 s1'' = if s1' < 0 then s1' + 2147483563 else s1'
 206
 207                 k'   = s2 `quot` 52774
 208                 s2'  = 40692 * (s2 - k' * 52774) - k' * 3791
 209                 s2'' = if s2' < 0 then s2' + 2147483399 else s2'
 210
 211 splitStdGen :: StdGen -> (StdGen, StdGen)
 212 splitStdGen std@(StdGen s1 s2) = (std, StdGen new_s1 new_s2)
 213    where
 214        -- simple in the extreme..
 215       new_s1
 216         | s1 == 2147483562 = 1
 217         | otherwise        = s1 + 1
 218
 219       new_s2
 220         | s2 == 1          = 2147483398
 221         | otherwise        = s2 - 1
 222
 223
 224
 225 \end{code}
 226
 227
 228 \begin{code}
 229 global_rng :: MutableVar RealWorld StdGen
 230 global_rng = unsafePerformIO $ do
 231    rng <- mkStdRNG
 232    stToIO (newVar rng)
 233
 234 setStdGen :: StdGen -> IO ()
 235 setStdGen sgen = stToIO (writeVar global_rng sgen)
 236
 237 getStdGen :: IO StdGen
 238 getStdGen = stToIO (readVar global_rng)
 239
 240 newStdGen :: IO StdGen
 241 newStdGen = do
 242   rng <- getStdGen
 243   let (a,b) = split rng
 244   setStdGen a
 245   return b
 246
 247 getStdRandom :: (StdGen -> (a,StdGen)) -> IO a
 248 getStdRandom f = do
 249    rng          <- getStdGen
 250    let (v, new_rng) = f rng
 251    setStdGen new_rng
 252    return v
 253 \end{code}