Arboreal Isomorphisms from Nuclear Pennies

The game of Nuclear Pennies is a game for one player played on a semi-infinite strip of sites, each of which may contain any non-negative number of pennies. The aim is to achieve a target configuration of pennies from a starting configuration by means of penny fusion and penny fission. In penny fission, a penny is split into a pair of pennies in the two neighbouring sites. In penny fusion, two pennies, separated by exactly one site, are fused into one penny in the intervening site. Obviously no fission or fusion take place at the very end site because there is no room for one of the fission products or fusion precursors. The following diagrams show some legal moves in Nuclear Pennies:

One using fission:



One using fusion:



Can you achieve this target configuration:



Starting from this configuration?



It's a surprisingly non-trivial task but I urge you to have a go. If you give up, here's a video solution:



So a natural question arises: what single penny configurations can be achieved from that starting position? As every move in Nuclear Pennies is reversible, if you solved the puzzle you know we can shift a single penny 6 places left or right. (But to move 6 spaces left we need one extra space on the left.) So labelling positions 0,1,2,... we know that starting from 7 we can place a penny on 6n+1 for all non-negative n. But what about other positions?

We can borrow a technique from M Reiss's seminal paper on peg solitaire: Beiträge zur Theorie des Solitär-Spiels, Crelles Journal 54 (1857) 344-379. The idea is to associate an algebraic value with each position and then ensure that each this is chosen so that each legal move leaves the value unchanged.

Label the sites as follows:



We consider the value of each penny to be the value of its site and to get the value of a position, add the values of the individual pennies. If we pick x to be one of the complex cube roots of -1 then we have x=x²+1, because x³+1=(1+x)(x²+1-x). We define a move to be paralegal if it doesn't change the value of a position. It should be clear that penny fusion and fission are both paralegal. In other words

legal ⇒ paralegal

and

paraillegal ⇒ illegal

Now consider the values of 1,x,x²,... and so on. It's not hard to show that the only time an element in this sequence takes the value 1 is when the exponent is 6. In fact, x is a sixth root of unity. So this shows we can never shift a single coin by anything other than a multiple of 6 sites. As the video shows how to shift by 6 sites we know we can shift a single coin by n sites if and only if n is a multiple of 6 (and there's one site available to the left of the leftmost penny). It's neat that we took an excursion through the field of complex numbers to prove something about a seemingly simple piece of discrete mathematics.

But who cares about nuclear coin reactions? That's not my real motivation here. Believe it or not, each coin shuffle above corresponds to an isomorphism between certain types and solving the puzzle above actually demonstrates a really neat isomorphism between tuples of data structures.

Go back to the algebra above. We have that legal ⇒ paralegal because coin reactions correspond to legal manipulations of x. Fissioning a coin represents replacing x by 1+x², for example. But the converse isn't true. We can manipulate complex numbers in all kinds of interesting ways that don't correspond to coin reactions. For example, we can write -x² but we aren't allowed negative numbers of coins in a site. But there's another algebraic structure which corresponds exactly to coin reactions. To see that, it's time to start writing Haskell code:

> import Prelude hiding (Left,Right)
> import Control.Monad
> import Test.QuickCheck

> data T = Leaf | Fork T T deriving (Eq,Show)

T is a simple binary tree type. That declaration simply says that in the algebra of types, T=T²+1. Now unlike complex numbers, when you form a polynomial in the type T you can't have anything other than a non-negative integer as a coefficient. In other words, positions in the game of Nuclear Pennies correspond precisely to polynomials, in the algebra of types, in T.

For example, consider the first move in the video solution I gave:



It corresponds to

T⁷->T⁸+T⁶

or in Haskell notation an isomorphism

(T,T,T,T,T,T,T) -> Either (T,T,T,T,T,T,T,T) (T,T,T,T,T,T)

So you know what I'm going to do next, I'm going to code up the entire solution as a Haskell program mapping seven trees to one and back again. You don't need to read all of the following, but I thought I'd put it here for completeness. So skip over the code to the end if you feel like it...

Firstly, Either is a bit tedious to write. So here's my own implementation that uses a slightly non-standard way to write a type constructor:

> data a :+ b = Left a | Right b

We know that the type algebra is commutative, so A+B=B+A. But that '=' sign is really an isomorphism, not equality, and we'll need that isomorphism explicitly:

> commute :: a :+ b -> b :+ a
> commute (Left a) = Right a
> commute (Right b) = Left b

Same goes for associativity:

> associate :: (a :+ b) :+ c -> a :+ (b :+ c)
> associate (Left (Left a)) = Left a
> associate (Left (Right b)) = Right (Left b)
> associate (Right c) = Right (Right c)

And associativity the other way. Here I'm starting a convention of using primes to represent the inverse of a function.

> associate' :: a :+ (b :+ c) -> (a :+ b) :+ c
> associate' (Left a) = Left (Left a)
> associate' (Right (Left b)) = Left (Right b)
> associate' (Right (Right c)) = Right c

For a fixed b, '+b' is a functor so f :: a -> b induces a map a+b -> c+b. Unfortunately I don't know an easy way to make '+b', as opposed to 'b+', a functor in Haskell. So I'll have to make do with this function that lifts f to a+b:

> liftLeft  :: (a -> c) -> a :+ b -> c :+ b
> liftLeft f (Left a) = Left (f a)
> liftLeft f (Right b) = Right b

I'm going to need all these.

> liftLeft2 = liftLeft . liftLeft
> liftLeft3 = liftLeft . liftLeft2
> liftLeft4 = liftLeft . liftLeft3
> liftLeft5 = liftLeft . liftLeft4

I'm going to represent tuples a little weirdly, though arguably this is much more natural than the usual way to write tuples. Tn corresponds to Tⁿ:

> type T0 = ()
> type T1 = (T,T0)
> type T2 = (T,T1)
> type T3 = (T,T2)
> type T4 = (T,T3)
> type T5 = (T,T4)
> type T6 = (T,T5)
> type T7 = (T,T6)
> type T8 = (T,T7)

And now we need functions to assemble and disassemble trees at the start of a tuple. These correspond to fusion and fission respectively in the game.

> assemble (Left x) = (Leaf,x)
> assemble (Right (a,(b,x))) = (Fork a b,x)

> assemble' (Leaf,x) = Left x
> assemble' (Fork a b,x) = Right (a,(b,x))

Here's the first step in the solution. That was easy.

> step1 :: T7 -> T6 :+ T8
> step1 = assemble'

It gets successively harder because we need to push the assemble' function down into the type additions:

> step2 :: T6 :+ T8 -> T5 :+ T7 :+ T8
> step2 = liftLeft assemble'

> step3 :: T5 :+ T7 :+ T8 -> T4 :+ T6 :+ T7 :+ T8
> step3 = liftLeft2 assemble'

> step4 :: T4 :+ T6 :+ T7 :+ T8 -> T3 :+ T5 :+ T6 :+ T7 :+ T8
> step4 = liftLeft3 assemble'

> step5 :: T3 :+ T5 :+ T6 :+ T7 :+ T8 -> T2 :+ T4 :+ T5 :+ T6 :+ T7 :+ T8
> step5 = liftLeft4 assemble'

> step6 :: T2 :+ T4 :+ T5 :+ T6 :+ T7 :+ T8 -> T1 :+ T3 :+ T4 :+ T5 :+ T6 :+ T7 :+ T8
> step6 = liftLeft5 assemble'

We're going to use a few functions whose sole purposes is to shuffle around type additions using commutativity and associativity.

> swap23 :: (a :+ b) :+ c -> (a :+ c) :+ b
> swap23 = commute . associate . liftLeft commute

> shuffle1 :: T1 :+ T3 :+ T4 :+ T5 :+ T6 :+ T7 :+ T8 -> T1 :+ T4 :+ T3 :+ T5 :+ T6 :+ T7 :+ T8
> shuffle1 = liftLeft4 swap23

> shuffle2 :: (T1 :+ T4) :+ T3 :+ T5 :+ T6 :+ T7 :+ T8 -> T3 :+ T5 :+ T6 :+ T7 :+ T8 :+ (T1 :+ T4)
> shuffle2 = swap23 . liftLeft swap23 . liftLeft2 swap23 . liftLeft3 swap23 . liftLeft4 commute

After the left-to-right sweep, here's the right-to-left sweep:

> step7 :: T3 :+ T5 :+ T6 :+ T7 :+ T8 :+ (T1 :+ T4) -> T4 :+ T6 :+ T7 :+ T8 :+ (T1 :+ T4)
> step7 = liftLeft4 assemble

> step8 :: T4 :+ T6 :+ T7 :+ T8 :+ (T1 :+ T4) -> T5 :+ T7 :+ T8 :+ (T1 :+ T4)
> step8 = liftLeft3 assemble

> step9 :: T5 :+ T7 :+ T8 :+ (T1 :+ T4) -> T6 :+ T8 :+ (T1 :+ T4)
> step9 = liftLeft2 assemble

> step10 :: T6 :+ T8 :+ (T1 :+ T4) -> T7 :+ (T1 :+ T4)
> step10 = liftLeft assemble

> shuffle3 :: T7 :+ (T1 :+ T4) -> T1 :+ T4 :+ T7
> shuffle3 = commute

And now we have another left-to-right sweep followed by a right-to-left sweep:

