packages feed

polynomial 0.6.5 → 0.7.1

raw patch · 6 files changed

+398/−135 lines, 6 filesdep +vector-th-unboxdep ~base

Dependencies added: vector-th-unbox

Dependency ranges changed: base

Files

polynomial.cabal view
@@ -1,5 +1,5 @@ name:                   polynomial-version:                0.6.5+version:                0.7.1 stability:              provisional  cabal-version:          >= 1.6@@ -22,9 +22,12 @@  Library   ghc-options:          -Wall -fno-warn-name-shadowing+  if impl(ghc >= 7.4)+    ghc-options:        -fwarn-unsafe   hs-source-dirs:       src   exposed-modules:      Math.Polynomial                         Math.Polynomial.Bernstein+                        Math.Polynomial.Bernoulli                         Math.Polynomial.Chebyshev                         Math.Polynomial.Hermite                         Math.Polynomial.Interpolation@@ -33,7 +36,10 @@                         Math.Polynomial.Newton                         Math.Polynomial.NumInstance                         Math.Polynomial.Type+                        Math.Polynomial.VectorSpace   other-modules:        Data.List.ZipSum+                        Data.VectorSpace.WrappedNum                         Math.Polynomial.Pretty                         -  build-depends:        base >= 3 && <5, deepseq, pretty, prettyclass, vector, vector-space+  build-depends:        base >= 3 && <5, deepseq, pretty, prettyclass, vector, vector-space,+                        vector-th-unbox >= 0.2.1
+ src/Data/VectorSpace/WrappedNum.hs view
@@ -0,0 +1,27 @@+{-# LANGUAGE TemplateHaskell, TypeFamilies, MultiParamTypeClasses, FlexibleInstances, FlexibleContexts #-}+{-# LANGUAGE GeneralizedNewtypeDeriving #-}+module Data.VectorSpace.WrappedNum+  (WrappedNum(..)) where++import Data.VectorSpace+import qualified Data.Vector.Unboxed         as U++import Data.Vector.Unboxed.Deriving++newtype WrappedNum a = WrapNum { unwrapNum :: a }+    deriving+        (Eq, Ord, Read, Show, Bounded+        , Enum, Num, Fractional, Real, RealFrac+        , Floating, RealFloat)++derivingUnbox "Wrapped"+    [t| (U.Unbox a) => WrappedNum a -> a |] [| unwrapNum |] [| \ a -> WrapNum a |]++instance Num a => AdditiveGroup (WrappedNum a) where+    zeroV = 0+    (^+^) = (+)+    negateV = negate++instance Num a => VectorSpace (WrappedNum a) where+    type Scalar (WrappedNum a) = WrappedNum a+    (*^) = (*)
src/Math/Polynomial.hs view
@@ -19,112 +19,67 @@ import Math.Polynomial.Type import Math.Polynomial.Pretty ({- instance -}) -import Data.List-import Data.List.ZipSum---- |The polynomial \"1\"-one :: (Num a, Eq a) => Poly a-one = constPoly 1---- |The polynomial (in x) \"x\"-x :: (Num a, Eq a) => Poly a-x = polyN 2 LE [0,1]+import Math.Polynomial.VectorSpace (one, x) -- to re-export+import qualified Math.Polynomial.VectorSpace as VS+import Data.VectorSpace.WrappedNum  -- |Given some constant 'k', construct the polynomial whose value is  -- constantly 'k'. constPoly :: (Num a, Eq a) => a -> Poly a-constPoly x = polyN 1 LE [x]+constPoly x = unwrapPoly (VS.constPoly (WrapNum x))  -- |Given some scalar 's' and a polynomial 'f', computes the polynomial 'g' -- such that: --  -- > evalPoly g x = s * evalPoly f x scalePoly :: (Num a, Eq a) => a -> Poly a -> Poly a-scalePoly 0 _ = zero-scalePoly s p = mapPoly (s*) p+scalePoly x f = unwrapPoly (VS.scalePoly (WrapNum x) (wrapPoly f))  -- |Given some polynomial 'f', computes the polynomial 'g' such that: --  -- > evalPoly g x = negate (evalPoly f x)-negatePoly :: Num a => Poly a -> Poly a-negatePoly = mapPoly negate+negatePoly :: (Num a, Eq a) => Poly a -> Poly a+negatePoly f = unwrapPoly (VS.negatePoly (wrapPoly f))  -- |Given polynomials 'f' and 'g', computes the polynomial 'h' such that: --  -- > evalPoly h x = evalPoly f x + evalPoly g x 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)+addPoly p q = unwrapPoly (VS.addPoly (wrapPoly p) (wrapPoly q))  {-# RULES   "sum Poly"    forall ps. foldl addPoly zero ps = sumPolys ps   #-} sumPolys :: (Num a, Eq a) => [Poly a] -> Poly a-sumPolys [] = zero-sumPolys ps = poly LE (foldl1 zipSum (map (polyCoeffs LE) ps))+sumPolys ps = unwrapPoly (VS.sumPolys (map wrapPoly ps))  -- |Given polynomials 'f' and 'g', computes the polynomial 'h' such that: --  -- > evalPoly h x = evalPoly f x * evalPoly g x 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---- |(Internal): multiply polynomials in LE order.  O(length xs * length ys).-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+multPoly p q = unwrapPoly (VS.multPolyWith (*) (wrapPoly p) (wrapPoly q))  -- |Given a polynomial 'f' and exponent 'n', computes the polynomial 'g' -- such that: --  -- > evalPoly g x = evalPoly f x ^ n 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 (n-1)-    | otherwise = (\x -> multPoly x x) (powPoly p (n`div`2))+powPoly p n = unwrapPoly (VS.powPolyWith 1 (*) (wrapPoly p) n)  -- |Given polynomials @a@ and @b@, with @b@ not 'zero', computes polynomials -- @q@ and @r@ such that: --  -- > addPoly (multPoly q b) r == a 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)+quotRemPoly u v = (unwrapPoly q, unwrapPoly r)     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' (n-1)-            where-                q0 = head u / v0-                u' = tail (zipSum u (map (negate q0 *) v))+        ~(q, r) = VS.quotRemPolyWith (*) (/) (wrapPoly u) (wrapPoly 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)+quotPoly u v = unwrapPoly (VS.quotPolyWith (*) (/) (wrapPoly u) (wrapPoly 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' (n-1)-            where-                q0 = head u / v0-                u' = tail (zipSum u (map (negate q0 *) v))+remPoly u v = unwrapPoly (VS.remPolyWith (*) (/) (wrapPoly u) (wrapPoly v))  -- |@composePoly f g@ constructs the polynomial 'h' such that: -- @@ -138,103 +93,57 @@ -- simply evaluate @f@ and @g@ and explicitly compose the resulting  -- functions.  This will usually be much more efficient. composePoly :: (Num a, Eq a) => Poly a -> Poly a -> Poly a-composePoly (polyCoeffs LE -> cs) (polyCoeffs LE -> ds) = poly LE (foldr mul [] 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 c acc = addScalarLE c (multPolyLE acc ds)---- |(internal) add a scalar to a list of polynomial coefficients in LE order-addScalarLE :: (Num a, Eq a) => a -> [a] -> [a]-addScalarLE 0 bs = bs-addScalarLE a [] = [a]-addScalarLE a (b:bs) = (a + b) : bs+composePoly p q = unwrapPoly (VS.composePolyWith (*) (wrapPoly p) (wrapPoly q))  -- |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 :: (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+evalPoly f x = unwrapNum (VS.evalPoly (wrapPoly f) (WrapNum x))  -- |Evaluate a polynomial and its derivative (respectively) at a point. evalPolyDeriv :: (Num a, Eq a) => Poly a -> a -> (a,a)-evalPolyDeriv (polyCoeffs LE -> cs) x = foldr mul (0,0) cs+evalPolyDeriv f x = (unwrapNum y, unwrapNum y')     where-        mul c (p, dp) = (p * x + c, dp * x + p)+        ~(y, y') = VS.evalPolyDeriv (wrapPoly f) (WrapNum x)  -- |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 :: (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]+evalPolyDerivs f x = map unwrapNum (VS.evalPolyDerivs (wrapPoly f) (WrapNum x))  -- |\"Contract\" a polynomial by attempting to divide out a root. -- -- @contractPoly p a@ returns @(q,r)@ such that @q*(x-a) + r == p@ contractPoly :: (Num a, Eq a) => Poly a -> a -> (Poly a, a)-contractPoly p@(polyCoeffs LE -> cs) a = (polyN n LE q, r)+contractPoly p a = (unwrapPoly q, unwrapNum r)     where-        n = rawPolyLength p-        cut remainder swap = (swap + remainder * a, remainder)-        (r,q) = mapAccumR cut 0 cs+        (q, r) = VS.contractPoly (wrapPoly p) (WrapNum a)  -- |@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. 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)+gcdPoly a b = unwrapPoly (VS.gcdPolyWith 1 (*) (/) (wrapPoly a) (wrapPoly b))  -- |Normalize a polynomial so that its highest-order coefficient is 1 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+monicPoly p = unwrapPoly (VS.monicPolyWith 1 (/) (wrapPoly p))  -- |Compute the derivative of a polynomial. 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-    ]+polyDeriv p = unwrapPoly (VS.polyDeriv (wrapPoly p))  -- |Compute all nonzero derivatives of a polynomial, starting with its  -- \"zero'th derivative\", the original polynomial itself. polyDerivs :: (Num a, Eq a) => Poly a -> [Poly a]-polyDerivs p = take (1 + polyDegree p) (iterate polyDeriv p)-+polyDerivs p = map unwrapPoly (VS.polyDerivs (wrapPoly p))  -- |Compute the definite integral (from 0 to x) of a polynomial. 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-    ]+polyIntegral p = unwrapPoly (VS.polyIntegral (wrapPoly p))  -- |Separate a nonzero polynomial into a set of factors none of which have -- multiple roots, and the product of which is the original polynomial.
+ src/Math/Polynomial/Bernoulli.hs view
@@ -0,0 +1,25 @@+module Math.Polynomial.Bernoulli (bernoulliPoly) where++import Math.Polynomial+import Data.VectorSpace++{- | Bernoulli polynomial with a nonstandard normalization++> b_i = bernoulliPoly !! i++Has the following generating function (C.2 in IH Sloan & S Joe+"Lattice Methods for multiple integration" 1994 page 227)++> t exp(x*t) / (exp(t) - 1) = sum_{i=0} b_i t^i++The standard normalization would have @= sum_{i=0} B_i t^i / i!@++-}+bernoulliPoly :: (Fractional a, Eq a) => [Poly a]+bernoulliPoly = map fst biIntegralBi++biIntegralBi :: (Fractional a, Eq a) => [(Poly a, Poly a)]+biIntegralBi = (constPoly 1, polyIntegral (constPoly 1)) : map f biIntegralBi+  where f (p, ip) = case polyIntegral ip of+                      ip2 -> case constPoly $ evalPoly ip2 0 - evalPoly ip2 1 of+                               c -> (c `addPoly` ip, polyIntegral c `addPoly` ip2)
src/Math/Polynomial/Type.hs view
@@ -1,4 +1,4 @@-{-# LANGUAGE ViewPatterns, TypeFamilies, GADTs #-}+{-# LANGUAGE ViewPatterns, TypeFamilies, GADTs, UndecidableInstances #-} -- |Low-level interface for the 'Poly' type. module Math.Polynomial.Type      ( Endianness(..)@@ -10,6 +10,9 @@     , unboxedPoly, unboxedPolyN          , mapPoly+    , rawMapPoly+    , wrapPoly+    , unwrapPoly          , unboxPoly     @@ -18,11 +21,13 @@     , rawVectorPoly     , rawUVectorPoly     , trim+    , vTrim          , polyIsZero     , polyIsOne          , polyCoeffs+    , vPolyCoeffs     , rawCoeffsOrder     , rawPolyCoeffs     , untrimmedPolyCoeffs@@ -37,10 +42,15 @@ -- import Data.List.Extras.LazyLength import Data.AdditiveGroup import Data.VectorSpace+import Data.VectorSpace.WrappedNum import Data.List.ZipSum import qualified Data.Vector as V import qualified Data.Vector.Unboxed as UV +-- 'unsafeCoerce' is only used in 'wrapPoly' and 'unwrapPoly', which are+-- type-safe alternatives to 'fmap'ing the 'WrappedNum' newtype constructor/projector+import Unsafe.Coerce (unsafeCoerce)+ data Endianness      = BE      -- ^ Big-Endian (head is highest-order term)@@ -84,14 +94,15 @@  -- TODO: specialize for case where one is a list and other is a vector; --  use native order of the list-instance (Num a, Eq a) => Eq (Poly a) where+-- TODO: think about plain Num support...+instance (AdditiveGroup a, Eq a) => Eq (Poly a) where     p == q           | rawCoeffsOrder p == rawCoeffsOrder q-        =  rawPolyCoeffs (trim (0==) p) -        == rawPolyCoeffs (trim (0==) q)+        =  rawPolyCoeffs (trim (zeroV==) p) +        == rawPolyCoeffs (trim (zeroV==) q)         | otherwise -        =  polyCoeffs LE p-        == polyCoeffs LE q+        =  vPolyCoeffs LE p+        == vPolyCoeffs LE q  -- -- Ord would be nice for some purposes, but it really just doesn't -- -- make sense (there is no natural order that is much better than any@@ -112,22 +123,38 @@     fmap f (UVectorPoly _ end cs) = VectorPoly False end (V.fromListN n . map f $ UV.toList cs)         where n = UV.length cs --- TODO: this needs to be renamed 'rawMapPoly' and wrapped with 'trim'. -- |Like fmap, but able to preserve unboxedness-mapPoly :: (a -> a) -> Poly a -> Poly a-mapPoly f (ListPoly    _ e cs) = ListPoly    False e (   map f cs)-mapPoly f (VectorPoly  _ e cs) = VectorPoly  False e ( V.map f cs)-mapPoly f (UVectorPoly _ e cs) = UVectorPoly False e (UV.map f cs)+mapPoly :: (Num a, Eq a) => (a -> a) -> Poly a -> Poly a+mapPoly f = trim (0==) . rawMapPoly f +rawMapPoly :: (a -> a) -> Poly a -> Poly a+rawMapPoly f (ListPoly    _ e cs) = ListPoly    False e (   map f cs)+rawMapPoly f (VectorPoly  _ e cs) = VectorPoly  False e ( V.map f cs)+rawMapPoly f (UVectorPoly _ e cs) = UVectorPoly False e (UV.map f cs)++{-# RULES "wrapPoly/unwrapPoly"   forall x. wrapPoly (unwrapPoly x) = x #-}+{-# RULES "unwrapPoly/wrapPoly"   forall x. unwrapPoly (wrapPoly x) = x #-}+{-# RULES "wrapPoly.unwrapPoly"   wrapPoly . unwrapPoly = id #-}+{-# RULES "unwrapPoly.wrapPoly"   unwrapPoly . wrapPoly = id #-}+-- |like @fmap WrapNum@ but using 'unsafeCoerce' to avoid a pointless traversal+wrapPoly :: Poly a -> Poly (WrappedNum a)+wrapPoly = unsafeCoerce++-- |like @fmap unwrapNum@ but using 'unsafeCoerce' to avoid a pointless traversal+unwrapPoly :: Poly (WrappedNum a) -> Poly a+unwrapPoly = unsafeCoerce+ instance AdditiveGroup a => AdditiveGroup (Poly a) where     zeroV = ListPoly True LE []     (untrimmedPolyCoeffs LE ->  a) ^+^ (untrimmedPolyCoeffs LE ->  b)          = ListPoly False LE (zipSumV a b)     negateV = fmap negateV -instance VectorSpace a => VectorSpace (Poly a) where+instance (Eq a, VectorSpace a, AdditiveGroup (Scalar a), Eq (Scalar a)) => VectorSpace (Poly a) where     type Scalar (Poly a) = Scalar a-    (*^) s = fmap (s *^)+    s *^ v+         | s == zeroV   = zeroV+         | otherwise    = vTrim (rawMapPoly (s *^) v)  -- |Trim zeroes from a polynomial (given a predicate for identifying zero). -- In particular, drops zeroes from the highest-order coefficients, so that@@ -145,6 +172,9 @@ trim isZero   (UVectorPoly _ LE cs) = UVectorPoly True LE (UV.reverse . UV.dropWhile isZero . UV.reverse $ cs) trim isZero   (UVectorPoly _ BE cs) = UVectorPoly True BE (UV.dropWhile isZero cs) +vTrim :: (Eq a, AdditiveGroup a) => Poly a -> Poly a+vTrim = trim (zeroV ==)+ -- |The polynomial \"0\" zero :: Poly a zero = ListPoly True LE []@@ -198,7 +228,11 @@  -- |Get the coefficients of a a 'Poly' in the specified order. polyCoeffs :: (Num a, Eq a) => Endianness -> Poly a -> [a]-polyCoeffs end p = untrimmedPolyCoeffs end (trim (0==) p)+polyCoeffs end p = untrimmedPolyCoeffs end (trim (0 ==) p)++-- |Get the coefficients of a a 'Poly' in the specified order.+vPolyCoeffs :: (Eq a, AdditiveGroup a) => Endianness -> Poly a -> [a]+vPolyCoeffs end p = untrimmedPolyCoeffs end (vTrim p)  polyIsZero :: (Num a, Eq a) => Poly a -> Bool polyIsZero = null . rawPolyCoeffs . trim (0==)
+ src/Math/Polynomial/VectorSpace.hs view
@@ -0,0 +1,262 @@+{-# 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 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++-- |The polynomial \"1\"+one :: (Num a, Eq a) => Poly a+one = polyN 1 LE [1]++-- |The polynomial (in x) \"x\"+x :: (Num a, Eq a) => Poly a+x = polyN 2 LE [0,1]++-- |Given some constant 'k', construct the polynomial whose value is +-- constantly 'k'.+constPoly :: (Eq a, AdditiveGroup a) => a -> Poly a+constPoly x = vPolyN 1 LE [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 = (*^)++-- |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 = vTrim . rawMapPoly 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 p@(vPolyCoeffs LE ->  a) q@(vPolyCoeffs LE ->  b) = vPolyN n LE (zipSumV a b)+    where n = max (rawPolyLength p) (rawPolyLength q)++{-# RULES+  "sum Poly"    forall ps. foldl addPoly zero ps = sumPolys ps+  #-}+sumPolys :: (AdditiveGroup a, Eq a) => [Poly a] -> Poly a+sumPolys [] = zero+sumPolys ps = poly LE (foldl1 zipSumV (map (vPolyCoeffs LE) 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 multiplyV p@(vPolyCoeffs LE -> xs) q@(vPolyCoeffs LE -> ys) = vPolyN n LE (multPolyWithLE multiplyV xs ys)+    where n = 1 + rawPolyDegree p + rawPolyDegree q++-- |(Internal): multiply polynomials in LE order.  O(length xs * length ys).+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++-- |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 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 (n-1)+            | otherwise = (\x -> multPoly x x) (powPoly p (n`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 _ _ _ 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' (n-1)+            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' (n-1)+            where+                q0 = divideV (head u) v0+                u' = tail (zipSumV u (map (multiplyV (negateV q0)) 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 multiplyV (vPolyCoeffs LE -> cs) (vPolyCoeffs LE -> ds) = poly LE (foldr mul [] 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 c acc = addScalarLE c (multPolyWithLE multiplyV acc ds)++-- |(internal) add a scalar to a list of polynomial coefficients in LE order+addScalarLE :: (AdditiveGroup a, Eq a) => a -> [a] -> [a]+addScalarLE a bs | a == zeroV = bs+addScalarLE a [] = [a]+addScalarLE a (b:bs) = (a ^+^ b) : 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 (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++-- |Evaluate a polynomial and its derivative (respectively) at a point.+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)++-- |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 (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]++-- |\"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 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++-- |@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 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++-- |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 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++-- |Compute the derivative of a polynomial.+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+    ]++-- |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 p = take (1 + polyDegree p) (iterate polyDeriv 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 p@(vPolyCoeffs LE -> cs) = vPolyN (1 + rawPolyLength p) LE $ zeroV :+    [ c ^/ n+    | c <- cs+    | n <- iterate (1+) 1+    ]