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