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