Small tweaks to fix Integer-defaulting warnings.
[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 instance Random Float where
349 randomR = randomIvalFrac
350 random rng =
351 case next rng of
352 (x,rng') ->
353 let
354 -- We use 24 bits of randomness corresponding to the 24 bit significand:
355 rand = fromIntegral (mask24 .&. x)
356 :: Float
357 in
358 -- (rand / 2^24, rng')
359 (rand / fromIntegral twoto24, rng')
360 -- Note, encodeFloat is another option, but I'm not seeing slightly
361 -- worse performance with the following [2011.06.25]:
362 -- (encodeFloat rand (-24), rng')
363 where
364 mask24 :: Int
365 -- mask24 = 2^24 - 1
366
367 mask24 = twoto24 - 1
368 -- RRN: Note, in my tests [2011.06.25] this worked as well as using Data.Bit:
369 twoto24 = (2::Int) ^ (24::Int)
370
371
372 instance Random CFloat where
373 randomR = randomIvalFrac
374 random rng = case random rng of
375 (x,rng') -> (realToFrac (x::Float), rng')
376
377 instance Random CDouble where
378 randomR = randomIvalFrac
379 random = randomFrac
380
381
382 mkStdRNG :: Integer -> IO StdGen
383 mkStdRNG o = do
384 ct <- getCPUTime
385 (sec, psec) <- getTime
386 return (createStdGen (sec * 12345 + psec + ct + o))
387
388 randomBounded :: (RandomGen g, Random a, Bounded a) => g -> (a, g)
389 randomBounded = randomR (minBound, maxBound)
390
391 -- The two integer functions below take an [inclusive,inclusive] range.
392 randomIvalIntegral :: (RandomGen g, Integral a) => (a, a) -> g -> (a, g)
393 randomIvalIntegral (l,h) = randomIvalInteger (toInteger l, toInteger h)
394
395 randomIvalInteger :: (RandomGen g, Num a) => (Integer, Integer) -> g -> (a, g)
396 randomIvalInteger (l,h) rng
397 | l > h = randomIvalInteger (h,l) rng
398 | otherwise = case (f n 1 rng) of (v, rng') -> (fromInteger (l + v `mod` k), rng')
399 where
400 k = h - l + 1
401 b = 2147483561
402 n = iLogBase b k
403
404 f 0 acc g = (acc, g)
405 f n' acc g =
406 let
407 (x,g') = next g
408 in
409 f (n' - 1) (fromIntegral x + acc * b) g'
410
411 -- The continuous functions on the other hand take an [inclusive,exclusive) range.
412 randomFrac :: (RandomGen g, Fractional a) => g -> (a, g)
413 randomFrac = randomIvalDouble (0::Double,1) realToFrac
414
415 randomIvalFrac :: (RandomGen g, Real a, Fractional b) => (a,a) -> g -> (b, g)
416 randomIvalFrac (a,b) = randomIvalDouble (realToFrac a, realToFrac b) realToFrac
417
418 randomIvalDouble :: (RandomGen g, Fractional a) => (Double, Double) -> (Double -> a) -> g -> (a, g)
419 randomIvalDouble (l,h) fromDouble rng
420 | l > h = randomIvalDouble (h,l) fromDouble rng
421 | otherwise =
422 case (randomIvalInteger (toInteger (minBound::Int32), toInteger (maxBound::Int32)) rng) of
423 (x, rng') ->
424 let
425 scaled_x =
426 fromDouble ((l+h)/2) +
427 fromDouble ((h-l) / realToFrac int32Count) *
428 fromIntegral (x::Int32)
429 in
430 (scaled_x, rng')
431
432 int32Count :: Integer
433 int32Count = toInteger (maxBound::Int32) - toInteger (minBound::Int32) + 1
434
435 iLogBase :: Integer -> Integer -> Integer
436 iLogBase b i = if i < b then 1 else 1 + iLogBase b (i `div` b)
437
438 stdRange :: (Int,Int)
439 stdRange = (0, 2147483562)
440
441 stdNext :: StdGen -> (Int, StdGen)
442 -- Returns values in the range stdRange
443 stdNext (StdGen s1 s2) = (fromIntegral z', StdGen s1'' s2'')
444 where z' = if z < 1 then z + 2147483562 else z
445 z = s1'' - s2''
446
447 k = s1 `quot` 53668
448 s1' = 40014 * (s1 - k * 53668) - k * 12211
449 s1'' = if s1' < 0 then s1' + 2147483563 else s1'
450
451 k' = s2 `quot` 52774
452 s2' = 40692 * (s2 - k' * 52774) - k' * 3791
453 s2'' = if s2' < 0 then s2' + 2147483399 else s2'
454
455 stdSplit :: StdGen -> (StdGen, StdGen)
456 stdSplit std@(StdGen s1 s2)
457 = (left, right)
458 where
459 -- no statistical foundation for this!
460 left = StdGen new_s1 t2
461 right = StdGen t1 new_s2
462
463 new_s1 | s1 == 2147483562 = 1
464 | otherwise = s1 + 1
465
466 new_s2 | s2 == 1 = 2147483398
467 | otherwise = s2 - 1
468
469 StdGen t1 t2 = snd (next std)
470
471 -- The global random number generator
472
473 {- $globalrng #globalrng#
474
475 There is a single, implicit, global random number generator of type
476 'StdGen', held in some global variable maintained by the 'IO' monad. It is
477 initialised automatically in some system-dependent fashion, for example, by
478 using the time of day, or Linux's kernel random number generator. To get
479 deterministic behaviour, use 'setStdGen'.
480 -}
481
482 -- |Sets the global random number generator.
483 setStdGen :: StdGen -> IO ()
484 setStdGen sgen = writeIORef theStdGen sgen
485
486 -- |Gets the global random number generator.
487 getStdGen :: IO StdGen
488 getStdGen = readIORef theStdGen
489
490 theStdGen :: IORef StdGen
491 theStdGen = unsafePerformIO $ do
492 rng <- mkStdRNG 0
493 newIORef rng
494
495 -- |Applies 'split' to the current global random generator,
496 -- updates it with one of the results, and returns the other.
497 newStdGen :: IO StdGen
498 newStdGen = atomicModifyIORef theStdGen split
499
500 {- |Uses the supplied function to get a value from the current global
501 random generator, and updates the global generator with the new generator
502 returned by the function. For example, @rollDice@ gets a random integer
503 between 1 and 6:
504
505 > rollDice :: IO Int
506 > rollDice = getStdRandom (randomR (1,6))
507
508 -}
509
510 getStdRandom :: (StdGen -> (a,StdGen)) -> IO a
511 getStdRandom f = atomicModifyIORef theStdGen (swap . f)
512 where swap (v,g) = (g,v)
513
514 {- $references
515
516 1. FW #Burton# Burton and RL Page, /Distributed random number generation/,
517 Journal of Functional Programming, 2(2):203-212, April 1992.
518
519 2. SK #Park# Park, and KW Miller, /Random number generators -
520 good ones are hard to find/, Comm ACM 31(10), Oct 1988, pp1192-1201.
521
522 3. DG #Carta# Carta, /Two fast implementations of the minimal standard
523 random number generator/, Comm ACM, 33(1), Jan 1990, pp87-88.
524
525 4. P #Hellekalek# Hellekalek, /Don\'t trust parallel Monte Carlo/,
526 Department of Mathematics, University of Salzburg,
527 <http://random.mat.sbg.ac.at/~peter/pads98.ps>, 1998.
528
529 5. Pierre #LEcuyer# L'Ecuyer, /Efficient and portable combined random
530 number generators/, Comm ACM, 31(6), Jun 1988, pp742-749.
531
532 The Web site <http://random.mat.sbg.ac.at/> is a great source of information.
533
534 -}