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