packages feed

polynomial 0.6 → 0.6.5

raw patch · 10 files changed

+305/−138 lines, 10 filesdep +vectordep ~base

Dependencies added: vector

Dependency ranges changed: base

Files

polynomial.cabal view
@@ -1,5 +1,5 @@ name:                   polynomial-version:                0.6+version:                0.6.5 stability:              provisional  cabal-version:          >= 1.6@@ -8,7 +8,7 @@ author:                 James Cook <mokus@deepbondi.net> maintainer:             James Cook <mokus@deepbondi.net> license:                PublicDomain-homepage:               /dev/null+homepage:               https://github.com/mokus0/polynomial  category:               Math, Numerical synopsis:               Polynomials@@ -17,8 +17,8 @@                         interesting polynomial sequences.  source-repository head-  type: darcs-  location: http://code.haskell.org/~mokus/polynomial+  type: git+  location: git://github.com/mokus0/polynomial.git  Library   ghc-options:          -Wall -fno-warn-name-shadowing@@ -26,13 +26,14 @@   exposed-modules:      Math.Polynomial                         Math.Polynomial.Bernstein                         Math.Polynomial.Chebyshev+                        Math.Polynomial.Hermite                         Math.Polynomial.Interpolation                         Math.Polynomial.Lagrange                         Math.Polynomial.Legendre                         Math.Polynomial.Newton                         Math.Polynomial.NumInstance-  other-modules:        Data.List.ZipSum                         Math.Polynomial.Type+  other-modules:        Data.List.ZipSum                         Math.Polynomial.Pretty                         -  build-depends:        base >= 3 && <5, deepseq, pretty, prettyclass, vector-space+  build-depends:        base >= 3 && <5, deepseq, pretty, prettyclass, vector, vector-space
src/Math/Polynomial.hs view
@@ -2,7 +2,8 @@ {-# OPTIONS_GHC -fno-warn-orphans #-} module Math.Polynomial     ( Endianness(..)-    , Poly, poly, polyCoeffs, polyIsZero, polyIsOne+    , Poly, poly, polyDegree+    , polyCoeffs, polyIsZero, polyIsOne     , zero, one, constPoly, x     , scalePoly, negatePoly     , composePoly@@ -10,8 +11,9 @@     , quotRemPoly, quotPoly, remPoly     , evalPoly, evalPolyDeriv, evalPolyDerivs     , contractPoly+    , monicPoly     , gcdPoly, separateRoots-    , polyDeriv, polyIntegral+    , polyDeriv, polyDerivs, polyIntegral     ) where  import Math.Polynomial.Type@@ -20,58 +22,56 @@ import Data.List import Data.List.ZipSum --- |The polynomial \"0\"-zero :: Num a => Poly a-zero = poly LE []- -- |The polynomial \"1\"-one :: Num a => Poly a+one :: (Num a, Eq a) => Poly a one = constPoly 1  -- |The polynomial (in x) \"x\"-x :: Num a => Poly a-x = poly LE [0,1]+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 :: Num a => a -> Poly a-constPoly x = poly LE [x]+constPoly :: (Num a, Eq a) => a -> Poly a+constPoly x = polyN 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 :: Num a => a -> Poly a -> Poly a+scalePoly :: (Num a, Eq a) => a -> Poly a -> Poly a scalePoly 0 _ = zero-scalePoly s p = fmap (s*) p+scalePoly s p = mapPoly (s*) p  -- |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 = fmap negate+negatePoly = mapPoly negate  -- |Given polynomials 'f' and 'g', computes the polynomial 'h' such that: --  -- > evalPoly h x = evalPoly f x + evalPoly g x-addPoly :: Num a => Poly a -> Poly a -> Poly a-addPoly (polyCoeffs LE ->  a) (polyCoeffs LE ->  b) = poly LE (zipSum a b)+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)  {-# RULES   "sum Poly"    forall ps. foldl addPoly zero ps = sumPolys ps   #-}-sumPolys :: Num a => [Poly a] -> Poly a+sumPolys :: (Num a, Eq a) => [Poly a] -> Poly a sumPolys [] = zero sumPolys ps = poly LE (foldl1 zipSum (map (polyCoeffs LE) ps))  -- |Given polynomials 'f' and 'g', computes the polynomial 'h' such that: --  -- > evalPoly h x = evalPoly f x * evalPoly g x-multPoly :: Num a => Poly a -> Poly a -> Poly a-multPoly (polyCoeffs LE -> xs) (polyCoeffs LE -> ys) = poly LE (multPolyLE xs ys)+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 => [a] -> [a] -> [a]+multPolyLE :: (Num a, Eq a) => [a] -> [a] -> [a] multPolyLE _  []     = [] multPolyLE xs (y:ys) = foldr mul [] xs     where@@ -82,8 +82,8 @@ -- such that: --  -- > evalPoly g x = evalPoly f x ^ n-powPoly :: (Num a, Integral b) => Poly a -> b -> Poly a-powPoly _ 0 = poly LE [1]+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"@@ -94,10 +94,10 @@ -- @q@ and @r@ such that: --  -- > addPoly (multPoly q b) r == a-quotRemPoly :: Fractional a => Poly a -> Poly a -> (Poly a, Poly a)+quotRemPoly :: (Fractional a, Eq a) => Poly a -> Poly a -> (Poly a, Poly a) quotRemPoly _ b | polyIsZero b = error "quotRemPoly: divide by zero"-quotRemPoly (polyCoeffs BE -> u) (polyCoeffs BE -> v)-    = go [] u (length u - length v)+quotRemPoly p@(polyCoeffs BE -> u) q@(polyCoeffs BE -> v)+    = go [] u (polyDegree p - polyDegree q)     where         v0  | null v    = 0             | otherwise = head v@@ -108,11 +108,11 @@                 q0 = head u / v0                 u' = tail (zipSum u (map (negate q0 *) v)) -quotPoly :: Fractional a => Poly a -> Poly a -> Poly a+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 => Poly a -> Poly a -> Poly a+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)@@ -137,7 +137,7 @@ -- 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.-composePoly :: Num a => Poly a -> Poly a -> Poly a+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@@ -153,7 +153,7 @@         mul c acc = addScalarLE c (multPolyLE acc ds)  -- |(internal) add a scalar to a list of polynomial coefficients in LE order-addScalarLE :: Num a => a -> [a] -> [a]+addScalarLE :: (Num a, Eq a) => a -> [a] -> [a] addScalarLE 0 bs = bs addScalarLE a [] = [a] addScalarLE a (b:bs) = (a + b) : bs@@ -161,7 +161,7 @@ -- |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 => Poly a -> a -> a+evalPoly :: (Num a, Eq a) => Poly a -> a -> a evalPoly (polyCoeffs LE -> cs) 0     | null cs   = 0     | otherwise = head cs@@ -170,7 +170,7 @@         mul c acc = c + acc * x  -- |Evaluate a polynomial and its derivative (respectively) at a point.-evalPolyDeriv :: Num a => Poly a -> a -> (a,a)+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)@@ -179,7 +179,7 @@ -- This is roughly equivalent to: --  -- > evalPolyDerivs p x = map (`evalPoly` x) (takeWhile (not . polyIsZero) (iterate polyDeriv p))-evalPolyDerivs :: Num a => Poly a -> a -> [a]+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@@ -190,39 +190,47 @@ -- |\"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 => Poly a -> a -> (Poly a, a)-contractPoly (polyCoeffs LE -> cs) a = (poly LE q, r)+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 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 => Poly a -> Poly a -> Poly a+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 monic a+        else monicPoly a     | otherwise     = gcdPoly b (a `remPoly` b) --- |(internal) Normalize a polynomial so that its highest-order coefficient is 1-monic :: Fractional a => Poly a -> Poly a-monic p = case polyCoeffs BE p of-    []      -> poly BE []-    (c:cs)  -> poly BE (1:map (/c) cs)+-- |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  -- |Compute the derivative of a polynomial.-polyDeriv :: Num a => Poly a -> Poly a-polyDeriv (polyCoeffs LE -> cs) = poly LE+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     ] +-- |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)++ -- |Compute the definite integral (from 0 to x) of a polynomial.-polyIntegral :: Fractional a => Poly a -> Poly a-polyIntegral (polyCoeffs LE -> cs) = poly LE $ 0 :+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@@ -234,7 +242,7 @@ -- Rational coefficients is a good idea. -- -- Useful when applicable as a way to simplify root-finding problems.-separateRoots :: Fractional a => Poly a -> [Poly a]+separateRoots :: (Fractional a, Eq a) => Poly a -> [Poly a] separateRoots p     | polyIsZero q  = error "separateRoots: zero polynomial"     | polyIsOne q   = [p]
src/Math/Polynomial/Chebyshev.hs view
@@ -22,12 +22,12 @@     ]  -- |Compute the coefficients of the n'th Chebyshev polynomial of the first kind.-t :: Num a => Int -> Poly a+t :: (Num a, Eq a) => Int -> Poly a t n | n >= 0    = poly LE . map fromInteger . polyCoeffs LE $ ts !! n     | otherwise = error "t: negative index"  -- |Compute the coefficients of the n'th Chebyshev polynomial of the second kind.-u :: Num a => Int -> Poly a+u :: (Num a, Eq a) => Int -> Poly a u n | n >= 0    = poly LE . map fromInteger . polyCoeffs LE $ us !! n     | otherwise = error "u: negative index"
+ src/Math/Polynomial/Hermite.hs view
@@ -0,0 +1,48 @@+module Math.Polynomial.Hermite where++import Math.Polynomial+import Data.VectorSpace++probHermite :: [Poly Integer]+probHermite +    = one+    : [ multPoly x h_n ^-^ polyDeriv h_n+      | h_n <- probHermite+      ]++physHermite :: [Poly Integer]+physHermite+    = one+    : [ scalePoly 2 (multPoly x h_n) ^-^ polyDeriv h_n+      | h_n <- physHermite+      ]++evalProbHermite :: (Integral a, Num b) => a -> b -> b+evalProbHermite n = fst . evalProbHermiteDeriv n++evalProbHermiteDeriv :: (Integral a, Num b) => a -> b -> (b, b)+evalProbHermiteDeriv 0 _ = (1, 0)+evalProbHermiteDeriv 1 x = (x, 1)+evalProbHermiteDeriv n x+    | n < 0     = error "evalProbHermite: n < 0"+    | otherwise = loop 1 x 1+    where+        loop k h_k h_km1+            | k == n    = (h_k, k' * h_km1)+            | otherwise = loop (k+1) (x * h_k - k' * h_km1) h_k+            where k' = fromIntegral k++evalPhysHermite :: (Integral a, Num b) => a -> b -> b+evalPhysHermite n = fst . evalPhysHermiteDeriv n++evalPhysHermiteDeriv :: (Integral a, Num b) => a -> b -> (b,b)+evalPhysHermiteDeriv 0 _ = (1,   0)+evalPhysHermiteDeriv 1 x = (2*x, 2)+evalPhysHermiteDeriv n x+    | n < 0     = error "evalProbHermite: n < 0"+    | otherwise = loop 1 (2*x) 1+    where+        loop k h_k h_km1+            | k == n    = (h_k, 2 * k' * h_km1)+            | otherwise = loop (k+1) (2 * (x * h_k - k' * h_km1)) h_k+            where k' = fromIntegral k
src/Math/Polynomial/Interpolation.hs view
@@ -67,7 +67,7 @@ -- inherently ill-conditioned problem.  In most cases it is both faster and  -- more accurate to use 'polyInterp' or 'nevilleDiffs' instead of evaluating -- a fitted polynomial.-iterativePolyFit :: Fractional a => [(a,a)] -> Poly a+iterativePolyFit :: (Fractional a, Eq a) => [(a,a)] -> Poly a iterativePolyFit = poly LE . loop     where         loop  [] = []@@ -85,7 +85,7 @@ -- inherently ill-conditioned problem.  In most cases it is both faster and  -- more accurate to use 'polyInterp' or 'nevilleDiffs' instead of evaluating -- a fitted polynomial.-lagrangePolyFit :: Fractional a => [(a,a)] -> Poly a+lagrangePolyFit :: (Fractional a, Eq a) => [(a,a)] -> Poly a lagrangePolyFit xys = sumPolys     [ scalePoly f (fst (contractPoly p x))     | f <- zipWith (/) ys phis
src/Math/Polynomial/Lagrange.hs view
@@ -28,12 +28,12 @@ -- 'Math.Polynomial.Interpolation.polyInterp' evaluates the interpolating -- polynomial directly, and is both quicker and more stable than any method -- I know of that computes the coefficients.-lagrangeBasis :: Fractional a => [a] -> [Poly a]+lagrangeBasis :: (Fractional a, Eq a) => [a] -> [Poly a] lagrangeBasis xs =     [ foldl1 multPoly         [ if q /= 0             then poly LE [negate x_j/q, 1/q]-            else error ("lagrangeBasis: duplicate root: " ++ show x_i)+            else error ("lagrangeBasis: duplicate root")         | x_j <- otherXs         , let q = x_i - x_j         ]@@ -42,7 +42,7 @@  -- |Construct the Lagrange "master polynomial" for the Lagrange barycentric form: -- That is, the monic polynomial with a root at each point in the input list.-lagrange :: Num a => [a] -> Poly a+lagrange :: (Num a, Eq a) => [a] -> Poly a lagrange [] = one lagrange xs = foldl1 multPoly     [ poly LE [negate x_i, 1]
src/Math/Polynomial/Legendre.hs view
@@ -20,7 +20,7 @@     ]  -- |Compute the coefficients of the n'th Legendre polynomial.-legendre :: Fractional a => Int -> Poly a+legendre :: (Fractional a, Eq a) => Int -> Poly a legendre n = poly LE . map fromRational . polyCoeffs LE $ legendres !! n  -- |Evaluate the n'th Legendre polynomial at a point X.  Both more efficient
src/Math/Polynomial/Newton.hs view
@@ -6,7 +6,7 @@ -- |Returns the Newton basis set of polynomials associated with a set of  -- abscissas.  This is the set of monic polynomials each of which is @0@  -- at all previous points in the input list.-newtonBasis :: Num a => [a] -> [Poly a]+newtonBasis :: (Num a, Eq a) => [a] -> [Poly a] newtonBasis xs =      [ foldl multPoly (poly LE [1])          [ poly LE [-x_i, 1]
src/Math/Polynomial/NumInstance.hs view
@@ -9,7 +9,7 @@  import Math.Polynomial -instance Num a => Num (Poly a) where+instance (Num a, Eq a) => Num (Poly a) where     fromInteger i = poly LE [fromInteger i]     (+) = addPoly     negate = negatePoly
src/Math/Polynomial/Type.hs view
@@ -1,10 +1,36 @@-{-# LANGUAGE ViewPatterns, TypeFamilies #-}+{-# LANGUAGE ViewPatterns, TypeFamilies, GADTs #-} -- |Low-level interface for the 'Poly' type. module Math.Polynomial.Type      ( Endianness(..)-    , Poly, poly, polyCoeffs-    , trim, rawPoly, rawPolyCoeffs-    , polyIsZero, polyIsOne+    , Poly+    +    , zero+    +    , poly, polyN+    , unboxedPoly, unboxedPolyN+    +    , mapPoly+    +    , unboxPoly+    +    , rawListPoly+    , rawListPolyN+    , rawVectorPoly+    , rawUVectorPoly+    , trim+    +    , polyIsZero+    , polyIsOne+    +    , polyCoeffs+    , rawCoeffsOrder+    , rawPolyCoeffs+    , untrimmedPolyCoeffs+    +    , polyDegree+    , rawPolyDegree+    , rawPolyLength+         ) where  import Control.DeepSeq@@ -12,60 +38,8 @@ import Data.AdditiveGroup import Data.VectorSpace import Data.List.ZipSum--dropEnd :: (a -> Bool) -> [a] -> [a]--- dropEnd p = reverse . dropWhile p . reverse-dropEnd p = go id-    where-        go t (x:xs)-            -- if p x, stash x (will only be used if 'not (any p xs)')-            | p x       =        go (t.(x:))  xs-            -- otherwise insert x and all stashed values in output and reset the stash-            | otherwise = t (x : go  id       xs)-        -- at end of string discard the stash-        go _ [] = []---- |Trim zeroes from a polynomial (given a predicate for identifying zero).--- In particular, drops zeroes from the highest-order coefficients, so that--- @0x^n + 0x^(n-1) + 0x^(n-2) + ... + ax^k + ...@, @a /= 0@--- is normalized to @ax^k + ...@.  --- --- The 'Eq' instance for 'Poly' and all the standard constructors / destructors--- are defined using @trim (0==)@.-trim :: (a -> Bool) -> Poly a -> Poly a-trim      _ p@(Poly _ True _) = p-trim isZero   (Poly LE _ cs) = Poly LE True (dropEnd   isZero cs)-trim isZero   (Poly BE _ cs) = Poly BE True (dropWhile isZero cs)---- |Make a 'Poly' from a list of coefficients using the specified coefficient order.-poly :: Num a => Endianness -> [a] -> Poly a-poly end cs = trim (0==) (rawPoly end cs)---- |Make a 'Poly' from a list of coefficients using the specified coefficient order,--- without the 'Num' context (and therefore without trimming zeroes from the --- coefficient list)-rawPoly :: Endianness -> [a] -> Poly a-rawPoly end cs = Poly end False cs ---- |Get the coefficients of a a 'Poly' in the specified order.-polyCoeffs :: Num a => Endianness -> Poly a -> [a]-polyCoeffs end p = rawPolyCoeffs end (trim (0==) p)---- |Get the coefficients of a a 'Poly' in the specified order, without the 'Num'--- constraint (and therefore without trimming zeroes).--- --- This function does not respect the 'Eq' instance:---   @x == y@ =/=> @rawPolyCoeffs e x == rawPolyCoeffs e y@.-rawPolyCoeffs :: Endianness -> Poly a -> [a]-rawPolyCoeffs end (Poly e _ cs)-    | e == end  = cs-    | otherwise = reverse cs--polyIsZero :: Num a => Poly a -> Bool-polyIsZero = null . coeffs . trim (0==)--polyIsOne :: Num a => Poly a -> Bool-polyIsOne = ([1]==) . coeffs . trim (0==)+import qualified Data.Vector as V+import qualified Data.Vector.Unboxed as UV  data Endianness      = BE @@ -77,31 +51,47 @@ instance NFData Endianness where     rnf x = seq x () -data Poly a = Poly -    { endianness :: !Endianness-    , _trimmed   :: !Bool-    , coeffs     :: ![a]-    }+data Poly a where+    ListPoly ::+        { trimmed    :: !Bool+        , endianness :: !Endianness+        , listCoeffs :: ![a]+        } -> Poly a+    VectorPoly ::+        { trimmed    :: !Bool+        , endianness :: !Endianness+        , vCoeffs    :: !(V.Vector a)+        } -> Poly a+    UVectorPoly :: UV.Unbox a => +        { trimmed    :: !Bool+        , endianness :: !Endianness+        , uvCoeffs   :: !(UV.Vector a)+        } -> Poly a  instance NFData a => NFData (Poly a) where-    rnf (Poly e t c) = rnf e `seq` rnf t `seq` rnf c+    rnf (ListPoly    _ _ c) = rnf c+    rnf (VectorPoly  _ _ c) = V.foldr' seq () c+    rnf (UVectorPoly _ _ _) = () -instance Num a => Show (Poly a) where-    showsPrec p (trim (0==) -> Poly end _ cs) +instance Show a => Show (Poly a) where+    showsPrec p f         = showParen (p > 10)              ( showString "poly "-            . showsPrec 11 end+            . showsPrec 11 (rawCoeffsOrder f)             . showChar ' '-            . showsPrec 11 cs+            . showsPrec 11 (rawPolyCoeffs f)             ) +-- 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     p == q  -        | endianness p == endianness q-        = coeffs (trim (0==) p) == coeffs (trim (0==) q)+        | rawCoeffsOrder p == rawCoeffsOrder q+        =  rawPolyCoeffs (trim (0==) p) +        == rawPolyCoeffs (trim (0==) q)         | otherwise -        = polyCoeffs BE p == polyCoeffs BE q-        +        =  polyCoeffs LE p+        == polyCoeffs 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@@ -116,14 +106,134 @@ --             qCoeffs = polyCoeffs BE q  instance Functor Poly where-    fmap f (Poly end _ cs) = Poly end False (map f cs)+    fmap f (ListPoly    _ end cs) = ListPoly   False end (map f cs)+    fmap f (VectorPoly  _ end cs) = VectorPoly False end (V.map f cs)+    -- TODO: make sure this gets fused+    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)+ instance AdditiveGroup a => AdditiveGroup (Poly a) where-    zeroV = Poly LE True []-    (rawPolyCoeffs LE ->  a) ^+^ (rawPolyCoeffs LE ->  b) -        = Poly LE False (zipSumV a b)+    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     type Scalar (Poly a) = Scalar a     (*^) s = fmap (s *^)++-- |Trim zeroes from a polynomial (given a predicate for identifying zero).+-- In particular, drops zeroes from the highest-order coefficients, so that+-- @0x^n + 0x^(n-1) + 0x^(n-2) + ... + ax^k + ...@, @a /= 0@+-- is normalized to @ax^k + ...@.  +-- +-- The 'Eq' instance for 'Poly' and all the standard constructors / destructors+-- are defined using @trim (0==)@.+trim :: (a -> Bool) -> Poly a -> Poly a+trim _ p | trimmed p = p+trim isZero   (ListPoly    _ LE cs) = ListPoly    True LE (dropEnd   isZero cs)+trim isZero   (ListPoly    _ BE cs) = ListPoly    True BE (dropWhile isZero cs)+trim isZero   (VectorPoly  _ LE cs) = VectorPoly  True LE (V.reverse . V.dropWhile isZero . V.reverse $ cs)+trim isZero   (VectorPoly  _ BE cs) = VectorPoly  True BE (V.dropWhile isZero cs)+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)++-- |The polynomial \"0\"+zero :: Poly a+zero = ListPoly True LE []++-- |Make a 'Poly' from a list of coefficients using the specified coefficient order.+poly :: (Num a, Eq a) => Endianness -> [a] -> Poly a+poly end = trim (0==) . rawListPoly end++-- |Make a 'Poly' from a list of coefficients, at most 'n' of which are significant.+polyN :: (Num a, Eq a) => Int -> Endianness -> [a] -> Poly a+polyN n end = trim (0==) . rawVectorPoly end . V.fromListN n++unboxedPoly :: (UV.Unbox a, Num a, Eq a) => Endianness -> [a] -> Poly a+unboxedPoly end = trim (0==) . rawUVectorPoly end . UV.fromList++unboxedPolyN :: (UV.Unbox a, Num a, Eq a) => Int -> Endianness -> [a] -> Poly a+unboxedPolyN n end = trim (0==) . rawUVectorPoly end . UV.fromListN n++unboxPoly :: UV.Unbox a => Poly a -> Poly a+unboxPoly (ListPoly   t e cs) = UVectorPoly t e (UV.fromList cs)+unboxPoly (VectorPoly t e cs) = UVectorPoly t e (UV.fromListN (V.length cs) (V.toList cs))+unboxPoly p@UVectorPoly{} = p++-- |Make a 'Poly' from a list of coefficients using the specified coefficient order,+-- without the 'Num' context (and therefore without trimming zeroes from the +-- coefficient list)+rawListPoly :: Endianness -> [a] -> Poly a+rawListPoly = ListPoly False++rawListPolyN :: Int -> Endianness -> [a] -> Poly a+rawListPolyN n e = rawVectorPoly e . V.fromListN n++rawVectorPoly :: Endianness -> V.Vector a -> Poly a+rawVectorPoly = VectorPoly False++rawUVectorPoly :: UV.Unbox a => Endianness -> UV.Vector a -> Poly a+rawUVectorPoly = UVectorPoly False++-- |Get the degree of a a 'Poly' (the highest exponent with nonzero coefficient)+polyDegree :: (Num a, Eq a) => Poly a -> Int+polyDegree p = rawPolyDegree (trim (0==) p)++rawPolyDegree :: Poly a -> Int+rawPolyDegree p = rawPolyLength p - 1++rawPolyLength :: Poly a -> Int+rawPolyLength (ListPoly    _ _ cs) =    length cs+rawPolyLength (VectorPoly  _ _ cs) =  V.length cs+rawPolyLength (UVectorPoly _ _ cs) = UV.length cs+++-- |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)++polyIsZero :: (Num a, Eq a) => Poly a -> Bool+polyIsZero = null . rawPolyCoeffs . trim (0==)++polyIsOne :: (Num a, Eq a) => Poly a -> Bool+polyIsOne = ([1]==) . rawPolyCoeffs . trim (0==)++rawCoeffsOrder :: Poly a -> Endianness+rawCoeffsOrder = endianness++rawPolyCoeffs :: Poly a -> [a]+rawPolyCoeffs p@ListPoly{}         = listCoeffs p+rawPolyCoeffs p@VectorPoly{}       = V.toList (vCoeffs p)+rawPolyCoeffs p@UVectorPoly{}      = UV.toList (uvCoeffs p)++-- TODO: make sure (V.toList . V.reverse) gets fused+untrimmedPolyCoeffs :: Endianness -> Poly a -> [a]+untrimmedPolyCoeffs e1 (VectorPoly  _ e2 cs)+    | e1 == e2  = V.toList cs+    | otherwise = V.toList  (V.reverse cs)+untrimmedPolyCoeffs e1 (UVectorPoly _ e2 cs)+    | e1 == e2  = UV.toList cs+    | otherwise = UV.toList (UV.reverse cs)+untrimmedPolyCoeffs e1 (ListPoly _ e2 cs)+    | e1 == e2  = cs+    | otherwise = reverse cs++dropEnd :: (a -> Bool) -> [a] -> [a]+-- dropEnd p = reverse . dropWhile p . reverse+dropEnd p = go id+    where+        go t (x:xs)+            -- if p x, stash x (will only be used if 'not (any p xs)')+            | p x       =        go (t.(x:))  xs+            -- otherwise insert x and all stashed values in output and reset the stash+            | otherwise = t (x : go  id       xs)+        -- at end of string discard the stash+        go _ [] = []