> step11 :: T1 :+ T4 :+ T7 -> T2 :+ T0 :+ T4 :+ T7
> step11 = liftLeft2 (commute . assemble')

> step12 :: T2 :+ T0 :+ T4 :+ T7 -> T3 :+ T1 :+ T0 :+ T4 :+ T7
> step12 = liftLeft3 (commute . assemble')

> step13 :: T3 :+ T1 :+ T0 :+ T4 :+ T7 -> T4 :+ T2 :+ T1 :+ T0 :+ T4 :+ T7
> step13 = liftLeft4 (commute . assemble')

> step14 :: T4 :+ T2 :+ T1 :+ T0 :+ T4 :+ T7-> T5 :+ T3 :+ T2 :+ T1 :+ T0 :+ T4 :+ T7
> step14 = liftLeft5 (commute . assemble')

> shuffle4 :: T5 :+ T3 :+ T2 :+ T1 :+ T0 :+ T4 :+ T7 -> T7 :+ T5 :+ T3 :+ T2 :+ T1 :+ T0 :+ T4
> shuffle4 = liftLeft5 commute . liftLeft4 swap23 . liftLeft3 swap23 . liftLeft2 swap23 . liftLeft swap23 . swap23

> shuffle5 :: T7 :+ T5 :+ T3 :+ T2 :+ T1 :+ T0 :+ T4 -> T7 :+ T5 :+ T4 :+ T3 :+ T2 :+ T1 :+ T0
> shuffle5 = liftLeft3 swap23 . liftLeft2 swap23 . liftLeft swap23 . swap23

> step15 :: T7 :+ T5 :+ T4 :+ T3 :+ T2 :+ T1 :+ T0 -> T6 :+ T4 :+ T3 :+ T2 :+ T1 :+ T0
> step15 = liftLeft5 (assemble . commute)

> step16 :: T6 :+ T4 :+ T3 :+ T2 :+ T1 :+ T0 -> T5 :+ T3 :+ T2 :+ T1 :+ T0
> step16 = liftLeft4 (assemble . commute)

> step17 :: T5 :+ T3 :+ T2 :+ T1 :+ T0 -> T4 :+ T2 :+ T1 :+ T0
> step17 = liftLeft3 (assemble . commute)

> step18 :: T4 :+ T2 :+ T1 :+ T0 -> T3 :+ T1 :+ T0
> step18 = liftLeft2 (assemble . commute)

> step19 :: T3 :+ T1 :+ T0 -> T2 :+ T0
> step19 = liftLeft (assemble . commute)

> step20 :: T2 :+ T0 -> T1
> step20 = assemble . commute

And put it all together:

> iso :: T7 -> T1
> iso = step20 . step19 . step18 . step17 . step16 . step15 . shuffle5 . shuffle4 . step14 .
>     step13 . step12 . step11 . shuffle3 . step10 . step9 . step8 . step7 .
>     shuffle2 . shuffle1 . step6 . step5 . step4 . step3 . step2 . step1 

At this point I might as well write the inverse. This was mostly derived from the forward function using vim macros:

> step20' :: T1 -> T2 :+ T0
> step20' = commute . assemble'

> step19' :: T2 :+ T0 -> T3 :+ T1 :+ T0
> step19' = liftLeft (commute . assemble')

> step18' :: T3 :+ T1 :+ T0 -> T4 :+ T2 :+ T1 :+ T0
> step18' = liftLeft2 (commute . assemble')

> step17' :: T4 :+ T2 :+ T1 :+ T0 -> T5 :+ T3 :+ T2 :+ T1 :+ T0
> step17' = liftLeft3 (commute . assemble')

> step16' :: T5 :+ T3 :+ T2 :+ T1 :+ T0 -> T6 :+ T4 :+ T3 :+ T2 :+ T1 :+ T0
> step16' = liftLeft4 (commute . assemble')

> step15' :: T6 :+ T4 :+ T3 :+ T2 :+ T1 :+ T0 -> T7 :+ T5 :+ T4 :+ T3 :+ T2 :+ T1 :+ T0
> step15' = liftLeft5 (commute . assemble')

> shuffle5' :: T7 :+ T5 :+ T4 :+ T3 :+ T2 :+ T1 :+ T0 -> T7 :+ T5 :+ T3 :+ T2 :+ T1 :+ T0 :+ T4
> shuffle5' = swap23 . liftLeft swap23 . liftLeft2 swap23 . liftLeft3 swap23

> shuffle4' :: T7 :+ T5 :+ T3 :+ T2 :+ T1 :+ T0 :+ T4 -> T5 :+ T3 :+ T2 :+ T1 :+ T0 :+ T4 :+ T7
> shuffle4' = swap23 . liftLeft swap23 . liftLeft2 swap23 . liftLeft3 swap23 . liftLeft4 swap23 . liftLeft5 commute

> step14' :: T5 :+ T3 :+ T2 :+ T1 :+ T0 :+ T4 :+ T7 -> T4 :+ T2 :+ T1 :+ T0 :+ T4 :+ T7
> step14' = liftLeft5 (assemble . commute)

> step13' :: T4 :+ T2 :+ T1 :+ T0 :+ T4 :+ T7 -> T3 :+ T1 :+ T0 :+ T4 :+ T7 
> step13' = liftLeft4 (assemble .commute)

> step12' :: T3 :+ T1 :+ T0 :+ T4 :+ T7 -> T2 :+ T0 :+ T4 :+ T7 
> step12' = liftLeft3 (assemble . commute)

> step11' :: T2 :+ T0 :+ T4 :+ T7 -> T1 :+ T4 :+ T7 
> step11' = liftLeft2 (assemble . commute)

> shuffle3' :: T1 :+ T4 :+ T7 -> T7 :+ (T1 :+ T4) 
> shuffle3' = commute

> step10' :: T7 :+ (T1 :+ T4) -> T6 :+ T8 :+ (T1 :+ T4) 
> step10' = liftLeft assemble'

> step9' :: T6 :+ T8 :+ (T1 :+ T4) -> T5 :+ T7 :+ T8 :+ (T1 :+ T4) 
> step9' = liftLeft2 assemble'

> step8' :: T5 :+ T7 :+ T8 :+ (T1 :+ T4) -> T4 :+ T6 :+ T7 :+ T8 :+ (T1 :+ T4) 
> step8' = liftLeft3 assemble'

> step7' :: T4 :+ T6 :+ T7 :+ T8 :+ (T1 :+ T4) -> T3 :+ T5 :+ T6 :+ T7 :+ T8 :+ (T1 :+ T4) 
> step7' = liftLeft4 assemble'

> shuffle2' :: T3 :+ T5 :+ T6 :+ T7 :+ T8 :+ (T1 :+ T4) -> (T1 :+ T4) :+ T3 :+ T5 :+ T6 :+ T7 :+ T8 
> shuffle2' = liftLeft4 commute . liftLeft3 swap23 . liftLeft2 swap23 . liftLeft swap23 . swap23

> shuffle1' :: T1 :+ T4 :+ T3 :+ T5 :+ T6 :+ T7 :+ T8 -> T1 :+ T3 :+ T4 :+ T5 :+ T6 :+ T7 :+ T8 
> shuffle1' = liftLeft4 swap23

> step6' :: T1 :+ T3 :+ T4 :+ T5 :+ T6 :+ T7 :+ T8 -> T2 :+ T4 :+ T5 :+ T6 :+ T7 :+ T8 
> step6' = liftLeft5 assemble

> step5' :: T2 :+ T4 :+ T5 :+ T6 :+ T7 :+ T8 -> T3 :+ T5 :+ T6 :+ T7 :+ T8 
> step5' = liftLeft4 assemble

> step4' :: T3 :+ T5 :+ T6 :+ T7 :+ T8 -> T4 :+ T6 :+ T7 :+ T8 
> step4' = liftLeft3 assemble

> step3' :: T4 :+ T6 :+ T7 :+ T8 -> T5 :+ T7 :+ T8
> step3' = liftLeft2 assemble

> step2' :: T5 :+ T7 :+ T8 -> T6 :+ T8
> step2' = liftLeft assemble

> step1' :: T6 :+ T8 -> T7
> step1' = assemble

> iso' :: T1 -> T7
> iso' = step1' . step2' . step3' . step4' .
>     step5' . step6' . shuffle1' . shuffle2' . step7' .
>     step8' . step9' . step10' . shuffle3' . step11' .
>     step12' . step13' . step14' . shuffle4' . shuffle5' .
>     step15' . step16' . step17' . step18' . step19' . step20'

Enough already! Time to test the code. It's already pretty obvious that each of the steps above is invertible but let's use QuickCheck anyway. And here's a small puzzle: why to I have return Leaf twice?

> instance Arbitrary T where
>     arbitrary = oneof [return Leaf,return Leaf,liftM2 Fork arbitrary arbitrary]

