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   genRange _ = stdRange
 166
 167 instance Show StdGen where
 168   showsPrec p (StdGen s1 s2) =
 169      showsPrec p s1 .
 170      showChar ' ' .
 171      showsPrec p s2
 172
 173 instance Read StdGen where
 174   readsPrec _p = \ r ->
 175      case try_read r of
 176        r@[_] -> r
 177        _   -> [stdFromString r] -- because it shouldn't ever fail.
 178     where
 179       try_read r = do
 180          (s1, r1) <- readDec (dropWhile isSpace r)
 181          (s2, r2) <- readDec (dropWhile isSpace r1)
 182          return (StdGen s1 s2, r2)
 183
 184 {-
 185  If we cannot unravel the StdGen from a string, create
 186  one based on the string given.
 187 -}
 188 stdFromString         :: String -> (StdGen, String)
 189 stdFromString s        = (mkStdGen num, rest)
 190         where (cs, rest) = splitAt 6 s
 191               num        = foldl (\a x -> x + 3 * a) 1 (map ord cs)
 192
 193
 194 {- |
 195 The function 'mkStdGen' provides an alternative way of producing an initial
 196 generator, by mapping an 'Int' into a generator. Again, distinct arguments
 197 should be likely to produce distinct generators.
 198
 199 Programmers may, of course, supply their own instances of 'RandomGen'.
 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 -}