Posts Tagged ‘generator’

Programming Praxis – Two Random Exercises

August 17, 2012

In today’s Programming Praxis exercise, our goal is to use convert two random number generators to a different range. Specifically, we have to convert an rng that produces numbers in the closed interval [1,3] to one that has the interval [1,9] and one with interval [1,5] to [1,7]. In both cases the resulting rngs must maintain an equal distribution across the interval. Let’s get started, shall we?

Some imports:

import Control.Applicative
import Control.Monad
import Data.List
import System.Random

These are the two functions we can use as imports.

rand3, rand5 :: IO Int
rand3 = randomRIO (1,3)
rand5 = randomRIO (1,5)

To test, we use a distribution function to see how many of each element occur in a sequence.

dist :: (Eq a, Ord a) => [a] -> [(a, Int)]
dist = map (\x -> (head x, length x)) . group . sort

Let’s start with the easiest one. Since 9 = 3 * 3, all we need to do scale up one input by a factor of 3 and “interpolate” using a second one. The reason for taking a source parameter is to facilitate testing: you can either feed it an IO Int for a random result, or you can give it a [Int] to see all possible results. We’ll get to why the testing is important in a bit.

rand9 :: (Applicative f, Num a, Enum a) => f a -> f a
rand9 src = (+) . (*3) . pred <$> src <*> src

The second one is more complicated, since 5*5 is not a multiple of 7. The correct solution is to do the same thing as above using a subrange that is a multiple of 7, i.e. 21 of the 25 numbers and rerolling when one of the four values outside the accepted range is chosen. The problem with this approach is that it isn’t guaranteed to terminate, since there is nothing preventing a true random number generator from generating the number 23 a few trillion times in a row. Granted, we’re only using pseudo-random number generators, so it probably will never come up. Here’s the implementation:

rand7 :: (Applicative m, Monad m, Integral a) => m a -> m a
rand7 src = (+) . (*5) . pred <$> src <*> src >>= \n -> 
    if n > 21 then rand7 src else return (mod n 7 + 1)

Since I didn’t like the chance of a non-terminating algorithm I came up with the following algorithm: sum seven random numbers in the range 1 to 5, and give the result modulo 7. Or in code:

rand7approx :: (Functor f, Monad f, Integral a) => f a -> f a
rand7approx src = succ . (`mod` 7) . sum <$> replicateM 7 src

Printing the distribution looked perfectly fine, with all of them getting a little over 14%. But something felt off to me. The code was shorter than the version with rerolling and didn’t suffer the inifinite loop problem, so there had to be reason it wasn’t the provided solution. I changed all the functions to take a source argument rather than calling rand3 and rand5 directly and printed the theoretical distributions as well (except for rand7, since the distribution for that one is non-terminating). This revealed that there is a 0.3% difference between the most and least common values. In retrospect, this should’ve been obvious, as 5^7 mod 7 isn’t equal to 0, which means it’s impossible for all 7 values to get the same distribution. Some more reading up revealed that there is indeed no solution that both provides an equal distribution and is guaranteed to terminate, since 7 is an infinitely repeating fraction in base 5.

Here are the tests I used:

main :: IO ()
main = do print $ dist (rand9 [1..3]) == zip [1..9] [1,1..]
          print . dist =<< replicateM 100000 (rand9 rand3)
          print . dist =<< replicateM 100000 (rand7 rand5)
          print $ dist (rand7approx [1..5])
          print . dist =<< replicateM 100000 (rand7approx rand5)

Programming Praxis – Rule 30 RNG

April 29, 2011

In today’s Programming Praxis exercise, our goal is to write a random number generator based on the Rule 30 cellular automaton. Let’s get started, shall we?

Some imports:

import Control.Monad.State
import Data.Bits
import Data.Word

First we define the amount of bits in our RNG state.

size :: Int
size = 7 * 43

We’re going to be needing a function to convert a list of bits to a number.

fromBits :: (Bits a, Integral a) => [Bool] -> a
fromBits = foldl (\a x -> shift a 1 .|. fromIntegral (fromEnum x)) 0

The step function replaces the state with the next generation of cells for a given rule and returns the middle bit.

step :: Monad m => Int -> StateT Integer m Bool
step r = modify (\s -> let b i = testBit s (mod i size) in
    fromBits [ testBit r $ fromBits [b (n+1), b n, b (n - 1)]
             | n <- [size - 1, size - 2..0]]) >> gets (`testBit` div size 2)

To produce a random number, we take the middle bit of the desired number of generations and convert that to an unsigned integer.

randomBits :: Monad m => Int -> Int -> StateT Integer m Word
randomBits n r = return . fromBits =<< replicateM n (step r)

To test the algorithm, we see if we get the correct result when we only set the bit in the middle. Afterwards, we produce a bunch of 8-bit numbers.

main :: IO ()
main = evalStateT (do r <- randomBits 32 30
                      liftIO $ print (r == 3112904540)
                      liftIO . print =<< replicateM 100 (randomBits 8 30)
                  ) $ bit (div size 2)

Looks like everything is working properly.