2 % (c) The GRASP/AQUA Project, Glasgow University, 1995-99
6 The June 1988 (v31 #6) issue of the Communications of the ACM has an
7 article by Pierre L'Ecuyer called, "Efficient and Portable Combined
8 Random Number Generators". Here is the Portable Combined Generator of
9 L'Ecuyer for 32-bit computers. It has a period of roughly 2.30584e18.
11 Transliterator: Lennart Augustsson
13 sof 1/99 - code brought (kicking and screaming) into the new Random
19 RandomGen(next, split)
22 , Random ( random, randomR,
32 import PrelGHC ( RealWorld )
33 import PrelNum ( fromInt )
34 import PrelShow ( showSignedInt, showSpace )
35 import PrelRead ( readDec )
36 import PrelIOBase ( unsafePerformIO, stToIO )
37 import PrelArr ( MutableVar, newVar, readVar, writeVar )
38 import PrelReal ( toInt )
39 import PrelFloat ( float2Double, double2Float )
40 import Time ( getClockTime, ClockTime(..) )
43 import privileged Prelude
52 import CPUTime ( getCPUTime )
53 import Char ( isSpace, chr, ord )
57 class RandomGen g where
67 instance RandomGen StdGen where
72 instance Show StdGen where
73 showsPrec p (StdGen s1 s2) =
78 instance Show StdGen where
79 showsPrec p (StdGen s1 s2) =
85 instance Read StdGen where
89 _ -> [stdFromString r] -- because it shouldn't ever fail.
92 (s1, r1) <- readDec (dropWhile isSpace r)
93 (s2, r2) <- readDec (dropWhile isSpace r1)
94 return (StdGen s1 s2, r2)
97 If we cannot unravel the StdGen from a string, create
98 one based on the string given.
100 stdFromString :: String -> (StdGen, String)
101 stdFromString s = (mkStdGen num, rest)
102 where (cs, rest) = splitAt 6 s
103 num = foldl (\a x -> x + 3 * a) 1 (map ord cs)
107 mkStdGen :: Int -> StdGen -- why not Integer ?
109 | s < 0 = mkStdGen (-s)
110 | otherwise = StdGen (s1+1) (s2+1)
112 (q, s1) = s `divMod` 2147483562
113 s2 = q `mod` 2147483398
115 createStdGen :: Integer -> StdGen
117 | s < 0 = createStdGen (-s)
118 | otherwise = StdGen (toInt (s1+1)) (toInt (s2+1))
120 (q, s1) = s `divMod` 2147483562
121 s2 = q `mod` 2147483398
125 The class definition - see library report for details.
129 -- Minimal complete definition: random and randomR
130 random :: RandomGen g => g -> (a, g)
131 randomR :: RandomGen g => (a,a) -> g -> (a,g)
133 randoms :: RandomGen g => g -> [a]
134 randoms g = x : randoms g' where (x,g') = random g
136 randomRs :: RandomGen g => (a,a) -> g -> [a]
137 randomRs ival g = x : randomRs ival g' where (x,g') = randomR ival g
140 randomIO = getStdRandom random
142 randomRIO :: (a,a) -> IO a
143 randomRIO range = getStdRandom (randomR range)
147 instance Random Int where
148 randomR (a,b) g = randomIvalInteger (toInteger a, toInteger b) g
149 random g = randomR (minBound,maxBound) g
151 instance Random Char where
153 case (randomIvalInteger (toInteger (ord a), toInteger (ord b)) g) of
155 random g = randomR (minBound,maxBound) g
157 instance Random Bool where
159 case (randomIvalInteger (toInteger (bool2Int a), toInteger (bool2Int b)) g) of
160 (x, g) -> (int2Bool x, g)
168 random g = randomR (minBound,maxBound) g
170 instance Random Integer where
171 randomR ival g = randomIvalInteger ival g
172 random g = randomR (toInteger (minBound::Int), toInteger (maxBound::Int)) g
174 instance Random Double where
175 randomR ival g = randomIvalDouble ival id g
176 random g = randomR (0::Double,1) g
178 -- hah, so you thought you were saving cycles by using Float?
181 instance Random Float where
182 random g = randomIvalDouble (0::Double,1) realToFrac g
183 randomR (a,b) g = randomIvalDouble (realToFrac a, realToFrac b) realToFrac g
185 instance Random Float where
186 randomR (a,b) g = randomIvalDouble (float2Double a, float2Double b) double2Float g
187 random g = randomIvalDouble (0::Double,1) double2Float g
195 mkStdRNG :: Integer -> IO StdGen
198 return (createStdGen (ct + o))
200 mkStdRNG :: Integer -> IO StdGen
203 (TOD sec _) <- getClockTime
204 return (createStdGen (sec * 12345 + ct + o))
207 randomIvalInteger :: (RandomGen g, Num a) => (Integer, Integer) -> g -> (a, g)
208 randomIvalInteger (l,h) rng
209 | l > h = randomIvalInteger (h,l) rng
210 | otherwise = case (f n 1 rng) of (v, rng') -> (fromInteger (l + v `mod` k), rng')
221 f (n-1) (fromInt x + acc * b) g'
223 randomIvalDouble :: (RandomGen g, Fractional a) => (Double, Double) -> (Double -> a) -> g -> (a, g)
224 randomIvalDouble (l,h) fromDouble rng
225 | l > h = randomIvalDouble (h,l) fromDouble rng
227 case (randomIvalInteger (toInteger (minBound::Int), toInteger (maxBound::Int)) rng) of
231 fromDouble ((l+h)/2) +
232 fromDouble ((h-l) / realToFrac intRange) *
233 fromIntegral (x::Int)
238 intRange = toInteger (maxBound::Int) - toInteger (minBound::Int)
240 iLogBase :: Integer -> Integer -> Integer
241 iLogBase b i = if i < b then 1 else 1 + iLogBase b (i `div` b)
243 stdNext :: StdGen -> (Int, StdGen)
244 stdNext (StdGen s1 s2) = (z', StdGen s1'' s2'')
245 where z' = if z < 1 then z + 2147483562 else z
249 s1' = 40014 * (s1 - k * 53668) - k * 12211
250 s1'' = if s1' < 0 then s1' + 2147483563 else s1'
253 s2' = 40692 * (s2 - k' * 52774) - k' * 3791
254 s2'' = if s2' < 0 then s2' + 2147483399 else s2'
256 stdSplit :: StdGen -> (StdGen, StdGen)
257 stdSplit std@(StdGen s1 s2)
260 -- no statistical foundation for this!
261 left = StdGen new_s1 t2
262 right = StdGen t1 new_s2
264 new_s1 | s1 == 2147483562 = 1
267 new_s2 | s2 == 1 = 2147483398
270 StdGen t1 t2 = snd (next std)
277 setStdGen :: StdGen -> IO ()
278 setStdGen sgen = writeIORef theStdGen sgen
280 getStdGen :: IO StdGen
281 getStdGen = readIORef theStdGen
283 theStdGen :: IORef StdGen
284 theStdGen = unsafePerformIO (newIORef (createStdGen 0))
288 global_rng :: MutableVar RealWorld StdGen
289 global_rng = unsafePerformIO $ do
293 setStdGen :: StdGen -> IO ()
294 setStdGen sgen = stToIO (writeVar global_rng sgen)
296 getStdGen :: IO StdGen
297 getStdGen = stToIO (readVar global_rng)
302 newStdGen :: IO StdGen
305 let (a,b) = split rng
309 getStdRandom :: (StdGen -> (a,StdGen)) -> IO a
312 let (v, new_rng) = f rng