Add a comment about why we can't use abs to get a non-negative number
[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 -- We want a non-negative number, but we can't just take the abs
225 -- of sMaybeNegative as -minBound == minBound.
226 s = sMaybeNegative .&. maxBound
227 (q, s1) = s `divMod` 2147483562
228 s2 = q `mod` 2147483398
229
230 createStdGen :: Integer -> StdGen
231 createStdGen s = mkStdGen32 $ fromIntegral s
232
233 -- FIXME: 1/2/3 below should be ** (vs@30082002) XXX
234
235 {- |
236 With a source of random number supply in hand, the 'Random' class allows the
237 programmer to extract random values of a variety of types.
238
239 Minimal complete definition: 'randomR' and 'random'.
240
241 -}
242
243 class Random a where
244 -- | Takes a range /(lo,hi)/ and a random number generator
245 -- /g/, and returns a random value uniformly distributed in the closed
246 -- interval /[lo,hi]/, together with a new generator. It is unspecified
247 -- what happens if /lo>hi/. For continuous types there is no requirement
248 -- that the values /lo/ and /hi/ are ever produced, but they may be,
249 -- depending on the implementation and the interval.
250 randomR :: RandomGen g => (a,a) -> g -> (a,g)
251
252 -- | The same as 'randomR', but using a default range determined by the type:
253 --
254 -- * For bounded types (instances of 'Bounded', such as 'Char'),
255 -- the range is normally the whole type.
256 --
257 -- * For fractional types, the range is normally the semi-closed interval
258 -- @[0,1)@.
259 --
260 -- * For 'Integer', the range is (arbitrarily) the range of 'Int'.
261 random :: RandomGen g => g -> (a, g)
262
263 -- | Plural variant of 'randomR', producing an infinite list of
264 -- random values instead of returning a new generator.
265 randomRs :: RandomGen g => (a,a) -> g -> [a]
266 randomRs ival g = x : randomRs ival g' where (x,g') = randomR ival g
267
268 -- | Plural variant of 'random', producing an infinite list of
269 -- random values instead of returning a new generator.
270 randoms :: RandomGen g => g -> [a]
271 randoms g = (\(x,g') -> x : randoms g') (random g)
272
273 -- | A variant of 'randomR' that uses the global random number generator
274 -- (see "System.Random#globalrng").
275 randomRIO :: (a,a) -> IO a
276 randomRIO range = getStdRandom (randomR range)
277
278 -- | A variant of 'random' that uses the global random number generator
279 -- (see "System.Random#globalrng").
280 randomIO :: IO a
281 randomIO = getStdRandom random
282
283
284 instance Random Integer where
285 randomR ival g = randomIvalInteger ival g
286 random g = randomR (toInteger (minBound::Int), toInteger (maxBound::Int)) g
287
288 instance Random Int where randomR = randomIvalIntegral; random = randomBounded
289 instance Random Int8 where randomR = randomIvalIntegral; random = randomBounded
290 instance Random Int16 where randomR = randomIvalIntegral; random = randomBounded
291 instance Random Int32 where randomR = randomIvalIntegral; random = randomBounded
292 instance Random Int64 where randomR = randomIvalIntegral; random = randomBounded
293
294 #ifndef __NHC__
295 -- Word is a type synonym in nhc98.
296 instance Random Word where randomR = randomIvalIntegral; random = randomBounded
297 #endif
298 instance Random Word8 where randomR = randomIvalIntegral; random = randomBounded
299 instance Random Word16 where randomR = randomIvalIntegral; random = randomBounded
300 instance Random Word32 where randomR = randomIvalIntegral; random = randomBounded
301 instance Random Word64 where randomR = randomIvalIntegral; random = randomBounded
302
303 instance Random CChar where randomR = randomIvalIntegral; random = randomBounded
304 instance Random CSChar where randomR = randomIvalIntegral; random = randomBounded
305 instance Random CUChar where randomR = randomIvalIntegral; random = randomBounded
306 instance Random CShort where randomR = randomIvalIntegral; random = randomBounded
307 instance Random CUShort where randomR = randomIvalIntegral; random = randomBounded
308 instance Random CInt where randomR = randomIvalIntegral; random = randomBounded
309 instance Random CUInt where randomR = randomIvalIntegral; random = randomBounded
310 instance Random CLong where randomR = randomIvalIntegral; random = randomBounded
311 instance Random CULong where randomR = randomIvalIntegral; random = randomBounded
312 instance Random CPtrdiff where randomR = randomIvalIntegral; random = randomBounded
313 instance Random CSize where randomR = randomIvalIntegral; random = randomBounded
314 instance Random CWchar where randomR = randomIvalIntegral; random = randomBounded
315 instance Random CSigAtomic where randomR = randomIvalIntegral; random = randomBounded
316 instance Random CLLong where randomR = randomIvalIntegral; random = randomBounded
317 instance Random CULLong where randomR = randomIvalIntegral; random = randomBounded
318 instance Random CIntPtr where randomR = randomIvalIntegral; random = randomBounded
319 instance Random CUIntPtr where randomR = randomIvalIntegral; random = randomBounded
320 instance Random CIntMax where randomR = randomIvalIntegral; random = randomBounded
321 instance Random CUIntMax where randomR = randomIvalIntegral; random = randomBounded
322
323 instance Random Char where
324 randomR (a,b) g =
325 case (randomIvalInteger (toInteger (ord a), toInteger (ord b)) g) of
326 (x,g') -> (chr x, g')
327 random g = randomR (minBound,maxBound) g
328
329 instance Random Bool where
330 randomR (a,b) g =
331 case (randomIvalInteger (bool2Int a, bool2Int b) g) of
332 (x, g') -> (int2Bool x, g')
333 where
334 bool2Int :: Bool -> Integer
335 bool2Int False = 0
336 bool2Int True = 1
337
338 int2Bool :: Int -> Bool
339 int2Bool 0 = False
340 int2Bool _ = True
341
342 random g = randomR (minBound,maxBound) g
343
344 instance Random Double where
345 randomR ival g = randomIvalDouble ival id g
346 random g = randomR (0::Double,1) g
347
348 -- hah, so you thought you were saving cycles by using Float?
349 instance Random Float where
350 randomR = randomIvalFrac
351 random = randomFrac
352
353 instance Random CFloat where
354 randomR = randomIvalFrac
355 random = randomFrac
356
357 instance Random CDouble where
358 randomR = randomIvalFrac
359 random = randomFrac
360
361
362 mkStdRNG :: Integer -> IO StdGen
363 mkStdRNG o = do
364 ct <- getCPUTime
365 (sec, psec) <- getTime
366 return (createStdGen (sec * 12345 + psec + ct + o))
367
368 randomBounded :: (RandomGen g, Random a, Bounded a) => g -> (a, g)
369 randomBounded = randomR (minBound, maxBound)
370
371 randomIvalIntegral :: (RandomGen g, Integral a) => (a, a) -> g -> (a, g)
372 randomIvalIntegral (l,h) = randomIvalInteger (toInteger l, toInteger h)
373
374 randomIvalInteger :: (RandomGen g, Num a) => (Integer, Integer) -> g -> (a, g)
375 randomIvalInteger (l,h) rng
376 | l > h = randomIvalInteger (h,l) rng
377 | otherwise = case (f n 1 rng) of (v, rng') -> (fromInteger (l + v `mod` k), rng')
378 where
379 k = h - l + 1
380 b = 2147483561
381 n = iLogBase b k
382
383 f 0 acc g = (acc, g)
384 f n' acc g =
385 let
386 (x,g') = next g
387 in
388 f (n' - 1) (fromIntegral x + acc * b) g'
389
390 randomFrac :: (RandomGen g, Fractional a) => g -> (a, g)
391 randomFrac = randomIvalDouble (0::Double,1) realToFrac
392
393 randomIvalFrac :: (RandomGen g, Real a, Fractional b) => (a,a) -> g -> (b, g)
394 randomIvalFrac (a,b) = randomIvalDouble (realToFrac a, realToFrac b) realToFrac
395
396 randomIvalDouble :: (RandomGen g, Fractional a) => (Double, Double) -> (Double -> a) -> g -> (a, g)
397 randomIvalDouble (l,h) fromDouble rng
398 | l > h = randomIvalDouble (h,l) fromDouble rng
399 | otherwise =
400 case (randomIvalInteger (toInteger (minBound::Int32), toInteger (maxBound::Int32)) rng) of
401 (x, rng') ->
402 let
403 scaled_x =
404 fromDouble ((l+h)/2) +
405 fromDouble ((h-l) / realToFrac int32Count) *
406 fromIntegral (x::Int32)
407 in
408 (scaled_x, rng')
409
410 int32Count :: Integer
411 int32Count = toInteger (maxBound::Int32) - toInteger (minBound::Int32) + 1
412
413 iLogBase :: Integer -> Integer -> Integer
414 iLogBase b i = if i < b then 1 else 1 + iLogBase b (i `div` b)
415
416 stdRange :: (Int,Int)
417 stdRange = (0, 2147483562)
418
419 stdNext :: StdGen -> (Int, StdGen)
420 -- Returns values in the range stdRange
421 stdNext (StdGen s1 s2) = (fromIntegral z', StdGen s1'' s2'')
422 where z' = if z < 1 then z + 2147483562 else z
423 z = s1'' - s2''
424
425 k = s1 `quot` 53668
426 s1' = 40014 * (s1 - k * 53668) - k * 12211
427 s1'' = if s1' < 0 then s1' + 2147483563 else s1'
428
429 k' = s2 `quot` 52774
430 s2' = 40692 * (s2 - k' * 52774) - k' * 3791
431 s2'' = if s2' < 0 then s2' + 2147483399 else s2'
432
433 stdSplit :: StdGen -> (StdGen, StdGen)
434 stdSplit std@(StdGen s1 s2)
435 = (left, right)
436 where
437 -- no statistical foundation for this!
438 left = StdGen new_s1 t2
439 right = StdGen t1 new_s2
440
441 new_s1 | s1 == 2147483562 = 1
442 | otherwise = s1 + 1
443
444 new_s2 | s2 == 1 = 2147483398
445 | otherwise = s2 - 1
446
447 StdGen t1 t2 = snd (next std)
448
449 -- The global random number generator
450
451 {- $globalrng #globalrng#
452
453 There is a single, implicit, global random number generator of type
454 'StdGen', held in some global variable maintained by the 'IO' monad. It is
455 initialised automatically in some system-dependent fashion, for example, by
456 using the time of day, or Linux's kernel random number generator. To get
457 deterministic behaviour, use 'setStdGen'.
458 -}
459
460 -- |Sets the global random number generator.
461 setStdGen :: StdGen -> IO ()
462 setStdGen sgen = writeIORef theStdGen sgen
463
464 -- |Gets the global random number generator.
465 getStdGen :: IO StdGen
466 getStdGen = readIORef theStdGen
467
468 theStdGen :: IORef StdGen
469 theStdGen = unsafePerformIO $ do
470 rng <- mkStdRNG 0
471 newIORef rng
472
473 -- |Applies 'split' to the current global random generator,
474 -- updates it with one of the results, and returns the other.
475 newStdGen :: IO StdGen
476 newStdGen = atomicModifyIORef theStdGen split
477
478 {- |Uses the supplied function to get a value from the current global
479 random generator, and updates the global generator with the new generator
480 returned by the function. For example, @rollDice@ gets a random integer
481 between 1 and 6:
482
483 > rollDice :: IO Int
484 > rollDice = getStdRandom (randomR (1,6))
485
486 -}
487
488 getStdRandom :: (StdGen -> (a,StdGen)) -> IO a
489 getStdRandom f = atomicModifyIORef theStdGen (swap . f)
490 where swap (v,g) = (g,v)
491
492 {- $references
493
494 1. FW #Burton# Burton and RL Page, /Distributed random number generation/,
495 Journal of Functional Programming, 2(2):203-212, April 1992.
496
497 2. SK #Park# Park, and KW Miller, /Random number generators -
498 good ones are hard to find/, Comm ACM 31(10), Oct 1988, pp1192-1201.
499
500 3. DG #Carta# Carta, /Two fast implementations of the minimal standard
501 random number generator/, Comm ACM, 33(1), Jan 1990, pp87-88.
502
503 4. P #Hellekalek# Hellekalek, /Don\'t trust parallel Monte Carlo/,
504 Department of Mathematics, University of Salzburg,
505 <http://random.mat.sbg.ac.at/~peter/pads98.ps>, 1998.
506
507 5. Pierre #LEcuyer# L'Ecuyer, /Efficient and portable combined random
508 number generators/, Comm ACM, 31(6), Jun 1988, pp742-749.
509
510 The Web site <http://random.mat.sbg.ac.at/> is a great source of information.
511
512 -}