Norsetto's Journal

In the last post I mentioned how we could use the matMult function to compute the powers of a given matrix.
I think it is worth expanding some of the concepts I briefly touched upon: the first thing we could consider, is that the power of two is simply the multiplication of our matrix by itself:

let a2 = matMult a a

We can then continue and define different powers, all to an exponent which is itself an exponent of 2:

let a4 = matMult a2 a2
let a8 = matMult a4 a4
.
.
.
let a256 = matMult a128 a128

We can do this by simply defining a matPow2 function:

matPow2 :: Array (Int, Int) Double -> Int -> Array (Int,Int) Double
matPow2 a n
  | n == 0 = a
  | otherwise = matMult (matPow2 a (n-1)) (matPow2 a (n-1))

I think that the interest of doing this should be evident: I’m trying to find a way to reduce the order of the power function, from O(n) (the n^th power of a matrix requires n multiplications) to O(log₂n) (the n^th power of a matrix requires in the order of log₂n multiplications).
The function we just defined is apparently working but only if n is exactly a power of 2. What if I want to have my matrix to a power which is not itself a power of 2?
Could I express my exponent as a sum of powers of 2? Say e = a₁*2ⁿ+a₂*2^(n-1)+…+a_(n-1)*2+a_n, does it ring somehow?
It really seems we can: its exactly as a change of base, we are expressing our exponent in base 2.

For instance, suppose we want to compute a¹⁵. 15 can be expressed as 8+4+2+1 (say 0xF if you wish), so we can simply compute the power as:

a¹⁵=a²³*a^22*a2*a

Or in code:

matMult (matMult (matPow2 a 3) (matPow2 a 2)) (matMult (matPow2 a 1) a)

So, if we have a list with all the required exponents of 2 (say [3,2,0]) we can compute our power simply as:

foldr1 matMult (map (matPow2 a) [3,2,0])

The map will expand to:

foldr1 matMult [matPow2 a 3, matPow2 a 2, matPow2 a 0]

And the foldr will expand to:

(matMult (matPow2 a 3) (matMult (matPow2 a 2) (matPow2 a 0)))

Now we just have to find a way to decompose our exponent into a list of powers of 2 (or, if you wish, find out its binary representation).
Haskell has a nice builtin function we can use, properFraction. It returns a pair whose first element is the integer part of the argument, and the second is its fractional part.
Hmmm, if the argument is log₂x we can check immediately if it is a power of 2, since its fractional part will be zero. If its not zero, the integer part will give us the exponent of the most significant power of 2 and then we can use the fractional part again to compute the remaining bits (pun intended). In code:

factorise :: Int -> [Int]
factorise n
  | b == 0 = [a]
  | otherwise = a:(factorise (n-2^a))
      where (a,b) = properFraction $ log (fromIntegral n) / log(2)

We can spot check if this works. For instance:

*Main> sum $ map (2^) (factorise 196)
196

So, putting it all together:

matPow :: Array (Int,Int) Double -> Int -> Array (Int,Int) Double
matPow a n = foldr1 matMult (map (matPow2 a) (factorise n))

Now, these functions do indeed do their job but I’m not very happy by their efficiencies. It just takes forever to simply compute anything with an exponent greater than a mere 4. Lets expand our matPow2 function and lets check what is going on:

matPow2 a 0
a
 
matPow2 a 1
matMult (matPow2 a 0) (matPow2 a 0)
matMult a a
 
matPow2 a 2
matMult (matPow2 a 1) (matPow2 a 1)
matMult (matMult (matPow2 a 0) (matPow2 a 0)) (matMult (matPow2 a 0) (matPow2 a 0))
matMult (matMult a a) (matMult a a)
 
matPow2 a 3
matMult (matPow2 a 2) (matPow2 a 2)
matMult (matMult (matPow2 a 1) (matPow2 a 1)) (matMult (matPow2 a 1) (matPow2 a 1))
matMult (matMult (matMult (matPow2 a 0)  (matPow2 a 0)) (matMult (matPow2 a 0) (matPow2 a 0))) (matMult (matMult (matPow2 a 0) (matPow2 a 0)) (matMult (matPow2 a 0) (matPow2 a 0)))
matMult (matMult (matMult a a) (matMult a a)) (matMult (matMult a a) (matMult a a))

