844dea8d5cc4519505caf16dd8f3d0017a54f4f5
[packages/random.git] / System / Random.hs
1 #if __GLASGOW_HASKELL__ >= 701
2 {-# LANGUAGE Trustworthy #-}
3 #endif
4
5 -----------------------------------------------------------------------------
6 -- |
7 -- Module : System.Random
8 -- Copyright : (c) The University of Glasgow 2001
9 -- License : BSD-style (see the file LICENSE in the 'random' repository)
10 --
11 -- Maintainer : libraries@haskell.org
12 -- Stability : stable
13 -- Portability : portable
14 --
15 -- This library deals with the common task of pseudo-random number
16 -- generation. The library makes it possible to generate repeatable
17 -- results, by starting with a specified initial random number generator,
18 -- or to get different results on each run by using the system-initialised
19 -- generator or by supplying a seed from some other source.
20 --
21 -- The library is split into two layers:
22 --
23 -- * A core /random number generator/ provides a supply of bits.
24 -- The class 'RandomGen' provides a common interface to such generators.
25 -- The library provides one instance of 'RandomGen', the abstract
26 -- data type 'StdGen'. Programmers may, of course, supply their own
27 -- instances of 'RandomGen'.
28 --
29 -- * The class 'Random' provides a way to extract values of a particular
30 -- type from a random number generator. For example, the 'Float'
31 -- instance of 'Random' allows one to generate random values of type
32 -- 'Float'.
33 --
34 -- This implementation uses the Portable Combined Generator of L'Ecuyer
35 -- ["System.Random\#LEcuyer"] for 32-bit computers, transliterated by
36 -- Lennart Augustsson. It has a period of roughly 2.30584e18.
37 --
38 -----------------------------------------------------------------------------
39
40 #include "MachDeps.h"
41
42 module System.Random
43 (
44
45 -- $intro
46
47 -- * Random number generators
48
49 #ifdef ENABLE_SPLITTABLEGEN
50 RandomGen(next, genRange)
51 , SplittableGen(split)
52 #else
53 RandomGen(next, genRange, split)
54 #endif
55 -- ** Standard random number generators
56 , StdGen
57 , mkStdGen
58
59 -- ** The global random number generator
60
61 -- $globalrng
62
63 , getStdRandom
64 , getStdGen
65 , setStdGen
66 , newStdGen
67
68 -- * Random values of various types
69 , Random ( random, randomR,
70 randoms, randomRs,
71 randomIO, randomRIO )
72
73 -- * References
74 -- $references
75
76 ) where
77
78 import Prelude
79
80 import Data.Bits
81 import Data.Int
82 import Data.Word
83 import Foreign.C.Types
84
85 #ifdef __NHC__
86 import CPUTime ( getCPUTime )
87 import Foreign.Ptr ( Ptr, nullPtr )
88 import Foreign.C ( CTime, CUInt )
89 #else
90 import System.CPUTime ( getCPUTime )
91 import Data.Time ( getCurrentTime, UTCTime(..) )
92 import Data.Ratio ( numerator, denominator )
93 #endif
94 import Data.Char ( isSpace, chr, ord )
95 import System.IO.Unsafe ( unsafePerformIO )
96 import Data.IORef
97 import Numeric ( readDec )
98
99 #ifdef __GLASGOW_HASKELL__
100 import GHC.Exts ( build )
101 #else
102 -- | A dummy variant of build without fusion.
103 {-# INLINE build #-}
104 build :: ((a -> [a] -> [a]) -> [a] -> [a]) -> [a]
105 build g = g (:) []
106 #endif
107
108 -- The standard nhc98 implementation of Time.ClockTime does not match
109 -- the extended one expected in this module, so we lash-up a quick
110 -- replacement here.
111 #ifdef __NHC__
112 foreign import ccall "time.h time" readtime :: Ptr CTime -> IO CTime
113 getTime :: IO (Integer, Integer)
114 getTime = do CTime t <- readtime nullPtr; return (toInteger t, 0)
115 #else
116 getTime :: IO (Integer, Integer)
117 getTime = do
118 utc <- getCurrentTime
119 let daytime = toRational $ utctDayTime utc
120 return $ quotRem (numerator daytime) (denominator daytime)
121 #endif
122
123 -- | The class 'RandomGen' provides a common interface to random number
124 -- generators.
125 --
126 #ifdef ENABLE_SPLITTABLEGEN
127 -- Minimal complete definition: 'next'.
128 #else
129 -- Minimal complete definition: 'next' and 'split'.
130 #endif
131
132 class RandomGen g where
133
134 -- |The 'next' operation returns an 'Int' that is uniformly distributed
135 -- in the range returned by 'genRange' (including both end points),
136 -- and a new generator.
137 next :: g -> (Int, g)
138
139 -- |The 'genRange' operation yields the range of values returned by
140 -- the generator.
141 --
142 -- It is required that:
143 --
144 -- * If @(a,b) = 'genRange' g@, then @a < b@.
145 --
146 -- * 'genRange' always returns a pair of defined 'Int's.
147 --
148 -- The second condition ensures that 'genRange' cannot examine its
149 -- argument, and hence the value it returns can be determined only by the
150 -- instance of 'RandomGen'. That in turn allows an implementation to make
151 -- a single call to 'genRange' to establish a generator's range, without
152 -- being concerned that the generator returned by (say) 'next' might have
153 -- a different range to the generator passed to 'next'.
154 --
155 -- The default definition spans the full range of 'Int'.
156 genRange :: g -> (Int,Int)
157
158 -- default method
159 genRange _ = (minBound, maxBound)
160
161 #ifdef ENABLE_SPLITTABLEGEN
162 -- | The class 'SplittableGen' proivides a way to specify a random number
163 -- generator that can be split into two new generators.
164 class SplittableGen g where
165 #endif
166 -- |The 'split' operation allows one to obtain two distinct random number
167 -- generators. This is very useful in functional programs (for example, when
168 -- passing a random number generator down to recursive calls), but very
169 -- little work has been done on statistically robust implementations of
170 -- 'split' (["System.Random\#Burton", "System.Random\#Hellekalek"]
171 -- are the only examples we know of).
172 split :: g -> (g, g)
173
174 {- |
175 The 'StdGen' instance of 'RandomGen' has a 'genRange' of at least 30 bits.
176
177 The result of repeatedly using 'next' should be at least as statistically
178 robust as the /Minimal Standard Random Number Generator/ described by
179 ["System.Random\#Park", "System.Random\#Carta"].
180 Until more is known about implementations of 'split', all we require is
181 that 'split' deliver generators that are (a) not identical and
182 (b) independently robust in the sense just given.
183
184 The 'Show' and 'Read' instances of 'StdGen' provide a primitive way to save the
185 state of a random number generator.
186 It is required that @'read' ('show' g) == g@.
187
188 In addition, 'reads' may be used to map an arbitrary string (not necessarily one
189 produced by 'show') onto a value of type 'StdGen'. In general, the 'Read'
190 instance of 'StdGen' has the following properties:
191
192 * It guarantees to succeed on any string.
193
194 * It guarantees to consume only a finite portion of the string.
195
196 * Different argument strings are likely to result in different results.
197
198 -}
199
200 data StdGen
201 = StdGen Int32 Int32
202
203 instance RandomGen StdGen where
204 next = stdNext
205 genRange _ = stdRange
206
207 #ifdef ENABLE_SPLITTABLEGEN
208 instance SplittableGen StdGen where
209 #endif
210 split = stdSplit
211
212 instance Show StdGen where
213 showsPrec p (StdGen s1 s2) =
214 showsPrec p s1 .
215 showChar ' ' .
216 showsPrec p s2
217
218 instance Read StdGen where
219 readsPrec _p = \ r ->
220 case try_read r of
221 r'@[_] -> r'
222 _ -> [stdFromString r] -- because it shouldn't ever fail.
223 where
224 try_read r = do
225 (s1, r1) <- readDec (dropWhile isSpace r)
226 (s2, r2) <- readDec (dropWhile isSpace r1)
227 return (StdGen s1 s2, r2)
228
229 {-
230 If we cannot unravel the StdGen from a string, create
231 one based on the string given.
232 -}
233 stdFromString :: String -> (StdGen, String)
234 stdFromString s = (mkStdGen num, rest)
235 where (cs, rest) = splitAt 6 s
236 num = foldl (\a x -> x + 3 * a) 1 (map ord cs)
237
238
239 {- |
240 The function 'mkStdGen' provides an alternative way of producing an initial
241 generator, by mapping an 'Int' into a generator. Again, distinct arguments
242 should be likely to produce distinct generators.
243 -}
244 mkStdGen :: Int -> StdGen -- why not Integer ?
245 mkStdGen s = mkStdGen32 $ fromIntegral s
246
247 mkStdGen32 :: Int32 -> StdGen
248 mkStdGen32 sMaybeNegative = StdGen (s1+1) (s2+1)
249 where
250 -- We want a non-negative number, but we can't just take the abs
251 -- of sMaybeNegative as -minBound == minBound.
252 s = sMaybeNegative .&. maxBound
253 (q, s1) = s `divMod` 2147483562
254 s2 = q `mod` 2147483398
255
256 createStdGen :: Integer -> StdGen
257 createStdGen s = mkStdGen32 $ fromIntegral s
258
259 -- FIXME: 1/2/3 below should be ** (vs@30082002) XXX
260
261 {- |
262 With a source of random number supply in hand, the 'Random' class allows the
263 programmer to extract random values of a variety of types.
264
265 Minimal complete definition: 'randomR' and 'random'.
266
267 -}
268
269 class Random a where
270 -- | Takes a range /(lo,hi)/ and a random number generator
271 -- /g/, and returns a random value uniformly distributed in the closed
272 -- interval /[lo,hi]/, together with a new generator. It is unspecified
273 -- what happens if /lo>hi/. For continuous types there is no requirement
274 -- that the values /lo/ and /hi/ are ever produced, but they may be,
275 -- depending on the implementation and the interval.
276 randomR :: RandomGen g => (a,a) -> g -> (a,g)
277
278 -- | The same as 'randomR', but using a default range determined by the type:
279 --
280 -- * For bounded types (instances of 'Bounded', such as 'Char'),
281 -- the range is normally the whole type.
282 --
283 -- * For fractional types, the range is normally the semi-closed interval
284 -- @[0,1)@.
285 --
286 -- * For 'Integer', the range is (arbitrarily) the range of 'Int'.
287 random :: RandomGen g => g -> (a, g)
288
289 -- | Plural variant of 'randomR', producing an infinite list of
290 -- random values instead of returning a new generator.
291 {-# INLINE randomRs #-}
292 randomRs :: RandomGen g => (a,a) -> g -> [a]
293 randomRs ival g = build (\cons _nil -> buildRandoms cons (randomR ival) g)
294
295 -- | Plural variant of 'random', producing an infinite list of
296 -- random values instead of returning a new generator.
297 {-# INLINE randoms #-}
298 randoms :: RandomGen g => g -> [a]
299 randoms g = build (\cons _nil -> buildRandoms cons random g)
300
301 -- | A variant of 'randomR' that uses the global random number generator
302 -- (see "System.Random#globalrng").
303 randomRIO :: (a,a) -> IO a
304 randomRIO range = getStdRandom (randomR range)
305
306 -- | A variant of 'random' that uses the global random number generator
307 -- (see "System.Random#globalrng").
308 randomIO :: IO a
309 randomIO = getStdRandom random
310
311 -- | Produce an infinite list-equivalent of random values.
312 {-# INLINE buildRandoms #-}
313 buildRandoms :: RandomGen g
314 => (a -> as -> as) -- ^ E.g. '(:)' but subject to fusion
315 -> (g -> (a,g)) -- ^ E.g. 'random'
316 -> g -- ^ A 'RandomGen' instance
317 -> as
318 buildRandoms cons rand = go
319 where
320 -- The seq fixes part of #4218 and also makes fused Core simpler.
321 go g = x `seq` (x `cons` go g') where (x,g') = rand g
322
323
324 instance Random Integer where
325 randomR ival g = randomIvalInteger ival g
326 random g = randomR (toInteger (minBound::Int), toInteger (maxBound::Int)) g
327
328 instance Random Int where randomR = randomIvalIntegral; random = randomBounded
329 instance Random Int8 where randomR = randomIvalIntegral; random = randomBounded
330 instance Random Int16 where randomR = randomIvalIntegral; random = randomBounded
331 instance Random Int32 where randomR = randomIvalIntegral; random = randomBounded
332 instance Random Int64 where randomR = randomIvalIntegral; random = randomBounded
333
334 #ifndef __NHC__
335 -- Word is a type synonym in nhc98.
336 instance Random Word where randomR = randomIvalIntegral; random = randomBounded
337 #endif
338 instance Random Word8 where randomR = randomIvalIntegral; random = randomBounded
339 instance Random Word16 where randomR = randomIvalIntegral; random = randomBounded
340 instance Random Word32 where randomR = randomIvalIntegral; random = randomBounded
341 instance Random Word64 where randomR = randomIvalIntegral; random = randomBounded
342
343 instance Random CChar where randomR = randomIvalIntegral; random = randomBounded
344 instance Random CSChar where randomR = randomIvalIntegral; random = randomBounded
345 instance Random CUChar where randomR = randomIvalIntegral; random = randomBounded
346 instance Random CShort where randomR = randomIvalIntegral; random = randomBounded
347 instance Random CUShort where randomR = randomIvalIntegral; random = randomBounded
348 instance Random CInt where randomR = randomIvalIntegral; random = randomBounded
349 instance Random CUInt where randomR = randomIvalIntegral; random = randomBounded
350 instance Random CLong where randomR = randomIvalIntegral; random = randomBounded
351 instance Random CULong where randomR = randomIvalIntegral; random = randomBounded
352 instance Random CPtrdiff where randomR = randomIvalIntegral; random = randomBounded
353 instance Random CSize where randomR = randomIvalIntegral; random = randomBounded
354 instance Random CWchar where randomR = randomIvalIntegral; random = randomBounded
355 instance Random CSigAtomic where randomR = randomIvalIntegral; random = randomBounded
356 instance Random CLLong where randomR = randomIvalIntegral; random = randomBounded
357 instance Random CULLong where randomR = randomIvalIntegral; random = randomBounded
358 instance Random CIntPtr where randomR = randomIvalIntegral; random = randomBounded
359 instance Random CUIntPtr where randomR = randomIvalIntegral; random = randomBounded
360 instance Random CIntMax where randomR = randomIvalIntegral; random = randomBounded
361 instance Random CUIntMax where randomR = randomIvalIntegral; random = randomBounded
362
363 instance Random Char where
364 randomR (a,b) g =
365 case (randomIvalInteger (toInteger (ord a), toInteger (ord b)) g) of
366 (x,g') -> (chr x, g')
367 random g = randomR (minBound,maxBound) g
368
369 instance Random Bool where
370 randomR (a,b) g =
371 case (randomIvalInteger (bool2Int a, bool2Int b) g) of
372 (x, g') -> (int2Bool x, g')
373 where
374 bool2Int :: Bool -> Integer
375 bool2Int False = 0
376 bool2Int True = 1
377
378 int2Bool :: Int -> Bool
379 int2Bool 0 = False
380 int2Bool _ = True
381
382 random g = randomR (minBound,maxBound) g
383
384 {-# INLINE randomRFloating #-}
385 randomRFloating :: (Fractional a, Num a, Ord a, Random a, RandomGen g) => (a, a) -> g -> (a, g)
386 randomRFloating (l,h) g
387 | l>h = randomRFloating (h,l) g
388 | otherwise = let (coef,g') = random g in
389 (2.0 * (0.5*l + coef * (0.5*h - 0.5*l)), g') -- avoid overflow
390
391 instance Random Double where
392 randomR = randomRFloating
393 random rng =
394 case random rng of
395 (x,rng') ->
396 -- We use 53 bits of randomness corresponding to the 53 bit significand:
397 ((fromIntegral (mask53 .&. (x::Int64)) :: Double)
398 / fromIntegral twoto53, rng')
399 where
400 twoto53 = (2::Int64) ^ (53::Int64)
401 mask53 = twoto53 - 1
402
403 instance Random Float where
404 randomR = randomRFloating
405 random rng =
406 -- TODO: Faster to just use 'next' IF it generates enough bits of randomness.
407 case random rng of
408 (x,rng') ->
409 -- We use 24 bits of randomness corresponding to the 24 bit significand:
410 ((fromIntegral (mask24 .&. (x::Int32)) :: Float)
411 / fromIntegral twoto24, rng')
412 -- Note, encodeFloat is another option, but I'm not seeing slightly
413 -- worse performance with the following [2011.06.25]:
414 -- (encodeFloat rand (-24), rng')
415 where
416 mask24 = twoto24 - 1
417 twoto24 = (2::Int32) ^ (24::Int32)
418
419 -- CFloat/CDouble are basically the same as a Float/Double:
420 instance Random CFloat where
421 randomR = randomRFloating
422 random rng = case random rng of
423 (x,rng') -> (realToFrac (x::Float), rng')
424
425 instance Random CDouble where
426 randomR = randomRFloating
427 -- A MYSTERY:
428 -- Presently, this is showing better performance than the Double instance:
429 -- (And yet, if the Double instance uses randomFrac then its performance is much worse!)
430 random = randomFrac
431 -- random rng = case random rng of
432 -- (x,rng') -> (realToFrac (x::Double), rng')
433
434 mkStdRNG :: Integer -> IO StdGen
435 mkStdRNG o = do
436 ct <- getCPUTime
437 (sec, psec) <- getTime
438 return (createStdGen (sec * 12345 + psec + ct + o))
439
440 randomBounded :: (RandomGen g, Random a, Bounded a) => g -> (a, g)
441 randomBounded = randomR (minBound, maxBound)
442
443 -- The two integer functions below take an [inclusive,inclusive] range.
444 randomIvalIntegral :: (RandomGen g, Integral a) => (a, a) -> g -> (a, g)
445 randomIvalIntegral (l,h) = randomIvalInteger (toInteger l, toInteger h)
446
447 randomIvalInteger :: (RandomGen g, Num a) => (Integer, Integer) -> g -> (a, g)
448 randomIvalInteger (l,h) rng
449 | l > h = randomIvalInteger (h,l) rng
450 | otherwise = case (f n 1 rng) of (v, rng') -> (fromInteger (l + v `mod` k), rng')
451 where
452 k = h - l + 1
453 -- ERROR: b here (2^31-87) represents a baked-in assumption about genRange:
454 b = 2147483561
455 n = iLogBase b k
456
457 -- Here we loop until we've generated enough randomness to cover the range:
458 f 0 acc g = (acc, g)
459 f n' acc g =
460 let
461 (x,g') = next g
462 in
463 -- We shift over the random bits generated thusfar (* b) and add in the new ones.
464 f (n' - 1) (fromIntegral x + acc * b) g'
465
466 -- The continuous functions on the other hand take an [inclusive,exclusive) range.
467 randomFrac :: (RandomGen g, Fractional a) => g -> (a, g)
468 randomFrac = randomIvalDouble (0::Double,1) realToFrac
469
470 randomIvalDouble :: (RandomGen g, Fractional a) => (Double, Double) -> (Double -> a) -> g -> (a, g)
471 randomIvalDouble (l,h) fromDouble rng
472 | l > h = randomIvalDouble (h,l) fromDouble rng
473 | otherwise =
474 case (randomIvalInteger (toInteger (minBound::Int32), toInteger (maxBound::Int32)) rng) of
475 (x, rng') ->
476 let
477 scaled_x =
478 fromDouble (0.5*l + 0.5*h) + -- previously (l+h)/2, overflowed
479 fromDouble ((0.5*h - 0.5*l) / (0.5 * realToFrac int32Count)) * -- avoid overflow
480 fromIntegral (x::Int32)
481 in
482 (scaled_x, rng')
483
484 int32Count :: Integer
485 int32Count = toInteger (maxBound::Int32) - toInteger (minBound::Int32) + 1
486
487 -- Perform an expensive logarithm on arbitrary-size integers by repeated division.
488 --
489 -- (NOTE: This actually returns ceiling(log(i) base b) except with an
490 -- incorrect result at iLogBase b b = 2.)
491 iLogBase :: Integer -> Integer -> Integer
492 iLogBase b i = if i < b then 1 else 1 + iLogBase b (i `div` b)
493
494 stdRange :: (Int,Int)
495 stdRange = (0, 2147483562)
496
497 stdNext :: StdGen -> (Int, StdGen)
498 -- Returns values in the range stdRange
499 stdNext (StdGen s1 s2) = (fromIntegral z', StdGen s1'' s2'')
500 where z' = if z < 1 then z + 2147483562 else z
501 z = s1'' - s2''
502
503 k = s1 `quot` 53668
504 s1' = 40014 * (s1 - k * 53668) - k * 12211
505 s1'' = if s1' < 0 then s1' + 2147483563 else s1'
506
507 k' = s2 `quot` 52774
508 s2' = 40692 * (s2 - k' * 52774) - k' * 3791
509 s2'' = if s2' < 0 then s2' + 2147483399 else s2'
510
511 stdSplit :: StdGen -> (StdGen, StdGen)
512 stdSplit std@(StdGen s1 s2)
513 = (left, right)
514 where
515 -- no statistical foundation for this!
516 left = StdGen new_s1 t2
517 right = StdGen t1 new_s2
518
519 new_s1 | s1 == 2147483562 = 1
520 | otherwise = s1 + 1
521
522 new_s2 | s2 == 1 = 2147483398
523 | otherwise = s2 - 1
524
525 StdGen t1 t2 = snd (next std)
526
527 -- The global random number generator
528
529 {- $globalrng #globalrng#
530
531 There is a single, implicit, global random number generator of type
532 'StdGen', held in some global variable maintained by the 'IO' monad. It is
533 initialised automatically in some system-dependent fashion, for example, by
534 using the time of day, or Linux's kernel random number generator. To get
535 deterministic behaviour, use 'setStdGen'.
536 -}
537
538 -- |Sets the global random number generator.
539 setStdGen :: StdGen -> IO ()
540 setStdGen sgen = writeIORef theStdGen sgen
541
542 -- |Gets the global random number generator.
543 getStdGen :: IO StdGen
544 getStdGen = readIORef theStdGen
545
546 theStdGen :: IORef StdGen
547 theStdGen = unsafePerformIO $ do
548 rng <- mkStdRNG 0
549 newIORef rng
550
551 -- |Applies 'split' to the current global random generator,
552 -- updates it with one of the results, and returns the other.
553 newStdGen :: IO StdGen
554 newStdGen = atomicModifyIORef theStdGen split
555
556 {- |Uses the supplied function to get a value from the current global
557 random generator, and updates the global generator with the new generator
558 returned by the function. For example, @rollDice@ gets a random integer
559 between 1 and 6:
560
561 > rollDice :: IO Int
562 > rollDice = getStdRandom (randomR (1,6))
563
564 -}
565
566 getStdRandom :: (StdGen -> (a,StdGen)) -> IO a
567 getStdRandom f = atomicModifyIORef theStdGen (swap . f)
568 where swap (v,g) = (g,v)
569
570 {- $references
571
572 1. FW #Burton# Burton and RL Page, /Distributed random number generation/,
573 Journal of Functional Programming, 2(2):203-212, April 1992.
574
575 2. SK #Park# Park, and KW Miller, /Random number generators -
576 good ones are hard to find/, Comm ACM 31(10), Oct 1988, pp1192-1201.
577
578 3. DG #Carta# Carta, /Two fast implementations of the minimal standard
579 random number generator/, Comm ACM, 33(1), Jan 1990, pp87-88.
580
581 4. P #Hellekalek# Hellekalek, /Don\'t trust parallel Monte Carlo/,
582 Department of Mathematics, University of Salzburg,
583 <http://random.mat.sbg.ac.at/~peter/pads98.ps>, 1998.
584
585 5. Pierre #LEcuyer# L'Ecuyer, /Efficient and portable combined random
586 number generators/, Comm ACM, 31(6), Jun 1988, pp742-749.
587
588 The Web site <http://random.mat.sbg.ac.at/> is a great source of information.
589
590 -}