Melding Monads

2009 April 7

Polynomial multiplication

Filed under: algorithms, corecursion, haskell — lpsmith @ 6:27 am

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 cxk.

> 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!

Blog at