f5b1360d0589167b3c94ecc42bac6619782c8dae
[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.Int
72 import Data.Word
73 import Foreign.C.Types
74
75 #ifdef __NHC__
76 import CPUTime ( getCPUTime )
77 import Foreign.Ptr ( Ptr, nullPtr )
78 import Foreign.C ( CTime, CUInt )
79 #else
80 import System.CPUTime ( getCPUTime )
81 import Data.Time ( getCurrentTime, UTCTime(..) )
82 import Data.Ratio ( numerator, denominator )
83 #endif
84 import Data.Char ( isSpace, chr, ord )
85 import System.IO.Unsafe ( unsafePerformIO )
86 import Data.IORef
87 import Numeric ( readDec )
88
89 -- The standard nhc98 implementation of Time.ClockTime does not match
90 -- the extended one expected in this module, so we lash-up a quick
91 -- replacement here.
92 #ifdef __NHC__
93 foreign import ccall "time.h time" readtime :: Ptr CTime -> IO CTime
94 getTime :: IO (Integer, Integer)
95 getTime = do CTime t <- readtime nullPtr; return (toInteger t, 0)
96 #else
97 getTime :: IO (Integer, Integer)
98 getTime = do
99 utc <- getCurrentTime
100 let daytime = toRational $ utctDayTime utc
101 return $ quotRem (numerator daytime) (denominator daytime)
102 #endif
103
104 -- | The class 'RandomGen' provides a common interface to random number
105 -- generators.
106 --
107 -- Minimal complete definition: 'next'.
108
109 class RandomGen g where
110
111 -- |The 'next' operation returns an 'Int' that is uniformly distributed
112 -- in the range returned by 'genRange' (including both end points),
113 -- and a new generator.
114 next :: g -> (Int, g)
115
116 -- |The 'genRange' operation yields the range of values returned by
117 -- the generator.
118 --
119 -- It is required that:
120 --
121 -- * If @(a,b) = 'genRange' g@, then @a < b@.
122 --
123 -- * 'genRange' always returns a pair of defined 'Int's.
124 --
125 -- The second condition ensures that 'genRange' cannot examine its
126 -- argument, and hence the value it returns can be determined only by the
127 -- instance of 'RandomGen'. That in turn allows an implementation to make
128 -- a single call to 'genRange' to establish a generator's range, without
129 -- being concerned that the generator returned by (say) 'next' might have
130 -- a different range to the generator passed to 'next'.
131 --
132 -- The default definition spans the full range of 'Int'.
133 genRange :: g -> (Int,Int)
134
135 -- default method
136 genRange _ = (minBound, maxBound)
137
138 -- | The class 'SplittableGen' proivides a way to specify a random number
139 -- generator that can be split into two new generators.
140 class SplittableGen g where
141 -- |The 'split' operation allows one to obtain two distinct random number
142 -- generators. This is very useful in functional programs (for example, when
143 -- passing a random number generator down to recursive calls), but very
144 -- little work has been done on statistically robust implementations of
145 -- 'split' (["System.Random\#Burton", "System.Random\#Hellekalek"]
146 -- are the only examples we know of).
147 split :: g -> (g, g)
148
149 {- |
150 The 'StdGen' instance of 'RandomGen' has a 'genRange' of at least 30 bits.
151
152 The result of repeatedly using 'next' should be at least as statistically
153 robust as the /Minimal Standard Random Number Generator/ described by
154 ["System.Random\#Park", "System.Random\#Carta"].
155 Until more is known about implementations of 'split', all we require is
156 that 'split' deliver generators that are (a) not identical and
157 (b) independently robust in the sense just given.
158
159 The 'Show' and 'Read' instances of 'StdGen' provide a primitive way to save the
160 state of a random number generator.
161 It is required that @'read' ('show' g) == g@.
162
163 In addition, 'reads' may be used to map an arbitrary string (not necessarily one
164 produced by 'show') onto a value of type 'StdGen'. In general, the 'Read'
165 instance of 'StdGen' has the following properties:
166
167 * It guarantees to succeed on any string.
168
169 * It guarantees to consume only a finite portion of the string.
170
171 * Different argument strings are likely to result in different results.
172
173 -}
174
175 data StdGen
176 = StdGen Int32 Int32
177
178 instance RandomGen StdGen where
179 next = stdNext
180 genRange _ = stdRange
181
182 instance SplittableGen StdGen where
183 split = stdSplit
184
185 instance Show StdGen where
186 showsPrec p (StdGen s1 s2) =
187 showsPrec p s1 .
188 showChar ' ' .
189 showsPrec p s2
190
191 instance Read StdGen where
192 readsPrec _p = \ r ->
193 case try_read r of
194 r'@[_] -> r'
195 _ -> [stdFromString r] -- because it shouldn't ever fail.
196 where
197 try_read r = do
198 (s1, r1) <- readDec (dropWhile isSpace r)
199 (s2, r2) <- readDec (dropWhile isSpace r1)
200 return (StdGen s1 s2, r2)
201
202 {-
203 If we cannot unravel the StdGen from a string, create
204 one based on the string given.
205 -}
206 stdFromString :: String -> (StdGen, String)
207 stdFromString s = (mkStdGen num, rest)
208 where (cs, rest) = splitAt 6 s
209 num = foldl (\a x -> x + 3 * a) 1 (map ord cs)
210
211
212 {- |
213 The function 'mkStdGen' provides an alternative way of producing an initial
214 generator, by mapping an 'Int' into a generator. Again, distinct arguments
215 should be likely to produce distinct generators.
216 -}
217 mkStdGen :: Int -> StdGen -- why not Integer ?
218 mkStdGen s = mkStdGen32 $ fromIntegral s
219
220 mkStdGen32 :: Int32 -> StdGen
221 mkStdGen32 s
222 | s < 0 = mkStdGen32 (-s)
223 | otherwise = StdGen (s1+1) (s2+1)
224 where
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 -}