Previously in this series:

In my previous post I introduced the 'knight moves problem': How many ways are there for a chess knight to reach cell (i,j) in exactly n moves? The recursive solution from last time is horribly inefficient for larger values of n. Today I will show some more efficient solutions.

If the knight first makes a move (-1,2) and then a move (2,1) it will end up at (1,3). If it first moves (2,1) and then (-1,2) it will also end up at (1,3). So, the order in which the moves happen does not matter for the final position! We can exploit this fact to make a faster program. Instead of determining what move to make at each step, we can count how many moves we make of each type and then determine in how many different orders these moves can be performed.

Denote by `n _{1}` the number of moves of the first type,

split n = [ (i,n-i) | i <- [0..n] ]

So for example, `split 3 = [(0,3),(1,2),(2,1),(3,0)]`.

By repeatedly splitting `n` we arrive at:

paths_{split}n (i,j) = sum $ do let n_{12345678}= n (n_{1},n_{2345678}) <- split n_{12345678}(n_{2},n_{345678}) <- split n_{2345678}(n_{3},n_{45678}) <- split n_{345678}(n_{4},n_{5678}) <- split n_{45678}(n_{5},n_{678}) <- split n_{5678}(n_{6},n_{78}) <- split n_{678}(n_{7},n_{8}) <- split n_{78}let counts = [n_{1},n_{2},n_{3},n_{4},n_{5},n_{6},n_{7},n_{8}] guard $ (i,j) == destination counts return $ multinomial counts

Here we only keep sequences of moves that end up in (i,j), as determined by the `destination` function:

destination counts = (sum hs, sum vs) where (hs,vs) = unzip [ (n*δ_{i},n*δ_{j}) | (n,(δ_{i},δ_{j})) <- zip counts moves ]

Next, we need to know how many different paths can be formed with a particular set of moves. You might remember binomial coefficients from high school, which give the number of ways to pick k items from a set of size n:

If we take n equal to m+k we get the number of different lists containing exactly k red balls and m green balls. Or put differently, the number of different paths containing k moves of the first type and m moves of the second type. This interpretation of binomial coefficients can be generalized two more than two types, giving *multinomial coefficients*. These are exactly what we need to determine the number of paths given the counts of each type of move:

multinomial xs | any (< 0) xs = 0 multinomial xs = factorial (sum xs) `div` product (map factorial xs)

This multinomial function requires calculating a lot of factorials, to make this as fast as possible they should be stored in an 'array':

factorial :: Int -> Integer factorial = unboundedArray $ scanl (*) 1 [1..]

