poly 0.5.0.0 → 0.5.1.0

raw patch · 34 files changed

+863/−319 lines, 34 filesdep +quickcheck-classes-basedep +tasty-benchdep −doctestdep −gaugedep ~basedep ~deepseqsetup-changed

Dependencies added: quickcheck-classes-base, tasty-bench

Dependencies removed: doctest, gauge

Dependency ranges changed: base, deepseq

Files

README.md +44/−24
Setup.hs +0/−2
bench/Bench.hs +5/−0
changelog.md +10/−1
poly.cabal +66/−43
src/Data/Poly.hs +7/−0
src/Data/Poly/Internal/Convert.hs +4/−0
src/Data/Poly/Internal/Dense.hs +127/−68
src/Data/Poly/Internal/Dense/DFT.hs +6/−0
src/Data/Poly/Internal/Dense/Field.hs +9/−6
src/Data/Poly/Internal/Dense/GcdDomain.hs +12/−9
src/Data/Poly/Internal/Dense/Laurent.hs +46/−8
src/Data/Poly/Internal/Multi.hs +76/−17
src/Data/Poly/Internal/Multi/Core.hs +23/−25
src/Data/Poly/Internal/Multi/Field.hs +3/−1
src/Data/Poly/Internal/Multi/GcdDomain.hs +2/−13
src/Data/Poly/Internal/Multi/Laurent.hs +66/−20
src/Data/Poly/Laurent.hs +2/−0
src/Data/Poly/Multi.hs +1/−0
src/Data/Poly/Multi/Semiring.hs +29/−8
src/Data/Poly/Orthogonal.hs +26/−5
src/Data/Poly/Semiring.hs +47/−12
src/Data/Poly/Sparse.hs +34/−7
src/Data/Poly/Sparse/Laurent.hs +32/−5
src/Data/Poly/Sparse/Semiring.hs +42/−9
test/Dense.hs +24/−5
test/DenseLaurent.hs +13/−0
test/Main.hs +8/−2
test/Multi.hs +13/−0
test/MultiLaurent.hs +13/−0
test/Sparse.hs +17/−0
test/SparseLaurent.hs +13/−0
test/TestUtils.hs +43/−25
test/doctests.hs +0/−4

README.md view

@@ -1,18 +1,30 @@-# poly [![Build Status](https://travis-ci.org/Bodigrim/poly.svg)](https://travis-ci.org/Bodigrim/poly) [![Hackage](http://img.shields.io/hackage/v/poly.svg)](https://hackage.haskell.org/package/poly) [![Hackage CI](https://matrix.hackage.haskell.org/api/v2/packages/poly/badge)](https://matrix.hackage.haskell.org/package/poly) [![Stackage LTS](http://stackage.org/package/poly/badge/lts)](http://stackage.org/lts/package/poly) [![Stackage Nightly](http://stackage.org/package/poly/badge/nightly)](http://stackage.org/nightly/package/poly) [![Coverage Status](https://coveralls.io/repos/github/Bodigrim/poly/badge.svg)](https://coveralls.io/github/Bodigrim/poly)+# poly [![Hackage](https://img.shields.io/hackage/v/poly.svg)](https://hackage.haskell.org/package/poly) [![Stackage LTS](https://www.stackage.org/package/poly/badge/lts)](https://www.stackage.org/lts/package/poly) [![Stackage Nightly](https://www.stackage.org/package/poly/badge/nightly)](https://www.stackage.org/nightly/package/poly) [![Coverage Status](https://coveralls.io/repos/github/Bodigrim/poly/badge.svg)](https://coveralls.io/github/Bodigrim/poly) -Haskell library for univariate and multivariate polynomials, backed by `Vector`.+Haskell library for univariate and multivariate polynomials, backed by `Vector`s.  ```haskell+> -- Univariate polynomials > (X + 1) + (X - 1) :: VPoly Integer-2 * X + 0-+2 * X > (X + 1) * (X - 1) :: UPoly Int-1 * X^2 + 0 * X + (-1)+1 * X^2 + (-1)++> -- Multivariate polynomials+> (X + Y) * (X - Y) :: VMultiPoly 2 Integer+1 * X^2 + (-1) * Y^2+> (X + Y + Z) ^ 2 :: UMultiPoly 3 Int+1 * X^2 + 2 * X * Y + 2 * X * Z + 1 * Y^2 + 2 * Y * Z + 1 * Z^2++> -- Laurent polynomials+> (X^-2 + 1) * (X - X^-1) :: VLaurent Integer+1 * X + (-1) * X^-3+> (X^-1 + Y) * (X + Y^-1) :: UMultiLaurent 2 Int+1 * X * Y + 2 + 1 * X^-1 * Y^-1 ```  ## Vectors -`Poly v a` is polymorphic over a container `v`, implementing `Vector` interface, and coefficients of type `a`. Usually `v` is either a boxed vector from `Data.Vector` or an unboxed vector from `Data.Vector.Unboxed`. Use unboxed vectors whenever possible, e. g., when coefficients are `Int` or `Double`.+`Poly v a` is polymorphic over a container `v`, implementing the `Vector` interface, and coefficients of type `a`. Usually `v` is either a boxed vector from [`Data.Vector`](https://hackage.haskell.org/package/vector/docs/Data-Vector.html) or an unboxed vector from [`Data.Vector.Unboxed`](https://hackage.haskell.org/package/vector/docs/Data-Vector-Unboxed.html). Use unboxed vectors whenever possible, e. g., when the coefficients are `Int`s or `Double`s.  There are handy type synonyms: @@ -62,7 +74,7 @@ 1 * X^4 + 0 * X^3 + 0 * X^2 + 0 * X + (-1) ``` -One can also find convenient to `scale` by monomial (cf. `monomial` above):+One can also find it convenient to `scale` by a monomial (cf. `monomial` above):  ```haskell > scale 2 3.5 (X^2 + 1) :: UPoly Double@@ -70,8 +82,8 @@ ```  While `Poly` cannot be made an instance of `Integral` (because there is no meaningful `toInteger`),-it is an instance of `GcdDomain` and `Euclidean` from `semirings` package. These type classes-cover main functionality of `Integral`, providing division with remainder and `gcd` / `lcm`:+it is an instance of `GcdDomain` and `Euclidean` from the [`semirings`](https://hackage.haskell.org/package/semirings) package. These type classes+cover the main functionality of `Integral`, providing division with remainder and `gcd` / `lcm`:  ```haskell > Data.Euclidean.gcd (X^2 + 7 * X + 6) (X^2 - 5 * X - 6) :: UPoly Int@@ -81,8 +93,8 @@ (1.0 * X + 0.0,1.0 * X + 2.0) ``` -Miscellaneous utilities include `eval` for evaluation at a given value of indeterminate,-and reciprocals `deriv` / `integral`:+Miscellaneous utilities include `eval` for evaluation at a given point,+and `deriv` / `integral` for taking the derivative and an indefinite integral, respectively:  ```haskell > eval (X^2 + 1 :: UPoly Int) 3@@ -114,34 +126,34 @@  ## Flavours -* `Data.Poly` provides dense univariate polynomials with `Num`-based interface.+* `Data.Poly` provides dense univariate polynomials with a `Num`-based interface.   This is a default choice for most users. -* `Data.Poly.Semiring` provides dense univariate polynomials with `Semiring`-based interface.+* `Data.Poly.Semiring` provides dense univariate polynomials with a `Semiring`-based interface. -* `Data.Poly.Laurent` provides dense univariate Laurent polynomials with `Semiring`-based interface.+* `Data.Poly.Laurent` provides dense univariate Laurent polynomials with a `Semiring`-based interface. -* `Data.Poly.Sparse` provides sparse univariate polynomials with `Num`-based interface.-  Besides that, you may find it easier to use in REPL+* `Data.Poly.Sparse` provides sparse univariate polynomials with a `Num`-based interface.+  Besides that, you may find it easier to use in the REPL   because of a more readable `Show` instance, skipping zero coefficients. -* `Data.Poly.Sparse.Semiring` provides sparse univariate polynomials with `Semiring`-based interface.+* `Data.Poly.Sparse.Semiring` provides sparse univariate polynomials with a `Semiring`-based interface. -* `Data.Poly.Sparse.Laurent` provides sparse univariate Laurent polynomials with `Semiring`-based interface.+* `Data.Poly.Sparse.Laurent` provides sparse univariate Laurent polynomials with a `Semiring`-based interface. -* `Data.Poly.Multi` provides sparse multivariate polynomials with `Num`-based interface.+* `Data.Poly.Multi` provides sparse multivariate polynomials with a `Num`-based interface. -* `Data.Poly.Multi.Semiring` provides sparse multivariate polynomials with `Semiring`-based interface.+* `Data.Poly.Multi.Semiring` provides sparse multivariate polynomials with a `Semiring`-based interface. -* `Data.Poly.Multi.Laurent` provides sparse multivariate Laurent polynomials with `Semiring`-based interface.+* `Data.Poly.Multi.Laurent` provides sparse multivariate Laurent polynomials with a `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`.+As a rough guide, `poly` is at least 20x-40x faster than the [`polynomial`](http://hackage.haskell.org/package/polynomial) library.+Multiplication is implemented via the Karatsuba algorithm.+Here are a couple of benchmarks for `UPoly Int`:  | Benchmark                     | polynomial, μs  | poly, μs | speedup | :---------------------------- | --------------: | -------: | ------:@@ -150,3 +162,11 @@ | addition, 10000 coeffs.       |            6545 |     167  |  39x | multiplication, 100 coeffs.   |            1733 |      33  |  52x | multiplication, 1000 coeffs.  |          442000 |    1456  | 303x++Due to being polymorphic by multiple axis, the performance of `poly` crucially depends on specialisation of instances. Clients are strongly recommended to compile with `ghc-options: -fspecialise-aggressively` and suggested to enable `-O2`.++## Additional resources++* __Polynomials in Haskell__, MuniHac, 12.09.2020:+  [slides](https://github.com/Bodigrim/my-talks/raw/master/munihac2020/slides.pdf),+  [video](https://youtu.be/NAs3ExQZUjA).

− Setup.hs

@@ -1,2 +0,0 @@-import Distribution.Simple-main = defaultMain

bench/Bench.hs view

