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