Its not surprising that this function is inefficient, look at all the repeating computations it has to make, all of which it has somehow already computed! Its pretty evident that we have not achieved our goal at all, the order is still O(n) …
When we do it manually, we compute the intermediate results and then combine them. This just gives us a limited number of multiplications. For instance, by computing first a2 = matMult a a and then a4 = matMult a2 a2 we only do two multiplications, while my function does 7, and thats just with an exponent of 2³. The number of multiplications is indeed exponential, 0 for 0, 1 for 1, 3 for 2, 7 for 3 or in general 2^n-1.
In effect, if we could find a way to automatically memorise the intermediate results of our matPow2 function we could reduce its order from O(2ⁿ) to O(n); that would be quite a gain (and this will really bring the order of the matPow function to a nice O(log₂n)).

You may have heard about it, and if not I encourage you to read about it, but the answer to this problem is memoization.

To memoize our power function we just need a way to easily map from its input arguments to an internal cache that contains already computed powers. The logic would be:

1. Map the input arguments to the cache
2. If there is no corresponding entry in the cache then:
2.1. Evaluate the function
2.2. Insert into the cache the computed value and return it
3. If the cache returns a value, just return it

I don’t think you will find it surprising that we can implement this with:

Data.Map for the mapping function
Control.Monad.State for the cache (which ought to remain mappable for all computations)

First of all we should define some aliases to keep our code readable:

type TM = Array (Int,Int) Double
type PowMap = Map.Map Int TM
type PowMapS = State PowMap TM

TM is our transition matrix. PowMap is the dictionary, indexed by the power exponent and containing the results of our computations (matrices). Note that we do not use the matrix to be exponentiated as part of the key even though strictly speaking it is one of its arguments. Even though we explicitly use it as an argument, in our application the
matrix to be exponentiated is more part of the state, it is indeed invariant in all calls.
PowMapS is the State implementation of the dictionary.

Our memoized function just implements in code the above algorithm.
We first map the input arguments and see if the dictionary returns something:

