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