> main = do
>     quickCheck $ \x -> x==iso (iso' x)
>     quickCheck $ \x -> x==iso' (iso x)

All done.

So time to think about what exactly we have here. We have a way to pack 7-tuples of trees into a single tree exactly. It's easy to store two trees in a tree say, just by making them the left and right branch of another tree. But then we get some 'wastage' because the isolated leaf doesn't correspond to any pair of trees. We have a perfect packing. Of course the set of trees and the set of 7-tuples of trees have the same cardinality, and it's not hard to find a bijection by enumerating these sets and laying them out side by side. But such a bijection could get arbitrarily complicated. But the isomorphism above only looks at the tops of the trees and can be executed lazily. It's a particularly nice isomorphism. And from what I showed above, you can only pack (6n+1)-tuples into a single tree. Weird eh?

And everything I say here is derived from the amazing paper Seven Trees in One by Andreas Blass. I've mentioned it a few times before but I've been meaning to implement the isomorphism explicitly for ages. Blass's paper also shows that you can, to some extent, argue paralegal ⇒ legal - something that's far from obvious.

Tries and their Derivatives

This was going to be part two of three, but I'm now squashing two parts down to this one article. So this is the conclusion to what I started in my previous installment.

Previously I talked about certain things you could do with types built up recursively using addition and multiplication of types - the so-called regular tree types. But I didn't talk about types that could be built using ->, ie. function types. These types seem somehow different. The type Either a b contains data pertaining to an a or a b, and (a,b) contains data about an a and a b. But a -> b is something completely different, a function rather than a datastructure. What I want to show is that this distinction isn't so clear, and that in fact functions can often be replaced by an isomorphic datastructure. There's an obvious application for this: memoisation, and most of what I say can be found in a paper by Ralf Hinze on that subject. However, I have an ulterior motive - computing generalised antidiagonals, and my notation will be biased accordingly.

Firstly, I'll be using conventional mathematical exponentiation notation for function types, so I'll use A^B to mean the same thing as B -> A, ie. functions mapping from B to A. This will help to make some of the algebra seem more natural.

So let's start with a type like A -> (). We could also write this as 1^A. From ordinary algebra we might expect this to equal 1. Sure enough, there is (up to pointwise equality) only one function of type A -> (), and it's called const ().

Now consider A^B+C. Algebraically this is just A^BA^C. In fact, the isomorphism is given by the standard Haskell either function. So far we can take any type B, built from singletons and addition, and reduce A^B to an isomorphic type with no exponentiation. Our goal will be to generalise this - for any type B (well...as many as possible anyway) to find a type constructor T such that T[A] = A^B, and in such a way that T is defined without exponentiation. Such a T is called a (generalised) trie. Like we did for the antidiagonal we can express this through a multiparameter type class:

> class Trie a t | a -> t where
>   tabulate :: (a -> b) -> t b
>   apply :: t b -> (a -> b)

tabulate converts a function into a datastructure and apply performs the opposite.

Again we're going to have the same issues with the antidiagonal. Haskell can't build tries for us automatically, but various generic forms of Haskell can. I'll be using plain Haskell so we'll be doing a bit of unnecessary labour. (I did try Hinze's Generics for the Masses approach but couldn't get it to work in this context. I suspect it can be made to work with more effort.)

So let's use the same example as in my previous installment, Bool.

> data BoolTrie a = BoolTrie a a deriving Show

> instance Trie Bool BoolTrie where
>   tabulate f = BoolTrie (f False) (f True)
>   apply (BoolTrie f _) False = f
>   apply (BoolTrie _ t) True = t

It's an easy exercise to show that apply and tabulate are mutual inverses. The BoolTrie stores the two possible values that a function Bool -> a could take. But here's the cool thing: given an expensive to compute function, f, of type Bool -> X, tabulate f is a datastructure that tells us everything we need to know about f. apply (tabulate f) gives us back our function again, but note how it never needs to call f more than once for each argument. In other words, we can define

> memo :: Trie a t => (a -> b) -> (a -> b)
> memo = apply . tabulate

which automatically converts functions into a memoised form.

Again, mimicking antidiagonals, let's implement the trie of Either a b:

> data EitherTrie u v a = EitherTrie (u a) (v a)

> instance (Trie a u,Trie b v) => Trie (Either a b) (EitherTrie u v) where
>   apply (EitherTrie f _) (Left x) = apply f x
>   apply (EitherTrie _ g) (Right y) = apply g y
>   tabulate f = EitherTrie (tabulate (f . Left)) (tabulate (f . Right))

Next consider products. This is slightly subtler, but only slightly. We use A^BC=(A^C)^B. If we define U[X]=X^B and V[X]=X^C then A^BC=U[V[A]].

> data PairTrie u v a = PairTrie (u (v a))

> instance (Trie a u,Trie b v) => Trie (a,b) (PairTrie u v) where
>   apply (PairTrie f) (b,c) = apply (apply f b) c
>   tabulate f = PairTrie $ tabulate (\a -> tabulate (\b -> f (a,b)))

This all applies recursively. So let's try tackling boolean lists again. We have L = 1+2L. Define T_L[X]=X^L. So T_L[X]=X¹X^2L=X(X^L)²=X(T_L[X])². This gives a nice bona fide definition of T_L. Just as with antidiagonals, it's annoying to have to write this by hand. A good generic programming language should be able to build a BoolTrie from a [Bool] automatically:

> data BoolListTrie a = BoolListTrie a (BoolListTrie a) (BoolListTrie a)

> instance Trie [Bool] BoolListTrie where
>   apply (BoolListTrie n _ _) [] = n
>   apply (BoolListTrie _ f _) (False:bs) = apply f bs
>   apply (BoolListTrie _ _ t) (True:bs) = apply t bs
>   tabulate f = BoolListTrie (f []) (tabulate (f . (False:))) (tabulate (f. (True:)))

Here's a silly example putting this in action on a binary version of the Ackermann function:

> inc (False : bs) = True : bs
> inc (True : bs) = False : inc bs
> inc [] = [True]

> dec [True] = []
> dec (True : bs) = False : bs
> dec (False : bs) = True : dec bs

> ack [] n = inc n
> ack m [] | not (null m) = ack (dec m) [True]
> ack m n = ack (dec m) (ack m (dec n))

> ack' = curry (memo (uncurry ack))

Note how once you've used ack' on some arguments, it runs instantly when you reuse it on the same or smaller arguments.

---

Part 2

Tries are containers. So suppose T is the trie corresponding to the type X with T[A] isomorphic to A^X. What is T[1]? Well it's just 1^X which from basic algebra equals 1. Or do look at it another way, it's the type of all functions taking values in the type 1. There's only one such function and so clearly there is a unique object of type T[1] and T[1] must be isomorphic to 1. An element of '1' can be thought of as a type with no information in it. So if we have a container with elements of 1 in it, it's as if the slots for those elements have simply been closed off. So another way to look at T[1] is that it is a container with all of its slots closed. So for tries, T[1] must be equivalent to the trivial type 1.

At this point, derivatives of containers have become mainstream. So I don't need to talk about this and can direct you to the Haskell wikibook if you need to brush up.

So what do we get if we differentiate the trie T? T is a container that holds one value for each element of X. The derivative, T', is a trie with a 'hole' in it. In other words, it's a container that contains a value for all elements of X but one. So what's T'[1]? All of the slots in the container have been stuffed with 1 and so are effectively 'closed'. But the hole is still there and the whole hasn't been 'closed'. So a T'[1] is a trie where one slot has been singled out as a hole. But the slots are in one-to-one correspondence with elements of X, and so T'[1]=X. There's another way to see this. Define the function on the reals t(y)=y^x. Then t'(y)=xy^x-1 so t'(1)=x. So it should come as no surprise that T'[1]=X. See the Container Types blog for more dicussion of this. They call F'[1] by the name log F. Their Y can be seen as the operator that maps a type to its trie.

Now, let's go back to the antidiagonal again. I showed how to compute X², but you can see that it's tricky to use the same approach to extend this to Xⁿ, for arbitrary naturals n. Instead, consider differentiating T' to make another hole in the trie. T''[X] must be a trie of X's with two holes, but the second hole obviously can't be where the first hole was. So T''[1] is a datastructure that consists of nothing but two distinct holes in a T[1]. As if by magic, T''[1] must be X². There's another way to see this. If t(y)=y^x, then t''(y)=x(x-1)y^x-2 so t''(1)=x². And more generally, T⁽ⁿ⁾[1]=Xⁿ. And that's the solution!

So a programming language powerful enough to differentiate types automatically, and powerful enough to build tries automatically (and Type-Indexed Data Types describes implementations of both of these in Generic Haskell) allows us to automatically construct the type of n-tuples where all of the elements are distinct. I'll leave the code to you as an exercise. :-)

Let me at least do the algebra for simple binary trees.

B = 1+B²

Let T be the trie of B so T[X] is isomorphic to X^B.

X^B = X(X^B)^B

T[X] = XT[T[B]]

(Exercise: try to draw a picture of this kind of tree. Because of the nesting of T's it's not a 'regular' type and you'll find yourself rapidly needing a lot of paper!).

So now we get

T'[X] = T[T[X]]+XT'[T[X]]T'[X]
T'[1] = T[1]+XT'[1]T'[1] = 1+X(T'[1])²

Ie. T'[1]=B

T''[X] = T'[T[X]]T'[X]+T'[T[X]]T'[X]+XT''[T[X]]T'[X]²
+XT'[T[X]]T''[X]
T''[1] = 2B²+B²T''[1]+BT''[1]

And that's a perfectly respectable recursive definition for T''[1]. I'll leave you to write the code. (It does work - this is actually the point which I started from...)

---

Some random thoughts

You can think of differentiation as annihilating an element of a container, leaving only a hole as the 'trace' that it was there before. So it's pretty weird that annihilating a 1 from a T[1] actually creates an X. It's almost as if T[1] is some kind of vacuum and X is an anti-particle. Weirder still, when we do this twice we get two distinct X's. So it's as if T[1] is some kind of vacuum for fermionic X's. This isn't completely insane. In quantum field theory we use differential operators to create and annihilate particles.

I haven't written the code, but it ought to be possible to go further and define Xⁿ/n! for instances of Ord. This is the type of n-tuples where the elements are distinct and in order. I expect you can construct this using a similar recursive technique to the way I built up the antidiagonal.

I don't know how to get an overline with HTML so I can't use Knuth's notation for this. But I expect that for instances of Ord there's a similar method that can be used to define X(X+1)...(X+n-1)/n!. This is the type of non-decreasing n-tuples, or, equivalently, n-element multisets. I'm not sure, but I think this may conflict with what Abbott et al. say in their paper on Derivatives of Containers where they want to identify exp (X) with multisets. I think it should be the sum of X(X+1)...(X+n-1)/n!, not the sum of Xⁿ/n!. But as that's just a sketch in the final paragraph of the paper, maybe I shouldn't get too worried about it. (If you're still reading Conor...)

Also, many identities satisfied by falling factorial should correspond to isomorphisms implementable in Haskell. And I'm sure there must be an Umbral Calculus connection.

And an important final note: X -> X² is not a functor. So don't even think about differentiating it.

---

Appendix

Changed my mind. Here's some code for implementing X²/2!.