@@ -1,13 +1,18 @@+{-# LANGUAGE CPP        #-} {-# LANGUAGE RankNTypes #-}  module Main where  import Gauge.Main import qualified DenseBench as Dense+#ifdef SupportSparse import qualified SparseBench as Sparse+#endif  main :: IO () main = defaultMain   [ Dense.benchSuite+#ifdef SupportSparse   , Sparse.benchSuite+#endif   ]

changelog.md view

@@ -1,3 +1,8 @@+# 0.5.1.0++* Add function `timesRing`.+* Tweak inlining pragmas.+ # 0.5.0.0  * Change definition of `Data.Euclidean.degree`@@ -5,7 +10,7 @@ * Implement multivariate polynomials (usual and Laurent). * Reimplement sparse univariate polynomials as a special case of multivariate ones. * Speed up `gcd` calculations for all flavours of polynomials.-* Decomission `PolyOverField`: it does not improve performance any more.+* Decomission `PolyOverField` and `LaurentOverField`: they do not improve performance any more. * Add function `quotRemFractional`. * Add an experimental implementation of the discrete Fourier transform. * Add conversion functions between dense and sparse polynomials.@@ -43,9 +48,13 @@  # 0.2.0.0 +* Parametrize `Poly` by underlying vector type.+* Introduce `Data.Poly.Semiring` module. * Fix a bug in `Num.(-)`. * Add functions `constant`, `eval`, `deriv`, `integral`. * Add a handy pattern synonym `X`.+* Add type synonyms `VPoly` and `UPoly`.+* Remove function `toPoly'`.  # 0.1.0.0

poly.cabal view

@@ -1,8 +1,8 @@ name: poly-version: 0.5.0.0+version: 0.5.1.0 synopsis: Polynomials description:-  Polynomials backed by `Vector`.+  Polynomials backed by `Vector`s. homepage: https://github.com/Bodigrim/poly#readme license: BSD3 license-file: LICENSE@@ -11,8 +11,8 @@ copyright: 2019-2020 Andrew Lelechenko category: Math, Numerical build-type: Simple-cabal-version: >=1.10-tested-with: GHC ==8.2.2 GHC ==8.4.4 GHC ==8.6.5 GHC ==8.8.4 GHC ==8.10.2+cabal-version: 2.0+tested-with: GHC ==8.6.5 GHC ==8.8.4 GHC ==8.10.7 GHC ==9.0.2 GHC ==9.2.5 GHC ==9.4.4 extra-source-files:   changelog.md   README.md@@ -21,6 +21,12 @@   type: git   location: https://github.com/Bodigrim/poly +flag sparse+  description:+    Enable sparse and multivariate polynomials, incurring a larger dependency footprint.+  default: True+  manual: True+ library   hs-source-dirs: src   exposed-modules:@@ -29,39 +35,52 @@     Data.Poly.Semiring     Data.Poly.Orthogonal -    Data.Poly.Sparse-    Data.Poly.Sparse.Laurent-    Data.Poly.Sparse.Semiring+  if flag(sparse)+    exposed-modules:+      Data.Poly.Sparse+      Data.Poly.Sparse.Laurent+      Data.Poly.Sparse.Semiring -    Data.Poly.Multi-    Data.Poly.Multi.Laurent-    Data.Poly.Multi.Semiring-  other-modules:-    Data.Poly.Internal.Convert+      Data.Poly.Multi+      Data.Poly.Multi.Laurent+      Data.Poly.Multi.Semiring +  other-modules:     Data.Poly.Internal.Dense     Data.Poly.Internal.Dense.Field     Data.Poly.Internal.Dense.DFT     Data.Poly.Internal.Dense.GcdDomain     Data.Poly.Internal.Dense.Laurent -    Data.Poly.Internal.Multi-    Data.Poly.Internal.Multi.Core-    Data.Poly.Internal.Multi.Field-    Data.Poly.Internal.Multi.GcdDomain-    Data.Poly.Internal.Multi.Laurent+  if flag(sparse)+    other-modules:+      Data.Poly.Internal.Convert+      Data.Poly.Internal.Multi+      Data.Poly.Internal.Multi.Core+      Data.Poly.Internal.Multi.Field+      Data.Poly.Internal.Multi.GcdDomain+      Data.Poly.Internal.Multi.Laurent+   build-depends:-    base >= 4.10 && < 5,+    base >= 4.12 && < 5,     deepseq >= 1.1 && < 1.5,-    finite-typelits >= 0.1,     primitive >= 0.6,     semirings >= 0.5.2,-    vector >= 0.12.0.2,-    vector-algorithms >= 0.8.0.3,-    vector-sized >= 1.1+    vector >= 0.12.0.2++  if flag(sparse)+    build-depends:+      finite-typelits >= 0.1,+      vector-algorithms >= 0.8.0.3,+      vector-sized >= 1.1+   default-language: Haskell2010+  other-extensions: QuantifiedConstraints   ghc-options: -Wall -Wcompat -Wredundant-constraints +  if flag(sparse)+    cpp-options: -DSupportSparse+ test-suite poly-tests   type: exitcode-stdio-1.0   main-is: Main.hs@@ -69,53 +88,57 @@     Dense     DenseLaurent     DFT-    Multi-    MultiLaurent     Orthogonal     Quaternion-    Sparse-    SparseLaurent     TestUtils+  if flag(sparse)+    other-modules:+      Multi+      MultiLaurent+      Sparse+      SparseLaurent   build-depends:     base >=4.10 && <5,-    finite-typelits,     mod >=0.1.2,     poly,     QuickCheck >=2.12,+    quickcheck-classes-base,     quickcheck-classes >=0.6.3,     semirings >= 0.5.2,     tasty >= 0.11,     tasty-quickcheck >= 0.8,-    vector >= 0.12.0.2,-    vector-sized >= 1.4.2+    vector >= 0.12.0.2+  if flag(sparse)+    build-depends:+      finite-typelits,+      vector-sized >= 1.4.2   default-language: Haskell2010   hs-source-dirs: test   ghc-options: -Wall -Wcompat -threaded -rtsopts--test-suite poly-doctests-  type:             exitcode-stdio-1.0-  main-is:          doctests.hs-  hs-source-dirs:   test-  default-language: Haskell2010-  build-depends:-      base,-      poly,-      doctest+  if flag(sparse)+    cpp-options: -DSupportSparse -benchmark poly-gauge+benchmark poly-bench   build-depends:     base >=4.10 && <5,     deepseq >= 1.1 && < 1.5,-    gauge >= 0.1,     mod >=0.1.2,     poly,     semirings >= 0.2,     vector >= 0.12.0.2+  build-depends:+    tasty-bench+  mixins:+    tasty-bench (Test.Tasty.Bench as Gauge.Main)   type: exitcode-stdio-1.0   main-is: Bench.hs   other-modules:     DenseBench-    SparseBench+  if flag(sparse)+    other-modules:+      SparseBench   default-language: Haskell2010   hs-source-dirs: bench-  ghc-options: -Wall -Wcompat+  ghc-options: -Wall -Wcompat -O2 -fspecialise-aggressively+  if flag(sparse)+    cpp-options: -DSupportSparse

src/Data/Poly.hs view

@@ -6,7 +6,10 @@ -- -- Dense polynomials and a 'Num'-based interface. --+-- @since 0.1.0.0+-- +{-# LANGUAGE CPP             #-} {-# LANGUAGE PatternSynonyms #-}  module Data.Poly@@ -24,11 +27,15 @@   , deriv   , integral   , quotRemFractional+#ifdef SupportSparse   , denseToSparse   , sparseToDense+#endif   ) where +#ifdef SupportSparse import Data.Poly.Internal.Convert+#endif import Data.Poly.Internal.Dense import Data.Poly.Internal.Dense.Field (quotRemFractional) import Data.Poly.Internal.Dense.GcdDomain ()

src/Data/Poly/Internal/Convert.hs view

@@ -31,6 +31,8 @@ -- >>> :set -XFlexibleContexts -- >>> denseToSparse (1 + Data.Poly.X^2) :: Data.Poly.Sparse.UPoly Int -- 1 * X^2 + 1+--+-- @since 0.5.0.0 denseToSparse   :: (Eq a, Num a, G.Vector v a, G.Vector v (SU.Vector 1 Word, a))   => Dense.Poly v a@@ -55,6 +57,8 @@ -- >>> :set -XFlexibleContexts -- >>> sparseToDense (1 + Data.Poly.Sparse.X^2) :: Data.Poly.UPoly Int -- 1 * X^2 + 0 * X + 1+--+-- @since 0.5.0.0 sparseToDense   :: (Num a, G.Vector v a, G.Vector v (SU.Vector 1 Word, a))   => Sparse.Poly v a

src/Data/Poly/Internal/Dense.hs view

@@ -11,6 +11,7 @@ {-# LANGUAGE GeneralizedNewtypeDeriving #-} {-# LANGUAGE PatternSynonyms            #-} {-# LANGUAGE ScopedTypeVariables        #-}+{-# LANGUAGE TypeApplications           #-} {-# LANGUAGE TypeFamilies               #-} {-# LANGUAGE ViewPatterns               #-} @@ -40,18 +41,18 @@   , deriv'   , unscale'   , integral'+  , timesRing   ) where  import Prelude hiding (quotRem, quot, rem, gcd, lcm) import Control.DeepSeq (NFData) import Control.Monad-import Control.Monad.Primitive import Control.Monad.ST import Data.Bits import Data.Euclidean (Euclidean, Field, quot) import Data.Kind import Data.List (foldl', intersperse)-import Data.Semiring (Semiring(..), Ring())+import Data.Semiring (Semiring(..), Ring(), minus) import qualified Data.Semiring as Semiring import qualified Data.Vector as V import qualified Data.Vector.Generic as G@@ -62,7 +63,7 @@ -- | Polynomials of one variable with coefficients from @a@, -- backed by a 'G.Vector' @v@ (boxed, unboxed, storable, etc.). ----- Use pattern 'X' for construction:+-- Use the pattern 'X' for construction: -- -- >>> (X + 1) + (X - 1) :: VPoly Integer -- 2 * X + 0@@ -72,16 +73,27 @@ -- Polynomials are stored normalized, without leading -- zero coefficients, so 0 * 'X' + 1 equals to 1. ----- 'Ord' instance does not make much sense mathematically,+-- The 'Ord' instance does not make much sense mathematically, -- it is defined only for the sake of 'Data.Set.Set', 'Data.Map.Map', etc. --+-- Due to being polymorphic by multiple axis, the performance of `Poly` crucially+-- depends on specialisation of instances. Clients are strongly recommended+-- to compile with @ghc-options:@ @-fspecialise-aggressively@ and suggested to enable @-O2@.+--+-- @since 0.1.0.0 newtype Poly (v :: Type -> Type) (a :: Type) = Poly   { unPoly :: v a-  -- ^ Convert 'Poly' to a vector of coefficients-  -- (first element corresponds to a constant term).+  -- ^ Convert a 'Poly' to a vector of coefficients+  -- (first element corresponds to the constant term).+  --+  -- @since 0.1.0.0   }-  deriving (Eq, NFData, Ord)+  deriving+  ( Eq, Ord+  , NFData -- ^ @since 0.3.2.0+  ) +-- | @since 0.3.1.0 instance (Eq a, Semiring a, G.Vector v a) => IsList (Poly v a) where   type Item (Poly v a) = a   fromList = toPoly' . G.fromList@@ -107,60 +119,78 @@       showCoeff i c = showsPrec 7 c . showString (" * X^" ++ show i)  -- | Polynomials backed by boxed vectors.+--+-- @since 0.2.0.0 type VPoly = Poly V.Vector  -- | Polynomials backed by unboxed vectors.+--+-- @since 0.2.0.0 type UPoly = Poly U.Vector --- | Make 'Poly' from a list of coefficients--- (first element corresponds to a constant term).+-- | Make a 'Poly' from a list of coefficients+-- (first element corresponds to the constant term). -- -- >>> :set -XOverloadedLists -- >>> toPoly [1,2,3] :: VPoly Integer -- 3 * X^2 + 2 * X + 1 -- >>> toPoly [0,0,0] :: UPoly Int -- 0+--+-- @since 0.1.0.0 toPoly :: (Eq a, Num a, G.Vector v a) => v a -> Poly v a toPoly = Poly . dropWhileEnd (== 0)+{-# INLINABLE toPoly #-}  toPoly' :: (Eq a, Semiring a, G.Vector v a) => v a -> Poly v a toPoly' = Poly . dropWhileEnd (== zero)+{-# INLINABLE toPoly' #-} --- | Return a leading power and coefficient of a non-zero polynomial.+-- | Return the leading power and coefficient of a non-zero polynomial. -- -- >>> leading ((2 * X + 1) * (2 * X^2 - 1) :: UPoly Int) -- Just (3,4) -- >>> leading (0 :: UPoly Int) -- Nothing+--+-- @since 0.3.0.0 leading :: G.Vector v a => Poly v a -> Maybe (Word, a) leading (Poly v)   | G.null v  = Nothing   | otherwise = Just (fromIntegral (G.length v - 1), G.last v)+{-# INLINABLE leading #-}  -- | Note that 'abs' = 'id' and 'signum' = 'const' 1. 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+  (+) = (toPoly .) . coerce (plusPoly @v @a (+))+  (-) = (toPoly .) . coerce (minusPoly @v @a negate (-))+  (*) = (toPoly .) . coerce (inline (karatsuba @v @a 0 (+) (-) (*)))+   negate (Poly xs) = Poly $ G.map negate xs   abs = id   signum = const 1   fromInteger n = case fromInteger n of     0 -> Poly G.empty     m -> Poly $ G.singleton m-  Poly xs * Poly ys = toPoly $ karatsuba xs ys+   {-# INLINE (+) #-}   {-# INLINE (-) #-}   {-# INLINE negate #-}   {-# INLINE fromInteger #-}   {-# INLINE (*) #-} +-- | Note that 'times' is significantly slower than '(*)' for large polynomials,+-- because Karatsuba multiplication algorithm requires subtraction, which is not+-- provided by 'Semiring' class. Use 'timesRing' instead. 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 $ 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++  plus  = (toPoly' .) . coerce (plusPoly @v @a plus)+  times = (toPoly' .) . coerce (inline (convolution @v @a zero plus times))+   {-# INLINE zero #-}   {-# INLINE one #-}   {-# INLINE plus #-}@@ -174,7 +204,13 @@  instance (Eq a, Ring a, G.Vector v a) => Ring (Poly v a) where   negate (Poly xs) = Poly $ G.map Semiring.negate xs+  {-# INLINABLE negate #-} +-- | Karatsuba multiplication algorithm for polynomials over rings.+timesRing :: forall v a. (Eq a, Ring a, G.Vector v a) => Poly v a -> Poly v a -> Poly v a+timesRing = (toPoly' .) . coerce (inline (karatsuba @v @a zero plus minus times))+{-# INLINE timesRing #-}+ dropWhileEnd   :: G.Vector v a   => (a -> Bool)@@ -187,10 +223,10 @@ {-# INLINE dropWhileEnd #-}  dropWhileEndM-  :: (PrimMonad m, G.Vector v a)+  :: G.Vector v a   => (a -> Bool)-  -> G.Mutable v (PrimState m) a-  -> m (G.Mutable v (PrimState m) a)+  -> G.Mutable v s a+  -> ST s (G.Mutable v s a) dropWhileEndM p xs = go (MG.length xs)   where     go 0 = pure $ MG.unsafeSlice 0 0 xs@@ -252,50 +288,57 @@ karatsubaThreshold = 32  karatsuba-  :: (Eq a, Num a, G.Vector v a)-  => v a+  :: G.Vector v a+  => a+  -> (a -> a -> a)+  -> (a -> a -> a)+  -> (a -> a -> a)   -> v a   -> v a-karatsuba xs ys-  | lenXs <= karatsubaThreshold || lenYs <= karatsubaThreshold-  = convolution 0 (+) (*) xs ys-  | otherwise = runST $ do-    zs <- MG.unsafeNew lenZs-    forM_ [0 .. lenZs - 1] $ \k -> do-      let z0 = if k < G.length zs0-               then G.unsafeIndex zs0 k-               else 0-          z11 = if k - m >= 0 && k - m < G.length zs11-               then G.unsafeIndex zs11 (k - m)-               else 0-          z10 = if k - m >= 0 && k - m < G.length zs0-               then G.unsafeIndex zs0 (k - m)-               else 0-          z12 = if k - m >= 0 && k - m < G.length zs2-               then G.unsafeIndex zs2 (k - m)-               else 0-          z2 = if k - 2 * m >= 0 && k - 2 * m < G.length zs2-               then G.unsafeIndex zs2 (k - 2 * m)-               else 0-      MG.unsafeWrite zs k (z0 + (z11 - z10 - z12) + z2)-    G.unsafeFreeze zs+  -> v a+karatsuba zer add sub mul = go   where-    lenXs = G.length xs-    lenYs = G.length ys-    lenZs = lenXs + lenYs - 1+    conv = inline convolution zer add mul+    go xs ys+      | lenXs <= karatsubaThreshold || lenYs <= karatsubaThreshold+      = conv xs ys+      | otherwise = runST $ do+        zs <- MG.unsafeNew lenZs+        forM_ [0 .. lenZs - 1] $ \k -> do+          let z0 = if k < G.length zs0+                   then G.unsafeIndex zs0 k+                   else zer+              z11 = if k - m >= 0 && k - m < G.length zs11+                   then G.unsafeIndex zs11 (k - m)+                   else zer+              z10 = if k - m >= 0 && k - m < G.length zs0+                   then G.unsafeIndex zs0 (k - m)+                   else zer+              z12 = if k - m >= 0 && k - m < G.length zs2+                   then G.unsafeIndex zs2 (k - m)+                   else zer+              z2 = if k - 2 * m >= 0 && k - 2 * m < G.length zs2+                   then G.unsafeIndex zs2 (k - 2 * m)+                   else zer+          MG.unsafeWrite zs k (z0 `add` (z11 `sub` (z10 `add` z12)) `add` z2)+        G.unsafeFreeze zs+      where+        lenXs = G.length xs+        lenYs = G.length ys+        lenZs = lenXs + lenYs - 1 -    m    = ((lenXs `min` lenYs) + 1) `shiftR` 1+        m    = ((lenXs `min` lenYs) + 1) `shiftR` 1 -    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+        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+        xs01 = plusPoly add xs0 xs1+        ys01 = plusPoly add ys0 ys1+        zs0  = go xs0 ys0+        zs2  = go xs1 ys1+        zs11 = go xs01 ys01 {-# INLINABLE karatsuba #-}  convolution@@ -306,19 +349,21 @@   -> v a   -> v a   -> v a-convolution zer add mul xs ys-  | lenXs == 0 || lenYs == 0 = G.empty-  | otherwise = G.generate lenZs $ \k -> foldl'+convolution zer add mul = \xs ys ->+  let lenXs = G.length xs+      lenYs = G.length ys+      lenZs = lenXs + lenYs - 1 in+  if lenXs == 0 || lenYs == 0+  then G.empty+  else G.generate lenZs $ \k -> foldl'     (\acc i -> acc `add` mul (G.unsafeIndex xs i) (G.unsafeIndex ys (k - i)))     zer     [max (k - lenYs + 1) 0 .. min k (lenXs - 1)]-  where-    lenXs = G.length xs-    lenYs = G.length ys-    lenZs = lenXs + lenYs - 1 {-# INLINABLE convolution #-}  -- | Create a monomial from a power and a coefficient.+--+-- @since 0.3.0.0 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)@@ -352,6 +397,8 @@ -- -- >>> scale 2 3 (X^2 + 1) :: UPoly Int -- 3 * X^4 + 0 * X^3 + 3 * X^2 + 0 * X + 0+--+-- @since 0.3.0.0 scale :: (Eq a, Num a, G.Vector v a) => Word -> a -> Poly v a -> Poly v a scale yp yc (Poly xs) = toPoly $ scaleInternal 0 (*) yp yc xs @@ -374,10 +421,12 @@ fst' :: StrictPair a b -> a fst' (a :*: _) = a --- | Evaluate at a given point.+-- | Evaluate the polynomial at a given point. -- -- >>> eval (X^2 + 1 :: UPoly Int) 3 -- 10+--+-- @since 0.2.0.0 eval :: (Num a, G.Vector v a) => Poly v a -> a -> a eval = substitute (*) {-# INLINE eval #-}@@ -390,6 +439,8 @@ -- -- >>> subst (X^2 + 1 :: UPoly Int) (X + 1 :: UPoly Int) -- 1 * X^2 + 2 * X + 2+--+-- @since 0.3.3.0 subst :: (Eq a, Num a, G.Vector v a, G.Vector w a) => Poly v a -> Poly w a -> Poly w a subst = substitute (scale 0) {-# INLINE subst #-}@@ -408,10 +459,12 @@   G.foldl' (\(acc :*: xn) cn -> acc `plus` f cn xn :*: x `times` xn) (zero :*: one) cs {-# INLINE substitute' #-} --- | Take a derivative.+-- | Take the derivative of the polynomial. -- -- >>> deriv (X^3 + 3 * X) :: UPoly Int -- 3 * X^2 + 0 * X + 3+--+-- @since 0.2.0.0 deriv :: (Eq a, Num a, G.Vector v a) => Poly v a -> Poly v a deriv (Poly xs)   | G.null xs = Poly G.empty@@ -424,11 +477,13 @@   | otherwise = toPoly' $ G.imap (\i x -> fromNatural (fromIntegral (i + 1)) `times` x) $ G.tail xs {-# INLINE deriv' #-} --- | Compute an indefinite integral of a polynomial,--- setting constant term to zero.+-- | Compute an indefinite integral of the polynomial,+-- setting the constant term to zero. -- -- >>> integral (3 * X^2 + 3) :: UPoly Double -- 1.0 * X^3 + 0.0 * X^2 + 3.0 * X + 0.0+--+-- @since 0.2.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@@ -455,7 +510,11 @@       lenXs = G.length xs {-# INLINABLE integral' #-} --- | Create an identity polynomial.+-- | The polynomial 'X'.+--+-- > X == monomial 1 1+--+-- @since 0.2.0.0 pattern X :: (Eq a, Num a, G.Vector v a) => Poly v a pattern X <- (isVar -> True)   where X = var

src/Data/Poly/Internal/Dense/DFT.hs view

@@ -26,6 +26,8 @@  -- | <https://en.wikipedia.org/wiki/Fast_Fourier_transform Discrete Fourier transform> -- \( y_k = \sum_{j=0}^{N-1} x_j \sqrt[N]{1}^{jk} \).+--+-- @since 0.5.0.0 dft   :: (Ring a, G.Vector v a)   => a   -- ^ primitive root \( \sqrt[N]{1} \), otherwise behaviour is undefined@@ -67,9 +69,12 @@           MG.unsafeWrite ys k             $! y0 `plus`  y1           MG.unsafeWrite ys (k + halfLen) $! y0 `minus` y1         G.unsafeFreeze ys+{-# INLINABLE dft #-}  -- | Inverse <https://en.wikipedia.org/wiki/Fast_Fourier_transform discrete Fourier transform> -- \( x_k = {1\over N} \sum_{j=0}^{N-1} y_j \sqrt[N]{1}^{-jk} \).+--+-- @since 0.5.0.0 inverseDft   :: (Field a, G.Vector v a)   => a   -- ^ primitive root \( \sqrt[N]{1} \), otherwise behaviour is undefined@@ -78,3 +83,4 @@ inverseDft primRoot ys = G.map (`times` invN) $ dft (recip primRoot) ys   where     invN = recip $ fromIntegral $ G.length ys+{-# INLINABLE inverseDft #-}

src/Data/Poly/Internal/Dense/Field.hs view

@@ -4,7 +4,7 @@ -- Licence:     BSD3 -- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com> ----- GcdDomain for Field underlying+-- 'Euclidean' instance with a 'Field' constraint on the coefficient type. --  {-# LANGUAGE ConstraintKinds            #-}@@ -21,7 +21,6 @@ import Prelude hiding (quotRem, quot, rem, gcd) import Control.Exception import Control.Monad-import Control.Monad.Primitive import Control.Monad.ST import Data.Euclidean (Euclidean(..), Field) import Data.Semiring (times, minus, zero, one)@@ -32,6 +31,8 @@ import Data.Poly.Internal.Dense.GcdDomain ()  -- | Note that 'degree' 0 = 0.+--+-- @since 0.3.0.0 instance (Eq a, Field a, G.Vector v a) => Euclidean (Poly v a) where   degree (Poly xs)     | G.null xs = 0@@ -49,6 +50,8 @@ -- -- >>> quotRemFractional (X^3 + 2) (X^2 - 1 :: UPoly Double) -- (1.0 * X + 0.0,1.0 * X + 2.0)+--+-- @since 0.5.0.0 quotRemFractional :: (Eq a, Fractional a, G.Vector v a) => Poly v a -> Poly v a -> (Poly v a, Poly v a) quotRemFractional (Poly xs) (Poly ys) = (toPoly qs, toPoly rs)   where@@ -108,10 +111,10 @@ {-# INLINABLE remainder #-}  remainderM-  :: (PrimMonad m, Eq a, Field a, G.Vector v a)-  => G.Mutable v (PrimState m) a-  -> G.Mutable v (PrimState m) a-  -> m ()+  :: (Eq a, Field a, G.Vector v a)+  => G.Mutable v s a+  -> G.Mutable v s a+  -> ST s () remainderM xs ys   | lenXs < lenYs = pure ()   | lenYs == 0 = throw DivideByZero

src/Data/Poly/Internal/Dense/GcdDomain.hs view

@@ -4,7 +4,7 @@ -- Licence:     BSD3 -- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com> ----- GcdDomain for GcdDomain underlying+-- 'GcdDomain' instance with a 'GcdDomain' constraint on the coefficient type. --  {-# LANGUAGE FlexibleInstances          #-}@@ -20,7 +20,6 @@ import Prelude hiding (gcd, lcm, (^)) import Control.Exception import Control.Monad-import Control.Monad.Primitive import Control.Monad.ST import Data.Euclidean import Data.Maybe@@ -30,9 +29,11 @@  import Data.Poly.Internal.Dense +-- | @since 0.3.0.0 instance (Eq a, Ring a, GcdDomain a, G.Vector v a) => GcdDomain (Poly v a) where   divide (Poly xs) (Poly ys) =     toPoly' <$> quotient xs ys+  {-# INLINABLE divide #-}    gcd (Poly xs) (Poly ys)     | G.null xs = Poly ys@@ -45,8 +46,10 @@   lcm x@(Poly xs) y@(Poly ys)     | G.null xs || G.null ys = zero     | otherwise = (x `divide'` gcd x y) `times` y+  {-# INLINABLE lcm #-}    coprime x y = isJust (one `divide` gcd x y)+  {-# INLINABLE coprime #-}  gcdNonEmpty   :: (Eq a, Ring a, GcdDomain a, G.Vector v a)@@ -76,10 +79,10 @@ {-# INLINABLE gcdNonEmpty #-}  gcdM-  :: (PrimMonad m, Eq a, Ring a, GcdDomain a, G.Vector v a)-  => G.Mutable v (PrimState m) a-  -> G.Mutable v (PrimState m) a-  -> m (G.Mutable v (PrimState m) a)+  :: (Eq a, Ring a, GcdDomain a, G.Vector v a)+  => G.Mutable v s a+  -> G.Mutable v s a+  -> ST s (G.Mutable v s a) gcdM xs ys   | MG.null xs = pure ys   | MG.null ys = pure xs@@ -143,9 +146,9 @@ divide' = (fromMaybe (error "gcd: violated internal invariant") .) . divide  isZeroM-  :: (Eq a, Semiring a, PrimMonad m, G.Vector v a)-  => G.Mutable v (PrimState m) a-  -> m Bool+  :: (Eq a, Semiring a, G.Vector v a)+  => G.Mutable v s a+  -> ST s Bool isZeroM xs = go (MG.length xs)   where     go 0 = pure True

src/Data/Poly/Internal/Dense/Laurent.hs view

@@ -51,7 +51,7 @@ -- of one variable with coefficients from @a@, -- backed by a 'G.Vector' @v@ (boxed, unboxed, storable, etc.). ----- Use pattern 'X' and operator '^-' for construction:+-- Use the pattern 'X' and the '^-' operator for construction: -- -- >>> (X + 1) + (X^-1 - 1) :: VLaurent Integer -- 1 * X + 0 + 1 * X^-1@@ -62,9 +62,15 @@ -- and trailing -- zero coefficients, so 0 * X + 1 + 0 * X^-1 equals to 1. ----- 'Ord' instance does not make much sense mathematically,+-- The 'Ord' instance does not make much sense mathematically, -- it is defined only for the sake of 'Data.Set.Set', 'Data.Map.Map', etc. --+-- Due to being polymorphic by multiple axis, the performance of `Laurent` crucially+-- depends on specialisation of instances. Clients are strongly recommended+-- to compile with @ghc-options:@ @-fspecialise-aggressively@ and suggested to enable @-O2@.+--+-- @since 0.4.0.0+-- data Laurent (v :: Type -> Type) (a :: Type) = Laurent !Int !(Poly v a)   deriving (Eq, Ord) @@ -79,6 +85,8 @@ -- (1,2 * X + 1) -- >>> unLaurent (0 :: ULaurent Int) -- (0,0)+--+-- @since 0.4.0.0 unLaurent :: Laurent v a -> (Int, Poly v a) unLaurent (Laurent off poly) = (off, poly) @@ -89,6 +97,8 @@ -- 2 * X^3 + 1 * X^2 -- >>> toLaurent (-2) (2 * Data.Poly.X + 1) :: ULaurent Int -- 2 * X^-1 + 1 * X^-2+--+-- @since 0.4.0.0 toLaurent   :: (Eq a, Semiring a, G.Vector v a)   => Int@@ -141,19 +151,26 @@       showCoeff i c = showsPrec 7 c . showString (" * X^" ++ show i)  -- | Laurent polynomials backed by boxed vectors.+--+-- @since 0.4.0.0 type VLaurent = Laurent V.Vector  -- | Laurent polynomials backed by unboxed vectors.+--+-- @since 0.4.0.0 type ULaurent = Laurent U.Vector --- | Return a leading power and coefficient of a non-zero polynomial.+-- | Return the leading power and coefficient of a non-zero polynomial. -- -- >>> leading ((2 * X + 1) * (2 * X^2 - 1) :: ULaurent Int) -- Just (3,4) -- >>> leading (0 :: ULaurent Int) -- Nothing+--+-- @since 0.4.0.0 leading :: G.Vector v a => Laurent v a -> Maybe (Int, a) leading (Laurent off poly) = first ((+ off) . fromIntegral) <$> Dense.leading poly+{-# INLINABLE leading #-}  -- | Note that 'abs' = 'id' and 'signum' = 'const' 1. instance (Eq a, Num a, G.Vector v a) => Num (Laurent v a) where@@ -196,6 +213,8 @@   negate (Laurent off poly) = Laurent off (Semiring.negate poly)  -- | Create a monomial from a power and a coefficient.+--+-- @since 0.4.0.0 monomial :: (Eq a, Semiring a, G.Vector v a) => Int -> a -> Laurent v a monomial p c   | c == zero = Laurent 0 zero@@ -206,13 +225,18 @@ -- -- >>> scale 2 3 (X^-2 + 1) :: ULaurent Int -- 3 * X^2 + 0 * X + 3+--+-- @since 0.4.0.0 scale :: (Eq a, Semiring a, G.Vector v a) => Int -> a -> Laurent v a -> Laurent v a scale yp yc (Laurent off poly) = toLaurent (off + yp) (Dense.scale' 0 yc poly)+{-# INLINABLE scale #-} --- | Evaluate at a given point.+-- | Evaluate the polynomial at a given point. -- -- >>> eval (X^-2 + 1 :: ULaurent Double) 2 -- 1.25+--+-- @since 0.4.0.0 eval :: (Field a, G.Vector v a) => Laurent v a -> a -> a eval (Laurent off poly) x = Dense.eval' poly x `times`   (if off >= 0 then x Semiring.^ off else quot one x Semiring.^ (- off))@@ -223,20 +247,28 @@ -- >>> import Data.Poly (UPoly) -- >>> subst (Data.Poly.X^2 + 1 :: UPoly Int) (X^-1 + 1 :: ULaurent Int) -- 2 + 2 * X^-1 + 1 * X^-2+--+-- @since 0.4.0.0 subst :: (Eq a, Semiring a, G.Vector v a, G.Vector w a) => Poly v a -> Laurent w a -> Laurent w a subst = Dense.substitute' (scale 0) {-# INLINE subst #-} --- | Take a derivative.+-- | Take the derivative of the polynomial. -- -- >>> deriv (X^-1 + 3 * X) :: ULaurent Int -- 3 + 0 * X^-1 + (-1) * X^-2+--+-- @since 0.4.0.0 deriv :: (Eq a, Ring a, G.Vector v a) => Laurent v a -> Laurent v a deriv (Laurent off (Poly xs)) =   toLaurent (off - 1) $ Dense.toPoly' $ G.imap (times . Semiring.fromIntegral . (+ off)) xs {-# INLINE deriv #-} --- | Create an identity polynomial.+-- | The polynomial 'X'.+--+-- > X == monomial 1 one+--+-- @since 0.4.0.0 pattern X :: (Eq a, Semiring a, G.Vector v a) => Laurent v a pattern X <- (isVar -> True)   where X = var@@ -253,12 +285,18 @@   | otherwise          = off == 1 && G.length xs == 1 && G.unsafeHead xs == one {-# INLINE isVar #-} --- | This operator can be applied only to monomials with unit coefficients,+-- | Used to construct monomials with negative powers.+--+-- This operator can be applied only to monomials with unit coefficients, -- but is instrumental to express Laurent polynomials--- in mathematical fashion:+-- in a mathematical fashion: --+-- >>> X^-3 :: ULaurent Int+-- 1 * X^-3 -- >>> X + 2 + 3 * (X^2)^-1 :: ULaurent Int -- 1 * X + 2 + 0 * X^-1 + 3 * X^-2+--+-- @since 0.4.0.0 (^-)   :: (Eq a, Num a, G.Vector v a)   => Laurent v a

src/Data/Poly/Internal/Multi.hs view

@@ -15,6 +15,7 @@ {-# LANGUAGE PolyKinds                  #-} {-# LANGUAGE ScopedTypeVariables        #-} {-# LANGUAGE StandaloneDeriving         #-}+{-# LANGUAGE TypeApplications           #-} {-# LANGUAGE TypeFamilies               #-} {-# LANGUAGE TypeOperators              #-} {-# LANGUAGE UndecidableInstances       #-}@@ -62,6 +63,7 @@ import Prelude hiding (quot, gcd) import Control.Arrow import Control.DeepSeq+import Data.Coerce import Data.Euclidean (Field, quot) import Data.Finite import Data.Kind@@ -82,7 +84,7 @@ -- | Sparse polynomials of @n@ variables with coefficients from @a@, -- backed by a 'G.Vector' @v@ (boxed, unboxed, storable, etc.). ----- Use patterns 'Data.Poly.Multi.X',+-- Use the patterns 'Data.Poly.Multi.X', -- 'Data.Poly.Multi.Y' and -- 'Data.Poly.Multi.Z' for construction: --@@ -95,12 +97,19 @@ -- Polynomials are stored normalized, without -- zero coefficients, so 0 * 'Data.Poly.Multi.X' + 1 equals to 1. ----- 'Ord' instance does not make much sense mathematically,+-- The 'Ord' instance does not make much sense mathematically, -- it is defined only for the sake of 'Data.Set.Set', 'Data.Map.Map', etc. --+-- Due to being polymorphic by multiple axis, the performance of `MultiPoly` crucially+-- depends on specialisation of instances. Clients are strongly recommended+-- to compile with @ghc-options:@ @-fspecialise-aggressively@ and suggested to enable @-O2@.+--+-- @since 0.5.0.0 newtype MultiPoly (v :: Type -> Type) (n :: Nat) (a :: Type) = MultiPoly   { unMultiPoly :: v (SU.Vector n Word, a)-  -- ^ Convert 'MultiPoly' to a vector of (powers, coefficient) pairs.+  -- ^ Convert a 'MultiPoly' to a vector of (powers, coefficient) pairs.+  --+  -- @since 0.5.0.0   }  deriving instance Eq     (v (SU.Vector n Word, a)) => Eq     (MultiPoly v n a)@@ -108,9 +117,13 @@ deriving instance NFData (v (SU.Vector n Word, a)) => NFData (MultiPoly v n a)  -- | Multivariate polynomials backed by boxed vectors.+--+-- @since 0.5.0.0 type VMultiPoly (n :: Nat) (a :: Type) = MultiPoly V.Vector n a  -- | Multivariate polynomials backed by unboxed vectors.+--+-- @since 0.5.0.0 type UMultiPoly (n :: Nat) (a :: Type) = MultiPoly U.Vector n a  -- | Sparse univariate polynomials with coefficients from @a@,@@ -129,15 +142,26 @@ -- 'Ord' instance does not make much sense mathematically, -- it is defined only for the sake of 'Data.Set.Set', 'Data.Map.Map', etc. --+-- Due to being polymorphic by multiple axis, the performance of `Poly` crucially+-- depends on specialisation of instances. Clients are strongly recommended+-- to compile with @ghc-options:@ @-fspecialise-aggressively@ and suggested to enable @-O2@.+--+-- @since 0.3.0.0 type Poly (v :: Type -> Type) (a :: Type) = MultiPoly v 1 a  -- | Polynomials backed by boxed vectors.+--+-- @since 0.3.0.0 type VPoly (a :: Type) = Poly V.Vector a  -- | Polynomials backed by unboxed vectors.+--+-- @since 0.3.0.0 type UPoly (a :: Type) = Poly U.Vector a --- | Convert 'Poly' to a vector of coefficients.+-- | Convert a 'Poly' to a vector of coefficients.+--+-- @since 0.3.0.0 unPoly   :: (G.Vector v (Word, a), G.Vector v (SU.Vector 1 Word, a))   => Poly v a@@ -178,7 +202,7 @@         2 -> "Z"         k -> "X" ++ show k --- | Make 'MultiPoly' from a list of (powers, coefficient) pairs.+-- | Make a 'MultiPoly' from a list of (powers, coefficient) pairs. -- -- >>> :set -XOverloadedLists -XDataKinds -- >>> import Data.Vector.Generic.Sized (fromTuple)@@ -186,29 +210,36 @@ -- 3 * X + 2 * Y + 1 -- >>> toMultiPoly [(fromTuple (0,0),0),(fromTuple (0,1),0),(fromTuple (1,0),0)] :: UMultiPoly 2 Int -- 0+--+-- @since 0.5.0.0 toMultiPoly   :: (Eq a, Num a, G.Vector v (SU.Vector n Word, a))   => v (SU.Vector n Word, a)   -> MultiPoly v n a toMultiPoly = MultiPoly . normalize (/= 0) (+)+{-# INLINABLE toMultiPoly #-}  toMultiPoly'   :: (Eq a, Semiring a, G.Vector v (SU.Vector n Word, a))   => v (SU.Vector n Word, a)   -> MultiPoly v n a toMultiPoly' = MultiPoly . normalize (/= zero) plus+{-# INLINABLE toMultiPoly' #-}  -- | Note that 'abs' = 'id' and 'signum' = 'const' 1. instance (Eq a, Num a, KnownNat n, G.Vector v (SU.Vector n Word, a)) => Num (MultiPoly v n a) where-  MultiPoly xs + MultiPoly ys = MultiPoly $ plusPoly (/= 0) (+) xs ys-  MultiPoly xs - MultiPoly ys = MultiPoly $ minusPoly (/= 0) negate (-) xs ys++  (+) = coerce (plusPoly    @v @(SU.Vector n Word) @a (/= 0) (+))+  (-) = coerce (minusPoly   @v @(SU.Vector n Word) @a (/= 0) negate (-))+  (*) = coerce (convolution @v @(SU.Vector n Word) @a (/= 0) (+) (*))+   negate (MultiPoly xs) = MultiPoly $ G.map (fmap negate) xs   abs = id   signum = const 1   fromInteger n = case fromInteger n of     0 -> MultiPoly G.empty     m -> MultiPoly $ G.singleton (0, m)-  MultiPoly xs * MultiPoly ys = MultiPoly $ convolution (/= 0) (+) (*) xs ys+   {-# INLINE (+) #-}   {-# INLINE (-) #-}   {-# INLINE negate #-}@@ -220,8 +251,10 @@   one     | (one :: a) == zero = zero     | otherwise = MultiPoly $ G.singleton (0, one)-  plus (MultiPoly xs) (MultiPoly ys) = MultiPoly $ plusPoly (/= zero) plus xs ys-  times (MultiPoly xs) (MultiPoly ys) = MultiPoly $ convolution (/= zero) plus times xs ys++  plus  = coerce (plusPoly    @v @(SU.Vector n Word) @a (/= zero) plus)+  times = coerce (convolution @v @(SU.Vector n Word) @a (/= zero) plus times)+   {-# INLINE zero #-}   {-# INLINE one #-}   {-# INLINE plus #-}@@ -236,13 +269,15 @@ instance (Eq a, Ring a, KnownNat n, G.Vector v (SU.Vector n Word, a)) => Ring (MultiPoly v n a) where   negate (MultiPoly xs) = MultiPoly $ G.map (fmap Semiring.negate) xs --- | Return a leading power and coefficient of a non-zero polynomial.+-- | Return the leading power and coefficient of a non-zero polynomial. -- -- >>> import Data.Poly.Sparse (UPoly) -- >>> leading ((2 * X + 1) * (2 * X^2 - 1) :: UPoly Int) -- Just (3,4) -- >>> leading (0 :: UPoly Int) -- Nothing+--+-- @since 0.3.0.0 leading :: G.Vector v (SU.Vector 1 Word, a) => Poly v a -> Maybe (Word, a) leading (MultiPoly v)   | G.null v  = Nothing@@ -254,6 +289,8 @@ -- >>> import Data.Vector.Generic.Sized (fromTuple) -- >>> scale (fromTuple (1, 1)) 3 (X^2 + Y) :: UMultiPoly 2 Int -- 3 * X^3 * Y + 3 * X * Y^2+--+-- @since 0.5.0.0 scale   :: (Eq a, Num a, KnownNat n, G.Vector v (SU.Vector n Word, a))   => SU.Vector n Word@@ -271,6 +308,8 @@ scale' yp yc = MultiPoly . scaleInternal (/= zero) times yp yc . unMultiPoly  -- | Create a monomial from powers and a coefficient.+--+-- @since 0.5.0.0 monomial   :: (Eq a, Num a, G.Vector v (SU.Vector n Word, a))   => SU.Vector n Word@@ -279,6 +318,7 @@ monomial p c   | c == 0    = MultiPoly G.empty   | otherwise = MultiPoly $ G.singleton (p, c)+{-# INLINABLE monomial #-}  monomial'   :: (Eq a, Semiring a, G.Vector v (SU.Vector n Word, a))@@ -288,13 +328,16 @@ monomial' p c   | c == zero = MultiPoly G.empty   | otherwise = MultiPoly $ G.singleton (p, c)+{-# INLINABLE monomial' #-} --- | Evaluate at a given point.+-- | Evaluate the polynomial at a given point. -- -- >>> :set -XDataKinds -- >>> import Data.Vector.Generic.Sized (fromTuple) -- >>> eval (X^2 + Y^2 :: UMultiPoly 2 Int) (fromTuple (3, 4) :: Data.Vector.Sized.Vector 2 Int) -- 25+--+-- @since 0.5.0.0 eval   :: (Num a, G.Vector v (SU.Vector n Word, a), G.Vector u a)   => MultiPoly v n a@@ -311,12 +354,14 @@ eval' = substitute' times {-# INLINE eval' #-} --- | Substitute another polynomials instead of variables.+-- | Substitute other polynomials instead of the variables. -- -- >>> :set -XDataKinds -- >>> import Data.Vector.Generic.Sized (fromTuple) -- >>> subst (X^2 + Y^2 + Z^2 :: UMultiPoly 3 Int) (fromTuple (X + 1, Y + 1, X + Y :: UMultiPoly 2 Int)) -- 2 * X^2 + 2 * X * Y + 2 * X + 2 * Y^2 + 2 * Y + 2+--+-- @since 0.5.0.0 subst   :: (Eq a, Num a, KnownNat m, G.Vector v (SU.Vector n Word, a), G.Vector w (SU.Vector m Word, a))   => MultiPoly v n a@@ -365,13 +410,15 @@     doMonom = SU.ifoldl' (\acc i p -> acc `times` ((xs `SG.index` i) Semiring.^ p)) one {-# INLINE substitute' #-} --- | Take a derivative with respect to the /i/-th variable.+-- | Take the derivative of the polynomial with respect to the /i/-th variable. -- -- >>> :set -XDataKinds -- >>> deriv 0 (X^3 + 3 * Y) :: UMultiPoly 2 Int -- 3 * X^2 -- >>> deriv 1 (X^3 + 3 * Y) :: UMultiPoly 2 Int -- 3+--+-- @since 0.5.0.0 deriv   :: (Eq a, Num a, G.Vector v (SU.Vector n Word, a))   => Finite n@@ -396,15 +443,17 @@   xs {-# INLINE deriv' #-} --- | Compute an indefinite integral of a polynomial--- by the /i/-th variable,--- setting constant term to zero.+-- | Compute an indefinite integral of the polynomial+-- with respect to the /i/-th variable,+-- setting the constant term to zero. -- -- >>> :set -XDataKinds -- >>> integral 0 (3 * X^2 + 2 * Y) :: UMultiPoly 2 Double -- 1.0 * X^3 + 2.0 * X * Y -- >>> integral 1 (3 * X^2 + 2 * Y) :: UMultiPoly 2 Double -- 3.0 * X^2 * Y + 1.0 * Y^2+--+-- @since 0.5.0.0 integral   :: (Fractional a, G.Vector v (SU.Vector n Word, a))   => Finite n@@ -428,6 +477,8 @@ {-# INLINE integral' #-}  -- | Create a polynomial equal to the first variable.+--+-- @since 0.5.0.0 pattern X   :: (Eq a, Num a, KnownNat n, 1 <= n, G.Vector v (SU.Vector n Word, a))   => MultiPoly v n a@@ -441,6 +492,8 @@   where X' = var' 0  -- | Create a polynomial equal to the second variable.+--+-- @since 0.5.0.0 pattern Y   :: (Eq a, Num a, KnownNat n, 2 <= n, G.Vector v (SU.Vector n Word, a))   => MultiPoly v n a@@ -454,6 +507,8 @@   where Y' = var' 1  -- | Create a polynomial equal to the third variable.+--+-- @since 0.5.0.0 pattern Z   :: (Eq a, Num a, KnownNat n, 3 <= n, G.Vector v (SU.Vector n Word, a))   => MultiPoly v n a@@ -527,6 +582,8 @@ -- | Interpret a multivariate polynomial over 1+/m/ variables -- as a univariate polynomial, whose coefficients are -- multivariate polynomials over the last /m/ variables.+--+-- @since 0.5.0.0 segregate   :: (G.Vector v (SU.Vector (1 + m) Word, a), G.Vector v (SU.Vector m Word, a))   => MultiPoly v (1 + m) a@@ -541,6 +598,8 @@ -- | Interpret a univariate polynomials, whose coefficients are -- multivariate polynomials over the first /m/ variables, -- as a multivariate polynomial over 1+/m/ variables.+--+-- @since 0.5.0.0 unsegregate   :: (G.Vector v (SU.Vector (1 + m) Word, a), G.Vector v (SU.Vector m Word, a))   => VPoly (MultiPoly v m a)

src/Data/Poly/Internal/Multi/Core.hs view

@@ -23,7 +23,6 @@   ) where  import Control.Monad-import Control.Monad.Primitive import Control.Monad.ST import Data.Bits import Data.Ord@@ -44,13 +43,14 @@     ws <- G.thaw vs     l' <- normalizeM p add ws     G.unsafeFreeze $ MG.unsafeSlice 0 l' ws+{-# INLINABLE normalize #-}  normalizeM-  :: (PrimMonad m, G.Vector v (t, a), Ord t)+  :: (G.Vector v (t, a), Ord t)   => (a -> Bool)   -> (a -> a -> a)-  -> G.Mutable v (PrimState m) (t, a)-  -> m Int+  -> G.Mutable v s (t, a)+  -> ST s Int normalizeM p add ws = do     let l = MG.length ws     let go i j acc@(accP, accC)@@ -72,6 +72,7 @@     Tim.sortBy (comparing fst) ws     wsHead <- MG.unsafeRead ws 0     go 0 1 wsHead+{-# INLINABLE normalizeM #-}  plusPoly   :: (G.Vector v (t, a), Ord t)@@ -80,20 +81,20 @@   -> v (t, a)   -> v (t, a)   -> v (t, a)-plusPoly p add xs ys = runST $ do+plusPoly p add = \xs ys -> runST $ do   zs <- MG.unsafeNew (G.length xs + G.length ys)   lenZs <- plusPolyM p add xs ys zs   G.unsafeFreeze $ MG.unsafeSlice 0 lenZs zs {-# INLINABLE plusPoly #-}  plusPolyM-  :: (PrimMonad m, G.Vector v (t, a), Ord t)+  :: (G.Vector v (t, a), Ord t)   => (a -> Bool)   -> (a -> a -> a)   -> v (t, a)   -> v (t, a)-  -> G.Mutable v (PrimState m) (t, a)-  -> m Int+  -> G.Mutable v s (t, a)+  -> ST s Int plusPolyM p add xs ys zs = go 0 0 0   where     lenXs = G.length xs@@ -127,7 +128,7 @@         GT -> do           MG.unsafeWrite zs iz (yp, yc)           go ix (iy + 1) (iz + 1)-{-# INLINABLE plusPolyM #-}+{-# INLINE plusPolyM #-}  minusPoly   :: (G.Vector v (t, a), Ord t)@@ -137,7 +138,9 @@   -> v (t, a)   -> v (t, a)   -> v (t, a)-minusPoly p neg sub xs ys = runST $ do+minusPoly p neg sub = \xs ys -> runST $ do+  let lenXs = G.length xs+      lenYs = G.length ys   zs <- MG.unsafeNew (lenXs + lenYs)   let go ix iy iz         | ix == lenXs, iy == lenYs = pure iz@@ -169,19 +172,16 @@             go ix (iy + 1) (iz + 1)   lenZs <- go 0 0 0   G.unsafeFreeze $ MG.unsafeSlice 0 lenZs zs-  where-    lenXs = G.length xs-    lenYs = G.length ys {-# INLINABLE minusPoly #-}  scaleM-  :: (PrimMonad m, G.Vector v (t, a), Num t)+  :: (G.Vector v (t, a), Num t)   => (a -> Bool)   -> (a -> a -> a)   -> v (t, a)   -> (t, a)-  -> G.Mutable v (PrimState m) (t, a)-  -> m Int+  -> G.Mutable v s (t, a)+  -> ST s Int scaleM p mul xs (yp, yc) zs = go 0 0   where     lenXs = G.length xs@@ -221,11 +221,10 @@   -> v (t, a)   -> v (t, a)   -> v (t, a)-convolution p add mult xs ys-  | G.length xs >= G.length ys-  = go mult xs ys-  | otherwise-  = go (flip mult) ys xs+convolution p add mult = \xs ys ->+  if G.length xs >= G.length ys+  then go mult xs ys+  else go (flip mult) ys xs   where     go :: (a -> a -> a) -> v (t, a) -> v (t, a) -> v (t, a)     go mul long short = runST $ do@@ -248,11 +247,10 @@       gogo slices' buffer' bufferNew      gogo-      :: PrimMonad m-      => U.Vector (Int, Int)+      :: U.Vector (Int, Int)       -> v (t, a)-      -> G.Mutable v (PrimState m) (t, a)-      -> m (v (t, a))+      -> G.Mutable v s (t, a)+      -> ST s (v (t, a))     gogo slices buffer bufferNew       | G.length slices == 0       = pure G.empty

src/Data/Poly/Internal/Multi/Field.hs view

@@ -4,7 +4,7 @@ -- Licence:     BSD3 -- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com> ----- Euclidean for Field underlying+-- 'Euclidean' instance with a 'Field' constraint on the coefficient type. --  {-# LANGUAGE ConstraintKinds            #-}@@ -44,6 +44,8 @@ -- -- >>> quotRemFractional (X^3 + 2) (X^2 - 1 :: UPoly Double) -- (1.0 * X,1.0 * X + 2.0)+--+-- @since 0.5.0.0 quotRemFractional :: (Eq a, Fractional a, G.Vector v (SU.Vector 1 Word, a)) => Poly v a -> Poly v a -> (Poly v a, Poly v a) quotRemFractional = quotientRemainder 0 (+) (-) (*) (/) {-# INLINE quotRemFractional #-}

src/Data/Poly/Internal/Multi/GcdDomain.hs view

@@ -4,7 +4,7 @@ -- Licence:     BSD3 -- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com> ----- GcdDomain for GcdDomain underlying+-- 'GcdDomain' instance with a 'GcdDomain' constraint on the coefficient type. --  {-# LANGUAGE CPP                        #-}@@ -12,15 +12,12 @@ {-# LANGUAGE FlexibleContexts           #-} {-# LANGUAGE FlexibleInstances          #-} {-# LANGUAGE GADTs                      #-}+{-# LANGUAGE QuantifiedConstraints      #-} {-# LANGUAGE ScopedTypeVariables        #-} {-# LANGUAGE TypeOperators              #-} {-# LANGUAGE TypeFamilies               #-} {-# LANGUAGE UndecidableInstances       #-} -#if __GLASGOW_HASKELL__ >= 806-{-# LANGUAGE QuantifiedConstraints      #-}-#endif- {-# OPTIONS_GHC -fno-warn-orphans #-}  module Data.Poly.Internal.Multi.GcdDomain@@ -40,10 +37,6 @@  import Data.Poly.Internal.Multi -#if __GLASGOW_HASKELL__ < 806-import qualified Data.Vector as V-#endif- instance {-# OVERLAPPING #-} (Eq a, Ring a, GcdDomain a, G.Vector v (SU.Vector 1 Word, a)) => GcdDomain (Poly v a) where   divide xs ys     | G.null (unMultiPoly ys) = throw DivideByZero@@ -71,11 +64,7 @@ isSucc = case someNatVal (natVal (Proxy :: Proxy n) - 1) of   SomeNat (_ :: Proxy m) -> IsSucc (unsafeCoerce Refl :: n :~: 1 + m) -#if __GLASGOW_HASKELL__ >= 806 instance (Eq a, Ring a, GcdDomain a, KnownNat n, forall m. KnownNat m => G.Vector v (SU.Vector m Word, a), forall m. KnownNat m => Eq (v (SU.Vector m Word, a))) => GcdDomain (MultiPoly v n a) where-#else-instance (Eq a, Ring a, GcdDomain a, KnownNat n, v ~ V.Vector) => GcdDomain (MultiPoly v n a) where-#endif   divide xs ys     | G.null (unMultiPoly ys) = throw DivideByZero     | G.length (unMultiPoly ys) == 1 = divideSingleton xs (G.unsafeHead (unMultiPoly ys))

src/Data/Poly/Internal/Multi/Laurent.hs view

@@ -14,6 +14,7 @@ {-# LANGUAGE FlexibleInstances          #-} {-# LANGUAGE LambdaCase                 #-} {-# LANGUAGE PatternSynonyms            #-}+{-# LANGUAGE QuantifiedConstraints      #-} {-# LANGUAGE ScopedTypeVariables        #-} {-# LANGUAGE StandaloneDeriving         #-} {-# LANGUAGE TypeFamilies               #-}@@ -21,10 +22,6 @@ {-# LANGUAGE UndecidableInstances       #-} {-# LANGUAGE ViewPatterns               #-} -#if __GLASGOW_HASKELL__ >= 806-{-# LANGUAGE QuantifiedConstraints      #-}-#endif- module Data.Poly.Internal.Multi.Laurent   ( MultiLaurent   , VMultiLaurent@@ -82,7 +79,7 @@ -- of @n@ variables with coefficients from @a@, -- backed by a 'G.Vector' @v@ (boxed, unboxed, storable, etc.). ----- Use patterns 'X', 'Y', 'Z' and operator '^-' for construction:+-- Use the patterns 'X', 'Y', 'Z' and the '^-' operator for construction: -- -- >>> (X + 1) + (Y^-1 - 1) :: VMultiLaurent 2 Integer -- 1 * X + 1 * Y^-1@@ -92,9 +89,14 @@ -- Polynomials are stored normalized, without -- zero coefficients, so 0 * X + 1 + 0 * X^-1 equals to 1. ----- 'Ord' instance does not make much sense mathematically,+-- The 'Ord' instance does not make much sense mathematically, -- it is defined only for the sake of 'Data.Set.Set', 'Data.Map.Map', etc. --+-- Due to being polymorphic by multiple axis, the performance of `MultiLaurent` crucially+-- depends on specialisation of instances. Clients are strongly recommended+-- to compile with @ghc-options:@ @-fspecialise-aggressively@ and suggested to enable @-O2@.+--+-- @since 0.5.0.0 data MultiLaurent (v :: Type -> Type) (n :: Nat) (a :: Type) =   MultiLaurent !(SU.Vector n Int) !(MultiPoly v n a) @@ -102,16 +104,20 @@ deriving instance Ord (v (SU.Vector n Word, a)) => Ord (MultiLaurent v n a)  -- | Multivariate Laurent polynomials backed by boxed vectors.+--+-- @since 0.5.0.0 type VMultiLaurent (n :: Nat) (a :: Type) = MultiLaurent V.Vector n a  -- | Multivariate Laurent polynomials backed by unboxed vectors.+--+-- @since 0.5.0.0 type UMultiLaurent (n :: Nat) (a :: Type) = MultiLaurent U.Vector n a  -- | <https://en.wikipedia.org/wiki/Laurent_polynomial Laurent polynomials> -- of one variable with coefficients from @a@, -- backed by a 'G.Vector' @v@ (boxed, unboxed, storable, etc.). ----- Use pattern 'X' and operator '^-' for construction:+-- Use the pattern 'X' and the '^-' operator for construction: -- -- >>> (X + 1) + (X^-1 - 1) :: VLaurent Integer -- 1 * X + 1 * X^-1@@ -121,15 +127,24 @@ -- Polynomials are stored normalized, without -- zero coefficients, so 0 * X + 1 + 0 * X^-1 equals to 1. ----- 'Ord' instance does not make much sense mathematically,+-- The 'Ord' instance does not make much sense mathematically, -- it is defined only for the sake of 'Data.Set.Set', 'Data.Map.Map', etc. --+-- Due to being polymorphic by multiple axis, the performance of `Laurent` crucially+-- depends on specialisation of instances. Clients are strongly recommended+-- to compile with @ghc-options:@ @-fspecialise-aggressively@ and suggested to enable @-O2@.+--+-- @since 0.4.0.0 type Laurent (v :: Type -> Type) (a :: Type) = MultiLaurent v 1 a  -- | Laurent polynomials backed by boxed vectors.+--+-- @since 0.4.0.0 type VLaurent (a :: Type) = Laurent V.Vector a  -- | Laurent polynomials backed by unboxed vectors.+--+-- @since 0.4.0.0 type ULaurent (a :: Type) = Laurent U.Vector a  instance (Eq a, Semiring a, KnownNat n, G.Vector v (SU.Vector n Int, a), G.Vector v (SU.Vector n Word, a)) => IsList (MultiLaurent v n a) where@@ -155,6 +170,8 @@ -- (Vector [1,0],2 * X + 1) -- >>> unMultiLaurent (0 :: UMultiLaurent 2 Int) -- (Vector [0,0],0)+--+-- @since 0.5.0.0 unMultiLaurent :: MultiLaurent v n a -> (SU.Vector n Int, MultiPoly v n a) unMultiLaurent (MultiLaurent off poly) = (off, poly) @@ -169,10 +186,12 @@ -- (1,2 * X + 1) -- >>> unLaurent (0 :: ULaurent Int) -- (0,0)+--+-- @since 0.4.0.0 unLaurent :: Laurent v a -> (Int, Poly v a) unLaurent = first SU.head . unMultiLaurent --- | Construct 'MultiLaurent' polynomial from an offset and a regular polynomial.+-- | Construct a 'MultiLaurent' polynomial from an offset and a regular polynomial. -- One can imagine it as 'Data.Poly.Multi.Semiring.scale', but allowing negative offsets. -- -- >>> :set -XDataKinds@@ -196,19 +215,22 @@         | otherwise = G.map (first (SU.zipWith subtract minPow)) xs {-# INLINE toMultiLaurent #-} --- | Construct 'Laurent' polynomial from an offset and a regular polynomial.+-- | Construct a 'Laurent' polynomial from an offset and a regular polynomial. -- One can imagine it as 'Data.Poly.Sparse.Semiring.scale', but allowing negative offsets. -- -- >>> toLaurent 2 (2 * Data.Poly.Sparse.X + 1) :: ULaurent Int -- 2 * X^3 + 1 * X^2 -- >>> toLaurent (-2) (2 * Data.Poly.Sparse.X + 1) :: ULaurent Int -- 2 * X^-1 + 1 * X^-2+--+-- @since 0.4.0.0 toLaurent   :: G.Vector v (SU.Vector 1 Word, a)   => Int   -> Poly v a   -> Laurent v a toLaurent = toMultiLaurent . SU.singleton+{-# INLINABLE toLaurent #-}  instance NFData (v (SU.Vector n Word, a)) => NFData (MultiLaurent v n a) where   rnf (MultiLaurent off poly) = rnf off `seq` rnf poly@@ -241,12 +263,14 @@         2 -> "Z"         k -> "X" ++ show k --- | Return a leading power and coefficient of a non-zero polynomial.+-- | Return the leading power and coefficient of a non-zero polynomial. -- -- >>> leading ((2 * X + 1) * (2 * X^2 - 1) :: ULaurent Int) -- Just (3,4) -- >>> leading (0 :: ULaurent Int) -- Nothing+--+-- @since 0.4.0.0 leading :: G.Vector v (SU.Vector 1 Word, a) => Laurent v a -> Maybe (Int, a) leading (MultiLaurent off poly) = first ((+ SU.head off) . fromIntegral) <$> Multi.leading poly @@ -294,6 +318,8 @@   negate (MultiLaurent off poly) = MultiLaurent off (Semiring.negate poly)  -- | Create a monomial from a power and a coefficient.+--+-- @since 0.5.0.0 monomial   :: (Eq a, Semiring a, KnownNat n, G.Vector v (SU.Vector n Word, a))   => SU.Vector n Int@@ -310,6 +336,8 @@ -- >>> import Data.Vector.Generic.Sized (fromTuple) -- >>> scale (fromTuple (1, 1)) 3 (X^-2 + Y) :: UMultiLaurent 2 Int -- 3 * X * Y^2 + 3 * X^-1 * Y+--+-- @since 0.5.0.0 scale   :: (Eq a, Semiring a, KnownNat n, G.Vector v (SU.Vector n Word, a))   => SU.Vector n Int@@ -318,12 +346,14 @@   -> MultiLaurent v n a scale yp yc (MultiLaurent off poly) = toMultiLaurent (off + yp) (Multi.scale' 0 yc poly) --- | Evaluate at a given point.+-- | Evaluate the polynomial at a given point. -- -- >>> :set -XDataKinds -- >>> import Data.Vector.Generic.Sized (fromTuple) -- >>> eval (X^2 + Y^-1 :: UMultiLaurent 2 Double) (fromTuple (3, 4) :: Data.Vector.Sized.Vector 2 Double) -- 9.25+--+-- @since 0.5.0.0 eval   :: (Field a, G.Vector v (SU.Vector n Word, a), G.Vector u a)   => MultiLaurent v n a@@ -340,6 +370,8 @@ -- >>> import Data.Poly.Multi (UMultiPoly) -- >>> subst (Data.Poly.Multi.X * Data.Poly.Multi.Y :: UMultiPoly 2 Int) (fromTuple (X + Y^-1, Y + X^-1 :: UMultiLaurent 2 Int)) -- 1 * X * Y + 2 + 1 * X^-1 * Y^-1+--+-- @since 0.5.0.0 subst   :: (Eq a, Semiring a, KnownNat n, G.Vector v (SU.Vector n Word, a), G.Vector w (SU.Vector n Word, a))   => MultiPoly v n a@@ -348,13 +380,15 @@ subst = Multi.substitute' (scale 0) {-# INLINE subst #-} --- | Take a derivative with respect to the /i/-th variable.+-- | Take the derivative of the polynomial with respect to the /i/-th variable. -- -- >>> :set -XDataKinds -- >>> deriv 0 (X^3 + 3 * Y) :: UMultiLaurent 2 Int -- 3 * X^2 -- >>> deriv 1 (X^3 + 3 * Y) :: UMultiLaurent 2 Int -- 3+--+-- @since 0.5.0.0 deriv   :: (Eq a, Ring a, KnownNat n, G.Vector v (SU.Vector n Word, a))   => Finite n@@ -369,6 +403,8 @@ {-# INLINE deriv #-}  -- | Create a polynomial equal to the first variable.+--+-- @since 0.5.0.0 pattern X   :: (Eq a, Semiring a, KnownNat n, 1 <= n, G.Vector v (SU.Vector n Word, a))   => MultiLaurent v n a@@ -376,6 +412,8 @@   where X = var 0  -- | Create a polynomial equal to the second variable.+--+-- @since 0.5.0.0 pattern Y   :: (Eq a, Semiring a, KnownNat n, 2 <= n, G.Vector v (SU.Vector n Word, a))   => MultiLaurent v n a@@ -383,6 +421,8 @@   where Y = var 1  -- | Create a polynomial equal to the third variable.+--+-- @since 0.5.0.0 pattern Z   :: (Eq a, Semiring a, KnownNat n, 3 <= n, G.Vector v (SU.Vector n Word, a))   => MultiLaurent v n a@@ -414,12 +454,18 @@   && G.length xs == 1 && G.unsafeHead xs == (0, one) {-# INLINE isVar #-} --- | This operator can be applied only to monomials with unit coefficients,--- but is still instrumental to express Laurent polynomials--- in mathematical fashion:+-- | Used to construct monomials with negative powers. --+-- This operator can be applied only to monomials with unit coefficients,+-- but is instrumental to express Laurent polynomials+-- in a mathematical fashion:+--+-- >>> X^-3 * Y^-1 :: UMultiLaurent 2 Int+-- 1 * X^-3 * Y^-1 -- >>> 3 * X^-1 + 2 * (Y^2)^-2 :: UMultiLaurent 2 Int -- 2 * Y^-4 + 3 * X^-1+--+-- @since 0.5.0.0 (^-)   :: (Eq a, Semiring a, KnownNat n, G.Vector v (SU.Vector n Word, a))   => MultiLaurent v n a@@ -448,11 +494,7 @@     coprime poly1 poly2   {-# INLINE coprime #-} -#if __GLASGOW_HASKELL__ >= 806 instance (Eq a, Ring a, GcdDomain a, KnownNat n, forall m. KnownNat m => G.Vector v (SU.Vector m Word, a), forall m. KnownNat m => Eq (v (SU.Vector m Word, a))) => GcdDomain (MultiLaurent v n a) where-#else-instance (Eq a, Ring a, GcdDomain a, KnownNat n, v ~ V.Vector) => GcdDomain (MultiLaurent v n a) where-#endif   divide (MultiLaurent off1 poly1) (MultiLaurent off2 poly2) =     toMultiLaurent (off1 - off2) <$> divide poly1 poly2   {-# INLINE divide #-}@@ -474,6 +516,8 @@ -- | Interpret a multivariate Laurent polynomial over 1+/m/ variables -- as a univariate Laurent polynomial, whose coefficients are -- multivariate Laurent polynomials over the last /m/ variables.+--+-- @since 0.5.0.0 segregate   :: (KnownNat m, G.Vector v (SU.Vector (1 + m) Word, a), G.Vector v (SU.Vector m Word, a))   => MultiLaurent v (1 + m) a@@ -488,6 +532,8 @@ -- | Interpret a univariate Laurent polynomials, whose coefficients are -- multivariate Laurent polynomials over the first /m/ variables, -- as a multivariate polynomial over 1+/m/ variables.+--+-- @since 0.5.0.0 unsegregate   :: forall v m a.      (KnownNat m, KnownNat (1 + m), G.Vector v (SU.Vector (1 + m) Word, a), G.Vector v (SU.Vector m Word, a))

src/Data/Poly/Laurent.hs view

@@ -6,6 +6,8 @@ -- -- <https://en.wikipedia.org/wiki/Laurent_polynomial Laurent polynomials>. --+-- @since 0.4.0.0+--  {-# LANGUAGE PatternSynonyms            #-}

src/Data/Poly/Multi.hs view

@@ -6,6 +6,7 @@ -- -- Sparse multivariate polynomials with 'Num' instance. --+-- @since 0.5.0.0  {-# LANGUAGE DataKinds        #-} {-# LANGUAGE FlexibleContexts #-}

src/Data/Poly/Multi/Semiring.hs view

@@ -4,8 +4,9 @@ -- Licence:     BSD3 -- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com> ----- Sparse multivariate polynomials with 'Semiring' instance.+-- Sparse multivariate polynomials with a 'Semiring' instance. --+-- @since 0.5.0.0  {-# LANGUAGE DataKinds        #-} {-# LANGUAGE FlexibleContexts #-}@@ -46,7 +47,7 @@ import Data.Poly.Internal.Multi.Field () import Data.Poly.Internal.Multi.GcdDomain () --- | Make 'MultiPoly' from a list of (powers, coefficient) pairs.+-- | Make a 'MultiPoly' from a list of (powers, coefficient) pairs. -- -- >>> :set -XOverloadedLists -XDataKinds -- >>> import Data.Vector.Generic.Sized (fromTuple)@@ -54,6 +55,8 @@ -- 3 * X + 2 * Y + 1 -- >>> toMultiPoly [(fromTuple (0,0),0),(fromTuple (0,1),0),(fromTuple (1,0),0)] :: UMultiPoly 2 Int -- 0+--+-- @since 0.5.0.0 toMultiPoly   :: (Eq a, Semiring a, G.Vector v (SU.Vector n Word, a))   => v (SU.Vector n Word, a)@@ -61,6 +64,8 @@ toMultiPoly = Multi.toMultiPoly'  -- | Create a monomial from powers and a coefficient.+--+-- @since 0.5.0.0 monomial   :: (Eq a, Semiring a, G.Vector v (SU.Vector n Word, a))   => SU.Vector n Word@@ -74,6 +79,8 @@ -- >>> import Data.Vector.Generic.Sized (fromTuple) -- >>> scale (fromTuple (1, 1)) 3 (X^2 + Y) :: UMultiPoly 2 Int -- 3 * X^3 * Y + 3 * X * Y^2+--+-- @since 0.5.0.0 scale   :: (Eq a, Semiring a, KnownNat n, G.Vector v (SU.Vector n Word, a))   => SU.Vector n Word@@ -83,29 +90,37 @@ scale = Multi.scale'  -- | Create a polynomial equal to the first variable.+--+-- @since 0.5.0.0 pattern X   :: (Eq a, Semiring a, KnownNat n, 1 <= n, G.Vector v (SU.Vector n Word, a))   => MultiPoly v n a pattern X = Multi.X'  -- | Create a polynomial equal to the second variable.+--+-- @since 0.5.0.0 pattern Y   :: (Eq a, Semiring a, KnownNat n, 2 <= n, G.Vector v (SU.Vector n Word, a))   => MultiPoly v n a pattern Y = Multi.Y'  -- | Create a polynomial equal to the third variable.+--+-- @since 0.5.0.0 pattern Z   :: (Eq a, Semiring a, KnownNat n, 3 <= n, G.Vector v (SU.Vector n Word, a))   => MultiPoly v n a pattern Z = Multi.Z' --- | Evaluate at a given point.+-- | Evaluate the polynomial at a given point. -- -- >>> :set -XDataKinds -- >>> import Data.Vector.Generic.Sized (fromTuple) -- >>> eval (X^2 + Y^2 :: UMultiPoly 2 Int) (fromTuple (3, 4) :: Data.Vector.Sized.Vector 2 Int) -- 25+--+-- @since 0.5.0.0 eval   :: (Semiring a, G.Vector v (SU.Vector n Word, a), G.Vector u a)   => MultiPoly v n a@@ -113,12 +128,14 @@   -> a eval = Multi.eval' --- | Substitute another polynomials instead of variables.+-- | Substitute other polynomials instead of the variables. -- -- >>> :set -XDataKinds -- >>> import Data.Vector.Generic.Sized (fromTuple) -- >>> subst (X^2 + Y^2 + Z^2 :: UMultiPoly 3 Int) (fromTuple (X + 1, Y + 1, X + Y :: UMultiPoly 2 Int)) -- 2 * X^2 + 2 * X * Y + 2 * X + 2 * Y^2 + 2 * Y + 2+--+-- @since 0.5.0.0 subst   :: (Eq a, Semiring a, KnownNat m, G.Vector v (SU.Vector n Word, a), G.Vector w (SU.Vector m Word, a))   => MultiPoly v n a@@ -126,13 +143,15 @@   -> MultiPoly w m a subst = Multi.subst' --- | Take a derivative with respect to the /i/-th variable.+-- | Take the derivative of the polynomial with respect to the /i/-th variable. -- -- >>> :set -XDataKinds -- >>> deriv 0 (X^3 + 3 * Y) :: UMultiPoly 2 Int -- 3 * X^2 -- >>> deriv 1 (X^3 + 3 * Y) :: UMultiPoly 2 Int -- 3+--+-- @since 0.5.0.0 deriv   :: (Eq a, Semiring a, G.Vector v (SU.Vector n Word, a))   => Finite n@@ -140,15 +159,17 @@   -> MultiPoly v n a deriv = Multi.deriv' --- | Compute an indefinite integral of a polynomial--- by the /i/-th variable,--- setting constant term to zero.+-- | Compute an indefinite integral of the polynomial+-- with respect to the /i/-th variable,+-- setting the constant term to zero. -- -- >>> :set -XDataKinds -- >>> integral 0 (3 * X^2 + 2 * Y) :: UMultiPoly 2 Double -- 1.0 * X^3 + 2.0 * X * Y -- >>> integral 1 (3 * X^2 + 2 * Y) :: UMultiPoly 2 Double -- 3.0 * X^2 * Y + 1.0 * Y^2+--+-- @since 0.5.0.0 integral   :: (Field a, G.Vector v (SU.Vector n Word, a))   => Finite n

src/Data/Poly/Orthogonal.hs view

@@ -6,6 +6,7 @@ -- -- Classical orthogonal polynomials. --+-- @since 0.4.0.0  {-# LANGUAGE OverloadedLists     #-} {-# LANGUAGE RebindableSyntax    #-}@@ -34,6 +35,8 @@ -- -- >>> take 3 legendre :: [Data.Poly.VPoly Double] -- [1.0,1.0 * X + 0.0,1.5 * X^2 + 0.0 * X + (-0.5)]+--+-- @since 0.4.0.0 legendre :: (Eq a, Field a, Vector v a) => [Poly v a] legendre = map (`subst'` toPoly [1 `quot` 2, 1 `quot` 2]) legendreShifted   where@@ -44,6 +47,8 @@ -- -- >>> take 3 legendreShifted :: [Data.Poly.VPoly Integer] -- [1,2 * X + (-1),6 * X^2 + (-6) * X + 1]+--+-- @since 0.4.0.0 legendreShifted :: (Eq a, Euclidean a, Ring a, Vector v a) => [Poly v a] legendreShifted = xs   where@@ -51,12 +56,16 @@     rec n pm1 p = unscale' 0 (n + 1) (toPoly [-1 - 2 * n, 2 + 4 * n] * p - scale 0 n pm1)  -- | <https://en.wikipedia.org/wiki/Gegenbauer_polynomials Gegenbauer polynomials>.+--+-- @since 0.4.0.0 gegenbauer :: (Eq a, Field a, Vector v a) => a -> [Poly v a] gegenbauer g = jacobi a a   where     a = g - 1 `quot` 2  -- | <https://en.wikipedia.org/wiki/Jacobi_polynomials Jacobi polynomials>.+--+-- @since 0.4.0.0 jacobi :: (Eq a, Field a, Vector v a) => a -> a -> [Poly v a] jacobi a b = xs   where@@ -74,8 +83,10 @@ -- | <https://en.wikipedia.org/wiki/Chebyshev_polynomials Chebyshev polynomials> -- of the first kind. ----- >>> take 3 chebyshev1 :: [VPoly Integer]+-- >>> take 3 chebyshev1 :: [Data.Poly.VPoly Integer] -- [1,1 * X + 0,2 * X^2 + 0 * X + (-1)]+--+-- @since 0.4.0.0 chebyshev1 :: (Eq a, Ring a, Vector v a) => [Poly v a] chebyshev1 = xs   where@@ -84,8 +95,10 @@ -- | <https://en.wikipedia.org/wiki/Chebyshev_polynomials Chebyshev polynomials> -- of the second kind. ----- >>> take 3 chebyshev2 :: [VPoly Integer]+-- >>> take 3 chebyshev2 :: [Data.Poly.VPoly Integer] -- [1,2 * X + 0,4 * X^2 + 0 * X + (-1)]+--+-- @since 0.4.0.0 chebyshev2 :: (Eq a, Ring a, Vector v a) => [Poly v a] chebyshev2 = xs   where@@ -93,8 +106,10 @@  -- | Probabilists' <https://en.wikipedia.org/wiki/Hermite_polynomials Hermite polynomials>. ----- >>> take 3 hermiteProb :: [VPoly Integer]+-- >>> take 3 hermiteProb :: [Data.Poly.VPoly Integer] -- [1,1 * X + 0,1 * X^2 + 0 * X + (-1)]+--+-- @since 0.4.0.0 hermiteProb :: (Eq a, Ring a, Vector v a) => [Poly v a] hermiteProb = xs   where@@ -103,8 +118,10 @@  -- | Physicists' <https://en.wikipedia.org/wiki/Hermite_polynomials Hermite polynomials>. ----- >>> take 3 hermitePhys :: [VPoly Double]+-- >>> take 3 hermitePhys :: [Data.Poly.VPoly Double] -- [1.0,2.0 * X + 0.0,4.0 * X^2 + 0.0 * X + (-2.0)]+--+-- @since 0.4.0.0 hermitePhys :: (Eq a, Ring a, Vector v a) => [Poly v a] hermitePhys = xs   where@@ -113,12 +130,16 @@  -- | <https://en.wikipedia.org/wiki/Laguerre_polynomials Laguerre polynomials>. ----- >>> take 3 laguerre :: [VPoly Double]+-- >>> take 3 laguerre :: [Data.Poly.VPoly Double] -- [1.0,(-1.0) * X + 1.0,0.5 * X^2 + (-2.0) * X + 1.0]+--+-- @since 0.4.0.0 laguerre :: (Eq a, Field a, Vector v a) => [Poly v a] laguerre = laguerreGen 0  -- | <https://en.wikipedia.org/wiki/Laguerre_polynomials#Generalized_Laguerre_polynomials Generalized Laguerre polynomials>+--+-- @since 0.4.0.0 laguerreGen :: (Eq a, Field a, Vector v a) => a -> [Poly v a] laguerreGen a = xs   where

src/Data/Poly/Semiring.hs view

@@ -6,7 +6,9 @@ -- -- Dense polynomials and a 'Semiring'-based interface. --+-- @since 0.2.0.0 +{-# LANGUAGE CPP              #-} {-# LANGUAGE DataKinds        #-} {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE PatternSynonyms  #-}@@ -25,8 +27,11 @@   , subst   , deriv   , integral+  , timesRing+#ifdef SupportSparse   , denseToSparse   , sparseToDense+#endif   , dft   , inverseDft   , dftMult@@ -36,28 +41,35 @@ import Data.Euclidean (Field) import Data.Semiring (Semiring(..)) import qualified Data.Vector.Generic as G-import qualified Data.Vector.Unboxed.Sized as SU -import qualified Data.Poly.Internal.Convert as Convert-import Data.Poly.Internal.Dense (Poly(..), VPoly, UPoly, leading)+import Data.Poly.Internal.Dense (Poly(..), VPoly, UPoly, leading, timesRing) import qualified Data.Poly.Internal.Dense as Dense import Data.Poly.Internal.Dense.Field () import Data.Poly.Internal.Dense.DFT import Data.Poly.Internal.Dense.GcdDomain ()++#ifdef SupportSparse+import qualified Data.Vector.Unboxed.Sized as SU import qualified Data.Poly.Internal.Multi as Sparse+import qualified Data.Poly.Internal.Convert as Convert+#endif --- | Make 'Poly' from a vector of coefficients--- (first element corresponds to a constant term).+-- | Make a 'Poly' from a vector of coefficients+-- (first element corresponds to the constant term). -- -- >>> :set -XOverloadedLists -- >>> toPoly [1,2,3] :: VPoly Integer -- 3 * X^2 + 2 * X + 1 -- >>> toPoly [0,0,0] :: UPoly Int -- 0+--+-- @since 0.2.0.0 toPoly :: (Eq a, Semiring a, G.Vector v a) => v a -> Poly v a toPoly = Dense.toPoly'  -- | Create a monomial from a power and a coefficient.+--+-- @since 0.3.0.0 monomial :: (Eq a, Semiring a, G.Vector v a) => Word -> a -> Poly v a monomial = Dense.monomial' @@ -65,17 +77,25 @@ -- -- >>> scale 2 3 (X^2 + 1) :: UPoly Int -- 3 * X^4 + 0 * X^3 + 3 * X^2 + 0 * X + 0+--+-- @since 0.3.0.0 scale :: (Eq a, Semiring a, G.Vector v a) => Word -> a -> Poly v a -> Poly v a scale = Dense.scale' --- | Create an identity polynomial.+-- | The polynomial 'X'.+--+-- > X == monomial 1 one+--+-- @since 0.2.0.0 pattern X :: (Eq a, Semiring a, G.Vector v a) => Poly v a pattern X = Dense.X' --- | Evaluate at a given point.+-- | Evaluate the polynomial at a given point. -- -- >>> eval (X^2 + 1 :: UPoly Int) 3 -- 10+--+-- @since 0.2.0.0 eval :: (Semiring a, G.Vector v a) => Poly v a -> a -> a eval = Dense.eval' @@ -83,21 +103,27 @@ -- -- >>> subst (X^2 + 1 :: UPoly Int) (X + 1 :: UPoly Int) -- 1 * X^2 + 2 * X + 2+--+-- @since 0.3.3.0 subst :: (Eq a, Semiring a, G.Vector v a, G.Vector w a) => Poly v a -> Poly w a -> Poly w a subst = Dense.subst' --- | Take a derivative.+-- | Take the derivative of the polynomial. -- -- >>> deriv (X^3 + 3 * X) :: UPoly Int -- 3 * X^2 + 0 * X + 3+--+-- @since 0.2.0.0 deriv :: (Eq a, Semiring a, G.Vector v a) => Poly v a -> Poly v a deriv = Dense.deriv' --- | Compute an indefinite integral of a polynomial,--- setting constant term to zero.+-- | Compute an indefinite integral of the polynomial,+-- setting the constant term to zero. -- -- >>> integral (3 * X^2 + 3) :: UPoly Double -- 1.0 * X^3 + 0.0 * X^2 + 3.0 * X + 0.0+--+-- @since 0.3.2.0 integral :: (Eq a, Field a, G.Vector v a) => Poly v a -> Poly v a integral = Dense.integral' @@ -105,6 +131,8 @@ -- <https://en.wikipedia.org/wiki/Fast_Fourier_transform discrete Fourier transform>. -- It could be faster than '(*)' for large polynomials -- if multiplication of coefficients is particularly expensive.+--+-- @since 0.5.0.0 dftMult   :: (Eq a, Field a, G.Vector v a)   => (Int -> a) -- ^ mapping from \( N = 2^n \) to a primitive root \( \sqrt[N]{1} \)@@ -125,19 +153,26 @@     xs' = padTo zl xs     ys' = padTo zl ys     primRoot = getPrimRoot zl+{-# INLINABLE dftMult #-} +#ifdef SupportSparse -- | Convert from dense to sparse polynomials. -- -- >>> :set -XFlexibleContexts--- >>> denseToSparse (1 `plus` Data.Poly.X^2) :: Data.Poly.Sparse.UPoly Int+-- >>> denseToSparse (1 `Data.Semiring.plus` Data.Poly.X^2) :: Data.Poly.Sparse.UPoly Int -- 1 * X^2 + 1+--+-- @since 0.5.0.0 denseToSparse :: (Eq a, Semiring a, G.Vector v a, G.Vector v (SU.Vector 1 Word, a)) => Dense.Poly v a -> Sparse.Poly v a denseToSparse = Convert.denseToSparse'  -- | Convert from sparse to dense polynomials. -- -- >>> :set -XFlexibleContexts--- >>> sparseToDense (1 `plus` Data.Poly.Sparse.X^2) :: Data.Poly.UPoly Int+-- >>> sparseToDense (1 `Data.Semiring.plus` Data.Poly.Sparse.X^2) :: Data.Poly.UPoly Int -- 1 * X^2 + 0 * X + 1+--+-- @since 0.5.0.0 sparseToDense :: (Semiring a, G.Vector v a, G.Vector v (SU.Vector 1 Word, a)) => Sparse.Poly v a -> Dense.Poly v a sparseToDense = Convert.sparseToDense'+#endif

src/Data/Poly/Sparse.hs view

@@ -4,8 +4,10 @@ -- Licence:     BSD3 -- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com> ----- Sparse polynomials with 'Num' instance.+-- Sparse polynomials with a 'Num' instance. --+-- @since 0.3.0.0+--  {-# LANGUAGE DataKinds        #-} {-# LANGUAGE FlexibleContexts #-}@@ -41,31 +43,39 @@ import Data.Poly.Internal.Multi.Field (quotRemFractional) import Data.Poly.Internal.Multi.GcdDomain () --- | Make 'Poly' from a list of (power, coefficient) pairs.+-- | Make a 'Poly' from a list of (power, coefficient) pairs. -- -- >>> :set -XOverloadedLists -- >>> toPoly [(0,1),(1,2),(2,3)] :: VPoly Integer -- 3 * X^2 + 2 * X + 1 -- >>> toPoly [(0,0),(1,0),(2,0)] :: UPoly Int -- 0+--+-- @since 0.3.0.0 toPoly   :: (Eq a, Num a, G.Vector v (Word, a), G.Vector v (SU.Vector 1 Word, a))   => v (Word, a)   -> Poly v a toPoly = Multi.toMultiPoly . G.map (first SU.singleton)+{-# INLINABLE toPoly #-}  -- | Create a monomial from a power and a coefficient.+--+-- @since 0.3.0.0 monomial   :: (Eq a, Num a, G.Vector v (SU.Vector 1 Word, a))   => Word   -> a   -> Poly v a monomial = Multi.monomial . SU.singleton+{-# INLINABLE monomial #-}  -- | Multiply a polynomial by a monomial, expressed as a power and a coefficient. -- -- >>> scale 2 3 (X^2 + 1) :: UPoly Int -- 3 * X^4 + 3 * X^2+--+-- @since 0.3.0.0 scale   :: (Eq a, Num a, G.Vector v (SU.Vector 1 Word, a))   => Word@@ -73,52 +83,69 @@   -> Poly v a   -> Poly v a scale = Multi.scale . SU.singleton+{-# INLINABLE scale #-} --- | Create an identity polynomial.+-- | The polynomial 'X'.+--+-- > X == monomial 1 1+--+-- @since 0.3.0.0 pattern X   :: (Eq a, Num a, G.Vector v (SU.Vector 1 Word, a))   => Poly v a pattern X = Multi.X --- | Evaluate at a given point.+-- | Evaluate the polynomial at a given point. -- -- >>> eval (X^2 + 1 :: UPoly Int) 3 -- 10+--+-- @since 0.3.0.0 eval   :: (Num a, G.Vector v (SU.Vector 1 Word, a))   => Poly v a   -> a   -> a eval p = Multi.eval p . SV.singleton+{-# INLINABLE eval #-}  -- | Substitute another polynomial instead of 'X'. -- -- >>> subst (X^2 + 1 :: UPoly Int) (X + 1 :: UPoly Int) -- 1 * X^2 + 2 * X + 2+--+-- @since 0.3.3.0 subst   :: (Eq a, Num a, G.Vector v (SU.Vector 1 Word, a), G.Vector w (SU.Vector 1 Word, a))   => Poly v a   -> Poly w a   -> Poly w a subst p = Multi.subst p . SV.singleton+{-# INLINABLE subst #-} --- | Take a derivative.+-- | Take the derivative of the polynomial. -- -- >>> deriv (X^3 + 3 * X) :: UPoly Int -- 3 * X^2 + 3+--+-- @since 0.3.0.0 deriv   :: (Eq a, Num a, G.Vector v (SU.Vector 1 Word, a))   => Poly v a   -> Poly v a deriv = Multi.deriv 0+{-# INLINABLE deriv #-} --- | Compute an indefinite integral of a polynomial,--- setting constant term to zero.+-- | Compute an indefinite integral of the polynomial,+-- setting the constant term to zero. -- -- >>> integral (3 * X^2 + 3) :: UPoly Double -- 1.0 * X^3 + 3.0 * X+--+-- @since 0.3.0.0 integral   :: (Fractional a, G.Vector v (SU.Vector 1 Word, a))   => Poly v a   -> Poly v a integral = Multi.integral 0+{-# INLINABLE integral #-}

src/Data/Poly/Sparse/Laurent.hs view

@@ -7,6 +7,7 @@ -- Sparse -- <https://en.wikipedia.org/wiki/Laurent_polynomial Laurent polynomials>. --+-- @since 0.4.0.0  {-# LANGUAGE DataKinds                  #-} {-# LANGUAGE FlexibleContexts           #-}@@ -39,17 +40,22 @@ import Data.Poly.Internal.Multi (Poly)  -- | Create a monomial from a power and a coefficient.+--+-- @since 0.4.0.0 monomial   :: (Eq a, Semiring a, G.Vector v (SU.Vector 1 Word, a))   => Int   -> a   -> Laurent v a monomial = Multi.monomial . SU.singleton+{-# INLINABLE monomial #-}  -- | Multiply a polynomial by a monomial, expressed as a power and a coefficient. -- -- >>> scale 2 3 (X^-2 + 1) :: ULaurent Int -- 3 * X^2 + 3+--+-- @since 0.4.0.0 scale   :: (Eq a, Semiring a, G.Vector v (SU.Vector 1 Word, a))   => Int@@ -57,18 +63,30 @@   -> Laurent v a   -> Laurent v a scale = Multi.scale . SU.singleton+{-# INLINABLE scale #-} --- | Create an identity polynomial.+-- | The polynomial 'X'.+--+-- > X == monomial 1 one+--+-- @since 0.4.0.0 pattern X   :: (Eq a, Semiring a, G.Vector v (SU.Vector 1 Word, a))   => Laurent v a pattern X = Multi.X --- | This operator can be applied only to monomials with unit coefficients,--- but is instrumental to express Laurent polynomials in mathematical fashion:+-- | Used to construct monomials with negative powers. --+-- This operator can be applied only to monomials with unit coefficients,+-- but is instrumental to express Laurent polynomials+-- in a mathematical fashion:+--+-- >>> X^-3 :: ULaurent Int+-- 1 * X^-3 -- >>> X + 2 + 3 * (X^2)^-1 :: ULaurent Int -- 1 * X + 2 + 3 * X^-2+--+-- @since 0.4.0.0 (^-)   :: (Eq a, Semiring a, G.Vector v (SU.Vector 1 Word, a))   => Laurent v a@@ -76,35 +94,44 @@   -> Laurent v a (^-) = (Multi.^-) --- | Evaluate at a given point.+-- | Evaluate the polynomial at a given point. -- -- >>> eval (X^-2 + 1 :: ULaurent Double) 2 -- 1.25+--+-- @since 0.4.0.0 eval   :: (Field a, G.Vector v (SU.Vector 1 Word, a))   => Laurent v a   -> a   -> a eval p = Multi.eval p . SV.singleton+{-# INLINABLE eval #-}  -- | Substitute another polynomial instead of 'X'. -- -- >>> import Data.Poly.Sparse (UPoly) -- >>> subst (Data.Poly.Sparse.X^2 + 1 :: UPoly Int) (X^-1 + 1 :: ULaurent Int) -- 2 + 2 * X^-1 + 1 * X^-2+--+-- @since 0.4.0.0 subst   :: (Eq a, Semiring a, G.Vector v (SU.Vector 1 Word, a), G.Vector w (SU.Vector 1 Word, a))   => Poly v a   -> Laurent w a   -> Laurent w a subst p = Multi.subst p . SV.singleton+{-# INLINABLE subst #-} --- | Take a derivative.+-- | Take the derivative of the polynomial. -- -- >>> deriv (X^-3 + 3 * X) :: ULaurent Int -- 3 + (-3) * X^-4+--+-- @since 0.4.0.0 deriv   :: (Eq a, Ring a, G.Vector v (SU.Vector 1 Word, a))   => Laurent v a   -> Laurent v a deriv = Multi.deriv 0+{-# INLINABLE deriv #-}

src/Data/Poly/Sparse/Semiring.hs view

@@ -4,8 +4,10 @@ -- Licence:     BSD3 -- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com> ----- Sparse polynomials with 'Semiring' instance.+-- Sparse polynomials with a 'Semiring' instance. --+-- @since 0.3.0.0+--  {-# LANGUAGE DataKinds        #-} {-# LANGUAGE FlexibleContexts #-}@@ -43,31 +45,39 @@ import Data.Poly.Internal.Multi.Field () import Data.Poly.Internal.Multi.GcdDomain () --- | Make 'Poly' from a list of (power, coefficient) pairs.+-- | Make a 'Poly' from a list of (power, coefficient) pairs. -- -- >>> :set -XOverloadedLists -- >>> toPoly [(0,1),(1,2),(2,3)] :: VPoly Integer -- 3 * X^2 + 2 * X + 1 -- >>> toPoly [(0,0),(1,0),(2,0)] :: UPoly Int -- 0+--+-- @since 0.3.0.0 toPoly   :: (Eq a, Semiring a, G.Vector v (Word, a), G.Vector v (SU.Vector 1 Word, a))   => v (Word, a)   -> Poly v a toPoly = Multi.toMultiPoly' . G.map (first SU.singleton)+{-# INLINABLE toPoly #-}  -- | Create a monomial from a power and a coefficient.+--+-- @since 0.3.0.0 monomial   :: (Eq a, Semiring a, G.Vector v (SU.Vector 1 Word, a))   => Word   -> a   -> Poly v a monomial = Multi.monomial' . SU.singleton+{-# INLINABLE monomial #-}  -- | Multiply a polynomial by a monomial, expressed as a power and a coefficient. -- -- >>> scale 2 3 (X^2 + 1) :: UPoly Int -- 3 * X^4 + 3 * X^2+--+-- @since 0.3.0.0 scale   :: (Eq a, Semiring a, G.Vector v (SU.Vector 1 Word, a))   => Word@@ -75,74 +85,97 @@   -> Poly v a   -> Poly v a scale = Multi.scale' . SU.singleton+{-# INLINABLE scale #-} --- | Create an identity polynomial.+-- | The polynomial 'X'.+--+-- > X == monomial 1 one+--+-- @since 0.3.0.0 pattern X   :: (Eq a, Semiring a, G.Vector v (SU.Vector 1 Word, a))   => Poly v a pattern X = Multi.X' --- | Evaluate at a given point.+-- | Evaluate the polynomial at a given point. -- -- >>> eval (X^2 + 1 :: UPoly Int) 3 -- 10+--+-- @since 0.3.0.0 eval   :: (Semiring a, G.Vector v (SU.Vector 1 Word, a))   => Poly v a   -> a   -> a eval p = Multi.eval' p . SV.singleton+{-# INLINABLE eval #-}  -- | Substitute another polynomial instead of 'X'. -- -- >>> subst (X^2 + 1 :: UPoly Int) (X + 1 :: UPoly Int) -- 1 * X^2 + 2 * X + 2+--+-- @since 0.3.3.0 subst   :: (Eq a, Semiring a, G.Vector v (SU.Vector 1 Word, a), G.Vector w (SU.Vector 1 Word, a))   => Poly v a   -> Poly w a   -> Poly w a subst p = Multi.subst' p . SV.singleton+{-# INLINABLE subst #-} --- | Take a derivative.+-- | Take the derivative of the polynomial. -- -- >>> deriv (X^3 + 3 * X) :: UPoly Int -- 3 * X^2 + 3+--+-- @since 0.3.0.0 deriv   :: (Eq a, Semiring a, G.Vector v (SU.Vector 1 Word, a))   => Poly v a   -> Poly v a deriv = Multi.deriv' 0+{-# INLINABLE deriv #-} --- | Compute an indefinite integral of a polynomial,--- setting constant term to zero.+-- | Compute an indefinite integral of the polynomial,+-- setting the constant term to zero. -- -- >>> integral (3 * X^2 + 3) :: UPoly Double -- 1.0 * X^3 + 3.0 * X+--+-- @since 0.3.2.0 integral   :: (Field a, G.Vector v (SU.Vector 1 Word, a))   => Poly v a   -> Poly v a integral = Multi.integral' 0+{-# INLINABLE integral #-}  -- | Convert from dense to sparse polynomials. -- -- >>> :set -XFlexibleContexts--- >>> denseToSparse (1 `plus` Data.Poly.X^2) :: Data.Poly.Sparse.UPoly Int+-- >>> denseToSparse (1 `Data.Semiring.plus` Data.Poly.X^2) :: Data.Poly.Sparse.UPoly Int -- 1 * X^2 + 1+--+-- @since 0.5.0.0 denseToSparse   :: (Eq a, Semiring a, G.Vector v a, G.Vector v (SU.Vector 1 Word, a))   => Dense.Poly v a   -> Multi.Poly v a denseToSparse = Convert.denseToSparse'+{-# INLINABLE denseToSparse #-}  -- | Convert from sparse to dense polynomials. -- -- >>> :set -XFlexibleContexts--- >>> sparseToDense (1 `plus` Data.Poly.Sparse.X^2) :: Data.Poly.UPoly Int+-- >>> sparseToDense (1 `Data.Semiring.plus` Data.Poly.Sparse.X^2) :: Data.Poly.UPoly Int -- 1 * X^2 + 0 * X + 1+--+-- @since 0.5.0.0 sparseToDense   :: (Semiring a, G.Vector v a, G.Vector v (SU.Vector 1 Word, a))   => Multi.Poly v a   -> Dense.Poly v a sparseToDense = Convert.sparseToDense'+{-# INLINABLE sparseToDense #-}

test/Dense.hs view

@@ -1,3 +1,4 @@+{-# LANGUAGE CPP                        #-} {-# LANGUAGE DataKinds                  #-} {-# LANGUAGE FlexibleContexts           #-} {-# LANGUAGE FlexibleInstances          #-}@@ -43,20 +44,28 @@   $ semiringTests ++ ringTests ++ numTests ++ euclideanTests ++ gcdDomainTests ++ isListTests ++ showTests  semiringTests :: [TestTree]+#ifdef MIN_VERSION_quickcheck_classes semiringTests =   [ mySemiringLaws (Proxy :: Proxy (UPoly ()))   , mySemiringLaws (Proxy :: Proxy (UPoly Int8))   , mySemiringLaws (Proxy :: Proxy (VPoly Integer))   , mySemiringLaws (Proxy :: Proxy (UPoly (Quaternion Int)))   ]+#else+semiringTests = []+#endif  ringTests :: [TestTree]+#ifdef MIN_VERSION_quickcheck_classes ringTests =   [ myRingLaws (Proxy :: Proxy (UPoly ()))   , myRingLaws (Proxy :: Proxy (UPoly Int8))   , myRingLaws (Proxy :: Proxy (VPoly Integer))   , myRingLaws (Proxy :: Proxy (UPoly (Quaternion Int)))   ]+#else+ringTests = []+#endif  numTests :: [TestTree] numTests =@@ -66,17 +75,25 @@   ]  gcdDomainTests :: [TestTree]+#ifdef MIN_VERSION_quickcheck_classes gcdDomainTests =   [ myGcdDomainLaws (Proxy :: Proxy (ShortPoly (VPoly Integer)))   , myGcdDomainLaws (Proxy :: Proxy (ShortPoly (UPoly (Mod 3))))   , myGcdDomainLaws (Proxy :: Proxy (ShortPoly (VPoly Rational)))   ]+#else+gcdDomainTests = []+#endif  euclideanTests :: [TestTree]+#ifdef MIN_VERSION_quickcheck_classes euclideanTests =   [ myEuclideanLaws (Proxy :: Proxy (ShortPoly (UPoly (Mod 3))))   , myEuclideanLaws (Proxy :: Proxy (ShortPoly (VPoly Rational)))   ]+#else+euclideanTests = []+#endif  isListTests :: [TestTree] isListTests =@@ -271,12 +288,14 @@  conversionTests :: TestTree conversionTests = testGroup "conversions"-  [ testProperty "sparseToDense . denseToSparse = id" $-    \(xs :: UPoly Int8) -> xs === sparseToDense (denseToSparse xs)-  , testProperty "sparseToDense' . denseToSparse' = id" $-    \(xs :: UPoly Int8) -> xs === S.sparseToDense (S.denseToSparse xs)-  , testProperty "toPoly . unPoly = id" $+  [ testProperty "toPoly . unPoly = id" $     \(xs :: UPoly Int8) -> xs === toPoly (unPoly xs)   , testProperty "S.toPoly . S.unPoly = id" $     \(xs :: UPoly Int8) -> xs === S.toPoly (S.unPoly xs)+#ifdef SupportSparse+  , testProperty "sparseToDense . denseToSparse = id" $+    \(xs :: UPoly Int8) -> xs === sparseToDense (denseToSparse xs)+  , testProperty "sparseToDense' . denseToSparse' = id" $+    \(xs :: UPoly Int8) -> xs === S.sparseToDense (S.denseToSparse xs)+#endif   ]

test/DenseLaurent.hs view

@@ -1,3 +1,4 @@+{-# LANGUAGE CPP                        #-} {-# LANGUAGE FlexibleContexts           #-} {-# LANGUAGE FlexibleInstances          #-} {-# LANGUAGE ScopedTypeVariables        #-}@@ -37,20 +38,28 @@   $ semiringTests ++ ringTests ++ numTests ++ gcdDomainTests ++ showTests  semiringTests :: [TestTree]+#ifdef MIN_VERSION_quickcheck_classes semiringTests =   [ mySemiringLaws (Proxy :: Proxy (ULaurent ()))   , mySemiringLaws (Proxy :: Proxy (ULaurent Int8))   , mySemiringLaws (Proxy :: Proxy (VLaurent Integer))   , mySemiringLaws (Proxy :: Proxy (ULaurent (Quaternion Int)))   ]+#else+semiringTests = []+#endif  ringTests :: [TestTree]+#ifdef MIN_VERSION_quickcheck_classes ringTests =   [ myRingLaws (Proxy :: Proxy (ULaurent ()))   , myRingLaws (Proxy :: Proxy (ULaurent Int8))   , myRingLaws (Proxy :: Proxy (VLaurent Integer))   , myRingLaws (Proxy :: Proxy (ULaurent (Quaternion Int)))   ]+#else+ringTests = []+#endif  numTests :: [TestTree] numTests =@@ -60,10 +69,14 @@   ]  gcdDomainTests :: [TestTree]+#ifdef MIN_VERSION_quickcheck_classes gcdDomainTests =   [ myGcdDomainLaws (Proxy :: Proxy (ShortPoly (VLaurent Integer)))   , myGcdDomainLaws (Proxy :: Proxy (ShortPoly (VLaurent Rational)))   ]+#else+gcdDomainTests = []+#endif  showTests :: [TestTree] showTests =

test/Main.hs view

@@ -1,3 +1,5 @@+{-# LANGUAGE CPP #-}+ module Main where  import Test.Tasty@@ -5,20 +7,24 @@ import qualified Dense import qualified DenseLaurent import qualified DFT+import qualified Orthogonal+#ifdef SupportSparse import qualified Multi import qualified MultiLaurent-import qualified Orthogonal import qualified Sparse import qualified SparseLaurent+#endif  main :: IO () main = defaultMain $ testGroup "All"     [ Dense.testSuite     , DenseLaurent.testSuite     , DFT.testSuite+    , Orthogonal.testSuite+#ifdef SupportSparse     , Sparse.testSuite     , SparseLaurent.testSuite     , Multi.testSuite     , MultiLaurent.testSuite-    , Orthogonal.testSuite+#endif     ]

test/Multi.hs view

@@ -1,3 +1,4 @@+{-# LANGUAGE CPP                        #-} {-# LANGUAGE DataKinds                  #-} {-# LANGUAGE FlexibleContexts           #-} {-# LANGUAGE FlexibleInstances          #-}@@ -49,6 +50,7 @@   $ semiringTests ++ ringTests ++ numTests ++ gcdDomainTests ++ isListTests ++ showTests  semiringTests :: [TestTree]+#ifdef MIN_VERSION_quickcheck_classes semiringTests =   [ mySemiringLaws (Proxy :: Proxy (UMultiPoly 3 ()))   , mySemiringLaws (Proxy :: Proxy (ShortPoly (UMultiPoly 2 Int8)))@@ -56,14 +58,21 @@   , tenTimesLess   $ mySemiringLaws (Proxy :: Proxy (ShortPoly (UMultiPoly 2 (Quaternion Int))))   ]+#else+semiringTests = []+#endif  ringTests :: [TestTree]+#ifdef MIN_VERSION_quickcheck_classes ringTests =   [ myRingLaws (Proxy :: Proxy (UMultiPoly 3 ()))   , myRingLaws (Proxy :: Proxy (UMultiPoly 3 Int8))   , myRingLaws (Proxy :: Proxy (VMultiPoly 3 Integer))   , myRingLaws (Proxy :: Proxy (UMultiPoly 3 (Quaternion Int)))   ]+#else+ringTests = []+#endif  numTests :: [TestTree] numTests =@@ -74,6 +83,7 @@   ]  gcdDomainTests :: [TestTree]+#ifdef MIN_VERSION_quickcheck_classes gcdDomainTests =   [ myGcdDomainLaws (Proxy :: Proxy (ShortPoly (VMultiPoly 3 Integer)))   , tenTimesLess@@ -81,6 +91,9 @@   , tenTimesLess   $ myGcdDomainLaws (Proxy :: Proxy (ShortPoly (VMultiPoly 3 Rational)))   ]+#else+gcdDomainTests = []+#endif  isListTests :: [TestTree] isListTests =

test/MultiLaurent.hs view

@@ -1,3 +1,4 @@+{-# LANGUAGE CPP                        #-} {-# LANGUAGE DataKinds                  #-} {-# LANGUAGE FlexibleContexts           #-} {-# LANGUAGE FlexibleInstances          #-}@@ -43,6 +44,7 @@   $ semiringTests ++ ringTests ++ numTests ++ gcdDomainTests ++ isListTests ++ showTests  semiringTests :: [TestTree]+#ifdef MIN_VERSION_quickcheck_classes semiringTests =   [ mySemiringLaws (Proxy :: Proxy (UMultiLaurent 3 ()))   , mySemiringLaws (Proxy :: Proxy (ShortPoly (UMultiLaurent 2 Int8)))@@ -50,14 +52,21 @@   , tenTimesLess   $ mySemiringLaws (Proxy :: Proxy (ShortPoly (UMultiLaurent 2 (Quaternion Int))))   ]+#else+semiringTests = []+#endif  ringTests :: [TestTree]+#ifdef MIN_VERSION_quickcheck_classes ringTests =   [ myRingLaws (Proxy :: Proxy (UMultiLaurent 3 ()))   , myRingLaws (Proxy :: Proxy (UMultiLaurent 3 Int8))   , myRingLaws (Proxy :: Proxy (VMultiLaurent 3 Integer))   , myRingLaws (Proxy :: Proxy (UMultiLaurent 3 (Quaternion Int)))   ]+#else+ringTests = []+#endif  numTests :: [TestTree] numTests =@@ -68,11 +77,15 @@   ]  gcdDomainTests :: [TestTree]+#ifdef MIN_VERSION_quickcheck_classes gcdDomainTests =   [ myGcdDomainLaws (Proxy :: Proxy (ShortPoly (VMultiLaurent 3 Integer)))   , tenTimesLess   $ myGcdDomainLaws (Proxy :: Proxy (ShortPoly (VMultiLaurent 3 Rational)))   ]+#else+gcdDomainTests = []+#endif  isListTests :: [TestTree] isListTests =

test/Sparse.hs view

@@ -1,3 +1,4 @@+{-# LANGUAGE CPP                        #-} {-# LANGUAGE DataKinds                  #-} {-# LANGUAGE FlexibleContexts           #-} {-# LANGUAGE FlexibleInstances          #-}@@ -47,6 +48,7 @@   $ semiringTests ++ ringTests ++ numTests ++ euclideanTests ++ gcdDomainTests ++ isListTests ++ showTests  semiringTests :: [TestTree]+#ifdef MIN_VERSION_quickcheck_classes semiringTests =   [ mySemiringLaws (Proxy :: Proxy (UPoly ()))   , mySemiringLaws (Proxy :: Proxy (UPoly Int8))@@ -54,14 +56,21 @@   , tenTimesLess   $ mySemiringLaws (Proxy :: Proxy (UPoly (Quaternion Int)))   ]+#else+semiringTests = []+#endif  ringTests :: [TestTree]+#ifdef MIN_VERSION_quickcheck_classes ringTests =   [ myRingLaws (Proxy :: Proxy (UPoly ()))   , myRingLaws (Proxy :: Proxy (UPoly Int8))   , myRingLaws (Proxy :: Proxy (VPoly Integer))   , myRingLaws (Proxy :: Proxy (UPoly (Quaternion Int)))   ]+#else+ringTests = []+#endif  numTests :: [TestTree] numTests =@@ -72,6 +81,7 @@   ]  gcdDomainTests :: [TestTree]+#ifdef MIN_VERSION_quickcheck_classes gcdDomainTests =   [ myGcdDomainLaws (Proxy :: Proxy (ShortPoly (VPoly Integer)))   , tenTimesLess@@ -79,12 +89,19 @@   , tenTimesLess   $ myGcdDomainLaws (Proxy :: Proxy (ShortPoly (VPoly Rational)))   ]+#else+gcdDomainTests = []+#endif  euclideanTests :: [TestTree]+#ifdef MIN_VERSION_quickcheck_classes euclideanTests =   [ myEuclideanLaws (Proxy :: Proxy (ShortPoly (UPoly (Mod 3))))   , myEuclideanLaws (Proxy :: Proxy (ShortPoly (VPoly Rational)))   ]+#else+euclideanTests = []+#endif  isListTests :: [TestTree] isListTests =

test/SparseLaurent.hs view

@@ -1,3 +1,4 @@+{-# LANGUAGE CPP                        #-} {-# LANGUAGE DataKinds                  #-} {-# LANGUAGE FlexibleContexts           #-} {-# LANGUAGE FlexibleInstances          #-}@@ -40,6 +41,7 @@   $ semiringTests ++ ringTests ++ numTests ++ gcdDomainTests ++ isListTests ++ showTests  semiringTests :: [TestTree]+#ifdef MIN_VERSION_quickcheck_classes semiringTests =   [ mySemiringLaws (Proxy :: Proxy (ULaurent ()))   , mySemiringLaws (Proxy :: Proxy (ULaurent Int8))@@ -47,14 +49,21 @@   , tenTimesLess   $ mySemiringLaws (Proxy :: Proxy (ULaurent (Quaternion Int)))   ]+#else+semiringTests = []+#endif  ringTests :: [TestTree]+#ifdef MIN_VERSION_quickcheck_classes ringTests =   [ myRingLaws (Proxy :: Proxy (ULaurent ()))   , myRingLaws (Proxy :: Proxy (ULaurent Int8))   , myRingLaws (Proxy :: Proxy (VLaurent Integer))   , myRingLaws (Proxy :: Proxy (ULaurent (Quaternion Int)))   ]+#else+ringTests = []+#endif  numTests :: [TestTree] numTests =@@ -65,11 +74,15 @@   ]  gcdDomainTests :: [TestTree]+#ifdef MIN_VERSION_quickcheck_classes gcdDomainTests =   [ myGcdDomainLaws (Proxy :: Proxy (ShortPoly (VLaurent Integer)))   , tenTimesLess   $ myGcdDomainLaws (Proxy :: Proxy (ShortPoly (VLaurent Rational)))   ]+#else+gcdDomainTests = []+#endif  isListTests :: [TestTree] isListTests =

test/TestUtils.hs view

@@ -1,3 +1,4 @@+{-# LANGUAGE CPP                        #-} {-# LANGUAGE DataKinds                  #-} {-# LANGUAGE FlexibleContexts           #-} {-# LANGUAGE FlexibleInstances          #-}@@ -10,35 +11,44 @@ module TestUtils   ( ShortPoly(..)   , tenTimesLess+  , myNumLaws+#ifdef MIN_VERSION_quickcheck_classes   , mySemiringLaws   , myRingLaws-  , myNumLaws   , myGcdDomainLaws   , myEuclideanLaws+#endif   , myIsListLaws   , myShowLaws   ) where  import Prelude hiding (lcm, rem)-import Control.Arrow import Data.Euclidean-import Data.Finite import Data.Mod.Word import Data.Proxy import Data.Semiring (Semiring(..), Ring) import qualified Data.Vector.Generic as G-import qualified Data.Vector.Generic.Sized as SG-import qualified Data.Vector.Unboxed.Sized as SU import GHC.Exts import GHC.TypeNats (KnownNat)-import Test.QuickCheck.Classes+import Test.QuickCheck.Classes.Base import Test.Tasty import Test.Tasty.QuickCheck +#ifdef MIN_VERSION_quickcheck_classes+import Test.QuickCheck.Classes+#endif+ import qualified Data.Poly.Semiring as Dense import qualified Data.Poly.Laurent as DenseLaurent++#ifdef SupportSparse+import Control.Arrow+import Data.Finite+import qualified Data.Vector.Generic.Sized as SG+import qualified Data.Vector.Unboxed.Sized as SU import Data.Poly.Multi.Semiring import qualified Data.Poly.Multi.Laurent as MultiLaurent+#endif  newtype ShortPoly a = ShortPoly { unShortPoly :: a }   deriving (Eq, Show, Semiring, GcdDomain, Euclidean, Num)@@ -47,13 +57,6 @@   arbitrary = oneof [arbitraryBoundedEnum, fromInteger <$> arbitrary]   shrink = map fromInteger . shrink . toInteger . unMod -instance KnownNat n => Arbitrary (Finite n) where-  arbitrary = elements finites--instance (Arbitrary a, KnownNat n, G.Vector v a) => Arbitrary (SG.Vector v n a) where-  arbitrary = SG.replicateM arbitrary-  shrink vs = [ vs SG.// [(i, x)] | i <- finites, x <- shrink (SG.index vs i) ]- instance (Eq a, Semiring a, Arbitrary a, G.Vector v a) => Arbitrary (Dense.Poly v a) where   arbitrary = Dense.toPoly . G.fromList <$> arbitrary   shrink = fmap (Dense.toPoly . G.fromList) . shrink . G.toList . Dense.unPoly@@ -70,6 +73,15 @@   arbitrary = (ShortPoly .) . DenseLaurent.toLaurent <$> ((`rem` 10) <$> arbitrary) <*> (unShortPoly <$> arbitrary)   shrink = fmap (ShortPoly . uncurry DenseLaurent.toLaurent . fmap unShortPoly) . shrink . fmap ShortPoly . DenseLaurent.unLaurent . unShortPoly +#ifdef SupportSparse++instance KnownNat n => Arbitrary (Finite n) where+  arbitrary = elements finites++instance (Arbitrary a, KnownNat n, G.Vector v a) => Arbitrary (SG.Vector v n a) where+  arbitrary = SG.replicateM arbitrary+  shrink vs = [ vs SG.// [(i, x)] | i <- finites, x <- shrink (SG.index vs i) ]+ instance (Eq a, Semiring a, Arbitrary a, KnownNat n, G.Vector v (SU.Vector n Word, a)) => Arbitrary (MultiPoly v n a) where   arbitrary = toMultiPoly . G.fromList <$> arbitrary   shrink = fmap (toMultiPoly . G.fromList) . shrink . G.toList . unMultiPoly@@ -86,16 +98,18 @@   arbitrary = (ShortPoly .) . MultiLaurent.toMultiLaurent <$> (SU.map (`rem` 10) <$> arbitrary) <*> (unShortPoly <$> arbitrary)   shrink = fmap (ShortPoly . uncurry MultiLaurent.toMultiLaurent . fmap unShortPoly) . shrink . fmap ShortPoly . MultiLaurent.unMultiLaurent . unShortPoly +#endif+ -------------------------------------------------------------------------------  tenTimesLess :: TestTree -> TestTree tenTimesLess = adjustOption $   \(QuickCheckTests n) -> QuickCheckTests (max 100 (n `div` 10)) -mySemiringLaws :: (Eq a, Semiring a, Arbitrary a, Show a) => Proxy a -> TestTree-mySemiringLaws proxy = testGroup tpclss $ map tune props+myNumLaws :: (Eq a, Num a, Arbitrary a, Show a) => Proxy a -> TestTree+myNumLaws proxy = testGroup tpclss $ map tune props   where-    Laws tpclss props = semiringLaws proxy+    Laws tpclss props = numLaws proxy      tune pair = case fst pair of       "Multiplicative Associativity" ->@@ -104,19 +118,18 @@         tenTimesLess test       "Multiplication Right Distributes Over Addition" ->         tenTimesLess test+      "Subtraction" ->+        tenTimesLess test       _ -> test       where         test = uncurry testProperty pair -myRingLaws :: (Eq a, Ring a, Arbitrary a, Show a) => Proxy a -> TestTree-myRingLaws proxy = testGroup tpclss $ map (uncurry testProperty) props-  where-    Laws tpclss props = ringLaws proxy+#ifdef MIN_VERSION_quickcheck_classes -myNumLaws :: (Eq a, Num a, Arbitrary a, Show a) => Proxy a -> TestTree-myNumLaws proxy = testGroup tpclss $ map tune props+mySemiringLaws :: (Eq a, Semiring a, Arbitrary a, Show a) => Proxy a -> TestTree+mySemiringLaws proxy = testGroup tpclss $ map tune props   where-    Laws tpclss props = numLaws proxy+    Laws tpclss props = semiringLaws proxy      tune pair = case fst pair of       "Multiplicative Associativity" ->@@ -125,12 +138,15 @@         tenTimesLess test       "Multiplication Right Distributes Over Addition" ->         tenTimesLess test-      "Subtraction" ->-        tenTimesLess test       _ -> test       where         test = uncurry testProperty pair +myRingLaws :: (Eq a, Ring a, Arbitrary a, Show a) => Proxy a -> TestTree+myRingLaws proxy = testGroup tpclss $ map (uncurry testProperty) props+  where+    Laws tpclss props = ringLaws proxy+ myGcdDomainLaws :: forall a. (Eq a, GcdDomain a, Arbitrary a, Show a) => Proxy a -> TestTree myGcdDomainLaws proxy = testGroup tpclss $ map tune $ lcm0 : props   where@@ -152,6 +168,8 @@ myEuclideanLaws proxy = testGroup tpclss $ map (uncurry testProperty) props   where     Laws tpclss props = euclideanLaws proxy++#endif  myIsListLaws :: (Eq a, IsList a, Arbitrary a, Show a, Show (Item a), Arbitrary (Item a)) => Proxy a -> TestTree myIsListLaws proxy = testGroup tpclss $ map (uncurry testProperty) props

− test/doctests.hs

@@ -1,4 +0,0 @@-import Test.DocTest (doctest)--main :: IO ()-main = doctest ["src"]