memoizedPow2 :: TM -> Int -> PowMapS
memoizedPow2 a n = do
             powmap <- get -- We fetch the dictionary
             maybe (matPow2' a n >>= update a n) -- combines Step 2.1 and 2.2
                    return
                    (Map.lookup n powmap) -- if it returns nothing we go to Step 2.1

You may want to check the syntax of the maybe function.
Step 2.1 is to compute the function:

matPow2' :: TM -> Int -> PowMapS
matPow2' a n
  | n == 0 = return a
  | otherwise = do
                b <- matPow2' a (n-1)
                return (matMult b b)

As you can see very similar to the unmemoized version, just adapted for the State Monad.
And then Step 2.2, insert it into the dictionary and return it:

update :: TM -> Int -> TM-> PowMapS
update a n ans = do
  powmap <- get
  put (Map.insert n ans powmap)
  return ans

Lets check if this works:

let a2 = matMult a a
let a4 = matMult a2 a2
let a8 = matMult a4 a4
let a16 = matMult a8 a8
let a32 = matMult a16 a16
let a64 = matMult a32 a32
let a128 = matMult a64 a64
let b = evalState (memoizedPow2 a 7) Map.empty
*Main Control.Monad.State> b == a128
True

Cool! It does its job.

We can now wrap it up and are ready to use it:

matPow2 :: TM -> Int -> TM
matPow2 a n = evalState (memoizedPow2 a n) Map.empty

And what about the time it takes? Well, feel free to time it on your machine, on mine it takes 1m30s to compute the unmemoized matPow2 with an exponent of 7 and just a couple of seconds to compute the memoized version, yuppie!
The difference will become even more evident with higher exponents (remember the order of these functions…). (matPow a 256)!(0,63) its really really fast now.
You may also want to monitor your memory usage (otherwise, why do you have that conky on your desktop?); we are trading speed vs. memory here, so it is not surprising that memory usage gets higher. Remember, we are memorizing the whole matrices as values of our dictionary.

Tags: haskell, programming, ubuntu

We have finally our transition matrix, found out that the game always succeeds, and now we should try to find how long it will take. We can’t obviously give a precise answer (512 steps or whatever) but we can estimate what is the probability that it will be within a certain amount of steps.

If you read the wikipedia article, you should know that being A a transition matrix, A² will also be a transition matrix. Its (i,j) element will give the probabiliy to go from state i to state j in two steps. In general, A^n (A to the power of n) is a transition matrix that will give the probabilities of going from a state to another in n steps.

We can therefore compute the probabilities that it will take to go from any given state to the end state (square 63) in n steps simply by taking the last column of A^n. The first element of this array in particular will give the probability of going from the start (square 0) to the end in n steps.

Lets take our A matrix, this element can be simply computed as:

countFreq 63 0 10000

(10000 being arbitrarily chosen, pick any number that won’t cause your computer to bend to its knees), which gives 0.1099 (very close to the theoretical probability of 0.111). It is clear that we are not very near to the end of the game yet.

To go further we need to be able to compute the nth power of our matrix A. Instead of reinventing the wheel we will just use the Array class, available from Data.Array. We need to modify our buildTM function so that the matrix will be constructed in accordance with the class definition:

buildTM' :: Double -> Array (Int, Int) Double
buildTM' n = array ((0,0),(63,63)) [((st,sq), countFreq sq st n) | st <- [0..63], sq <- [0..63]]

We are also not going to reinvent the matrix multiplication function, there is already a nice one given here:

matMult :: Array (Int,Int) Double -> Array (Int,Int) Double
	   -> Array (Int,Int) Double
matMult x y =  array ((0,0),(63,63))
                     [((i,j), sum [x!(i,k) * y!(k,j) |k <- range (0,63)])
                                       | i <- range (0,63),
                                         j <- range (0,63)]

Now we can compute A² simply by doing this in ghci:

let a = buildTM' 10000 -- (again feel free to experiment, this takes about 15 minutes on my PC)
let a2 = matMult a a

And the last column of it will be [a2!(i,63)|i<-[0..63]]. The first element of this column will be the probability of going from the first square to the last square in 2 steps:

a2!(0,63) -- 0.1099

Well, same result as before. Were would we expecting anything different!?
Now, increasing the exponent should start giving some results.

let a4 = matMult a2 a2
let a8 = matMult a4 a4
let a16 = matMult a8 a8
let a32 = matMult a16 a16

And now:

a32!(0,63) -- 0.7189049811868098

So the probability of ending the game in less than 32 steps is already quite high, about 72%.

let a64   = matMult a32 a32
let a128 = matMult a64 a64
let a256 = matMult a128 a128
 
a64!(0,63)  -- 0.9292194901730458
a128!(0,63) -- 0.9955122418509581
a256!(0,63) -- 0.9999819589853614

Remember that previously we found that the longest sequence in 10000 simulations was 219? As we see, the probability of ever finding a sequence longer than 256 steps is pretty low (I don’t know what you think, but a probability of 0.0018% is pretty near to “impossible” in my vocabulary). Just for curiosity, what is the probability of finding a sequence longer than 219 steps? We can write 219 as 128+64+16+8+2+1 so using our already computed matrices we easily find that:

a219!(0,63) -- 0.9999111233274545

which is indeed very close to 9999/10000 …

Do you remember the clumsy countFreq function I wrote? It turns out that the Array class provides a built-in function to do exactly this (accumArray). We can therefore cut the middle man and build our transition matrix directly using it:

countFreq :: Int -> Int -> Int -> Array Int Int
countFreq st n s = accumArray (+) 0 (0,63) [(i,1)| i <- runSim (n-1) st s]

buildTM :: Int -> Int -> Array (Int,Int) Double
buildTM n s = array ((0,0),(63,63)) [((i,j),
	fromIntegral ((countFreq i n s)!j) / fromIntegral n) |
	i <- range (0,63),
	j <- range (0,63)]

I don’t know how smart ghc is, but I expect that the countFreq in (countFreq i n s)!j) will be evaluated only once (somehow I think there should be a simpler method to combine arrays into a matrix). Both methods seem to be equivalent in terms of timing.

