



> import Ratio





z

count

count 1 z = 1

count n z = 0





> sample a = map (flip count a) [0..10]

> ex0 :: Fractional a => [a]

> ex0 = z





z

sample ex0

list

list a

a





> ex1 :: Fractional a => [a]

> ex1 = list z





count n (list z)

pair a

a's





> ex2 :: Fractional a => [a]

> ex2 = pair (list z)









> ex3 :: Fractional a => [a]

> ex3 = list (pair z)









> ex4 :: Fractional a => [a]

> ex4 = pair (nonEmpty (list z))





//





> ex5 :: Fractional a => [a]

> ex5 = pair ((list z) // (>=3))





ring





> ex6 :: Fractional a => [a]

> ex6 = ring z





sample ex6

set

set a

a

count n (set z)





> ex7 :: Fractional a => [a]

> ex7 = set (nonEmpty (ring z))





sample

ex7





> ex8 :: Fractional a => [a]

> ex8 = set (nonEmpty (set z))





oneOf

oneOf a b

a

b

count 3 ex9





> ex9 :: Fractional a => [a]

> ex9 = set (oneOf z (pair z))









> ex10 :: Fractional a => [a]

> ex10 = set ((list z) // (>=2))









> ex11 :: Fractional a => [a]

> ex11 = list (nonEmpty (ring z))









> ex12 :: Fractional a => [a]

> ex12 = set (nonEmpty (set z)) // even









> ex13 :: Fractional a => [a]

> ex13 = set (nonEmpty (set (nonEmpty (set z))))





0

1

2



a 0 +a 1 z+a 2 z2/2!+a 3 z3/3!+...



0

1

2





> (*!) _ 0 = 0

> (*!) a b = a*b



> (!*) 0 _ = 0

> (!*) a b = a*b









> (^+) a b = zipWith (+) a b

> (^-) a b = zipWith (-) a b









> ~(a:as) `convolve` (b:bs) = (a *! b):

> ((map (a !*) bs) ^+ (as `convolve` (b:bs)))









> compose (f:fs) (0:gs) = f:(gs `convolve` (compose fs (0:gs)))









> inverse (0:f:fs) = x where x = map (recip f *) (0:1:g)

> _:_:g = map negate (compose (0:0:fs) x)









> invert x = r where r = map (/x0) ((1:repeat 0) ^- (r `convolve` (0:xs)))

> x0:xs = x









> (^/) (0:a) (0:b) = a ^/ b

> (^/) a b = a `convolve` (invert b)





z





> z :: Fractional a => [a]

> z = 0:1:repeat 0









> d (_:x) = zipWith (*) (map fromInteger [1..]) x









> integrate x = 0 : zipWith (/) x (map fromInteger [1..])









> square x = x `convolve` x



> instance (Num r) => Num [r] where

> x+y = zipWith (+) x y

> x-y = zipWith (-) x y

> ~x*y = x `convolve` y

> fromInteger x = fromInteger x:repeat 0

> negate x = map negate x

> signum (x:_) = signum x:repeat 0

> abs (x:xs) = error "Can't form abs of a power series"



> instance (Fractional r) => Fractional [r] where

> x/y = x ^/ y

> fromRational x = fromRational x:repeat 0



> sqrt' x = 1:rs where rs = map (/2) (xs ^- (rs `convolve` (0:rs)))

> _:xs = x



> instance (Fractional r) => Floating [r] where

> sqrt (1:x) = sqrt' (1:x)

> sqrt _ = error "Can only find sqrt when leading term is 1"









> exp x = e where e = 1+integrate (e * d x)









> log x = integrate (d x/x)









> sin x = integrate ((cos x)*(d x))

> cos x = [1] ... negate (integrate ((sin x)*(d x)))









> asin x = integrate (d x/sqrt(1-x*x))

> atan x = integrate (d x/(1+x*x))

> acos x = error "Unable to form power series for acos"









> sinh x = integrate ((cosh x)*(d x))

> cosh x = [1] ... integrate ((sinh x)*(d x))



> asinh x = integrate (d x/sqrt(1+x*x))

> atanh x = integrate (d x/(1-x*x))

> acosh x = error "Unable to form power series for acosh"



> pi = error "There is no formal power series for pi"









> t :: (Num a) => ([a])

> t = 0:1:repeat 0



> t' :: [Rational]

> t' = t









> lead [] x = x

> lead (a:as) x = a : (lead as (tail x))



> a ... x = lead a x









> one = t'

> list x = 1/(1-x)

> set x = exp x

> ring x = -log(1-x)

> pair x = x*x

> oneOf a b = a+b









> necklace x = -log(1-x)/2+x/2+x*x/4





union a b





> union a b = a*b



> (//) :: Fractional a => [a] -> (Integer -> Bool) -> [a]

> (//) a c = zipWith (\a-> \b->(if (c a :: Bool) then b else 0)) [(0::Integer)..] a

> nonEmpty a = a // (/= 0)



> factorial 0 = 1

> factorial n = n*factorial (n-1)



> count n a = numerator ((a!!(fromInteger n)) * (factorial (fromInteger n)))









> tree x = p where p = [0] ... union (set p) x





sample (tree z)

graph





> graph = [2^((n*(n-1) `div` 2)) / product (map fromInteger [1..n]) | n







> connectedGraph = 1 + log graph





I thought it'd be fun to dig up my first ever Haskell project from about 2 years ago. I wrote a little about part of it before but I didn't write about what I did with it. Apologies for some of the formatting, I haven't yet got around to updating every line of this code. But it definitely still works. (As usual with old code, I had to fix it in places. How is it that old code ceases to work? Is there some kind of bit rot that sets in if you leave code unused for more than a couple of years?)And as usual, this post is written in literate Haskell so you can just run it interacively in ghci.The idea here was to define a little DSL in Haskell for counting the size of a variety of combinatorial objects. Rather than start with the theory, I'm going to look at some examples in action first.The idea is that I'll be considering different ways to structure the set {1,2,...,n} (with n=0 meaning the empty set). Let's start with the simplest example: the singleton set with one element. I'll denote this by the symbol. Consider how many ways we can structure the set {1,2,...,n} as a one element set. It's pretty obvious that you can only structure this set as a one element set if it has one element. So if we useto count the number of ways of doing things, we expectandfor n equal to anything else.In fact, let's defineI haven't told you what valuehas, just how it responds to the count function. You can now try evaluatingNow consider the next function in our DSL,gives the number of ways to form a list of's from {1,...,n}. A good example is simplyie. the number of ways to form a list of n http://sigfpe.blogspot.com/elements. For example, for n=3 we have [1,2,3],[1,3,2],[2,1,3],[2,3,1],[3,1,2],[3,2,1]. Obviouslyis n!.Many of the examples I'm going to look at are listed in Sloane's famous encyclopedia of integer sequences. n! is listed as A000142 Next considerwhich is the number of ways of forming an ordered pair of. For example, the pairs of lists from {1,2,3} are ([],[1,2,3]),([1],[2,3]),([2],[1,3]),([3],[1,2]),([1,2],[3]) and so on. With a little work you should count 24 of those. This is given byConversely we can consider lists of pairs. Clearly this is zero for n odd. For n=4 we get terms like [(1,2),(3,4)], [(1,3),(2,4)] and so on. Again there are 24. The code is:This DSL also supports constraints. For example we can redo the pair of lists example with the proviso that the lists are both non-empty:More generally theoperator can provide any numerical constraint. For example:Gives the ways of forming pairs of lists of length >=3. Clearly this is zero for n<6 and for n=6 we just get 6!At this point I need to point out that these examples all run petty quickly so we're not simply generating and counting.The next function is. A ring is like a list except that cyclic permutations are considered equal. So [1,2,3], [2,3,1] and [3,1,2] are the same ring. How many rings are there on {1,2,...,n}?If you evaluateyou may recognise that it's simply (n-1)!.The next useful function iscounts the number of sets of's you can form from {1,...,n}. Clearlyis 1 as there's only one set you can form from {1,...,n}. So here's a nice example:That's counting the ways of forming {1,...,n} into sets of non-empty rings. If you applyyou may spot that we have n! again. There's a really neat explanation why. Every permutation of {1,...,n} is a product of disjoint cycles. The order of the product doesn't matter so we have a set (rather than a list) of disjoint cycles. A cycle is obviously non-empty and cycles are just rings. Socounts permutations. And the number of permutations on {1,...,n} is n!.Here's an old classic ( A000110 ):This is the number of ways of partitioning {1,...,n} into sets of non-empty sets. For n=3 we get {{1,2,3}},{{1},{2,3}},{{2},{1,3}},{{3},{1,2}} and {{1},{2},{3}}. In other words, 5.Another useful function ismeans forming {1,...,n} into either anor aFor example we can count the ways to split a set into a set of singletons and pairs. For n=3 we have {1,(2,3)},{2,(1,3)},{3,(1,2)}, the previous three with the pair flipped and {1,2,3}, in other wordsis 7:By the way, that's Sloane's A047974 Now some more examples just for fun:Partitions into lists of size >=2. (Sloane's A052845 ):Partitions into lists of cycles (Sloane's A007840):Partitions into an even number of non-empty sets (Sloane's A020557):Partitions into sets of non-empty sets of non-empty sets (Sloane's A000258):Thirteen examples is enough for anyone. So how does it work? I haven't time to give a crash course on generating functions but I can at least sketch the idea. In particular, exponential generating functions. The idea is this: if I have a sequence a,a,a,... then its exponential generating function is:So we can represent any sequence of integers as a power series. For some sequences this series doesn't converge. But that's not a problem. We can still formally manipulate the series as long as we don't try to substitute an actual real number for z. So we call these formal power series. The amazing fact is this: the sorts of operations we might perform on power series when performing ordinary arithmetic (adding hem, multiplying them, finding their logarithms and so on) turn out to be exactly the kind of operations you would perform on them to solve combinatorial problems. So many combinatorial problems can be reduced to finding power series for various arithmetical combinations of other power series.So now we need some code to manipulate power series. This is such a popular Haskell example that I'll breeze through the code and give some links to other people's documentation. I represent power series simply as a list of coefficients.When computing with formal power series we frequently require products of the form (a+a*t+a*t+...)*0 We know the result is zero no matter what the ai are and by immediately returning zero we can eliminate bottomless recursions. We use the operator *! for this purpose. !* is the symmetric counterpart of *!We represent formal power series over an algebraic structure a as infinite list of type a.Addition and subtraction are straightforward pointwise applications of addition and subtraction in R.Multiplication itself is achieved through the operation of convolution:If f(z) and g(z) are two power series, this computes f(g(z)):Note how that function fails if the series for g starts with anything other than zero. A few of my operations will be like that. They correspond to situations where the operation doesn't make sense on a power series (at least no without some additional structure or assumptions).The (functional) inverse of a function:Reciprocal of power series:Division of power series is defined in terms of invert though there is probably a direct definition that is as simple. Note the special case handling of power series that both start with zero:Here comes the definition ofthat we've been using:z = [0,1,0,0,...] in other words the generator of formal power series.one = [1,0,0,...], ie. the multiplicative unit.d computes the derivative of formal power series usingd (t^n) = n*t^(n-1)Integration is simply the inverse of differentiation. We assume the integral has a zero leading term.Compute the square of a formal power series:We apply transcendental functions to formal power series by finding suitable integro-differential representations of the function. For exp consider:d/dt exp(f(t)) = f'(t)*exp(f(t))We can immediately writeexp(f(t)) = integral f'(t)*exp(f(t))This gives a recursive definition of the exponential of a power series.d/dt log(f(t)) = f'(t)/f(t) which gives the following definition:We give mutually recursive definitions for sine and cosine:The inverse trignometric functions are similarly straightforward except for the inverse cosine. The inverse cosine at zero is transcendental and so might not exist in R. On the other hand there is no power series for acos(1+x). So we can give no formal power series for this function.The hyperbolic trigonometric functions are similar to the ordinary trigonometric functions.t is the generator of the formal power series. t' is the generator for formal power series over the rationals.The ... operator is used to specify a formal power series with known leading terms. If l is a list of n elements then the first n elements of l ... f are given by the elements of l. The remaining terms are those of l.And now we can get back to combinatorics. The definitions of the functions I gave above are stunningly simple:How does that work? Just about any book on combinatorics will tell you. I quite like Generatingfunctionology though it does make a mountain out of a molehill in places.is the number of ways of splitting the set {1,...,n} into two so that one half is a combinatorial structure given by a and the other half is given by b. It's definition is amazingly simple:Now for some more examples. A tree is, by definition, a single parent node and a set of child trees. So we have the recursive definition:The [0] ... is to get the recursion started. Tryto count the number of trees.The last example is neat.is the exponential generating function of the number of ways of making a graph from {1,...,n}. ( A006125 ) It's not too hard to see what's going on: a graph on a set S is given by a set of edges. In other words it's a subset of the set of unordered pairs of elements of S. If S has size n, there are n(n-1)/2 such pairs.exp has a really amazing property. If f is the expnential power series corresponding to counting some structure, say F, on {1,...,n} then exp(f) is the power series of the structure consisting of the union of a bunch of non-zero sized F's. We can put this into practice: all graphs can be decomposed as the union of connected graphs. So the generating function of connected graphs is given by(The 1+ is because of the non-zero I just mentioned.)That was Sloane's A001187 By the way, this plays a big role in quantum field theory where we have to enumerate Feynman diagrams Anyway, if you didn't previously know about generating functions I hope this whets your appetite to read up on them some more. They are an almost magical part of mathematics that makes the solution of some very difficult seeming combinatorial problems almost trivial. They also have countless real world applications ranging from statistics to digital signal processing.