poly 0.1.0.0 → 0.2.0.0
raw patch · 7 files changed
+518/−96 lines, 7 filesdep +primitivedep ~basedep ~semirings
Dependencies added: primitive
Dependency ranges changed: base, semirings
Files
- README.md +19/−2
- changelog.md +6/−0
- poly.cabal +7/−5
- src/Data/Poly.hs +15/−3
- src/Data/Poly/Semiring.hs +64/−0
- src/Data/Poly/Uni/Dense.hs +300/−74
- test/Main.hs +107/−12
README.md view
@@ -1,3 +1,20 @@-# poly+# poly [](https://travis-ci.org/Bodigrim/poly) [](https://hackage.haskell.org/package/poly) -A type to represent polynomials with Num and Semiring instances.+Polynomials with `Num` and `Semiring` instances, backed by `Vector`.++```haskell+> (X + 1) + (X - 1) :: VPoly Integer+2 * X + 0++> (X + 1) * (X - 1) :: UPoly Int+1 * X^2 + 0 * X + (-1)++> eval (X^2 + 1 :: UPoly Int) 3+10++> eval (X^2 + 1 :: VPoly (UPoly Int)) (X + 1)+1 * X^2 + 2 * X + 2++> deriv (X^3 + 3 * X) :: UPoly Int+3 * X^2 + 0 * X + 3+```
changelog.md view
@@ -1,3 +1,9 @@+# 0.2.0.0++* Fix a bug in `Num.(-)`.+* Add functions `constant`, `eval`, `deriv`, `integral`.+* Add a handy pattern synonym `X`.+ # 0.1.0.0 * Initial release.
poly.cabal view
@@ -1,8 +1,8 @@ name: poly-version: 0.1.0.0+version: 0.2.0.0 synopsis: Polynomials description:- A type to represent polynomials with Num and Semiring instances.+ Polynomials with `Num` and `Semiring` instances, backed by `Vector`. homepage: https://github.com/Bodigrim/poly#readme license: BSD3 license-file: LICENSE@@ -13,7 +13,7 @@ build-type: Simple extra-source-files: README.md cabal-version: >=1.10-tested-with: GHC ==7.10.3 GHC ==8.0.2 GHC ==8.2.2 GHC ==8.4.4 GHC ==8.6.4+tested-with: GHC ==8.0.2 GHC ==8.2.2 GHC ==8.4.4 GHC ==8.6.5 GHC ==8.8.1 extra-source-files: changelog.md @@ -25,10 +25,12 @@ hs-source-dirs: src exposed-modules: Data.Poly+ Data.Poly.Semiring other-modules: Data.Poly.Uni.Dense build-depends:- base >= 4.8 && < 5,+ base >= 4.9 && < 5,+ primitive, semirings, vector default-language: Haskell2010@@ -38,7 +40,7 @@ type: exitcode-stdio-1.0 main-is: Main.hs build-depends:- base >=4.8 && <5,+ base >=4.9 && <5, poly, QuickCheck >=2.10, quickcheck-classes >=0.6.1,
src/Data/Poly.hs view
@@ -4,11 +4,23 @@ -- Licence: BSD3 -- Maintainer: Andrew Lelechenko <andrew.lelechenko@gmail.com> ----- Polynomials.+-- Dense polynomials and a 'Num'-based interface. -- +{-# LANGUAGE PatternSynonyms #-}+ module Data.Poly- ( module Data.Poly.Uni.Dense+ ( Poly+ , VPoly+ , UPoly+ , unPoly+ -- * Num interface+ , toPoly+ , constant+ , pattern X+ , eval+ , deriv+ , integral ) where -import Data.Poly.Uni.Dense+import Data.Poly.Uni.Dense hiding (quotRem)
+ src/Data/Poly/Semiring.hs view
@@ -0,0 +1,64 @@+-- |+-- Module: Data.Poly.Semiring+-- Copyright: (c) 2019 Andrew Lelechenko+-- Licence: BSD3+-- Maintainer: Andrew Lelechenko <andrew.lelechenko@gmail.com>+--+-- Dense polynomials and a 'Semiring'-based interface.+--++{-# LANGUAGE PatternSynonyms #-}++module Data.Poly.Semiring+ ( Poly+ , VPoly+ , UPoly+ , unPoly+ -- * Semiring interface+ , toPoly+ , constant+ , pattern X+ , eval+ , deriv+ ) where++import Data.Semiring (Semiring)+import qualified Data.Vector.Generic as G++import Data.Poly.Uni.Dense (Poly(..), VPoly, UPoly)+import qualified Data.Poly.Uni.Dense as Dense++-- | Make 'Poly' from a vector of coefficients+-- (first element corresponds to a constant term).+--+-- >>> :set -XOverloadedLists+-- >>> toPoly [1,2,3] :: VPoly Integer+-- 3 * X^2 + 2 * X + 1+-- >>> toPoly [0,0,0] :: UPoly Int+-- 0+toPoly :: (Eq a, Semiring a, G.Vector v a) => v a -> Poly v a+toPoly = Dense.toPoly'++-- | Create a polynomial from a constant term.+constant :: (Eq a, Semiring a, G.Vector v a) => a -> Poly v a+constant = Dense.constant'++-- | Create an identity polynomial.+pattern X :: (Eq a, Semiring a, G.Vector v a, Eq (v a)) => Poly v a+pattern X = Dense.X'++-- | Evaluate at a given point.+--+-- >>> eval (X^2 + 1 :: UPoly Int) 3+-- 10+-- >>> eval (X^2 + 1 :: VPoly (UPoly Int)) (X + 1)+-- 1 * X^2 + 2 * X + 2+eval :: (Semiring a, G.Vector v a) => Poly v a -> a -> a+eval = Dense.eval'++-- | Take a derivative.+--+-- >>> deriv (X^3 + 3 * X) :: UPoly Int+-- 3 * X^2 + 0 * X + 3+deriv :: (Eq a, Semiring a, G.Vector v a) => Poly v a -> Poly v a+deriv = Dense.deriv'
src/Data/Poly/Uni/Dense.hs view
@@ -7,124 +7,350 @@ -- Dense polynomials of one variable. -- +{-# LANGUAGE PatternSynonyms #-} {-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE ViewPatterns #-} module Data.Poly.Uni.Dense ( Poly+ , VPoly+ , UPoly , unPoly+ -- * Num interface , toPoly+ , constant+ , pattern X+ , eval+ , deriv+ , integral+ , quotRem+ -- * Semiring interface , toPoly'+ , constant'+ , pattern X'+ , eval'+ , deriv' ) where -import Prelude hiding (negate)+import Prelude hiding (quotRem)+import Control.Exception import Control.Monad+import Control.Monad.Primitive import Control.Monad.ST-import Data.List (foldl')-import Data.Semiring (Semiring(..), Ring(..))-import Data.Vector (Vector)+import Data.List (foldl', intersperse)+import Data.Semigroup (stimes)+import Data.Semiring (Semiring(..), Add(..))+import qualified Data.Semiring as Semiring import qualified Data.Vector as V-import qualified Data.Vector.Mutable as MV+import qualified Data.Vector.Generic as G+import qualified Data.Vector.Generic.Mutable as MG+import qualified Data.Vector.Unboxed as U --- | Polynomials of one variable.+-- | Polynomials of one variable with coefficients from @a@,+-- backed by a 'G.Vector' @v@ (boxed, unboxed, storable, etc.). ----- >>> :set -XOverloadedLists--- >>> -- (1 + x) * (-1 + x) = (-1 + x^2)--- >>> toPoly [1,1] * toPoly [-1,1]--- Poly {unPoly = [-1,0,1]}+-- Use pattern 'X' for construction: ----- >>> :set -XOverloadedLists--- >>> -- (1 + x) + (1 - x) = 2--- >>> toPoly [1,1] + toPoly [1,-1]--- Poly {unPoly = [2]}-newtype Poly a = Poly- { unPoly :: Vector a+-- >>> (X + 1) + (X - 1) :: VPoly Integer+-- 2 * X + 0+-- >>> (X + 1) * (X - 1) :: UPoly Int+-- 1 * X^2 + 0 * X + (-1)+--+-- Polynomials are stored normalized, without leading+-- zero coefficients, so 0 * 'X' + 1 equals to 1.+--+-- 'Ord' instance does not make much sense mathematically,+-- it is defined only for the sake of 'Data.Set.Set', 'Data.Map.Map', etc.+--+newtype Poly v a = Poly+ { unPoly :: v a -- ^ Convert 'Poly' to a vector of coefficients -- (first element corresponds to a constant term). }- deriving (Eq, Ord, Show)+ deriving (Eq, Ord) +instance (Show a, G.Vector v a) => Show (Poly v a) where+ showsPrec d (Poly xs)+ | G.null xs+ = showString "0"+ | G.length xs == 1+ = showsPrec d (G.head xs)+ | otherwise+ = showParen (d > 0)+ $ foldl (.) id+ $ intersperse (showString " + ")+ $ G.ifoldl (\acc i c -> showCoeff i c : acc) [] xs+ where+ showCoeff 0 c = showsPrec 7 c+ showCoeff 1 c = showsPrec 7 c . showString " * X"+ showCoeff i c = showsPrec 7 c . showString " * X^" . showsPrec 7 i++-- | Polynomials backed by boxed vectors.+type VPoly = Poly V.Vector++-- | Polynomials backed by unboxed vectors.+type UPoly = Poly U.Vector+ -- | Make 'Poly' from a list of coefficients -- (first element corresponds to a constant term). -- -- >>> :set -XOverloadedLists--- >>> toPoly [1,2,3]--- Poly {unPoly = [1,2,3]}------ >>> :set -XOverloadedLists--- >>> toPoly [0,0,0]--- Poly {unPoly = []}-toPoly :: (Eq a, Num a) => Vector a -> Poly a+-- >>> toPoly [1,2,3] :: VPoly Integer+-- 3 * X^2 + 2 * X + 1+-- >>> toPoly [0,0,0] :: UPoly Int+-- 0+toPoly :: (Eq a, Num a, G.Vector v a) => v a -> Poly v a toPoly = Poly . dropWhileEnd (== 0) --- | Make 'Poly' from a vector of coefficients--- (first element corresponds to a constant term).------ >>> :set -XOverloadedLists--- >>> toPoly' [1,2,3]--- Poly {unPoly = [1,2,3]}------ >>> :set -XOverloadedLists--- >>> toPoly' [0,0,0]--- Poly {unPoly = []}-toPoly' :: (Eq a, Semiring a) => Vector a -> Poly a+toPoly' :: (Eq a, Semiring a, G.Vector v a) => v a -> Poly v a toPoly' = Poly . dropWhileEnd (== zero) -instance (Eq a, Num a) => Num (Poly a) where- Poly xs + Poly ys = toPoly $ zipOrCopy (+) xs ys- Poly xs - Poly ys = toPoly $ zipOrCopy (-) xs ys+instance (Eq a, Num a, G.Vector v a) => Num (Poly v a) where+ Poly xs + Poly ys = toPoly $ plusPoly (+) xs ys+ Poly xs - Poly ys = toPoly $ minusPoly negate (-) xs ys+ negate (Poly xs) = Poly $ G.map negate xs abs = id signum = const 1 fromInteger n = case fromInteger n of- 0 -> Poly $ V.empty- m -> Poly $ V.singleton m+ 0 -> Poly $ G.empty+ m -> Poly $ G.singleton m Poly xs * Poly ys = toPoly $ convolution 0 (+) (*) xs ys -instance (Eq a, Semiring a) => Semiring (Poly a) where- zero = Poly V.empty+instance (Eq a, Semiring a, G.Vector v a) => Semiring (Poly v a) where+ zero = Poly G.empty one | (one :: a) == zero = zero- | otherwise = Poly $ V.singleton one- plus (Poly xs) (Poly ys) = toPoly' $ zipOrCopy plus xs ys+ | otherwise = Poly $ G.singleton one+ plus (Poly xs) (Poly ys) = toPoly' $ plusPoly plus xs ys times (Poly xs) (Poly ys) = toPoly' $ convolution zero plus times xs ys -instance (Eq a, Ring a) => Ring (Poly a) where- negate (Poly xs) = Poly $ V.map negate xs+instance (Eq a, Semiring.Ring a, G.Vector v a) => Semiring.Ring (Poly v a) where+ negate (Poly xs) = Poly $ G.map Semiring.negate xs -dropWhileEnd :: (a -> Bool) -> Vector a -> Vector a-dropWhileEnd p xs = V.slice 0 (go (V.length xs)) xs+dropWhileEnd+ :: G.Vector v a+ => (a -> Bool)+ -> v a+ -> v a+dropWhileEnd p xs = G.basicUnsafeSlice 0 (go (G.basicLength xs)) xs where go 0 = 0- go n = if p (xs V.! (n - 1)) then go (n - 1) else n+ go n = if p (G.unsafeIndex xs (n - 1)) then go (n - 1) else n -zipOrCopy :: (a -> a -> a) -> Vector a -> Vector a -> Vector a-zipOrCopy f xs ys = runST $ do- zs <- MV.new lenZs- forM_ [0 .. lenZs - 1] $ \i ->- MV.write zs i (f (xs V.! i) (ys V.! i))- when (lenXs < lenYs) $- forM_ [lenXs .. lenYs - 1] $ \i ->- MV.write zs i (ys V.! i)- when (lenYs < lenXs) $- forM_ [lenYs .. lenXs - 1] $ \i ->- MV.write zs i (xs V.! i)- V.unsafeFreeze zs- where- lenXs = V.length xs- lenYs = V.length ys- lenZs = lenXs `max` lenYs+plusPoly+ :: G.Vector v a+ => (a -> a -> a)+ -> v a+ -> v a+ -> v a+plusPoly add xs ys = runST $ do+ zs <- MG.new (G.basicLength xs `max` G.basicLength ys)+ plusPolyM add xs ys zs+ G.unsafeFreeze zs -convolution :: a -> (a -> a -> a) -> (a -> a -> a) -> Vector a -> Vector a -> Vector a+plusPolyM+ :: (PrimMonad m, G.Vector v a)+ => (a -> a -> a)+ -> v a+ -> v a+ -> G.Mutable v (PrimState m) a+ -> m ()+plusPolyM add xs ys zs = do+ let lenXs = G.basicLength xs+ lenYs = G.basicLength ys+ case lenXs `compare` lenYs of+ LT -> do+ forM_ [0 .. lenXs - 1] $ \i ->+ MG.unsafeWrite zs i (add (G.unsafeIndex xs i) (G.unsafeIndex ys i))+ G.unsafeCopy+ (MG.basicUnsafeSlice lenXs (lenYs - lenXs) zs)+ (G.basicUnsafeSlice lenXs (lenYs - lenXs) ys)+ EQ -> do+ forM_ [0 .. lenXs - 1] $ \i ->+ MG.unsafeWrite zs i (add (G.unsafeIndex xs i) (G.unsafeIndex ys i))+ GT -> do+ forM_ [0 .. lenYs - 1] $ \i ->+ MG.unsafeWrite zs i (add (G.unsafeIndex xs i) (G.unsafeIndex ys i))+ G.unsafeCopy+ (MG.basicUnsafeSlice lenYs (lenXs - lenYs) zs)+ (G.basicUnsafeSlice lenYs (lenXs - lenYs) xs)++minusPoly+ :: G.Vector v a+ => (a -> a)+ -> (a -> a -> a)+ -> v a+ -> v a+ -> v a+minusPoly neg sub xs ys = runST $ do+ zs <- MG.new (G.basicLength xs `max` G.basicLength ys)+ minusPolyM neg sub xs ys zs+ G.unsafeFreeze zs++minusPolyM+ :: (PrimMonad m, G.Vector v a)+ => (a -> a)+ -> (a -> a -> a)+ -> v a+ -> v a+ -> G.Mutable v (PrimState m) a+ -> m ()+minusPolyM neg sub xs ys zs = do+ let lenXs = G.basicLength xs+ lenYs = G.basicLength ys+ case lenXs `compare` lenYs of+ LT -> do+ forM_ [0 .. lenXs - 1] $ \i ->+ MG.unsafeWrite zs i (sub (G.unsafeIndex xs i) (G.unsafeIndex ys i))+ forM_ [lenXs .. lenYs - 1] $ \i ->+ MG.unsafeWrite zs i (neg (G.unsafeIndex ys i))+ EQ -> do+ forM_ [0 .. lenXs - 1] $ \i ->+ MG.unsafeWrite zs i (sub (G.unsafeIndex xs i) (G.unsafeIndex ys i))+ GT -> do+ forM_ [0 .. lenYs - 1] $ \i ->+ MG.unsafeWrite zs i (sub (G.unsafeIndex xs i) (G.unsafeIndex ys i))+ G.unsafeCopy+ (MG.basicUnsafeSlice lenYs (lenXs - lenYs) zs)+ (G.basicUnsafeSlice lenYs (lenXs - lenYs) xs)++convolution+ :: G.Vector v a+ => a+ -> (a -> a -> a)+ -> (a -> a -> a)+ -> v a+ -> v a+ -> v a convolution zer add mul xs ys- | V.null xs || V.null ys = V.empty+ | G.null xs || G.null ys = G.empty | otherwise = runST $ do- zs <- MV.new lenZs+ zs <- MG.new lenZs forM_ [0 .. lenZs - 1] $ \k -> do let is = [max (k - lenYs + 1) 0 .. min k (lenXs - 1)]- -- js = reverse [max (k - lenXs) 0 .. min k lenYs]- let acc = foldl' add zer $ flip map is $ \i -> mul (xs V.! i) (ys V.! (k - i))- MV.write zs k acc- V.unsafeFreeze zs+ acc = foldl' add zer $ flip map is $ \i ->+ mul (G.unsafeIndex xs i) (G.unsafeIndex ys (k - i))+ MG.unsafeWrite zs k acc+ G.unsafeFreeze zs where- lenXs = V.length xs- lenYs = V.length ys+ lenXs = G.basicLength xs+ lenYs = G.basicLength ys lenZs = lenXs + lenYs - 1++-- | This is just a proof of concept,+-- which should be replaced by a proper 'Euclidean' interface.+quotRem+ :: (Integral a, G.Vector v a)+ => Poly v a+ -> Poly v a+ -> (Poly v a, Poly v a)+quotRem (Poly xs) (Poly ys) = (toPoly qs, toPoly rs)+ where+ (qs, rs) = quotRem' xs ys++quotRem'+ :: (Integral a, G.Vector v a)+ => v a+ -> v a+ -> (v a, v a)+quotRem' xs ys+ | G.null ys = throw DivideByZero+ | G.basicLength xs < G.basicLength ys = (G.empty, xs)+ | otherwise = runST $ do+ let lenXs = G.basicLength xs+ lenYs = G.basicLength ys+ lenQs = lenXs - lenYs + 1+ qs <- MG.new lenQs+ rs <- MG.new lenXs+ G.unsafeCopy rs xs+ forM_ [lenQs - 1, lenQs - 2 .. 0] $ \i -> do+ let j = lenXs - 1 + i - (lenQs - 1)+ r <- MG.unsafeRead rs j+ let q = r `quot` G.unsafeLast ys+ MG.unsafeWrite qs i q+ forM_ [0 .. lenYs - 1] $ \k -> do+ MG.unsafeModify rs (\c -> c - q * G.unsafeIndex ys k) (j + k - lenYs + 1)+ (,) <$> G.unsafeFreeze qs <*> G.unsafeFreeze rs+++-- | Create a polynomial from a constant term.+constant :: (Eq a, Num a, G.Vector v a) => a -> Poly v a+constant 0 = Poly G.empty+constant c = Poly $ G.singleton c++constant' :: (Eq a, Semiring a, G.Vector v a) => a -> Poly v a+constant' c+ | c == zero = Poly G.empty+ | otherwise = Poly $ G.singleton c++data StrictPair a b = !a :*: !b++infixr 1 :*:++fst' :: StrictPair a b -> a+fst' (a :*: _) = a++-- | Evaluate at a given point.+--+-- >>> eval (X^2 + 1 :: UPoly Int) 3+-- 10+-- >>> eval (X^2 + 1 :: VPoly (UPoly Int)) (X + 1)+-- 1 * X^2 + 2 * X + 2+eval :: (Num a, G.Vector v a) => Poly v a -> a -> a+eval (Poly cs) x = fst' $+ G.foldl' (\(acc :*: xn) cn -> (acc + cn * xn :*: x * xn)) (0 :*: 1) cs++eval' :: (Semiring a, G.Vector v a) => Poly v a -> a -> a+eval' (Poly cs) x = fst' $+ G.foldl' (\(acc :*: xn) cn -> (acc `plus` cn `times` xn :*: x `times` xn)) (zero :*: one) cs++-- | Take a derivative.+--+-- >>> deriv (X^3 + 3 * X) :: UPoly Int+-- 3 * X^2 + 0 * X + 3+deriv :: (Eq a, Num a, G.Vector v a) => Poly v a -> Poly v a+deriv (Poly xs)+ | G.null xs = Poly G.empty+ | otherwise = toPoly $ G.imap (\i x -> fromIntegral (i + 1) * x) $ G.tail xs++deriv' :: (Eq a, Semiring a, G.Vector v a) => Poly v a -> Poly v a+deriv' (Poly xs)+ | G.null xs = Poly G.empty+ | otherwise = toPoly' $ G.imap (\i x -> getAdd (stimes (i + 1) (Add x))) $ G.tail xs++-- | Compute an indefinite integral of a polynomial,+-- setting constant term to zero.+--+-- >>> integral (constant 3.0 * X^2 + constant 3.0) :: UPoly Double+-- 1.0 * X^3 + 0.0 * X^2 + 3.0 * X + 0.0+integral :: (Eq a, Fractional a, G.Vector v a) => Poly v a -> Poly v a+integral (Poly xs)+ | G.null xs = Poly G.empty+ | otherwise = toPoly $ runST $ do+ zs <- MG.new (lenXs + 1)+ MG.unsafeWrite zs 0 0+ forM_ [0 .. lenXs - 1] $ \i ->+ MG.unsafeWrite zs (i + 1) (G.unsafeIndex xs i * recip (fromIntegral i + 1))+ G.unsafeFreeze zs+ where+ lenXs = G.basicLength xs++-- | Create an identity polynomial.+pattern X :: (Eq a, Num a, G.Vector v a, Eq (v a)) => Poly v a+pattern X <- ((==) var -> True)+ where X = var++var :: forall a v. (Eq a, Num a, G.Vector v a, Eq (v a)) => Poly v a+var+ | (1 :: a) == 0 = Poly G.empty+ | otherwise = Poly $ G.fromList [0, 1]++-- | Create an identity polynomial.+pattern X' :: (Eq a, Semiring a, G.Vector v a, Eq (v a)) => Poly v a+pattern X' <- ((==) var' -> True)+ where X' = var'++var' :: forall a v. (Eq a, Semiring a, G.Vector v a, Eq (v a)) => Poly v a+var'+ | (one :: a) == zero = Poly G.empty+ | otherwise = Poly $ G.fromList [zero, one]
test/Main.hs view
@@ -1,19 +1,33 @@+{-# LANGUAGE ScopedTypeVariables #-}+ {-# OPTIONS_GHC -fno-warn-orphans #-} module Main where +import Prelude hiding (quotRem) import Data.Int import Data.Poly+import qualified Data.Poly.Semiring as S import Data.Proxy-import Data.Semiring+import Data.Semiring (Semiring) import qualified Data.Vector as V+import qualified Data.Vector.Generic as G+import qualified Data.Vector.Unboxed as U import Test.Tasty import Test.Tasty.QuickCheck-import Test.QuickCheck.Classes+import Test.QuickCheck.Classes (lawsProperties, semiringLaws, ringLaws) +instance (Eq a, Semiring a, Arbitrary a, G.Vector v a) => Arbitrary (Poly v a) where+ arbitrary = S.toPoly . G.fromList <$> arbitrary+ shrink = fmap (S.toPoly . G.fromList) . shrink . G.toList . unPoly+ main :: IO () main = defaultMain $ testGroup "All"- [ semiringTests+ [ arithmeticTests+ , semiringTests+ , evalTests+ , derivTests+ , quotRemTests ] semiringTests :: TestTree@@ -21,14 +35,95 @@ = testGroup "Semiring" $ map (uncurry testProperty) $ concatMap lawsProperties- [ semiringLaws (Proxy :: Proxy (Poly ()))- , ringLaws (Proxy :: Proxy (Poly ()))- , semiringLaws (Proxy :: Proxy (Poly Int8))- , ringLaws (Proxy :: Proxy (Poly Int8))- , semiringLaws (Proxy :: Proxy (Poly Integer))- , ringLaws (Proxy :: Proxy (Poly Integer))+ [ semiringLaws (Proxy :: Proxy (Poly U.Vector ()))+ , ringLaws (Proxy :: Proxy (Poly U.Vector ()))+ , semiringLaws (Proxy :: Proxy (Poly U.Vector Int8))+ , ringLaws (Proxy :: Proxy (Poly U.Vector Int8))+ , semiringLaws (Proxy :: Proxy (Poly V.Vector Integer))+ , ringLaws (Proxy :: Proxy (Poly V.Vector Integer)) ] -instance (Eq a, Semiring a, Arbitrary a) => Arbitrary (Poly a) where- arbitrary = toPoly' . V.fromList <$> arbitrary- shrink = fmap (toPoly' . V.fromList) . shrink . V.toList . unPoly+arithmeticTests :: TestTree+arithmeticTests = testGroup "Arithmetic"+ [ testProperty "addition matches reference" $+ \(xs :: [Int]) ys -> toPoly (V.fromList (addRef xs ys)) ===+ toPoly (V.fromList xs) + toPoly (V.fromList ys)+ , testProperty "subtraction matches reference" $+ \(xs :: [Int]) ys -> toPoly (V.fromList (subRef xs ys)) ===+ toPoly (V.fromList xs) - toPoly (V.fromList ys)+ ]++addRef :: Num a => [a] -> [a] -> [a]+addRef [] ys = ys+addRef xs [] = xs+addRef (x : xs) (y : ys) = (x + y) : addRef xs ys++subRef :: Num a => [a] -> [a] -> [a]+subRef [] ys = map negate ys+subRef xs [] = xs+subRef (x : xs) (y : ys) = (x - y) : subRef xs ys++evalTests :: TestTree+evalTests = testGroup "eval" $ concat+ [ evalTestGroup (Proxy :: Proxy (Poly U.Vector Int8))+ , evalTestGroup (Proxy :: Proxy (Poly V.Vector Integer))+ ]++evalTestGroup+ :: forall v a.+ (Eq a, Num a, Semiring a, Arbitrary a, Show a, Eq (v a), Show (v a), G.Vector v a)+ => Proxy (Poly v a)+ -> [TestTree]+evalTestGroup _ =+ [ testProperty "eval (p + q) r = eval p r + eval q r" $+ \p q r -> e (p + q) r === e p r + e q r+ , testProperty "eval (p * q) r = eval p r * eval q r" $+ \p q r -> e (p * q) r === e p r * e q r+ , testProperty "eval x p = p" $+ \p -> e X p === p+ , testProperty "eval (constant c) p = c" $+ \c p -> e (constant c) p === c++ , testProperty "eval' (p + q) r = eval' p r + eval' q r" $+ \p q r -> e' (p + q) r === e' p r + e' q r+ , testProperty "eval' (p * q) r = eval' p r * eval' q r" $+ \p q r -> e' (p * q) r === e' p r * e' q r+ , testProperty "eval' x p = p" $+ \p -> e' S.X p === p+ , testProperty "eval' (S.constant c) p = c" $+ \c p -> e' (S.constant c) p === c+ ]++ where+ e :: Poly v a -> a -> a+ e = eval+ e' :: Poly v a -> a -> a+ e' = S.eval++derivTests :: TestTree+derivTests = testGroup "deriv"+ [ testProperty "deriv = S.deriv" $+ \(p :: Poly V.Vector Integer) -> deriv p === S.deriv p+ , testProperty "deriv . integral = id" $+ \(p :: Poly V.Vector Rational) -> deriv (integral p) === p+ , testProperty "deriv c = 0" $+ \c -> deriv (constant c :: Poly V.Vector Int) === 0+ , testProperty "deriv cX = c" $+ \c -> deriv (constant c * X :: Poly V.Vector Int) === constant c+ , testProperty "deriv (p + q) = deriv p + deriv q" $+ \p q -> deriv (p + q) === (deriv p + deriv q :: Poly V.Vector Int)+ , testProperty "deriv (p * q) = p * deriv q + q * deriv p" $+ \p q -> deriv (p * q) === (p * deriv q + q * deriv p :: Poly V.Vector Int)+ -- The following property takes too long for a regular test-suite+ -- , testProperty "deriv (eval p q) = deriv q * eval (deriv p) q" $+ -- \(p :: Poly V.Vector Int) (q :: Poly U.Vector Int) ->+ -- deriv (eval (toPoly $ fmap constant $ unPoly p) q) ===+ -- deriv q * eval (toPoly $ fmap constant $ unPoly $ deriv p) q+ ]++quotRemTests :: TestTree+quotRemTests = testGroup "quotRem" []+ -- [ testProperty "(q, r) = x `quotRem` y ==> q * y + r == x" $+ -- \(x :: Poly U.Vector Int) y -> let (q, r) = x `quotRem` y in+ -- y === 0 .||. q * y + r === x+ -- ]