[project @ 2004-01-05 10:07:52 by ross]
[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 -- Random numbers.
12 --
13 -----------------------------------------------------------------------------
14
15 module System.Random
16 (
17
18 -- $intro
19
20 -- * The 'RandomGen' class, and the 'StdGen' generator
21
22 RandomGen(next, split, genRange)
23 , StdGen
24 , mkStdGen
25
26 -- * The 'Random' class
27 , Random ( random, randomR,
28 randoms, randomRs,
29 randomIO, randomRIO )
30
31 -- * The global random number generator
32
33 -- $globalrng
34
35 , getStdRandom
36 , getStdGen
37 , setStdGen
38 , newStdGen
39
40 -- * References
41 -- $references
42
43 ) where
44
45 import Prelude
46
47 #ifdef __NHC__
48 import CPUTime ( getCPUTime )
49 import Foreign.Ptr ( Ptr, nullPtr )
50 #else
51 import System.CPUTime ( getCPUTime )
52 import System.Time ( getClockTime, ClockTime(..) )
53 #endif
54 import Data.Char ( isSpace, chr, ord )
55 import System.IO.Unsafe ( unsafePerformIO )
56 import Data.IORef
57 import Numeric ( readDec )
58
59 -- The standard nhc98 implementation of Time.ClockTime does not match
60 -- the extended one expected in this module, so we lash-up a quick
61 -- replacement here.
62 #ifdef __NHC__
63 data ClockTime = TOD Integer ()
64 foreign import ccall "time.h time" readtime :: Ptr () -> IO Int
65 getClockTime :: IO ClockTime
66 getClockTime = do t <- readtime nullPtr; return (TOD (toInteger t) ())
67 #endif
68
69 {- $intro
70
71 This library deals with the common task of pseudo-random
72 number generation. The library makes it possible to generate
73 repeatable results, by starting with a specified initial random
74 number generator; or to get different results on each run by using the
75 system-initialised generator, or by supplying a seed from some other
76 source.
77
78 The library is split into two layers:
79
80 * A core /random number generator/ provides a supply of bits. The class
81 'RandomGen' provides a common interface to such generators.
82
83 * The class 'Random' provides a way to extract particular values from
84 a random number generator. For example, the 'Float' instance of 'Random'
85 allows one to generate random values of type 'Float'.
86
87 [Comment found in this file when merging with Library Report:]
88
89 The June 1988 (v31 \#6) issue of the Communications of the ACM has an
90 article by Pierre L'Ecuyer called, /Efficient and Portable Combined
91 Random Number Generators/. Here is the Portable Combined Generator of
92 L'Ecuyer for 32-bit computers. It has a period of roughly 2.30584e18.
93
94 Transliterator: Lennart Augustsson
95
96 -}
97
98 -- | The class 'RandomGen' provides a common interface to random number
99 -- generators.
100
101 class RandomGen g where
102
103 -- |The 'next' operation allows one to extract at least 30 bits (one 'Int''s
104 -- worth) from the generator, returning a new generator as well. The
105 -- integer returned may be positive or negative.
106 next :: g -> (Int, g)
107
108 -- |The 'split' operation allows one to obtain two distinct random number
109 -- generators. This is very useful in functional programs (for example, when
110 -- passing a random number generator down to recursive calls), but very
111 -- little work has been done on statistically robust implementations of
112 -- 'split' ([1,4] are the only examples we know of).
113 split :: g -> (g, g)
114
115 -- |The 'genRange' operation yields the range of values returned by
116 -- the generator.
117 --
118 -- It is required that:
119 --
120 -- * If @(a,b) = 'genRange' g@, then @a < b@.
121 --
122 -- * 'genRange' is not strict.
123 --
124 -- The second condition ensures that 'genRange' cannot examine its
125 -- argument, and hence the value it returns can be determined only by the
126 -- instance of 'RandomGen'. That in turn allows an implementation to make
127 -- a single call to 'genRange' to establish a generator's range, without
128 -- being concerned that the generator returned by (say) 'next' might have
129 -- a different range to the generator passed to 'next'.
130 genRange :: g -> (Int,Int)
131
132 -- default method
133 genRange g = (minBound,maxBound)
134
135 {- |The "System.Random" library provides one instance of 'RandomGen', the
136 abstract data type 'StdGen'.
137
138 The 'StdGen' instance of 'RandomGen' has a 'genRange' of at least 30 bits.
139
140 The result of repeatedly using 'next' should be at least as statistically
141 robust as the /Minimal Standard Random Number Generator/ described by
142 ["System.Random\#Park", "System.Random\#Carta"].
143 Until more is known about implementations of 'split', all we require is
144 that 'split' deliver generators that are (a) not identical and
145 (b) independently robust in the sense just given.
146
147 The 'Show'\/'Read' instances of 'StdGen' provide a primitive way to save the
148 state of a random number generator.
149 It is required that @'read' ('show' g) == g@.
150
151 In addition, 'read' may be used to map an arbitrary string (not necessarily one
152 produced by 'show') onto a value of type 'StdGen'. In general, the 'read'
153 instance of 'StdGen' has the following properties:
154
155 * It guarantees to succeed on any string.
156
157 * It guarantees to consume only a finite portion of the string.
158
159 * Different argument strings are likely to result in different results.
160
161 -}
162
163 data StdGen
164 = StdGen Int Int
165
166 instance RandomGen StdGen where
167 next = stdNext
168 split = stdSplit
169
170 instance Show StdGen where
171 showsPrec p (StdGen s1 s2) =
172 showsPrec p s1 .
173 showChar ' ' .
174 showsPrec p s2
175
176 instance Read StdGen where
177 readsPrec _p = \ r ->
178 case try_read r of
179 r@[_] -> r
180 _ -> [stdFromString r] -- because it shouldn't ever fail.
181 where
182 try_read r = do
183 (s1, r1) <- readDec (dropWhile isSpace r)
184 (s2, r2) <- readDec (dropWhile isSpace r1)
185 return (StdGen s1 s2, r2)
186
187 {-
188 If we cannot unravel the StdGen from a string, create
189 one based on the string given.
190 -}
191 stdFromString :: String -> (StdGen, String)
192 stdFromString s = (mkStdGen num, rest)
193 where (cs, rest) = splitAt 6 s
194 num = foldl (\a x -> x + 3 * a) 1 (map ord cs)
195
196
197 {- |
198 The function 'mkStdGen' provides an alternative way of producing an initial
199 generator, by mapping an 'Int' into a generator. Again, distinct arguments
200 should be likely to produce distinct generators.
201
202 Programmers may, of course, supply their own instances of 'RandomGen'.
203 -}
204 mkStdGen :: Int -> StdGen -- why not Integer ?
205 mkStdGen s
206 | s < 0 = mkStdGen (-s)
207 | otherwise = StdGen (s1+1) (s2+1)
208 where
209 (q, s1) = s `divMod` 2147483562
210 s2 = q `mod` 2147483398
211
212 createStdGen :: Integer -> StdGen
213 createStdGen s
214 | s < 0 = createStdGen (-s)
215 | otherwise = StdGen (fromInteger (s1+1)) (fromInteger (s2+1))
216 where
217 (q, s1) = s `divMod` 2147483562
218 s2 = q `mod` 2147483398
219
220 -- FIXME: 1/2/3 below should be ** (vs@30082002) XXX
221
222 {- |
223 With a source of random number supply in hand, the 'Random' class allows the
224 programmer to extract random values of a variety of types.
225
226 Minimal complete definition: 'randomR' and 'random'.
227
228 -}
229
230 class Random a where
231 -- | Takes a range /(lo,hi)/ and a random number generator
232 -- /g/, and returns a random value uniformly distributed in the closed
233 -- interval /[lo,hi]/, together with a new generator. It is unspecified
234 -- what happens if /lo>hi/. For continuous types there is no requirement
235 -- that the values /lo/ and /hi/ are ever produced, but they may be,
236 -- depending on the implementation and the interval.
237 randomR :: RandomGen g => (a,a) -> g -> (a,g)
238
239 -- | The same as 'randomR', but using a default range determined by the type:
240 --
241 -- * For bounded types (instances of 'Bounded', such as 'Char'),
242 -- the range is normally the whole type.
243 --
244 -- * For fractional types, the range is normally the semi-closed interval
245 -- @[0,1)@.
246 --
247 -- * For 'Integer', the range is (arbitrarily) the range of 'Int'.
248 random :: RandomGen g => g -> (a, g)
249
250 -- | Plural variant of 'randomR', producing an infinite list of
251 -- random values instead of returning a new generator.
252 randomRs :: RandomGen g => (a,a) -> g -> [a]
253 randomRs ival g = x : randomRs ival g' where (x,g') = randomR ival g
254
255 -- | Plural variant of 'random', producing an infinite list of
256 -- random values instead of returning a new generator.
257 randoms :: RandomGen g => g -> [a]
258 randoms g = (\(x,g') -> x : randoms g') (random g)
259
260 -- | A variant of 'randomR' that uses the global random number generator
261 -- (see "System.Random#globalrng").
262 randomRIO :: (a,a) -> IO a
263 randomRIO range = getStdRandom (randomR range)
264
265 -- | A variant of 'random' that uses the global random number generator
266 -- (see "System.Random#globalrng").
267 randomIO :: IO a
268 randomIO = getStdRandom random
269
270
271 instance Random Int where
272 randomR (a,b) g = randomIvalInteger (toInteger a, toInteger b) g
273 random g = randomR (minBound,maxBound) g
274
275 instance Random Char where
276 randomR (a,b) g =
277 case (randomIvalInteger (toInteger (ord a), toInteger (ord b)) g) of
278 (x,g) -> (chr x, g)
279 random g = randomR (minBound,maxBound) g
280
281 instance Random Bool where
282 randomR (a,b) g =
283 case (randomIvalInteger (toInteger (bool2Int a), toInteger (bool2Int b)) g) of
284 (x, g) -> (int2Bool x, g)
285 where
286 bool2Int False = 0
287 bool2Int True = 1
288
289 int2Bool 0 = False
290 int2Bool _ = True
291
292 random g = randomR (minBound,maxBound) g
293
294 instance Random Integer where
295 randomR ival g = randomIvalInteger ival g
296 random g = randomR (toInteger (minBound::Int), toInteger (maxBound::Int)) g
297
298 instance Random Double where
299 randomR ival g = randomIvalDouble ival id g
300 random g = randomR (0::Double,1) g
301
302 -- hah, so you thought you were saving cycles by using Float?
303 instance Random Float where
304 random g = randomIvalDouble (0::Double,1) realToFrac g
305 randomR (a,b) g = randomIvalDouble (realToFrac a, realToFrac b) realToFrac g
306
307 mkStdRNG :: Integer -> IO StdGen
308 mkStdRNG o = do
309 ct <- getCPUTime
310 (TOD sec _) <- getClockTime
311 return (createStdGen (sec * 12345 + ct + o))
312
313 randomIvalInteger :: (RandomGen g, Num a) => (Integer, Integer) -> g -> (a, g)
314 randomIvalInteger (l,h) rng
315 | l > h = randomIvalInteger (h,l) rng
316 | otherwise = case (f n 1 rng) of (v, rng') -> (fromInteger (l + v `mod` k), rng')
317 where
318 k = h - l + 1
319 b = 2147483561
320 n = iLogBase b k
321
322 f 0 acc g = (acc, g)
323 f n acc g =
324 let
325 (x,g') = next g
326 in
327 f (n-1) (fromIntegral x + acc * b) g'
328
329 randomIvalDouble :: (RandomGen g, Fractional a) => (Double, Double) -> (Double -> a) -> g -> (a, g)
330 randomIvalDouble (l,h) fromDouble rng
331 | l > h = randomIvalDouble (h,l) fromDouble rng
332 | otherwise =
333 case (randomIvalInteger (toInteger (minBound::Int), toInteger (maxBound::Int)) rng) of
334 (x, rng') ->
335 let
336 scaled_x =
337 fromDouble ((l+h)/2) +
338 fromDouble ((h-l) / realToFrac intRange) *
339 fromIntegral (x::Int)
340 in
341 (scaled_x, rng')
342
343 intRange :: Integer
344 intRange = toInteger (maxBound::Int) - toInteger (minBound::Int)
345
346 iLogBase :: Integer -> Integer -> Integer
347 iLogBase b i = if i < b then 1 else 1 + iLogBase b (i `div` b)
348
349 stdNext :: StdGen -> (Int, StdGen)
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 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, /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 The Web site <http://random.mat.sbg.ac.at/> is a great source of information.
444
445 -}