Random number quality with given seed in haskell

Question

I'm playing with procedural genaration in a small ascii game and I encountered this issue with random numbers in haskell. The basic idea is to provide a random number seeded with (x,y) of some part of game world to for example decide if there's a tree there or not (this guy explains it nicely)

This is what I get when trying a different seed for each generation:

randomFromSeed :: Int -> Int -> Int
randomFromSeed max seed = fst (randomR (0, max - 1) (mkStdGen seed))

Prelude> map (randomFromSeed 10) [1..20]
[5,9,3,7,1,5,9,3,7,1,5,9,3,7,1,5,9,3,7,1]

It clearly has a period of 5, but on the other hand on the mkStdGen docs it says:

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.

So how come, there seem to be only 5 distinct generators coming?

How can I get those to be truly random when given different seeds?

Edit For some weird reason using bigger numbers makes it better:

Prelude> let mult = 1000000 in map (randomFromSeed 10) [0,mult .. 20*mult]
[3,7,0,6,9,2,8,1,4,0,3,9,2,5,1,4,7,3,6,9,5]

Answer 1

So how come, there seems to be only 5 distinct generators coming ?

It is an illusion to think there are just 5 generators. If you print the second numbers of each sequence instead of the first ones, you get this:

random2ndFromSeed :: Int -> Int -> Int
random2ndFromSeed max seed =
    let   g0       = mkStdGen seed
          (v1, g1) = randomR (0, max - 1) g0
          (v2, g2) = randomR (0, max - 1) g1
    in   v2

 λ> 
 λ> map  (random2ndFromSeed 10)  [1..40]
[6,9,3,8,1,4,8,3,6,9,3,8,1,4,8,3,6,9,3,8,1,4,8,3,6,9,3,8,1,4,8,3,6,9,3,8,1,4,8,3]
 λ> 

So the periodicity appears then to be 8 instead of 5 !

A way to get rid of the apparent problem is to replace the standard generator by the Threefish one, which is of more recent design and has better statistical properties. Or you can also use pcg-random as mentioned by Dave Compton.

import  System.Random.TF

tfRandomFromSeed :: Int -> Int -> Int
tfRandomFromSeed max seed = let   g0 = mkTFGen seed
                            in    fst $ randomR (0, max - 1) g0

 λ> 
 λ> map  (tfRandomFromSeed 10)  [1..40]
[4,5,6,7,5,3,3,0,0,4,2,8,0,4,1,0,0,1,3,5,6,4,3,6,4,0,3,6,4,0,2,4,5,9,7,3,8,5,2,4]
 λ> 

More generally, the appearance of randomness is supposed to come from repeated applications of the generator next function. Here, the function is only applied once per seed/sequence, so there is no claim to randomness.

How to create a persistent 2D random field

From the comments, the actual need is of a “random” function of a point in 2D space. If the player, after some random walk, returns to some already visited point, it is expected to find the same random value as before, and this without memorizing previous random values.

And to achieve this in a way that gets us some guarantees about the statistical properties of the random values, we need to do it with a single seed and a single random sequence; for thats' what our applied mathematicians are testing.

We need two things to produce such a persistent two-dimensional random field:

1. a way to map a 2D point to a single natural number, used as an offset into a random sequence
2. a random number generator with efficient access to arbitrary points of each sequence

Mapping 2D points to natural numbers

This can be done for example by leveraging the Cantor Pairing Function from elementary set theory.

We can use this code:

-- limited to first quadrant, x >= 0 and y >= 0:
cantor1 :: Int -> Int -> Int
cantor1 x y = y + (let s = x + y  in  div  (s * (s+1))  2)

-- for all 4 quadrants:
cantor :: (Int, Int) -> Int
cantor (x,y) =
    let quadrant
          | x >= 0  &&  y >= 0   =  0
          | x <  0  &&  y >= 0   =  1
          | x <  0  &&  y  < 0   =  2  
          | x >= 0  &&  y <  0   =  3
          | otherwise            =  error  "cantor: internal error #1"
        cant1
          | x >= 0  &&  y >= 0   =  cantor1     x      y
          | x <  0  &&  y >= 0   =  cantor1  (-1-x)    y
          | x <  0  &&  y  < 0   =  cantor1  (-1-x)  (-1-y)
          | x >= 0  &&  y <  0   =  cantor1     x    (-1-y)
          | otherwise            =  error  "cantor: internal error #2"
    in
         4*cant1 + quadrant

Arranging arbitrary access

With this preliminary step out of the way, we have to recognize that the regular Haskell random number generation API is not well suited to the task at hand.

The API provides sequential access to the random sequence thru the next function. But there is no arbitrary access, such as provided in the C++ random library by the discard function. And the classic monadic style using the MonadRandom interface is all about sequential access. It's basically like a state monad.

Furthermore, with some random number generators, efficient access to an arbitrary point of the sequence is simply impossible. In such a case, the C++ discard function just uses costly single stepping to get to the wanted point.

