-----------------------------------------------------------------------------
-- |
-- Module : Polynomial.Basic
-- Copyright : (c) Matthew Donadio 2002
-- License : GPL
--
-- Maintainer : m.p.donadio@ieee.org
-- 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
-- | Add two polynomials
polyadd :: Num a => [a] -> [a] -> [a]
polyadd [] ys = ys
polyadd xs [] = xs
polyadd (x:xs) (y:ys) = (x+y) : polyadd xs ys
polyAddScalar :: Num a => a -> [a] -> [a]
polyAddScalar c [] = [c]
polyAddScalar c (x:xs) = (c+x):xs
-- | 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 [0] (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 [] = [1]
roots2poly (r:rs) = polymult [-r, 1] (roots2poly rs)