diff --git a/README.md b/README.md --- a/README.md +++ b/README.md @@ -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). diff --git a/Setup.hs b/Setup.hs deleted file mode 100644 --- a/Setup.hs +++ /dev/null @@ -1,2 +0,0 @@ -import Distribution.Simple -main = defaultMain diff --git a/bench/Bench.hs b/bench/Bench.hs --- a/bench/Bench.hs +++ b/bench/Bench.hs @@ -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 ] diff --git a/changelog.md b/changelog.md --- a/changelog.md +++ b/changelog.md @@ -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 diff --git a/poly.cabal b/poly.cabal --- a/poly.cabal +++ b/poly.cabal @@ -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 diff --git a/src/Data/Poly.hs b/src/Data/Poly.hs --- a/src/Data/Poly.hs +++ b/src/Data/Poly.hs @@ -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 () diff --git a/src/Data/Poly/Internal/Convert.hs b/src/Data/Poly/Internal/Convert.hs --- a/src/Data/Poly/Internal/Convert.hs +++ b/src/Data/Poly/Internal/Convert.hs @@ -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 diff --git a/src/Data/Poly/Internal/Dense.hs b/src/Data/Poly/Internal/Dense.hs --- a/src/Data/Poly/Internal/Dense.hs +++ b/src/Data/Poly/Internal/Dense.hs @@ -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 diff --git a/src/Data/Poly/Internal/Dense/DFT.hs b/src/Data/Poly/Internal/Dense/DFT.hs --- a/src/Data/Poly/Internal/Dense/DFT.hs +++ b/src/Data/Poly/Internal/Dense/DFT.hs @@ -26,6 +26,8 @@ -- | -- \( 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 -- \( 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 #-} diff --git a/src/Data/Poly/Internal/Dense/Field.hs b/src/Data/Poly/Internal/Dense/Field.hs --- a/src/Data/Poly/Internal/Dense/Field.hs +++ b/src/Data/Poly/Internal/Dense/Field.hs @@ -4,7 +4,7 @@ -- Licence: BSD3 -- Maintainer: Andrew Lelechenko -- --- 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 diff --git a/src/Data/Poly/Internal/Dense/GcdDomain.hs b/src/Data/Poly/Internal/Dense/GcdDomain.hs --- a/src/Data/Poly/Internal/Dense/GcdDomain.hs +++ b/src/Data/Poly/Internal/Dense/GcdDomain.hs @@ -4,7 +4,7 @@ -- Licence: BSD3 -- Maintainer: Andrew Lelechenko -- --- 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 diff --git a/src/Data/Poly/Internal/Dense/Laurent.hs b/src/Data/Poly/Internal/Dense/Laurent.hs --- a/src/Data/Poly/Internal/Dense/Laurent.hs +++ b/src/Data/Poly/Internal/Dense/Laurent.hs @@ -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 diff --git a/src/Data/Poly/Internal/Multi.hs b/src/Data/Poly/Internal/Multi.hs --- a/src/Data/Poly/Internal/Multi.hs +++ b/src/Data/Poly/Internal/Multi.hs @@ -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) diff --git a/src/Data/Poly/Internal/Multi/Core.hs b/src/Data/Poly/Internal/Multi/Core.hs --- a/src/Data/Poly/Internal/Multi/Core.hs +++ b/src/Data/Poly/Internal/Multi/Core.hs @@ -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 diff --git a/src/Data/Poly/Internal/Multi/Field.hs b/src/Data/Poly/Internal/Multi/Field.hs --- a/src/Data/Poly/Internal/Multi/Field.hs +++ b/src/Data/Poly/Internal/Multi/Field.hs @@ -4,7 +4,7 @@ -- Licence: BSD3 -- Maintainer: Andrew Lelechenko -- --- 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 #-} diff --git a/src/Data/Poly/Internal/Multi/GcdDomain.hs b/src/Data/Poly/Internal/Multi/GcdDomain.hs --- a/src/Data/Poly/Internal/Multi/GcdDomain.hs +++ b/src/Data/Poly/Internal/Multi/GcdDomain.hs @@ -4,7 +4,7 @@ -- Licence: BSD3 -- Maintainer: Andrew Lelechenko -- --- 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)) diff --git a/src/Data/Poly/Internal/Multi/Laurent.hs b/src/Data/Poly/Internal/Multi/Laurent.hs --- a/src/Data/Poly/Internal/Multi/Laurent.hs +++ b/src/Data/Poly/Internal/Multi/Laurent.hs @@ -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 -- | -- 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)) diff --git a/src/Data/Poly/Laurent.hs b/src/Data/Poly/Laurent.hs --- a/src/Data/Poly/Laurent.hs +++ b/src/Data/Poly/Laurent.hs @@ -6,6 +6,8 @@ -- -- . -- +-- @since 0.4.0.0 +-- {-# LANGUAGE PatternSynonyms #-} diff --git a/src/Data/Poly/Multi.hs b/src/Data/Poly/Multi.hs --- a/src/Data/Poly/Multi.hs +++ b/src/Data/Poly/Multi.hs @@ -6,6 +6,7 @@ -- -- Sparse multivariate polynomials with 'Num' instance. -- +-- @since 0.5.0.0 {-# LANGUAGE DataKinds #-} {-# LANGUAGE FlexibleContexts #-} diff --git a/src/Data/Poly/Multi/Semiring.hs b/src/Data/Poly/Multi/Semiring.hs --- a/src/Data/Poly/Multi/Semiring.hs +++ b/src/Data/Poly/Multi/Semiring.hs @@ -4,8 +4,9 @@ -- Licence: BSD3 -- Maintainer: Andrew Lelechenko -- --- 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 diff --git a/src/Data/Poly/Orthogonal.hs b/src/Data/Poly/Orthogonal.hs --- a/src/Data/Poly/Orthogonal.hs +++ b/src/Data/Poly/Orthogonal.hs @@ -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) -- | . +-- +-- @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 -- | . +-- +-- @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 @@ -- | -- 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 @@ -- | -- 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' . -- --- >>> 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' . -- --- >>> 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 @@ -- | . -- --- >>> 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 -- | +-- +-- @since 0.4.0.0 laguerreGen :: (Eq a, Field a, Vector v a) => a -> [Poly v a] laguerreGen a = xs where diff --git a/src/Data/Poly/Semiring.hs b/src/Data/Poly/Semiring.hs --- a/src/Data/Poly/Semiring.hs +++ b/src/Data/Poly/Semiring.hs @@ -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 @@ -- . -- 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 diff --git a/src/Data/Poly/Sparse.hs b/src/Data/Poly/Sparse.hs --- a/src/Data/Poly/Sparse.hs +++ b/src/Data/Poly/Sparse.hs @@ -4,8 +4,10 @@ -- Licence: BSD3 -- Maintainer: Andrew Lelechenko -- --- 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 #-} diff --git a/src/Data/Poly/Sparse/Laurent.hs b/src/Data/Poly/Sparse/Laurent.hs --- a/src/Data/Poly/Sparse/Laurent.hs +++ b/src/Data/Poly/Sparse/Laurent.hs @@ -7,6 +7,7 @@ -- Sparse -- . -- +-- @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 #-} diff --git a/src/Data/Poly/Sparse/Semiring.hs b/src/Data/Poly/Sparse/Semiring.hs --- a/src/Data/Poly/Sparse/Semiring.hs +++ b/src/Data/Poly/Sparse/Semiring.hs @@ -4,8 +4,10 @@ -- Licence: BSD3 -- Maintainer: Andrew Lelechenko -- --- 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 #-} diff --git a/test/Dense.hs b/test/Dense.hs --- a/test/Dense.hs +++ b/test/Dense.hs @@ -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 ] diff --git a/test/DenseLaurent.hs b/test/DenseLaurent.hs --- a/test/DenseLaurent.hs +++ b/test/DenseLaurent.hs @@ -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 = diff --git a/test/Main.hs b/test/Main.hs --- a/test/Main.hs +++ b/test/Main.hs @@ -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 ] diff --git a/test/Multi.hs b/test/Multi.hs --- a/test/Multi.hs +++ b/test/Multi.hs @@ -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 = diff --git a/test/MultiLaurent.hs b/test/MultiLaurent.hs --- a/test/MultiLaurent.hs +++ b/test/MultiLaurent.hs @@ -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 = diff --git a/test/Sparse.hs b/test/Sparse.hs --- a/test/Sparse.hs +++ b/test/Sparse.hs @@ -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 = diff --git a/test/SparseLaurent.hs b/test/SparseLaurent.hs --- a/test/SparseLaurent.hs +++ b/test/SparseLaurent.hs @@ -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 = diff --git a/test/TestUtils.hs b/test/TestUtils.hs --- a/test/TestUtils.hs +++ b/test/TestUtils.hs @@ -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 diff --git a/test/doctests.hs b/test/doctests.hs deleted file mode 100644 --- a/test/doctests.hs +++ /dev/null @@ -1,4 +0,0 @@ -import Test.DocTest (doctest) - -main :: IO () -main = doctest ["src"]