Call the set of pairs (a,a') with a<a' the subdiagonal. An element of (a,a) either has a lower first element, a lower second element, or the two elements are equal. This is expressed in the type Pair below.

I'll sketch some algebra. Write s(X) for the subdiagonal. So s(X) is a little like X(X-1)/2. It's not hard to see that

s(1) = 0
s(a+b) = s(a)+ab+s(b)

s(ab) = s(a)b²+as(b)

If L[X] is a list of X's, then

L[X]=1+X L[X]

s(L[X]) = XL[X]+s(X)L[X]+Xs(L[X])

> import Test.QuickCheck

> data Pair a u = Lower u | Diagonal a | Upper u deriving Show

> class Ord a => SubDiagonal a u | a -> u where
>     twine :: (a,a) -> Pair a u
>     untwine :: Pair a u -> (a,a)
>     twine' :: (a,a) -> u
>     untwine' :: u -> (a,a)
>     twine (a,b) = case compare a b of
>         LT -> Lower (twine' (a,b))
>         EQ -> Diagonal a
>         GT -> Upper (twine' (b,a))
>     untwine (Lower u) = untwine' u
>     untwine (Diagonal a) = (a,a)
>     untwine (Upper u) = uncurry (flip (,)) $ untwine' u

> data SubEither a b u v = BothLeft u | Diff a b | BothRight v deriving Show

> instance (SubDiagonal a u,SubDiagonal b v) => SubDiagonal (Either a b) (SubEither a b u v) where
>     twine' (Left a,Left a') = BothLeft (twine' (a,a'))
>     twine' (Left a,Right b) = Diff a b
>     twine' (Right b,Right b') = BothRight (twine' (b,b'))
>     untwine' (BothLeft u) = let (a,a') = untwine' u in (Left a,Left a')
>     untwine' (Diff a b) = (Left a,Right b)
>     untwine' (BothRight u) = let (b,b') = untwine' u in (Right b,Right b')

> data SubPair a b u v = LeftDiffers u b b | LeftSame a v deriving Show

> instance (SubDiagonal a u,SubDiagonal b v) => SubDiagonal (a,b) (SubPair a b u v) where
>     twine' ((a,b),(a',b')) | a/=a'      = LeftDiffers (twine' (a,a')) b b'
>                            | otherwise  = LeftSame a (twine' (b,b'))
>     untwine' (LeftDiffers u b b') = let (a,a') = untwine' u in ((a,b),(a',b'))
>     untwine' (LeftSame a v) = let (b,b') = untwine' v in ((a,b),(a,b'))

> instance SubDiagonal Bool () where
>     twine' (False,True) = ()
>     untwine' () = (False,True)

> data SubList a u = FirstDiffers u (Pair [a] (SubList a u))| FirstSame a (SubList a u) | LeftNil a [a] deriving Show

> instance (SubDiagonal a u) => SubDiagonal [a] (SubList a u) where
>     twine' ([],(a:as)) = LeftNil a as
>     twine' ((a:as),(b:bs)) | a/=b       = FirstDiffers (twine' (a,b)) (twine (as,bs))
>                            | otherwise  = FirstSame a (twine' (as,bs))
>     untwine' (FirstDiffers u v) = let (a,a') = untwine' u in let (as,bs) = untwine v in (a:as,a':bs)
>     untwine' (FirstSame a u) = let (bs,bs') = untwine' u in (a:bs,a:bs')
>     untwine' (LeftNil b bs) = ([],b:bs)

> type Natural = Integer

> instance SubDiagonal Natural (Natural,Natural) where
>     twine' (a,b) = (a,b-a)
>     untwine' (a,b) = (a,a+b)

> main = do
>     quickCheck (\(x,y) -> untwine (twine (x,y))==(x::Bool,y))
>     quickCheck (\(x,y) -> untwine (twine (x,y))==(x::[Bool],y))
>     quickCheck (\(x,y) -> untwine (twine (x,y))==(x::(Bool,[Bool]),y))
>     quickCheck (\(x,y) -> untwine (twine (x,y))==(x::([Bool],Natural),y))
>     quickCheck (\(x,y) -> untwine (twine (x,y))==(x::[[(Bool,Natural)]],y))

It's interesting to play with this a little and see what actually gets constructed. For example

*Main> twine ([1..10],[1,2,3,6,5,6,6,8,9,10,11::Integer])
Lower (FirstSame 1 (FirstSame 2 (FirstSame 3 (FirstDiffers (4,2) (Upper
(FirstSame 5 (FirstSame 6 (FirstDiffers (6,1) (Upper (FirstSame 8
(FirstSame 9 (FirstSame 10 (LeftNil 11 [])))))))))))))

Spans of values that are the same in both lists are represented with just one set of values. Individual differences are flagged as such. And of course the whole thing works recursively if you have lists of lists as in the last quickCheck example above.

---

A Closing Puzzle

Can you define the type of unordered binary trees? An unordered binary tree is either the empty tree, or an unordered pair of unordered binary trees. I'd like to say T=1+T(T+1)/2 = 1+T+T², but those can't be turned directly into a valid definition is Haskell.

I don't know if there's a solution.

The Antidiagonal

In a programming language with constrained types we can construct a type like "the type of pairs of X's where the two X's are distinct". But can we make such a type in Haskell? Answering this question will take us on a long journey which I think I'm going to split over three parts. And as usual, just put this blog post into a file called main.lhs and you can run it with ghci -fglasgow-exts -fallow-undecidable-instances.

We can write our question a little more formally. Given a type X, can we form a type U with the property that

X² = U+X

The idea is that the = sign is an isomorphism with the property that the diagonal in X², ie. elements of the form (x,x), get mapped to the right component of U+X. When we are able to do this, we'll call U the antidiagonal of X, and say that X is splittable.

We can express the relationship between U and X through a multiparameter type class

> diagonal x = (x,x)

> class Eq x => AntiDiagonal u x | x -> u where
>   twine :: (x,x) -> Either x u
>   untwine :: Either x u -> (x,x)

>   twine' :: (x,x) -> u
>   untwine' :: u -> (x,x)

>   twine (x,y) | x==y      = Left x
>               | otherwise = Right $ twine' (x,y)

>   untwine (Left x)  = (x,x)
>   untwine (Right y) = untwine' y
    

The isomorphism between X² and X+U is given by twine and untwine. But to save writing similar code over and over again, and to ensure we really are mapping the diagonal of X² correctly, we define these in terms of twine' and untwine'. (Note that twine' is partial, it's only guaranteed to take a value off the diagonal.)

Just to get into the swing of things, here's a simple example:

> instance AntiDiagonal Bool Bool where
>   twine' (a,b) = a
>   untwine' a = (a,not a)

It's not hard to check that twine and untwine are mutual inverses and that twine . diagonal=Left. You can view this as a special case of 2²=2+2 as 2 is essentially a synonym for Bool.

It looks a lot like what we're trying to do is subtract types by forming X²-X or X(X-1). Much as I've enjoyed trying to find interpretations of conventional algebra and calculus in the context of the algebra of types, subtraction of types, in general, really doesn't make much sense. Consider what we might mean by X-1. If X=Bool then we could simply define

data Bool' = False'

This certainly has some of the properties you might expect of Bool-1, such as having only one instance. But it's not very natural. How should we embed this type back in Bool? There are two obvious ways of doing it and neither stands out as better than the other, and neither is a natural choice. But what we've shown above is that there is a natural way to subtract X from X² because a copy of X appears naturally in X² as the diagonal. So the question is, can we extend this notion to types beyond Bool?

How about a theorem. (I'll give a more intuitive explanation below.)

Theorem

In any commutative semiring fix a and b. If the equation a²=u+a has a solution, and the equation b²=v+b has a solution, then the equations (a+b)²=x+(a+b) and (ab)²=y+ab also have solutions.

Proof

Simply define

x = u+v+2ab

and

y = au+bv+uv.

QED

So if some types are splittable, so are their sums and products. As types form a commutative semiring we see that this theorem, twined with 1²=1+0 allows us to form "X(X-1)" for any type X built purely using non-recursive Haskell data declarations. In fact, we can use the above theorem to define "X(X-1)" for types, and use the notation X² for this. (I hope your browser shows that the exponent here is underlined.) There's a reason I use this notation which I'll get to later. So (a+b)²=a²+b²+2ab and (ab)²=ab²+ba²+a²b².

Here's a more intuitive explanation of that theorem above. Suppose u and v are both in Either a b. Then there are several ways they could be different to each other. For example, they could both be of the form Left _ in which case u = Left u' and v = Left v' and u' and v' must be distinct. They could be of the form Left u' and Right v' in which case it doesn't matter what u' and v' are. A little thought shows there are four distinct cases in total. These can be written in algebraic notation as u+v+ab+ba=u+v+2ab. To make this clearer, here's an implementation:

> data Either' a b u v = BothLeft u | BothRight v | LeftAndRight Bool a b

> instance (AntiDiagonal u a,AntiDiagonal v b) => AntiDiagonal (Either' a b u v) (Either a b) where
>   twine' (Left x,Left y)   = BothLeft (twine' (x,y))
>   twine' (Right x,Right y) = BothRight (twine' (x,y))
>   twine' (Left x,Right y)  = LeftAndRight False x y
>   twine' (Right x,Left y)  = LeftAndRight True y x

>   untwine' (BothLeft u) = let (a,b) = untwine' u in (Left a,Left b)
>   untwine' (BothRight v) = let (a,b) = untwine' v in (Right a,Right b)
>   untwine' (LeftAndRight False x y) = (Left x,Right y)
>   untwine' (LeftAndRight True x y) = (Right y,Left x)

A similar argument can be carried through for (a,b) leading to:

> data Pair' a b u v = LeftSame a v | RightSame u b | BothDiffer u v deriving (Eq,Show)

> instance (AntiDiagonal u a,AntiDiagonal v b) => AntiDiagonal (Pair' a b u v) (a,b) where
>   twine' ((a,b),(a',b')) | a==a'     = LeftSame a (twine' (b,b'))
>                          | b==b'     = RightSame (twine' (a,a')) b
>                          | otherwise = BothDiffer (twine' (a,a')) (twine' (b,b'))

