System/Random.hs

   1 -----------------------------------------------------------------------------
   2 -- |
   3 -- Module      :  System.Random
   4 -- Copyright   :  (c) The University of Glasgow 2001
   5 -- License     :  BSD-style (see the file libraries/base/LICENSE)
   6 --
   7 -- Maintainer  :  libraries@haskell.org
   8 -- Stability   :  provisional
   9 -- Portability :  portable
  10 --
  11 -- Random numbers.
  12 --
  13 -----------------------------------------------------------------------------
  14
  15 module System.Random
  16         (
  17
  18         -- $intro
  19
  20         -- * The 'RandomGen' class, and the 'StdGen' generator
  21
  22           RandomGen(next, split, genRange)
  23         , StdGen
  24         , mkStdGen
  25
  26         -- * The 'Random' class
  27         , Random ( random,   randomR,
  28                    randoms,  randomRs,
  29                    randomIO, randomRIO )
  30
  31         -- * The global random number generator
  32
  33         -- $globalrng
  34
  35         , getStdRandom
  36         , getStdGen
  37         , setStdGen
  38         , newStdGen
  39
  40         -- * References
  41         -- $references
  42
  43         ) where
  44
  45 import Prelude
  46
  47 #ifdef __NHC__
  48 import CPUTime          ( getCPUTime )
  49 #else
  50 import System.CPUTime   ( getCPUTime )
  51 import System.Time      ( getClockTime, ClockTime(..) )
  52 #endif
  53 import Data.Char        ( isSpace, chr, ord )
  54 import System.IO.Unsafe ( unsafePerformIO )
  55 import Data.IORef
  56 import Numeric          ( readDec )
  57
  58 #ifdef __GLASGOW_HASKELL__
  59 import GHC.IOBase       ( stToIO )
  60 #endif
  61
  62 -- The standard nhc98 implementation of Time.ClockTime does not match
  63 -- the extended one expected in this module, so we lash-up a quick
  64 -- replacement here.
  65 #ifdef __NHC__
  66 data ClockTime = TOD Integer ()
  67 foreign import ccall "time.h time" readtime :: Int -> IO Int
  68 getClockTime :: IO ClockTime
  69 getClockTime = do t <- readtime 0;  return (TOD (toInteger t) ())
  70 #endif
  71
  72 {- $intro
  73
  74 The "Random" library deals with the common task of pseudo-random
  75 number generation. The library makes it possible to generate
  76 repeatable results, by starting with a specified initial random
  77 number generator; or to get different results on each run by using the
  78 system-initialised generator, or by supplying a seed from some other
  79 source.
  80
  81 The library is split into two layers:
  82
  83 * A core /random number generator/ provides a supply of bits. The class
  84 'RandomGen' provides a common interface to such generators.
  85
  86 * The class 'Random' provides a way to extract particular values from
  87 a random number generator. For example, the 'Float' instance of 'Random'
  88 allows one to generate random values of type 'Float'.
  89
  90 [Comment found in this file when merging with Library Report:]
  91
  92 The June 1988 (v31 \#6) issue of the Communications of the ACM has an
  93 article by Pierre L'Ecuyer called, /Efficient and Portable Combined
  94 Random Number Generators/.  Here is the Portable Combined Generator of
  95 L'Ecuyer for 32-bit computers.  It has a period of roughly 2.30584e18.
  96
  97 Transliterator: Lennart Augustsson
  98
  99 -}
 100
 101 -- |RandomGen
 102 -- The class 'RandomGen' provides a common interface to random number generators.
 103
 104 class RandomGen g where
 105
 106    -- |The 'next' operation allows one to extract at least 30 bits (one 'Int''s
 107    -- worth) from the generator, returning a new generator as well.  The
 108    -- integer returned may be positive or negative.
 109    next     :: g -> (Int, g)
 110
 111    -- |The 'split' operation allows one to obtain two distinct random number
 112    -- generators. This is very useful in functional programs (for example, when
 113    -- passing a random number generator down to recursive calls), but very
 114    -- little work has been done on statistically robust implementations of
 115    -- @split ([1,4]@ are the only examples we know of).
 116    split    :: g -> (g, g)
 117
 118    genRange :: g -> (Int,Int)
 119
 120    -- default mathod
 121    genRange g = (minBound,maxBound)
 122
 123 {- |The "Random" library provides one instance of 'RandomGen', the abstract data
 124 type 'StdGen'.
 125
 126 The result of repeatedly using next should be at least as statistically robust
 127 as the /Minimal Standard Random Number Generator/ described by
 128 ["System.Random\#Park", "System.Random\#Carta"]. Until more
 129 is known about implementations of 'split', all we require is that 'split' deliver
 130 generators that are (a) not identical and (b) independently robust in the sense
 131 just given.
 132
 133 The 'show'\/'Read' instances of 'StdGen' provide a primitive way to save the
 134 state of a random number generator. It is required that @read (show g) == g@.
 135
 136 In addition, 'read' may be used to map an arbitrary string (not necessarily one
 137 produced by 'show') onto a value of type 'StdGen'. In general, the 'read'
 138 instance of 'StdGen' has the following properties:
 139
 140 * It guarantees to succeed on any string.
 141
 142 *It guarantees to consume only a finite portion of the string.
 143
 144 * Different argument strings are likely to result in different results.
 145
 146 The function 'mkStdGen' provides an alternative way of producing an initial
 147 generator, by mapping an 'Int' into a generator. Again, distinct arguments
 148 should be likely to produce distinct generators.
 149
 150 Programmers may, of course, supply their own instances of 'RandomGen'.
 151
 152 -}
 153
 154 data StdGen
 155  = StdGen Int Int
 156
 157 instance RandomGen StdGen where
 158   next  = stdNext
 159   split = stdSplit
 160
 161 instance Show StdGen where
 162   showsPrec p (StdGen s1 s2) =
 163      showsPrec p s1 .
 164      showChar ' ' .
 165      showsPrec p s2
 166
 167 instance Read StdGen where
 168   readsPrec _p = \ r ->
 169      case try_read r of
 170        r@[_] -> r
 171        _   -> [stdFromString r] -- because it shouldn't ever fail.
 172     where
 173       try_read r = do
 174          (s1, r1) <- readDec (dropWhile isSpace r)
 175          (s2, r2) <- readDec (dropWhile isSpace r1)
 176          return (StdGen s1 s2, r2)
 177
 178 {-
 179  If we cannot unravel the StdGen from a string, create
 180  one based on the string given.
 181 -}
 182 stdFromString         :: String -> (StdGen, String)
 183 stdFromString s        = (mkStdGen num, rest)
 184         where (cs, rest) = splitAt 6 s
 185               num        = foldl (\a x -> x + 3 * a) 1 (map ord cs)
 186
 187
 188 mkStdGen :: Int -> StdGen -- why not Integer ?
 189 mkStdGen s
 190  | s < 0     = mkStdGen (-s)
 191  | otherwise = StdGen (s1+1) (s2+1)
 192       where
 193         (q, s1) = s `divMod` 2147483562
 194         s2      = q `mod` 2147483398
 195
 196 createStdGen :: Integer -> StdGen
 197 createStdGen s
 198  | s < 0     = createStdGen (-s)
 199  | otherwise = StdGen (fromInteger (s1+1)) (fromInteger (s2+1))
 200       where
 201         (q, s1) = s `divMod` 2147483562
 202         s2      = q `mod` 2147483398
 203
 204 -- FIXME: 1/2/3 below should be ** (vs@30082002) XXX
 205
 206 {- |The 'Random' class
 207 With a source of random number supply in hand, the 'Random' class allows the
 208 programmer to extract random values of a variety of types.
 209
 210 * 'randomR' takes a range /(lo,hi)/ and a random number generator /g/, and returns
 211 a random value uniformly distributed in the closed interval /[lo,hi]/, together
 212 with a new generator. It is unspecified what happens if /lo>hi/. For continuous
 213 types there is no requirement that the values /lo/ and /hi/ are ever produced,
 214 but they may be, depending on the implementation and the interval.
 215
 216 * 'random' does the same as 'randomR', but does not take a range.
 217
 218 (1) For bounded types (instances of 'Bounded', such as 'Char'), the range is
 219 normally the whole type.
 220
 221 (2) For fractional types, the range is normally the semi-closed interval @[0,1)@.
 222
 223 (3) For 'Integer', the range is (arbitrarily) the range of 'Int'.
 224
 225 * The plural versions, 'randomRs' and 'randoms', produce an infinite list of
 226 random values, and do not return a new generator.
 227
 228 * The 'IO' versions, 'randomRIO' and 'randomIO', use the global random number
 229 generator (see Section 17.3
 230 <http://www.haskell.org/onlinelibrary/random.html#global-rng>).
 231 -}
 232
 233 class Random a where
 234   -- |Minimal complete definition: 'random' and 'randomR'
 235   random  :: RandomGen g => g -> (a, g)
 236   randomR :: RandomGen g => (a,a) -> g -> (a,g)
 237
 238   -- |Default methods
 239   randoms  :: RandomGen g => g -> [a]
 240   randoms  g      = (\(x,g') -> x : randoms g') (random g)
 241
 242   randomRs :: RandomGen g => (a,a) -> g -> [a]
 243   randomRs ival g = x : randomRs ival g' where (x,g') = randomR ival g
 244
 245   randomIO  :: IO a
 246   randomIO         = getStdRandom random
 247
 248   randomRIO :: (a,a) -> IO a
 249   randomRIO range  = getStdRandom (randomR range)
 250
 251
 252 instance Random Int where
 253   randomR (a,b) g = randomIvalInteger (toInteger a, toInteger b) g
 254   random g        = randomR (minBound,maxBound) g
 255
 256 instance Random Char where
 257   randomR (a,b) g =
 258       case (randomIvalInteger (toInteger (ord a), toInteger (ord b)) g) of
 259         (x,g) -> (chr x, g)
 260   random g        = randomR (minBound,maxBound) g
 261
 262 instance Random Bool where
 263   randomR (a,b) g =
 264       case (randomIvalInteger (toInteger (bool2Int a), toInteger (bool2Int b)) g) of
 265         (x, g) -> (int2Bool x, g)
 266        where
 267          bool2Int False = 0
 268          bool2Int True  = 1
 269
 270          int2Bool 0     = False
 271          int2Bool _     = True
 272
 273   random g        = randomR (minBound,maxBound) g
 274
 275 instance Random Integer where
 276   randomR ival g = randomIvalInteger ival g
 277   random g       = randomR (toInteger (minBound::Int), toInteger (maxBound::Int)) g
 278
 279 instance Random Double where
 280   randomR ival g = randomIvalDouble ival id g
 281   random g       = randomR (0::Double,1) g
 282
 283 -- hah, so you thought you were saving cycles by using Float?
 284 instance Random Float where
 285   random g        = randomIvalDouble (0::Double,1) realToFrac g
 286   randomR (a,b) g = randomIvalDouble (realToFrac a, realToFrac b) realToFrac g
 287
 288 mkStdRNG :: Integer -> IO StdGen
 289 mkStdRNG o = do
 290     ct          <- getCPUTime
 291     (TOD sec _) <- getClockTime
 292     return (createStdGen (sec * 12345 + ct + o))
 293
 294 randomIvalInteger :: (RandomGen g, Num a) => (Integer, Integer) -> g -> (a, g)
 295 randomIvalInteger (l,h) rng
 296  | l > h     = randomIvalInteger (h,l) rng
 297  | otherwise = case (f n 1 rng) of (v, rng') -> (fromInteger (l + v `mod` k), rng')
 298      where
 299        k = h - l + 1
 300        b = 2147483561
 301        n = iLogBase b k
 302
 303        f 0 acc g = (acc, g)
 304        f n acc g =
 305           let
 306            (x,g')   = next g
 307           in
 308           f (n-1) (fromIntegral x + acc * b) g'
 309
 310 randomIvalDouble :: (RandomGen g, Fractional a) => (Double, Double) -> (Double -> a) -> g -> (a, g)
 311 randomIvalDouble (l,h) fromDouble rng
 312   | l > h     = randomIvalDouble (h,l) fromDouble rng
 313   | otherwise =
 314        case (randomIvalInteger (toInteger (minBound::Int), toInteger (maxBound::Int)) rng) of
 315          (x, rng') ->
 316             let
 317              scaled_x =
 318                 fromDouble ((l+h)/2) +
 319                 fromDouble ((h-l) / realToFrac intRange) *
 320                 fromIntegral (x::Int)
 321             in
 322             (scaled_x, rng')
 323
 324 intRange :: Integer
 325 intRange  = toInteger (maxBound::Int) - toInteger (minBound::Int)
 326
 327 iLogBase :: Integer -> Integer -> Integer
 328 iLogBase b i = if i < b then 1 else 1 + iLogBase b (i `div` b)
 329
 330 stdNext :: StdGen -> (Int, StdGen)
 331 stdNext (StdGen s1 s2) = (z', StdGen s1'' s2'')
 332         where   z'   = if z < 1 then z + 2147483562 else z
 333                 z    = s1'' - s2''
 334
 335                 k    = s1 `quot` 53668
 336                 s1'  = 40014 * (s1 - k * 53668) - k * 12211
 337                 s1'' = if s1' < 0 then s1' + 2147483563 else s1'
 338
 339                 k'   = s2 `quot` 52774
 340                 s2'  = 40692 * (s2 - k' * 52774) - k' * 3791
 341                 s2'' = if s2' < 0 then s2' + 2147483399 else s2'
 342
 343 stdSplit            :: StdGen -> (StdGen, StdGen)
 344 stdSplit std@(StdGen s1 s2)
 345                      = (left, right)
 346                        where
 347                         -- no statistical foundation for this!
 348                         left    = StdGen new_s1 t2
 349                         right   = StdGen t1 new_s2
 350
 351                         new_s1 | s1 == 2147483562 = 1
 352                                | otherwise        = s1 + 1
 353
 354                         new_s2 | s2 == 1          = 2147483398
 355                                | otherwise        = s2 - 1
 356
 357                         StdGen t1 t2 = snd (next std)
 358
 359 -- The global random number generator
 360
 361 {- $globalrng
 362
 363 There is a single, implicit, global random number generator of type
 364 'StdGen', held in some global variable maintained by the 'IO' monad. It is
 365 initialised automatically in some system-dependent fashion, for example, by
 366 using the time of day, or Linux's kernel random number generator. To get
 367 deterministic behaviour, use 'setStdGen'.
 368 -}
 369
 370 -- |'setStdGen' sets the global random number generator.
 371 setStdGen :: StdGen -> IO ()
 372 setStdGen sgen = writeIORef theStdGen sgen
 373
 374 -- |'getStdGen' gets the global random number generator.
 375 getStdGen :: IO StdGen
 376 getStdGen  = readIORef theStdGen
 377
 378 -- |'newStdGen' applies 'split' to the current global random generator, updates it
 379 -- with one of the results, and returns the other.
 380 theStdGen :: IORef StdGen
 381 theStdGen  = unsafePerformIO $ do
 382    rng <- mkStdRNG 0
 383    newIORef rng
 384
 385 newStdGen :: IO StdGen
 386 newStdGen = do
 387   rng <- getStdGen
 388   let (a,b) = split rng
 389   setStdGen a
 390   return b
 391
 392 {- |'getStdRandom' uses the supplied function to get a value from the current
 393 global random generator, and updates the global generator with the new generator
 394 returned by the function. For example, 'rollDice' gets a random integer between 1 and 6:
 395
 396 >  rollDice :: IO Int
 397 >  rollDice = getStdRandom (randomR (1,6))
 398
 399 -}
 400
 401 getStdRandom :: (StdGen -> (a,StdGen)) -> IO a
 402 getStdRandom f = do
 403    rng          <- getStdGen
 404    let (v, new_rng) = f rng
 405    setStdGen new_rng
 406    return v
 407
 408 {- $references
 409
 410 * [1] FW Burton and RL Page, /Distributed random number generation/,
 411 Journal of Functional Programming, 2(2):203-212, April 1992.
 412
 413 * [2] SK #Park# Park, and KW Miller, /Random number generators -
 414 good ones are hard to find/, Comm ACM 31(10), Oct 1988, pp1192-1201.
 415
 416 * [3] DG #Carta# Carta, /Two fast implementations of the minimal standard
 417 random number generator/, Comm ACM, 33(1), Jan 1990, pp87-88.
 418
 419 * [4] P Hellekalek, /Don\'t trust parallel Monte Carlo/,
 420 Department of Mathematics, University of Salzburg,
 421 <http://random.mat.sbg.ac.at/~peter/pads98.ps>, 1998.
 422
 423 The Web site <http://random.mat.sbg.ac.at/> is a great source of information.
 424
 425 -}