diff --git a/polynomial.cabal b/polynomial.cabal --- a/polynomial.cabal +++ b/polynomial.cabal @@ -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 diff --git a/src/Data/VectorSpace/WrappedNum.hs b/src/Data/VectorSpace/WrappedNum.hs new file mode 100644 --- /dev/null +++ b/src/Data/VectorSpace/WrappedNum.hs @@ -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 + (*^) = (*) diff --git a/src/Math/Polynomial.hs b/src/Math/Polynomial.hs --- a/src/Math/Polynomial.hs +++ b/src/Math/Polynomial.hs @@ -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. diff --git a/src/Math/Polynomial/Bernoulli.hs b/src/Math/Polynomial/Bernoulli.hs new file mode 100644 --- /dev/null +++ b/src/Math/Polynomial/Bernoulli.hs @@ -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) 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,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==) diff --git a/src/Math/Polynomial/VectorSpace.hs b/src/Math/Polynomial/VectorSpace.hs new file mode 100644 --- /dev/null +++ b/src/Math/Polynomial/VectorSpace.hs @@ -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 + ]