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