-----------------------------------------------------------------------------
-- |
-- 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 :: forall a. Num a => [a] -> a -> a
polyeval []     a
_ = a
0
polyeval (a
p:[a]
ps) a
x = a
p forall a. Num a => a -> a -> a
+ a
x forall a. Num a => a -> a -> a
* forall a. Num a => [a] -> a -> a
polyeval [a]
ps a
x

-- | Add two polynomials

polyadd :: Num a => [a] -> [a] -> [a]
polyadd :: forall a. Num a => [a] -> [a] -> [a]
polyadd [] [a]
ys          = [a]
ys
polyadd [a]
xs []          = [a]
xs
polyadd (a
x:[a]
xs) (a
y:[a]
ys)  = (a
xforall a. Num a => a -> a -> a
+a
y) forall a. a -> [a] -> [a]
: forall a. Num a => [a] -> [a] -> [a]
polyadd [a]
xs [a]
ys

polyAddScalar :: Num a => a -> [a] -> [a]
polyAddScalar :: forall a. Num a => a -> [a] -> [a]
polyAddScalar a
c [] = [a
c]
polyAddScalar a
c (a
x:[a]
xs) = (a
cforall a. Num a => a -> a -> a
+a
x)forall a. a -> [a] -> [a]
:[a]
xs

-- | Subtract two polynomials

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

-- | Scale a polynomial

polyscale :: Num a => a -> [a] -> [a]
polyscale :: forall a. Num a => a -> [a] -> [a]
polyscale a
a [a]
x = forall a b. (a -> b) -> [a] -> [b]
map (a
aforall a. Num a => a -> a -> a
*) [a]
x

-- | Multiply two polynomials

polymult :: Num a => [a] -> [a] -> [a]
polymult :: forall a. Num a => [a] -> [a] -> [a]
polymult [a]
ys =
   forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr (\a
x [a]
acc -> forall a. Num a => [a] -> [a] -> [a]
polyadd (forall a. Num a => a -> [a] -> [a]
polyscale a
x [a]
ys) (a
0 forall a. a -> [a] -> [a]
: [a]
acc)) []

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

-- | Divide two polynomials

polydiv :: Fractional a => [a] -> [a] -> [a]
polydiv :: forall a. Fractional a => [a] -> [a] -> [a]
polydiv [a]
x0 [a]
y0 = forall a. [a] -> [a]
reverse forall a b. (a -> b) -> a -> b
$ forall a. Fractional a => [a] -> [a] -> [a]
polydiv' (forall a. [a] -> [a]
reverse [a]
x0) (forall a. [a] -> [a]
reverse [a]
y0)
    where polydiv' :: [a] -> [a] -> [a]
polydiv' (a
x:[a]
xs) [a]
y
             | forall (t :: * -> *) a. Foldable t => t a -> Int
length (a
xforall a. a -> [a] -> [a]
:[a]
xs) forall a. Ord a => a -> a -> Bool
< forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
y = []
             | Bool
otherwise = a
z forall a. a -> [a] -> [a]
: ([a] -> [a] -> [a]
polydiv' (forall a. [a] -> [a]
tail (forall a. Num a => [a] -> [a] -> [a]
polysub (a
xforall a. a -> [a] -> [a]
:[a]
xs) (forall a. Num a => [a] -> [a] -> [a]
polymult [a
z] [a]
y))) [a]
y)
                where z :: a
z = a
x forall a. Fractional a => a -> a -> a
/ forall a. [a] -> a
head [a]
y
          polydiv' [] [a]
_ = []

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

polymod :: Fractional a => [a] -> [a] -> [a]
polymod :: forall a. Fractional a => [a] -> [a] -> [a]
polymod [a]
x0 [a]
y0 = forall a. [a] -> [a]
reverse forall a b. (a -> b) -> a -> b
$ forall a. Fractional a => [a] -> [a] -> [a]
polymod' (forall a. [a] -> [a]
reverse [a]
x0) (forall a. [a] -> [a]
reverse [a]
y0)
    where polymod' :: [a] -> [a] -> [a]
polymod' (a
x:[a]
xs) [a]
y
             | forall (t :: * -> *) a. Foldable t => t a -> Int
length (a
xforall a. a -> [a] -> [a]
:[a]
xs) forall a. Ord a => a -> a -> Bool
< forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
y = (a
xforall a. a -> [a] -> [a]
:[a]
xs)
             | Bool
otherwise = [a] -> [a] -> [a]
polymod' (forall a. [a] -> [a]
tail (forall a. Num a => [a] -> [a] -> [a]
polysub (a
xforall a. a -> [a] -> [a]
:[a]
xs) (forall a. Num a => [a] -> [a] -> [a]
polymult [a
z] [a]
y))) [a]
y
                where z :: a
z = a
x forall a. Fractional a => a -> a -> a
/ forall a. [a] -> a
head [a]
y
          polymod' [] [a]
_ = []

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

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

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

polysubst :: Num a => [a] -> [a] -> [a]
polysubst :: forall a. Num a => [a] -> [a] -> [a]
polysubst [a]
ws = forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr (\a
x -> forall a. Num a => a -> [a] -> [a]
polyAddScalar a
x forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. Num a => [a] -> [a] -> [a]
polymult [a]
ws) []

polysubstAlt :: Num a => [a] -> [a] -> [a]
polysubstAlt :: forall a. Num a => [a] -> [a] -> [a]
polysubstAlt [a]
w0 [a]
x0 = forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr forall a. Num a => [a] -> [a] -> [a]
polyadd [a
0] (forall {a}. Num a => Int -> [a] -> [a] -> [[a]]
polysubst' Int
0 [a]
w0 [a]
x0)
    where polysubst' :: Int -> [a] -> [a] -> [[a]]
polysubst' Int
_ [a]
_ []     = []
          polysubst' Int
n [a]
w (a
x:[a]
xs) = forall a. Num a => a -> [a] -> [a]
polyscale a
x (forall a b. (Num a, Integral b) => [a] -> b -> [a]
polypow [a]
w (Int
n::Int)) forall a. a -> [a] -> [a]
: Int -> [a] -> [a] -> [[a]]
polysubst' (Int
nforall a. Num a => a -> a -> a
+Int
1) [a]
w [a]
xs

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

-- | Polynomial derivative

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

-- | Polynomial integration

polyinteg :: Fractional a => [a] -> a -> [a]
polyinteg :: forall a. Fractional a => [a] -> a -> [a]
polyinteg [a]
x0 a
c = a
c forall a. a -> [a] -> [a]
: forall {t}. Fractional t => t -> [t] -> [t]
polyinteg' a
1 [a]
x0
    where polyinteg' :: t -> [t] -> [t]
polyinteg' t
_ []     = []
          polyinteg' t
n (t
x:[t]
xs) = t
x forall a. Fractional a => a -> a -> a
/ t
n forall a. a -> [a] -> [a]
: t -> [t] -> [t]
polyinteg' (t
nforall a. Num a => a -> a -> a
+t
1) [t]
xs

-- | Convert roots to a polynomial

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