Assuming that 4 players are playing and that each take 15 sec per turn, it seems pretty likely that few games (if any) will last longer than half hour (since we have 4 kids, the single probabilities have to be combined). Since the probability of having 4 kids playing together in a room for half an hour without a fight breaking is pretty low, I think we can safely enunciate the Goose theorem:

It is very probable that 4 kids playing the Game of the Goose will have a fight

And this i can definitively prove from personal experience.

Tags: haskell, programming, ubuntu

I left the previous post with a couple of interesting questions. We can answer them using Markov chains (I’m not going to explain what they are, the linked wikipedia article is more than clear on the subject).

It is clear that the Game of the Goose is a Markov process: the state at any given step is only a function of the last state. What is interesting is that we can use the code of the last post to estimate the transition matrix.
The element i,j of the matrix is the probability of going from the state i to the state j. But the state in this example is just our position on the board. We can simulate, for each square, where we end up to. If we do that for a reasonable amount of time we end up having a frequency distribution, from which we can easily compute the probabilities simply dividing by the number of tries.

To do this, we need to sligthly modify our turn function, so that it takes as an input the current square we are investigating, and doesn’t loop but yield the resulting (simulated) destination square:

turn :: Int -> State ([Int], StdGen) ()
turn cur = do
     (xs,gen) <- get
     (d1,gen) <- return $ randomR (1,6) gen
     (d2,gen) <- return $ randomR (1,6) gen
     int <- inPrison cur (d1+d2)
     fin <- goTo int (d1+d2)
     put(fin:xs,gen)

We also need to modify the runSim function, so that it returns a list with all the destination squares:

runSim :: Int -> Int -> Int -> [Int]
runSim n sq s = fst $ execState (foldr (>>) (turn sq) (replicate n (turn sq))) ([], mkStdGen s)

For instance running it with runSim 100 0 0 we get as a result:

[12,7,12,7,8,10,4,63,8,7,12,7,10,2,8,10,8,8,8,7,8,7,7,7,10,7,63,63,4,7,8,7,4,63,8,10,3,63,3,7,
10,10,8,12,10,10,8,10,10,63,8,10,3,63,12,12,6,4,4,63,8,63,11,10,10,10,8,63,8,8,8,12,10,3,2,
8,63,4,12,10,63,8,8,10,63,63,11,7,11,7,12,10,7,10,12,11,3,12,7,10,6]

To excerpt the frequency we can use the following clumsy function (I guess it could be done more efficiently, but right now I can’t be bothered to think about it …):

countFreq :: Int -> Int -> Double -> Double
countFreq sq st n = fromInteger ( sum $ map (\x -> if x == sq then 1 else 0) (runSim (round (n-1)) st 0) ) / n

Now we can check if we get some meaningful results. For instance:

Main> [countFreq n 42 10000 | n <- [0..63]]
[0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.719,0.0,0.0,0.0,0.0,0.1138,0.0,0.1672,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0]

You can see this already as a single row of the transition matrix, the 42th row. The first number is the probability of going to the start from square 42 (0), the second one is the probability of going to the first square from square 42 (again 0) and so on.
We can check the correctness of this result since we can compute the probabilities of exiting from the prison (square 42) pretty easily.
We can exit throwing a 5 or a 7. We can get a 5 with 4 possible dice combination, out of 36 ((1,4),(2,3),(3,2),(4,1)), so the probability is 4/36=0.111. We can get a 7 with 6 possible dice combination, so the probability is 6/36=0.167. We are stuck in the prison in all other cases (26/36=0.722). Indeed, the probabilities computed with haskell match these pretty well.
We can spot check other rows. For row 0:

*Main> [countFreq n 0 10000 | n <- [0..63]]
[0.0,0.0,2.97e-2,5.42e-2,8.21e-2,0.0,2.5e-2,0.1672,0.1346,0.0,0.2007,5.53e-2,0.1383,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.1129]