>   untwine' (LeftSame a v)   = let (b,b') = untwine' v in ((a,b),(a,b'))
>   untwine' (RightSame u b)  = let (a,a') = untwine' u in ((a,b),(a',b))
>   untwine' (BothDiffer u v) = let (a,a') = untwine' u
>                                   (b,b') = untwine' v
>                                   in ((a,b),(a,b'))        

This all works very well, but at this point it becomes clear that Haskell has some weaknesses. The type Either () () is isomorphic to Bool so we should be able to use the above code to construct the antidiagonal of Bool automatically. But Haskell doesn't give us access to that information. We can't ask, at runtime, if Bool is the sum of more primitive types. There are a number of solutions to this problem - we can use various more generic types of Haskell, or use Template Haskell. But I'm just going to stick with Haskell and manually construct the antidiagonal. (I wonder what Coq has to offer here.)

So I've solved the problem for 'finite' types built from singles, addition and multiplication. But what about recursive types. Before doing that, let's consider an approach to forming the antidiagonal of the naturals. Haskell has no natural type but let's pretend anyway

> type Natural = Integer

There's an obvious packing of a distinct pair of naturals as a pair of naturals:

> instance AntiDiagonal (Bool,Natural,Natural) Natural where
>   twine' (a,b) | compare a b == GT = (False,a-b,b)
>                | compare a b == LT = (True,b-a,a)
>   untwine' (False,d,b) = (b+d,b)
>   untwine' (True,d,a) = (a,a+d)

(I had to use compare to work around a blogger.com HTML bug!) It'd be cool if the code above could have been derived automatically, but alas it's not to be. But we can get something formally similar. Define the natural numbers like this:

> data N = Zero | S N deriving (Eq,Show)

A natural is zero or the successor of a natural. Algebraically this is just N=1+N. Now we wish to find N². Using 1²=0 and the earlier theorem we get M=1²+N²+2N, ie. M=M+2N. Let's code this up:

> data M = Loop M | Finish Bool N deriving Show

> instance AntiDiagonal M N where
>   twine' (Zero,S x) = Finish False x
>   twine' (S x,Zero) = Finish True x
>   twine' (S x,S y) = Loop (twine' (x,y))
>   untwine' (Finish False x) = (Zero,S x)
>   untwine' (Finish True x) = (S x,Zero)
>   untwine' (Loop m) = let (a,b) = untwine' m in (S a,S b)

Note that I've more or less just coded up the same thing as what I did for Either above. Can you see that that this is a disguised version of the code for the antidiagonal of Natural above? Think about the type N=1+N. We can view this as the set of paths through this finite state diagram, starting at N and ending at 1:



Essentially there's just one path for each natural number, and this number counts how many times you loop around. Now consider the same sort of thing with the type M=M+2N, starting at M and ending at 1:



We can describe such paths by the number of loops taken at M, the number of loops taken at N, and a Bool specifying whether we took path 0 or 1 from state 2 to state N. In other words, M is much the same thing as (Bool,Natural,Natural) above! M is a kind of 'compressed' version of a pair of N's. Suppose we want to twine S (S (S Zero)) and S (S (S (S Zero))). Both of these share a S (S (S _)) part. What the type M does is allow you to factor out this part (that's the part that goes into Loop) and the remainder is stored in the Finish part, with a boolean specifying whether it was the first or second natural that terminated first.

Let's step back or a second. Earlier I showed how for any type X, built from addition, multiplication and 1, we could form X². Now we've gone better, we can now form X² even for recursive types. (At least for data, not codata.) We haven't defined subtraction in general, but we have shown how to form X(X-1) in a meaningful way.

Let's try lists of booleans next, [Bool]. We can write this type as L=1+2L. Let's do the algebra above to find P = L²:

P = L²
= (1+2L)²
= 1²+(2L)²+4L
= (2L²+2²L+2²L²)+4L
= 2P+2P+2L+4L
= 4P+6L

In other words P = 4P+6L. The easiest thing is to code this up. Remember that each '=' sign in the above derivation corresponds to an isomorphism defined by twine and untwine so after lots of unpacking, and rewriting 4P+6L as 2L+2L+2P+2L+2P we get

> data SharedList = LeftNil Bool [Bool] | RightNil Bool [Bool]
>                   | HeadSame Bool SharedList | TailSame Bool [Bool]
>                   | Diff Bool SharedList deriving Show

> instance AntiDiagonal SharedList [Bool] where
>   twine' ([],b:t) = LeftNil b t
>   twine' (b:t,[]) = RightNil b t
>   twine' (a:b,a':b') | a==a' = HeadSame a (twine' (b,b'))
>                   | b==b' = TailSame a b
>                   | otherwise = Diff a (twine' (b,b'))

>   untwine' (LeftNil b t)  = ([],b:t)
>   untwine' (RightNil b t) = (b:t,[])
>   untwine' (HeadSame a b) = let (t1,t2) = untwine' b in (a:t1,a:t2)
>   untwine' (TailSame a b) = (a:b,not a:b)
>   untwine' (Diff a b)     = let (t1,t2) = untwine' b in (a:t1,not a:t2)

This looks pretty hairy, but it's really just a slight extension of the M=M+2N example. What's happening is that if two lists have the same prefix, then SharedList makes that sharing explicit. In other words this type implements a form of compression by factoring out shared prefixes. Unfortunately it took a bit of work to code that up. However, if we were programming in generic Haskell, the above would come absolutely for free once we'd defined how to handle addition and multiplication of types. What's more, it doesn't stop with lists. If you try it with trees you automatically get factoring of common subtrees and it works with any other datatype you can build from a Haskell data declaration (that doesn't use builtin primitive types like Int or Double).

So now I can say why I used the underlined superscript notation. It's the falling factorial. More generally Xⁿ=X(X-1)...(X-n+1). You may be able to guess what this actually means - it's an n-tuple with n distinct elements. Unfortunately, you can probably see that generalising the above theorem to n from 2 gets a bit messy. But (and this is what I really want to talk about) there's an amazing bit of calculus that allows you to define a generating function (more like generating functor) that gives you all of the Xⁿ in one go. But before I can talk about that I need to write a blog about generalised tries...

---

Here are some exercises.

(1) Can you code up the antidiagonal of binary boolean trees, T = 2+T²:

data BoolTree = Leaf Bool | Fork BoolTree BoolTree

(2) There are more efficient ways to define the naturals than through the successor function. Can you come up with a more efficient binary scheme and then code up its antidiagonal?

(3) The antidiagonal of the integers can be approximated by (Integer,Integer). This seems a bit useless - after all, the whole point of what I've written above is to split this up. But we can use this approximation to construct approximations of other types where you do get a payoff. Implement an approximation to [Integer]² this way so that you still get the benefit of prefix sharing. This looks a lot like traditional tries.

I'll have a Buchburger with fries

This is, as usual, a literate Haskell program. So let me start with

> import Data.Map as M
> import Data.List as L

Everyone else is at it, so why not me? Here's my approach to solving this problem. Let's warm up with a simpler menu. Consider:

Fries $3
Burger $5

What can I get that costs exactly $13?

My approach is going to revolve around rewrite rules. I'm going to find a set of rewrite rules that converts money into fries and burgers, and does so in such a way that it always eliminates the money, if it is possible to do so.

Let's use D to represent a dollar, F to represent fries and B to represent a burger. Here's our first attempt at some rewrite rules:

DDD->F
DDDDD->B

Now we can take the string DDDDDDDDDDDDD and try to rewrite it. We could start with the first rule:
DDDDDDDDDDDDD
->DDDDDDDDDDF
->DDDDDDDFF
->DDDDFFF
->DFFFF
But now we're stuck with a D and no way to get rid of it.

Before we go much further, however, it's time to simplify our notation. Firstly, unlike normal rewrite rules, we don't care about the order of our symbols. So let's assume that we can reorder symbols any time we want in our strings and rules. So FDD=DFD=DDF and so on. Let's also use Xⁿ to mean n copies of X in a row. Now we can write our strings as mononomials in D, F and B. That makes things much easier. We can now write our rules as

D³->F
D⁵->B

We've seen that just picking the first rule each time can get us stuck. It seems that the solution should be to come up with an algorithm to generate the correct order. But there's another strategy - it's to modify the rules so that the order doesn't matter. A set of rules with this property is said to be confluent. Amazingly you can turn any reasonable set of rewriting rules into a confluent one merely by adding in elements.

So let's define a monomial type:

> newtype Monomial a i = Monomial (Map a i) deriving (Eq,Show)
> unM (Monomial x) = x

> divides (Monomial a) (Monomial b) = isSubmapOfBy (<=) a b
> lcm (Monomial a) (Monomial b) = Monomial $ unionWith max a b

> instance (Show a, Show i, Num i,Ord a) => Num (Monomial a i) where
>     Monomial a*Monomial b = Monomial $ unionWith (+) a b

> instance (Show a, Show i, Num i,Ord a) => Fractional (Monomial a i) where
>     Monomial a/Monomial b = Monomial $ M.filter (/=0) $ unionWith (+) a (M.map negate b)

Rewrite rules are simply pairs of monomials:

> data Rule a i = Monomial a i :-> Monomial a i deriving (Eq,Show,Ord)
> lm (a :-> _) = a

Now we can write code to apply these rules to a string. The first function simply applies the rule to a term:

> apply (p :-> q) r = r/p*q

applicable searches for the first applicable rule in a list of rules:

> applicable r rules = find (\(p :-> _) -> divides p r) rules

And reduceStep applies it. Note that it returns the result twice in a pair. This is simply so that it works with the next function. In the event that no rule is applicable it returns Nothing.

> reduceStep rules r = applicable r rules >>= return . dup . flip apply r
>     where dup x = (x,x)

reduce repeatedly applies reduceStep until it can be applied no more. Repeatedly doing something until Nothing is returned is best handled by unfoldr:

> reduce rules r = last $ (r:) $ unfoldr (reduceStep rules) r

Now consider our example rewrite rules. Let's find a simple example of where confluence fails to hold. Consider two different ways to rewrite D⁵. We could use the first rule and get D⁵->D²F or use the second and get D⁵->B. Either way we terminate - and get different results. How can we get out of this impasse? There's a really obvious solution, add in a new rule D²F->B. Now if we start with D⁵ and apply the first rule then we can get out of it by applying this new third rule. Either way, we get B.