Fortunately, there is a Haskell implementation of Pierre L'Ecuyer et al MRG32k3a random number generator.

With MRG32k3a, arbitrary access into the random sequence boils down to exponentiation of small matrices in 2 Galois fields. Thanks to the ancient and revered Indian exponentiation algorithm, this can be done in O(log n) time.

The MRG32k3a code in github does not provide a full Haskell style interface, such as a RandomGen instance, so we have to add a bit of wrapper code around it.

First, we need some import clauses:

import  System.Random
import  System.Random.TF
import qualified  Data.List           as  L
import qualified  Text.Printf         as  TP
import qualified  Data.Text           as  TL
import qualified  Data.ByteString     as  BS
import qualified  Data.Text.Encoding  as  TSE
import qualified  Crypto.Hash.SHA256  as  SHA
import qualified  System.Random.MRG32K3A.Simple as MRG

and then the wrapper code itself:

newtype MRGen = MRGen MRG.State  -- wrapper type for MRG32k3a generator
                deriving  Show

instance RandomGen  MRGen  where
    genRange = let  mrg32k3a_m1 = ((2::Integer)^32 - 209)
               in   const  (0::Int, fromIntegral (mrg32k3a_m1 - 1))

    next (MRGen g0) = let  (v, g1) = MRG.next g0
                      in   ((fromIntegral v)::Int, MRGen g1)

    split (MRGen g0) = let  g1 = MRG.advance ((2::Integer)^96) g0
                       in   (MRGen g0, MRGen g1) 

mkMRGen :: Int -> MRGen
mkMRGen userSeed = let  longSeed = hashSeed userSeed
                        g0       =  MRG.seed longSeed
                   in   MRGen g0

ranSeek :: MRGen -> Integer -> MRGen
ranSeek (MRGen g0) count =  let  g1 = (MRG.advance count g0)  in   MRGen g1

hashSeed :: Int -> Integer
hashSeed userSeed =
    let str   = "MRG32k3a:" ++ (TP.printf "0x%x" userSeed)
        bytes =  (TSE.encodeUtf8 . TL.pack) $ str
        ints  = (map (fromIntegral) $ BS.unpack (SHA.hash bytes)) :: [Integer]
    in
        L.foldl'  (\acc d -> acc*256 + d)  0  (take 20 ints)

Function mkMRGen is similar to mkStdGen. Arbitrary access into the random sequence is provided by function ranSeek :: MRGen -> Integer -> MRGen in O(log n) time.

Side note: I am rehashing the user-provided seed in mkMRGen. this is because the github package uses its seed as just an offset into the random sequence. So in order to avoid the risk of sequence overlap for small user seeds, I need to generate a large number from the user seed.

Thanks to our RandomGen instance, we have access to usual functions such as random :: RandomGen g => g -> (a, g). For example, we can generate a 2D random field of type Double from a simple Int seed like this:

randomDoubleField :: Int -> (Int, Int) -> Double
randomDoubleField userSeed (x,y) =
    let  k  = 1  -- number of needed random values per plane point
         g0 = mkMRGen userSeed
         g1 = ranSeek  g0  (fromIntegral (k * cantor (x,y)))
    in   fst (random g1)

Now that we have that little toolkit, we can write a small test program, drawing some random landscape for a neighboorhood of point zero, with one character per 2D point.

Say, character 't' stands for one type of tree and 'T' for another type of tree. The absence of tree is denoted by a minus sign.

Main program:

randomCharField :: Int -> (Int, Int) -> Char
randomCharField  userSeed  (x,y) =
    let  n = floor (8.0 * randomDoubleField userSeed (x,y) )
    in   "------tT"  !!  n


rowString :: Int -> Int -> Int -> String
rowString userSeed size y =
               let  xRange = [(-size) .. size]
               in   map  (randomCharField userSeed)  [ (x,y) | x <- xRange ]


main = do
    let  userSeed = 42
         size     = 6
         yRange   = [(-size) .. size]
    mapM_  (putStrLn . (rowString userSeed size))  yRange

Program output:

--t-T----TT-t
------t-----T
-T--T--T-----
--t-T--tTTT--
--T--t---T---
t-Tt------t--
-T-----t-T---
-T-t-t----T--
tT-tT---tT--t
---TTt---t---
-------T---t-
--t---------t
-tT-t---t----

Optimization note: If performance is a concern, you probably want to move the (mkMRGen userSeed) computation out of the loops.

Answer 2

You can avoid the unexpected behavior you see by using pcg-random instead of random :

import System.Random.PCG 
import Control.Monad.ST

randomFromSeed :: Int -> Int -> Int
randomFromSeed max seed = runST $ do
  g <- initialize (fromIntegral seed) 0
  uniformR (0, max - 1) g

main :: IO ()
main = print $ map (randomFromSeed 10) [1..20]

pcg-random has other nice properties too.