Again we can verify it easily. The probability of going from the start to 63 is by throwing a 9, that can be done with 4 possible dice combination -> 0.111. You can check yourself that the results are indeed reasonable, just mind the special squares. For instance the probability of going from the start to square 6 is 1/36=2.78e-2 (we can only go to this square passing through the well at square 12, which is only reachable with a double 6 from the start).

If you check other rows you will see a couple of interesting results. First for square 63, I’d expect to see all zeroes and a final 1. I don’t get this, and thats because inPrison needs to be modified to not advance if our starting square is indeed 63 (just add a | square == 63 = return 63 as a first guard in there).

I’m also getting an endless loop for row 51 and row 59. Ah ah! This must be the source
for the strange results we were getting in the previous post. Indeed, throwing an 8 at square 51 leads to square 59, but square 59 requires you to advance by 8, which brings you back to 59 and so on. The same if you throw an 8 directly in square 59.

Now, I don’t have a cardboard version available and frankly I don’t remember, but by checking some pictures on the net confirms that square 59 is a Goose square. So, either the rules I use are wrong, or really this is a specific design of the original game.
Anyhow, I don’t like it, so I got this square out from my board and rerun the findLongest computation:

*Main> findLongest 10000
[63,60,60,60,60,60,61,57,62,60,61,56,61,62,57,49,42,42,42,42,42,28,8,0,48,43,39,46,26,19,16,10,2,0,59,57,61,60,60,62,57,
48,34,20,10,6,6,7,0,62,57,48,39,47,42,42,42,42,22,10,0,55,61,55,35,29,22,6,7,0,55,61,47,42,42,42,42,37,22,10,0,62,60,57,
62,57,53,62,59,56,62,55,62,56,61,59,53,61,61,62,55,62,57,61,53,37,17,15,7,0,61,55,61,61,59,61,53,62,62,62,56,62,59,49,
43,40,37,33,21,10,0,62,55,62,59,59,49,42,42,24,15,10,0,56,49,33,26,19,10,0,62,56,62,56,61,61,56,61,57,59,55,61,57,61,57,
60,55,61,61,61,55,60,61,62,59,57,46,43,29,22,15,3,0,61,57,55,61,57,51,43,37,26,2,0,53,47,35,19,12,0,61,59,53,62,55,60,59,
59,60,62,55,47,38,35,26,22,13,8,0]
*Main> length it
219

Oh well, a partial answer already Having removed the Goose from square 59 gives us a game that is likely to always converge.

Just for curiosity, take a look at rows 51 to 62, and see what are the probabilities of going back to the start or winning from those. Amazingly (or perhaps not), there are a number of squares where you are more likely to have to start again than win. The highest probability of winning is a mere 16% from square 56.

Finally, we can compute our transition matrix:

buildTM :: Double -> [Double]
buildTM n = [countFreq sq st n | st <- [0..63], sq <- [0..63]]

We can now fully answer the first question:

Are there really periodic solutions which will never converge, and how probable are they?

Having got rid of square 59, we can easily verify that the array p = [0,0,0 ...., 0,1] (which is row 63 of the TM, and an absorbing state) satisfies the stationary distribution constraint:

Sum (i=0 ..63) p(i)TM(i,j) = p(j)

so the answer to the question is that the game always succeed.
To answer the second question we need a way to treat the TM matrix in haskell, which we will see in part 3.

Tags: haskell, programming, ubuntu

Now that I know how to add random numbers to a state-full sequence I can start doing some funky stuff.
Meet the “Game of the Goose”, a very simple board game suitable for young kids (of any age ); if I remember it correctly the first board game I ever played.
The rules are very simple, there is a board with 64 squares, numbered 0 to 63. All players start at square 0 (also surprisingly called the start square). At each turn a player throws two dice and advances by the given amount. Wins the first player to reach EXACTLY square number 63. If a player advances further than this, it will then proceed backwards by the remaining amount and continue playing.
There are a number of special squares too, like 9 and 14 (the goose squares, and every 9 after those), which allow the player to advance by the same amount once again, or square 58 (the death) which brings the player back to the start square, the prison, the wells, the bridge, the labyrinth and so on (details in the code, I simulated all of them but square 19, which makes the player loose 3 turns, not very relevant for a simulation and pretty annoying since it requires to augment the state).

