From the GHC 6.6 implementation of the module System.Random

instance Read StdGen where readsPrec _p = \ r -> case try_read r of r@[_] -> r _ -> [stdFromString r] -- because it shouldn't ever fail. where try_read r = do

One can inject invalid state here (out of range) or 0 0.

(s1, r1) <- readDec (dropWhile isSpace r) (s2, r2) <- readDec (dropWhile isSpace r1) return (StdGen s1 s2, r2) {- If we cannot unravel the StdGen from a string, create one based on the string given. -}

I wish this would use a real hash function, for example, SHA-256. It consumes only 6 characters! Since it has about 62 bits of state, it ought to in principle consume 8, though problems with overflow and the 3 multiplier might cause complications.

stdFromString :: String -> (StdGen, String) stdFromString s = (mkStdGen num, rest) where (cs, rest) = splitAt 6 s num = foldl (\a x -> x + 3 * a) 1 (map ord cs) {- | The function 'mkStdGen' provides an alternative way of producing an initial generator, by mapping an 'Int' into a generator. Again, distinct arguments should be likely to produce distinct generators. Programmers may, of course, supply their own instances of 'RandomGen'. -} mkStdGen :: Int -> StdGen -- why not Integer ?

Yes, why not? Then the user could specify lots of initial entropy.

mkStdGen s | s < 0 = mkStdGen (-s) | otherwise = StdGen (s1+1) (s2+1) where (q, s1) = s `divMod` 2147483562 s2 = q `mod` 2147483398

It should be commented that we use one less than the moduli to avoid the initial state being zero.

Look at all these hard coded numbers! Ideally they should be specified exactly once, in factored form, as in L'Ecuyer's paper. The factorization (showing no common factors between the moduli) demonstrates the long period.

Sadly, the function below is not exported.

createStdGen :: Integer -> StdGen createStdGen s | s < 0 = createStdGen (-s) | otherwise = StdGen (fromInteger (s1+1)) (fromInteger (s2+1)) where (q, s1) = s `divMod` 2147483562 s2 = q `mod` 2147483398 -- FIXME: 1/2/3 below should be ** (vs@30082002) XXX

FIXME huh?

...The lower 8 bits are used to create a random Char. (See the definition, especially the `mod` call in randomIvalInteger below.) The lower bits of RNG are suspect.

Even worse the lower 1 bits are used to create a random Bool.

instance Random Char where randomR (a,b) g = case (randomIvalInteger (toInteger (ord a), toInteger (ord b)) g) of (x,g) -> (chr x, g) random g = randomR (minBound,maxBound) g instance Random Bool where randomR (a,b) g = case (randomIvalInteger (toInteger (bool2Int a), toInteger (bool2Int b)) g) of (x, g) -> (int2Bool x, g) where bool2Int False = 0 bool2Int True = 1 int2Bool 0 = False int2Bool _ = True random g = randomR (minBound,maxBound) g -- hah, so you thought you were saving cycles by using Float?

Yes, it sure would be nice.

instance Random Float where random g = randomIvalDouble (0::Double,1) realToFrac g randomR (a,b) g = randomIvalDouble (realToFrac a, realToFrac b) realToFrac g mkStdRNG :: Integer -> IO StdGen mkStdRNG o = do ct <- getCPUTime (TOD sec _) <- getClockTime return (createStdGen (sec * 12345 + ct + o))

Is 12345 the optimal constant? Does it interact badly with the modulus 2147483562?

randomIvalInteger :: (RandomGen g, Num a) => (Integer, Integer) -> g -> (a, g) randomIvalInteger (l,h) rng | l > h = randomIvalInteger (h,l) rng | otherwise = case (f n 1 rng) of (v, rng') -> (fromInteger (l + v `mod` k), rng')

So iLogBase overshoots, then `mod` is used to bring it back into range. The wraparound overshoot makes it no longer uniform.

where k = h - l + 1 b = 2147483561

Shouldn't this be 2147483562? We are creating a big number in base 2147483562.

n = iLogBase b k f 0 acc g = (acc, g) f n acc g = let (x,g') = next g

Might need to subtract 1 from x.

in f (n-1) (fromIntegral x + acc * b) g' randomIvalDouble :: (RandomGen g, Fractional a) => (Double, Double) -> (Double -> a) -> g -> (a, g) randomIvalDouble (l,h) fromDouble rng | l > h = randomIvalDouble (h,l) fromDouble rng | otherwise =

Instead of choosing the pre-scaled range to be the same as the modulus (about 2 billion), we instead do from -2G to 2G, or 4 billion. This will require two calls to the RNG, then throwing away almost all the bits of the second call (which ought to be used for double-precision mantissa!)

case (randomIvalInteger (toInteger (minBound::Int), toInteger (maxBound::Int)) rng) of (x, rng') ->

Consumes only 32-ish bits of randomness to produce 52 bits of double-precision mantissa!

let scaled_x = fromDouble ((l+h)/2) +

The divide by 2 factor is because intRange goes from -2G to 2G.

fromDouble ((h-l) / realToFrac intRange) * fromIntegral (x::Int) in (scaled_x, rng')

intRange returns approximately 4 * 10^{9}.

intRange :: Integer intRange = toInteger (maxBound::Int) - toInteger (minBound::Int)

This function iLogBase has quadratic running time.

iLogBase :: Integer -> Integer -> Integer iLogBase b i = if i < b then 1 else 1 + iLogBase b (i `div` b) stdRange :: (Int,Int)

The range is actually from 1.

stdRange = (0, 2147483562) stdNext :: StdGen -> (Int, StdGen) -- Returns values in the range stdRange stdNext (StdGen s1 s2) = (z', StdGen s1'' s2'') where z' = if z < 1 then z + 2147483562 else z

The value z' = 0 cannot be produced! The statistical independence between s1'' and s2'' is violated, so I'm not sure about the uniformity.

z = s1'' - s2'' k = s1 `quot` 53668 s1' = 40014 * (s1 - k * 53668) - k * 12211 s1'' = if s1' < 0 then s1' + 2147483563 else s1' k' = s2 `quot` 52774 s2' = 40692 * (s2 - k' * 52774) - k' * 3791 s2'' = if s2' < 0 then s2' + 2147483399 else s2' stdSplit :: StdGen -> (StdGen, StdGen) stdSplit std@(StdGen s1 s2) = (left, right) where -- no statistical foundation for this! left = StdGen new_s1 t2 right = StdGen t1 new_s2 new_s1 | s1 == 2147483562 = 1 | otherwise = s1 + 1 new_s2 | s2 == 1 = 2147483398 | otherwise = s2 - 1 StdGen t1 t2 = snd (next std)