polynomial 0.5 → 0.6

raw patch · 5 files changed

+163/−53 lines, 5 filesdep +deepseqPVP ok

version bump matches the API change (PVP)

Dependencies added: deepseq

API changes (from Hackage documentation)

Files

• polynomial.cabal +2/−2
• src/Math/Polynomial.hs +97/−13
• src/Math/Polynomial/Bernstein.hs +3/−8
• src/Math/Polynomial/Chebyshev.hs +17/−10
• src/Math/Polynomial/Type.hs +44/−20

@@ -1,5 +1,5 @@ name:                   polynomial-version:                0.5+version:                0.6 stability:              provisional  cabal-version:          >= 1.6@@ -35,4 +35,4 @@                         Math.Polynomial.Type                         Math.Polynomial.Pretty                         -  build-depends:        base >= 3 && <5, pretty, prettyclass, vector-space+  build-depends:        base >= 3 && <5, deepseq, pretty, prettyclass, vector-space

@@ -3,8 +3,9 @@ module Math.Polynomial     ( Endianness(..)     , Poly, poly, polyCoeffs, polyIsZero, polyIsOne-    , zero, one, x+    , zero, one, constPoly, x     , scalePoly, negatePoly+    , composePoly     , addPoly, sumPolys, multPoly, powPoly     , quotRemPoly, quotPoly, remPoly     , evalPoly, evalPolyDeriv, evalPolyDerivs@@ -19,21 +20,40 @@ 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 = poly LE [1]+one = constPoly 1 +-- |The polynomial (in x) \"x\" x :: Num a => Poly a x = poly 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]++-- |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 0 _ = zero scalePoly s p = fmap (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 +-- |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) @@ -44,24 +64,38 @@ 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 $ multX ys+multPoly (polyCoeffs LE -> xs) (polyCoeffs LE -> ys) = poly LE (multPolyLE xs ys)++-- |(Internal): multiply polynomials in LE order.  O(length xs * length ys).+multPolyLE :: Num a => [a] -> [a] -> [a]+multPolyLE _  []     = []+multPolyLE xs (y:ys) = foldr mul [] xs     where-        multX (0:ys) = 0:multX ys-        multX ys = foldl zipSum []-            [ shift ++ map (x *) ys-            | (x, shift) <- zip xs (inits (repeat 0))-            , x /= 0-            ]+        mul 0 bs = 0 : bs+        mul x bs = (x*y) : zipSum (map (x*) ys) bs +-- |Given a polynomial 'f' and exponent 'n', computes the polynomial 'g'+-- 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 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)) +-- |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 => 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)     where@@ -75,8 +109,11 @@                 u' = tail (zipSum u (map (negate q0 *) v))  quotPoly :: Fractional a => Poly a -> Poly a -> Poly a-quotPoly u v = fst (quotRemPoly u v)+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 _ b | polyIsZero b = error "remPoly: divide by zero" remPoly (polyCoeffs BE -> u) (polyCoeffs BE -> v)     = go u (length u - length v)     where@@ -89,17 +126,59 @@                 q0 = head u / v0                 u' = tail (zipSum u (map (negate 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.+composePoly :: Num 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 => a -> [a] -> [a]+addScalarLE 0 bs = 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 :: Num 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 +-- |Evaluate a polynomial and its derivative (respectively) at a point. evalPolyDeriv :: Num 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) +-- |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 => Poly a -> a -> [a] evalPolyDerivs (polyCoeffs LE -> cs) x = trunc . zipWith (*) factorials $ foldr mul [] cs     where@@ -117,6 +196,9 @@         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 a b      | polyIsZero b  = if polyIsZero a@@ -124,13 +206,13 @@         else monic a     | otherwise     = gcdPoly b (a `remPoly` b) --- |Normalize a polynomial so that its highest-order coefficient is 1+-- |(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) -+-- |Compute the derivative of a polynomial. polyDeriv :: Num a => Poly a -> Poly a polyDeriv (polyCoeffs LE -> cs) = poly LE     [ c * n@@ -138,6 +220,7 @@     | n <- iterate (1+) 1     ] +-- |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 :     [ c / n@@ -145,7 +228,7 @@     | n <- iterate (1+) 1     ] --- |Separate a polynomial into a set of factors none of which have+-- |Separate a nonzero polynomial into a set of factors none of which have -- multiple roots, and the product of which is the original polynomial. -- Note that if division is not exact, it may fail to separate roots.   -- Rational coefficients is a good idea.@@ -153,6 +236,7 @@ -- Useful when applicable as a way to simplify root-finding problems. separateRoots :: Fractional a => Poly a -> [Poly a] separateRoots p+    | polyIsZero q  = error "separateRoots: zero polynomial"     | polyIsOne q   = [p]     | otherwise     = p `quotPoly` q : separateRoots q     where

@@ -58,14 +58,9 @@ -- |de Casteljau's algorithm, returning the whole tableau.  Used both for -- evaluating and splitting polynomials in Bernstein form. deCasteljau :: Num a => [a] -> a -> [[a]]-deCasteljau cs t = takeWhile (not.null) table-    where-        table = cs : -            [ [ b_i * (1-t) + b_ip1 * t-              | b_i:b_ip1:_ <- tails row-              ]-            | row <- table-            ]+deCasteljau [] _ = []+deCasteljau cs t = cs : deCasteljau (zipWith (interp t) cs (tail cs)) t+    where interp t x0 x1 = (1-t)*x0 + t*x1  -- |Given a polynomial in Bernstein form (that is, a list of coefficients -- for a basis set from 'bernstein', such as is returned by 'bernsteinFit')

@@ -1,4 +1,4 @@-{-# LANGUAGE ParallelListComp #-}+{-# LANGUAGE ParallelListComp, BangPatterns #-} module Math.Polynomial.Chebyshev where  import Math.Polynomial@@ -7,7 +7,7 @@ -- |The Chebyshev polynomials of the first kind with 'Integer' coefficients. ts :: [Poly Integer] ts = poly LE [1] : -    [ addPoly (poly LE [0, 1]    `multPoly` t_n)+    [ addPoly (x                `multPoly` t_n)               (poly LE [-1,0,1] `multPoly` u_n)     | t_n <- ts     | u_n <- poly LE [0] : us@@ -16,24 +16,26 @@ -- The Chebyshev polynomials of the second kind with 'Integer' coefficients. us :: [Poly Integer] us = -    [ addPoly t_n (multPoly u_n (poly LE [0,1]))+    [ addPoly t_n (multPoly x u_n)     | t_n <- ts     | u_n <- poly LE [0] : us     ]  -- |Compute the coefficients of the n'th Chebyshev polynomial of the first kind. t :: Num a => Int -> Poly a-t n = poly LE . map fromInteger . polyCoeffs LE $ ts !! n+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 n = poly LE . map fromInteger . polyCoeffs LE $ us !! n+u n | n >= 0    = poly LE . map fromInteger . polyCoeffs LE $ us !! n+    | otherwise = error "u: negative index"  -- |Evaluate the n'th Chebyshev polynomial of the first kind at a point X.   -- Both more efficient and more numerically stable than computing the  -- coefficients and evaluating the polynomial. evalT :: Num a => Int -> a -> a-evalT n x = evalTs x !! n+evalT n x = fst (evalTU n x)  -- |Evaluate all the Chebyshev polynomials of the first kind at a point X. evalTs :: Num a => a -> [a]@@ -43,7 +45,7 @@ -- Both more efficient and more numerically stable than computing the  -- coefficients and evaluating the polynomial. evalU :: Num a => Int -> a -> a-evalU n x = evalUs x !! n+evalU n x = snd (evalTU n x)  -- |Evaluate all the Chebyshev polynomials of the second kind at a point X. evalUs :: Num a => a -> [a]@@ -51,10 +53,15 @@  -- |Evaluate the n'th Chebyshev polynomials of both kinds at a point X. evalTU :: Num a => Int -> a -> (a,a)-evalTU n x = (ts!!n, us!!n)-    where (ts,us) = evalTsUs x+evalTU n x = go n 1 0+    where+        go !0 !t_n !u_n = (t_n, u_n)+        go !n !t_n !u_n = go (n-1) t_np1 u_np1+            where+                t_np1 = x * t_n - (1-x*x)*u_n+                u_np1 = x * u_n + t_n --- |Evaluate all the Chebyshev polynomials of the both kinds at a point X.+-- |Evaluate all the Chebyshev polynomials of both kinds at a point X. evalTsUs :: Num a => a -> ([a], [a]) evalTsUs x = (ts, tail us)     where

@@ -1,10 +1,13 @@ {-# LANGUAGE ViewPatterns, TypeFamilies #-}+-- |Low-level interface for the 'Poly' type. module Math.Polynomial.Type      ( Endianness(..)     , Poly, poly, polyCoeffs+    , trim, rawPoly, rawPolyCoeffs     , polyIsZero, polyIsOne     ) where +import Control.DeepSeq -- import Data.List.Extras.LazyLength import Data.AdditiveGroup import Data.VectorSpace@@ -22,26 +25,47 @@         -- at end of string discard the stash         go _ [] = [] -trim :: Num a => Poly a -> Poly a-trim p@(Poly _ True _) = p-trim   (Poly LE _ cs) = Poly LE True (dropEnd   (==0) cs)-trim   (Poly BE _ cs) = Poly BE True (dropWhile (==0) cs)+-- |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 (Poly end False cs)+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 = case trim p of-    Poly e _ cs | e == end  -> cs-                | otherwise -> reverse cs+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+polyIsZero = null . coeffs . trim (0==)  polyIsOne :: Num a => Poly a -> Bool-polyIsOne = ([1]==) . coeffs . trim+polyIsOne = ([1]==) . coeffs . trim (0==)  data Endianness      = BE @@ -50,14 +74,20 @@     -- ^ Little-Endian (head is const term)     deriving (Eq, Ord, Enum, Bounded, Show) +instance NFData Endianness where+    rnf x = seq x ()+ data Poly a = Poly      { endianness :: !Endianness     , _trimmed   :: !Bool     , coeffs     :: ![a]     } +instance NFData a => NFData (Poly a) where+    rnf (Poly e t c) = rnf e `seq` rnf t `seq` rnf c+ instance Num a => Show (Poly a) where-    showsPrec p (trim -> Poly end _ cs) +    showsPrec p (trim (0==) -> Poly end _ cs)          = showParen (p > 10)              ( showString "poly "             . showsPrec 11 end@@ -68,7 +98,7 @@ instance (Num a, Eq a) => Eq (Poly a) where     p == q           | endianness p == endianness q-        = coeffs (trim p) == coeffs (trim q)+        = coeffs (trim (0==) p) == coeffs (trim (0==) q)         | otherwise          = polyCoeffs BE p == polyCoeffs BE q         @@ -88,16 +118,10 @@ instance Functor Poly where     fmap f (Poly end _ cs) = Poly end False (map f cs) ---- Local-use-only: extract coefficients in LE order, without Num constraint--- (and therefore without trimming)-le :: Poly a -> [a]-le p@(endianness -> LE) = coeffs p-le p                    = reverse (coeffs p)- instance AdditiveGroup a => AdditiveGroup (Poly a) where     zeroV = Poly LE True []-    (le ->  a) ^+^ (le ->  b) = Poly LE False (zipSumV a b)+    (rawPolyCoeffs LE ->  a) ^+^ (rawPolyCoeffs LE ->  b) +        = Poly LE False (zipSumV a b)     negateV = fmap negateV  instance VectorSpace a => VectorSpace (Poly a) where