Before we generate new rules, we need to say something about how we can ensure termination. If we have one rule that says X->Y and another that says Y->X then we're going to get stuck. So what we can do is insist on an ordering on our monomials. We'll restrict ourselves to rules where the left hand side is larger than the right hand side. So here's an ordering known to algebraists as lexicographic ordering. The idea is that the symbols in the monomial have an ordering on them. One monomial, a, is considered bigger than another, b, if the power of the least symbol (in the symbol ordering) in a/b is negative. (BTW See postscript for a comment on this.) What this means is that the rewrite rules will try to eliminate symbols earlier in the ordering and replace them with symbols that come later. By making our currency symbol the first in the ordering the rewrite rules will tend to act so as to eliminate the currency. And that's what the original problem asked us to do. Note that there are other orderings we could use so for generality Monomial ought to be a typeclass.

> instance (Show a,Ord a,Ord i,Num i) => Ord (Monomial a i) where
>     compare x y = let
>         Monomial u = x/y
>         in
>         if M.null u
>             then EQ
>             else if snd (M.findMin u)>0
>                 then LT
>                 else GT

We also use a function to ensure that any rule is in canonical form with the greater term on the right:

> order (a :-> b) = if a>b then b :-> a else a :-> b

If we find a complete confluent set of rules that respect the ordering in this way, how do we know the rules terminate? Suppose the complete ordered list of symbols in our rules is (a₁,a₂,...,a_n) then any string a₁^k₁a₂^k₂...a_n^k_n can be written as a vector (k₁,k₂,...,k_n). Each application of a rewrite rule replaces the corresponding vector with one that is lexicographically earlier. It's pretty clear that such a process must terminate. (Aside: You can think of this as using the fact that ωⁿ is well-ordered. More general kinds of rewrite rules make use of more complex ordinals to prove termination. This goes to show that ordinal arithmetic is actually a useful tool.)

So back to what we were doing: finding how we can add new rules to ensure confluence. Suppose we have two rules whose LHS 'overlap'. If we think in terms of monomials then we mean two monomials that have a non-trivial common factor. If we form the lowest common multiple of these monomials then we now have a 'minimal' term to which both rules can be applied. In fact, that's how I got the example above. One rule has P³ on the LHS and one has P⁵. The lowest common multiple of these is P⁵ and I investigated how the two rules apply to this one term. When we did this, we ended up with two different terms which we simply made into a new rule. Here it is in Haskell:

> newRule pq@(p :-> _) pq'@(p' :-> _) = let
>     j = Main.lcm p p'
>     in order $ apply pq j :-> apply pq' j

Given a rule, we also simplify the left and right hand sides with normal:

> normal rules (p :-> q) = let
>     p' = reduce rules p
>     q' = reduce rules q
>     in if p'/=q'
>         then Just $ order $ reduce rules p :-> reduce rules q
>         else Nothing

The problem is that when we generate a new rule there now appear a whole bunch of new pairs that could be considered because you need to consider pairs made up of the new rule and all of the old rules. Our strategy for dealing with this will be simple: take the list of all pairs of rules and generate new rules from those. Every time we generate a new rule we need to add in the new rule pairs that are generated. We maintain a list of pairs found so far, chewing pairs for newRule from the front, and adding newly generated pairs to the back.

The first line is easy:

> complete incomplete = let

We make all possible pairs. We don't want to consider pairs twice so we use the inherited ordering on pairs:

>     pairs = [(x,y) | x <- incomplete, y <- incomplete, y>x]

Here's a little heuristic. Things go faster if we use monomials that come lower in the ordering before higher ranking monomials. So the completed set of rules is going to be sorted:

>     rules = sort incomplete

Now here's a weird bit to explain. This isn't essential to the mathematics but is is essential to this piece of code. The local function iterate is going to return a pair that consists of the remaining list of rule pairs and the completed rules. But iterate also consumes rule pairs. So this function is actually going to consume its own output. There's an important issue with this: eventually the rate of new rule creation will peter out so that it will consume its way to the end. At this point iterate will try to consume the very element it's expected to generate and spin in an infinite loop. (The dangers of unguarded recursion!) So we also maintain another list whose length is the same as the distance between the front and back of our pair list. When this first list is empty, we know we can return [] for the second list. It doesn't matter what this list is, as long as it has the correct length. This 'phantom' list is the first argument to iterate. But for now I suggest ignoring it and coming back to it if you ever get to a second reading.

Here's the base case I mentioned:

>     iterate [] pairs rules = ([],rules)

And here is the guts of the loop:

>     iterate (_:n) (p:pairs) rules = 

p is the first pair chewed off the front. We generate a new rule from it and normalise it:

>         case normal rules $ uncurry newRule p of

If we get a new rule we now need to add it to our list of rules, and add in all the new pairs it generates:

>             Just pq -> let u = [(x,pq) | x <- rules, x /= pq]

Note the use of insert to maintain the order.

>                            (pairs', rules') = iterate (u ++ n) pairs $ L.insert pq rules
>                            in (u ++ pairs', rules')

Otherwise we don't need to generate anything:

>             Nothing -> iterate n pairs rules

Here's where we tie the knot by feeding our output back in as input:

>     (pairs', rules') = iterate pairs (pairs ++ pairs') rules
>     in rules'

It's not obvious at this point that the algorithm terminates. It does, but you'll have to trust me.

Anyway, that's not all. We also do a bit of tidying up by eliminating redundant rules. In a confluent system, if there are two rules we can apply, it doesn't matter which one we choose. Now consider a pair of rules looking something like: AB²C -> ... and AB³C² -> ... . If we can apply the second rule, then we could also apply the first rule. But in a confluent system, it doesn't matter which rule we choose. So we can just throw out the second rule here. So we tidy up our rules by eliminating rules made redundant by another. Turns out that at this stage we have a canonical form. For any given set of rules, any completion with the reducant elements deleted is the same as any other (modulo reordering). So our completion process is itself confluent! This isn't hard to prove but I'm not going to do it now.

> nonredundant rules r = L.null [s | s <- rules, s/=r, lm s `divides` lm r]
> tidy rules = L.filter (nonredundant rules) rules

And now canonical puts our rules into a complete canonical form:

> canonical rules = tidy $ complete rules

Back to the original problem:

> data Menu = P 
>     | MixedFruit 
>     | FrenchFries 
>     | SideSalad 
>     | HotWings 
>     | MozzarellaSticks 
>     | SamplerPlate deriving (Eq,Ord,Show)

Note that deriving Ord. This is so that the lexicographic ordering above can work.

> makeRule a b = (Monomial $ M.fromList a) :-> (Monomial $ M.fromList b)
> rules = [
>         makeRule [(P,215)] [(MixedFruit,1)],
>         makeRule [(P,275)] [(FrenchFries,1)],
>         makeRule [(P,335)] [(SideSalad,1)],
>         makeRule [(P,355)] [(HotWings,1)],
>         makeRule [(P,420)] [(MozzarellaSticks,1)],
>         makeRule [(P,580)] [(SamplerPlate,1)]     
>     ]

> crules = canonical rules

> ex n = Monomial $ fromList [(P,n)]
> example0 = reduce crules (ex 1505)

> main = print (reduce crules (ex 1505))

Now, some caveats. I don't claim that this is a fast method, just a different approach. Completing a set of rules is an expensive operation and if you start with N rules, the worst case running time is something horrendous like O(exp(exp(N)). On the other hand, once you've found a set of canonical rules, reduction with them is very fast. This makes it ideal for problems involving repeated use of the same menu. Or for when you have a lot of spare cash:

> example = print (reduce crules (ex 10000115))

And now I have a confession to make. I've really just described the Buchberger algorithm for Gröbner basis computation. It's not quite the Buchberger algorithm because I'm only considering binomials, not polynomials - but you'll notice my code is remarkably similar to the code at the other end of that link. The rule p :-> q represents the binomial p-q. Unfortunately, the Buchberger algorithm is usually presented in its more general form - but understanding that requires a background in commutative algebra and algebraic geometry. I hope I've managed to present a form of it in a way that's digestible by people without abstract algebra experience. The binomial Buchberger algorithm has many applications, including a variety of applications in integer programming.

Note that this is also very close to the Knuth-Bendix algorithm. However, the latter algorithm is for the noncommutative case, and it isn't guaranteed to terminate.

BTW I recommend playing with the completed rules by hand - at least for simple examples. It's interesting to see what steps are taken to reach the solution.

PS My monomial ordering is the opposite of the usual lexicographic monomial ordering. This is because I want my rules sorted with what is usually considered the largest monomial at the front.

PPS Can I get rid of that weird phantom list thing? You can replace the list with an integer for its length, but it's not clear to me that this is an actual improvement.

Update: I forgot to mention this. Check out David Amos's implementation of the full Buchberger algorithm which you can use to solve many more problems than my simplistic monomial version.

Data and Codata

In Gödel, Escher, Bach, Hofstadter introduces the programming language Bloop. Bloop is a little like BASIC except that it forces programmers to specify in advance how many times each loop will iterate. As a result, Bloop programs have the handy property that they are guaranteed to terminate. Unfortunately this property also makes it impossible to write something like an operating system in Bloop. To write an OS we need open-ended loops that keep running until the user explicitly chooses to shut the OS down. One solution to this is to write code in Floop. Floop allows us to write unbounded loops, the equivalent of C's while (1) { ... }. The problem with that, however, is that we can write runaway infinite loops that never terminate and never give us any output. Is there some language that lies between Bloop and Floop that can give us unbounded looping when we need it, but which never allows us to hoist ourselves by our petards by writing runaway loops?

Wishing no disrepsect to OS writers, at first blush it might seem that the distinction between a runaway loop and an idle OS is too fine - if we can write an infinite loop that does something useful, then surely we can write a useless one too. But it turns out that there is a very elegant and well-principled way to distinguish between these kinds of loops, and this allows us to write open-ended interactive software in a programming language that nonetheless always produces a well-defined output, no matter what the input. In order to do this we need to distinguish between two kinds of data: data and codata. By ensuring that a function expecting codata never receives data, and vice versa, we can ensure that even programs with open-ended loops always produce a well defined output.

The concepts I want to talk about are very general and can apply to whatever programming language you use. I'm going to use some simple Haskell examples but most of these will translate to other languages. So consider something like

sum [] = 0
sum (a:as) = a + sum as

This sums the elements of a list. Note how it's well behaved as long as we give it a finite list as input. On the other hand, consider

sum' [] = 0
sum' a = sum' (1:a) - 1

This isn't well behaved at all. Except when you input empty lists, it never gives a result. From a mathematical perspective it's not a good definition either, there are many functions that satisfy these two properties. Is there some general principle at work here that allows us to tell immediately that one of these terminates and the other doesn't? We know from Turing that there is no procedure that guarantees we can always tell such programs apart, but in this case there is something that we can easily point to. In the first program, the right hand side of the second line uses the sum function recursively but we only apply it to a strict subpart of the input, its tail in fact. In the second example we apply sum' to something that contains the argument. The former function is using what is known as structural recursion, and it's not hard to convince yourself that structural recursion, applied to finite datastructures, always terminates.

So if we limit ourselves to structural recursion we can guarantee our programs will always terminate. But what about a definition like:

fact 0 = 1
fact n = n * fact (n-1)

That doesn't appear to use structural recursion. But we can view it as such like this. Define the natural numbers as follows:

data Nat = Zero | S Nat

0 is represented as Zero, 1 is represented as S Zero and so on. We can represent every natural this way. Here's the important thing: if n>0 then n-1 is simply a subpart of n. So we can view this kind of recursion as a type of structural recursion.

(In fact, by a curious quirk of the Haskell 98 standard we can rewrite our definition to look more like a structural recursion:

fact' 0 = 1
fact' (n+1) = (n+1) * fact' n

I'm guessing this feature is in Haskell precisely so that people can 'pretend' they are using structural recursion with the + in n+1 serving a role as a kind of pseudo-constructor.)

So we have a nice rule for ensuring our code terminates. But sum fails when applied to infinite lists. Should we simply rule out infinite datastructures? That seems a bit drastic. The example that convinced me to look into Haskell was

fib = 1 : 1 : zipWith (+) fib (tail fib)

We really don't want to rule out such a succinct definition of the Fibonacci numbers. But how can we allow such structures when we have functions like sum sitting around. Applying sum to fibs will obviously cause non-termination.

Let's consider another example:

sumSoFar x [] = [x]
sumSoFar x (y:ys) = x : sumSoFar (x+y) ys

Like sum, this fails to terminate for an infinite input. But unlike sum, it's possible to make sense of it. If the inputs were 0 and the infinite list [1,1,1,1,...] then the result would be [0,1,2,3,...]. The program might not terminate, but from a mathematical perspective this is a completely well defined function. What's more, suppose the input list represented a stream of data being input at a keyboard, and that the output was displayed on a screen one element at a time, then we'd have a simple cash register. This program might not terminate, but it's completely well behaved. Note that this could only work in a lazy language. A strict language would want to evaluate the entire list before outputting anything. But in a lazy language we can start outputting the beginning of the list before the rest of it is computed.

From the point of view of using infinite lists, it's sum that's badly behaved, and sumSoFar that's well behaved. Again, can we point to a simple distinction between these two programs that explains this? Turns out we can, and in some sense it's dual to what we said above. sumSoFar is well behaved because when we recursively call sumSoFar on the right hand side we do so from inside a list constructor. (Remember that : is the list constructor.) This is known as guarded recursion and it guarantees that even though our programs don't terminate, they still define a unique mathematical function and result in a well behaved program. In the case of sumSoFar, each time we look at another element of the result, we trigger another lazy evaluation that terminates. But the right hand side won't just run off and compute endlessly until we do that triggering because the recursion is 'guarded' within a constructor. (That, by the way, was the main point of this article, so you can probably relax now.)

Note the duality: in structural recursion we 'deconstruct' the argument and then we're allowed to recurse. In guarded recursion we recurse first, and then we're allowed to use the constructor.

So we've almost achieved our goal of describing rules to allow open-ended loops because we've managed to give a rule for writing functions that are guaranteed to be well-behaved even though they manipulate infinite data structures. But we're not quite home yet - we still can't have functions like sum coexist with infinite lists. How can we ensure that an infinite list is never handed to sum?

Consider the definition of Haskell lists. It's something like this:

data [a] = [] | a : [a]

I.e. a list of a's is either the empty list [] or it's made from an a and a list of a's.
You can think of this as an equation in [a]. In Haskell we take this as uniquely defining what [a] is, but in reality there is more than one solution to this equation. Consider the type consisting of only finite lists. That satisfies this equation. A finite list is either an empty list, or an element followed by a finite list. Similarly a possibly infinite list is either an empty list, or an element followed by a possibly infinite list. There is an ambiguity. Finite lists form, what is in some sense, the smallest possible solution to this equation. The possibly infinite lists form the largest possible solution. Haskell takes the largest possible solution.

But suppose we were to distinguish between these two different solutions. We could use the keyword data to mean the smallest solution and codata to mean the largest solution. The former represents data, and it's always finite. The latter represents what we call codata, and it's possibly infinite. And now we can refine our rules for well-behavedness. Consider data and codata to be distinct types. In a strongly typed language this immediately allows us to restrict sum to data, not codata. The rule is: you're only allowed to use structural recursion with data and guarded recursion with codata. With that rule, we're guaranteed that our recursions will always be safe, and yet that we can still have open-ended loops in our code. Sometimes these are called recursion and corecursion respectively.

And now we can go a little crazy with language. When we want to prove that a structurally recursive program terminates we use induction. This doesn't work straightforwardly for corecursion, so instead we use a principle called coinduction. Recursive programs typically terminate. Corecursive programs don't necessarily terminate, but they're still well-behaved as long as they keep on going whenever we give them input. We can call this cotermination. And so on... I'm not going to say what coinduction is because I'd have to talk about bisimulation, and this post would get way too long.

So now a mathematical aside. The reason for all the "co"s is that data and codata have categorical descriptions and they turn out to be dual to each other. But you don't hear mathematicians talking about corecursion and coinduction much. Why not? Well one of the axioms of set theory is the Axiom of Foundation. One way of interpreting this is the statement that there is no infinite sequence

...a₃∈a₂∈a₁∈a₀.

So even though we can construct infinite lists in mathematics, we can't construct 'infinitely deep' lists. This means that in mathematics we can use a form of structural recursion. In fact, the familiar principle of induction follows from this. So for many of the things that mathematicians do, induction works fine. But if we decide to use a non-standard variation of set theory where the axiom of foundation doesn't hold we can no longer use recursion. For example the Axiom of Extension says that two sets are equal if their elements are equal. This is a recursive definition, but it's useless in the presence of a set a such that a∈a. At this point mathematicians need a principle of coinduction. And for more on that, I refer you to Vicious Circles.

Oh...time for a quick rant. Over the years I've seen a few people argue that there's something fundamentally wrong with the notion of the algorithm because it doesn't apply to the kind of open-ended loop we see in operating systems and interactive applications. Some have even gone further to suggest that somehow mathematics and computer science are fundamentally different because mathematics can't seek to describe these kinds of open-ended phenomena. As I've tried to show above, not only are there nice ways to look at open-ended computations, but from a mathematical perspective they are precisely dual, in the category theoretical sense, to terminating computations. It may be true that mathematicians sometimes spend more time with things, and computer scientists with cothings. But this really isn't such a big difference and the same langiage can be used to talk about both.

I learnt all this from a variety of sources including Turner's paper on Total Functional Programming. My original motivation for trying to understand this stuff came from this post on Ars Mathematica. Curiously, that post had a significant effect on me. It made me realise there was this entire co-universe out there that I knew nothing about, and I was forced to make a real effort with learning Haskell because Haskell was the programming language closest to the mathematical notation used to talk about this co-universe. (Thanks Walt, but sorry to everyone who wants to see less Haskell here!)

I didn't really get why anyone would want to distinguish between data and codata until I saw some clues in some slides by Altenkirch and some comments by Conor McBride. I'm still not sure anyone says it quite as clearly as this: distinguishing between data and codata means we can allow the coexistence of infinite lists, structural recursion and open-ended loops, without risk of causing bad behaviour.

Death>Dishonour =/=> ~Dishonour>~Death

That you rate death before dishonour should not imply that you rate honour above life. Or more precisely: if A and B are (not necessarily disjoint) outcomes, define A>B to mean that the expected utility of outcome A is higher than the expected utility of outcome B. Use the symbol '~' to mean the negation of an outcome. Then Death>Dishonour does not imply that ~Dishonour>~Death, or more generally, A>B =/=> ~B>~A.

This example from the book Subjective Probability: The Real Thing by Richard Jeffrey was pointed out to me by a friend of mine yesterday. Let's draw a little diagram of possible outcomes and their utility:

B\A	death	life
dishonour	1	*
honour	2	3

The '*' indicates an outcome with probability zero and the other outcomes are considered to be equally likely. A is the variable describing your state of liveness and B is the variable representing your honour.

We have:

E(utility|A=death) = 1.5

E(utility|B=dishonour) = 1



E(utility|B=honour) = 2.5

E(utility|A=life) = 3

With this set of utilities one does indeed rate death before dishonour but not honour before life.

This gets mixed reactions. Some people think that it's a nice counterintuitive result. Some people think it's so trivial they don't see why anyone would even mention it. If you're in the latter group, sorry for wasting your time. :-)

(I leave the Haskell program using a probability monad as an exercise...)

Monads from Algebra and the the Gray Code from Groups

There is a close association between algebraic structures and monads. I've mentioned this in passing before but I think it's interesting to work out some details. Among other things, it gives us a way of converting your favourite algebraic structures into domain specific languages embedded in Haskell.

First an aside. Haskell monads are in some sense only an approximation to mathematical monads. In fact, Haskell is only an approximation to mathematics. It's easy to define a Haskell function f, say, such that x == y but f x /= f y. Once you do such a thing you can no longer reason about Haskell functions with the assumption that == represents equality in the mathematical sense. (For an example, see Saizan's comment here.) So in the following I'm going to assume that we've limited ourselves to functions such that x == y implies that f x == f y, x < y and y < z imply x < z and so on. I'd love to see some way, in Haskell, to make explicit 'promises' about the properties of functions. (I guess that's what they call programming by contract.) But for the following we'll just assume it holds.

