```-----------------------------------------------------------------------------
-- |
-- Module      :  Polynomial.Basic
--
-- Stability   :  experimental
-- Portability :  portable
--
-- Simple module for handling polynomials.
--
-----------------------------------------------------------------------------

-- TODO: We should really create a datatype for polynomials...

-- TODO: Should polydiv return the quotient and the remainder as a tuple?

module Polynomial.Basic where

-- * Types

-- | Polynomials are lists of numbers:
-- [ a0, a1, ... , an ] == an*x^n + ... + a1*x + a0
-- and negative exponents are currently verboten.

-- * Functions

-- | Evaluate a polynomial using Horner's method.

polyeval :: Num a => [a] -> a -> a
polyeval []     _ = 0
polyeval (p:ps) x = p + x * polyeval ps x

polyadd :: Num a => [a] -> [a] -> [a]

polyAddScalar :: Num a => a -> [a] -> [a]

-- | Subtract two polynomials

polysub :: Num a => [a] -> [a] -> [a]
polysub [] ys          = map negate ys
polysub xs []          = xs
polysub (x:xs) (y:ys)  = (x-y) : polysub xs ys

-- | Scale a polynomial

polyscale :: Num a => a -> [a] -> [a]
polyscale a x = map (a*) x

-- | Multiply two polynomials

polymult :: Num a => [a] -> [a] -> [a]
polymult ys =
foldr (\x acc -> polyadd (polyscale x ys) (0 : acc)) []

polymultAlt :: Num a => [a] -> [a] -> [a]
polymultAlt _  []     = []
polymultAlt ys (x:xs) = polyadd (polyscale x ys) (0 : polymult ys xs)

-- | Divide two polynomials

polydiv :: Fractional a => [a] -> [a] -> [a]
polydiv x0 y0 = reverse \$ polydiv' (reverse x0) (reverse y0)
where polydiv' (x:xs) y
| length (x:xs) < length y = []
| otherwise = z : (polydiv' (tail (polysub (x:xs) (polymult [z] y))) y)
where z = x / head y
polydiv' [] _ = []

-- | Modulus of two polynomials (remainder of division)

polymod :: Fractional a => [a] -> [a] -> [a]
polymod x0 y0 = reverse \$ polymod' (reverse x0) (reverse y0)
where polymod' (x:xs) y
| length (x:xs) < length y = (x:xs)
| otherwise = polymod' (tail (polysub (x:xs) (polymult [z] y))) y
where z = x / head y
polymod' [] _ = []

-- | Raise a polynomial to a non-negative integer power

polypow :: (Num a, Integral b) => [a] -> b -> [a]
polypow _ 0 = [ 1 ]
polypow x 1 = x
polypow x n | even n    = xSqr
| otherwise = polymult x xSqr
where xSqr = polymult x2 x2
x2   = polypow x (n `div` 2)

-- | Polynomial substitution y(n) = x(w(n))

polysubst :: Num a => [a] -> [a] -> [a]
polysubst ws = foldr (\x -> polyAddScalar x . polymult ws) []

polysubstAlt :: Num a => [a] -> [a] -> [a]
polysubstAlt w0 x0 = foldr polyadd  (polysubst' 0 w0 x0)
where polysubst' _ _ []     = []
polysubst' n w (x:xs) = polyscale x (polypow w (n::Int)) : polysubst' (n+1) w xs

{- |
Polynomial substitution @y(n) = x(w(n))@
where the coefficients of @x@ are also polynomials.
-}
polyPolySubst :: Num a => [a] -> [[a]] -> [a]
polyPolySubst ws = foldr (\x -> polyadd x . polymult ws) []

-- | Polynomial derivative

polyderiv :: Num a => [a] -> [a]
polyderiv [] = []
polyderiv (_:xs0) = polyderiv' 1 xs0
where polyderiv' _ []     = []
polyderiv' n (x:xs) = n * x : polyderiv' (n+1) xs

-- | Polynomial integration

polyinteg :: Fractional a => [a] -> a -> [a]
polyinteg x0 c = c : polyinteg' 1 x0
where polyinteg' _ []     = []
polyinteg' n (x:xs) = x / n : polyinteg' (n+1) xs

-- | Convert roots to a polynomial

roots2poly :: Num a => [a] -> [a]
roots2poly []     = 
roots2poly (r:rs) = polymult [-r, 1] (roots2poly rs)
```