module Math.Polynomial.VectorSpace
( Endianness(..)
, Poly, poly, polyDegree
, vPolyCoeffs, polyIsZero, polyIsOne
, zero, one, constPoly, x
, scalePoly, negatePoly
, composePolyWith
, addPoly, sumPolys, multPolyWith, powPolyWith
, quotRemPolyWith, quotPolyWith, remPolyWith
, evalPoly, evalPolyDeriv, evalPolyDerivs
, contractPoly
, monicPolyWith
, gcdPolyWith
, polyDeriv, polyDerivs, polyIntegral
) where
import Math.Polynomial.Type hiding (poly, polyDegree, polyIsZero)
import Math.Polynomial.Pretty ()
import Data.List
import Data.List.ZipSum
import Data.VectorSpace
vPolyN :: (Eq a, AdditiveGroup a) => Int -> Endianness -> [a] -> Poly a
vPolyN n e = vTrim . rawListPolyN n e
poly :: (Eq a, AdditiveGroup a) => Endianness -> [a] -> Poly a
poly e = vTrim . rawListPoly e
polyDegree :: (Eq a, AdditiveGroup a) => Poly a -> Int
polyDegree p = rawPolyDegree (vTrim p)
polyIsZero :: (Eq a, AdditiveGroup a) => Poly a -> Bool
polyIsZero = null . rawPolyCoeffs . vTrim
one :: (Num a, Eq a) => Poly a
one = polyN 1 LE [1]
x :: (Num a, Eq a) => Poly a
x = polyN 2 LE [0,1]
constPoly :: (Eq a, AdditiveGroup a) => a -> Poly a
constPoly x = vPolyN 1 LE [x]
scalePoly :: (Eq a, VectorSpace a, AdditiveGroup (Scalar a), Eq (Scalar a))
=> Scalar a -> Poly a -> Poly a
scalePoly = (*^)
negatePoly :: (AdditiveGroup a, Eq a) => Poly a -> Poly a
negatePoly = vTrim . rawMapPoly negateV
addPoly :: (AdditiveGroup a, Eq a) => Poly a -> Poly a -> Poly a
addPoly p@(vPolyCoeffs LE -> a) q@(vPolyCoeffs LE -> b) = vPolyN n LE (zipSumV a b)
where n = max (rawPolyLength p) (rawPolyLength q)
sumPolys :: (AdditiveGroup a, Eq a) => [Poly a] -> Poly a
sumPolys [] = zero
sumPolys ps = poly LE (foldl1 zipSumV (map (vPolyCoeffs LE) ps))
multPolyWith :: (AdditiveGroup a, Eq a) => (a -> a -> a) -> Poly a -> Poly a -> Poly a
multPolyWith multiplyV p@(vPolyCoeffs LE -> xs) q@(vPolyCoeffs LE -> ys) = vPolyN n LE (multPolyWithLE multiplyV xs ys)
where n = 1 + rawPolyDegree p + rawPolyDegree q
multPolyWithLE :: (AdditiveGroup a, Eq a) => (a -> a -> a) -> [a] -> [a] -> [a]
multPolyWithLE _ _ [] = []
multPolyWithLE multiplyV xs (y:ys) = foldr mul [] xs
where
mul x bs
| x == zeroV = zeroV : bs
| otherwise = (multiplyV x y) : zipSumV (map (multiplyV x) ys) bs
powPolyWith :: (AdditiveGroup a, Eq a, Integral b) => a -> (a -> a -> a) -> Poly a -> b -> Poly a
powPolyWith one multiplyV p n
| n < 0 = error "powPolyWith: negative exponent"
| otherwise = powPoly p n
where
multPoly = multPolyWith multiplyV
powPoly p 0 = constPoly one
powPoly p 1 = p
powPoly p n
| odd n = p `multPoly` powPoly p (n1)
| otherwise = (\x -> multPoly x x) (powPoly p (n`div`2))
quotRemPolyWith :: (AdditiveGroup a, Eq a) => (a -> a -> a) -> (a -> a -> a) -> Poly a -> Poly a -> (Poly a, Poly a)
quotRemPolyWith _ _ _ b | polyIsZero b = error "quotRemPoly: divide by zero"
quotRemPolyWith multiplyV divideV p@(vPolyCoeffs BE -> u) q@(vPolyCoeffs BE -> v)
= go [] u (polyDegree p polyDegree q)
where
v0 | null v = zeroV
| 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 = divideV (head u) v0
u' = tail (zipSumV u (map (multiplyV (negateV q0)) v))
quotPolyWith :: (AdditiveGroup a, Eq a) => (a -> a -> a) -> (a -> a -> a) -> Poly a -> Poly a -> Poly a
quotPolyWith multiplyV divideV u v
| polyIsZero v = error "quotPoly: divide by zero"
| otherwise = fst (quotRemPolyWith multiplyV divideV u v)
remPolyWith :: (AdditiveGroup a, Eq a) => (a -> a -> a) -> (a -> a -> a) -> Poly a -> Poly a -> Poly a
remPolyWith _ _ _ b | polyIsZero b = error "remPoly: divide by zero"
remPolyWith multiplyV divideV (vPolyCoeffs BE -> u) (vPolyCoeffs BE -> v)
= go u (length u length v)
where
v0 | null v = zeroV
| otherwise = head v
go u n
| null u || n < 0 = poly BE u
| otherwise = go u' (n1)
where
q0 = divideV (head u) v0
u' = tail (zipSumV u (map (multiplyV (negateV q0)) v))
composePolyWith :: (AdditiveGroup a, Eq a) => (a -> a -> a) -> Poly a -> Poly a -> Poly a
composePolyWith multiplyV (vPolyCoeffs LE -> cs) (vPolyCoeffs LE -> ds) = poly LE (foldr mul [] cs)
where
mul c acc = addScalarLE c (multPolyWithLE multiplyV acc ds)
addScalarLE :: (AdditiveGroup a, Eq a) => a -> [a] -> [a]
addScalarLE a bs | a == zeroV = bs
addScalarLE a [] = [a]
addScalarLE a (b:bs) = (a ^+^ b) : bs
evalPoly :: (VectorSpace a, Eq a, AdditiveGroup (Scalar a), Eq (Scalar a)) => Poly a -> Scalar a -> a
evalPoly (vPolyCoeffs LE -> cs) x
| x == zeroV =
if null cs
then zeroV
else head cs
| otherwise = foldr mul zeroV cs
where
mul c acc = c ^+^ acc ^* x
evalPolyDeriv :: (VectorSpace a, Eq a) => Poly a -> Scalar a -> (a,a)
evalPolyDeriv (vPolyCoeffs LE -> cs) x = foldr mul (zeroV, zeroV) cs
where
mul c (p, dp) = ((x *^ p) ^+^ c, (x *^ dp) ^+^ p)
evalPolyDerivs :: (VectorSpace a, Eq a, Num (Scalar a)) => Poly a -> Scalar a -> [a]
evalPolyDerivs (vPolyCoeffs 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) = (x *^ p ^+^ c) : map (x *^) pd `zipSumV` pds
mul c [] = [c]
contractPoly :: (VectorSpace a, Eq a) => Poly a -> Scalar a -> (Poly a, a)
contractPoly p@(vPolyCoeffs LE -> cs) a = (vPolyN n LE q, r)
where
n = rawPolyLength p
cut remainder swap = (swap ^+^ (a *^ remainder), remainder)
(r,q) = mapAccumR cut zeroV cs
gcdPolyWith :: (AdditiveGroup a, Eq a) => a -> (a -> a -> a) -> (a -> a -> a) -> Poly a -> Poly a -> Poly a
gcdPolyWith oneV multiplyV divideV a b
| polyIsZero b = if polyIsZero a
then error "gcdPolyWith: gcdPoly zero zero is undefined"
else monicPolyWith oneV divideV a
| otherwise = gcdPolyWith oneV multiplyV divideV b (a `remPoly` b)
where remPoly = remPolyWith multiplyV divideV
monicPolyWith :: (AdditiveGroup a, Eq a) => a -> (a -> a -> a) -> Poly a -> Poly a
monicPolyWith oneV divideV p = case vPolyCoeffs BE p of
[] -> vPolyN n BE []
(c:cs) -> vPolyN n BE (oneV : map (`divideV` c) cs)
where n = rawPolyLength p
polyDeriv :: (VectorSpace a, Eq a, Num (Scalar a)) => Poly a -> Poly a
polyDeriv p@(vPolyCoeffs LE -> cs) = vPolyN (rawPolyDegree p) LE
[ n *^ c
| c <- drop 1 cs
| n <- iterate (1+) 1
]
polyDerivs :: (VectorSpace a, Eq a, Num (Scalar a)) => Poly a -> [Poly a]
polyDerivs p = take (1 + polyDegree p) (iterate polyDeriv p)
polyIntegral :: (VectorSpace a, Eq a, Fractional (Scalar a)) => Poly a -> Poly a
polyIntegral p@(vPolyCoeffs LE -> cs) = vPolyN (1 + rawPolyLength p) LE $ zeroV :
[ c ^/ n
| c <- cs
| n <- iterate (1+) 1
]