The first abstract algebraic structure that people study is usually the group, but I'm going to start with something simpler, the monoid. A monoid is a set M (called the underlying set of the monoid) with a binary operator, ·:M×M→M and identity element e such that e·x = x·e = x and x ·(y·z) = (x·y)·z for all x,y,z in M. By abuse of notation we'll often use M for the name of the monoid as well as the underlying set. Monoids are easy to come by and there are lots of examples. Any time you have a sequence of operations that can be applied, or not, in any order, to some system, you have a monoid. Sequences of operations form the underlying set of the monad and the binary operator means "do what's on the right" followed by "do what's on the left". For example the set of uniform scaling operations we can perform on a circle form a monoid. Write a scaling by a factor of x as s(x). If we double the size of a circle, and then triple it we have s(2) followed by s(3), written s(3)·s(2). Note that we have s(m)·s(n) = s(mn) and that s(1)·s(m)=s(m)·s(1)=s(m) so s(1) plays the role of the identity, e. Also note that for all m apart from 0 we have s(m)·s(1/m)=e so every scaling apart from s(0) has in inverse. Note however that this monoid is special because for all x and y in it, x·y=y·x. This doesn't always hold in monoids and when it does hold our monoid is said to be commutative.

If we have a type of algebraic structure then we can often form "free" versions of that structure. For monoids we get free monoids. There are many different ways of looking at "freeness" and I'm going to go through some of them informally:

1. Given any monoid there are equations satisfied by its elements. From the above example we have that x·y=y·x and also that x·e·x=x·x. But notice how there is a big difference between these equations. The former doesn't hold in all monads, but the latter does. A monoid is said to be free when the only equations that hold follow from equations that hold in every monoid. You can think of a free monoid as being a generic monoid. It has no special properties above and beyond being a monoid. Given a set S, the free monoid generated by S is the smallest free monoid containing the set S and with no equations relating the elements of S. Write this monoid as FS. For example, suppose S = {x,y}. Then we know that e, x and y are all in FS. We also know that x·x, x·y, y·x and y·y are all in S. Importantly we know that all 4 of these elements are distinct because we know there can be no equations connecting them beyond those that define a monoid. In fact, it's not hard to see that the elements of FS are essentially just the (possibly zero length) strings of x's and y's.
   
2. Given a set S, consider its elements to be 'labels' for unknown variables. In fact, just use the elements of S as variables. Then we can also consider the free monoid generated by S to be the set of "monoid expressions" in these unknowns. By "monoid expression" we just mean functions of these variables that can be written using the language of monoids, ie. the symbols e and ·. If S = {x,y} then examples of such expressions are e, x, y, x·y and so on. I hope it's not hard to see that this is simply another description of the same thing as in the previous paragraph.
   
3. Another way to think of free monoids requires a tiny bit more algebra. Define a monoid homomorphism from one monoid, M, to another, N, to be a function f:M→N on the underlying sets such that f(e)=e and f(x·y)=f(x)·f(y). A bijective monoid homomorphism is called an isomorphism. If there is an isomorphism between two monoids then in some sense they are the same monoid. Note that e and · are being overloaded here - we have two different monoids and these symbols have different meanings in each one. Now we can define the free monoid generated by S to be the monoid, M, such that (1) there is a function i:S→M (2) given any monoid N, and any function f:S→N, then f can be factored as f' o i where f' is a monoid homomorphism.
   
4. Very vaguely, the set S lives in the category of sets. But FS lives in the category of monoids. FS is the closest approximation to S in the category of monoids. The previous property gives a sense of what this means: anything that can be said about a function from S to a monoid N can be said in the language of monoids using a homomorphism from FS to N. Don't worry if this is too vague as I won't be using it below - but it may click in some people's minds.
   

The important thing here is that the operator F, ie. the "free monoid generated by" operator, forms a monad. It's tempting to think that monads got their name from monoids for this reason, but I don't think it is because just about any "free X generated by" operator forms a monad. Haskell monads give rise to DSLs via do-notation. So this means that algebraic structures give rise to monads, which may give rise to Haskell monads, and hence DSLs.

So now to explain why F forms a monad. Note that just about everything I say here works with other kinds of algebraic structures ranging from monoids through groups to vector spaces. Consider our set S={x,y} above. For convenience, let's drop the · and write the binary operator as multiplication in the usual way. Given an element of S we can easily get an element of FS. In fact, we have an embedding i:S→FS with i(x) = x. This is an abuse of notation. The element on the LHS is an element of S and the x on the RHS is a variable labeled with the symbol x, but we also write this as x because v_x, or something like that, would be a pain in the ass to keep writing.

Now think about the elements of F(FS). These are strings of elements of FS, ie. strings of monoid expressions in x and y. We could write a typical member as (xyx)(e)(xxy) where the xyx, e and xxy are each elements of FS and I'm using parentheses to make clear which bits they are. It's tempting to say that this is simply xyxxxy, but that would be wrong. The former is an element of F(FS) and the latter is an element of FS. It would be clearer that these weren't equal if we used the v notation to write v_xyxv_ev_xxy. But in this case the oscurity is to our advantage. Even though (xyx)(e)(xxy) doesn't equal xyxxxy, the notation strongly suggests that we define a function that maps the first to the second. In fact, we can define a function m:F(FS)→FS that essentially just erases parentheses.

So we have functions i:S→FS and m:F(FS)→FS. Looks suspiciously like a monad. In fact, i and m satisfy the monad laws (exercise!) and make F into a monad.

So what monad is this? It's probably not hard to guess: elements of FS are finite strings of elements of S. So it's essentially the List monad. Unfortunately, Haskell allows you to form infinite lists and so the correspondence isn't 100% precise. Nonetheless, it's good enough that if you hadn't already invented the List monad (as in Haskell monad), you probably would if you had considered the free monoid monad (as in mathematical monad). i is just the embedding \x -> [x] and m is concat, which essentially just erases brackets. So you can think of [x,y,z] has the Haskell way to write the monoid element xyz.

In fact, if you repeat the above discussion with other algebraic structures that I and others have discussed you'll get other familiar monads (modulo details like being able to create infinite lists). Here's a table (with apologies for bad formatting caused by blogger.com):

Algebraic Structure	DSL
Monoid	Combinatorial search
M-Set	"Write-only" side effects
Vector Space	Probability theory/Quantum Mechanics
Semirings	Two player game tree search
Modules over tropical semiring*	Min-sum algorithm
Group	???

(Note that the above table is approximate in the sense that sometimes you may need to restrict to instances of Eq and tweak the code slightly to make the Haskell monad behave like the mathematical one.)

Notice the gap in the "group" row. A group is a monoid in which every element has a left and right inverse. The monad is the "free group generated by" monad and I'll call it G. It doesn't correspond to any of the standard Haskell monads. So firstly - if free monoids give lists, what does the underlying datatype for free groups look like? Go back to our set S = {x,y}. GS contains e, x and y as well as all of the products of x and y that appear in FS. But additionally GS contains the inverses of x and y, x^-1 and y^-1. And of course you need all strings of x, y and their inverses. But do you need more? What about (xy)^-1y? Well we can expand out the inverse using (xy)^-1=y^-1x^-1. The net effect is that the free group contains precisely all (possibly empty) strings of x, y, x^-1 and y^-1, where substrings like xx^-1 and x^-1x are removed. We can model this in Haskell using the following type:

> import Control.Monad
> data Group a = G [Either a a] deriving Show

We're using Left x to represent x and Right x to represent x^-1. We implement >>= so that it uses the inverse of product formula above. It also ought to cancel out terms like xx^-1 but to do that requires that we restrict the monad to instances of Eq and if I do that I risk castigation from Saizan. So I'll leave out the cancellation to get the following Haskell monad:

> instance Monad Group where
>     return x = G [Left x]
>     G x >>= f = G $ concatMap g x where
>         g (Left y) = let G u = f y in u
>         g (Right z) = let G u = f z in reverse (map (either Right Left) u)

But what purpose might this serve? Try these:

> test1 = sequence $ replicate 4 (G [Left 0,Left 1])
> test2 = sequence $ replicate 4 (G [Left 0,Right 1])

The first does something almost identical to what you might expect from the ordinary list monad. But note the ordering in the second. This is the Gray code monad! By using Left and Right we can control the direction in which future combinatorial searches are carried out. Admittedly not the most useful monad, but it's curious that we do get something recognisable from this abstract nonsense. And maybe you spotted the "beautiful" Haskell implementation of the power set function on reddit recently. Here's a cute twist on that:

> powerset1 = filterM (const $ G [Left False,Left True])
> powerset2 = filterM (const $ G [Left False,Right True])

The first gives a result similar to the original, but the second lists the subsets in such a way that each element differs from the previous one by one element.

I wonder what other algebarically motivated monads are waiting to be discovered.

And I hope to write something about deriving comonads from coalgebras, as soon as I've read this. It looks really intriguing. Problem is, I'm having trouble making sense of it.

And sorry, this is a bit incoherent today. I only get a short period of time to write this stuff up and I blew most of mine this week in traffic jams. As usual, if the text makes no sense, you can be sure the code above works as I just tested it. But ask questions and complain about errors...

^* The tropical semiring is the proper name for the (R,min,+) semiring I talked about earlier.

Update: Saizan's doing a great job of keeping me honest. He noticed that I'd omitted a step. I was originally finding the inverse of a free group element by simply reversing the elements in the product. But I was failing to flip the polarities of each of the elements (so to speak) despite having spelled out in the text exactly what I needed to do. The code is now fixed.