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!
