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 CTime -> IO CTime
  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 -- Minimal complete definition: 'next' and 'split'.
  96
  97 class RandomGen g where
  98
  99    -- |The 'next' operation returns an 'Int' that is uniformly distributed
 100    -- in the range returned by 'genRange' (including both end points),
 101    -- and a new generator.
 102    next     :: g -> (Int, g)
 103
 104    -- |The 'split' operation allows one to obtain two distinct random number
 105    -- generators. This is very useful in functional programs (for example, when
 106    -- passing a random number generator down to recursive calls), but very
 107    -- little work has been done on statistically robust implementations of
 108    -- 'split' (["System.Random\#Burton", "System.Random\#Hellekalek"]
 109    -- are the only examples we know of).
 110    split    :: g -> (g, g)
 111
 112    -- |The 'genRange' operation yields the range of values returned by
 113    -- the generator.
 114    --
 115    -- It is required that:
 116    --
 117    -- * If @(a,b) = 'genRange' g@, then @a < b@.
 118    --
 119    -- * 'genRange' always returns a pair of defined 'Int's.
 120    --
 121    -- The second condition ensures that 'genRange' cannot examine its
 122    -- argument, and hence the value it returns can be determined only by the
 123    -- instance of 'RandomGen'.  That in turn allows an implementation to make
 124    -- a single call to 'genRange' to establish a generator's range, without
 125    -- being concerned that the generator returned by (say) 'next' might have
 126    -- a different range to the generator passed to 'next'.
 127    --
 128    -- The default definition spans the full range of 'Int'.
 129    genRange :: g -> (Int,Int)
 130
 131    -- default method
 132    genRange g = (minBound,maxBound)
 133
 134 {- |
 135 The 'StdGen' instance of 'RandomGen' has a 'genRange' of at least 30 bits.
 136
 137 The result of repeatedly using 'next' should be at least as statistically
 138 robust as the /Minimal Standard Random Number Generator/ described by
 139 ["System.Random\#Park", "System.Random\#Carta"].
 140 Until more is known about implementations of 'split', all we require is
 141 that 'split' deliver generators that are (a) not identical and
 142 (b) independently robust in the sense just given.
 143
 144 The 'Show' and 'Read' instances of 'StdGen' provide a primitive way to save the
 145 state of a random number generator.
 146 It is required that @'read' ('show' g) == g@.
 147
 148 In addition, 'read' may be used to map an arbitrary string (not necessarily one
 149 produced by 'show') onto a value of type 'StdGen'. In general, the 'read'
 150 instance of 'StdGen' has the following properties:
 151
 152 * It guarantees to succeed on any string.
 153
 154 * It guarantees to consume only a finite portion of the string.
 155
 156 * Different argument strings are likely to result in different results.
 157
 158 -}
 159
 160 data StdGen
 161  = StdGen Int Int
 162
 163 instance RandomGen StdGen where
 164   next  = stdNext
 165   split = stdSplit
 166   genRange _ = stdRange
 167
 168 instance Show StdGen where
 169   showsPrec p (StdGen s1 s2) =
 170      showsPrec p s1 .
 171      showChar ' ' .
 172      showsPrec p s2
 173
 174 instance Read StdGen where
 175   readsPrec _p = \ r ->
 176      case try_read r of
 177        r@[_] -> r
 178        _   -> [stdFromString r] -- because it shouldn't ever fail.
 179     where
 180       try_read r = do
 181          (s1, r1) <- readDec (dropWhile isSpace r)
 182          (s2, r2) <- readDec (dropWhile isSpace r1)
 183          return (StdGen s1 s2, r2)
 184
 185 {-
 186  If we cannot unravel the StdGen from a string, create
 187  one based on the string given.
 188 -}
 189 stdFromString         :: String -> (StdGen, String)
 190 stdFromString s        = (mkStdGen num, rest)
 191         where (cs, rest) = splitAt 6 s
 192               num        = foldl (\a x -> x + 3 * a) 1 (map ord cs)
 193
 194
 195 {- |
 196 The function 'mkStdGen' provides an alternative way of producing an initial
 197 generator, by mapping an 'Int' into a generator. Again, distinct arguments
 198 should be likely to produce distinct generators.
 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 stdRange :: (Int,Int)
 346 stdRange = (0, 2147483562)
 347
 348 stdNext :: StdGen -> (Int, StdGen)
 349 -- Returns values in the range stdRange
 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# 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# 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 5. Pierre #LEcuyer# L'Ecuyer, /Efficient and portable combined random
 444 number generators/, Comm ACM, 31(6), Jun 1988, pp742-749.
 445
 446 The Web site <http://random.mat.sbg.ac.at/> is a great source of information.
 447
 448 -}