Calculating `paths _{split}` only takes O(n

The above function uses a "generate and test" approach: Generate all possibilities and test which ones reach the destination. It would be more efficient to generate *only* those possibilities.

Algebraic reasoning can help us here. Let's start by expanding the condition in the `guard` statement:

(i,j) == destination counts = {- by definition of destination -} (i,j) == (sum hs, sum vs) where (hs,vs) = unzip [ (n*δ_{i},n*δ_{j}) | (n,(δ_{i},δ_{j})) <- zip counts moves ] = {- expand unzip and simplify -} i == sum (zipWith (*) counts (map fst moves) && j == sum (zipWith (*) counts (map snd moves) = {- by definition of moves (see previous post) -} i == sum (zipWith (*) counts [2,2,-2,-2,1,-1,1,-1] && j == sum (zipWith (*) counts [1,-1,1,-1,2,2,-2,-2] = {- expanding the sum and product, remember n_{12}= n_{1}+n_{2}, etc. -} i == 2*n_{12}- 2*n_{34}+ n_{57}- n_{68}&& j == 2*n_{56}- 2*n_{78}+ n_{13}- n_{24}= {- reordering -} n_{57}- n_{68}== i - 2*n_{12}+ 2*n_{34}&& n_{13}- n_{24}== j - 2*n_{56}+ 2*n_{78}

These are equations we can work with.
Take the equation involving `i`. We know that `n _{57} + n_{68} = n_{5678}`, and that

-- | find a and b such that a+b == c, a-b == d, a,b >= 0 solve_{pm}c d | ok == 0 && a >= 0 && a <= c = return (a,c-a) | otherwise = mzero where (a,ok) = (c + d) `divMod` 2

This gives an O(n^{5}) algorithm:

paths_{pm}n (i,j) = sum $ do let n_{12345678}= n (n_{1234},n_{5678}) <- split n_{12345678}(n_{12},n_{34}) <- split n_{1234}(n_{56},n_{78}) <- split n_{5678}(n_{57},n_{68}) <- solve_{pm}n_{5678}(i - 2*n_{12}+ 2*n_{34}) (n_{13},n_{24}) <- solve_{pm}n_{1234}(j - 2*n_{56}+ 2*n_{78}) (n_{1},n_{2}) <- split n_{12}let n_{3}= n_{13}- n_{1}let n_{4}= n_{24}- n_{2}(n_{5},n_{6}) <- split n_{56}let n_{7}= n_{57}- n_{5}let n_{8}= n_{68}- n_{6}return $ multinomial [n_{1},n_{2},n_{3},n_{4},n_{5},n_{6},n_{7},n_{8}]

It turns out that we don't actually need to know `n _{1}`,

But we could also first decide where the brightly colored balls (red and green) go and where the dark collored ones (blue) go. Now there are only two types of balls, so this is a binomial coefficient, or in terms of a multinomial, `multinomial [n _{rg},n_{b}]`. Then for the positions with brightly colored balls, we need to determine which ones are which color, which can be done in

This same arguments also holds when there are eight types of balls (or moves), so

multinomial [n_{1},n_{2},n_{3},n_{4},n_{5},n_{6},n_{7},n_{8}] == multinomial [n_{12},n_{34},n_{56},n_{78}] * multinomial [n_{1},n_{2}] * multinomial[n_{3},n_{4}] * multinomial [n_{5},n_{6}] * multinomial[n_{7},n_{8}]

If you plug this into the `paths _{pm}` function, you might notice that the last part of the function is calculating the product of two independent things. One part is about

paths_{O4}n (i,j) = sum $ do let n_{12345678}= n (n_{1234},n_{5678}) <- split n_{12345678}(n_{12},n_{34}) <- split n_{1234}(n_{56},n_{78}) <- split n_{5678}(n_{57},n_{68}) <- solve_{pm}n_{5678}(i - 2*n_{12}+ 2*n_{34}) (n_{13},n_{24}) <- solve_{pm}n_{1234}(j - 2*n_{56}+ 2*n_{78}) let result_{1234}= sum $ do (n_{1},n_{2}) <- split n_{12}let n_{3}= n_{13}- n_{1}let n_{4}= n_{24}- n_{2}return $ multinomial [n_{1},n_{2}] * multinomial[n_{3},n_{4}] let result_{5678}= sum $ do (n_{5},n_{6}) <- split n_{56}let n_{7}= n_{57}- n_{5}let n_{8}= n_{68}- n_{6}return $ multinomial [n_{5},n_{6}] * multinomial[n_{7},n_{8}] return $ multinomial [n_{12},n_{34},n_{56},n_{78}] * result_{1234}* result_{5678}

Here both of the `result` parts are of the form

sum [ multinomial [a,b] * multinomial[x-a,y-b] | (a,b) <- split n ]

which just so happens to be equivalent to just `multinomial [x,y]` (a proof of this statement is left as an exercise, i.e. I am too lazy to write it out).
This equation immediately leads to a (much simpler) O(n^{3}) algorithm:

paths_{O3}n (i,j) = sum $ do let n_{12345678}= n (n_{1234},n_{5678}) <- split n_{12345678}(n_{12},n_{34}) <- split n_{1234}(n_{56},n_{78}) <- split n_{5678}(n_{57},n_{68}) <- solve_{pm}n_{5678}(i - 2*n_{12}+ 2*n_{34}) (n_{13},n_{24}) <- solve_{pm}n_{1234}(j - 2*n_{56}+ 2*n_{78}) return $ multinomial [n_{12},n_{34},n_{56},n_{78}] * multinomial [n_{57},n_{68}] * multinomial [n_{13},n_{24}]

After all this manipulation it is a good idea to check whether the program still does the right thing. We can either manually compare the path matrices:

check paths = and [ pathMatrix paths_{rec}n == pathMatrix paths n | n <- [0..3] ]

Or use QuickCheck or SmallCheck:

Knight2> smallCheck 5 (\(N n) ij -> paths_{O3}n ij == paths_{rec}n ij) ... Depth 5: Completed 726 test(s) without failure.

Finally, to contrast with the first part of this series, here is the time it takes to calculate the number of paths in 100 steps:

Knight2> paths_{O3}100 (4,4) 2422219241769802380469882122062019059350760968380804461263234408581143863208781993964800 (4.75 secs, 270708940 bytes)

The recursive algorithm would need in the order of 10^{77} years to arrive at this answer.

Still, `paths _{O3}` is not the fastest possible algorithm.
Next time I will look at a completely different approach, but further improvements to the solution in this post are possible as well.
As an exercise for the reader, you should try transforming

## Comments

I'd be interested to see this done with sigfpe's combinatorial library:

http://sigfpe.blogspot.com/2007/11/small-combinatorial-library.html

Thanks for the link. For the third post in this series I have been working on a similar idea with a Num instance using convolution. It was developed independently, so the presentation will be a bit different.

Why the *2 after n12 and n56?

i == 2*n12*2 - 2*n34 + n57 - n68 && j == 2*n56*2 - 2*n78 + n13 - n24

Herbert: you are right, that is was error. Obviously there is only a single factor 2, from the 2 in

moves. I updated the article, it should be correct now.This is so related!

## Reply