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