To simulate the game in haskell we need first of all a function that is computed every turn. I wrote the following:

turn :: State ([Int], StdGen) ()
turn = do
     (cur:xs,gen) <- get
     (d1,gen) <- return $ randomR (1,6) gen
     (d2,gen) <- return $ randomR (1,6) gen
     int <- inPrison cur (d1+d2)
     fin <- goTo int (d1+d2)
     put(fin:cur:xs,gen)
     if fin /= 63 then turn else return ()

For each turn, we throw the dice, we check if we are in prison and if we are allowed to exit (forget for the time being that you can also be freed if another player ends up in it).
We then go to the new square and we do whatever action the new square obliges us to do and loop back until we are on square 63.
The two functions this uses are very simple to write:

inPrison :: (Num a, Monad m) => a -> a -> m a
inPrison square dice
      | square == 42 = if (dice == 5) || (dice == 7) then
      	       return (42+dice) else return 42
      | otherwise = return (square+dice)

goTo :: (Num a, Ord a, Monad m) => a -> a -> m a
goTo square dice
     | square == 63 = return 63
     | square > 63 = goTo (126-square) dice
     | square == 5  || square == 9  || square == 14 || square == 18 ||
       square == 23 || square == 27 || square == 32 || square == 36 ||
       square == 41 || square == 45 || square == 50 || square == 54 ||
       square == 59 = goTo (square+dice) dice
     | square == 6  = return 12
     | square == 12 = return 6
     | square == 31 = return 42
     | square == 52 = return 39
     | square == 58 = return 0
     | otherwise = return square

Note that I do not take a record of each visited square, only of the starting and ending one. So, for instance, if my first throw is a 9 I win immediately and the result of the simulation will be [63,0].

Now, we only need some logic to start the simulation:

runSim :: Int -> [Int]
runSim s = fst $ execState turn ([0], mkStdGen s)

I don’t know how the pseudo-random number generator works for different haskell implementation (and different arches). On my machine (The Glorious Glasgow Haskell Compilation System, version 6.12.1, i686) runSim 0 gives the following result:

[63,55,62,53,48,43,33,25,16,8,0,60,60,61,62,62,56,61,61,61,56,61,57,61,55,62,53,62,55,60,60,60,61,47,42,42,42,42,25,11,6,0]

This is one of my typical games, end up in prison and be there for few turns, approach quickly the end and fail it several times, being thrown back to square 0 and finally getting to the end square. In the mean time my opponents are likely to play like:

[63,56,60,61,61,60,55,61,56,62,57,47,43,21,17,10,0]

or even the infamous:

[63,0]

Needless to say, this gave us hours of endless fun, punctuated by accusations of cheating, heavy cheating, stealing and finally a big fight. Ah, those were the great days ….

Anyhow, this just wet my appetite. What could it be the longest possible game?
What I really like about haskell is that its pretty easy to add computations to some early result, and greatly expand the capability of your code with few additional lines.

In this case, I just add the following:

longest :: [a] -> [a] -> [a]
longest a b = if length a > length b then a else b

findLongest :: Int -> [Int]
findLongest n = foldl1 longest [runSim s | s <- [0..n]]

and I’m set.

Running it with findLongest 100 however, gives me the first surprise. The code apparently loops forever. Is this a bug in my code? A bug (feature!?) of the game design? A problem with the pseudo-random number generator?
Checking it in details I uncover that there is a number of seeds that cause this, surprisingly I find at least 9 in the first 100 seeds I try. Fishing those out like this:

findLongest n = foldl1 longest [runSim s | s <- [0..n], s `notElem` [22,30,31,33,58,75,79,88,94]]

and findLongest 100 gives me a list of 127 squares.

Being the curious twit I am, I immediately pose myself few questions:

1. Are there really periodic solutions which will never converge, and how probable are they?
2. What is the probability of finding sequences longer than 127 squares?
3. Is yahoo/bing being given prime spot in Ubuntu the final drop that will make me stop using it?

We can answer the first two of these questions with Markov chains and some more bits of haskell in Part 2.

Tags: haskell, programming, ubuntu

It took me a full two hours to come up with this, and most probably the smarter way to do it was using Monad transformers and Control.Monad.Random (or something as contrived as that). But I still have to understand monads so I’m quite a way off to understand transformers and the significance of Control.Monad.Random really eludes me.
Anyhow, the problem is how to add random numbers to a state-full computation.
Anyone which played with Haskell (or I guess any other functional language) should know that random numbers are kind of special, as you are practically forced to carry the generator around.
Now, suppose you are doing something using the state monad, and you need to use random
numbers in there (Markov knocks at the door).
A very trivial example just to get to the point could be the following code:

import Control.Monad.State

tick :: State Int ()
tick = do
     a <- get
     put(a+1)

Now, executing it in ghci with:

execState tick 0

gives 1 (doh!).
Now, suppose that instead of just adding 1 we want to add a random (integer) number in the range 1 to 6. To be able to use randomR we need to somehow carry the random number generator with us.
The trivial choice would be to augment the state. Instead of having a state of one number, we have a state with a pair (number, generator).
For instance:

import Control.Monad.State
import Random

tick:: State (Int, StdGen) ()
tick = do
     (a,gen) <- get
     (r,gen) <- return $ randomR (1,6) gen
     put(a+r,gen)

We can test it using the following simple function:

addRand :: Int -> Int -> Int
addRand m n = fst $ execState (sequence $ replicate m tick) (n, mkStdGen 0)

And testing it with:

addRand 10 0

gives 33 (pretty near to the expected average of 35).
Now, if only some smarty-pants could point me to a kindergarten level tutorial about monad transformers I would really be obliged.

Tags: haskell, programming, ubuntu

Still on the previous example, I wanted to see if I could find a monad to solve the problem. I thought that the state one
could be usefull, and indeed this is the code I came up with (using Control.Monad.State, which is in libghc6-mtl-dev if you are on Debian or derivatives):

import System.Environment
import Control.Monad.State

taylor :: Double -> Double -> Double
taylor x tol = reduce1 (x - fromIntegral(truncate(x / (2*pi)))*(2*pi)) tol

reduce1 :: Double-> Double -> Double
reduce1 x tol
   | x < 0 = (-reduce2 (-x) tol)
   | otherwise = reduce2 x tol

reduce2 :: Double -> Double -> Double
reduce2 x tol
   | x <= pi/2 = sumTaylor x tol
   | x <= pi = sumTaylor (pi - x) tol
   | x <= 3*pi/2 = (-sumTaylor (x - pi) tol)
   | otherwise = (-sumTaylor (2*pi - x) tol)

sumTaylor :: Double->Double->Double
sumTaylor x tol = taylorInternal 1.0 1.0 1.0 1 x tol

taylorInternal :: Double->Double->Double->Int->Double->Double->Double
taylorInternal num den sign index x tol =
  execState (loop num den sign index) 0
    where
      loop num den sign index
        | index > 100 = return ()
        | abs(num/den) < tol = return ()
        | even index = do
                        modify (sign*num/den +)
                        loop (num * x)
                             (den * fromIntegral index)
                             (-sign)
                             (index + 1)
        | otherwise = loop (num * x)
                           (den * fromIntegral index)
                           (sign)
                           (index + 1)

main :: IO ()
main = do
    (x:tol:[]) <- getArgs
    print $ taylor (read x) (read tol)

It looks indeed more “iterative” than the previous version using scanl. Performance is about the same with the executable size sligthly larger. One of the advantages is that one can check if the result is within the given tolerance before computing the new state, which spares one iteration cycle.

Now, to really enter into monadoom I should really try to define an ad-hoc monad to solve this particular problem.

Tags: haskell, programming, ubuntu

