fragmentary ideas ䷿ intellectual what-nots ䷷ and haskell programming ䷴

The Algorithm

We can use Brent’s Algorithm to detect infinitely-repeating sequences in a list of values generated by some iterated function: that is, any list in which the next value in the sequence is generated from the previous value alone; if we find duplicate values in the list, we know we have a cycle.

The implementation below isn’t particularly elegant, and since I want to use it as a stand-alone tool I’m having it output strings:

module Main 
    where

cycling :: (Show a, Eq a) => Int -> [a] -> String
cycling k []     = "Empty list"
cycling k (a:as) = find 0 a 1 2 as
    where find _ _ c _  [] = "reached end at " ++ show c ++": no cycles"
          find i x c p (x':xs) 
            | c >  k    =  "no cycles after " ++ show k
            | x == x'   =  "cycle at "++ show c ++": "++(show$take (c-i) xs)
            | c == p    = find c x' (c+1) (p*2) xs 
            | otherwise = find i x  (c+1) p xs

  --- SOME RANDOM NUMBER GENERATORS TO TEST ---
-- bad generators?:
g1 = 0 :      [ (g*7 + 1)   `mod`   32       | g <- g1]
g2 = 17 :     [ (g*22 + 221)   `mod`  2^32   | g <- g2]
g3 = 3249 :   [ (g*22695477 + 1)  `mod` 300  | g <- g3]
g4 = 234587 : [ (g*22695476 + 1) `mod` 2^32  | g <- g4]

-- good generators:
g5 = 294587 : [ (g*22695477 + 1) `mod` 2^32  | g <- g5]
g6 = 0 :      [ (g*22695477 + 1) `mod` 2^32  | g <- g6]

Randomness is hard…

Linear Congruential Generator’s are notoriously easy to screw up with the wrong parameters. Let’s use our function to test just how finicky this RNG algorithm can be. random stream g4 is almost identical to g5, a known-good generator. Let’s compare the two with a cut off limit of 999,999. First the good:

*Main> cycling 999999 g5
"no cycles after 999999"

Great, now let’s see how our typo-ed generator fares:

*Main> cycling 999999 g4
"cycle at 17: [180790021]"

… yikes.

GHC has an extension to the list comprehensions syntax that replicates the functionality of zipWith, by allowing you to have two generators running in parallel. Turn on the extension in GHCi like so:

Prelude> :set -XParallelListComp

I use the extension below to generate an infinite list of random coordinates (scaled to a  1000×1000 grid) using two different  Linear Congruential Generators running in parallel. It should be simple to modify the code below to actually  run the generators concurrently using Data Parallel Haskell (although I haven’t had a chance to play with that yet).

randPoints :: (Integral i) => i -> [(i,i)]
randPoints s = drop 1 [ (scale x, scale y) | x <- g1 | y <- g2 ]
    where -- two Linear Congruential Generators (from Borland C, glibc):
          g1 = s : [ (g*22695477 + 1)   `mod`   2^32   | g <- g1]
          g2 = s : [ (g*1103515245 + 12345) `mod` 2^32 | g <- g2]
          -- adjust the random output to fit a 1000x1000 plane:
          scale v = (v `mod` 1000) - 500

Then we can use another list comprehension, this time with boolean guards, to generate random points that lie within an outer circle (and outside an inner circle):

monteDonut = [p | p@(x,y) <- randPoints 42, let b= x^2+y^2, b< 500^2, b> 10000]

Finally from GHCi we can use the Gnuplot library to easily plot the points we generated:

*Main> :m + Graphics.Gnuplot.Simple
*Main> plotDots [Aspect(Ratio 1), XTicks Nothing, YTicks Nothing] (take 10000 monteDonut)


…and we get a nice Monte Carlo donut with sprinkles!: