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