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