Alright, it’s been almost a month since I’ve last posted, and I wanted to post weekly. Perhaps bi-weekly would be a better goal, as long as I really can stick to it.

So, a few months ago, while working on a combinatorial problem for Project Euler, I was re-implementing polynomial arithmetic, and made an interesting observation. I was using a sparse list representation, with each element of the list being `(M c k)`

, representing the monomial `cx`

.^{k}

> data Mono = M !Integer !Integer deriving (Eq, Ord, Show)

A polynomial is represented by a list of monomials, which is assumed to be sorted in ascending order of the power `k`

above, and every power is assumed to be unique in our representation. Thus we can easily define addition and subtraction using standard techniques for dealing with ordered lists:

> mergeBy :: (Integer -> Integer -> Integer) -> [Mono] -> [Mono] -> [Mono] > mergeBy f xs ys = loop xs ys > where > loop [] ys = ys > loop xs [] = xs > loop (M c0 k0:xs) (M c1 k1:ys) > = case compare k0 k1 of > LT -> M c0 k0 : loop xs (M c1 k1:ys) > EQ -> let c' = f c0 c1 > in if c' == 0 > then loop xs ys > else M c' k0 : loop xs ys > GT -> M c1 k1 : loop (M c0 k0:xs) ys > add = mergeBy (+) > sub = mergeBy (-)

It’s also very straightforward to multiply two monomials, and to perform a “scalar” multiplication of a single monomial with a polynomial. Note that the degree of each term is changed by a constant amount, so we don’t have to worry about re-ordering the list, and that the first case of `smul`

is a completely optional minor optimization.

> mmul (M c0 k0) (M c1 k1) = M (c0 * c1) (k0 + k1) > smul (M 1 0) ys = ys > smul x ys = map (mmul x) ys

Taking into account the possibility of polynomials with an infinite number of terms, it’s now fairly easy to implement a polynomial multiplication program in terms of `add`

and `mmul`

:

> pmul _ [] = [] > pmul [] (_:_) = [] > pmul (x:xs) (y:ys) = mmul x y : add (smul x ys) (pmul xs (y:ys))

This is based around the algebraic identity `(x + xs) * (y + ys) == x * y + x * ys + xs * (y + ys)`

. It works on infinite polynomials because `x`

and `y`

are the terms of least degree to be multiplied, and every other product will have a higher degree. Unfortunately, this code turns out to be hopelessly inefficient.

The second thing I tried was simply changing the order in which arguments are passed to the recursive call to `pmul`

:

> pmul2 _ [] = [] > pmul2 [] (_:_) = [] > pmul2 (x:xs) (y:ys) = mmul x y : add (smul x ys) (pmul2 (y:ys) xs)

Empirically, this is asymptotically faster than the first attempt. I do not know why. Any insight would be fun to read about!