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