Knight in n, part 2: combinatorics

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.

Ignoring the order of moves

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 n1 the number of moves of the first type, n123 the number of moves of type 1, 2 or 3, etc. So n1234 = n1+n2+n3+n4, and since there are eight different moves, n12345678 = n. A count nab can be split into na+nb in several ways, for now we will consider all possibilities:

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:

pathssplit n (i,j) = sum $ do
    let n12345678 = n
    (n1,n2345678) <- split n12345678
    (n2,n345678) <- split n2345678
    (n3,n45678) <- split n345678
    (n4,n5678) <- split n45678
    (n5,n678) <- split n5678
    (n6,n78) <- split n678
    (n7,n8) <- split n78
    let counts = [n1,n2,n3,n4,n5,n6,n7,n8]
    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: multinomial [2,1,1]

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 pathssplit only takes O(n7) integer operations, since each split effectively costs a factor n. While this is better than the previous result, it is still not satisfactory.

Solving the guard condition

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 n12 = n1+n2, etc. -}
   i == 2*n12 - 2*n34 + n57 - n68 &&
   j == 2*n56 - 2*n78 + n13 - n24
= {- reordering -}
   n57 - n68 == i - 2*n12 + 2*n34 &&
   n13 - n24 == j - 2*n56 + 2*n78

These are equations we can work with. Take the equation involving i. We know that n57 + n68 = n5678, and that n57 - n68 == i - 2*n12 + 2*n34. From these two equations, we can solve for n57 and n68, without needing an expensive split:

-- | find a and b such that a+b == c, a-b == d, a,b >= 0
solvepm 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(n5) algorithm:

pathspm n (i,j) = sum $ do
    let n12345678 = n
    (n1234,n5678) <- split n12345678
    (n12,n34) <- split n1234
    (n56,n78) <- split n5678
    (n57,n68) <- solvepm n5678 (i - 2*n12 + 2*n34)
    (n13,n24) <- solvepm n1234 (j - 2*n56 + 2*n78)
    (n1,n2) <- split n12
    let n3 = n13 - n1
    let n4 = n24 - n2
    (n5,n6) <- split n56
    let n7 = n57 - n5
    let n8 = n68 - n6
    return $ multinomial [n1,n2,n3,n4,n5,n6,n7,n8]

Multinomial laws

It turns out that we don't actually need to know n1, n2, etc. If you think about it, the multinomial coefficient [n1,n2,n3,n4,n5,n6,n7,n8] means: "The number of different lists with n1 red balls, n2 of green balls, etc.". To make such a list we can first pick where to put the red balls, then where to put the blue balls, then the green balls and so on.

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 [nrg,nb]. Then for the positions with brightly colored balls, we need to determine which ones are which color, which can be done in multinomial [nr,ng] ways. In a picture:

multinomial [2,1,1] = multinomial [3,1] * multinomial [2,1]

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

multinomial [n1,n2,n3,n4,n5,n6,n7,n8]
 == multinomial [n12,n34,n56,n78]
  * multinomial [n1,n2] * multinomial[n3,n4]
  * multinomial [n5,n6] * multinomial[n7,n8]

If you plug this into the pathspm function, you might notice that the last part of the function is calculating the product of two independent things. One part is about n1..n4 and the other about n5..n8. Now remember that the function paths takes the sum of all possibilities, and that products distributes over sums. This means that the two loops for n1234 and n5678 can be performed independently, giving us an O(n4) algorithm:

pathsO4 n (i,j) = sum $ do
    let n12345678 = n
    (n1234,n5678) <- split n12345678
    (n12,n34) <- split n1234
    (n56,n78) <- split n5678
    (n57,n68) <- solvepm n5678 (i - 2*n12 + 2*n34)
    (n13,n24) <- solvepm n1234 (j - 2*n56 + 2*n78)
    let result1234 = sum $ do
         (n1,n2) <- split n12
         let n3 = n13 - n1
         let n4 = n24 - n2
         return $ multinomial [n1,n2] * multinomial[n3,n4]
    let result5678 = sum $ do
         (n5,n6) <- split n56
         let n7 = n57 - n5
         let n8 = n68 - n6
         return $ multinomial [n5,n6] * multinomial[n7,n8]
    return $ multinomial [n12,n34,n56,n78] * result1234 * result5678

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(n3) algorithm:

pathsO3 n (i,j) = sum $ do
    let n12345678 = n
    (n1234,n5678) <- split n12345678
    (n12,n34) <- split n1234
    (n56,n78) <- split n5678
    (n57,n68) <- solvepm n5678 (i - 2*n12 + 2*n34)
    (n13,n24) <- solvepm n1234 (j - 2*n56 + 2*n78)
    return $ multinomial [n12,n34,n56,n78]
           * multinomial [n57,n68]
           * multinomial [n13,n24]

Verifying the results

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 pathsrec n == pathMatrix paths n | n <- [0..3] ]

Or use QuickCheck or SmallCheck:

Knight2> smallCheck 5 (\(N n) ij -> pathsO3 n ij == pathsrec 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> pathsO3 100 (4,4)
2422219241769802380469882122062019059350760968380804461263234408581143863208781993964800
(4.75 secs, 270708940 bytes)

The recursive algorithm would need in the order of 1077 years to arrive at this answer.

Still, pathsO3 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 pathsO3 into an O(n2) solution. Hint: there are more sums-of-products of independent values.

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.

HerbertDate: 2011-04-16T15:21Zx

Why the *2 after n12 and n56?

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

Twan van LaarhovenDate: 2011-04-18T19:26Zx

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.

HerbertDate: 2011-05-18T06:19Zx

This is so related!

Reply

(optional)
(optional, will not be revealed)
Name a function of type [[a]] -> [a]:
Use > code for code blocks, @code@ for inline code. Some html is also allowed.