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 #else
  74 import System.CPUTime   ( getCPUTime )
  75 import System.Time      ( getClockTime, ClockTime(..) )
  76 #endif
  77 import Data.Char        ( isSpace, chr, ord )
  78 import System.IO.Unsafe ( unsafePerformIO )
  79 import Data.IORef
  80 import Numeric          ( readDec )
  81
  82 -- The standard nhc98 implementation of Time.ClockTime does not match
  83 -- the extended one expected in this module, so we lash-up a quick
  84 -- replacement here.
  85 #ifdef __NHC__
  86 data ClockTime = TOD Integer ()
  87 foreign import ccall "time.h time" readtime :: Ptr () -> IO Int
  88 getClockTime :: IO ClockTime
  89 getClockTime = do t <- readtime nullPtr;  return (TOD (toInteger t) ())
  90 #endif
  91
  92 -- | The class 'RandomGen' provides a common interface to random number
  93 -- generators.
  94
  95 class RandomGen g where
  96
  97    -- |The 'next' operation returns an 'Int' that is uniformly distributed
  98    -- in the range returned by 'genRange' (including both end points),
  99    -- and a new generator.
 100    next     :: g -> (Int, g)
 101
 102    -- |The 'split' operation allows one to obtain two distinct random number
 103    -- generators. This is very useful in functional programs (for example, when
 104    -- passing a random number generator down to recursive calls), but very
 105    -- little work has been done on statistically robust implementations of
 106    -- 'split' (["System.Random\#Burton", "System.Random\#Hellekalek"]
 107    -- are the only examples we know of).
 108    split    :: g -> (g, g)
 109
 110    -- |The 'genRange' operation yields the range of values returned by
 111    -- the generator.
 112    --
 113    -- It is required that:
 114    --
 115    -- * If @(a,b) = 'genRange' g@, then @a < b@.
 116    --
 117    -- * 'genRange' is not strict.
 118    --
 119    -- The second condition ensures that 'genRange' cannot examine its
 120    -- argument, and hence the value it returns can be determined only by the
 121    -- instance of 'RandomGen'.  That in turn allows an implementation to make
 122    -- a single call to 'genRange' to establish a generator's range, without
 123    -- being concerned that the generator returned by (say) 'next' might have
 124    -- a different range to the generator passed to 'next'.
 125    genRange :: g -> (Int,Int)
 126
 127    -- default method
 128    genRange g = (minBound,maxBound)
 129
 130 {- |
 131 The 'StdGen' instance of 'RandomGen' has a 'genRange' of at least 30 bits.
 132
 133 The result of repeatedly using 'next' should be at least as statistically
 134 robust as the /Minimal Standard Random Number Generator/ described by
 135 ["System.Random\#Park", "System.Random\#Carta"].
 136 Until more is known about implementations of 'split', all we require is
 137 that 'split' deliver generators that are (a) not identical and
 138 (b) independently robust in the sense just given.
 139
 140 The 'Show' and 'Read' instances of 'StdGen' provide a primitive way to save the
 141 state of a random number generator.
 142 It is required that @'read' ('show' g) == g@.
 143
 144 In addition, 'read' may be used to map an arbitrary string (not necessarily one
 145 produced by 'show') onto a value of type 'StdGen'. In general, the 'read'
 146 instance of 'StdGen' has the following properties:
 147
 148 * It guarantees to succeed on any string.
 149
 150 * It guarantees to consume only a finite portion of the string.
 151
 152 * Different argument strings are likely to result in different results.
 153
 154 -}
 155
 156 data StdGen
 157  = StdGen Int Int
 158
 159 instance RandomGen StdGen where
 160   next  = stdNext
 161   split = stdSplit
 162   genRange _ = stdRange
 163
 164 instance Show StdGen where
 165   showsPrec p (StdGen s1 s2) =
 166      showsPrec p s1 .
 167      showChar ' ' .
 168      showsPrec p s2
 169
 170 instance Read StdGen where
 171   readsPrec _p = \ r ->
 172      case try_read r of
 173        r@[_] -> r
 174        _   -> [stdFromString r] -- because it shouldn't ever fail.
 175     where
 176       try_read r = do
 177          (s1, r1) <- readDec (dropWhile isSpace r)
 178          (s2, r2) <- readDec (dropWhile isSpace r1)
 179          return (StdGen s1 s2, r2)
 180
 181 {-
 182  If we cannot unravel the StdGen from a string, create
 183  one based on the string given.
 184 -}
 185 stdFromString         :: String -> (StdGen, String)
 186 stdFromString s        = (mkStdGen num, rest)
 187         where (cs, rest) = splitAt 6 s
 188               num        = foldl (\a x -> x + 3 * a) 1 (map ord cs)
 189
 190
 191 {- |
 192 The function 'mkStdGen' provides an alternative way of producing an initial
 193 generator, by mapping an 'Int' into a generator. Again, distinct arguments
 194 should be likely to produce distinct generators.
 195 -}
 196 mkStdGen :: Int -> StdGen -- why not Integer ?
 197 mkStdGen s
 198  | s < 0     = mkStdGen (-s)
 199  | otherwise = StdGen (s1+1) (s2+1)
 200       where
 201         (q, s1) = s `divMod` 2147483562
 202         s2      = q `mod` 2147483398
 203
 204 createStdGen :: Integer -> StdGen
 205 createStdGen s
 206  | s < 0     = createStdGen (-s)
 207  | otherwise = StdGen (fromInteger (s1+1)) (fromInteger (s2+1))
 208       where
 209         (q, s1) = s `divMod` 2147483562
 210         s2      = q `mod` 2147483398
 211
 212 -- FIXME: 1/2/3 below should be ** (vs@30082002) XXX
 213
 214 {- |
 215 With a source of random number supply in hand, the 'Random' class allows the
 216 programmer to extract random values of a variety of types.
 217
 218 Minimal complete definition: 'randomR' and 'random'.
 219
 220 -}
 221
 222 class Random a where
 223   -- | Takes a range /(lo,hi)/ and a random number generator
 224   -- /g/, and returns a random value uniformly distributed in the closed
 225   -- interval /[lo,hi]/, together with a new generator. It is unspecified
 226   -- what happens if /lo>hi/. For continuous types there is no requirement
 227   -- that the values /lo/ and /hi/ are ever produced, but they may be,
 228   -- depending on the implementation and the interval.
 229   randomR :: RandomGen g => (a,a) -> g -> (a,g)
 230
 231   -- | The same as 'randomR', but using a default range determined by the type:
 232   --
 233   -- * For bounded types (instances of 'Bounded', such as 'Char'),
 234   --   the range is normally the whole type.
 235   --
 236   -- * For fractional types, the range is normally the semi-closed interval
 237   -- @[0,1)@.
 238   --
 239   -- * For 'Integer', the range is (arbitrarily) the range of 'Int'.
 240   random  :: RandomGen g => g -> (a, g)
 241
 242   -- | Plural variant of 'randomR', producing an infinite list of
 243   -- random values instead of returning a new generator.
 244   randomRs :: RandomGen g => (a,a) -> g -> [a]
 245   randomRs ival g = x : randomRs ival g' where (x,g') = randomR ival g
 246
 247   -- | Plural variant of 'random', producing an infinite list of
 248   -- random values instead of returning a new generator.
 249   randoms  :: RandomGen g => g -> [a]
 250   randoms  g      = (\(x,g') -> x : randoms g') (random g)
 251
 252   -- | A variant of 'randomR' that uses the global random number generator
 253   -- (see "System.Random#globalrng").
 254   randomRIO :: (a,a) -> IO a
 255   randomRIO range  = getStdRandom (randomR range)
 256
 257   -- | A variant of 'random' that uses the global random number generator
 258   -- (see "System.Random#globalrng").
 259   randomIO  :: IO a
 260   randomIO         = getStdRandom random
 261
 262
 263 instance Random Int where
 264   randomR (a,b) g = randomIvalInteger (toInteger a, toInteger b) g
 265   random g        = randomR (minBound,maxBound) g
 266
 267 instance Random Char where
 268   randomR (a,b) g =
 269       case (randomIvalInteger (toInteger (ord a), toInteger (ord b)) g) of
 270         (x,g) -> (chr x, g)
 271   random g        = randomR (minBound,maxBound) g
 272
 273 instance Random Bool where
 274   randomR (a,b) g =
 275       case (randomIvalInteger (toInteger (bool2Int a), toInteger (bool2Int b)) g) of
 276         (x, g) -> (int2Bool x, g)
 277        where
 278          bool2Int False = 0
 279          bool2Int True  = 1
 280
 281          int2Bool 0     = False
 282          int2Bool _     = True
 283
 284   random g        = randomR (minBound,maxBound) g
 285
 286 instance Random Integer where
 287   randomR ival g = randomIvalInteger ival g
 288   random g       = randomR (toInteger (minBound::Int), toInteger (maxBound::Int)) g
 289
 290 instance Random Double where
 291   randomR ival g = randomIvalDouble ival id g
 292   random g       = randomR (0::Double,1) g
 293
 294 -- hah, so you thought you were saving cycles by using Float?
 295 instance Random Float where
 296   random g        = randomIvalDouble (0::Double,1) realToFrac g
 297   randomR (a,b) g = randomIvalDouble (realToFrac a, realToFrac b) realToFrac g
 298
 299 mkStdRNG :: Integer -> IO StdGen
 300 mkStdRNG o = do
 301     ct          <- getCPUTime
 302     (TOD sec _) <- getClockTime
 303     return (createStdGen (sec * 12345 + ct + o))
 304
 305 randomIvalInteger :: (RandomGen g, Num a) => (Integer, Integer) -> g -> (a, g)
 306 randomIvalInteger (l,h) rng
 307  | l > h     = randomIvalInteger (h,l) rng
 308  | otherwise = case (f n 1 rng) of (v, rng') -> (fromInteger (l + v `mod` k), rng')
 309      where
 310        k = h - l + 1
 311        b = 2147483561
 312        n = iLogBase b k
 313
 314        f 0 acc g = (acc, g)
 315        f n acc g =
 316           let
 317            (x,g')   = next g
 318           in
 319           f (n-1) (fromIntegral x + acc * b) g'
 320
 321 randomIvalDouble :: (RandomGen g, Fractional a) => (Double, Double) -> (Double -> a) -> g -> (a, g)
 322 randomIvalDouble (l,h) fromDouble rng
 323   | l > h     = randomIvalDouble (h,l) fromDouble rng
 324   | otherwise =
 325        case (randomIvalInteger (toInteger (minBound::Int), toInteger (maxBound::Int)) rng) of
 326          (x, rng') ->
 327             let
 328              scaled_x =
 329                 fromDouble ((l+h)/2) +
 330                 fromDouble ((h-l) / realToFrac intRange) *
 331                 fromIntegral (x::Int)
 332             in
 333             (scaled_x, rng')
 334
 335 intRange :: Integer
 336 intRange  = toInteger (maxBound::Int) - toInteger (minBound::Int)
 337
 338 iLogBase :: Integer -> Integer -> Integer
 339 iLogBase b i = if i < b then 1 else 1 + iLogBase b (i `div` b)
 340
 341 stdRange :: (Int,Int)
 342 stdRange = (0, 2147483562)
 343
 344 stdNext :: StdGen -> (Int, StdGen)
 345 -- Returns values in the range stdRange
 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 -}