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