diff --git a/polynomial.cabal b/polynomial.cabal --- a/polynomial.cabal +++ b/polynomial.cabal @@ -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 maintainer: James Cook 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 diff --git a/src/Math/Polynomial.hs b/src/Math/Polynomial.hs --- a/src/Math/Polynomial.hs +++ b/src/Math/Polynomial.hs @@ -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] diff --git a/src/Math/Polynomial/Chebyshev.hs b/src/Math/Polynomial/Chebyshev.hs --- a/src/Math/Polynomial/Chebyshev.hs +++ b/src/Math/Polynomial/Chebyshev.hs @@ -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" diff --git a/src/Math/Polynomial/Hermite.hs b/src/Math/Polynomial/Hermite.hs new file mode 100644 --- /dev/null +++ b/src/Math/Polynomial/Hermite.hs @@ -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 diff --git a/src/Math/Polynomial/Interpolation.hs b/src/Math/Polynomial/Interpolation.hs --- a/src/Math/Polynomial/Interpolation.hs +++ b/src/Math/Polynomial/Interpolation.hs @@ -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 diff --git a/src/Math/Polynomial/Lagrange.hs b/src/Math/Polynomial/Lagrange.hs --- a/src/Math/Polynomial/Lagrange.hs +++ b/src/Math/Polynomial/Lagrange.hs @@ -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] diff --git a/src/Math/Polynomial/Legendre.hs b/src/Math/Polynomial/Legendre.hs --- a/src/Math/Polynomial/Legendre.hs +++ b/src/Math/Polynomial/Legendre.hs @@ -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 diff --git a/src/Math/Polynomial/Newton.hs b/src/Math/Polynomial/Newton.hs --- a/src/Math/Polynomial/Newton.hs +++ b/src/Math/Polynomial/Newton.hs @@ -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] diff --git a/src/Math/Polynomial/NumInstance.hs b/src/Math/Polynomial/NumInstance.hs --- a/src/Math/Polynomial/NumInstance.hs +++ b/src/Math/Polynomial/NumInstance.hs @@ -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 diff --git a/src/Math/Polynomial/Type.hs b/src/Math/Polynomial/Type.hs --- a/src/Math/Polynomial/Type.hs +++ b/src/Math/Polynomial/Type.hs @@ -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 _ [] = []