poly 0.3.0.0 → 0.3.1.0
raw patch · 16 files changed
+432/−59 lines, 16 filesdep ~QuickCheckdep ~basedep ~gauge
Dependency ranges changed: QuickCheck, base, gauge, primitive, quickcheck-classes, semirings, tasty, tasty-quickcheck, vector
Files
- README.md +35/−13
- bench/DenseBench.hs +16/−3
- changelog.md +5/−0
- poly.cabal +16/−15
- src/Data/Poly.hs +5/−0
- src/Data/Poly/Internal/Dense.hs +90/−6
- src/Data/Poly/Internal/Dense/Fractional.hs +9/−0
- src/Data/Poly/Internal/Dense/GcdDomain.hs +5/−0
- src/Data/Poly/Internal/PolyOverFractional.hs +9/−0
- src/Data/Poly/Internal/Sparse.hs +31/−4
- src/Data/Poly/Internal/Sparse/Fractional.hs +9/−0
- src/Data/Poly/Internal/Sparse/GcdDomain.hs +5/−0
- src/Data/Poly/Semiring.hs +5/−0
- test/Dense.hs +43/−9
- test/Quaternion.hs +113/−0
- test/Sparse.hs +36/−9
README.md view
@@ -1,7 +1,9 @@-# poly [](https://travis-ci.org/Bodigrim/poly) [](https://hackage.haskell.org/package/poly)+# poly [](https://travis-ci.org/Bodigrim/poly) [](https://hackage.haskell.org/package/poly) [](https://matrix.hackage.haskell.org/package/poly) [](http://stackage.org/lts/package/poly) [](http://stackage.org/nightly/package/poly) -Univariate polynomials, backed by `Vector`. ++Haskell library for univariate polynomials, backed by `Vector`.+ ```haskell > (X + 1) + (X - 1) :: VPoly Integer 2 * X + 0@@ -26,7 +28,7 @@ The simplest way to construct a polynomial is using the pattern `X`: ```haskell-> X^2 - 3*X + 2 :: UPoly Int+> X^2 - 3 * X + 2 :: UPoly Int 1 * X^2 + (-3) * X + 2 ``` @@ -39,11 +41,18 @@ 1 * X^2 + (-3) * X + 2 ``` +Alternatively one can enable `{-# LANGUAGE OverloadedLists #-}` and simply write++```haskell+> [2, -3, 1] :: UPoly Int+1 * X^2 + (-3) * X + 2+```+ There is a shortcut to construct a monomial: ```haskell-> monomial 2 3 :: UPoly Int-3 * X^2 + 0 * X + 0+> monomial 2 3.5 :: UPoly Double+3.5 * X^2 + 0.0 * X + 0.0 ``` ## Operations@@ -58,8 +67,8 @@ One can also find convenient to `scale` by monomial (cf. `monomial` above): ```haskell-> scale 2 3 (X^2 + 1) :: UPoly Int-3 * X^4 + 0 * X^3 + 3 * X^2 + 0 * X + 0+> scale 2 3.5 (X^2 + 1) :: UPoly Double+3.5 * X^4 + 0.0 * X^3 + 3.5 * X^2 + 0.0 * X + 0.0 ``` While `Poly` cannot be made an instance of `Integral` (because there is no meaningful `toInteger`),@@ -81,9 +90,6 @@ > 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 Double 3.0 * X^2 + 0.0 * X + 3.0 @@ -100,12 +106,12 @@ [2,-3,1] ``` -Further, `leading` is a shortcut to to obtain the leading term of a non-zero polynomial,+Further, `leading` is a shortcut to obtain the leading term of a non-zero polynomial, expressed as a power and a coefficient: ```haskell-> leading (X^2 - 3 * X + 2 :: UPoly Int)-Just (2,1)+> leading (X^2 - 3 * X + 2 :: UPoly Double)+Just (2,1.0) ``` ## Flavours@@ -122,3 +128,19 @@ because of a more readable `Show` instance, skipping zero coefficients. * `Data.Poly.Sparse.Semiring` provides sparse polynomials with `Semiring`-based interface.++All flavours are available backed by boxed or unboxed vectors.++## Performance++As a rough guide, `poly` is at least 20x-40x faster than [`polynomial`](http://hackage.haskell.org/package/polynomial) library.+Multiplication is implemented via Karatsuba algorithm.+Here is a couple of benchmarks for `UPoly Int`.++| Benchmark | polynomial, μs | poly, μs | speedup+| :---------------------------- | --------------: | -------: | ------:+| addition, 100 coeffs. | 45 | 2 | 22x+| addition, 1000 coeffs. | 441 | 17 | 25x+| addition, 10000 coeffs. | 6545 | 167 | 39x+| multiplication, 100 coeffs. | 1733 | 33 | 52x+| multiplication, 1000 coeffs. | 442000 | 1456 | 303x
bench/DenseBench.hs view
@@ -1,3 +1,4 @@+{-# LANGUAGE CPP #-} {-# LANGUAGE RankNTypes #-} module DenseBench@@ -6,21 +7,25 @@ import Prelude hiding (quotRem, gcd) import Gauge.Main-import Data.Euclidean import Data.Poly-import qualified Data.Vector as V import qualified Data.Vector.Unboxed as U+#if MIN_VERSION_semirings(0,4,2)+import Data.Euclidean+import qualified Data.Vector as V+#endif benchSuite :: Benchmark benchSuite = bgroup "dense" $ concat [ map benchAdd [100, 1000, 10000]- , map benchMul [10, 100]+ , map benchMul [100, 1000, 10000] , map benchEval [100, 1000, 10000] , map benchDeriv [100, 1000, 10000] , map benchIntegral [100, 1000, 10000]+#if MIN_VERSION_semirings(0,4,2) , map benchQuotRem [10, 100] , map benchGcdFrac [10, 100] , map benchGcd [10, 100]+#endif ] benchAdd :: Int -> Benchmark@@ -38,6 +43,8 @@ benchIntegral :: Int -> Benchmark benchIntegral k = bench ("integral/" ++ show k) $ nf doIntegral k +#if MIN_VERSION_semirings(0,4,2)+ benchQuotRem :: Int -> Benchmark benchQuotRem k = bench ("quotRem/" ++ show k) $ nf doQuotRem k @@ -47,6 +54,8 @@ benchGcdFrac :: Int -> Benchmark benchGcdFrac k = bench ("gcdFrac/" ++ show k) $ nf doGcdFrac k +#endif+ doBinOp :: (forall a. Num a => a -> a -> a) -> Int -> Int doBinOp op n = U.sum zs where@@ -72,6 +81,8 @@ xs = toPoly $ U.generate n ((* 2) . fromIntegral) zs = unPoly $ integral xs +#if MIN_VERSION_semirings(0,4,2)+ doQuotRem :: Int -> Double doQuotRem n = U.sum (unPoly qs) + U.sum (unPoly rs) where@@ -92,3 +103,5 @@ xs = PolyOverFractional $ toPoly $ V.generate n ((+ 1) . (* 2) . fromIntegral) ys = PolyOverFractional $ toPoly $ V.generate n ((+ 2) . (* 3) . fromIntegral) gs = unPoly $ unPolyOverFractional $ xs `gcd` ys++#endif
changelog.md view
@@ -1,3 +1,8 @@+# 0.3.1.0++* Implement Karatsuba multiplication.+* Add `IsList` instance.+ # 0.3.0.0 * Implement sparse polynomials.
poly.cabal view
@@ -1,5 +1,5 @@ name: poly-version: 0.3.0.0+version: 0.3.1.0 synopsis: Polynomials description: Polynomials backed by `Vector`.@@ -38,10 +38,10 @@ Data.Poly.Internal.Sparse.GcdDomain build-depends: base >= 4.9 && < 5,- primitive,- semirings >= 0.4,- vector,- vector-algorithms+ primitive >= 0.6,+ semirings >= 0.2,+ vector >= 0.12.0.2,+ vector-algorithms >= 0.7 default-language: Haskell2010 ghc-options: -Wall @@ -50,27 +50,28 @@ main-is: Main.hs other-modules: Dense+ Quaternion Sparse build-depends: base >=4.9 && <5, poly,- QuickCheck >=2.10,- quickcheck-classes >=0.6.1,- semirings,- tasty,- tasty-quickcheck,- vector+ QuickCheck >=2.12,+ quickcheck-classes >=0.5,+ semirings >= 0.2,+ tasty >= 0.11,+ tasty-quickcheck >= 0.8,+ vector >= 0.12.0.2 default-language: Haskell2010 hs-source-dirs: test ghc-options: -Wall benchmark poly-gauge build-depends:- base,- gauge,+ base >=4.9 && <5,+ gauge >= 0.1, poly,- semirings,- vector+ semirings >= 0.2,+ vector >= 0.12.0.2 type: exitcode-stdio-1.0 main-is: Bench.hs other-modules:
src/Data/Poly.hs view
@@ -7,6 +7,7 @@ -- Dense polynomials and a 'Num'-based interface. -- +{-# LANGUAGE CPP #-} {-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE GeneralizedNewtypeDeriving #-} {-# LANGUAGE PatternSynonyms #-}@@ -25,11 +26,15 @@ , eval , deriv , integral+#if MIN_VERSION_semirings(0,4,2) -- * Fractional coefficients , PolyOverFractional(..)+#endif ) where import Data.Poly.Internal.Dense+#if MIN_VERSION_semirings(0,4,2) import Data.Poly.Internal.Dense.Fractional () import Data.Poly.Internal.Dense.GcdDomain () import Data.Poly.Internal.PolyOverFractional+#endif
src/Data/Poly/Internal/Dense.hs view
@@ -7,6 +7,7 @@ -- Dense polynomials of one variable. -- +{-# LANGUAGE CPP #-} {-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE GeneralizedNewtypeDeriving #-} {-# LANGUAGE PatternSynonyms #-}@@ -37,7 +38,7 @@ , deriv' ) where -import Prelude hiding (quotRem, quot, rem, gcd, lcm, (^))+import Prelude hiding (quotRem, rem, gcd, lcm, (^)) import Control.Monad import Control.Monad.Primitive import Control.Monad.ST@@ -48,6 +49,11 @@ import qualified Data.Vector.Generic as G import qualified Data.Vector.Generic.Mutable as MG import qualified Data.Vector.Unboxed as U+import GHC.Exts+#if !MIN_VERSION_semirings(0,4,0)+import Data.Semigroup+import Numeric.Natural+#endif -- | Polynomials of one variable with coefficients from @a@, -- backed by a 'G.Vector' @v@ (boxed, unboxed, storable, etc.).@@ -72,6 +78,12 @@ } deriving (Eq, Ord) +instance (Eq a, Semiring a, G.Vector v a) => IsList (Poly v a) where+ type Item (Poly v a) = a+ fromList = toPoly' . G.fromList+ fromListN = (toPoly' .) . G.fromListN+ toList = G.toList . unPoly+ instance (Show a, G.Vector v a) => Show (Poly v a) where showsPrec d (Poly xs) | G.null xs@@ -129,7 +141,7 @@ fromInteger n = case fromInteger n of 0 -> Poly $ G.empty m -> Poly $ G.singleton m- Poly xs * Poly ys = toPoly $ convolution 0 (+) (*) xs ys+ Poly xs * Poly ys = toPoly $ karatsuba xs ys {-# INLINE (+) #-} {-# INLINE (-) #-} {-# INLINE negate #-}@@ -148,6 +160,14 @@ {-# INLINE plus #-} {-# INLINE times #-} +#if MIN_VERSION_semirings(0,4,0)+ fromNatural n = if n' == zero then zero else Poly $ G.singleton n'+ where+ n' :: a+ n' = fromNatural n+ {-# INLINE fromNatural #-}+#endif+ 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 @@ -224,6 +244,56 @@ G.unsafeFreeze zs {-# INLINE minusPoly #-} +karatsubaThreshold :: Int+karatsubaThreshold = 32++karatsuba+ :: (Eq a, Num a, G.Vector v a)+ => v a+ -> v a+ -> v a+karatsuba xs ys+ | lenXs <= karatsubaThreshold || lenYs <= karatsubaThreshold+ = convolution 0 (+) (*) xs ys+ | otherwise = runST $ do+ zs <- MG.basicUnsafeNew lenZs+ forM_ [0 .. lenZs - 1] $ \k -> do+ let z0 = if k < G.basicLength zs0+ then G.unsafeIndex zs0 k+ else 0+ z11 = if k - m >= 0 && k - m < G.basicLength zs11+ then G.unsafeIndex zs11 (k - m)+ else 0+ z10 = if k - m >= 0 && k - m < G.basicLength zs0+ then G.unsafeIndex zs0 (k - m)+ else 0+ z12 = if k - m >= 0 && k - m < G.basicLength zs2+ then G.unsafeIndex zs2 (k - m)+ else 0+ z2 = if k - 2 * m >= 0 && k - 2 * m < G.basicLength zs2+ then G.unsafeIndex zs2 (k - 2 * m)+ else 0+ MG.unsafeWrite zs k (z0 + (z11 - z10 - z12) + z2)+ G.unsafeFreeze zs+ where+ lenXs = G.basicLength xs+ lenYs = G.basicLength ys+ lenZs = lenXs + lenYs - 1++ m = ((lenXs `min` lenYs) + 1) `quot` 2++ xs0 = G.slice 0 m xs+ xs1 = G.slice m (lenXs - m) xs+ ys0 = G.slice 0 m ys+ ys1 = G.slice m (lenYs - m) ys++ xs01 = plusPoly (+) xs0 xs1+ ys01 = plusPoly (+) ys0 ys1+ zs0 = karatsuba xs0 ys0+ zs2 = karatsuba xs1 ys1+ zs11 = karatsuba xs01 ys01+{-# INLINE karatsuba #-}+ convolution :: G.Vector v a => a@@ -251,11 +321,13 @@ monomial :: (Eq a, Num a, G.Vector v a) => Word -> a -> Poly v a monomial _ 0 = Poly G.empty monomial p c = Poly $ G.generate (fromIntegral p + 1) (\i -> if i == fromIntegral p then c else 0)+{-# INLINE monomial #-} monomial' :: (Eq a, Semiring a, G.Vector v a) => Word -> a -> Poly v a monomial' p c | c == zero = Poly G.empty | otherwise = Poly $ G.generate (fromIntegral p + 1) (\i -> if i == fromIntegral p then c else zero)+{-# INLINE monomial' #-} scaleInternal :: (Eq a, G.Vector v a)@@ -263,9 +335,9 @@ -> (a -> a -> a) -> Word -> a- -> Poly v a -> v a-scaleInternal zer mul yp yc (Poly xs) = runST $ do+ -> v a+scaleInternal zer mul yp yc xs = runST $ do let lenXs = G.basicLength xs zs <- MG.basicUnsafeNew (fromIntegral yp + lenXs) forM_ [0 .. fromIntegral yp - 1] $ \k ->@@ -273,16 +345,17 @@ forM_ [0 .. lenXs - 1] $ \k -> MG.unsafeWrite zs (fromIntegral yp + k) (mul yc $ G.unsafeIndex xs k) G.unsafeFreeze zs+{-# INLINE scaleInternal #-} -- | Multiply a polynomial by a monomial, expressed as a power and a coefficient. -- -- >>> scale 2 3 (X^2 + 1) :: UPoly Int -- 3 * X^4 + 0 * X^3 + 3 * X^2 + 0 * X + 0 scale :: (Eq a, Num a, G.Vector v a) => Word -> a -> Poly v a -> Poly v a-scale yp yc xs = toPoly $ scaleInternal 0 (*) yp yc xs+scale yp yc (Poly xs) = toPoly $ scaleInternal 0 (*) yp yc xs scale' :: (Eq a, Semiring a, G.Vector v a) => Word -> a -> Poly v a -> Poly v a-scale' yp yc xs = toPoly' $ scaleInternal zero times yp yc xs+scale' yp yc (Poly xs) = toPoly' $ scaleInternal zero times yp yc xs data StrictPair a b = !a :*: !b @@ -322,6 +395,17 @@ | G.null xs = Poly G.empty | otherwise = toPoly' $ G.imap (\i x -> fromNatural (fromIntegral (i + 1)) `times` x) $ G.tail xs {-# INLINE deriv' #-}++#if !MIN_VERSION_semirings(0,4,0)+fromNatural :: Semiring a => Natural -> a+fromNatural 0 = zero+fromNatural n = getAdd' (stimes n (Add' one))++newtype Add' a = Add' { getAdd' :: a }++instance Semiring a => Semigroup (Add' a) where+ Add' a <> Add' b = Add' (a `plus` b)+#endif -- | Compute an indefinite integral of a polynomial, -- setting constant term to zero.
src/Data/Poly/Internal/Dense/Fractional.hs view
@@ -7,6 +7,7 @@ -- GcdDomain for Fractional underlying -- +{-# LANGUAGE CPP #-} {-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE GeneralizedNewtypeDeriving #-} {-# LANGUAGE PatternSynonyms #-}@@ -16,6 +17,8 @@ {-# OPTIONS_GHC -fno-warn-orphans #-} +#if MIN_VERSION_semirings(0,4,2)+ module Data.Poly.Internal.Dense.Fractional ( fractionalGcd ) where@@ -127,3 +130,9 @@ remainderM xs ys' gcdM ys' xs {-# INLINE gcdM #-}++#else++module Data.Poly.Internal.Dense.Fractional () where++#endif
src/Data/Poly/Internal/Dense/GcdDomain.hs view
@@ -7,6 +7,7 @@ -- GcdDomain for GcdDomain underlying -- +{-# LANGUAGE CPP #-} {-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE GeneralizedNewtypeDeriving #-} {-# LANGUAGE PatternSynonyms #-}@@ -19,6 +20,8 @@ module Data.Poly.Internal.Dense.GcdDomain () where +#if MIN_VERSION_semirings(0,4,2)+ import Prelude hiding (gcd, lcm, (^)) import Control.Exception import Control.Monad@@ -171,3 +174,5 @@ go (lenQs - 1) {-# INLINE quotient #-}++#endif
src/Data/Poly/Internal/PolyOverFractional.hs view
@@ -7,10 +7,13 @@ -- Wrapper with a more efficient 'Euclidean' instance. -- +{-# LANGUAGE CPP #-} {-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE GeneralizedNewtypeDeriving #-} {-# LANGUAGE UndecidableInstances #-} +#if MIN_VERSION_semirings(0,4,2)+ module Data.Poly.Internal.PolyOverFractional ( PolyOverFractional(..) ) where@@ -44,3 +47,9 @@ rem (PolyOverFractional x) (PolyOverFractional y) = PolyOverFractional (rem x y) {-# INLINE rem #-}++#else++module Data.Poly.Internal.PolyOverFractional () where++#endif
src/Data/Poly/Internal/Sparse.hs view
@@ -7,10 +7,12 @@ -- Sparse polynomials of one variable. -- +{-# LANGUAGE CPP #-} {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE PatternSynonyms #-} {-# LANGUAGE ScopedTypeVariables #-} {-# LANGUAGE StandaloneDeriving #-}+{-# LANGUAGE TypeFamilies #-} {-# LANGUAGE UndecidableInstances #-} {-# LANGUAGE ViewPatterns #-} @@ -48,6 +50,11 @@ import qualified Data.Vector.Generic.Mutable as MG import qualified Data.Vector.Unboxed as U import qualified Data.Vector.Algorithms.Tim as Tim+import GHC.Exts+#if !MIN_VERSION_semirings(0,4,0)+import Data.Semigroup+import Numeric.Natural+#endif -- | Polynomials of one variable with coefficients from @a@, -- backed by a 'G.Vector' @v@ (boxed, unboxed, storable, etc.).@@ -74,6 +81,12 @@ deriving instance Eq (v (Word, a)) => Eq (Poly v a) deriving instance Ord (v (Word, a)) => Ord (Poly v a) +instance (Eq a, Semiring a, G.Vector v (Word, a)) => IsList (Poly v a) where+ type Item (Poly v a) = (Word, a)+ fromList = toPoly' . G.fromList+ fromListN = (toPoly' .) . G.fromListN+ toList = G.toList . unPoly+ instance (Show a, G.Vector v (Word, a)) => Show (Poly v a) where showsPrec d (Poly xs) | G.null xs@@ -184,15 +197,18 @@ | otherwise = Poly $ G.singleton (0, one) plus (Poly xs) (Poly ys) = Poly $ plusPoly (/= zero) plus xs ys times (Poly xs) (Poly ys) = Poly $ convolution (/= zero) plus times xs ys- fromNatural n = if n' == zero then zero else Poly $ G.singleton (0, n')- where- n' :: a- n' = fromNatural n {-# INLINE zero #-} {-# INLINE one #-} {-# INLINE plus #-} {-# INLINE times #-}++#if MIN_VERSION_semirings(0,4,0)+ fromNatural n = if n' == zero then zero else Poly $ G.singleton (0, n')+ where+ n' :: a+ n' = fromNatural n {-# INLINE fromNatural #-}+#endif instance (Eq a, Semiring.Ring a, G.Vector v (Word, a)) => Semiring.Ring (Poly v a) where negate (Poly xs) = Poly $ G.map (fmap Semiring.negate) xs@@ -472,6 +488,17 @@ (\p c -> fromNatural (fromIntegral p) `times` c) xs {-# INLINE deriv' #-}++#if !MIN_VERSION_semirings(0,4,0)+fromNatural :: Semiring a => Natural -> a+fromNatural 0 = zero+fromNatural n = getAdd' (stimes n (Add' one))++newtype Add' a = Add' { getAdd' :: a }++instance Semiring a => Semigroup (Add' a) where+ Add' a <> Add' b = Add' (a `plus` b)+#endif derivPoly :: G.Vector v (Word, a)
src/Data/Poly/Internal/Sparse/Fractional.hs view
@@ -7,6 +7,7 @@ -- GcdDomain for Fractional underlying -- +{-# LANGUAGE CPP #-} {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE GeneralizedNewtypeDeriving #-}@@ -18,6 +19,8 @@ {-# OPTIONS_GHC -fno-warn-orphans #-} +#if MIN_VERSION_semirings(0,4,2)+ module Data.Poly.Internal.Sparse.Fractional ( fractionalGcd ) where@@ -67,3 +70,9 @@ | G.null (unPoly ys) = xs | otherwise = fractionalGcd ys $ snd $ quotientRemainder xs ys {-# INLINE fractionalGcd #-}++#else++module Data.Poly.Internal.Sparse.Fractional () where++#endif
src/Data/Poly/Internal/Sparse/GcdDomain.hs view
@@ -7,6 +7,7 @@ -- GcdDomain for GcdDomain underlying -- +{-# LANGUAGE CPP #-} {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE GeneralizedNewtypeDeriving #-}@@ -21,6 +22,8 @@ module Data.Poly.Internal.Sparse.GcdDomain () where +#if MIN_VERSION_semirings(0,4,2)+ import Prelude hiding (gcd, lcm, (^)) import Control.Exception import Data.Euclidean@@ -76,3 +79,5 @@ gx = fromMaybe err $ divide g xc gy = fromMaybe err $ divide g yc err = error "gcd: violated internal invariant"++#endif
src/Data/Poly/Semiring.hs view
@@ -7,6 +7,7 @@ -- Dense polynomials and a 'Semiring'-based interface. -- +{-# LANGUAGE CPP #-} {-# LANGUAGE PatternSynonyms #-} module Data.Poly.Semiring@@ -22,8 +23,10 @@ , pattern X , eval , deriv+#if MIN_VERSION_semirings(0,4,2) -- * Fractional coefficients , PolyOverFractional(..)+#endif ) where import Data.Semiring (Semiring)@@ -31,9 +34,11 @@ import Data.Poly.Internal.Dense (Poly(..), VPoly, UPoly, leading) import qualified Data.Poly.Internal.Dense as Dense+#if MIN_VERSION_semirings(0,4,2) import Data.Poly.Internal.Dense.Fractional () import Data.Poly.Internal.Dense.GcdDomain () import Data.Poly.Internal.PolyOverFractional+#endif -- | Make 'Poly' from a vector of coefficients -- (first element corresponds to a constant term).
test/Dense.hs view
@@ -1,3 +1,5 @@+{-# LANGUAGE CPP #-}+{-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE GeneralizedNewtypeDeriving #-} {-# LANGUAGE ScopedTypeVariables #-}@@ -9,7 +11,9 @@ ) where import Prelude hiding (quotRem)+#if MIN_VERSION_semirings(0,4,2) import Data.Euclidean+#endif import Data.Int import Data.Poly import qualified Data.Poly.Semiring as S@@ -22,16 +26,28 @@ import Test.Tasty.QuickCheck hiding (scale) import Test.QuickCheck.Classes +import Quaternion+ 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 +#if MIN_VERSION_semirings(0,4,2) instance (Eq a, Semiring a, Arbitrary a, G.Vector v a) => Arbitrary (PolyOverFractional (Poly v a)) where arbitrary = PolyOverFractional . S.toPoly . G.fromList . (\xs -> take (length xs `mod` 10) xs) <$> arbitrary shrink = fmap (PolyOverFractional . S.toPoly . G.fromList) . shrink . G.toList . unPoly . unPolyOverFractional+#endif newtype ShortPoly a = ShortPoly { unShortPoly :: a }- deriving (Eq, Show, Semiring, GcdDomain, Euclidean)+ deriving+ ( Eq+ , Show+ , Semiring+#if MIN_VERSION_semirings(0,4,2)+ , GcdDomain+ , Euclidean+#endif+ ) instance (Eq a, Semiring a, Arbitrary a, G.Vector v a) => Arbitrary (ShortPoly (Poly v a)) where arbitrary = ShortPoly . S.toPoly . G.fromList . (\xs -> take (length xs `mod` 10) xs) <$> arbitrary@@ -44,7 +60,9 @@ , semiringTests , evalTests , derivTests+#if MIN_VERSION_semirings(0,4,2) -- , euclideanTests+#endif ] semiringTests :: TestTree@@ -53,13 +71,18 @@ $ map (uncurry testProperty) $ concatMap lawsProperties [ 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))+ , semiringLaws (Proxy :: Proxy (Poly U.Vector (Quaternion Int)))+#if MIN_VERSION_quickcheck_classes(0,6,1)+ , ringLaws (Proxy :: Proxy (Poly U.Vector ()))+ , ringLaws (Proxy :: Proxy (Poly U.Vector Int8))+ , ringLaws (Proxy :: Proxy (Poly V.Vector Integer))+ , ringLaws (Proxy :: Proxy (Poly U.Vector (Quaternion Int)))+#endif ] +#if MIN_VERSION_semirings(0,4,2) -- euclideanTests :: TestTree -- euclideanTests -- = testGroup "Euclidean"@@ -69,6 +92,7 @@ -- , gcdDomainLaws (Proxy :: Proxy (PolyOverFractional (Poly V.Vector Rational))) -- , euclideanLaws (Proxy :: Proxy (ShortPoly (Poly V.Vector Rational))) -- ]+#endif arithmeticTests :: TestTree arithmeticTests = testGroup "Arithmetic"@@ -100,15 +124,25 @@ $ iterate (0 :) ys otherTests :: TestTree-otherTests = testGroup "Other"+otherTests = testGroup "other" $ concat+ [ otherTestGroup (Proxy :: Proxy Int8)+ , otherTestGroup (Proxy :: Proxy (Quaternion Int))+ ]++otherTestGroup+ :: forall a.+ (Eq a, Show a, Semiring a, Num a, Arbitrary a, U.Unbox a, G.Vector U.Vector a)+ => Proxy a+ -> [TestTree]+otherTestGroup _ = [ testProperty "leading p 0 == Nothing" $- \p -> leading (monomial p 0 :: UPoly Int) === Nothing+ \p -> leading (monomial p 0 :: UPoly a) === Nothing , testProperty "leading . monomial = id" $- \p c -> c /= 0 ==> leading (monomial p c :: UPoly Int) === Just (p, c)+ \p c -> c /= 0 ==> leading (monomial p c :: UPoly a) === Just (p, c) , testProperty "monomial matches reference" $- \p (c :: Int) -> monomial p c === toPoly (V.fromList (monomialRef p c))+ \p (c :: a) -> monomial p c === toPoly (V.fromList (monomialRef p c)) , testProperty "scale matches multiplication by monomial" $- \p c (xs :: UPoly Int) -> scale p c xs === monomial p c * xs+ \p c (xs :: UPoly a) -> scale p c xs === monomial p c * xs ] monomialRef :: Num a => Word -> a -> [a]
+ test/Quaternion.hs view
@@ -0,0 +1,113 @@+-- |+-- Module: Quaternion+-- Copyright: (c) 2019 Andrew Lelechenko+-- Licence: BSD3+-- Maintainer: Andrew Lelechenko <andrew.lelechenko@gmail.com>+--+-- This is a toy implementtion of quaternions,+-- serving solely to test polynomials+-- over non-commutative rings.+--++{-# LANGUAGE DeriveGeneric #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE TypeFamilies #-}++module Quaternion+ ( Quaternion(..)+ ) where++import Prelude hiding (negate)+import Control.Monad+import Data.Semiring (Semiring(..), Ring(..), minus)+import GHC.Generics+import Test.Tasty.QuickCheck hiding (scale)++import Data.Vector.Unboxed (Vector)+import qualified Data.Vector.Generic as G+import Data.Vector.Unboxed.Mutable (MVector)+import qualified Data.Vector.Generic.Mutable as M+import Data.Vector.Unboxed (Unbox)++data Quaternion a = Quaternion a a a a+ deriving (Eq, Ord, Show, Generic)++instance Ring a => Semiring (Quaternion a) where+ zero = Quaternion zero zero zero zero+ one = Quaternion one zero zero zero+ plus (Quaternion a1 b1 c1 d1) (Quaternion a2 b2 c2 d2) =+ Quaternion (a1 `plus` a2) (b1 `plus` b2) (c1 `plus` c2) (d1 `plus` d2)+ times (Quaternion a1 b1 c1 d1) (Quaternion a2 b2 c2 d2) =+ Quaternion+ (a1 `times` a2 `minus` b1 `times` b2 `minus` c1 `times` c2 `minus` d1 `times` d2)+ (a1 `times` b2 `plus` b1 `times` a2 `plus` c1 `times` d2 `minus` d1 `times` c2)+ (a1 `times` c2 `minus` b1 `times` d2 `plus` c1 `times` a2 `plus` d1 `times` b2)+ (a1 `times` d2 `plus` b1 `times` c2 `minus` c1 `times` b2 `plus` d1 `times` a2)++instance Ring a => Ring (Quaternion a) where+ negate (Quaternion a b c d) =+ Quaternion (negate a) (negate b) (negate c) (negate d)++instance (Ring a, Num a) => Num (Quaternion a) where+ (+) = plus+ (-) = minus+ (*) = times+ abs = id+ signum = const one+ fromInteger n = Quaternion (fromInteger n) zero zero zero++instance Arbitrary a => Arbitrary (Quaternion a) where+ arbitrary = Quaternion <$> arbitrary <*> arbitrary <*> arbitrary <*> arbitrary+ shrink = genericShrink++newtype instance MVector s (Quaternion a) = MV_Quaternion (MVector s (a, a, a, a))+newtype instance Vector (Quaternion a) = V_Quaternion (Vector (a, a, a, a))++instance (Unbox a) => Unbox (Quaternion a)++instance (Unbox a) => M.MVector MVector (Quaternion a) where+ {-# INLINE basicLength #-}+ {-# INLINE basicUnsafeSlice #-}+ {-# INLINE basicOverlaps #-}+ {-# INLINE basicUnsafeNew #-}+ {-# INLINE basicInitialize #-}+ {-# INLINE basicUnsafeReplicate #-}+ {-# INLINE basicUnsafeRead #-}+ {-# INLINE basicUnsafeWrite #-}+ {-# INLINE basicClear #-}+ {-# INLINE basicSet #-}+ {-# INLINE basicUnsafeCopy #-}+ {-# INLINE basicUnsafeGrow #-}+ basicLength (MV_Quaternion v) = M.basicLength v+ basicUnsafeSlice i n (MV_Quaternion v) = MV_Quaternion $ M.basicUnsafeSlice i n v+ basicOverlaps (MV_Quaternion v1) (MV_Quaternion v2) = M.basicOverlaps v1 v2+ basicUnsafeNew n = MV_Quaternion `liftM` M.basicUnsafeNew n+ basicInitialize (MV_Quaternion v) = M.basicInitialize v+ basicUnsafeReplicate n (Quaternion a b c d) = MV_Quaternion `liftM` M.basicUnsafeReplicate n (a, b, c, d)+ basicUnsafeRead (MV_Quaternion v) i = (\(a, b, c, d) -> Quaternion a b c d) `liftM` M.basicUnsafeRead v i+ basicUnsafeWrite (MV_Quaternion v) i (Quaternion a b c d) = M.basicUnsafeWrite v i (a, b, c, d)+ basicClear (MV_Quaternion v) = M.basicClear v+ basicSet (MV_Quaternion v) (Quaternion a b c d) = M.basicSet v (a, b, c, d)+ basicUnsafeCopy (MV_Quaternion v1) (MV_Quaternion v2) = M.basicUnsafeCopy v1 v2+ basicUnsafeMove (MV_Quaternion v1) (MV_Quaternion v2) = M.basicUnsafeMove v1 v2+ basicUnsafeGrow (MV_Quaternion v) n = MV_Quaternion `liftM` M.basicUnsafeGrow v n++instance (Unbox a) => G.Vector Vector (Quaternion a) where+ {-# INLINE basicUnsafeFreeze #-}+ {-# INLINE basicUnsafeThaw #-}+ {-# INLINE basicLength #-}+ {-# INLINE basicUnsafeSlice #-}+ {-# INLINE basicUnsafeIndexM #-}+ {-# INLINE elemseq #-}+ basicUnsafeFreeze (MV_Quaternion v) = V_Quaternion `liftM` G.basicUnsafeFreeze v+ basicUnsafeThaw (V_Quaternion v) = MV_Quaternion `liftM` G.basicUnsafeThaw v+ basicLength (V_Quaternion v) = G.basicLength v+ basicUnsafeSlice i n (V_Quaternion v) = V_Quaternion $ G.basicUnsafeSlice i n v+ basicUnsafeIndexM (V_Quaternion v) i+ = (\(a, b, c, d) -> Quaternion a b c d) `liftM` G.basicUnsafeIndexM v i+ basicUnsafeCopy (MV_Quaternion mv) (V_Quaternion v)+ = G.basicUnsafeCopy mv v+ elemseq _ (Quaternion a b c d) z = G.elemseq (undefined :: Vector a) a+ $ G.elemseq (undefined :: Vector a) b+ $ G.elemseq (undefined :: Vector a) c+ $ G.elemseq (undefined :: Vector a) d z
test/Sparse.hs view
@@ -1,3 +1,4 @@+{-# LANGUAGE CPP #-} {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE GeneralizedNewtypeDeriving #-}@@ -11,7 +12,9 @@ ) where import Prelude hiding (quotRem)+#if MIN_VERSION_semirings(0,4,2) import Data.Euclidean+#endif import Data.Function import Data.Int import Data.List@@ -26,12 +29,22 @@ import Test.Tasty.QuickCheck hiding (scale) import Test.QuickCheck.Classes +import Quaternion+ instance (Eq a, Semiring a, Arbitrary a, G.Vector v (Word, a)) => Arbitrary (Poly v a) where arbitrary = S.toPoly . G.fromList <$> arbitrary shrink = fmap (S.toPoly . G.fromList) . shrink . G.toList . unPoly newtype ShortPoly a = ShortPoly { unShortPoly :: a }- deriving (Eq, Show, Semiring, GcdDomain, Euclidean)+ deriving+ ( Eq+ , Show+ , Semiring+#if MIN_VERSION_semirings(0,4,2)+ , GcdDomain+ , Euclidean+#endif+ ) instance (Eq a, Semiring a, Arbitrary a, G.Vector v (Word, a)) => Arbitrary (ShortPoly (Poly v a)) where arbitrary = ShortPoly . S.toPoly . G.fromList . (\xs -> take (length xs `mod` 5) xs) <$> arbitrary@@ -52,11 +65,15 @@ $ map (uncurry testProperty) $ concatMap lawsProperties [ 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))+ , semiringLaws (Proxy :: Proxy (Poly U.Vector (Quaternion Int)))+#if MIN_VERSION_quickcheck_classes(0,6,1)+ , ringLaws (Proxy :: Proxy (Poly U.Vector ()))+ , ringLaws (Proxy :: Proxy (Poly U.Vector Int8))+ , ringLaws (Proxy :: Proxy (Poly V.Vector Integer))+ , ringLaws (Proxy :: Proxy (Poly U.Vector (Quaternion Int)))+#endif ] arithmeticTests :: TestTree@@ -98,15 +115,25 @@ $ [ (xp + yp, xc * yc) | (xp, xc) <- xs, (yp, yc) <- ys ] otherTests :: TestTree-otherTests = testGroup "Other"+otherTests = testGroup "other" $ concat+ [ otherTestGroup (Proxy :: Proxy Int8)+ , otherTestGroup (Proxy :: Proxy (Quaternion Int))+ ]++otherTestGroup+ :: forall a.+ (Eq a, Show a, Semiring a, Num a, Arbitrary a, U.Unbox a, G.Vector U.Vector a)+ => Proxy a+ -> [TestTree]+otherTestGroup _ = [ testProperty "leading p 0 == Nothing" $- \p -> leading (monomial p 0 :: UPoly Int) === Nothing+ \p -> leading (monomial p 0 :: UPoly a) === Nothing , testProperty "leading . monomial = id" $- \p c -> c /= 0 ==> leading (monomial p c :: UPoly Int) === Just (p, c)+ \p c -> c /= 0 ==> leading (monomial p c :: UPoly a) === Just (p, c) , testProperty "monomial matches reference" $- \p (c :: Int) -> monomial p c === toPoly (V.fromList (monomialRef p c))+ \p (c :: a) -> monomial p c === toPoly (V.fromList (monomialRef p c)) , testProperty "scale matches multiplication by monomial" $- \p c (xs :: UPoly Int) -> scale p c xs === monomial p c * xs+ \p c (xs :: UPoly a) -> scale p c xs === monomial p c * xs ] monomialRef :: Num a => Word -> a -> [(Word, a)]