{-# LANGUAGE ParallelListComp, ViewPatterns, FlexibleContexts #-}
{-# OPTIONS_GHC -fno-warn-orphans #-}
-- TODO: update all haddock comments
-- |Same general interface as Math.Polynomial, but using AdditiveGroup, 
-- VectorSpace, etc., instead of Num where sensible.
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 ({- instance -})

import Data.List
import Data.List.ZipSum

import Data.VectorSpace

vPolyN :: (Eq a, AdditiveGroup a) => Int -> Endianness -> [a] -> Poly a
vPolyN :: Int -> Endianness -> [a] -> Poly a
vPolyN n :: Int
n e :: Endianness
e = Poly a -> Poly a
forall a. (Eq a, AdditiveGroup a) => Poly a -> Poly a
vTrim (Poly a -> Poly a) -> ([a] -> Poly a) -> [a] -> Poly a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> Endianness -> [a] -> Poly a
forall a. Int -> Endianness -> [a] -> Poly a
rawListPolyN Int
n Endianness
e

poly :: (Eq a, AdditiveGroup a) => Endianness -> [a] -> Poly a
poly :: Endianness -> [a] -> Poly a
poly e :: Endianness
e = Poly a -> Poly a
forall a. (Eq a, AdditiveGroup a) => Poly a -> Poly a
vTrim (Poly a -> Poly a) -> ([a] -> Poly a) -> [a] -> Poly a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Endianness -> [a] -> Poly a
forall a. Endianness -> [a] -> Poly a
rawListPoly Endianness
e

polyDegree :: (Eq a, AdditiveGroup a) => Poly a -> Int
polyDegree :: Poly a -> Int
polyDegree p :: Poly a
p = Poly a -> Int
forall a. Poly a -> Int
rawPolyDegree (Poly a -> Poly a
forall a. (Eq a, AdditiveGroup a) => Poly a -> Poly a
vTrim Poly a
p)

polyIsZero :: (Eq a, AdditiveGroup a) => Poly a -> Bool
polyIsZero :: Poly a -> Bool
polyIsZero = [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null ([a] -> Bool) -> (Poly a -> [a]) -> Poly a -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Poly a -> [a]
forall a. Poly a -> [a]
rawPolyCoeffs (Poly a -> [a]) -> (Poly a -> Poly a) -> Poly a -> [a]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Poly a -> Poly a
forall a. (Eq a, AdditiveGroup a) => Poly a -> Poly a
vTrim

-- |The polynomial \"1\"
one :: (Num a, Eq a) => Poly a
one :: Poly a
one = Int -> Endianness -> [a] -> Poly a
forall a. (Num a, Eq a) => Int -> Endianness -> [a] -> Poly a
polyN 1 Endianness
LE [1]

-- |The polynomial (in x) \"x\"
x :: (Num a, Eq a) => Poly a
x :: Poly a
x = Int -> Endianness -> [a] -> Poly a
forall a. (Num a, Eq a) => Int -> Endianness -> [a] -> Poly a
polyN 2 Endianness
LE [0,1]

-- |Given some constant 'k', construct the polynomial whose value is 
-- constantly 'k'.
constPoly :: (Eq a, AdditiveGroup a) => a -> Poly a
constPoly :: a -> Poly a
constPoly x :: a
x = Int -> Endianness -> [a] -> Poly a
forall a.
(Eq a, AdditiveGroup a) =>
Int -> Endianness -> [a] -> Poly a
vPolyN 1 Endianness
LE [a
x]

-- |Given some scalar 's' and a polynomial 'f', computes the polynomial 'g'
-- such that:
-- 
-- > evalPoly g x = s * evalPoly f x
scalePoly :: (Eq a, VectorSpace a, AdditiveGroup (Scalar a), Eq (Scalar a)) 
    => Scalar a -> Poly a -> Poly a
scalePoly :: Scalar a -> Poly a -> Poly a
scalePoly = Scalar a -> Poly a -> Poly a
forall v. VectorSpace v => Scalar v -> v -> v
(*^)

-- |Given some polynomial 'f', computes the polynomial 'g' such that:
-- 
-- > evalPoly g x = negate (evalPoly f x)
negatePoly :: (AdditiveGroup a, Eq a) => Poly a -> Poly a
negatePoly :: Poly a -> Poly a
negatePoly = Poly a -> Poly a
forall a. (Eq a, AdditiveGroup a) => Poly a -> Poly a
vTrim (Poly a -> Poly a) -> (Poly a -> Poly a) -> Poly a -> Poly a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (a -> a) -> Poly a -> Poly a
forall a. (a -> a) -> Poly a -> Poly a
rawMapPoly a -> a
forall v. AdditiveGroup v => v -> v
negateV

-- |Given polynomials 'f' and 'g', computes the polynomial 'h' such that:
-- 
-- > evalPoly h x = evalPoly f x + evalPoly g x
addPoly :: (AdditiveGroup a, Eq a) => Poly a -> Poly a -> Poly a
addPoly :: Poly a -> Poly a -> Poly a
addPoly p :: Poly a
p@(Endianness -> Poly a -> [a]
forall a. (Eq a, AdditiveGroup a) => Endianness -> Poly a -> [a]
vPolyCoeffs Endianness
LE ->  [a]
a) q :: Poly a
q@(Endianness -> Poly a -> [a]
forall a. (Eq a, AdditiveGroup a) => Endianness -> Poly a -> [a]
vPolyCoeffs Endianness
LE ->  [a]
b) = Int -> Endianness -> [a] -> Poly a
forall a.
(Eq a, AdditiveGroup a) =>
Int -> Endianness -> [a] -> Poly a
vPolyN Int
n Endianness
LE ([a] -> [a] -> [a]
forall t. AdditiveGroup t => [t] -> [t] -> [t]
zipSumV [a]
a [a]
b)
    where n :: Int
n = Int -> Int -> Int
forall a. Ord a => a -> a -> a
max (Poly a -> Int
forall a. Poly a -> Int
rawPolyLength Poly a
p) (Poly a -> Int
forall a. Poly a -> Int
rawPolyLength Poly a
q)

{-# RULES
  "sum Poly"    forall ps. foldl addPoly zero ps = sumPolys ps
  #-}
sumPolys :: (AdditiveGroup a, Eq a) => [Poly a] -> Poly a
sumPolys :: [Poly a] -> Poly a
sumPolys [] = Poly a
forall a. Poly a
zero
sumPolys ps :: [Poly a]
ps = Endianness -> [a] -> Poly a
forall a. (Eq a, AdditiveGroup a) => Endianness -> [a] -> Poly a
poly Endianness
LE (([a] -> [a] -> [a]) -> [[a]] -> [a]
forall (t :: * -> *) a. Foldable t => (a -> a -> a) -> t a -> a
foldl1 [a] -> [a] -> [a]
forall t. AdditiveGroup t => [t] -> [t] -> [t]
zipSumV ((Poly a -> [a]) -> [Poly a] -> [[a]]
forall a b. (a -> b) -> [a] -> [b]
map (Endianness -> Poly a -> [a]
forall a. (Eq a, AdditiveGroup a) => Endianness -> Poly a -> [a]
vPolyCoeffs Endianness
LE) [Poly a]
ps))

-- |Given polynomials 'f' and 'g', computes the polynomial 'h' such that:
-- 
-- > evalPoly h x = evalPoly f x * evalPoly g x
multPolyWith :: (AdditiveGroup a, Eq a) => (a -> a -> a) -> Poly a -> Poly a -> Poly a
multPolyWith :: (a -> a -> a) -> Poly a -> Poly a -> Poly a
multPolyWith multiplyV :: a -> a -> a
multiplyV p :: Poly a
p@(Endianness -> Poly a -> [a]
forall a. (Eq a, AdditiveGroup a) => Endianness -> Poly a -> [a]
vPolyCoeffs Endianness
LE -> [a]
xs) q :: Poly a
q@(Endianness -> Poly a -> [a]
forall a. (Eq a, AdditiveGroup a) => Endianness -> Poly a -> [a]
vPolyCoeffs Endianness
LE -> [a]
ys) = Int -> Endianness -> [a] -> Poly a
forall a.
(Eq a, AdditiveGroup a) =>
Int -> Endianness -> [a] -> Poly a
vPolyN Int
n Endianness
LE ((a -> a -> a) -> [a] -> [a] -> [a]
forall a.
(AdditiveGroup a, Eq a) =>
(a -> a -> a) -> [a] -> [a] -> [a]
multPolyWithLE a -> a -> a
multiplyV [a]
xs [a]
ys)
    where n :: Int
n = 1 Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Poly a -> Int
forall a. Poly a -> Int
rawPolyDegree Poly a
p Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Poly a -> Int
forall a. Poly a -> Int
rawPolyDegree Poly a
q

-- |(Internal): multiply polynomials in LE order.  O(length xs * length ys).
multPolyWithLE :: (AdditiveGroup a, Eq a) => (a -> a -> a) -> [a] -> [a] -> [a]
multPolyWithLE :: (a -> a -> a) -> [a] -> [a] -> [a]
multPolyWithLE _         _  []     = []
multPolyWithLE multiplyV :: a -> a -> a
multiplyV xs :: [a]
xs (y :: a
y:ys :: [a]
ys) = (a -> [a] -> [a]) -> [a] -> [a] -> [a]
forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr a -> [a] -> [a]
mul [] [a]
xs
    where
        mul :: a -> [a] -> [a]
mul x :: a
x bs :: [a]
bs
            | a
x a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
forall v. AdditiveGroup v => v
zeroV    = a
forall v. AdditiveGroup v => v
zeroV a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [a]
bs
            | Bool
otherwise     = (a -> a -> a
multiplyV a
x a
y) a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [a] -> [a] -> [a]
forall t. AdditiveGroup t => [t] -> [t] -> [t]
zipSumV ((a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (a -> a -> a
multiplyV a
x) [a]
ys) [a]
bs

-- |Given a polynomial 'f' and exponent 'n', computes the polynomial 'g'
-- such that:
-- 
-- > evalPoly g x = evalPoly f x ^ n
powPolyWith :: (AdditiveGroup a, Eq a, Integral b) => a -> (a -> a -> a) -> Poly a -> b -> Poly a
powPolyWith :: a -> (a -> a -> a) -> Poly a -> b -> Poly a
powPolyWith one :: a
one multiplyV :: a -> a -> a
multiplyV p :: Poly a
p n :: b
n
    | b
n b -> b -> Bool
forall a. Ord a => a -> a -> Bool
< 0     = [Char] -> Poly a
forall a. HasCallStack => [Char] -> a
error "powPolyWith: negative exponent"
    | Bool
otherwise = Poly a -> b -> Poly a
forall a. Integral a => Poly a -> a -> Poly a
powPoly Poly a
p b
n
    where
        multPoly :: Poly a -> Poly a -> Poly a
multPoly = (a -> a -> a) -> Poly a -> Poly a -> Poly a
forall a.
(AdditiveGroup a, Eq a) =>
(a -> a -> a) -> Poly a -> Poly a -> Poly a
multPolyWith a -> a -> a
multiplyV
        powPoly :: Poly a -> a -> Poly a
powPoly p :: Poly a
p 0 = a -> Poly a
forall a. (Eq a, AdditiveGroup a) => a -> Poly a
constPoly a
one
        powPoly p :: Poly a
p 1 = Poly a
p
        powPoly p :: Poly a
p n :: a
n 
            | a -> Bool
forall a. Integral a => a -> Bool
odd a
n     = Poly a
p Poly a -> Poly a -> Poly a
`multPoly` Poly a -> a -> Poly a
powPoly Poly a
p (a
na -> a -> a
forall a. Num a => a -> a -> a
-1)
            | Bool
otherwise = (\x :: Poly a
x -> Poly a -> Poly a -> Poly a
multPoly Poly a
x Poly a
x) (Poly a -> a -> Poly a
powPoly Poly a
p (a
na -> a -> a
forall a. Integral a => a -> a -> a
`div`2))

-- |Given polynomials @a@ and @b@, with @b@ not 'zero', computes polynomials
-- @q@ and @r@ such that:
-- 
-- > addPoly (multPoly q b) r == a
quotRemPolyWith :: (AdditiveGroup a, Eq a) => (a -> a -> a) -> (a -> a -> a) -> Poly a -> Poly a -> (Poly a, Poly a)
quotRemPolyWith :: (a -> a -> a)
-> (a -> a -> a) -> Poly a -> Poly a -> (Poly a, Poly a)
quotRemPolyWith _ _ _ b :: Poly a
b | Poly a -> Bool
forall a. (Eq a, AdditiveGroup a) => Poly a -> Bool
polyIsZero Poly a
b = [Char] -> (Poly a, Poly a)
forall a. HasCallStack => [Char] -> a
error "quotRemPoly: divide by zero"
quotRemPolyWith multiplyV :: a -> a -> a
multiplyV divideV :: a -> a -> a
divideV p :: Poly a
p@(Endianness -> Poly a -> [a]
forall a. (Eq a, AdditiveGroup a) => Endianness -> Poly a -> [a]
vPolyCoeffs Endianness
BE -> [a]
u) q :: Poly a
q@(Endianness -> Poly a -> [a]
forall a. (Eq a, AdditiveGroup a) => Endianness -> Poly a -> [a]
vPolyCoeffs Endianness
BE -> [a]
v)
    = [a] -> [a] -> Int -> (Poly a, Poly a)
forall t. (Ord t, Num t) => [a] -> [a] -> t -> (Poly a, Poly a)
go [] [a]
u (Poly a -> Int
forall a. (Eq a, AdditiveGroup a) => Poly a -> Int
polyDegree Poly a
p Int -> Int -> Int
forall a. Num a => a -> a -> a
- Poly a -> Int
forall a. (Eq a, AdditiveGroup a) => Poly a -> Int
polyDegree Poly a
q)
    where
        v0 :: a
v0  | [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [a]
v    = a
forall v. AdditiveGroup v => v
zeroV
            | Bool
otherwise = [a] -> a
forall a. [a] -> a
head [a]
v
        go :: [a] -> [a] -> t -> (Poly a, Poly a)
go q :: [a]
q u :: [a]
u n :: t
n
            | [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [a]
u Bool -> Bool -> Bool
|| t
n t -> t -> Bool
forall a. Ord a => a -> a -> Bool
< 0   = (Endianness -> [a] -> Poly a
forall a. (Eq a, AdditiveGroup a) => Endianness -> [a] -> Poly a
poly Endianness
LE [a]
q, Endianness -> [a] -> Poly a
forall a. (Eq a, AdditiveGroup a) => Endianness -> [a] -> Poly a
poly Endianness
BE [a]
u)
            | Bool
otherwise         = [a] -> [a] -> t -> (Poly a, Poly a)
go (a
q0a -> [a] -> [a]
forall a. a -> [a] -> [a]
:[a]
q) [a]
u' (t
nt -> t -> t
forall a. Num a => a -> a -> a
-1)
            where
                q0 :: a
q0 = a -> a -> a
divideV ([a] -> a
forall a. [a] -> a
head [a]
u) a
v0
                u' :: [a]
u' = [a] -> [a]
forall a. [a] -> [a]
tail ([a] -> [a] -> [a]
forall t. AdditiveGroup t => [t] -> [t] -> [t]
zipSumV [a]
u ((a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (a -> a -> a
multiplyV (a -> a
forall v. AdditiveGroup v => v -> v
negateV a
q0)) [a]
v))

quotPolyWith :: (AdditiveGroup a, Eq a) => (a -> a -> a) -> (a -> a -> a) -> Poly a -> Poly a -> Poly a
quotPolyWith :: (a -> a -> a) -> (a -> a -> a) -> Poly a -> Poly a -> Poly a
quotPolyWith multiplyV :: a -> a -> a
multiplyV divideV :: a -> a -> a
divideV u :: Poly a
u v :: Poly a
v
    | Poly a -> Bool
forall a. (Eq a, AdditiveGroup a) => Poly a -> Bool
polyIsZero Poly a
v  = [Char] -> Poly a
forall a. HasCallStack => [Char] -> a
error "quotPoly: divide by zero"
    | Bool
otherwise     = (Poly a, Poly a) -> Poly a
forall a b. (a, b) -> a
fst ((a -> a -> a)
-> (a -> a -> a) -> Poly a -> Poly a -> (Poly a, Poly a)
forall a.
(AdditiveGroup a, Eq a) =>
(a -> a -> a)
-> (a -> a -> a) -> Poly a -> Poly a -> (Poly a, Poly a)
quotRemPolyWith a -> a -> a
multiplyV a -> a -> a
divideV Poly a
u Poly a
v)
remPolyWith :: (AdditiveGroup a, Eq a) => (a -> a -> a) -> (a -> a -> a) -> Poly a -> Poly a -> Poly a
remPolyWith :: (a -> a -> a) -> (a -> a -> a) -> Poly a -> Poly a -> Poly a
remPolyWith _ _ _ b :: Poly a
b | Poly a -> Bool
forall a. (Eq a, AdditiveGroup a) => Poly a -> Bool
polyIsZero Poly a
b = [Char] -> Poly a
forall a. HasCallStack => [Char] -> a
error "remPoly: divide by zero"
remPolyWith multiplyV :: a -> a -> a
multiplyV divideV :: a -> a -> a
divideV (Endianness -> Poly a -> [a]
forall a. (Eq a, AdditiveGroup a) => Endianness -> Poly a -> [a]
vPolyCoeffs Endianness
BE -> [a]
u) (Endianness -> Poly a -> [a]
forall a. (Eq a, AdditiveGroup a) => Endianness -> Poly a -> [a]
vPolyCoeffs Endianness
BE -> [a]
v)
    = [a] -> Int -> Poly a
forall t. (Ord t, Num t) => [a] -> t -> Poly a
go [a]
u ([a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
u Int -> Int -> Int
forall a. Num a => a -> a -> a
- [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
v)
    where
        v0 :: a
v0  | [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [a]
v    = a
forall v. AdditiveGroup v => v
zeroV
            | Bool
otherwise = [a] -> a
forall a. [a] -> a
head [a]
v
        go :: [a] -> t -> Poly a
go u :: [a]
u n :: t
n
            | [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [a]
u Bool -> Bool -> Bool
|| t
n t -> t -> Bool
forall a. Ord a => a -> a -> Bool
< 0   = Endianness -> [a] -> Poly a
forall a. (Eq a, AdditiveGroup a) => Endianness -> [a] -> Poly a
poly Endianness
BE [a]
u
            | Bool
otherwise         = [a] -> t -> Poly a
go [a]
u' (t
nt -> t -> t
forall a. Num a => a -> a -> a
-1)
            where
                q0 :: a
q0 = a -> a -> a
divideV ([a] -> a
forall a. [a] -> a
head [a]
u) a
v0
                u' :: [a]
u' = [a] -> [a]
forall a. [a] -> [a]
tail ([a] -> [a] -> [a]
forall t. AdditiveGroup t => [t] -> [t] -> [t]
zipSumV [a]
u ((a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (a -> a -> a
multiplyV (a -> a
forall v. AdditiveGroup v => v -> v
negateV a
q0)) [a]
v))

-- |@composePoly f g@ constructs the polynomial 'h' such that:
-- 
-- > evalPoly h = evalPoly f . evalPoly g
-- 
-- This is a very expensive operation and, in general, returns a polynomial 
-- that is quite a bit more expensive to evaluate than @f@ and @g@ together
-- (because it is of a much higher order than either).  Unless your 
-- polynomials are quite small or you are quite certain you need the
-- coefficients of the composed polynomial, it is recommended that you 
-- simply evaluate @f@ and @g@ and explicitly compose the resulting 
-- functions.  This will usually be much more efficient.
composePolyWith :: (AdditiveGroup a, Eq a) => (a -> a -> a) -> Poly a -> Poly a -> Poly a
composePolyWith :: (a -> a -> a) -> Poly a -> Poly a -> Poly a
composePolyWith multiplyV :: a -> a -> a
multiplyV (Endianness -> Poly a -> [a]
forall a. (Eq a, AdditiveGroup a) => Endianness -> Poly a -> [a]
vPolyCoeffs Endianness
LE -> [a]
cs) (Endianness -> Poly a -> [a]
forall a. (Eq a, AdditiveGroup a) => Endianness -> Poly a -> [a]
vPolyCoeffs Endianness
LE -> [a]
ds) = Endianness -> [a] -> Poly a
forall a. (Eq a, AdditiveGroup a) => Endianness -> [a] -> Poly a
poly Endianness
LE ((a -> [a] -> [a]) -> [a] -> [a] -> [a]
forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr a -> [a] -> [a]
mul [] [a]
cs)
    where
        -- Implementation note: this is a hand-inlining of the following
        -- (with the 'Num' instance in "Math.Polynomial.NumInstance"):
        -- > composePoly f g = evalPoly (fmap constPoly f) g
        -- 
        -- This is a very expensive operation, something like
        -- O(length cs ^ 2 * length ds) I believe. There may be some more 
        -- tricks to improve that, but I suspect there isn't much room for 
        -- improvement. The number of terms in the resulting polynomial is 
        -- O(length cs * length ds) already, and each one is the sum of 
        -- quite a few terms.
        mul :: a -> [a] -> [a]
mul c :: a
c acc :: [a]
acc = a -> [a] -> [a]
forall a. (AdditiveGroup a, Eq a) => a -> [a] -> [a]
addScalarLE a
c ((a -> a -> a) -> [a] -> [a] -> [a]
forall a.
(AdditiveGroup a, Eq a) =>
(a -> a -> a) -> [a] -> [a] -> [a]
multPolyWithLE a -> a -> a
multiplyV [a]
acc [a]
ds)

-- |(internal) add a scalar to a list of polynomial coefficients in LE order
addScalarLE :: (AdditiveGroup a, Eq a) => a -> [a] -> [a]
addScalarLE :: a -> [a] -> [a]
addScalarLE a :: a
a bs :: [a]
bs | a
a a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
forall v. AdditiveGroup v => v
zeroV = [a]
bs
addScalarLE a :: a
a [] = [a
a]
addScalarLE a :: a
a (b :: a
b:bs :: [a]
bs) = (a
a a -> a -> a
forall v. AdditiveGroup v => v -> v -> v
^+^ a
b) a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [a]
bs

-- |Evaluate a polynomial at a point or, equivalently, convert a polynomial
-- to the function it represents.  For example, @evalPoly 'x' = 'id'@ and 
-- @evalPoly ('constPoly' k) = 'const' k.@
evalPoly :: (VectorSpace a, Eq a, AdditiveGroup (Scalar a), Eq (Scalar a)) => Poly a -> Scalar a -> a
evalPoly :: Poly a -> Scalar a -> a
evalPoly (Endianness -> Poly a -> [a]
forall a. (Eq a, AdditiveGroup a) => Endianness -> Poly a -> [a]
vPolyCoeffs Endianness
LE -> [a]
cs) x :: Scalar a
x
    | Scalar a
x Scalar a -> Scalar a -> Bool
forall a. Eq a => a -> a -> Bool
== Scalar a
forall v. AdditiveGroup v => v
zeroV =
        if [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [a]
cs
            then a
forall v. AdditiveGroup v => v
zeroV
            else [a] -> a
forall a. [a] -> a
head [a]
cs
    | Bool
otherwise = (a -> a -> a) -> a -> [a] -> a
forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr a -> a -> a
mul a
forall v. AdditiveGroup v => v
zeroV [a]
cs
    where
        mul :: a -> a -> a
mul c :: a
c acc :: a
acc = a
c a -> a -> a
forall v. AdditiveGroup v => v -> v -> v
^+^ a
acc a -> Scalar a -> a
forall v s. (VectorSpace v, s ~ Scalar v) => v -> s -> v
^* Scalar a
x

-- |Evaluate a polynomial and its derivative (respectively) at a point.
evalPolyDeriv :: (VectorSpace a, Eq a) => Poly a -> Scalar a -> (a,a)
evalPolyDeriv :: Poly a -> Scalar a -> (a, a)
evalPolyDeriv (Endianness -> Poly a -> [a]
forall a. (Eq a, AdditiveGroup a) => Endianness -> Poly a -> [a]
vPolyCoeffs Endianness
LE -> [a]
cs) x :: Scalar a
x = (a -> (a, a) -> (a, a)) -> (a, a) -> [a] -> (a, a)
forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr a -> (a, a) -> (a, a)
mul (a
forall v. AdditiveGroup v => v
zeroV, a
forall v. AdditiveGroup v => v
zeroV) [a]
cs
    where
        mul :: a -> (a, a) -> (a, a)
mul c :: a
c (p :: a
p, dp :: a
dp) = ((Scalar a
x Scalar a -> a -> a
forall v. VectorSpace v => Scalar v -> v -> v
*^ a
p) a -> a -> a
forall v. AdditiveGroup v => v -> v -> v
^+^ a
c, (Scalar a
x Scalar a -> a -> a
forall v. VectorSpace v => Scalar v -> v -> v
*^ a
dp) a -> a -> a
forall v. AdditiveGroup v => v -> v -> v
^+^ a
p)

-- |Evaluate a polynomial and all of its nonzero derivatives at a point.  
-- This is roughly equivalent to:
-- 
-- > evalPolyDerivs p x = map (`evalPoly` x) (takeWhile (not . polyIsZero) (iterate polyDeriv p))
evalPolyDerivs :: (VectorSpace a, Eq a, Num (Scalar a)) => Poly a -> Scalar a -> [a]
evalPolyDerivs :: Poly a -> Scalar a -> [a]
evalPolyDerivs (Endianness -> Poly a -> [a]
forall a. (Eq a, AdditiveGroup a) => Endianness -> Poly a -> [a]
vPolyCoeffs Endianness
LE -> [a]
cs) x :: Scalar a
x = [a] -> [a]
forall a. [a] -> [a]
trunc ([a] -> [a]) -> ([a] -> [a]) -> [a] -> [a]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Scalar a -> a -> a) -> [Scalar a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Scalar a -> a -> a
forall v. VectorSpace v => Scalar v -> v -> v
(*^) [Scalar a]
factorials ([a] -> [a]) -> [a] -> [a]
forall a b. (a -> b) -> a -> b
$ (a -> [a] -> [a]) -> [a] -> [a] -> [a]
forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr a -> [a] -> [a]
mul [] [a]
cs
    where
        trunc :: [c] -> [c]
trunc list :: [c]
list = (c -> a -> c) -> [c] -> [a] -> [c]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith c -> a -> c
forall a b. a -> b -> a
const [c]
list [a]
cs
        factorials :: [Scalar a]
factorials = (Scalar a -> Scalar a -> Scalar a)
-> Scalar a -> [Scalar a] -> [Scalar a]
forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl Scalar a -> Scalar a -> Scalar a
forall a. Num a => a -> a -> a
(*) 1 ((Scalar a -> Scalar a) -> Scalar a -> [Scalar a]
forall a. (a -> a) -> a -> [a]
iterate (Scalar a -> Scalar a -> Scalar a
forall a. Num a => a -> a -> a
+1) 1)
        mul :: a -> [a] -> [a]
mul c :: a
c pds :: [a]
pds@(p :: a
p:pd :: [a]
pd) = (Scalar a
x Scalar a -> a -> a
forall v. VectorSpace v => Scalar v -> v -> v
*^ a
p a -> a -> a
forall v. AdditiveGroup v => v -> v -> v
^+^ a
c) a -> [a] -> [a]
forall a. a -> [a] -> [a]
: (a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (Scalar a
x Scalar a -> a -> a
forall v. VectorSpace v => Scalar v -> v -> v
*^) [a]
pd [a] -> [a] -> [a]
forall t. AdditiveGroup t => [t] -> [t] -> [t]
`zipSumV` [a]
pds
        mul c :: a
c [] = [a
c]

-- |\"Contract\" a polynomial by attempting to divide out a root.
--
-- @contractPoly p a@ returns @(q,r)@ such that @q*(x-a) + r == p@
contractPoly :: (VectorSpace a, Eq a) => Poly a -> Scalar a -> (Poly a, a)
contractPoly :: Poly a -> Scalar a -> (Poly a, a)
contractPoly p :: Poly a
p@(Endianness -> Poly a -> [a]
forall a. (Eq a, AdditiveGroup a) => Endianness -> Poly a -> [a]
vPolyCoeffs Endianness
LE -> [a]
cs) a :: Scalar a
a = (Int -> Endianness -> [a] -> Poly a
forall a.
(Eq a, AdditiveGroup a) =>
Int -> Endianness -> [a] -> Poly a
vPolyN Int
n Endianness
LE [a]
q, a
r)
    where
        n :: Int
n = Poly a -> Int
forall a. Poly a -> Int
rawPolyLength Poly a
p
        cut :: a -> a -> (a, a)
cut remainder :: a
remainder swap :: a
swap = (a
swap a -> a -> a
forall v. AdditiveGroup v => v -> v -> v
^+^ (Scalar a
a Scalar a -> a -> a
forall v. VectorSpace v => Scalar v -> v -> v
*^ a
remainder), a
remainder)
        (r :: a
r,q :: [a]
q) = (a -> a -> (a, a)) -> a -> [a] -> (a, [a])
forall (t :: * -> *) a b c.
Traversable t =>
(a -> b -> (a, c)) -> a -> t b -> (a, t c)
mapAccumR a -> a -> (a, a)
cut a
forall v. AdditiveGroup v => v
zeroV [a]
cs

-- |@gcdPoly a b@ computes the highest order monic polynomial that is a 
-- divisor of both @a@ and @b@.  If both @a@ and @b@ are 'zero', the 
-- result is undefined.
gcdPolyWith :: (AdditiveGroup a, Eq a) => a -> (a -> a -> a) -> (a -> a -> a) -> Poly a -> Poly a -> Poly a
gcdPolyWith :: a -> (a -> a -> a) -> (a -> a -> a) -> Poly a -> Poly a -> Poly a
gcdPolyWith oneV :: a
oneV multiplyV :: a -> a -> a
multiplyV divideV :: a -> a -> a
divideV a :: Poly a
a b :: Poly a
b 
    | Poly a -> Bool
forall a. (Eq a, AdditiveGroup a) => Poly a -> Bool
polyIsZero Poly a
b  = if Poly a -> Bool
forall a. (Eq a, AdditiveGroup a) => Poly a -> Bool
polyIsZero Poly a
a
        then [Char] -> Poly a
forall a. HasCallStack => [Char] -> a
error "gcdPolyWith: gcdPoly zero zero is undefined"
        else a -> (a -> a -> a) -> Poly a -> Poly a
forall a.
(AdditiveGroup a, Eq a) =>
a -> (a -> a -> a) -> Poly a -> Poly a
monicPolyWith a
oneV a -> a -> a
divideV Poly a
a
    | Bool
otherwise     = a -> (a -> a -> a) -> (a -> a -> a) -> Poly a -> Poly a -> Poly a
forall a.
(AdditiveGroup a, Eq a) =>
a -> (a -> a -> a) -> (a -> a -> a) -> Poly a -> Poly a -> Poly a
gcdPolyWith a
oneV a -> a -> a
multiplyV a -> a -> a
divideV Poly a
b (Poly a
a Poly a -> Poly a -> Poly a
`remPoly` Poly a
b)
    where remPoly :: Poly a -> Poly a -> Poly a
remPoly = (a -> a -> a) -> (a -> a -> a) -> Poly a -> Poly a -> Poly a
forall a.
(AdditiveGroup a, Eq a) =>
(a -> a -> a) -> (a -> a -> a) -> Poly a -> Poly a -> Poly a
remPolyWith a -> a -> a
multiplyV a -> a -> a
divideV

-- |Normalize a polynomial so that its highest-order coefficient is 1
monicPolyWith :: (AdditiveGroup a, Eq a) => a -> (a -> a -> a) -> Poly a -> Poly a
monicPolyWith :: a -> (a -> a -> a) -> Poly a -> Poly a
monicPolyWith oneV :: a
oneV divideV :: a -> a -> a
divideV p :: Poly a
p = case Endianness -> Poly a -> [a]
forall a. (Eq a, AdditiveGroup a) => Endianness -> Poly a -> [a]
vPolyCoeffs Endianness
BE Poly a
p of
    []      -> Int -> Endianness -> [a] -> Poly a
forall a.
(Eq a, AdditiveGroup a) =>
Int -> Endianness -> [a] -> Poly a
vPolyN Int
n Endianness
BE []
    (c :: a
c:cs :: [a]
cs)  -> Int -> Endianness -> [a] -> Poly a
forall a.
(Eq a, AdditiveGroup a) =>
Int -> Endianness -> [a] -> Poly a
vPolyN Int
n Endianness
BE (a
oneV a -> [a] -> [a]
forall a. a -> [a] -> [a]
: (a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (a -> a -> a
`divideV` a
c) [a]
cs)
    where n :: Int
n = Poly a -> Int
forall a. Poly a -> Int
rawPolyLength Poly a
p

-- |Compute the derivative of a polynomial.
polyDeriv :: (VectorSpace a, Eq a, Num (Scalar a)) => Poly a -> Poly a
polyDeriv :: Poly a -> Poly a
polyDeriv p :: Poly a
p@(Endianness -> Poly a -> [a]
forall a. (Eq a, AdditiveGroup a) => Endianness -> Poly a -> [a]
vPolyCoeffs Endianness
LE -> [a]
cs) = Int -> Endianness -> [a] -> Poly a
forall a.
(Eq a, AdditiveGroup a) =>
Int -> Endianness -> [a] -> Poly a
vPolyN (Poly a -> Int
forall a. Poly a -> Int
rawPolyDegree Poly a
p) Endianness
LE
    [ Scalar a
n Scalar a -> a -> a
forall v. VectorSpace v => Scalar v -> v -> v
*^ a
c
    | a
c <- Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
drop 1 [a]
cs
    | Scalar a
n <- (Scalar a -> Scalar a) -> Scalar a -> [Scalar a]
forall a. (a -> a) -> a -> [a]
iterate (1Scalar a -> Scalar a -> Scalar a
forall a. Num a => a -> a -> a
+) 1
    ]

-- |Compute all nonzero derivatives of a polynomial, starting with its 
-- \"zero'th derivative\", the original polynomial itself.
polyDerivs :: (VectorSpace a, Eq a, Num (Scalar a)) => Poly a -> [Poly a]
polyDerivs :: Poly a -> [Poly a]
polyDerivs p :: Poly a
p = Int -> [Poly a] -> [Poly a]
forall a. Int -> [a] -> [a]
take (1 Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Poly a -> Int
forall a. (Eq a, AdditiveGroup a) => Poly a -> Int
polyDegree Poly a
p) ((Poly a -> Poly a) -> Poly a -> [Poly a]
forall a. (a -> a) -> a -> [a]
iterate Poly a -> Poly a
forall a. (VectorSpace a, Eq a, Num (Scalar a)) => Poly a -> Poly a
polyDeriv Poly a
p)


-- |Compute the definite integral (from 0 to x) of a polynomial.
polyIntegral :: (VectorSpace a, Eq a, Fractional (Scalar a)) => Poly a -> Poly a
polyIntegral :: Poly a -> Poly a
polyIntegral p :: Poly a
p@(Endianness -> Poly a -> [a]
forall a. (Eq a, AdditiveGroup a) => Endianness -> Poly a -> [a]
vPolyCoeffs Endianness
LE -> [a]
cs) = Int -> Endianness -> [a] -> Poly a
forall a.
(Eq a, AdditiveGroup a) =>
Int -> Endianness -> [a] -> Poly a
vPolyN (1 Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Poly a -> Int
forall a. Poly a -> Int
rawPolyLength Poly a
p) Endianness
LE ([a] -> Poly a) -> [a] -> Poly a
forall a b. (a -> b) -> a -> b
$ a
forall v. AdditiveGroup v => v
zeroV a -> [a] -> [a]
forall a. a -> [a] -> [a]
:
    [ a
c a -> Scalar a -> a
forall v s.
(VectorSpace v, s ~ Scalar v, Fractional s) =>
v -> s -> v
^/ Scalar a
n
    | a
c <- [a]
cs
    | Scalar a
n <- (Scalar a -> Scalar a) -> Scalar a -> [Scalar a]
forall a. (a -> a) -> a -> [a]
iterate (1Scalar a -> Scalar a -> Scalar a
forall a. Num a => a -> a -> a
+) 1
    ]