I’ve been dabbling with haskell in the recent days. My first semi-serious attempt to some code was inspired by a prompt at programming praxis.
What I concocted should ashame any veteran haskell programmer (if you are and have stumbled on this by accident feel free to point me in the right direction;-)):

import System.Environment

taylor :: Double -> Double -> Double
taylor x tol = reduce1 (x - fromIntegral(truncate(x / (2*pi)))*(2*pi)) tol

reduce1 :: Double-> Double -> Double
reduce1 x tol
   | x < 0 = (-reduce2 (-x) tol)
   | otherwise = reduce2 x tol

reduce2 :: Double -> Double -> Double
reduce2 x tol
   | x <= pi/2 = sumTaylor x tol
   | x <= pi = sumTaylor (pi - x) tol
   | x <= 3*pi/2 = (-sumTaylor (x - pi) tol)
   | otherwise = (-sumTaylor (2*pi - x) tol)

sumTaylor :: Double -> Double -> Double
sumTaylor x tol = checkdelta (scanl sumOdd (0,1,1,-1) [1..]) tol
                  where sumOdd (p,n,d,s) el
                          | even el = (p-s*n/d,n*x,d*fromIntegral el,-s)
                          | otherwise = (p,n*x,d*fromIntegral el,s)

checkdelta :: [(Double,Double,Double,Double)] -> Double -> Double
checkdelta (a:b:c:xs) tol = if abs ( first4 a - first4 c ) < tol then first4 c else checkdelta (c:xs) tol
checkdelta [] _ = 0

first4 :: (Double,Double,Double,Double) -> Double
first4 (a,b,c,d) = a

main :: IO ()
main = do
    (x:tol:[]) <- getArgs
    print $ taylor (read x) (read tol)

I was curious to see how this would fare against a C equivalent, so I quickly typed the following:

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define max_iter 100

/*
 * Compute the taylor sum x-x**3/3!+x**5/5!+ ...
 * Stop if the difference between 2 consecutive terms is minor than a given
 * tolerance or if it doesn't converge after 100 iterations.
 *
 */

double sum(const double x, const double tol) {

  double y, old_y, num, den, sign;
  int i;

  y = 0.0;
  num = 1.0;
  den = 1.0;
  sign = 1.0;
  old_y = 100.0;

  for ( i=1; i<max_iter; i++ ) {
    if (!(i & 1)) {
	y += sign*num/den;
	sign = -sign;
	if(fabs(y-old_y)<tol) return y;
	old_y = y;
    }
    num *= x;
    den *= (double)i;
  }
  return NAN;
}

int main(int argc, char **argv) {

  double x, tol, pi, y;
  int i, sign = 0;

  if (argc<2) {
    printf("Requires at least a numeric argument\n");
    return -1;
  }

  pi = atan(1.0)*4.0;
  x = strtod(argv[1], NULL);

  if (argc==3) {
     tol = strtod(argv[2], NULL);
  } else {
     tol = 1.0e-9;
  }

  //Only compute for positive x (since sin(x)=-sin(-x))
  if ( x < 0 ) {
     x = -x;
     sign = 1;
  }

  //Limit x to the 0-2pi range (sine is periodic with period=2pi
  x -= floor(x/(2*pi)) * (2*pi);

  //Only compute in the first quadrant and use symmetries
  if ( x <= pi/2 ) {
     y = sum(x, tol);
  } else if ( x <= pi) {
     y = sum(pi-x, tol);
  } else if ( x <= 3*pi/2 ) {
     y = -sum(x-pi, tol);
  } else {
     y = -sum(2*pi-x, tol);
  }

  if (sign) y = -y;

  printf("%.16f\n", y);

  return 0;
}

Code size wise the haskell one is indeed more compact (and I think there is still room for improvement), however the size of the executable is amazingly different. With default compiler settings for gcc and gch the C program is 8.2 kB and the haskell one 1007 KB!? I wonder what is all that overhead …

Also, performance wise the C code is about 100% faster than the haskell one (but that I expected).

Well, its not really a fair comparison, the problem is indeed more tailored to imperative programming.

Tags: haskell, programming, ubuntu