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