module Math.Polynomial
( Endianness(..)
, Poly, poly, polyDegree
, polyCoeffs, polyIsZero, polyIsOne
, zero, one, constPoly, x
, scalePoly, negatePoly
, composePoly
, addPoly, sumPolys, multPoly, powPoly
, quotRemPoly, quotPoly, remPoly
, evalPoly, evalPolyDeriv, evalPolyDerivs
, contractPoly
, monicPoly
, gcdPoly, separateRoots
, polyDeriv, polyDerivs, polyIntegral
) where
import Math.Polynomial.Type
import Math.Polynomial.Pretty ()
import Data.List
import Data.List.ZipSum
one :: (Num a, Eq a) => Poly a
one = constPoly 1
x :: (Num a, Eq a) => Poly a
x = polyN 2 LE [0,1]
constPoly :: (Num a, Eq a) => a -> Poly a
constPoly x = polyN 1 LE [x]
scalePoly :: (Num a, Eq a) => a -> Poly a -> Poly a
scalePoly 0 _ = zero
scalePoly s p = mapPoly (s*) p
negatePoly :: Num a => Poly a -> Poly a
negatePoly = mapPoly negate
addPoly :: (Num a, Eq a) => Poly a -> Poly a -> Poly a
addPoly p@(polyCoeffs LE -> a) q@(polyCoeffs LE -> b) = polyN n LE (zipSum a b)
where n = max (rawPolyLength p) (rawPolyLength q)
sumPolys :: (Num a, Eq a) => [Poly a] -> Poly a
sumPolys [] = zero
sumPolys ps = poly LE (foldl1 zipSum (map (polyCoeffs LE) ps))
multPoly :: (Num a, Eq a) => Poly a -> Poly a -> Poly a
multPoly p@(polyCoeffs LE -> xs) q@(polyCoeffs LE -> ys) = polyN n LE (multPolyLE xs ys)
where n = 1 + rawPolyDegree p + rawPolyDegree q
multPolyLE :: (Num a, Eq a) => [a] -> [a] -> [a]
multPolyLE _ [] = []
multPolyLE xs (y:ys) = foldr mul [] xs
where
mul 0 bs = 0 : bs
mul x bs = (x*y) : zipSum (map (x*) ys) bs
powPoly :: (Num a, Eq a, Integral b) => Poly a -> b -> Poly a
powPoly _ 0 = one
powPoly p 1 = p
powPoly p n
| n < 0 = error "powPoly: negative exponent"
| odd n = p `multPoly` powPoly p (n1)
| otherwise = (\x -> multPoly x x) (powPoly p (n`div`2))
quotRemPoly :: (Fractional a, Eq a) => Poly a -> Poly a -> (Poly a, Poly a)
quotRemPoly _ b | polyIsZero b = error "quotRemPoly: divide by zero"
quotRemPoly p@(polyCoeffs BE -> u) q@(polyCoeffs BE -> v)
= go [] u (polyDegree p polyDegree q)
where
v0 | null v = 0
| otherwise = head v
go q u n
| null u || n < 0 = (poly LE q, poly BE u)
| otherwise = go (q0:q) u' (n1)
where
q0 = head u / v0
u' = tail (zipSum u (map (negate q0 *) v))
quotPoly :: (Fractional a, Eq a) => Poly a -> Poly a -> Poly a
quotPoly u v
| polyIsZero v = error "quotPoly: divide by zero"
| otherwise = fst (quotRemPoly u v)
remPoly :: (Fractional a, Eq a) => Poly a -> Poly a -> Poly a
remPoly _ b | polyIsZero b = error "remPoly: divide by zero"
remPoly (polyCoeffs BE -> u) (polyCoeffs BE -> v)
= go u (length u length v)
where
v0 | null v = 0
| otherwise = head v
go u n
| null u || n < 0 = poly BE u
| otherwise = go u' (n1)
where
q0 = head u / v0
u' = tail (zipSum u (map (negate q0 *) v))
composePoly :: (Num a, Eq a) => Poly a -> Poly a -> Poly a
composePoly (polyCoeffs LE -> cs) (polyCoeffs LE -> ds) = poly LE (foldr mul [] cs)
where
mul c acc = addScalarLE c (multPolyLE acc ds)
addScalarLE :: (Num a, Eq a) => a -> [a] -> [a]
addScalarLE 0 bs = bs
addScalarLE a [] = [a]
addScalarLE a (b:bs) = (a + b) : bs
evalPoly :: (Num a, Eq a) => Poly a -> a -> a
evalPoly (polyCoeffs LE -> cs) 0
| null cs = 0
| otherwise = head cs
evalPoly (polyCoeffs LE -> cs) x = foldr mul 0 cs
where
mul c acc = c + acc * x
evalPolyDeriv :: (Num a, Eq a) => Poly a -> a -> (a,a)
evalPolyDeriv (polyCoeffs LE -> cs) x = foldr mul (0,0) cs
where
mul c (p, dp) = (p * x + c, dp * x + p)
evalPolyDerivs :: (Num a, Eq a) => Poly a -> a -> [a]
evalPolyDerivs (polyCoeffs LE -> cs) x = trunc . zipWith (*) factorials $ foldr mul [] cs
where
trunc list = zipWith const list cs
factorials = scanl (*) 1 (iterate (+1) 1)
mul c pds@(p:pd) = (p * x + c) : map (x *) pd `zipSum` pds
mul c [] = [c]
contractPoly :: (Num a, Eq a) => Poly a -> a -> (Poly a, a)
contractPoly p@(polyCoeffs LE -> cs) a = (polyN n LE q, r)
where
n = rawPolyLength p
cut remainder swap = (swap + remainder * a, remainder)
(r,q) = mapAccumR cut 0 cs
gcdPoly :: (Fractional a, Eq a) => Poly a -> Poly a -> Poly a
gcdPoly a b
| polyIsZero b = if polyIsZero a
then error "gcdPoly: gcdPoly zero zero is undefined"
else monicPoly a
| otherwise = gcdPoly b (a `remPoly` b)
monicPoly :: (Fractional a, Eq a) => Poly a -> Poly a
monicPoly p = case polyCoeffs BE p of
[] -> polyN n BE []
(c:cs) -> polyN n BE (1:map (/c) cs)
where n = rawPolyLength p
polyDeriv :: (Num a, Eq a) => Poly a -> Poly a
polyDeriv p@(polyCoeffs LE -> cs) = polyN (rawPolyDegree p) LE
[ c * n
| c <- drop 1 cs
| n <- iterate (1+) 1
]
polyDerivs :: (Num a, Eq a) => Poly a -> [Poly a]
polyDerivs p = take (1 + polyDegree p) (iterate polyDeriv p)
polyIntegral :: (Fractional a, Eq a) => Poly a -> Poly a
polyIntegral p@(polyCoeffs LE -> cs) = polyN (1 + rawPolyLength p) LE $ 0 :
[ c / n
| c <- cs
| n <- iterate (1+) 1
]
separateRoots :: (Fractional a, Eq a) => Poly a -> [Poly a]
separateRoots p
| polyIsZero q = error "separateRoots: zero polynomial"
| polyIsOne q = [p]
| otherwise = p `quotPoly` q : separateRoots q
where
q = gcdPoly p (polyDeriv p)