poly 0.4.0.0 → 0.5.0.0

raw patch · 42 files changed

+3863/−1715 lines, 42 filesdep +doctestdep +finite-typelitsdep +vector-sizeddep ~basedep ~modPVP ok

version bump matches the API change (PVP)

Dependencies added: doctest, finite-typelits, vector-sized

Dependency ranges changed: base, mod

API changes (from Hackage documentation)

Files

• README.md +16/−10
• bench/DenseBench.hs +14/−22
• bench/SparseBench.hs +16/−15
• changelog.md +12/−0
• poly.cabal +43/−14
• src/Data/Poly.hs +5/−3
• src/Data/Poly/Internal/Convert.hs +88/−0
• src/Data/Poly/Internal/Dense.hs +25/−10
• src/Data/Poly/Internal/Dense/DFT.hs +80/−0
• src/Data/Poly/Internal/Dense/Field.hs +35/−42
• src/Data/Poly/Internal/Dense/GcdDomain.hs +66/−41
• src/Data/Poly/Internal/Dense/Laurent.hs +288/−0
• src/Data/Poly/Internal/Multi.hs +553/−0
• src/Data/Poly/Internal/Multi/Core.hs +311/−0
• src/Data/Poly/Internal/Multi/Field.hs +73/−0
• src/Data/Poly/Internal/Multi/GcdDomain.hs +179/−0
• src/Data/Poly/Internal/Multi/Laurent.hs +507/−0
• src/Data/Poly/Internal/PolyOverField.hs +0/−46
• src/Data/Poly/Internal/Sparse.hs +0/−583
• src/Data/Poly/Internal/Sparse/Field.hs +0/−56
• src/Data/Poly/Internal/Sparse/GcdDomain.hs +0/−74
• src/Data/Poly/Laurent.hs +1/−257
• src/Data/Poly/Multi.hs +35/−0
• src/Data/Poly/Multi/Laurent.hs +32/−0
• src/Data/Poly/Multi/Semiring.hs +157/−0
• src/Data/Poly/Orthogonal.hs +2/−2
• src/Data/Poly/Semiring.hs +56/−5
• src/Data/Poly/Sparse.hs +99/−5
• src/Data/Poly/Sparse/Laurent.hs +62/−235
• src/Data/Poly/Sparse/Semiring.hs +81/−24
• test/DFT.hs +69/−0
• test/Dense.hs +99/−68
• test/DenseLaurent.hs +62/−49
• test/Main.hs +11/−5
• test/Multi.hs +307/−0
• test/MultiLaurent.hs +222/−0
• test/Orthogonal.hs +1/−1
• test/Quaternion.hs +11/−13
• test/Sparse.hs +99/−66
• test/SparseLaurent.hs +73/−57
• test/TestUtils.hs +69/−12
• test/doctests.hs +4/−0

@@ -1,8 +1,6 @@-# 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)--+# 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) -Haskell library for univariate polynomials, backed by `Vector`.+Haskell library for univariate and multivariate polynomials, backed by `Vector`.  ```haskell > (X + 1) + (X - 1) :: VPoly Integer@@ -116,18 +114,26 @@  ## Flavours -The same API is exposed in four flavours:--* `Data.Poly` provides dense polynomials with `Num`-based interface.+* `Data.Poly` provides dense univariate polynomials with `Num`-based interface.   This is a default choice for most users. -* `Data.Poly.Semiring` provides dense polynomials with `Semiring`-based interface.+* `Data.Poly.Semiring` provides dense univariate polynomials with `Semiring`-based interface. -* `Data.Poly.Sparse` provides sparse polynomials with `Num`-based interface.+* `Data.Poly.Laurent` provides dense univariate Laurent polynomials with `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   because of a more readable `Show` instance, skipping zero coefficients. -* `Data.Poly.Sparse.Semiring` provides sparse polynomials with `Semiring`-based interface.+* `Data.Poly.Sparse.Semiring` provides sparse univariate polynomials with `Semiring`-based interface.++* `Data.Poly.Sparse.Laurent` provides sparse univariate Laurent polynomials with `Semiring`-based interface.++* `Data.Poly.Multi` provides sparse multivariate polynomials with `Num`-based interface.++* `Data.Poly.Multi.Semiring` provides sparse multivariate polynomials with `Semiring`-based interface.++* `Data.Poly.Multi.Laurent` provides sparse multivariate Laurent polynomials with `Semiring`-based interface.  All flavours are available backed by boxed or unboxed vectors. 

@@ -1,4 +1,3 @@-{-# LANGUAGE GeneralizedNewtypeDeriving #-} {-# LANGUAGE RankNTypes                 #-} {-# LANGUAGE TypeApplications           #-} @@ -8,7 +7,7 @@  import Prelude hiding (quotRem, gcd) import Gauge.Main-import Data.Euclidean (Euclidean(..), GcdDomain(..), Field)+import Data.Euclidean (Euclidean(..), GcdDomain(..)) import Data.Poly import qualified Data.Poly.Semiring as S (toPoly) import Data.Semiring (Semiring(..), Ring, Mod2(..))@@ -23,10 +22,10 @@   , map benchEval     [100, 1000, 10000]   , map benchDeriv    [100, 1000, 10000]   , map benchIntegral [100, 1000, 10000]-  , map benchQuotRem    [10, 100]-  , map benchGcd        [10, 100]-  , map benchGcdFracRat [10, 20, 40]-  , map benchGcdFracM   [10, 100, 1000]+  , map benchQuotRem  [10, 100]+  , map benchGcd      [10, 100]+  , map benchGcdRat   [10, 20, 40]+  , map benchGcdM     [10, 100, 1000]   ]  benchAdd :: Int -> Benchmark@@ -48,13 +47,13 @@ benchQuotRem k = bench ("quotRem/" ++ show k) $ nf doQuotRem k  benchGcd :: Int -> Benchmark-benchGcd k = bench ("gcd/" ++ show k) $ nf doGcd k+benchGcd k = bench ("gcd/Integer/" ++ show k) $ nf (doGcd @Integer) k -benchGcdFracRat :: Int -> Benchmark-benchGcdFracRat k = bench ("gcdFrac/Rational/" ++ show k) $ nf (doGcdFrac @Rational) k+benchGcdRat :: Int -> Benchmark+benchGcdRat k = bench ("gcd/Rational/" ++ show k) $ nf (doGcd @Rational) k -benchGcdFracM :: Int -> Benchmark-benchGcdFracM k = bench ("gcdFrac/Mod2/" ++ show k) $ nf (getMod2 . doGcdFrac @Mod2) k+benchGcdM :: Int -> Benchmark+benchGcdM k = bench ("gcd/Mod2/" ++ show k) $ nf (getMod2 . doGcd @Mod2) k  doBinOp :: (forall a. Num a => a -> a -> a) -> Int -> Int doBinOp op n = U.sum zs@@ -94,16 +93,9 @@     ys = toPoly $ U.generate n       gen2     (qs, rs) = xs `quotRem` ys -doGcd :: Int -> Integer-doGcd n = V.sum gs+doGcd :: (Eq a, Ring a, GcdDomain a) => Int -> a+doGcd n = V.foldl' plus zero gs   where-    xs = toPoly $ V.generate n gen1-    ys = toPoly $ V.generate n gen2+    xs = S.toPoly $ V.generate n gen1+    ys = S.toPoly $ V.generate n gen2     gs = unPoly $ xs `gcd` ys--doGcdFrac :: (Eq a, Field a) => Int -> a-doGcdFrac n = V.foldl' plus zero gs-  where-    xs = PolyOverField $ S.toPoly $ V.generate n gen1-    ys = PolyOverField $ S.toPoly $ V.generate n gen2-    gs = unPoly $ unPolyOverField $ xs `gcd` ys

@@ -12,11 +12,12 @@  benchSuite :: Benchmark benchSuite = bgroup "sparse" $ concat-  [ map benchAdd      $ zip3 tabs vecs2 vecs3-  , map benchMul      $ take 2 $ zip3 tabs vecs2 vecs3-  , map benchEval     $ zip tabs vecs2-  , map benchDeriv    $ zip tabs vecs2-  , map benchIntegral $ zip tabs vecs2'+  [ zipWith3 benchAdd tabs vecs2 vecs3+  , take 2+  $ zipWith3 benchMul tabs vecs2 vecs3+  , zipWith benchEval tabs vecs2+  , zipWith benchDeriv tabs vecs2+  , zipWith benchIntegral tabs vecs2'   ]  tabs :: [Int]@@ -34,20 +35,20 @@ vecs3 = flip map tabs $   \n -> toPoly $ U.generate n (\i -> (fromIntegral i ^ 3, i * 3)) -benchAdd :: (Int, UPoly Int, UPoly Int) -> Benchmark-benchAdd (k, xs, ys) = bench ("add/" ++ show k) $ nf (doBinOp (+) xs) ys+benchAdd :: Int -> UPoly Int -> UPoly Int -> Benchmark+benchAdd k xs ys = bench ("add/" ++ show k) $ nf (doBinOp (+) xs) ys -benchMul :: (Int, UPoly Int, UPoly Int) -> Benchmark-benchMul (k, xs, ys) = bench ("mul/" ++ show k) $ nf (doBinOp (*) xs) ys+benchMul :: Int -> UPoly Int -> UPoly Int -> Benchmark+benchMul k xs ys = bench ("mul/" ++ show k) $ nf (doBinOp (*) xs) ys -benchEval :: (Int, UPoly Int) -> Benchmark-benchEval (k, xs) = bench ("eval/" ++ show k) $ nf doEval xs+benchEval :: Int -> UPoly Int -> Benchmark+benchEval k xs = bench ("eval/" ++ show k) $ nf doEval xs -benchDeriv :: (Int, UPoly Int) -> Benchmark-benchDeriv (k, xs) = bench ("deriv/" ++ show k) $ nf doDeriv xs+benchDeriv :: Int -> UPoly Int -> Benchmark+benchDeriv k xs = bench ("deriv/" ++ show k) $ nf doDeriv xs -benchIntegral :: (Int, UPoly Double) -> Benchmark-benchIntegral (k, xs) = bench ("integral/" ++ show k) $ nf doIntegral xs+benchIntegral :: Int -> UPoly Double -> Benchmark+benchIntegral k xs = bench ("integral/" ++ show k) $ nf doIntegral xs  doBinOp :: (forall a. Num a => a -> a -> a) -> UPoly Int -> UPoly Int -> Int doBinOp op xs ys = U.foldl' (\acc (_, x) -> acc + x) 0 zs

@@ -1,3 +1,15 @@+# 0.5.0.0++* Change definition of `Data.Euclidean.degree`+  to coincide with the degree of polynomial.+* 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.+* Add function `quotRemFractional`.+* Add an experimental implementation of the discrete Fourier transform.+* Add conversion functions between dense and sparse polynomials.+ # 0.4.0.0  * Implement Laurent polynomials.

@@ -1,5 +1,5 @@ name: poly-version: 0.4.0.0+version: 0.5.0.0 synopsis: Polynomials description:   Polynomials backed by `Vector`.@@ -12,7 +12,7 @@ category: Math, Numerical build-type: Simple cabal-version: >=1.10-tested-with: GHC ==8.0.2 GHC ==8.2.2 GHC ==8.4.4 GHC ==8.6.5 GHC ==8.8.3 GHC ==8.10.1+tested-with: GHC ==8.2.2 GHC ==8.4.4 GHC ==8.6.5 GHC ==8.8.4 GHC ==8.10.2 extra-source-files:   changelog.md   README.md@@ -26,28 +26,41 @@   exposed-modules:     Data.Poly     Data.Poly.Laurent-    Data.Poly.Orthogonal     Data.Poly.Semiring+    Data.Poly.Orthogonal+     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.Internal.Dense     Data.Poly.Internal.Dense.Field+    Data.Poly.Internal.Dense.DFT     Data.Poly.Internal.Dense.GcdDomain-    Data.Poly.Internal.PolyOverField-    Data.Poly.Internal.Sparse-    Data.Poly.Internal.Sparse.Field-    Data.Poly.Internal.Sparse.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   build-depends:-    base >= 4.9 && < 5,+    base >= 4.10 && < 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-algorithms >= 0.8.0.3,+    vector-sized >= 1.1   default-language: Haskell2010-  ghc-options: -Wall -Wcompat+  ghc-options: -Wall -Wcompat -Wredundant-constraints  test-suite poly-tests   type: exitcode-stdio-1.0@@ -55,30 +68,46 @@   other-modules:     Dense     DenseLaurent+    DFT+    Multi+    MultiLaurent     Orthogonal     Quaternion     Sparse     SparseLaurent     TestUtils   build-depends:-    base >=4.9 && <5,-    mod,+    base >=4.10 && <5,+    finite-typelits,+    mod >=0.1.2,     poly,     QuickCheck >=2.12,     quickcheck-classes >=0.6.3,     semirings >= 0.5.2,     tasty >= 0.11,     tasty-quickcheck >= 0.8,-    vector >= 0.12.0.2+    vector >= 0.12.0.2,+    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+ benchmark poly-gauge   build-depends:-    base >=4.9 && <5,+    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

@@ -23,10 +23,12 @@   , subst   , deriv   , integral-  , PolyOverField(..)+  , quotRemFractional+  , denseToSparse+  , sparseToDense   ) where +import Data.Poly.Internal.Convert import Data.Poly.Internal.Dense-import Data.Poly.Internal.Dense.Field ()+import Data.Poly.Internal.Dense.Field (quotRemFractional) import Data.Poly.Internal.Dense.GcdDomain ()-import Data.Poly.Internal.PolyOverField

@@ -0,0 +1,88 @@+-- |+-- Module:      Data.Poly.Internal.Convert+-- Copyright:   (c) 2020 Andrew Lelechenko+-- Licence:     BSD3+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+--+-- Conversions between polynomials.+--++{-# LANGUAGE DataKinds        #-}+{-# LANGUAGE FlexibleContexts #-}++module Data.Poly.Internal.Convert+  ( denseToSparse+  , denseToSparse'+  , sparseToDense+  , sparseToDense'+  ) where++import Control.Monad.ST+import Data.Semiring (Semiring(..))+import qualified Data.Vector.Generic as G+import qualified Data.Vector.Generic.Mutable as MG+import qualified Data.Vector.Unboxed.Sized as SU++import qualified Data.Poly.Internal.Dense as Dense+import qualified Data.Poly.Internal.Multi as Sparse++-- | Convert from dense to sparse polynomials.+--+-- >>> :set -XFlexibleContexts+-- >>> denseToSparse (1 + Data.Poly.X^2) :: Data.Poly.Sparse.UPoly Int+-- 1 * X^2 + 1+denseToSparse+  :: (Eq a, Num a, G.Vector v a, G.Vector v (SU.Vector 1 Word, a))+  => Dense.Poly v a+  -> Sparse.Poly v a+denseToSparse = denseToSparseInternal 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' = denseToSparseInternal zero++denseToSparseInternal+  :: (Eq a, G.Vector v a, G.Vector v (SU.Vector 1 Word, a))+  => a+  -> Dense.Poly v a+  -> Sparse.Poly v a+denseToSparseInternal z = Sparse.MultiPoly . G.imapMaybe (\i c -> if c == z then Nothing else Just (fromIntegral i, c)) . Dense.unPoly++-- | Convert from sparse to dense polynomials.+--+-- >>> :set -XFlexibleContexts+-- >>> sparseToDense (1 + Data.Poly.Sparse.X^2) :: Data.Poly.UPoly Int+-- 1 * X^2 + 0 * X + 1+sparseToDense+  :: (Num a, G.Vector v a, G.Vector v (SU.Vector 1 Word, a))+  => Sparse.Poly v a+  -> Dense.Poly v a+sparseToDense = sparseToDenseInternal 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' = sparseToDenseInternal zero++sparseToDenseInternal+  :: (G.Vector v a, G.Vector v (SU.Vector 1 Word, a))+  => a+  -> Sparse.Poly v a+  -> Dense.Poly v a+sparseToDenseInternal z (Sparse.MultiPoly xs)+  | G.null xs = Dense.Poly G.empty+  | otherwise = runST $ do+  let len = fromIntegral (SU.head (fst (G.unsafeLast xs)) + 1)+  ys <- MG.unsafeNew len+  MG.set ys z+  let go xi yi+        | xi >= G.length xs = pure ()+        | (yi', c) <- G.unsafeIndex xs xi+        , yi == fromIntegral (SU.head yi')+        = MG.unsafeWrite ys yi c >> go (xi + 1) (yi + 1)+        | otherwise = go xi (yi + 1)+  go 0 0+  Dense.Poly <$> G.unsafeFreeze ys

@@ -42,13 +42,14 @@   , integral'   ) where -import Prelude hiding (quotRem, quot, rem, gcd, lcm, (^))+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 qualified Data.Semiring as Semiring@@ -74,7 +75,7 @@ -- 'Ord' instance does not make much sense mathematically, -- it is defined only for the sake of 'Data.Set.Set', 'Data.Map.Map', etc. ---newtype Poly v a = Poly+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).@@ -99,9 +100,11 @@       $ intersperse (showString " + ")       $ G.ifoldl (\acc i c -> showCoeff i c : acc) [] xs     where+      -- Powers are guaranteed to be non-negative+      showCoeff :: Int -> a -> String -> String       showCoeff 0 c = showsPrec 7 c       showCoeff 1 c = showsPrec 7 c . showString " * X"-      showCoeff i c = showsPrec 7 c . showString " * X^" . showsPrec 7 i+      showCoeff i c = showsPrec 7 c . showString (" * X^" ++ show i)  -- | Polynomials backed by boxed vectors. type VPoly = Poly V.Vector@@ -328,7 +331,7 @@ {-# INLINE monomial' #-}  scaleInternal-  :: (Eq a, G.Vector v a)+  :: G.Vector v a   => a   -> (a -> a -> a)   -> Word@@ -453,23 +456,35 @@ {-# INLINABLE integral' #-}  -- | Create an identity polynomial.-pattern X :: (Eq a, Num a, G.Vector v a, Eq (v a)) => Poly v a-pattern X <- ((==) var -> True)+pattern X :: (Eq a, Num a, G.Vector v a) => Poly v a+pattern X <- (isVar -> True)   where X = var -var :: forall a v. (Eq a, Num a, G.Vector v a, Eq (v a)) => Poly v a+var :: forall a v. (Eq a, Num a, G.Vector v a) => Poly v a var   | (1 :: a) == 0 = Poly G.empty   | otherwise     = Poly $ G.fromList [0, 1] {-# INLINE var #-} +isVar :: forall v a. (Eq a, Num a, G.Vector v a) => Poly v a -> Bool+isVar (Poly xs)+  | (1 :: a) == 0 = G.null xs+  | otherwise     = G.length xs == 2 && xs G.! 0 == 0 && xs G.! 1 == 1+{-# INLINE isVar #-}+ -- | Create an identity polynomial.-pattern X' :: (Eq a, Semiring a, G.Vector v a, Eq (v a)) => Poly v a-pattern X' <- ((==) var' -> True)+pattern X' :: (Eq a, Semiring a, G.Vector v a) => Poly v a+pattern X' <- (isVar' -> True)   where X' = var' -var' :: forall a v. (Eq a, Semiring a, G.Vector v a, Eq (v a)) => Poly v a+var' :: forall a v. (Eq a, Semiring a, G.Vector v a) => Poly v a var'   | (one :: a) == zero = Poly G.empty   | otherwise          = Poly $ G.fromList [zero, one] {-# INLINE var' #-}++isVar' :: forall v a. (Eq a, Semiring a, G.Vector v a) => Poly v a -> Bool+isVar' (Poly xs)+  | (one :: a) == zero = G.null xs+  | otherwise          = G.length xs == 2 && xs G.! 0 == zero && xs G.! 1 == one+{-# INLINE isVar' #-}

@@ -0,0 +1,80 @@+-- |+-- Module:      Data.Poly.Internal.Dense.FFT+-- Copyright:   (c) 2020 Andrew Lelechenko+-- Licence:     BSD3+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+--+-- Discrete Fourier transform.+--++{-# LANGUAGE BangPatterns        #-}+{-# LANGUAGE ScopedTypeVariables #-}++module Data.Poly.Internal.Dense.DFT+  ( dft+  , inverseDft+  ) where++import Prelude hiding (recip, fromIntegral)+import Control.Monad.ST+import Data.Bits hiding (shift)+import Data.Foldable+import Data.Semiring (Semiring(..), Ring(..), minus, fromIntegral)+import Data.Field (Field, recip)+import qualified Data.Vector.Generic as G+import qualified Data.Vector.Generic.Mutable as MG++-- | <https://en.wikipedia.org/wiki/Fast_Fourier_transform Discrete Fourier transform>+-- \( y_k = \sum_{j=0}^{N-1} x_j \sqrt[N]{1}^{jk} \).+dft+  :: (Ring a, G.Vector v a)+  => a   -- ^ primitive root \( \sqrt[N]{1} \), otherwise behaviour is undefined+  -> v a -- ^ \( \{ x_k \}_{k=0}^{N-1} \) (currently only  \( N = 2^n \) is supported)+  -> v a -- ^ \( \{ y_k \}_{k=0}^{N-1} \)+dft primRoot (xs :: v a)+  | popCount nn /= 1 = error "dft: only vectors of length 2^n are supported"+  | otherwise = go 0 0+  where+    nn = G.length xs+    n = countTrailingZeros nn++    roots :: v a+    roots = G.iterateN+      (1 `unsafeShiftL` (n - 1))+      (\x -> x `seq` (x `times` primRoot))+      one++    go !offset !shift+      | shift >= n = G.unsafeSlice offset 1 xs+      | otherwise = runST $ do+        let halfLen = 1 `unsafeShiftL` (n - shift - 1)+            ys0 = go offset (shift + 1)+            ys1 = go (offset + 1 `unsafeShiftL` shift) (shift + 1)+        ys <- MG.new (halfLen `unsafeShiftL` 1)++        -- This corresponds to k = 0 in the loop below.+        -- It improves performance by avoiding multiplication+        -- by roots V.! 0 = 1.+        let y00 = G.unsafeIndex ys0 0+            y10 = G.unsafeIndex ys1 0+        MG.unsafeWrite ys 0       $! y00 `plus`  y10+        MG.unsafeWrite ys halfLen $! y00 `minus` y10++        forM_ [1..halfLen - 1] $ \k -> do+          let y0 = G.unsafeIndex ys0 k+              y1 = G.unsafeIndex ys1 k `times`+                   G.unsafeIndex roots (k `unsafeShiftL` shift)+          MG.unsafeWrite ys k             $! y0 `plus`  y1+          MG.unsafeWrite ys (k + halfLen) $! y0 `minus` y1+        G.unsafeFreeze ys++-- | Inverse <https://en.wikipedia.org/wiki/Fast_Fourier_transform discrete Fourier transform>+-- \( x_k = {1\over N} \sum_{j=0}^{N-1} y_j \sqrt[N]{1}^{-jk} \).+inverseDft+  :: (Field a, G.Vector v a)+  => a   -- ^ primitive root \( \sqrt[N]{1} \), otherwise behaviour is undefined+  -> v a -- ^ \( \{ y_k \}_{k=0}^{N-1} \) (currently only  \( N = 2^n \) is supported)+  -> v a -- ^ \( \{ x_k \}_{k=0}^{N-1} \)+inverseDft primRoot ys = G.map (`times` invN) $ dft (recip primRoot) ys+  where+    invN = recip $ fromIntegral $ G.length ys

@@ -9,23 +9,21 @@  {-# LANGUAGE ConstraintKinds            #-} {-# LANGUAGE FlexibleInstances          #-}-{-# LANGUAGE PatternSynonyms            #-} {-# LANGUAGE ScopedTypeVariables        #-} {-# LANGUAGE TypeFamilies               #-}  {-# OPTIONS_GHC -fno-warn-orphans #-}  module Data.Poly.Internal.Dense.Field-  ( fieldGcd+  ( quotRemFractional   ) where -import Prelude hiding (quotRem, quot, rem, gcd, recip)+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.Field (recip) import Data.Semiring (times, minus, zero, one) import qualified Data.Vector.Generic as G import qualified Data.Vector.Generic.Mutable as MG@@ -33,42 +31,60 @@ import Data.Poly.Internal.Dense import Data.Poly.Internal.Dense.GcdDomain () -instance (Eq a, Eq (v a), Field a, G.Vector v a) => Euclidean (Poly v a) where-  degree (Poly xs) = fromIntegral (G.length xs)+-- | Note that 'degree' 0 = 0.+instance (Eq a, Field a, G.Vector v a) => Euclidean (Poly v a) where+  degree (Poly xs)+    | G.null xs = 0+    | otherwise = fromIntegral (G.length xs - 1)    quotRem (Poly xs) (Poly ys) = (toPoly' qs, toPoly' rs)     where-      (qs, rs) = quotientAndRemainder xs ys+      (qs, rs) = quotientAndRemainder zero (== one) minus times (one `quot`) xs ys   {-# INLINE quotRem #-}    rem (Poly xs) (Poly ys) = toPoly' $ remainder xs ys   {-# INLINE rem #-} +-- | Polynomial division with remainder.+--+-- >>> quotRemFractional (X^3 + 2) (X^2 - 1 :: UPoly Double)+-- (1.0 * X + 0.0,1.0 * X + 2.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+    (qs, rs) = quotientAndRemainder 0 (== 1) (-) (*) recip xs ys+{-# INLINE quotRemFractional #-}+ quotientAndRemainder-  :: (Eq a, Field a, G.Vector v a)-  => v a-  -> v a+  :: (Eq a, G.Vector v a)+  => a             -- ^ zero+  -> (a -> Bool)   -- ^ is one?+  -> (a -> a -> a) -- ^ subtract+  -> (a -> a -> a) -- ^ multiply+  -> (a -> a)      -- ^ invert+  -> v a           -- ^ dividend+  -> v a           -- ^ divisor   -> (v a, v a)-quotientAndRemainder xs ys+quotientAndRemainder zer isOne sub mul inv xs ys   | lenXs < lenYs = (G.empty, xs)   | lenYs == 0 = throw DivideByZero-  | lenYs == 1 = let invY = recip (G.unsafeHead ys) in-                 (G.map (`times` invY) xs, G.empty)+  | lenYs == 1 = let invY = inv (G.unsafeHead ys) in+                 (G.map (`mul` invY) xs, G.empty)   | otherwise = runST $ do     qs <- MG.unsafeNew lenQs     rs <- MG.unsafeNew lenXs     G.unsafeCopy rs xs     let yLast = G.unsafeLast ys-        invYLast = recip yLast+        invYLast = inv yLast     forM_ [lenQs - 1, lenQs - 2 .. 0] $ \i -> do       r <- MG.unsafeRead rs (lenYs - 1 + i)-      let q = if yLast == one then r else r `times` invYLast+      let q = if isOne yLast then r else r `mul` invYLast       MG.unsafeWrite qs i q-      MG.unsafeWrite rs (lenYs - 1 + i) zero+      MG.unsafeWrite rs (lenYs - 1 + i) zer       forM_ [0 .. lenYs - 2] $ \k -> do         let y = G.unsafeIndex ys k-        when (y /= zero) $-          MG.unsafeModify rs (\c -> c `minus` q `times` y) (i + k)+        when (y /= zer) $+          MG.unsafeModify rs (\c -> c `sub` (q `mul` y)) (i + k)     let rs' = MG.unsafeSlice 0 lenYs rs     (,) <$> G.unsafeFreeze qs <*> G.unsafeFreeze rs'   where@@ -102,7 +118,7 @@   | lenYs == 1 = MG.set xs zero   | otherwise = do     yLast <- MG.unsafeRead ys (lenYs - 1)-    let invYLast = recip yLast+    let invYLast = one `quot` yLast     forM_ [lenQs - 1, lenQs - 2 .. 0] $ \i -> do       r <- MG.unsafeRead xs (lenYs - 1 + i)       MG.unsafeWrite xs (lenYs - 1 + i) zero@@ -116,26 +132,3 @@     lenYs = MG.length ys     lenQs = lenXs - lenYs + 1 {-# INLINABLE remainderM #-}--fieldGcd-  :: (Eq a, Field a, G.Vector v a)-  => Poly v a-  -> Poly v a-  -> Poly v a-fieldGcd (Poly xs) (Poly ys) = toPoly' $ runST $ do-  xs' <- G.thaw xs-  ys' <- G.thaw ys-  gcdM xs' ys'-{-# INLINE fieldGcd #-}--gcdM-  :: (PrimMonad m, Eq a, Field a, G.Vector v a)-  => G.Mutable v (PrimState m) a-  -> G.Mutable v (PrimState m) a-  -> m (v a)-gcdM xs ys = do-  ys' <- dropWhileEndM (== zero) ys-  if MG.null ys' then G.unsafeFreeze xs else do-    remainderM xs ys'-    gcdM ys' xs-{-# INLINE gcdM #-}

@@ -8,7 +8,7 @@ --  {-# LANGUAGE FlexibleInstances          #-}-{-# LANGUAGE PatternSynonyms            #-}+{-# LANGUAGE MultiWayIf                 #-} {-# LANGUAGE ScopedTypeVariables        #-} {-# LANGUAGE TypeFamilies               #-} @@ -23,25 +23,31 @@ import Control.Monad.Primitive import Control.Monad.ST import Data.Euclidean+import Data.Maybe import Data.Semiring (Semiring(..), Ring(), isZero, minus) import qualified Data.Vector.Generic as G import qualified Data.Vector.Generic.Mutable as MG  import Data.Poly.Internal.Dense --- | Consider using 'Data.Poly.Semiring.PolyOverField' wrapper,--- which provides a much faster implementation of--- 'Data.Euclidean.gcd' for polynomials over 'Field'.-instance (Eq a, Ring a, GcdDomain a, Eq (v a), G.Vector v a) => GcdDomain (Poly v a) where+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    gcd (Poly xs) (Poly ys)     | G.null xs = Poly ys     | G.null ys = Poly xs+    | G.length xs == 1 = Poly $ G.singleton $ G.foldl' gcd (G.unsafeHead xs) ys+    | G.length ys == 1 = Poly $ G.singleton $ G.foldl' gcd (G.unsafeHead ys) xs     | otherwise = toPoly' $ gcdNonEmpty xs ys   {-# INLINE gcd #-} +  lcm x@(Poly xs) y@(Poly ys)+    | G.null xs || G.null ys = zero+    | otherwise = (x `divide'` gcd x y) `times` y++  coprime x y = isJust (one `divide` gcd x y)+ gcdNonEmpty   :: (Eq a, Ring a, GcdDomain a, G.Vector v a)   => v a@@ -63,14 +69,11 @@     a <- MG.unsafeRead zs' (lenZs - 1)     z <- go a (lenZs - 1) -    let err = error "gcdNonEmpty: violated internal invariant"     forM_ [0 .. lenZs - 1] $ \i ->-      MG.unsafeModify-        zs'-        (\c -> maybe err (`times` xy) (c `divide` z))-        i+      MG.unsafeModify zs'((`times` xy) . (`divide'` z)) i      G.unsafeFreeze zs'+{-# INLINABLE gcdNonEmpty #-}  gcdM   :: (PrimMonad m, Eq a, Ring a, GcdDomain a, G.Vector v a)@@ -85,38 +88,60 @@       lenYs = MG.length ys   xLast <- MG.unsafeRead xs (lenXs - 1)   yLast <- MG.unsafeRead ys (lenYs - 1)-  let z = xLast `lcm` yLast-      zx = case z `divide` xLast of-        Nothing -> error "gcdM: highest coefficient is 0"-        Just t  -> t-      zy = case z `divide` yLast of-        Nothing -> error "gcdM: highest coefficient is 0"-        Just t  -> t+  let z  = xLast `lcm` yLast+      zx = z `divide'` xLast+      zy = z `divide'` yLast -  if lenXs <= lenYs then do-    forM_ [0 .. lenXs - 1] $ \i -> do-      x <- MG.unsafeRead xs i-      MG.unsafeModify-        ys-        (\y -> (y `times` zy) `minus` x `times` zx)-        (i + lenYs - lenXs)-    forM_ [0 .. lenYs - lenXs - 1] $-      MG.unsafeModify ys (`times` zy)-    ys' <- dropWhileEndM isZero ys-    gcdM xs ys'-  else do-    forM_ [0 .. lenYs - 1] $ \i -> do-      y <- MG.unsafeRead ys i-      MG.unsafeModify-        xs-        (\x -> (x `times` zx) `minus` y `times` zy)-        (i + lenXs - lenYs)-    forM_ [0 .. lenXs - lenYs - 1] $-      MG.unsafeModify xs (`times` zx)-    xs' <- dropWhileEndM isZero xs-    gcdM xs' ys+  if+    | lenYs <= lenXs+    , Just xy <- xLast `divide` yLast -> do+      forM_ [0 .. lenYs - 1] $ \i -> do+        y <- MG.unsafeRead ys i+        when (y /= zero) $+          MG.unsafeModify+            xs+            (\x -> x `minus` y `times` xy)+            (i + lenXs - lenYs)+      xs' <- dropWhileEndM isZero xs+      gcdM xs' ys+    | lenXs <= lenYs+    , Just yx <- yLast `divide` xLast -> do+      forM_ [0 .. lenXs - 1] $ \i -> do+        x <- MG.unsafeRead xs i+        when (x /= zero) $+          MG.unsafeModify+            ys+            (\y -> y `minus` x `times` yx)+            (i + lenYs - lenXs)+      ys' <- dropWhileEndM isZero ys+      gcdM xs ys'+    | lenYs <= lenXs -> do+      forM_ [0 .. lenYs - 1] $ \i -> do+        y <- MG.unsafeRead ys i+        MG.unsafeModify+          xs+          (\x -> x `times` zx `minus` y `times` zy)+          (i + lenXs - lenYs)+      forM_ [0 .. lenXs - lenYs - 1] $+        MG.unsafeModify xs (`times` zx)+      xs' <- dropWhileEndM isZero xs+      gcdM xs' ys+    | otherwise -> do+      forM_ [0 .. lenXs - 1] $ \i -> do+        x <- MG.unsafeRead xs i+        MG.unsafeModify+          ys+          (\y -> y `times` zy `minus` x `times` zx)+          (i + lenYs - lenXs)+      forM_ [0 .. lenYs - lenXs - 1] $+        MG.unsafeModify ys (`times` zy)+      ys' <- dropWhileEndM isZero ys+      gcdM xs ys' {-# INLINABLE gcdM #-} +divide' :: GcdDomain a => a -> a -> a+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@@ -130,7 +155,7 @@ {-# INLINE isZeroM #-}  quotient-  :: (Eq a, Eq (v a), Ring a, GcdDomain a, G.Vector v a)+  :: (Eq a, Ring a, GcdDomain a, G.Vector v a)   => v a   -> v a   -> Maybe (v a)@@ -158,7 +183,7 @@               Nothing -> pure Nothing               Just q -> do                 MG.unsafeWrite qs i q-                forM_ [0 .. lenYs - 1] $ \k -> do+                forM_ [0 .. lenYs - 1] $ \k ->                   MG.unsafeModify                     rs                     (\c -> c `minus` q `times` G.unsafeIndex ys k)

@@ -0,0 +1,288 @@+-- |+-- Module:      Data.Poly.Internal.Dense.Laurent+-- Copyright:   (c) 2020 Andrew Lelechenko+-- Licence:     BSD3+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+--+-- <https://en.wikipedia.org/wiki/Laurent_polynomial Laurent polynomials>.+--++{-# LANGUAGE FlexibleInstances          #-}+{-# LANGUAGE KindSignatures             #-}+{-# LANGUAGE PatternSynonyms            #-}+{-# LANGUAGE ScopedTypeVariables        #-}+{-# LANGUAGE ViewPatterns               #-}++module Data.Poly.Internal.Dense.Laurent+  ( Laurent+  , VLaurent+  , ULaurent+  , unLaurent+  , toLaurent+  , leading+  , monomial+  , scale+  , pattern X+  , (^-)+  , eval+  , subst+  , deriv+  ) where++import Prelude hiding (quotRem, quot, rem, gcd, lcm)+import Control.Arrow (first)+import Control.DeepSeq (NFData(..))+import Control.Exception+import Data.Euclidean (GcdDomain(..), Euclidean(..), Field)+import Data.Kind+import Data.List (intersperse)+import Data.Semiring (Semiring(..), Ring())+import qualified Data.Semiring as Semiring+import qualified Data.Vector as V+import qualified Data.Vector.Generic as G+import qualified Data.Vector.Unboxed as U++import Data.Poly.Internal.Dense (Poly(..))+import qualified Data.Poly.Internal.Dense as Dense+import Data.Poly.Internal.Dense.Field ()+import Data.Poly.Internal.Dense.GcdDomain ()++-- | <https://en.wikipedia.org/wiki/Laurent_polynomial Laurent polynomials>+-- of one variable with coefficients from @a@,+-- backed by a 'G.Vector' @v@ (boxed, unboxed, storable, etc.).+--+-- Use pattern 'X' and operator '^-' for construction:+--+-- >>> (X + 1) + (X^-1 - 1) :: VLaurent Integer+-- 1 * X + 0 + 1 * X^-1+-- >>> (X + 1) * (1 - X^-1) :: ULaurent Int+-- 1 * X + 0 + (-1) * X^-1+--+-- Polynomials are stored normalized, without leading+-- and trailing+-- zero coefficients, so 0 * X + 1 + 0 * X^-1 equals to 1.+--+-- 'Ord' instance does not make much sense mathematically,+-- it is defined only for the sake of 'Data.Set.Set', 'Data.Map.Map', etc.+--+data Laurent (v :: Type -> Type) (a :: Type) = Laurent !Int !(Poly v a)+  deriving (Eq, Ord)++-- | Deconstruct a 'Laurent' polynomial into an offset (largest possible)+-- and a regular polynomial.+--+-- >>> unLaurent (2 * X + 1 :: ULaurent Int)+-- (0,2 * X + 1)+-- >>> unLaurent (1 + 2 * X^-1 :: ULaurent Int)+-- (-1,1 * X + 2)+-- >>> unLaurent (2 * X^2 + X :: ULaurent Int)+-- (1,2 * X + 1)+-- >>> unLaurent (0 :: ULaurent Int)+-- (0,0)+unLaurent :: Laurent v a -> (Int, Poly v a)+unLaurent (Laurent off poly) = (off, poly)++-- | Construct 'Laurent' polynomial from an offset and a regular polynomial.+-- One can imagine it as 'Data.Poly.Semiring.scale', but allowing negative offsets.+--+-- >>> toLaurent 2 (2 * Data.Poly.X + 1) :: ULaurent Int+-- 2 * X^3 + 1 * X^2+-- >>> toLaurent (-2) (2 * Data.Poly.X + 1) :: ULaurent Int+-- 2 * X^-1 + 1 * X^-2+toLaurent+  :: (Eq a, Semiring a, G.Vector v a)+  => Int+  -> Poly v a+  -> Laurent v a+toLaurent off (Poly xs) = go 0+  where+    go k+      | k >= G.length xs+      = Laurent 0 zero+      | G.unsafeIndex xs k == zero+      = go (k + 1)+      | otherwise+      = Laurent (off + k) (Poly (G.unsafeDrop k xs))+{-# INLINE toLaurent #-}++toLaurentNum+  :: (Eq a, Num a, G.Vector v a)+  => Int+  -> Poly v a+  -> Laurent v a+toLaurentNum off (Poly xs) = go 0+  where+    go k+      | k >= G.length xs+      = Laurent 0 0+      | G.unsafeIndex xs k == 0+      = go (k + 1)+      | otherwise+      = Laurent (off + k) (Poly (G.unsafeDrop k xs))+{-# INLINE toLaurentNum #-}++instance NFData (v a) => NFData (Laurent v a) where+  rnf (Laurent off poly) = rnf off `seq` rnf poly++instance (Show a, G.Vector v a) => Show (Laurent v a) where+  showsPrec d (Laurent off poly)+    | G.null (unPoly poly)+      = showString "0"+    | otherwise+      = showParen (d > 0)+      $ foldl (.) id+      $ intersperse (showString " + ")+      $ G.ifoldl (\acc i c -> showCoeff (i + off) c : acc) []+      $ unPoly poly+    where+      -- Negative powers should be displayed without surrounding brackets+      showCoeff 0 c = showsPrec 7 c+      showCoeff 1 c = showsPrec 7 c . showString " * X"+      showCoeff i c = showsPrec 7 c . showString (" * X^" ++ show i)++-- | Laurent polynomials backed by boxed vectors.+type VLaurent = Laurent V.Vector++-- | Laurent polynomials backed by unboxed vectors.+type ULaurent = Laurent U.Vector++-- | Return a 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+leading :: G.Vector v a => Laurent v a -> Maybe (Int, a)+leading (Laurent off poly) = first ((+ off) . fromIntegral) <$> Dense.leading poly++-- | Note that 'abs' = 'id' and 'signum' = 'const' 1.+instance (Eq a, Num a, G.Vector v a) => Num (Laurent v a) where+  Laurent off1 poly1 * Laurent off2 poly2 = toLaurentNum (off1 + off2) (poly1 * poly2)+  Laurent off1 poly1 + Laurent off2 poly2 = case off1 `compare` off2 of+    LT -> toLaurentNum off1 (poly1 + Dense.scale (fromIntegral $ off2 - off1) 1 poly2)+    EQ -> toLaurentNum off1 (poly1 + poly2)+    GT -> toLaurentNum off2 (Dense.scale (fromIntegral $ off1 - off2) 1 poly1 + poly2)+  Laurent off1 poly1 - Laurent off2 poly2 = case off1 `compare` off2 of+    LT -> toLaurentNum off1 (poly1 - Dense.scale (fromIntegral $ off2 - off1) 1 poly2)+    EQ -> toLaurentNum off1 (poly1 - poly2)+    GT -> toLaurentNum off2 (Dense.scale (fromIntegral $ off1 - off2) 1 poly1 - poly2)+  negate (Laurent off poly) = Laurent off (negate poly)+  abs = id+  signum = const 1+  fromInteger n = Laurent 0 (fromInteger n)+  {-# INLINE (+) #-}+  {-# INLINE (-) #-}+  {-# INLINE negate #-}+  {-# INLINE fromInteger #-}+  {-# INLINE (*) #-}++instance (Eq a, Semiring a, G.Vector v a) => Semiring (Laurent v a) where+  zero = Laurent 0 zero+  one  = Laurent 0 one+  Laurent off1 poly1 `times` Laurent off2 poly2 =+    toLaurent (off1 + off2) (poly1 `times` poly2)+  Laurent off1 poly1 `plus` Laurent off2 poly2 = case off1 `compare` off2 of+    LT -> toLaurent off1 (poly1 `plus` Dense.scale' (fromIntegral $ off2 - off1) one poly2)+    EQ -> toLaurent off1 (poly1 `plus` poly2)+    GT -> toLaurent off2 (Dense.scale' (fromIntegral $ off1 - off2) one poly1 `plus` poly2)+  fromNatural n = Laurent 0 (fromNatural n)+  {-# INLINE zero #-}+  {-# INLINE one #-}+  {-# INLINE plus #-}+  {-# INLINE times #-}+  {-# INLINE fromNatural #-}++instance (Eq a, Ring a, G.Vector v a) => Ring (Laurent v a) where+  negate (Laurent off poly) = Laurent off (Semiring.negate poly)++-- | Create a monomial from a power and a coefficient.+monomial :: (Eq a, Semiring a, G.Vector v a) => Int -> a -> Laurent v a+monomial p c+  | c == zero = Laurent 0 zero+  | otherwise = Laurent p (Dense.monomial' 0 c)+{-# INLINE 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 + 0 * X + 3+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)++-- | Evaluate at a given point.+--+-- >>> eval (X^-2 + 1 :: ULaurent Double) 2+-- 1.25+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))+{-# INLINE eval #-}++-- | Substitute another polynomial instead of 'Data.Poly.X'.+--+-- >>> import Data.Poly (UPoly)+-- >>> subst (Data.Poly.X^2 + 1 :: UPoly Int) (X^-1 + 1 :: ULaurent Int)+-- 2 + 2 * X^-1 + 1 * X^-2+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.+--+-- >>> deriv (X^-1 + 3 * X) :: ULaurent Int+-- 3 + 0 * X^-1 + (-1) * X^-2+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.+pattern X :: (Eq a, Semiring a, G.Vector v a) => Laurent v a+pattern X <- (isVar -> True)+  where X = var++var :: forall a v. (Eq a, Semiring a, G.Vector v a) => Laurent v a+var+  | (one :: a) == zero = Laurent 0 zero+  | otherwise          = Laurent 1 one+{-# INLINE var #-}++isVar :: forall v a. (Eq a, Semiring a, G.Vector v a) => Laurent v a -> Bool+isVar (Laurent off (Poly xs))+  | (one :: a) == zero = off == 0 && G.null xs+  | otherwise          = off == 1 && G.length xs == 1 && G.unsafeHead xs == one+{-# INLINE isVar #-}++-- | This operator can be applied only to monomials with unit coefficients,+-- but is instrumental to express Laurent polynomials+-- in mathematical fashion:+--+-- >>> X + 2 + 3 * (X^2)^-1 :: ULaurent Int+-- 1 * X + 2 + 0 * X^-1 + 3 * X^-2+(^-)+  :: (Eq a, Num a, G.Vector v a)+  => Laurent v a+  -> Int+  -> Laurent v a+Laurent off (Poly xs) ^- n+  | G.length xs == 1, G.unsafeHead xs == 1+  = Laurent (off * (-n)) (Poly xs)+  | otherwise+  = throw $ PatternMatchFail "(^-) can be applied only to a monom with unit coefficient"++instance (Eq a, Ring a, GcdDomain a, G.Vector v a) => GcdDomain (Laurent v a) where+  divide (Laurent off1 poly1) (Laurent off2 poly2) =+    toLaurent (off1 - off2) <$> divide poly1 poly2+  {-# INLINE divide #-}++  gcd (Laurent _ poly1) (Laurent _ poly2) =+    toLaurent 0 (gcd poly1 poly2)+  {-# INLINE gcd #-}++  lcm (Laurent _ poly1) (Laurent _ poly2) =+    toLaurent 0 (lcm poly1 poly2)+  {-# INLINE lcm #-}++  coprime (Laurent _ poly1) (Laurent _ poly2) =+    coprime poly1 poly2+  {-# INLINE coprime #-}

@@ -0,0 +1,553 @@+-- |+-- Module:      Data.Poly.Internal.Multi+-- Copyright:   (c) 2020 Andrew Lelechenko+-- Licence:     BSD3+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+--++{-# LANGUAGE DataKinds                  #-}+{-# LANGUAGE FlexibleContexts           #-}+{-# LANGUAGE FlexibleInstances          #-}+{-# LANGUAGE GADTs                      #-}+{-# LANGUAGE GeneralizedNewtypeDeriving #-}+{-# LANGUAGE LambdaCase                 #-}+{-# LANGUAGE PatternSynonyms            #-}+{-# LANGUAGE PolyKinds                  #-}+{-# LANGUAGE ScopedTypeVariables        #-}+{-# LANGUAGE StandaloneDeriving         #-}+{-# LANGUAGE TypeFamilies               #-}+{-# LANGUAGE TypeOperators              #-}+{-# LANGUAGE UndecidableInstances       #-}+{-# LANGUAGE ViewPatterns               #-}++{-# OPTIONS_GHC -fno-warn-orphans #-}++module Data.Poly.Internal.Multi+  ( MultiPoly(..)+  , VMultiPoly+  , UMultiPoly+  , toMultiPoly+  , toMultiPoly'+  , leading+  , monomial+  , monomial'+  , scale+  , scale'+  , pattern X+  , pattern Y+  , pattern Z+  , pattern X'+  , pattern Y'+  , pattern Z'+  , eval+  , eval'+  , subst+  , subst'+  , substitute+  , substitute'+  , deriv+  , deriv'+  , integral+  , integral'+  -- * Univariate polynomials+  , Poly+  , VPoly+  , UPoly+  , unPoly+  -- * Conversions+  , segregate+  , unsegregate+  ) where++import Prelude hiding (quot, gcd)+import Control.Arrow+import Control.DeepSeq+import Data.Euclidean (Field, quot)+import Data.Finite+import Data.Kind+import Data.List (intersperse)+import Data.Semiring (Semiring(..), Ring())+import qualified Data.Semiring as Semiring+import qualified Data.Vector as V+import qualified Data.Vector.Generic as G+import qualified Data.Vector.Generic.Sized as SG+import qualified Data.Vector.Unboxed as U+import qualified Data.Vector.Unboxed.Sized as SU+import qualified Data.Vector.Sized as SV+import GHC.Exts (IsList(..))+import GHC.TypeNats (KnownNat, Nat, type (+), type (<=))++import Data.Poly.Internal.Multi.Core++-- | 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',+-- 'Data.Poly.Multi.Y' and+-- 'Data.Poly.Multi.Z' for construction:+--+-- >>> :set -XDataKinds+-- >>> (X + 1) + (Y - 1) + Z :: VMultiPoly 3 Integer+-- 1 * X + 1 * Y + 1 * Z+-- >>> (X + 1) * (Y - 1) :: UMultiPoly 2 Int+-- 1 * X * Y + (-1) * X + 1 * Y + (-1)+--+-- 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,+-- it is defined only for the sake of 'Data.Set.Set', 'Data.Map.Map', etc.+--+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.+  }++deriving instance Eq     (v (SU.Vector n Word, a)) => Eq     (MultiPoly v n a)+deriving instance Ord    (v (SU.Vector n Word, a)) => Ord    (MultiPoly v n a)+deriving instance NFData (v (SU.Vector n Word, a)) => NFData (MultiPoly v n a)++-- | Multivariate polynomials backed by boxed vectors.+type VMultiPoly (n :: Nat) (a :: Type) = MultiPoly V.Vector n a++-- | Multivariate polynomials backed by unboxed vectors.+type UMultiPoly (n :: Nat) (a :: Type) = MultiPoly U.Vector n a++-- | Sparse univariate polynomials with coefficients from @a@,+-- backed by a 'G.Vector' @v@ (boxed, unboxed, storable, etc.).+--+-- Use pattern 'Data.Poly.Multi.X' for construction:+--+-- >>> (X + 1) + (X - 1) :: VPoly Integer+-- 2 * X+-- >>> (X + 1) * (X - 1) :: UPoly Int+-- 1 * X^2 + (-1)+--+-- 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,+-- it is defined only for the sake of 'Data.Set.Set', 'Data.Map.Map', etc.+--+type Poly (v :: Type -> Type) (a :: Type) = MultiPoly v 1 a++-- | Polynomials backed by boxed vectors.+type VPoly (a :: Type) = Poly V.Vector a++-- | Polynomials backed by unboxed vectors.+type UPoly (a :: Type) = Poly U.Vector a++-- | Convert 'Poly' to a vector of coefficients.+unPoly+  :: (G.Vector v (Word, a), G.Vector v (SU.Vector 1 Word, a))+  => Poly v a+  -> v (Word, a)+unPoly = G.map (first SU.head) . unMultiPoly++instance (Eq a, Semiring a, G.Vector v (SU.Vector n Word, a)) => IsList (MultiPoly v n a) where+  type Item (MultiPoly v n a) = (SU.Vector n Word, a)+  fromList = toMultiPoly' . G.fromList+  fromListN = (toMultiPoly' .) . G.fromListN+  toList = G.toList . unMultiPoly++instance (Show a, KnownNat n, G.Vector v (SU.Vector n Word, a)) => Show (MultiPoly v n a) where+  showsPrec d (MultiPoly xs)+    | G.null xs+      = showString "0"+    | otherwise+      = showParen (d > 0)+      $ foldl (.) id+      $ intersperse (showString " + ")+      $ G.foldl (\acc (is, c) -> showCoeff is c : acc) [] xs+    where+      showCoeff is c+        = showsPrec 7 c . foldl (.) id+        ( map ((showString " * " .) . uncurry showPower)+        $ filter ((/= 0) . fst)+        $ zip (SU.toList is) (finites :: [Finite n]))++      -- Powers are guaranteed to be non-negative+      showPower :: Word -> Finite n -> String -> String+      showPower 1 n = showString (showVar n)+      showPower i n = showString (showVar n) . showString ("^" ++ show i)++      showVar :: Finite n -> String+      showVar = \case+        0 -> "X"+        1 -> "Y"+        2 -> "Z"+        k -> "X" ++ show k++-- | Make 'MultiPoly' from a list of (powers, coefficient) pairs.+--+-- >>> :set -XOverloadedLists -XDataKinds+-- >>> import Data.Vector.Generic.Sized (fromTuple)+-- >>> toMultiPoly [(fromTuple (0,0),1),(fromTuple (0,1),2),(fromTuple (1,0),3)] :: VMultiPoly 2 Integer+-- 3 * X + 2 * Y + 1+-- >>> toMultiPoly [(fromTuple (0,0),0),(fromTuple (0,1),0),(fromTuple (1,0),0)] :: UMultiPoly 2 Int+-- 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) (+)++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++-- | 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+  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 #-}+  {-# INLINE fromInteger #-}+  {-# INLINE (*) #-}++instance (Eq a, Semiring a, KnownNat n, G.Vector v (SU.Vector n Word, a)) => Semiring (MultiPoly v n a) where+  zero = MultiPoly G.empty+  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+  {-# INLINE zero #-}+  {-# INLINE one #-}+  {-# INLINE plus #-}+  {-# INLINE times #-}++  fromNatural n = if n' == zero then zero else MultiPoly $ G.singleton (0, n')+    where+      n' :: a+      n' = fromNatural n+  {-# INLINE fromNatural #-}++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.+--+-- >>> import Data.Poly.Sparse (UPoly)+-- >>> leading ((2 * X + 1) * (2 * X^2 - 1) :: UPoly Int)+-- Just (3,4)+-- >>> leading (0 :: UPoly Int)+-- Nothing+leading :: G.Vector v (SU.Vector 1 Word, a) => Poly v a -> Maybe (Word, a)+leading (MultiPoly v)+  | G.null v  = Nothing+  | otherwise = Just $ first SU.head $ G.last v++-- | Multiply a polynomial by a monomial, expressed as powers and a coefficient.+--+-- >>> :set -XDataKinds+-- >>> 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+scale+  :: (Eq a, Num a, KnownNat n, G.Vector v (SU.Vector n Word, a))+  => SU.Vector n Word+  -> a+  -> MultiPoly v n a+  -> MultiPoly v n a+scale yp yc = MultiPoly . scaleInternal (/= 0) (*) yp yc . unMultiPoly++scale'+  :: (Eq a, Semiring a, KnownNat n, G.Vector v (SU.Vector n Word, a))+  => SU.Vector n Word+  -> a+  -> MultiPoly v n a+  -> MultiPoly v n a+scale' yp yc = MultiPoly . scaleInternal (/= zero) times yp yc . unMultiPoly++-- | Create a monomial from powers and a coefficient.+monomial+  :: (Eq a, Num a, G.Vector v (SU.Vector n Word, a))+  => SU.Vector n Word+  -> a+  -> MultiPoly v n a+monomial p c+  | c == 0    = MultiPoly G.empty+  | otherwise = MultiPoly $ G.singleton (p, c)++monomial'+  :: (Eq a, Semiring a, G.Vector v (SU.Vector n Word, a))+  => SU.Vector n Word+  -> a+  -> MultiPoly v n a+monomial' p c+  | c == zero = MultiPoly G.empty+  | otherwise = MultiPoly $ G.singleton (p, c)++-- | Evaluate 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+eval+  :: (Num a, G.Vector v (SU.Vector n Word, a), G.Vector u a)+  => MultiPoly v n a+  -> SG.Vector u n a+  -> a+eval = substitute (*)+{-# INLINE eval #-}++eval'+  :: (Semiring a, G.Vector v (SU.Vector n Word, a), G.Vector u a)+  => MultiPoly v n a+  -> SG.Vector u n a+  -> a+eval' = substitute' times+{-# INLINE eval' #-}++-- | Substitute another polynomials instead of 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+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+  -> SV.Vector n (MultiPoly w m a)+  -> MultiPoly w m a+subst = substitute (scale 0)+{-# INLINE subst #-}++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+  -> SV.Vector n (MultiPoly w m a)+  -> MultiPoly w m a+subst' = substitute' (scale' 0)+{-# INLINE subst' #-}++substitute+  :: forall v u n a b.+     (G.Vector v (SU.Vector n Word, a), G.Vector u b, Num b)+  => (a -> b -> b)+  -> MultiPoly v n a+  -> SG.Vector u n b+  -> b+substitute f (MultiPoly cs) xs = G.foldl' go 0 cs+  where+    go :: b -> (SU.Vector n Word, a) -> b+    go acc (ps, c) = acc + f c (doMonom ps)++    doMonom :: SU.Vector n Word -> b+    doMonom = SU.ifoldl' (\acc i p -> acc * ((xs `SG.index` i) ^ p)) 1+{-# INLINE substitute #-}++substitute'+  :: forall v u n a b.+     (G.Vector v (SU.Vector n Word, a), G.Vector u b, Semiring b)+  => (a -> b -> b)+  -> MultiPoly v n a+  -> SG.Vector u n b+  -> b+substitute' f (MultiPoly cs) xs = G.foldl' go zero cs+  where+    go :: b -> (SU.Vector n Word, a) -> b+    go acc (ps, c) = acc `plus` f c (doMonom ps)++    doMonom :: SU.Vector n Word -> b+    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.+--+-- >>> :set -XDataKinds+-- >>> deriv 0 (X^3 + 3 * Y) :: UMultiPoly 2 Int+-- 3 * X^2+-- >>> deriv 1 (X^3 + 3 * Y) :: UMultiPoly 2 Int+-- 3+deriv+  :: (Eq a, Num a, G.Vector v (SU.Vector n Word, a))+  => Finite n+  -> MultiPoly v n a+  -> MultiPoly v n a+deriv i (MultiPoly xs) = MultiPoly $ derivPoly+  (/= 0)+  (\ps -> ps SU.// [(i, ps `SU.index` i - 1)])+  (\ps c -> fromIntegral (ps `SU.index` i) * c)+  xs+{-# INLINE deriv #-}++deriv'+  :: (Eq a, Semiring a, G.Vector v (SU.Vector n Word, a))+  => Finite n+  -> MultiPoly v n a+  -> MultiPoly v n a+deriv' i (MultiPoly xs) = MultiPoly $ derivPoly+  (/= zero)+  (\ps -> ps SU.// [(i, ps `SU.index` i - 1)])+  (\ps c -> fromNatural (fromIntegral (ps `SU.index` i)) `times` c)+  xs+{-# INLINE deriv' #-}++-- | Compute an indefinite integral of a polynomial+-- by the /i/-th variable,+-- setting 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+integral+  :: (Fractional a, G.Vector v (SU.Vector n Word, a))+  => Finite n+  -> MultiPoly v n a+  -> MultiPoly v n a+integral i (MultiPoly xs)+  = MultiPoly+  $ G.map (\(ps, c) -> let p = ps `SU.index` i in+    (ps SU.// [(i, p + 1)], c / fromIntegral (p + 1))) xs+{-# INLINE integral #-}++integral'+  :: (Field a, G.Vector v (SU.Vector n Word, a))+  => Finite n+  -> MultiPoly v n a+  -> MultiPoly v n a+integral' i (MultiPoly xs)+  = MultiPoly+  $ G.map (\(ps, c) -> let p = ps `SU.index` i in+    (ps SU.// [(i, p + 1)], c `quot` Semiring.fromIntegral (p + 1))) xs+{-# INLINE integral' #-}++-- | Create a polynomial equal to the first variable.+pattern X+  :: (Eq a, Num a, KnownNat n, 1 <= n, G.Vector v (SU.Vector n Word, a))+  => MultiPoly v n a+pattern X <- (isVar 0 -> True)+  where X = var 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' <- (isVar' 0 -> True)+  where X' = var' 0++-- | Create a polynomial equal to the second variable.+pattern Y+  :: (Eq a, Num a, KnownNat n, 2 <= n, G.Vector v (SU.Vector n Word, a))+  => MultiPoly v n a+pattern Y <- (isVar 1 -> True)+  where Y = var 1++pattern Y'+  :: (Eq a, Semiring a, KnownNat n, 2 <= n, G.Vector v (SU.Vector n Word, a))+  => MultiPoly v n a+pattern Y' <- (isVar' 1 -> True)+  where Y' = var' 1++-- | Create a polynomial equal to the third variable.+pattern Z+  :: (Eq a, Num a, KnownNat n, 3 <= n, G.Vector v (SU.Vector n Word, a))+  => MultiPoly v n a+pattern Z <- (isVar 2 -> True)+  where Z = var 2++pattern Z'+  :: (Eq a, Semiring a, KnownNat n, 3 <= n, G.Vector v (SU.Vector n Word, a))+  => MultiPoly v n a+pattern Z' <- (isVar' 2 -> True)+  where Z' = var' 2++var+  :: forall v n a.+     (Eq a, Num a, KnownNat n, G.Vector v (SU.Vector n Word, a))+  => Finite n+  -> MultiPoly v n a+var i+  | (1 :: a) == 0 = MultiPoly G.empty+  | otherwise     = MultiPoly $ G.singleton+    (SU.generate (\j -> if i == j then 1 else 0), 1)+{-# INLINE var #-}++var'+  :: forall v n a.+     (Eq a, Semiring a, KnownNat n, G.Vector v (SU.Vector n Word, a))+  => Finite n+  -> MultiPoly v n a+var' i+  | (one :: a) == zero = MultiPoly G.empty+  | otherwise          = MultiPoly $ G.singleton+    (SU.generate (\j -> if i == j then 1 else 0), one)+{-# INLINE var' #-}++isVar+  :: forall v n a.+     (Eq a, Num a, KnownNat n, G.Vector v (SU.Vector n Word, a))+  => Finite n+  -> MultiPoly v n a+  -> Bool+isVar i (MultiPoly xs)+  | (1 :: a) == 0 = G.null xs+  | otherwise     = G.length xs == 1 && G.unsafeHead xs == (SU.generate (\j -> if i == j then 1 else 0), 1)+{-# INLINE isVar #-}++isVar'+  :: forall v n a.+     (Eq a, Semiring a, KnownNat n, G.Vector v (SU.Vector n Word, a))+  => Finite n+  -> MultiPoly v n a+  -> Bool+isVar' i (MultiPoly xs)+  | (one :: a) == zero = G.null xs+  | otherwise          = G.length xs == 1 && G.unsafeHead xs == (SU.generate (\j -> if i == j then 1 else 0), one)+{-# INLINE isVar' #-}++-------------------------------------------------------------------------------++groupOn :: (G.Vector v a, Eq b) => (a -> b) -> v a -> [v a]+groupOn f = go+  where+    go xs+      | G.null xs = []+      | otherwise = case mk of+        Nothing -> [xs]+        Just k  -> G.unsafeTake (k + 1) xs : go (G.unsafeDrop (k + 1) xs)+        where+          fy = f (G.unsafeHead xs)+          mk = G.findIndex ((/= fy) . f) (G.unsafeTail xs)++-- | Interpret a multivariate polynomial over 1+/m/ variables+-- as a univariate polynomial, whose coefficients are+-- multivariate polynomials over the last /m/ variables.+segregate+  :: (G.Vector v (SU.Vector (1 + m) Word, a), G.Vector v (SU.Vector m Word, a))+  => MultiPoly v (1 + m) a+  -> VPoly (MultiPoly v m a)+segregate+  = MultiPoly+  . G.fromList+  . map (\vs -> (SU.take (fst (G.unsafeHead vs)), MultiPoly $ G.map (first SU.tail) vs))+  . groupOn (SU.head . fst)+  . unMultiPoly++-- | Interpret a univariate polynomials, whose coefficients are+-- multivariate polynomials over the first /m/ variables,+-- as a multivariate polynomial over 1+/m/ variables.+unsegregate+  :: (G.Vector v (SU.Vector (1 + m) Word, a), G.Vector v (SU.Vector m Word, a))+  => VPoly (MultiPoly v m a)+  -> MultiPoly v (1 + m) a+unsegregate+  = MultiPoly+  . G.concat+  . G.toList+  . G.map (\(v, MultiPoly vs) -> G.map (first (v SU.++)) vs)+  . unMultiPoly

@@ -0,0 +1,311 @@+-- |+-- Module:      Data.Poly.Internal.Multi.Core+-- Copyright:   (c) 2019 Andrew Lelechenko+-- Licence:     BSD3+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+--+-- Sparse polynomials of one variable.+--++{-# LANGUAGE DataKinds                  #-}+{-# LANGUAGE FlexibleContexts           #-}+{-# LANGUAGE ScopedTypeVariables        #-}+{-# LANGUAGE TypeFamilies               #-}+{-# LANGUAGE UndecidableInstances       #-}++module Data.Poly.Internal.Multi.Core+  ( normalize+  , plusPoly+  , minusPoly+  , convolution+  , scaleInternal+  , derivPoly+  ) where++import Control.Monad+import Control.Monad.Primitive+import Control.Monad.ST+import Data.Bits+import Data.Ord+import qualified Data.Vector.Algorithms.Tim as Tim+import qualified Data.Vector.Generic as G+import qualified Data.Vector.Generic.Mutable as MG+import qualified Data.Vector.Unboxed as U++normalize+  :: (G.Vector v (t, a), Ord t)+  => (a -> Bool)+  -> (a -> a -> a)+  -> v (t, a)+  -> v (t, a)+normalize p add vs+  | G.null vs = vs+  | otherwise = runST $ do+    ws <- G.thaw vs+    l' <- normalizeM p add ws+    G.unsafeFreeze $ MG.unsafeSlice 0 l' ws++normalizeM+  :: (PrimMonad m, G.Vector v (t, a), Ord t)+  => (a -> Bool)+  -> (a -> a -> a)+  -> G.Mutable v (PrimState m) (t, a)+  -> m Int+normalizeM p add ws = do+    let l = MG.length ws+    let go i j acc@(accP, accC)+          | j >= l =+            if p accC+              then do+                MG.write ws i acc+                pure $ i + 1+              else pure i+          | otherwise = do+            v@(vp, vc) <- MG.unsafeRead ws j+            if vp == accP+              then go i (j + 1) (accP, accC `add` vc)+              else if p accC+                then do+                  MG.write ws i acc+                  go (i + 1) (j + 1) v+                else go i (j + 1) v+    Tim.sortBy (comparing fst) ws+    wsHead <- MG.unsafeRead ws 0+    go 0 1 wsHead++plusPoly+  :: (G.Vector v (t, a), Ord t)+  => (a -> Bool)+  -> (a -> a -> a)+  -> v (t, a)+  -> v (t, a)+  -> v (t, a)+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)+  => (a -> Bool)+  -> (a -> a -> a)+  -> v (t, a)+  -> v (t, a)+  -> G.Mutable v (PrimState m) (t, a)+  -> m Int+plusPolyM p add xs ys zs = go 0 0 0+  where+    lenXs = G.length xs+    lenYs = G.length ys++    go ix iy iz+      | ix == lenXs, iy == lenYs = pure iz+      | ix == lenXs = do+        G.unsafeCopy+          (MG.unsafeSlice iz (lenYs - iy) zs)+          (G.unsafeSlice iy (lenYs - iy) ys)+        pure $ iz + lenYs - iy+      | iy == lenYs = do+        G.unsafeCopy+          (MG.unsafeSlice iz (lenXs - ix) zs)+          (G.unsafeSlice ix (lenXs - ix) xs)+        pure $ iz + lenXs - ix+      | (xp, xc) <- G.unsafeIndex xs ix+      , (yp, yc) <- G.unsafeIndex ys iy+      = case xp `compare` yp of+        LT -> do+          MG.unsafeWrite zs iz (xp, xc)+          go (ix + 1) iy (iz + 1)+        EQ -> do+          let zc = xc `add` yc+          if p zc then do+            MG.unsafeWrite zs iz (xp, zc)+            go (ix + 1) (iy + 1) (iz + 1)+          else+            go (ix + 1) (iy + 1) iz+        GT -> do+          MG.unsafeWrite zs iz (yp, yc)+          go ix (iy + 1) (iz + 1)+{-# INLINABLE plusPolyM #-}++minusPoly+  :: (G.Vector v (t, a), Ord t)+  => (a -> Bool)+  -> (a -> a)+  -> (a -> a -> a)+  -> v (t, a)+  -> v (t, a)+  -> v (t, a)+minusPoly p neg sub xs ys = runST $ do+  zs <- MG.unsafeNew (lenXs + lenYs)+  let go ix iy iz+        | ix == lenXs, iy == lenYs = pure iz+        | ix == lenXs = do+          forM_ [iy .. lenYs - 1] $ \i ->+            MG.unsafeWrite zs (iz + i - iy)+              (fmap neg (G.unsafeIndex ys i))+          pure $ iz + lenYs - iy+        | iy == lenYs = do+          G.unsafeCopy+            (MG.unsafeSlice iz (lenXs - ix) zs)+            (G.unsafeSlice ix (lenXs - ix) xs)+          pure $ iz + lenXs - ix+        | (xp, xc) <- G.unsafeIndex xs ix+        , (yp, yc) <- G.unsafeIndex ys iy+        = case xp `compare` yp of+          LT -> do+            MG.unsafeWrite zs iz (xp, xc)+            go (ix + 1) iy (iz + 1)+          EQ -> do+            let zc = xc `sub` yc+            if p zc then do+              MG.unsafeWrite zs iz (xp, zc)+              go (ix + 1) (iy + 1) (iz + 1)+            else+              go (ix + 1) (iy + 1) iz+          GT -> do+            MG.unsafeWrite zs iz (yp, neg yc)+            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)+  => (a -> Bool)+  -> (a -> a -> a)+  -> v (t, a)+  -> (t, a)+  -> G.Mutable v (PrimState m) (t, a)+  -> m Int+scaleM p mul xs (yp, yc) zs = go 0 0+  where+    lenXs = G.length xs++    go ix iz+      | ix == lenXs = pure iz+      | (xp, xc) <- G.unsafeIndex xs ix+      = do+        let zc = xc `mul` yc+        if p zc then do+          MG.unsafeWrite zs iz (xp + yp, zc)+          go (ix + 1) (iz + 1)+        else+          go (ix + 1) iz+{-# INLINABLE scaleM #-}++scaleInternal+  :: (G.Vector v (t, a), Num t)+  => (a -> Bool)+  -> (a -> a -> a)+  -> t+  -> a+  -> v (t, a)+  -> v (t, a)+scaleInternal p mul yp yc xs = runST $ do+  zs <- MG.unsafeNew (G.length xs)+  len <- scaleM p (flip mul) xs (yp, yc) zs+  G.unsafeFreeze $ MG.unsafeSlice 0 len zs+{-# INLINABLE scaleInternal #-}++convolution+  :: forall v t a.+     (G.Vector v (t, a), Ord t, Num t)+  => (a -> Bool)+  -> (a -> a -> a)+  -> (a -> a -> a)+  -> 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+  where+    go :: (a -> a -> a) -> v (t, a) -> v (t, a) -> v (t, a)+    go mul long short = runST $ do+      let lenLong   = G.length long+          lenShort  = G.length short+          lenBuffer = lenLong * lenShort+      slices <- MG.unsafeNew lenShort+      buffer <- MG.unsafeNew lenBuffer++      forM_ [0 .. lenShort - 1] $ \iShort -> do+        let (pShort, cShort) = G.unsafeIndex short iShort+            from = iShort * lenLong+            bufferSlice = MG.unsafeSlice from lenLong buffer+        len <- scaleM p mul long (pShort, cShort) bufferSlice+        MG.unsafeWrite slices iShort (from, len)++      slices' <- G.unsafeFreeze slices+      buffer' <- G.unsafeFreeze buffer+      bufferNew <- MG.unsafeNew lenBuffer+      gogo slices' buffer' bufferNew++    gogo+      :: PrimMonad m+      => U.Vector (Int, Int)+      -> v (t, a)+      -> G.Mutable v (PrimState m) (t, a)+      -> m (v (t, a))+    gogo slices buffer bufferNew+      | G.length slices == 0+      = pure G.empty+      | G.length slices == 1+      , (from, len) <- G.unsafeIndex slices 0+      = pure $ G.unsafeSlice from len buffer+      | otherwise = do+        let nSlices = G.length slices+        slicesNew <- MG.unsafeNew ((nSlices + 1) `shiftR` 1)+        forM_ [0 .. (nSlices - 2) `shiftR` 1] $ \i -> do+          let (from1, len1) = G.unsafeIndex slices (2 * i)+              (from2, len2) = G.unsafeIndex slices (2 * i + 1)+              slice1 = G.unsafeSlice from1 len1 buffer+              slice2 = G.unsafeSlice from2 len2 buffer+              slice3 = MG.unsafeSlice from1 (len1 + len2) bufferNew+          len3 <- plusPolyM p add slice1 slice2 slice3+          MG.unsafeWrite slicesNew i (from1, len3)++        when (odd nSlices) $ do+          let (from, len) = G.unsafeIndex slices (nSlices - 1)+              slice1 = G.unsafeSlice from len buffer+              slice3 = MG.unsafeSlice from len bufferNew+          G.unsafeCopy slice3 slice1+          MG.unsafeWrite slicesNew (nSlices `shiftR` 1) (from, len)++        slicesNew' <- G.unsafeFreeze slicesNew+        buffer'    <- G.unsafeThaw   buffer+        bufferNew' <- G.unsafeFreeze bufferNew+        gogo slicesNew' bufferNew' buffer'+{-# INLINABLE convolution #-}++derivPoly+  :: (G.Vector v (t, a))+  => (a -> Bool)   -- ^ is coefficient non-zero?+  -> (t -> t)      -- ^ how to modify powers?+  -> (t -> a -> a) -- ^ how to modify coefficient?+  -> v (t, a)+  -> v (t, a)+derivPoly p dec mul xs+  | G.null xs = G.empty+  | otherwise = runST $ do+    let lenXs = G.length xs+    zs <- MG.unsafeNew lenXs+    let go ix iz+          | ix == lenXs = pure iz+          | (xp, xc) <- G.unsafeIndex xs ix+          = do+            let zc = xp `mul` xc+            if p zc then do+              MG.unsafeWrite zs iz (dec xp, zc)+              go (ix + 1) (iz + 1)+            else+              go (ix + 1) iz+    lenZs <- go 0 0+    G.unsafeFreeze $ MG.unsafeSlice 0 lenZs zs+{-# INLINABLE derivPoly #-}

@@ -0,0 +1,73 @@+-- |+-- Module:      Data.Poly.Internal.Multi.Field+-- Copyright:   (c) 2019 Andrew Lelechenko+-- Licence:     BSD3+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+--+-- Euclidean for Field underlying+--++{-# LANGUAGE ConstraintKinds            #-}+{-# LANGUAGE DataKinds                  #-}+{-# LANGUAGE FlexibleContexts           #-}+{-# LANGUAGE FlexibleInstances          #-}+{-# LANGUAGE ScopedTypeVariables        #-}+{-# LANGUAGE TypeFamilies               #-}+{-# LANGUAGE UndecidableInstances       #-}++{-# OPTIONS_GHC -fno-warn-orphans #-}++module Data.Poly.Internal.Multi.Field+  ( quotRemFractional+  ) where++import Prelude hiding (quotRem, quot, rem, div, gcd)+import Control.Arrow+import Control.Exception+import Data.Euclidean (Euclidean(..), Field)+import Data.Semiring (Semiring(..), minus)+import qualified Data.Vector.Generic as G+import qualified Data.Vector.Unboxed.Sized as SU++import Data.Poly.Internal.Multi+import Data.Poly.Internal.Multi.GcdDomain ()++-- | Note that 'degree' 0 = 0.+instance (Eq a, Field a, G.Vector v (SU.Vector 1 Word, a)) => Euclidean (Poly v a) where+  degree (MultiPoly xs)+    | G.null xs = 0+    | otherwise = fromIntegral (SU.head (fst (G.unsafeLast xs)))++  quotRem = quotientRemainder zero plus minus times quot++-- | Polynomial division with remainder.+--+-- >>> quotRemFractional (X^3 + 2) (X^2 - 1 :: UPoly Double)+-- (1.0 * X,1.0 * X + 2.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 #-}++quotientRemainder+  :: G.Vector v (SU.Vector 1 Word, a)+  => Poly v a                           -- ^ zero+  -> (Poly v a -> Poly v a -> Poly v a) -- ^ add+  -> (Poly v a -> Poly v a -> Poly v a) -- ^ subtract+  -> (Poly v a -> Poly v a -> Poly v a) -- ^ multiply+  -> (a -> a -> a)                      -- ^ divide+  -> Poly v a                           -- ^ dividend+  -> Poly v a                           -- ^ divisor+  -> (Poly v a, Poly v a)+quotientRemainder zer add sub mul div ts ys = case leading ys of+  Nothing -> throw DivideByZero+  Just (yp, yc) -> go ts+    where+      go xs = case leading xs of+        Nothing -> (zer, zer)+        Just (xp, xc) -> case xp `compare` yp of+          LT -> (zer, xs)+          EQ -> (zs, xs')+          GT -> first (`add` zs) $ go xs'+          where+            zs = MultiPoly $ G.singleton (SU.singleton (xp - yp), xc `div` yc)+            xs' = xs `sub` (zs `mul` ys)

@@ -0,0 +1,179 @@+-- |+-- Module:      Data.Poly.Internal.Multi.GcdDomain+-- Copyright:   (c) 2019 Andrew Lelechenko+-- Licence:     BSD3+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+--+-- GcdDomain for GcdDomain underlying+--++{-# LANGUAGE CPP                        #-}+{-# LANGUAGE DataKinds                  #-}+{-# LANGUAGE FlexibleContexts           #-}+{-# LANGUAGE FlexibleInstances          #-}+{-# LANGUAGE GADTs                      #-}+{-# 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+  () where++import Prelude hiding (gcd, lcm, (^))+import Control.Exception+import Data.Euclidean+import Data.Maybe+import Data.Proxy+import Data.Semiring (Semiring(..), Ring(), minus)+import Data.Type.Equality+import qualified Data.Vector.Generic as G+import qualified Data.Vector.Unboxed.Sized as SU+import GHC.TypeNats (KnownNat, type (+), SomeNat(..), natVal, sameNat, someNatVal)+import Unsafe.Coerce++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+    | G.length (unMultiPoly ys) == 1 = divideSingleton xs (G.unsafeHead (unMultiPoly ys))+    | otherwise = divide1 xs ys++  gcd xs ys+    | G.null (unMultiPoly xs) = ys+    | G.null (unMultiPoly ys) = xs+    | G.length (unMultiPoly xs) == 1 = gcdSingleton (G.unsafeHead (unMultiPoly xs)) ys+    | G.length (unMultiPoly ys) == 1 = gcdSingleton (G.unsafeHead (unMultiPoly ys)) xs+    | otherwise = gcd1 xs ys++  lcm xs ys+    | G.null (unMultiPoly xs) || G.null (unMultiPoly ys) = zero+    | otherwise = (xs `divide'` gcd xs ys) `times` ys++  coprime x y = isJust (one `divide` gcd x y)++data IsSucc n where+  IsSucc :: KnownNat m => n :~: 1 + m -> IsSucc n++-- | This is unsafe when n ~ 0.+isSucc :: forall n. KnownNat n => IsSucc n+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))+    -- Polynomials of zero variables are necessarily constants,+    -- so they have been dealt with above.+    | Just Refl <- sameNat (Proxy :: Proxy n) (Proxy :: Proxy 1)+    = divide1 xs ys+    | otherwise = case isSucc :: IsSucc n of+      IsSucc Refl -> unsegregate <$> segregate xs `divide` segregate ys+  gcd xs ys+    | G.null (unMultiPoly xs) = ys+    | G.null (unMultiPoly ys) = xs+    | G.length (unMultiPoly xs) == 1 = gcdSingleton (G.unsafeHead (unMultiPoly xs)) ys+    | G.length (unMultiPoly ys) == 1 = gcdSingleton (G.unsafeHead (unMultiPoly ys)) xs+    -- Polynomials of zero variables are necessarily constants,+    -- so they have been dealt with above.+    | Just Refl <- sameNat (Proxy :: Proxy n) (Proxy :: Proxy 1)+    = gcd1 xs ys+    | otherwise = case isSucc :: IsSucc n of+      IsSucc Refl -> unsegregate $ segregate xs `gcd` segregate ys++divideSingleton+  :: (GcdDomain a, G.Vector v (SU.Vector n Word, a))+  => MultiPoly v n a+  -> (SU.Vector n Word, a)+  -> Maybe (MultiPoly v n a)+divideSingleton (MultiPoly pcs) (p, c) = MultiPoly <$> G.mapM divideMonomial pcs+  where+    divideMonomial (p', c')+      | SU.and (SU.zipWith (>=) p' p)+      , Just c'' <- c' `divide` c+      = Just (SU.zipWith (-) p' p, c'')+      | otherwise+      = Nothing++gcdSingleton+  :: (Eq a, GcdDomain a, G.Vector v (SU.Vector n Word, a))+  => (SU.Vector n Word, a)+  -> MultiPoly v n a+  -> MultiPoly v n a+gcdSingleton pc (MultiPoly pcs) = uncurry monomial' $+  G.foldl' (\(accP, accC) (p, c) -> (SU.zipWith min accP p, gcd accC c)) pc pcs++divide1+  :: (Eq a, GcdDomain a, Ring a, G.Vector v (SU.Vector 1 Word, a))+  => Poly v a+  -> Poly v a+  -> Maybe (Poly v a)+divide1 xs ys = case leading ys of+  Nothing -> throw DivideByZero+  Just (yp, yc) -> case leading xs of+    Nothing -> Just xs+    Just (xp, xc)+      | xp < yp -> Nothing+      | otherwise -> do+        zc <- divide xc yc+        let z = MultiPoly $ G.singleton (SU.singleton (xp - yp), zc)+        rest <- divide1 (xs `minus` z `times` ys) ys+        pure $ rest `plus` z++gcd1+  :: (Eq a, GcdDomain a, Ring a, G.Vector v (SU.Vector 1 Word, a))+  => Poly v a+  -> Poly v a+  -> Poly v a+gcd1 x@(MultiPoly xs) y@(MultiPoly ys) =+  times xy (divide1' z (monomial' 0 (content zs)))+  where+    z@(MultiPoly zs) = gcdHelper x y+    xy = monomial' 0 (gcd (content xs) (content ys))+    divide1' = (fromMaybe (error "gcd: violated internal invariant") .) . divide1++content :: (GcdDomain a, G.Vector v (t, a)) => v (t, a) -> a+content = G.foldl' (\acc (_, t) -> gcd acc t) zero++gcdHelper+  :: (Eq a, Ring a, GcdDomain a, G.Vector v (SU.Vector 1 Word, a))+  => Poly v a+  -> Poly v a+  -> Poly v a+gcdHelper xs ys = case (leading xs, leading ys) of+  (Nothing, _) -> ys+  (_, Nothing) -> xs+  (Just (xp, xc), Just (yp, yc))+    | yp <= xp+    , Just xy <- xc `divide` yc+    -> gcdHelper ys (xs `minus` ys `times` monomial' (SU.singleton (xp - yp)) xy)+    | xp <= yp+    , Just yx <- yc `divide` xc+    -> gcdHelper xs (ys `minus` xs `times` monomial' (SU.singleton (yp - xp)) yx)+    | yp <= xp+    -> gcdHelper ys (xs `times` monomial' 0 gx `minus` ys `times` monomial' (SU.singleton (xp - yp)) gy)+    | otherwise+    -> gcdHelper xs (ys `times` monomial' 0 gy `minus` xs `times` monomial' (SU.singleton (yp - xp)) gx)+    where+      g = lcm xc yc+      gx = divide' g xc+      gy = divide' g yc++divide' :: GcdDomain a => a -> a -> a+divide' = (fromMaybe (error "gcd: violated internal invariant") .) . divide

@@ -0,0 +1,507 @@+-- |+-- Module:      Data.Poly.Internal.Multi.Laurent+-- Copyright:   (c) 2020 Andrew Lelechenko+-- Licence:     BSD3+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+--+-- Sparse multivariate+-- <https://en.wikipedia.org/wiki/Laurent_polynomial Laurent polynomials>.+--++{-# LANGUAGE CPP                        #-}+{-# LANGUAGE DataKinds                  #-}+{-# LANGUAGE FlexibleContexts           #-}+{-# LANGUAGE FlexibleInstances          #-}+{-# LANGUAGE LambdaCase                 #-}+{-# LANGUAGE PatternSynonyms            #-}+{-# LANGUAGE ScopedTypeVariables        #-}+{-# LANGUAGE StandaloneDeriving         #-}+{-# LANGUAGE TypeFamilies               #-}+{-# LANGUAGE TypeOperators              #-}+{-# LANGUAGE UndecidableInstances       #-}+{-# LANGUAGE ViewPatterns               #-}++#if __GLASGOW_HASKELL__ >= 806+{-# LANGUAGE QuantifiedConstraints      #-}+#endif++module Data.Poly.Internal.Multi.Laurent+  ( MultiLaurent+  , VMultiLaurent+  , UMultiLaurent+  , unMultiLaurent+  , toMultiLaurent+  , leading+  , monomial+  , scale+  , pattern X+  , pattern Y+  , pattern Z+  , (^-)+  , eval+  , subst+  , deriv+  -- * Univariate polynomials+  , Laurent+  , VLaurent+  , ULaurent+  , unLaurent+  , toLaurent+  -- * Conversions+  , segregate+  , unsegregate+  ) where++import Prelude hiding (quotRem, quot, rem, gcd, lcm)+import Control.Arrow (first)+import Control.DeepSeq (NFData(..))+import Control.Exception+import Data.Euclidean (GcdDomain(..), Euclidean(..), Field)+import Data.Finite+import Data.Kind+import Data.List (intersperse, foldl1')+import Data.Semiring (Semiring(..), Ring())+import qualified Data.Semiring as Semiring+import qualified Data.Vector as V+import qualified Data.Vector.Generic as G+import qualified Data.Vector.Generic.Sized as SG+import qualified Data.Vector.Sized as SV+import qualified Data.Vector.Unboxed as U+import qualified Data.Vector.Unboxed.Sized as SU+import GHC.Exts+import GHC.TypeNats (KnownNat, Nat, type (+), type (<=))++import Data.Poly.Internal.Multi.Core (derivPoly)+import Data.Poly.Internal.Multi (Poly, MultiPoly(..))+import qualified Data.Poly.Internal.Multi as Multi+import Data.Poly.Internal.Multi.Field ()+import Data.Poly.Internal.Multi.GcdDomain ()++-- | Sparse+-- <https://en.wikipedia.org/wiki/Laurent_polynomial Laurent polynomials>+-- 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:+--+-- >>> (X + 1) + (Y^-1 - 1) :: VMultiLaurent 2 Integer+-- 1 * X + 1 * Y^-1+-- >>> (X + 1) * (Z - X^-1) :: UMultiLaurent 3 Int+-- 1 * X * Z + 1 * Z + (-1) + (-1) * X^-1+--+-- 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,+-- it is defined only for the sake of 'Data.Set.Set', 'Data.Map.Map', etc.+--+data MultiLaurent (v :: Type -> Type) (n :: Nat) (a :: Type) =+  MultiLaurent !(SU.Vector n Int) !(MultiPoly v n a)++deriving instance Eq  (v (SU.Vector n Word, a)) => Eq  (MultiLaurent v n a)+deriving instance Ord (v (SU.Vector n Word, a)) => Ord (MultiLaurent v n a)++-- | Multivariate Laurent polynomials backed by boxed vectors.+type VMultiLaurent (n :: Nat) (a :: Type) = MultiLaurent V.Vector n a++-- | Multivariate Laurent polynomials backed by unboxed vectors.+type UMultiLaurent (n :: Nat) (a :: Type) = MultiLaurent U.Vector n a++-- | <https://en.wikipedia.org/wiki/Laurent_polynomial Laurent polynomials>+-- of one variable with coefficients from @a@,+-- backed by a 'G.Vector' @v@ (boxed, unboxed, storable, etc.).+--+-- Use pattern 'X' and operator '^-' for construction:+--+-- >>> (X + 1) + (X^-1 - 1) :: VLaurent Integer+-- 1 * X + 1 * X^-1+-- >>> (X + 1) * (1 - X^-1) :: ULaurent Int+-- 1 * X + (-1) * X^-1+--+-- 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,+-- it is defined only for the sake of 'Data.Set.Set', 'Data.Map.Map', etc.+--+type Laurent (v :: Type -> Type) (a :: Type) = MultiLaurent v 1 a++-- | Laurent polynomials backed by boxed vectors.+type VLaurent (a :: Type) = Laurent V.Vector a++-- | Laurent polynomials backed by unboxed vectors.+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+  type Item (MultiLaurent v n a) = (SU.Vector n Int, a)++  fromList [] = MultiLaurent 0 zero+  fromList xs = toMultiLaurent minPow (fromList ys)+    where+      minPow = foldl1' (SU.zipWith min) (map fst xs)+      ys = map (first (SU.map fromIntegral . subtract minPow)) xs++  toList (MultiLaurent off (MultiPoly poly)) =+    map (first ((+ off) . SU.map fromIntegral)) $ G.toList poly++-- | Deconstruct a 'MultiLaurent' polynomial into an offset (largest possible)+-- and a regular polynomial.+--+-- >>> unMultiLaurent (2 * X + 1 :: UMultiLaurent 2 Int)+-- (Vector [0,0],2 * X + 1)+-- >>> unMultiLaurent (1 + 2 * X^-1 :: UMultiLaurent 2 Int)+-- (Vector [-1,0],1 * X + 2)+-- >>> unMultiLaurent (2 * X^2 + X :: UMultiLaurent 2 Int)+-- (Vector [1,0],2 * X + 1)+-- >>> unMultiLaurent (0 :: UMultiLaurent 2 Int)+-- (Vector [0,0],0)+unMultiLaurent :: MultiLaurent v n a -> (SU.Vector n Int, MultiPoly v n a)+unMultiLaurent (MultiLaurent off poly) = (off, poly)++-- | Deconstruct a 'Laurent' polynomial into an offset (largest possible)+-- and a regular polynomial.+--+-- >>> unLaurent (2 * X + 1 :: ULaurent Int)+-- (0,2 * X + 1)+-- >>> unLaurent (1 + 2 * X^-1 :: ULaurent Int)+-- (-1,1 * X + 2)+-- >>> unLaurent (2 * X^2 + X :: ULaurent Int)+-- (1,2 * X + 1)+-- >>> unLaurent (0 :: ULaurent Int)+-- (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.+-- One can imagine it as 'Data.Poly.Multi.Semiring.scale', but allowing negative offsets.+--+-- >>> :set -XDataKinds+-- >>> import Data.Vector.Generic.Sized (fromTuple)+-- >>> toMultiLaurent (fromTuple (2, 0)) (2 * Data.Poly.Multi.X + 1) :: UMultiLaurent 2 Int+-- 2 * X^3 + 1 * X^2+-- >>> toMultiLaurent (fromTuple (0, -2)) (2 * Data.Poly.Multi.X + 1) :: UMultiLaurent 2 Int+-- 2 * X * Y^-2 + 1 * Y^-2+toMultiLaurent+  :: (KnownNat n, G.Vector v (SU.Vector n Word, a))+  => SU.Vector n Int+  -> MultiPoly v n a+  -> MultiLaurent v n a+toMultiLaurent off (MultiPoly xs)+  | G.null xs = MultiLaurent 0 (MultiPoly G.empty)+  | otherwise = MultiLaurent (SU.zipWith (\o m -> o + fromIntegral m) off minPow) (MultiPoly ys)+    where+      minPow = G.foldl'(\acc (x, _) -> SU.zipWith min acc x) (SU.replicate maxBound) xs+      ys+        | SU.all (== 0) minPow = xs+        | otherwise = G.map (first (SU.zipWith subtract minPow)) xs+{-# INLINE toMultiLaurent #-}++-- | Construct '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+toLaurent+  :: G.Vector v (SU.Vector 1 Word, a)+  => Int+  -> Poly v a+  -> Laurent v a+toLaurent = toMultiLaurent . SU.singleton++instance NFData (v (SU.Vector n Word, a)) => NFData (MultiLaurent v n a) where+  rnf (MultiLaurent off poly) = rnf off `seq` rnf poly++instance (Show a, KnownNat n, G.Vector v (SU.Vector n Word, a)) => Show (MultiLaurent v n a) where+  showsPrec d (MultiLaurent off (MultiPoly xs))+    | G.null xs+      = showString "0"+    | otherwise+      = showParen (d > 0)+      $ foldl (.) id+      $ intersperse (showString " + ")+      $ G.foldl (\acc (is, c) -> showCoeff (SU.map fromIntegral is + off) c : acc) [] xs+    where+      showCoeff is c+        = showsPrec 7 c . foldl (.) id+        ( map ((showString " * " .) . uncurry showPower)+        $ filter ((/= 0) . fst)+        $ zip (SU.toList is) (finites :: [Finite n]))++      -- Negative powers should be displayed without surrounding brackets+      showPower :: Int -> Finite n -> String -> String+      showPower 1 n = showString (showVar n)+      showPower i n = showString (showVar n) . showString ("^" ++ show i)++      showVar :: Finite n -> String+      showVar = \case+        0 -> "X"+        1 -> "Y"+        2 -> "Z"+        k -> "X" ++ show k++-- | Return a 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+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++-- | Note that 'abs' = 'id' and 'signum' = 'const' 1.+instance (Eq a, Num a, KnownNat n, G.Vector v (SU.Vector n Word, a)) => Num (MultiLaurent v n a) where+  MultiLaurent off1 poly1 * MultiLaurent off2 poly2 = toMultiLaurent (off1 + off2) (poly1 * poly2)+  MultiLaurent off1 poly1 + MultiLaurent off2 poly2 = toMultiLaurent off (poly1' + poly2')+    where+      off    = SU.zipWith min off1 off2+      poly1' = Multi.scale (SU.zipWith (\x y -> fromIntegral (x - y)) off1 off) 1 poly1+      poly2' = Multi.scale (SU.zipWith (\x y -> fromIntegral (x - y)) off2 off) 1 poly2+  MultiLaurent off1 poly1 - MultiLaurent off2 poly2 = toMultiLaurent off (poly1' - poly2')+    where+      off    = SU.zipWith min off1 off2+      poly1' = Multi.scale (SU.zipWith (\x y -> fromIntegral (x - y)) off1 off) 1 poly1+      poly2' = Multi.scale (SU.zipWith (\x y -> fromIntegral (x - y)) off2 off) 1 poly2+  negate (MultiLaurent off poly) = MultiLaurent off (negate poly)+  abs = id+  signum = const 1+  fromInteger n = MultiLaurent 0 (fromInteger n)+  {-# INLINE (+) #-}+  {-# INLINE (-) #-}+  {-# INLINE negate #-}+  {-# INLINE fromInteger #-}+  {-# INLINE (*) #-}++instance (Eq a, Semiring a, KnownNat n, G.Vector v (SU.Vector n Word, a)) => Semiring (MultiLaurent v n a) where+  zero = MultiLaurent 0 zero+  one  = MultiLaurent 0 one+  MultiLaurent off1 poly1 `times` MultiLaurent off2 poly2 =+    toMultiLaurent (off1 + off2) (poly1 `times` poly2)+  MultiLaurent off1 poly1 `plus` MultiLaurent off2 poly2 = toMultiLaurent off (poly1' `plus` poly2')+    where+      off    = SU.zipWith min off1 off2+      poly1' = Multi.scale' (SU.zipWith (\x y -> fromIntegral (x - y)) off1 off) one poly1+      poly2' = Multi.scale' (SU.zipWith (\x y -> fromIntegral (x - y)) off2 off) one poly2+  fromNatural n = MultiLaurent 0 (fromNatural n)+  {-# INLINE zero #-}+  {-# INLINE one #-}+  {-# INLINE plus #-}+  {-# INLINE times #-}+  {-# INLINE fromNatural #-}++instance (Eq a, Ring a, KnownNat n, G.Vector v (SU.Vector n Word, a)) => Ring (MultiLaurent v n a) where+  negate (MultiLaurent off poly) = MultiLaurent off (Semiring.negate poly)++-- | Create a monomial from a power and a coefficient.+monomial+  :: (Eq a, Semiring a, KnownNat n, G.Vector v (SU.Vector n Word, a))+  => SU.Vector n Int+  -> a+  -> MultiLaurent v n a+monomial p c+  | c == zero = MultiLaurent 0 zero+  | otherwise = MultiLaurent p (Multi.monomial' 0 c)+{-# INLINE monomial #-}++-- | Multiply a polynomial by a monomial, expressed as a power and a coefficient.+--+-- >>> :set -XDataKinds+-- >>> 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+scale+  :: (Eq a, Semiring a, KnownNat n, G.Vector v (SU.Vector n Word, a))+  => SU.Vector n Int+  -> a+  -> MultiLaurent v n a+  -> MultiLaurent v n a+scale yp yc (MultiLaurent off poly) = toMultiLaurent (off + yp) (Multi.scale' 0 yc poly)++-- | Evaluate 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+eval+  :: (Field a, G.Vector v (SU.Vector n Word, a), G.Vector u a)+  => MultiLaurent v n a+  -> SG.Vector u n a+  -> a+eval (MultiLaurent off poly) xs = Multi.eval' poly xs `times`+  SU.ifoldl' (\acc i o -> acc `times` (let x = SG.index xs i in if o >= 0 then x Semiring.^ o else quot one x Semiring.^ (- o))) one off+{-# INLINE eval #-}++-- | Substitute another polynomial instead of 'Data.Poly.Multi.X'.+--+-- >>> :set -XDataKinds+-- >>> import Data.Vector.Generic.Sized (fromTuple)+-- >>> 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+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+  -> SV.Vector n (MultiLaurent w n a)+  -> MultiLaurent w n a+subst = Multi.substitute' (scale 0)+{-# INLINE subst #-}++-- | Take a derivative 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+deriv+  :: (Eq a, Ring a, KnownNat n, G.Vector v (SU.Vector n Word, a))+  => Finite n+  -> MultiLaurent v n a+  -> MultiLaurent v n a+deriv i (MultiLaurent off (MultiPoly xs)) =+  toMultiLaurent (off SU.// [(i, off `SU.index` i - 1)]) $ MultiPoly $ derivPoly+    (/= zero)+    id+    (\ps c -> Semiring.fromIntegral (fromIntegral (ps `SU.index` i) + off `SU.index` i) `times` c)+    xs+{-# INLINE deriv #-}++-- | Create a polynomial equal to the first variable.+pattern X+  :: (Eq a, Semiring a, KnownNat n, 1 <= n, G.Vector v (SU.Vector n Word, a))+  => MultiLaurent v n a+pattern X <- (isVar 0 -> True)+  where X = var 0++-- | Create a polynomial equal to the second variable.+pattern Y+  :: (Eq a, Semiring a, KnownNat n, 2 <= n, G.Vector v (SU.Vector n Word, a))+  => MultiLaurent v n a+pattern Y <- (isVar 1 -> True)+  where Y = var 1++-- | Create a polynomial equal to the third variable.+pattern Z+  :: (Eq a, Semiring a, KnownNat n, 3 <= n, G.Vector v (SU.Vector n Word, a))+  => MultiLaurent v n a+pattern Z <- (isVar 2 -> True)+  where Z = var 2++var+  :: forall v n a.+     (Eq a, Semiring a, KnownNat n, G.Vector v (SU.Vector n Word, a))+  => Finite n+  -> MultiLaurent v n a+var i+  | (one :: a) == zero = MultiLaurent 0 zero+  | otherwise          = MultiLaurent+      (SU.generate (\j -> if i == j then 1 else 0)) one+{-# INLINE var #-}++isVar+  :: forall v n a.+     (Eq a, Semiring a, KnownNat n, G.Vector v (SU.Vector n Word, a))+  => Finite n+  -> MultiLaurent v n a+  -> Bool+isVar i (MultiLaurent off (MultiPoly xs))+  | (one :: a) == zero+  = off == 0 && G.null xs+  | otherwise+  = off == SU.generate (\j -> if i == j then 1 else 0)+  && 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:+--+-- >>> 3 * X^-1 + 2 * (Y^2)^-2 :: UMultiLaurent 2 Int+-- 2 * Y^-4 + 3 * X^-1+(^-)+  :: (Eq a, Semiring a, KnownNat n, G.Vector v (SU.Vector n Word, a))+  => MultiLaurent v n a+  -> Int+  -> MultiLaurent v n a+MultiLaurent off (MultiPoly xs) ^- n+  | G.length xs == 1, G.unsafeHead xs == (0, one)+  = MultiLaurent (SU.map (* (-n)) off) (MultiPoly xs)+  | otherwise+  = throw $ PatternMatchFail "(^-) can be applied only to a monom with unit coefficient"++instance {-# OVERLAPPING #-} (Eq a, Ring a, GcdDomain a, G.Vector v (SU.Vector 1 Word, a)) => GcdDomain (Laurent v a) where+  divide (MultiLaurent off1 poly1) (MultiLaurent off2 poly2) =+    toMultiLaurent (off1 - off2) <$> divide poly1 poly2+  {-# INLINE divide #-}++  gcd (MultiLaurent _ poly1) (MultiLaurent _ poly2) =+    toMultiLaurent 0 (gcd poly1 poly2)+  {-# INLINE gcd #-}++  lcm (MultiLaurent _ poly1) (MultiLaurent _ poly2) =+    toMultiLaurent 0 (lcm poly1 poly2)+  {-# INLINE lcm #-}++  coprime (MultiLaurent _ poly1) (MultiLaurent _ poly2) =+    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 #-}++  gcd (MultiLaurent _ poly1) (MultiLaurent _ poly2) =+    toMultiLaurent 0 (gcd poly1 poly2)+  {-# INLINE gcd #-}++  lcm (MultiLaurent _ poly1) (MultiLaurent _ poly2) =+    toMultiLaurent 0 (lcm poly1 poly2)+  {-# INLINE lcm #-}++  coprime (MultiLaurent _ poly1) (MultiLaurent _ poly2) =+    coprime poly1 poly2+  {-# INLINE coprime #-}++-------------------------------------------------------------------------------++-- | 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.+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+  -> VLaurent (MultiLaurent v m a)+segregate (MultiLaurent off poly)+  = toMultiLaurent (SU.take off)+  $ MultiPoly+  $ G.map (fmap (toMultiLaurent (SU.tail off)))+  $ Multi.unMultiPoly+  $ Multi.segregate poly++-- | Interpret a univariate Laurent polynomials, whose coefficients are+-- multivariate Laurent polynomials over the first /m/ variables,+-- as a multivariate polynomial over 1+/m/ variables.+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))+  => VLaurent (MultiLaurent v m a)+  -> MultiLaurent v (1 + m) a+unsegregate (MultiLaurent off poly)+  | G.null (unMultiPoly poly)+  = MultiLaurent 0 (MultiPoly G.empty)+  | otherwise+  = toMultiLaurent (off SU.++ offs) (MultiPoly (G.concat (G.toList ys)))+  where+    xs :: V.Vector (SU.Vector 1 Word, (SU.Vector m Int, MultiPoly v m a))+    xs = G.map (fmap unMultiLaurent) $ Multi.unMultiPoly poly+    offs :: SU.Vector m Int+    offs = G.foldl' (\acc (_, (v, _)) -> SU.zipWith min acc v) (SU.replicate maxBound) xs+    ys :: V.Vector (v (SU.Vector (1 + m) Word, a))+    ys = G.map (\(v, (vs, p)) -> G.map (first ((v SU.++) . SU.zipWith3 (\a b c -> c + fromIntegral (b - a)) offs vs)) (unMultiPoly p)) xs

@@ -1,46 +0,0 @@--- |--- Module:      Data.Poly.Internal.PolyOverField--- Copyright:   (c) 2019 Andrew Lelechenko--- Licence:     BSD3--- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>------ Wrapper with a more efficient 'Euclidean' instance.-----{-# LANGUAGE ConstraintKinds            #-}-{-# LANGUAGE FlexibleInstances          #-}-{-# LANGUAGE GeneralizedNewtypeDeriving #-}-{-# LANGUAGE PatternSynonyms            #-}--module Data.Poly.Internal.PolyOverField-  ( PolyOverField(..)-  ) where--import Prelude hiding (quotRem, quot, rem, gcd, lcm, (^))-import Control.DeepSeq (NFData)-import Data.Euclidean-import Data.Semiring-import qualified Data.Vector.Generic as G--import qualified Data.Poly.Internal.Dense as Dense-import qualified Data.Poly.Internal.Dense.Field as Dense (fieldGcd)---- | Wrapper for polynomials over 'Field',--- providing a faster 'GcdDomain' instance.-newtype PolyOverField poly = PolyOverField { unPolyOverField :: poly }-  deriving (Eq, NFData, Num, Ord, Ring, Semiring, Show)--instance (Eq a, Eq (v a), Field a, G.Vector v a) => GcdDomain (PolyOverField (Dense.Poly v a)) where-  gcd (PolyOverField x) (PolyOverField y) = PolyOverField (Dense.fieldGcd x y)-  {-# INLINE gcd #-}--instance (Eq a, Eq (v a), Field a, G.Vector v a) => Euclidean (PolyOverField (Dense.Poly v a)) where-  degree (PolyOverField x) =-    degree x-  quotRem (PolyOverField x) (PolyOverField y) =-    let (q, r) = quotRem x y in-      (PolyOverField q, PolyOverField r)-  {-# INLINE quotRem #-}-  rem (PolyOverField x) (PolyOverField y) =-    PolyOverField (rem x y)-  {-# INLINE rem #-}

@@ -1,583 +0,0 @@--- |--- Module:      Data.Poly.Internal.Sparse--- Copyright:   (c) 2019 Andrew Lelechenko--- Licence:     BSD3--- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>------ Sparse polynomials of one variable.-----{-# LANGUAGE FlexibleContexts           #-}-{-# LANGUAGE GeneralizedNewtypeDeriving #-}-{-# LANGUAGE PatternSynonyms            #-}-{-# LANGUAGE ScopedTypeVariables        #-}-{-# LANGUAGE StandaloneDeriving         #-}-{-# LANGUAGE TypeFamilies               #-}-{-# LANGUAGE UndecidableInstances       #-}-{-# LANGUAGE ViewPatterns               #-}--module Data.Poly.Internal.Sparse-  ( Poly(..)-  , VPoly-  , UPoly-  , leading-  -- * Num interface-  , toPoly-  , monomial-  , scale-  , pattern X-  , eval-  , subst-  , deriv-  , integral-  -- * Semiring interface-  , toPoly'-  , monomial'-  , scale'-  , pattern X'-  , eval'-  , subst'-  , substitute'-  , deriv'-  , integral'-  ) where--import Prelude hiding (quot)-import Control.DeepSeq (NFData)-import Control.Monad-import Control.Monad.Primitive-import Control.Monad.ST-import Data.Bits-import Data.Euclidean (Field, quot)-import Data.List (intersperse)-import Data.Ord-import Data.Semiring (Semiring(..), Ring())-import qualified Data.Semiring as Semiring-import qualified Data.Vector as V-import qualified Data.Vector.Generic as G-import qualified Data.Vector.Generic.Mutable as MG-import qualified Data.Vector.Unboxed as U-import qualified Data.Vector.Algorithms.Tim as Tim-import GHC.Exts---- | Polynomials of one variable with coefficients from @a@,--- backed by a 'G.Vector' @v@ (boxed, unboxed, storable, etc.).------ Use pattern 'X' for construction:------ >>> (X + 1) + (X - 1) :: VPoly Integer--- 2 * X--- >>> (X + 1) * (X - 1) :: UPoly Int--- 1 * X^2 + (-1)------ Polynomials are stored normalized, without--- zero coefficients, so 0 * 'X' + 1 equals to 1.------ 'Ord' instance does not make much sense mathematically,--- it is defined only for the sake of 'Data.Set.Set', 'Data.Map.Map', etc.----newtype Poly v a = Poly-  { unPoly :: v (Word, a)-  -- ^ Convert 'Poly' to a vector of coefficients-  -- (first element corresponds to a constant term).-  }--deriving instance Eq     (v (Word, a)) => Eq     (Poly v a)-deriving instance Ord    (v (Word, a)) => Ord    (Poly v a)-deriving instance NFData (v (Word, a)) => NFData (Poly v a)--instance (Eq a, Semiring a, G.Vector v (Word, a)) => IsList (Poly v a) where-  type Item (Poly v a) = (Word, a)-  fromList = toPoly' . G.fromList-  fromListN = (toPoly' .) . G.fromListN-  toList = G.toList . unPoly--instance (Show a, G.Vector v (Word, a)) => Show (Poly v a) where-  showsPrec d (Poly xs)-    | G.null xs-      = showString "0"-    | otherwise-      = showParen (d > 0)-      $ foldl (.) id-      $ intersperse (showString " + ")-      $ G.foldl (\acc (i, c) -> showCoeff i c : acc) [] xs-    where-      showCoeff 0 c = showsPrec 7 c-      showCoeff 1 c = showsPrec 7 c . showString " * X"-      showCoeff i c = showsPrec 7 c . showString " * X^" . showsPrec 7 i---- | Polynomials backed by boxed vectors.-type VPoly = Poly V.Vector---- | Polynomials backed by unboxed vectors.-type UPoly = Poly U.Vector---- | Make 'Poly' from a list of (power, coefficient) pairs.--- (first element corresponds to a constant term).------ >>> :set -XOverloadedLists--- >>> toPoly [(0,1),(1,2),(2,3)] :: VPoly Integer--- 3 * X^2 + 2 * X + 1--- >>> S.toPoly [(0,0),(1,0),(2,0)] :: UPoly Int--- 0-toPoly :: (Eq a, Num a, G.Vector v (Word, a)) => v (Word, a) -> Poly v a-toPoly = Poly . normalize (/= 0) (+)--toPoly' :: (Eq a, Semiring a, G.Vector v (Word, a)) => v (Word, a) -> Poly v a-toPoly' = Poly . normalize (/= zero) plus---- | Return a 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-leading :: G.Vector v (Word, a) => Poly v a -> Maybe (Word, a)-leading (Poly v)-  | G.null v  = Nothing-  | otherwise = Just (G.last v)--normalize-  :: G.Vector v (Word, a)-  => (a -> Bool)-  -> (a -> a -> a)-  -> v (Word, a)-  -> v (Word, a)-normalize p add vs-  | G.null vs = vs-  | otherwise = runST $ do-    ws <- G.thaw vs-    l' <- normalizeM p add ws-    G.unsafeFreeze $ MG.unsafeSlice 0 l' ws--normalizeM-  :: (PrimMonad m, G.Vector v (Word, a))-  => (a -> Bool)-  -> (a -> a -> a)-  -> G.Mutable v (PrimState m) (Word, a)-  -> m Int-normalizeM p add ws = do-    let l = MG.length ws-    let go i j acc@(accP, accC)-          | j >= l =-            if p accC-              then do-                MG.write ws i acc-                pure $ i + 1-              else pure i-          | otherwise = do-            v@(vp, vc) <- MG.unsafeRead ws j-            if vp == accP-              then go i (j + 1) (accP, accC `add` vc)-              else if p accC-                then do-                  MG.write ws i acc-                  go (i + 1) (j + 1) v-                else go i (j + 1) v-    Tim.sortBy (comparing fst) ws-    wsHead <- MG.unsafeRead ws 0-    go 0 1 wsHead---- | Note that 'abs' = 'id' and 'signum' = 'const' 1.-instance (Eq a, Num a, G.Vector v (Word, a)) => Num (Poly v a) where-  Poly xs + Poly ys = Poly $ plusPoly (/= 0) (+) xs ys-  Poly xs - Poly ys = Poly $ minusPoly (/= 0) negate (-) xs ys-  negate (Poly xs) = Poly $ G.map (fmap negate) xs-  abs = id-  signum = const 1-  fromInteger n = case fromInteger n of-    0 -> Poly G.empty-    m -> Poly $ G.singleton (0, m)-  Poly xs * Poly ys = Poly $ convolution (/= 0) (+) (*) xs ys-  {-# INLINE (+) #-}-  {-# INLINE (-) #-}-  {-# INLINE negate #-}-  {-# INLINE fromInteger #-}-  {-# INLINE (*) #-}--instance (Eq a, Semiring a, G.Vector v (Word, a)) => Semiring (Poly v a) where-  zero = Poly G.empty-  one-    | (one :: a) == zero = zero-    | otherwise = Poly $ G.singleton (0, one)-  plus (Poly xs) (Poly ys) = Poly $ plusPoly (/= zero) plus xs ys-  times (Poly xs) (Poly ys) = Poly $ convolution (/= zero) plus times xs ys-  {-# INLINE zero #-}-  {-# INLINE one #-}-  {-# INLINE plus #-}-  {-# INLINE times #-}--  fromNatural n = if n' == zero then zero else Poly $ G.singleton (0, n')-    where-      n' :: a-      n' = fromNatural n-  {-# INLINE fromNatural #-}--instance (Eq a, Ring a, G.Vector v (Word, a)) => Ring (Poly v a) where-  negate (Poly xs) = Poly $ G.map (fmap Semiring.negate) xs--plusPoly-  :: G.Vector v (Word, a)-  => (a -> Bool)-  -> (a -> a -> a)-  -> v (Word, a)-  -> v (Word, a)-  -> v (Word, a)-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 (Word, a))-  => (a -> Bool)-  -> (a -> a -> a)-  -> v (Word, a)-  -> v (Word, a)-  -> G.Mutable v (PrimState m) (Word, a)-  -> m Int-plusPolyM p add xs ys zs = go 0 0 0-  where-    lenXs = G.length xs-    lenYs = G.length ys--    go ix iy iz-      | ix == lenXs, iy == lenYs = pure iz-      | ix == lenXs = do-        G.unsafeCopy-          (MG.unsafeSlice iz (lenYs - iy) zs)-          (G.unsafeSlice iy (lenYs - iy) ys)-        pure $ iz + lenYs - iy-      | iy == lenYs = do-        G.unsafeCopy-          (MG.unsafeSlice iz (lenXs - ix) zs)-          (G.unsafeSlice ix (lenXs - ix) xs)-        pure $ iz + lenXs - ix-      | (xp, xc) <- G.unsafeIndex xs ix-      , (yp, yc) <- G.unsafeIndex ys iy-      = case xp `compare` yp of-        LT -> do-          MG.unsafeWrite zs iz (xp, xc)-          go (ix + 1) iy (iz + 1)-        EQ -> do-          let zc = xc `add` yc-          if p zc then do-            MG.unsafeWrite zs iz (xp, zc)-            go (ix + 1) (iy + 1) (iz + 1)-          else-            go (ix + 1) (iy + 1) iz-        GT -> do-          MG.unsafeWrite zs iz (yp, yc)-          go ix (iy + 1) (iz + 1)-{-# INLINABLE plusPolyM #-}--minusPoly-  :: G.Vector v (Word, a)-  => (a -> Bool)-  -> (a -> a)-  -> (a -> a -> a)-  -> v (Word, a)-  -> v (Word, a)-  -> v (Word, a)-minusPoly p neg sub xs ys = runST $ do-  zs <- MG.unsafeNew (lenXs + lenYs)-  let go ix iy iz-        | ix == lenXs, iy == lenYs = pure iz-        | ix == lenXs = do-          forM_ [iy .. lenYs - 1] $ \i ->-            MG.unsafeWrite zs (iz + i - iy)-              (fmap neg (G.unsafeIndex ys i))-          pure $ iz + lenYs - iy-        | iy == lenYs = do-          G.unsafeCopy-            (MG.unsafeSlice iz (lenXs - ix) zs)-            (G.unsafeSlice ix (lenXs - ix) xs)-          pure $ iz + lenXs - ix-        | (xp, xc) <- G.unsafeIndex xs ix-        , (yp, yc) <- G.unsafeIndex ys iy-        = case xp `compare` yp of-          LT -> do-            MG.unsafeWrite zs iz (xp, xc)-            go (ix + 1) iy (iz + 1)-          EQ -> do-            let zc = xc `sub` yc-            if p zc then do-              MG.unsafeWrite zs iz (xp, zc)-              go (ix + 1) (iy + 1) (iz + 1)-            else-              go (ix + 1) (iy + 1) iz-          GT -> do-            MG.unsafeWrite zs iz (yp, neg yc)-            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 (Word, a))-  => (a -> Bool)-  -> (a -> a -> a)-  -> v (Word, a)-  -> (Word, a)-  -> G.Mutable v (PrimState m) (Word, a)-  -> m Int-scaleM p mul xs (yp, yc) zs = go 0 0-  where-    lenXs = G.length xs--    go ix iz-      | ix == lenXs = pure iz-      | (xp, xc) <- G.unsafeIndex xs ix-      = do-        let zc = xc `mul` yc-        if p zc then do-          MG.unsafeWrite zs iz (xp + yp, zc)-          go (ix + 1) (iz + 1)-        else-          go (ix + 1) iz-{-# INLINABLE scaleM #-}--scaleInternal-  :: G.Vector v (Word, a)-  => (a -> Bool)-  -> (a -> a -> a)-  -> Word-  -> a-  -> Poly v a-  -> Poly v a-scaleInternal p mul yp yc (Poly xs) = runST $ do-  zs <- MG.unsafeNew (G.length xs)-  len <- scaleM p (flip mul) xs (yp, yc) zs-  fmap Poly $ G.unsafeFreeze $ MG.unsafeSlice 0 len zs-{-# INLINABLE scaleInternal #-}---- | Multiply a polynomial by a monomial, expressed as a power and a coefficient.------ >>> scale 2 3 (X^2 + 1) :: UPoly Int--- 3 * X^4 + 3 * X^2-scale :: (Eq a, Num a, G.Vector v (Word, a)) => Word -> a -> Poly v a -> Poly v a-scale = scaleInternal (/= 0) (*)--scale' :: (Eq a, Semiring a, G.Vector v (Word, a)) => Word -> a -> Poly v a -> Poly v a-scale' = scaleInternal (/= zero) times--convolution-  :: forall v a.-     G.Vector v (Word, a)-  => (a -> Bool)-  -> (a -> a -> a)-  -> (a -> a -> a)-  -> v (Word, a)-  -> v (Word, a)-  -> v (Word, a)-convolution p add mult xs ys-  | G.length xs >= G.length ys-  = go mult xs ys-  | otherwise-  = go (flip mult) ys xs-  where-    go :: (a -> a -> a) -> v (Word, a) -> v (Word, a) -> v (Word, a)-    go mul long short = runST $ do-      let lenLong   = G.length long-          lenShort  = G.length short-          lenBuffer = lenLong * lenShort-      slices <- MG.unsafeNew lenShort-      buffer <- MG.unsafeNew lenBuffer--      forM_ [0 .. lenShort - 1] $ \iShort -> do-        let (pShort, cShort) = G.unsafeIndex short iShort-            from = iShort * lenLong-            bufferSlice = MG.unsafeSlice from lenLong buffer-        len <- scaleM p mul long (pShort, cShort) bufferSlice-        MG.unsafeWrite slices iShort (from, len)--      slices' <- G.unsafeFreeze slices-      buffer' <- G.unsafeFreeze buffer-      bufferNew <- MG.unsafeNew lenBuffer-      gogo slices' buffer' bufferNew--    gogo-      :: PrimMonad m-      => U.Vector (Int, Int)-      -> v (Word, a)-      -> G.Mutable v (PrimState m) (Word, a)-      -> m (v (Word, a))-    gogo slices buffer bufferNew-      | G.length slices == 0-      = pure G.empty-      | G.length slices == 1-      , (from, len) <- G.unsafeIndex slices 0-      = pure $ G.unsafeSlice from len buffer-      | otherwise = do-        let nSlices = G.length slices-        slicesNew <- MG.unsafeNew ((nSlices + 1) `shiftR` 1)-        forM_ [0 .. (nSlices - 2) `shiftR` 1] $ \i -> do-          let (from1, len1) = G.unsafeIndex slices (2 * i)-              (from2, len2) = G.unsafeIndex slices (2 * i + 1)-              slice1 = G.unsafeSlice from1 len1 buffer-              slice2 = G.unsafeSlice from2 len2 buffer-              slice3 = MG.unsafeSlice from1 (len1 + len2) bufferNew-          len3 <- plusPolyM p add slice1 slice2 slice3-          MG.unsafeWrite slicesNew i (from1, len3)--        when (odd nSlices) $ do-          let (from, len) = G.unsafeIndex slices (nSlices - 1)-              slice1 = G.unsafeSlice from len buffer-              slice3 = MG.unsafeSlice from len bufferNew-          G.unsafeCopy slice3 slice1-          MG.unsafeWrite slicesNew (nSlices `shiftR` 1) (from, len)--        slicesNew' <- G.unsafeFreeze slicesNew-        buffer'    <- G.unsafeThaw   buffer-        bufferNew' <- G.unsafeFreeze bufferNew-        gogo slicesNew' bufferNew' buffer'-{-# INLINABLE convolution #-}---- | Create a monomial from a power and a coefficient.-monomial :: (Eq a, Num a, G.Vector v (Word, a)) => Word -> a -> Poly v a-monomial _ 0 = Poly G.empty-monomial p c = Poly $ G.singleton (p, c)--monomial' :: (Eq a, Semiring a, G.Vector v (Word, a)) => Word -> a -> Poly v a-monomial' p c-  | c == zero = Poly G.empty-  | otherwise = Poly $ G.singleton (p, c)--data Strict3 a b c = Strict3 !a !b !c--fst3 :: Strict3 a b c -> a-fst3 (Strict3 a _ _) = a---- | Evaluate at a given point.------ >>> eval (X^2 + 1 :: UPoly Int) 3--- 10-eval :: (Num a, G.Vector v (Word, a)) => Poly v a -> a -> a-eval = substitute (*)-{-# INLINE eval #-}--eval' :: (Semiring a, G.Vector v (Word, a)) => Poly v a -> a -> a-eval' = substitute' times-{-# INLINE eval' #-}---- | Substitute another polynomial instead of 'X'.------ >>> subst (X^2 + 1 :: UPoly Int) (X + 1 :: UPoly Int)--- 1 * X^2 + 2 * X + 2-subst-  :: (Eq a, Num a, G.Vector v (Word, a), G.Vector w (Word, a))-  => Poly v a-  -> Poly w a-  -> Poly w a-subst = substitute (scale 0)-{-# INLINE subst #-}--subst'-  :: (Eq a, Semiring a, G.Vector v (Word, a), G.Vector w (Word, a))-  => Poly v a-  -> Poly w a-  -> Poly w a-subst' = substitute' (scale' 0)-{-# INLINE subst' #-}--substitute :: (G.Vector v (Word, a), Num b) => (a -> b -> b) -> Poly v a -> b -> b-substitute f (Poly cs) x = fst3 $ G.foldl' go (Strict3 0 0 1) cs-  where-    go (Strict3 acc q xq) (p, c) =-      let xp = xq * x ^ (p - q) in-        Strict3 (acc + f c xp) p xp-{-# INLINE substitute #-}--substitute' :: (G.Vector v (Word, a), Semiring b) => (a -> b -> b) -> Poly v a -> b -> b-substitute' f (Poly cs) x = fst3 $ G.foldl' go (Strict3 zero 0 one) cs-  where-    go (Strict3 acc q xq) (p, c) =-      let xp = xq `times` (if p == q then one else x Semiring.^ (p - q)) in-        Strict3 (acc `plus` f c xp) p xp-{-# INLINE substitute' #-}---- | Take a derivative.------ >>> deriv (X^3 + 3 * X) :: UPoly Int--- 3 * X^2 + 3-deriv :: (Eq a, Num a, G.Vector v (Word, a)) => Poly v a -> Poly v a-deriv (Poly xs) = Poly $ derivPoly-  (/= 0)-  (\p c -> fromIntegral p * c)-  xs-{-# INLINE deriv #-}--deriv' :: (Eq a, Semiring a, G.Vector v (Word, a)) => Poly v a -> Poly v a-deriv' (Poly xs) = Poly $ derivPoly-  (/= zero)-  (\p c -> fromNatural (fromIntegral p) `times` c)-  xs-{-# INLINE deriv' #-}--derivPoly-  :: G.Vector v (Word, a)-  => (a -> Bool)-  -> (Word -> a -> a)-  -> v (Word, a)-  -> v (Word, a)-derivPoly p mul xs-  | G.null xs = G.empty-  | otherwise = runST $ do-    let lenXs = G.length xs-    zs <- MG.unsafeNew lenXs-    let go ix iz-          | ix == lenXs = pure iz-          | (xp, xc) <- G.unsafeIndex xs ix-          = do-            let zc = xp `mul` xc-            if xp > 0 && p zc then do-              MG.unsafeWrite zs iz (xp - 1, zc)-              go (ix + 1) (iz + 1)-            else-              go (ix + 1) iz-    lenZs <- go 0 0-    G.unsafeFreeze $ MG.unsafeSlice 0 lenZs zs-{-# INLINABLE derivPoly #-}---- | Compute an indefinite integral of a polynomial,--- setting constant term to zero.------ >>> integral (3 * X^2 + 3) :: UPoly Double--- 1.0 * X^3 + 3.0 * X-integral :: (Eq a, Fractional a, G.Vector v (Word, a)) => Poly v a -> Poly v a-integral (Poly xs)-  = Poly-  $ G.map (\(p, c) -> (p + 1, c / (fromIntegral p + 1))) xs-{-# INLINE integral #-}--integral' :: (Eq a, Field a, G.Vector v (Word, a)) => Poly v a -> Poly v a-integral' (Poly xs)-  = Poly-  $ G.map (\(p, c) -> (p + 1, c `quot` Semiring.fromIntegral (p + 1))) xs-{-# INLINE integral' #-}---- | Create an identity polynomial.-pattern X :: (Eq a, Num a, G.Vector v (Word, a), Eq (v (Word, a))) => Poly v a-pattern X <- ((==) var -> True)-  where X = var--var :: forall a v. (Eq a, Num a, G.Vector v (Word, a), Eq (v (Word, a))) => Poly v a-var-  | (1 :: a) == 0 = Poly G.empty-  | otherwise     = Poly $ G.singleton (1, 1)-{-# INLINE var #-}---- | Create an identity polynomial.-pattern X' :: (Eq a, Semiring a, G.Vector v (Word, a), Eq (v (Word, a))) => Poly v a-pattern X' <- ((==) var' -> True)-  where X' = var'--var' :: forall a v. (Eq a, Semiring a, G.Vector v (Word, a), Eq (v (Word, a))) => Poly v a-var'-  | (one :: a) == zero = Poly G.empty-  | otherwise          = Poly $ G.singleton (1, one)-{-# INLINE var' #-}

@@ -1,56 +0,0 @@--- |--- Module:      Data.Poly.Internal.Sparse.Field--- Copyright:   (c) 2019 Andrew Lelechenko--- Licence:     BSD3--- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>------ GcdDomain for Field underlying-----{-# LANGUAGE ConstraintKinds            #-}-{-# LANGUAGE FlexibleContexts           #-}-{-# LANGUAGE FlexibleInstances          #-}-{-# LANGUAGE PatternSynonyms            #-}-{-# LANGUAGE ScopedTypeVariables        #-}-{-# LANGUAGE TypeFamilies               #-}-{-# LANGUAGE UndecidableInstances       #-}--{-# OPTIONS_GHC -fno-warn-orphans #-}--module Data.Poly.Internal.Sparse.Field () where--import Prelude hiding (quotRem, quot, rem, gcd)-import Control.Arrow-import Control.Exception-import Data.Euclidean (Euclidean(..), Field)-import Data.Semiring (minus, plus, times, zero)-import qualified Data.Vector.Generic as G--import Data.Poly.Internal.Sparse-import Data.Poly.Internal.Sparse.GcdDomain ()--instance (Eq a, Eq (v (Word, a)), Field a, G.Vector v (Word, a)) => Euclidean (Poly v a) where-  degree (Poly xs)-    | G.null xs = 0-    | otherwise = 1 + fromIntegral (fst (G.last xs))--  quotRem = quotientRemainder--quotientRemainder-  :: (Eq a, Field a, G.Vector v (Word, a))-  => Poly v a-  -> Poly v a-  -> (Poly v a, Poly v a)-quotientRemainder ts ys = case leading ys of-  Nothing -> throw DivideByZero-  Just (yp, yc) -> go ts-    where-      go xs = case leading xs of-        Nothing -> (zero, zero)-        Just (xp, xc) -> case xp `compare` yp of-          LT -> (zero, xs)-          EQ -> (zs, xs')-          GT -> first (`plus` zs) $ go xs'-          where-            zs = Poly $ G.singleton (xp `minus` yp, xc `quot` yc)-            xs' = xs `minus` zs `times` ys

@@ -1,74 +0,0 @@--- |--- Module:      Data.Poly.Internal.Sparse.GcdDomain--- Copyright:   (c) 2019 Andrew Lelechenko--- Licence:     BSD3--- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>------ GcdDomain for GcdDomain underlying-----{-# LANGUAGE FlexibleContexts           #-}-{-# LANGUAGE FlexibleInstances          #-}-{-# LANGUAGE PatternSynonyms            #-}-{-# LANGUAGE ScopedTypeVariables        #-}-{-# LANGUAGE TypeFamilies               #-}-{-# LANGUAGE UndecidableInstances       #-}--{-# OPTIONS_GHC -fno-warn-orphans #-}--module Data.Poly.Internal.Sparse.GcdDomain-  () where--import Prelude hiding (gcd, lcm, (^))-import Control.Exception-import Data.Euclidean-import Data.Maybe-import Data.Semiring (Semiring(..), Ring(), minus)-import qualified Data.Vector.Generic as G--import Data.Poly.Internal.Sparse---- | Consider using 'Data.Poly.Sparse.Semiring.PolyOverField' wrapper,--- which provides a much faster implementation of--- 'Data.Euclidean.gcd' for polynomials over 'Field'.-instance (Eq a, Ring a, GcdDomain a, Eq (v (Word, a)), G.Vector v (Word, a)) => GcdDomain (Poly v a) where-  divide xs ys = case leading ys of-    Nothing -> throw DivideByZero-    Just (yp, yc) -> case leading xs of-      Nothing -> Just xs-      Just (xp, xc)-        | xp < yp -> Nothing-        | otherwise -> do-          zc <- divide xc yc-          let z = Poly $ G.singleton (xp - yp, zc)-          rest <- divide (xs `minus` z `times` ys) ys-          pure $ rest `plus` z--  gcd xs ys-    | G.null (unPoly xs) = ys-    | G.null (unPoly ys) = xs-    | otherwise = maybe err (times xy) (divide zs (monomial' 0 (cont zs)))-      where-        err = error "gcd: violated internal invariant"-        zs = gcdHelper xs ys-        cont ts = G.foldl' (\acc (_, t) -> gcd acc t) zero (unPoly ts)-        xy = monomial' 0 (gcd (cont xs) (cont ys))--gcdHelper-  :: (Eq a, Ring a, GcdDomain a, G.Vector v (Word, a))-  => Poly v a-  -> Poly v a-  -> Poly v a-gcdHelper xs ys = case leading xs of-  Nothing -> ys-  Just (xp, xc) -> case leading ys of-    Nothing -> xs-    Just (yp, yc) -> case xp `compare` yp of-      LT -> gcdHelper xs (ys `times` monomial' 0 gy `minus` xs `times` monomial' (yp - xp) gx)-      EQ -> gcdHelper xs (ys `times` monomial' 0 gy `minus` xs `times` monomial' 0 gx)-      GT -> gcdHelper (xs `times` monomial' 0 gx `minus` ys `times` monomial' (xp - yp) gy) ys-      where-        g = lcm xc yc-        gx = fromMaybe err $ divide g xc-        gy = fromMaybe err $ divide g yc-        err = error "gcd: violated internal invariant"

@@ -7,11 +7,7 @@ -- <https://en.wikipedia.org/wiki/Laurent_polynomial Laurent polynomials>. -- -{-# LANGUAGE FlexibleInstances          #-}-{-# LANGUAGE GeneralizedNewtypeDeriving #-} {-# LANGUAGE PatternSynonyms            #-}-{-# LANGUAGE ScopedTypeVariables        #-}-{-# LANGUAGE ViewPatterns               #-}  module Data.Poly.Laurent   ( Laurent@@ -27,258 +23,6 @@   , eval   , subst   , deriv-  , LaurentOverField(..)   ) where -import Prelude hiding (quotRem, quot, rem, gcd)-import Control.Arrow (first)-import Control.DeepSeq (NFData(..))-import Data.Euclidean (GcdDomain(..), Euclidean(..), Field)-import Data.List (intersperse)-import Data.Semiring (Semiring(..), Ring())-import qualified Data.Semiring as Semiring-import qualified Data.Vector as V-import qualified Data.Vector.Generic as G-import qualified Data.Vector.Unboxed as U--import Data.Poly.Internal.Dense (Poly(..))-import qualified Data.Poly.Internal.Dense as Dense-import Data.Poly.Internal.Dense.Field ()-import Data.Poly.Internal.Dense.GcdDomain ()-import Data.Poly.Internal.PolyOverField---- | <https://en.wikipedia.org/wiki/Laurent_polynomial Laurent polynomials>--- of one variable with coefficients from @a@,--- backed by a 'G.Vector' @v@ (boxed, unboxed, storable, etc.).------ Use pattern 'X' and operator '^-' for construction:------ >>> (X + 1) + (X^-1 - 1) :: VLaurent Integer--- 1 * X + 0 + 1 * X^-1--- >>> (X + 1) * (1 - X^-1) :: ULaurent Int--- 1 * X + 0 + (-1) * X^-1------ Polynomials are stored normalized, without leading--- and trailing--- zero coefficients, so 0 * X + 1 + 0 * X^-1 equals to 1.------ 'Ord' instance does not make much sense mathematically,--- it is defined only for the sake of 'Data.Set.Set', 'Data.Map.Map', etc.----data Laurent v a = Laurent !Int !(Poly v a)-  deriving (Eq, Ord)---- | Deconstruct a 'Laurent' polynomial into an offset (largest possible)--- and a regular polynomial.------ >>> unLaurent (2 * X + 1 :: ULaurent Int)--- (0,2 * X + 1)--- >>> unLaurent (1 + 2 * X^-1 :: ULaurent Int)--- (-1,1 * X + 2)--- >>> unLaurent (2 * X^2 + X :: ULaurent Int)--- (1,2 * X + 1)--- >>> unLaurent (0 :: ULaurent Int)--- (0,0)-unLaurent :: Laurent v a -> (Int, Poly v a)-unLaurent (Laurent off poly) = (off, poly)---- | Construct 'Laurent' polynomial from an offset and a regular polynomial.--- One can imagine it as 'Data.Poly.scale'', but allowing negative offsets.------ >>> toLaurent 2 (2 * Data.Poly.X + 1) :: ULaurent Int--- 2 * X^3 + 1 * X^2--- >>> toLaurent (-2) (2 * Data.Poly.X + 1) :: ULaurent Int--- 2 * X^-1 + 1 * X^-2-toLaurent-  :: (Eq a, Semiring a, G.Vector v a)-  => Int-  -> Poly v a-  -> Laurent v a-toLaurent off (Poly xs) = go 0-  where-    go k-      | k >= G.length xs-      = Laurent 0 zero-      | G.unsafeIndex xs k == zero-      = go (k + 1)-      | otherwise-      = Laurent (off + k) (Poly (G.unsafeDrop k xs))-{-# INLINE toLaurent #-}--toLaurentNum-  :: (Eq a, Num a, G.Vector v a)-  => Int-  -> Poly v a-  -> Laurent v a-toLaurentNum off (Poly xs) = go 0-  where-    go k-      | k >= G.length xs-      = Laurent 0 0-      | G.unsafeIndex xs k == 0-      = go (k + 1)-      | otherwise-      = Laurent (off + k) (Poly (G.unsafeDrop k xs))-{-# INLINE toLaurentNum #-}--instance NFData (v a) => NFData (Laurent v a) where-  rnf (Laurent off poly) = rnf off `seq` rnf poly--instance (Show a, G.Vector v a) => Show (Laurent v a) where-  showsPrec d (Laurent off poly)-    | G.null (unPoly poly)-      = showString "0"-    | otherwise-      = showParen (d > 0)-      $ foldl (.) id-      $ intersperse (showString " + ")-      $ G.ifoldl (\acc i c -> showCoeff (i + off) c : acc) []-      $ unPoly poly-    where-      showCoeff 0 c = showsPrec 7 c-      showCoeff 1 c = showsPrec 7 c . showString " * X"-      showCoeff i c = showsPrec 7 c . showString (" * X^" ++ show i)---- | Laurent polynomials backed by boxed vectors.-type VLaurent = Laurent V.Vector---- | Laurent polynomials backed by unboxed vectors.-type ULaurent = Laurent U.Vector---- | Return a 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-leading :: G.Vector v a => Laurent v a -> Maybe (Int, a)-leading (Laurent off poly) = first ((+ off) . fromIntegral) <$> Dense.leading poly---- | Note that 'abs' = 'id' and 'signum' = 'const' 1.-instance (Eq a, Num a, G.Vector v a) => Num (Laurent v a) where-  Laurent off1 poly1 * Laurent off2 poly2 = toLaurentNum (off1 + off2) (poly1 * poly2)-  Laurent off1 poly1 + Laurent off2 poly2 = case off1 `compare` off2 of-    LT -> toLaurentNum off1 (poly1 + Dense.scale (fromIntegral $ off2 - off1) 1 poly2)-    EQ -> toLaurentNum off1 (poly1 + poly2)-    GT -> toLaurentNum off2 (Dense.scale (fromIntegral $ off1 - off2) 1 poly1 + poly2)-  Laurent off1 poly1 - Laurent off2 poly2 = case off1 `compare` off2 of-    LT -> toLaurentNum off1 (poly1 - Dense.scale (fromIntegral $ off2 - off1) 1 poly2)-    EQ -> toLaurentNum off1 (poly1 - poly2)-    GT -> toLaurentNum off2 (Dense.scale (fromIntegral $ off1 - off2) 1 poly1 - poly2)-  negate (Laurent off poly) = Laurent off (negate poly)-  abs = id-  signum = const 1-  fromInteger n = Laurent 0 (fromInteger n)-  {-# INLINE (+) #-}-  {-# INLINE (-) #-}-  {-# INLINE negate #-}-  {-# INLINE fromInteger #-}-  {-# INLINE (*) #-}--instance (Eq a, Semiring a, G.Vector v a) => Semiring (Laurent v a) where-  zero = Laurent 0 zero-  one  = Laurent 0 one-  Laurent off1 poly1 `times` Laurent off2 poly2 =-    toLaurent (off1 + off2) (poly1 `times` poly2)-  Laurent off1 poly1 `plus` Laurent off2 poly2 = case off1 `compare` off2 of-    LT -> toLaurent off1 (poly1 `plus` Dense.scale' (fromIntegral $ off2 - off1) one poly2)-    EQ -> toLaurent off1 (poly1 `plus` poly2)-    GT -> toLaurent off2 (Dense.scale' (fromIntegral $ off1 - off2) one poly1 `plus` poly2)-  fromNatural n = Laurent 0 (fromNatural n)-  {-# INLINE zero #-}-  {-# INLINE one #-}-  {-# INLINE plus #-}-  {-# INLINE times #-}-  {-# INLINE fromNatural #-}--instance (Eq a, Ring a, G.Vector v a) => Ring (Laurent v a) where-  negate (Laurent off poly) = Laurent off (Semiring.negate poly)---- | Create a monomial from a power and a coefficient.-monomial :: (Eq a, Semiring a, G.Vector v a) => Int -> a -> Laurent v a-monomial p c-  | c == zero = Laurent 0 zero-  | otherwise = Laurent p (Dense.monomial' 0 c)-{-# INLINE monomial #-}---- | Multiply a polynomial by a monomial, expressed as a power and a coefficient.------ >>> scale 2 3 (X^2 + 1) :: ULaurent Int--- 3 * X^4 + 0 * X^3 + 3 * X^2 + 0 * X + 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)---- | Evaluate at a given point.------ >>> eval (X^2 + 1 :: ULaurent Int) 3--- 10-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))-{-# INLINE eval #-}---- | Substitute another polynomial instead of 'Data.Poly.X'.------ >>> subst (X^2 + 1 :: UPoly Int) (X + 1 :: ULaurent Int)--- 1 * X^2 + 2 * X + 2-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.------ >>> deriv (X^3 + 3 * X) :: ULaurent Int--- 3 * X^2 + 0 * X + 3-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.-pattern X :: (Eq a, Semiring a, G.Vector v a, Eq (v a)) => Laurent v a-pattern X <- ((==) var -> True)-  where X = var--var :: forall a v. (Eq a, Semiring a, G.Vector v a, Eq (v a)) => Laurent v a-var-  | (one :: a) == zero = Laurent 0 zero-  | otherwise          = Laurent 1 one-{-# INLINE var #-}---- | This operator can be applied only to 'X',--- but is instrumental to express Laurent polynomials in mathematical fashion:------ >>> X + 2 + 3 * X^-1 :: ULaurent Int--- 1 * X + 2 + 3 * X^(-1)-(^-)-  :: (Eq a, Semiring a, G.Vector v a, Eq (v a))-  => Laurent v a-  -> Int-  -> Laurent v a-X^-n = monomial (negate n) one-_^-_ = error "(^-) can be applied only to X"---- | Consider using 'LaurentOverField' wrapper,--- which provides a much faster implementation of--- 'Data.Euclidean.gcd' for polynomials over 'Field'.-instance (Eq a, Ring a, GcdDomain a, Eq (v a), G.Vector v a) => GcdDomain (Laurent v a) where-  divide (Laurent off1 poly1) (Laurent off2 poly2) =-    toLaurent (off1 - off2) <$> divide poly1 poly2-  {-# INLINE divide #-}--  gcd (Laurent _ poly1) (Laurent _ poly2) =-    toLaurent 0 (gcd poly1 poly2)-  {-# INLINE gcd #-}---- | Wrapper for Laurent polynomials over 'Field',--- providing a faster 'GcdDomain' instance.-newtype LaurentOverField laurent = LaurentOverField { unLaurentOverField :: laurent }-  deriving (Eq, NFData, Num, Ord, Ring, Semiring, Show)--instance (Eq a, Eq (v a), Field a, G.Vector v a) => GcdDomain (LaurentOverField (Laurent v a)) where-  divide (LaurentOverField (Laurent off1 poly1)) (LaurentOverField (Laurent off2 poly2)) =-    LaurentOverField . toLaurent (off1 - off2) . unPolyOverField <$> divide (PolyOverField poly1) (PolyOverField poly2)--  gcd (LaurentOverField (Laurent _ poly1)) (LaurentOverField (Laurent _ poly2)) =-    LaurentOverField (toLaurent 0 (unPolyOverField (gcd (PolyOverField poly1) (PolyOverField poly2))))-  {-# INLINE gcd #-}+import Data.Poly.Internal.Dense.Laurent

@@ -0,0 +1,35 @@+-- |+-- Module:      Data.Poly.Multi+-- Copyright:   (c) 2020 Andrew Lelechenko+-- Licence:     BSD3+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+--+-- Sparse multivariate polynomials with 'Num' instance.+--++{-# LANGUAGE DataKinds        #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE PatternSynonyms  #-}++module Data.Poly.Multi+  ( MultiPoly+  , VMultiPoly+  , UMultiPoly+  , unMultiPoly+  , toMultiPoly+  , monomial+  , scale+  , pattern X+  , pattern Y+  , pattern Z+  , eval+  , subst+  , deriv+  , integral+  , segregate+  , unsegregate+  ) where++import Data.Poly.Internal.Multi+import Data.Poly.Internal.Multi.Field ()+import Data.Poly.Internal.Multi.GcdDomain ()

@@ -0,0 +1,32 @@+-- |+-- Module:      Data.Poly.Multi.Laurent+-- Copyright:   (c) 2020 Andrew Lelechenko+-- Licence:     BSD3+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+--+-- Sparse multivariate+-- <https://en.wikipedia.org/wiki/Laurent_polynomial Laurent polynomials>.+--++{-# LANGUAGE PatternSynonyms            #-}++module Data.Poly.Multi.Laurent+  ( MultiLaurent+  , VMultiLaurent+  , UMultiLaurent+  , unMultiLaurent+  , toMultiLaurent+  , monomial+  , scale+  , pattern X+  , pattern Y+  , pattern Z+  , (^-)+  , eval+  , subst+  , deriv+  , segregate+  , unsegregate+  ) where++import Data.Poly.Internal.Multi.Laurent

@@ -0,0 +1,157 @@+-- |+-- Module:      Data.Poly.Multi.Semiring+-- Copyright:   (c) 2020 Andrew Lelechenko+-- Licence:     BSD3+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+--+-- Sparse multivariate polynomials with 'Semiring' instance.+--++{-# LANGUAGE DataKinds        #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE PatternSynonyms  #-}+{-# LANGUAGE TypeFamilies     #-}+{-# LANGUAGE TypeOperators    #-}++module Data.Poly.Multi.Semiring+  ( MultiPoly+  , VMultiPoly+  , UMultiPoly+  , unMultiPoly+  , toMultiPoly+  , monomial+  , scale+  , pattern X+  , pattern Y+  , pattern Z+  , eval+  , subst+  , deriv+  , integral+  , segregate+  , unsegregate+  ) where++import Data.Finite+import Data.Euclidean (Field)+import Data.Semiring (Semiring(..))+import qualified Data.Vector.Generic as G+import qualified Data.Vector.Generic.Sized as SG+import qualified Data.Vector.Sized as SV+import qualified Data.Vector.Unboxed.Sized as SU+import GHC.TypeNats (KnownNat, type (<=))++import Data.Poly.Internal.Multi (MultiPoly, VMultiPoly, UMultiPoly, unMultiPoly, segregate, unsegregate)+import qualified Data.Poly.Internal.Multi as Multi+import Data.Poly.Internal.Multi.Field ()+import Data.Poly.Internal.Multi.GcdDomain ()++-- | Make 'MultiPoly' from a list of (powers, coefficient) pairs.+--+-- >>> :set -XOverloadedLists -XDataKinds+-- >>> import Data.Vector.Generic.Sized (fromTuple)+-- >>> toMultiPoly [(fromTuple (0,0),1),(fromTuple (0,1),2),(fromTuple (1,0),3)] :: VMultiPoly 2 Integer+-- 3 * X + 2 * Y + 1+-- >>> toMultiPoly [(fromTuple (0,0),0),(fromTuple (0,1),0),(fromTuple (1,0),0)] :: UMultiPoly 2 Int+-- 0+toMultiPoly+  :: (Eq a, Semiring a, G.Vector v (SU.Vector n Word, a))+  => v (SU.Vector n Word, a)+  -> MultiPoly v n a+toMultiPoly = Multi.toMultiPoly'++-- | Create a monomial from powers and a coefficient.+monomial+  :: (Eq a, Semiring a, G.Vector v (SU.Vector n Word, a))+  => SU.Vector n Word+  -> a+  -> MultiPoly v n a+monomial = Multi.monomial'++-- | Multiply a polynomial by a monomial, expressed as powers and a coefficient.+--+-- >>> :set -XDataKinds+-- >>> 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+scale+  :: (Eq a, Semiring a, KnownNat n, G.Vector v (SU.Vector n Word, a))+  => SU.Vector n Word+  -> a+  -> MultiPoly v n a+  -> MultiPoly v n a+scale = Multi.scale'++-- | Create a polynomial equal to the first variable.+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.+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.+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.+--+-- >>> :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+eval+  :: (Semiring a, G.Vector v (SU.Vector n Word, a), G.Vector u a)+  => MultiPoly v n a+  -> SG.Vector u n a+  -> a+eval = Multi.eval'++-- | Substitute another polynomials instead of 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+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+  -> SV.Vector n (MultiPoly w m a)+  -> MultiPoly w m a+subst = Multi.subst'++-- | Take a derivative 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+deriv+  :: (Eq a, Semiring a, G.Vector v (SU.Vector n Word, a))+  => Finite n+  -> MultiPoly v n a+  -> MultiPoly v n a+deriv = Multi.deriv'++-- | Compute an indefinite integral of a polynomial+-- by the /i/-th variable,+-- setting 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+integral+  :: (Field a, G.Vector v (SU.Vector n Word, a))+  => Finite n+  -> MultiPoly v n a+  -> MultiPoly v n a+integral = Multi.integral'

@@ -35,7 +35,7 @@ -- >>> take 3 legendre :: [Data.Poly.VPoly Double] -- [1.0,1.0 * X + 0.0,1.5 * X^2 + 0.0 * X + (-0.5)] legendre :: (Eq a, Field a, Vector v a) => [Poly v a]-legendre = map (flip subst' (toPoly [1 `quot` 2, 1 `quot` 2])) legendreShifted+legendre = map (`subst'` toPoly [1 `quot` 2, 1 `quot` 2]) legendreShifted   where     subst' :: (Eq a, Semiring a, Vector v a) => Poly v a -> Poly v a -> Poly v a     subst' = subst@@ -104,7 +104,7 @@ -- | Physicists' <https://en.wikipedia.org/wiki/Hermite_polynomials Hermite polynomials>. -- -- >>> take 3 hermitePhys :: [VPoly Double]--- [1,2 * X + 0,4 * X^2 + 0 * X + (-2)]+-- [1.0,2.0 * X + 0.0,4.0 * X^2 + 0.0 * X + (-2.0)] hermitePhys :: (Eq a, Ring a, Vector v a) => [Poly v a] hermitePhys = xs   where

@@ -7,7 +7,9 @@ -- Dense polynomials and a 'Semiring'-based interface. -- -{-# LANGUAGE PatternSynonyms #-}+{-# LANGUAGE DataKinds        #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE PatternSynonyms  #-}  module Data.Poly.Semiring   ( Poly@@ -23,18 +25,26 @@   , subst   , deriv   , integral-  , PolyOverField(..)+  , denseToSparse+  , sparseToDense+  , dft+  , inverseDft+  , dftMult   ) where +import Data.Bits import Data.Euclidean (Field)-import Data.Semiring (Semiring)+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 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 ()-import Data.Poly.Internal.PolyOverField+import qualified Data.Poly.Internal.Multi as Sparse  -- | Make 'Poly' from a vector of coefficients -- (first element corresponds to a constant term).@@ -59,7 +69,7 @@ scale = Dense.scale'  -- | Create an identity polynomial.-pattern X :: (Eq a, Semiring a, G.Vector v a, Eq (v a)) => Poly v a+pattern X :: (Eq a, Semiring a, G.Vector v a) => Poly v a pattern X = Dense.X'  -- | Evaluate at a given point.@@ -90,3 +100,44 @@ -- 1.0 * X^3 + 0.0 * X^2 + 3.0 * X + 0.0 integral :: (Eq a, Field a, G.Vector v a) => Poly v a -> Poly v a integral = Dense.integral'++-- | Multiplication of polynomials using+-- <https://en.wikipedia.org/wiki/Fast_Fourier_transform discrete Fourier transform>.+-- It could be faster than '(*)' for large polynomials+-- if multiplication of coefficients is particularly expensive.+dftMult+  :: (Eq a, Field a, G.Vector v a)+  => (Int -> a) -- ^ mapping from \( N = 2^n \) to a primitive root \( \sqrt[N]{1} \)+  -> Poly v a+  -> Poly v a+  -> Poly v a+dftMult getPrimRoot (Poly xs) (Poly ys) =+  toPoly $ inverseDft primRoot $ G.zipWith times (dft primRoot xs') (dft primRoot ys')+  where+    nextPowerOf2 :: Int -> Int+    nextPowerOf2 0 = 1+    nextPowerOf2 1 = 1+    nextPowerOf2 x = 1 `unsafeShiftL` (finiteBitSize (0 :: Int) - countLeadingZeros (x - 1))++    padTo l vs = G.generate l (\k -> if k < G.length vs then vs G.! k else zero)++    zl = nextPowerOf2 (G.length xs + G.length ys)+    xs' = padTo zl xs+    ys' = padTo zl ys+    primRoot = getPrimRoot zl++-- | Convert from dense to sparse polynomials.+--+-- >>> :set -XFlexibleContexts+-- >>> denseToSparse (1 `plus` Data.Poly.X^2) :: Data.Poly.Sparse.UPoly Int+-- 1 * X^2 + 1+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+-- 1 * X^2 + 0 * X + 1+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'

@@ -7,15 +7,17 @@ -- Sparse polynomials with 'Num' instance. -- -{-# LANGUAGE PatternSynonyms #-}+{-# LANGUAGE DataKinds        #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE PatternSynonyms  #-}  module Data.Poly.Sparse   ( Poly   , VPoly   , UPoly   , unPoly-  , leading   , toPoly+  , leading   , monomial   , scale   , pattern X@@ -23,8 +25,100 @@   , subst   , deriv   , integral+  , quotRemFractional+  , denseToSparse+  , sparseToDense   ) where -import Data.Poly.Internal.Sparse-import Data.Poly.Internal.Sparse.Field ()-import Data.Poly.Internal.Sparse.GcdDomain ()+import Control.Arrow+import qualified Data.Vector.Generic as G+import qualified Data.Vector.Unboxed.Sized as SU+import qualified Data.Vector.Sized as SV++import Data.Poly.Internal.Convert+import Data.Poly.Internal.Multi (Poly, VPoly, UPoly, unPoly, leading)+import qualified Data.Poly.Internal.Multi as Multi+import Data.Poly.Internal.Multi.Field (quotRemFractional)+import Data.Poly.Internal.Multi.GcdDomain ()++-- | Make '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+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)++-- | Create a monomial from a power and a coefficient.+monomial+  :: (Eq a, Num a, G.Vector v (SU.Vector 1 Word, a))+  => Word+  -> a+  -> Poly v a+monomial = Multi.monomial . SU.singleton++-- | 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+scale+  :: (Eq a, Num a, G.Vector v (SU.Vector 1 Word, a))+  => Word+  -> a+  -> Poly v a+  -> Poly v a+scale = Multi.scale . SU.singleton++-- | Create an identity polynomial.+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.+--+-- >>> eval (X^2 + 1 :: UPoly Int) 3+-- 10+eval+  :: (Num a, G.Vector v (SU.Vector 1 Word, a))+  => Poly v a+  -> a+  -> a+eval p = Multi.eval p . SV.singleton++-- | Substitute another polynomial instead of 'X'.+--+-- >>> subst (X^2 + 1 :: UPoly Int) (X + 1 :: UPoly Int)+-- 1 * X^2 + 2 * X + 2+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++-- | Take a derivative.+--+-- >>> deriv (X^3 + 3 * X) :: UPoly Int+-- 3 * X^2 + 3+deriv+  :: (Eq a, Num a, G.Vector v (SU.Vector 1 Word, a))+  => Poly v a+  -> Poly v a+deriv = Multi.deriv 0++-- | Compute an indefinite integral of a polynomial,+-- setting constant term to zero.+--+-- >>> integral (3 * X^2 + 3) :: UPoly Double+-- 1.0 * X^3 + 3.0 * X+integral+  :: (Fractional a, G.Vector v (SU.Vector 1 Word, a))+  => Poly v a+  -> Poly v a+integral = Multi.integral 0

@@ -4,18 +4,13 @@ -- Licence:     BSD3 -- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com> ----- Sparse <https://en.wikipedia.org/wiki/Laurent_polynomial Laurent polynomials>.+-- Sparse+-- <https://en.wikipedia.org/wiki/Laurent_polynomial Laurent polynomials>. -- +{-# LANGUAGE DataKinds                  #-} {-# LANGUAGE FlexibleContexts           #-}-{-# LANGUAGE FlexibleInstances          #-}-{-# LANGUAGE GeneralizedNewtypeDeriving #-} {-# LANGUAGE PatternSynonyms            #-}-{-# LANGUAGE ScopedTypeVariables        #-}-{-# LANGUAGE StandaloneDeriving         #-}-{-# LANGUAGE TypeFamilies               #-}-{-# LANGUAGE UndecidableInstances       #-}-{-# LANGUAGE ViewPatterns               #-}  module Data.Poly.Sparse.Laurent   ( Laurent@@ -33,251 +28,83 @@   , deriv   ) where -import Prelude hiding (quotRem, quot, rem, gcd)-import Control.Arrow (first)-import Control.DeepSeq (NFData(..))-import Data.Euclidean (GcdDomain(..), Euclidean(..), Field)-import Data.List (intersperse)-import Data.Ord-import Data.Semiring (Semiring(..), Ring())-import qualified Data.Semiring as Semiring-import qualified Data.Vector as V+import Data.Euclidean (Field)+import Data.Semiring (Semiring(..), Ring) import qualified Data.Vector.Generic as G-import qualified Data.Vector.Unboxed as U-import GHC.Exts--import Data.Poly.Internal.Sparse (Poly(..))-import qualified Data.Poly.Internal.Sparse as Sparse-import Data.Poly.Internal.Sparse.Field ()-import Data.Poly.Internal.Sparse.GcdDomain ()---- | <https://en.wikipedia.org/wiki/Laurent_polynomial Laurent polynomials>--- of one variable with coefficients from @a@,--- backed by a 'G.Vector' @v@ (boxed, unboxed, storable, etc.).------ Use pattern 'X' and operator '^-' for construction:------ >>> (X + 1) + (X^-1 - 1) :: VLaurent Integer--- 1 * X + 1 * X^-1--- >>> (X + 1) * (1 - X^-1) :: ULaurent Int--- 1 * X + (-1) * X^-1------ 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,--- it is defined only for the sake of 'Data.Set.Set', 'Data.Map.Map', etc.----data Laurent v a = Laurent !Int !(Poly v a)--deriving instance Eq  (v (Word, a)) => Eq  (Laurent v a)-deriving instance Ord (v (Word, a)) => Ord (Laurent v a)--instance (Eq a, Semiring a, G.Vector v (Word, a)) => IsList (Laurent v a) where-  type Item (Laurent v a) = (Int, a)--  fromList xs = toLaurent minPow (fromList ys)-    where-      minPow = minimum $ maxBound : map fst xs-      ys = map (first (fromIntegral . (subtract minPow))) xs--  fromListN n xs = toLaurent minPow (fromListN n ys)-    where-      minPow = minimum $ maxBound : map fst xs-      ys = map (first (fromIntegral . (subtract minPow))) xs--  toList (Laurent off poly) =-    map (first ((+ off) . fromIntegral)) $ G.toList $ unPoly poly+import qualified Data.Vector.Unboxed.Sized as SU+import qualified Data.Vector.Sized as SV --- | Deconstruct a 'Laurent' polynomial into an offset (largest possible)--- and a regular polynomial.------ >>> unLaurent (2 * X + 1 :: ULaurent Int)--- (0,2 * X + 1)--- >>> unLaurent (1 + 2 * X^-1 :: ULaurent Int)--- (-1,1 * X + 2)--- >>> unLaurent (2 * X^2 + X :: ULaurent Int)--- (1,2 * X + 1)--- >>> unLaurent (0 :: ULaurent Int)--- (0,0)-unLaurent :: Laurent v a -> (Int, Poly v a)-unLaurent (Laurent off poly) = (off, poly)+import Data.Poly.Internal.Multi.Laurent hiding (monomial, scale, pattern X, (^-), eval, subst, deriv)+import qualified Data.Poly.Internal.Multi.Laurent as Multi+import Data.Poly.Internal.Multi (Poly) --- | Construct 'Laurent' polynomial from an offset and a regular polynomial.--- One can imagine it as 'Data.Poly.Sparse.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-toLaurent-  :: (Eq a, Semiring a, G.Vector v (Word, a))+-- | Create a monomial from a power and a coefficient.+monomial+  :: (Eq a, Semiring a, G.Vector v (SU.Vector 1 Word, a))   => Int-  -> Poly v a+  -> a   -> Laurent v a-toLaurent off (Poly xs)-  | G.null xs = Laurent 0 zero-  | otherwise = Laurent (off + fromIntegral minPow) (Poly ys)-    where-      minPow = fst $ G.minimumBy (comparing fst) xs-      ys = if minPow == 0 then xs else G.map (first (subtract minPow)) xs-{-# INLINE toLaurent #-}+monomial = Multi.monomial . SU.singleton -toLaurentNum-  :: (Eq a, Num a, G.Vector v (Word, a))+-- | 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+scale+  :: (Eq a, Semiring a, G.Vector v (SU.Vector 1 Word, a))   => Int-  -> Poly v a+  -> a   -> Laurent v a-toLaurentNum off (Poly xs)-  | G.null xs = Laurent 0 0-  | otherwise = Laurent (off + fromIntegral minPow) (Poly ys)-    where-      minPow = fst $ G.minimumBy (comparing fst) xs-      ys = if minPow == 0 then xs else G.map (first (subtract minPow)) xs-{-# INLINE toLaurentNum #-}--instance NFData (v (Word, a)) => NFData (Laurent v a) where-  rnf (Laurent off poly) = rnf off `seq` rnf poly--instance (Show a, G.Vector v (Word, a)) => Show (Laurent v a) where-  showsPrec d (Laurent off poly)-    | G.null (unPoly poly)-      = showString "0"-    | otherwise-      = showParen (d > 0)-      $ foldl (.) id-      $ intersperse (showString " + ")-      $ G.ifoldl (\acc i c -> showCoeff (i + off) c : acc) []-      $ unPoly poly-    where-      showCoeff 0 c = showsPrec 7 c-      showCoeff 1 c = showsPrec 7 c . showString " * X"-      showCoeff i c = showsPrec 7 c . showString (" * X^" ++ show i)---- | Laurent polynomials backed by boxed vectors.-type VLaurent = Laurent V.Vector---- | Laurent polynomials backed by unboxed vectors.-type ULaurent = Laurent U.Vector---- | Return a 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-leading :: G.Vector v (Word, a) => Laurent v a -> Maybe (Int, a)-leading (Laurent off poly) = first ((+ off) . fromIntegral) <$> Sparse.leading poly---- | Note that 'abs' = 'id' and 'signum' = 'const' 1.-instance (Eq a, Num a, G.Vector v (Word, a)) => Num (Laurent v a) where-  Laurent off1 poly1 * Laurent off2 poly2 = toLaurentNum (off1 + off2) (poly1 * poly2)-  Laurent off1 poly1 + Laurent off2 poly2 = case off1 `compare` off2 of-    LT -> toLaurentNum off1 (poly1 + Sparse.scale (fromIntegral $ off2 - off1) 1 poly2)-    EQ -> toLaurentNum off1 (poly1 + poly2)-    GT -> toLaurentNum off2 (Sparse.scale (fromIntegral $ off1 - off2) 1 poly1 + poly2)-  Laurent off1 poly1 - Laurent off2 poly2 = case off1 `compare` off2 of-    LT -> toLaurentNum off1 (poly1 - Sparse.scale (fromIntegral $ off2 - off1) 1 poly2)-    EQ -> toLaurentNum off1 (poly1 - poly2)-    GT -> toLaurentNum off2 (Sparse.scale (fromIntegral $ off1 - off2) 1 poly1 - poly2)-  negate (Laurent off poly) = Laurent off (negate poly)-  abs = id-  signum = const 1-  fromInteger n = Laurent 0 (fromInteger n)-  {-# INLINE (+) #-}-  {-# INLINE (-) #-}-  {-# INLINE negate #-}-  {-# INLINE fromInteger #-}-  {-# INLINE (*) #-}--instance (Eq a, Semiring a, G.Vector v (Word, a)) => Semiring (Laurent v a) where-  zero = Laurent 0 zero-  one  = Laurent 0 one-  Laurent off1 poly1 `times` Laurent off2 poly2 =-    toLaurent (off1 + off2) (poly1 `times` poly2)-  Laurent off1 poly1 `plus` Laurent off2 poly2 = case off1 `compare` off2 of-    LT -> toLaurent off1 (poly1 `plus` Sparse.scale' (fromIntegral $ off2 - off1) one poly2)-    EQ -> toLaurent off1 (poly1 `plus` poly2)-    GT -> toLaurent off2 (Sparse.scale' (fromIntegral $ off1 - off2) one poly1 `plus` poly2)-  fromNatural n = Laurent 0 (fromNatural n)-  {-# INLINE zero #-}-  {-# INLINE one #-}-  {-# INLINE plus #-}-  {-# INLINE times #-}-  {-# INLINE fromNatural #-}--instance (Eq a, Ring a, G.Vector v (Word, a)) => Ring (Laurent v a) where-  negate (Laurent off poly) = Laurent off (Semiring.negate poly)+  -> Laurent v a+scale = Multi.scale . SU.singleton --- | Create a monomial from a power and a coefficient.-monomial :: (Eq a, Semiring a, G.Vector v (Word, a)) => Int -> a -> Laurent v a-monomial p c-  | c == zero = Laurent 0 zero-  | otherwise = Laurent p (Sparse.monomial' 0 c)-{-# INLINE monomial #-}+-- | Create an identity polynomial.+pattern X+  :: (Eq a, Semiring a, G.Vector v (SU.Vector 1 Word, a))+  => Laurent v a+pattern X = Multi.X --- | Multiply a polynomial by a monomial, expressed as a power and a coefficient.+-- | This operator can be applied only to monomials with unit coefficients,+-- but is instrumental to express Laurent polynomials in mathematical fashion: ----- >>> scale 2 3 (X^2 + 1) :: ULaurent Int--- 3 * X^4 + 3 * X^2-scale :: (Eq a, Semiring a, G.Vector v (Word, a)) => Int -> a -> Laurent v a -> Laurent v a-scale yp yc (Laurent off poly) = toLaurent (off + yp) (Sparse.scale' 0 yc poly)+-- >>> X + 2 + 3 * (X^2)^-1 :: ULaurent Int+-- 1 * X + 2 + 3 * X^-2+(^-)+  :: (Eq a, Semiring a, G.Vector v (SU.Vector 1 Word, a))+  => Laurent v a+  -> Int+  -> Laurent v a+(^-) = (Multi.^-)  -- | Evaluate at a given point. ----- >>> eval (X^2 + 1 :: ULaurent Int) 3--- 10-eval :: (Field a, G.Vector v (Word, a)) => Laurent v a -> a -> a-eval (Laurent off poly) x = Sparse.eval' poly x `times`-  (if off >= 0 then x Semiring.^ off else quot one x Semiring.^ (- off))-{-# INLINE eval #-}+-- >>> eval (X^-2 + 1 :: ULaurent Double) 2+-- 1.25+eval+  :: (Field a, G.Vector v (SU.Vector 1 Word, a))+  => Laurent v a+  -> a+  -> a+eval p = Multi.eval p . SV.singleton --- | Substitute another polynomial instead of 'Data.Poly.Sparse.X'.+-- | Substitute another polynomial instead of 'X'. ----- >>> subst (X^2 + 1 :: UPoly Int) (X + 1 :: ULaurent Int)--- 1 * X^2 + 2 * X + 2-subst :: (Eq a, Semiring a, G.Vector v (Word, a), G.Vector w (Word, a)) => Poly v a -> Laurent w a -> Laurent w a-subst = Sparse.substitute' (scale 0)-{-# INLINE subst #-}+-- >>> 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+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  -- | Take a derivative. ----- >>> deriv (X^3 + 3 * X) :: ULaurent Int--- 3 * X^2 + 3-deriv :: (Eq a, Ring a, G.Vector v (Word, a)) => Laurent v a -> Laurent v a-deriv (Laurent off (Poly xs)) =-  toLaurent (off - 1) $ Sparse.toPoly' $ G.map (\(i, x) -> (i, x `times` Semiring.fromIntegral (fromIntegral i + off))) xs-{-# INLINE deriv #-}---- | Create an identity polynomial.-pattern X :: (Eq a, Semiring a, G.Vector v (Word, a), Eq (v (Word, a))) => Laurent v a-pattern X <- ((==) var -> True)-  where X = var--var :: forall a v. (Eq a, Semiring a, G.Vector v (Word, a), Eq (v (Word, a))) => Laurent v a-var-  | (one :: a) == zero = Laurent 0 zero-  | otherwise          = Laurent 1 one-{-# INLINE var #-}---- | This operator can be applied only to 'X',--- but is instrumental to express Laurent polynomials in mathematical fashion:------ >>> X + 2 + 3 * X^-1 :: ULaurent Int--- 1 * X + 2 + 3 * X^(-1)-(^-)-  :: (Eq a, Semiring a, G.Vector v (Word, a), Eq (v (Word, a)))+-- >>> deriv (X^-3 + 3 * X) :: ULaurent Int+-- 3 + (-3) * X^-4+deriv+  :: (Eq a, Ring a, G.Vector v (SU.Vector 1 Word, a))   => Laurent v a-  -> Int   -> Laurent v a-X^-n = monomial (negate n) one-_^-_ = error "(^-) can be applied only to X"--instance (Eq a, Ring a, GcdDomain a, Eq (v (Word, a)), G.Vector v (Word, a)) => GcdDomain (Laurent v a) where-  divide (Laurent off1 poly1) (Laurent off2 poly2) =-    toLaurent (off1 - off2) <$> divide poly1 poly2-  {-# INLINE divide #-}--  gcd (Laurent _ poly1) (Laurent _ poly2) =-    toLaurent 0 (gcd poly1 poly2)-  {-# INLINE gcd #-}+deriv = Multi.deriv 0

@@ -7,6 +7,7 @@ -- Sparse polynomials with 'Semiring' instance. -- +{-# LANGUAGE DataKinds        #-} {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE PatternSynonyms  #-} @@ -15,8 +16,8 @@   , VPoly   , UPoly   , unPoly-  , leading   , toPoly+  , leading   , monomial   , scale   , pattern X@@ -24,68 +25,124 @@   , subst   , deriv   , integral+  , denseToSparse+  , sparseToDense   ) where +import Control.Arrow import Data.Euclidean (Field)-import Data.Semiring (Semiring)+import Data.Semiring (Semiring(..)) import qualified Data.Vector.Generic as G+import qualified Data.Vector.Unboxed.Sized as SU+import qualified Data.Vector.Sized as SV -import Data.Poly.Internal.Sparse (Poly(..), VPoly, UPoly, leading)-import qualified Data.Poly.Internal.Sparse as Sparse-import Data.Poly.Internal.Sparse.Field ()-import Data.Poly.Internal.Sparse.GcdDomain ()+import qualified Data.Poly.Internal.Convert as Convert+import qualified Data.Poly.Internal.Dense as Dense+import Data.Poly.Internal.Multi (Poly, VPoly, UPoly, unPoly, leading)+import qualified Data.Poly.Internal.Multi as Multi+import Data.Poly.Internal.Multi.Field ()+import Data.Poly.Internal.Multi.GcdDomain ()  -- | Make 'Poly' from a list of (power, coefficient) pairs.--- (first element corresponds to a constant term). -- -- >>> :set -XOverloadedLists -- >>> toPoly [(0,1),(1,2),(2,3)] :: VPoly Integer -- 3 * X^2 + 2 * X + 1--- >>> S.toPoly [(0,0),(1,0),(2,0)] :: UPoly Int+-- >>> toPoly [(0,0),(1,0),(2,0)] :: UPoly Int -- 0-toPoly :: (Eq a, Semiring a, G.Vector v (Word, a)) => v (Word, a) -> Poly v a-toPoly = Sparse.toPoly'+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)  -- | Create a monomial from a power and a coefficient.-monomial :: (Eq a, Semiring a, G.Vector v (Word, a)) => Word -> a -> Poly v a-monomial = Sparse.monomial'+monomial+  :: (Eq a, Semiring a, G.Vector v (SU.Vector 1 Word, a))+  => Word+  -> a+  -> Poly v a+monomial = Multi.monomial' . SU.singleton  -- | 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-scale :: (Eq a, Semiring a, G.Vector v (Word, a)) => Word -> a -> Poly v a -> Poly v a-scale = Sparse.scale'+scale+  :: (Eq a, Semiring a, G.Vector v (SU.Vector 1 Word, a))+  => Word+  -> a+  -> Poly v a+  -> Poly v a+scale = Multi.scale' . SU.singleton  -- | Create an identity polynomial.-pattern X :: (Eq a, Semiring a, G.Vector v (Word, a), Eq (v (Word, a))) => Poly v a-pattern X = Sparse.X'+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. -- -- >>> eval (X^2 + 1 :: UPoly Int) 3 -- 10-eval :: (Semiring a, G.Vector v (Word, a)) => Poly v a -> a -> a-eval = Sparse.eval'+eval+  :: (Semiring a, G.Vector v (SU.Vector 1 Word, a))+  => Poly v a+  -> a+  -> a+eval p = Multi.eval' p . SV.singleton  -- | Substitute another polynomial instead of 'X'. -- -- >>> subst (X^2 + 1 :: UPoly Int) (X + 1 :: UPoly Int) -- 1 * X^2 + 2 * X + 2-subst :: (Eq a, Semiring a, G.Vector v (Word, a), G.Vector w (Word, a)) => Poly v a -> Poly w a -> Poly w a-subst = Sparse.subst'+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  -- | Take a derivative. -- -- >>> deriv (X^3 + 3 * X) :: UPoly Int -- 3 * X^2 + 3-deriv :: (Eq a, Semiring a, G.Vector v (Word, a)) => Poly v a -> Poly v a-deriv = Sparse.deriv'+deriv+  :: (Eq a, Semiring a, G.Vector v (SU.Vector 1 Word, a))+  => Poly v a+  -> Poly v a+deriv = Multi.deriv' 0  -- | Compute an indefinite integral of a polynomial, -- setting constant term to zero. -- -- >>> integral (3 * X^2 + 3) :: UPoly Double -- 1.0 * X^3 + 3.0 * X-integral :: (Eq a, Field a, G.Vector v (Word, a)) => Poly v a -> Poly v a-integral = Sparse.integral'+integral+  :: (Field a, G.Vector v (SU.Vector 1 Word, a))+  => Poly v a+  -> Poly v a+integral = Multi.integral' 0++-- | Convert from dense to sparse polynomials.+--+-- >>> :set -XFlexibleContexts+-- >>> denseToSparse (1 `plus` Data.Poly.X^2) :: Data.Poly.Sparse.UPoly Int+-- 1 * X^2 + 1+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'++-- | Convert from sparse to dense polynomials.+--+-- >>> :set -XFlexibleContexts+-- >>> sparseToDense (1 `plus` Data.Poly.Sparse.X^2) :: Data.Poly.UPoly Int+-- 1 * X^2 + 0 * X + 1+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'

@@ -0,0 +1,69 @@+{-# LANGUAGE DataKinds                  #-}+{-# LANGUAGE TypeOperators              #-}++module DFT+  ( testSuite+  ) where++import Data.Complex+import Data.Mod.Word+import Data.Poly.Semiring (UPoly, unPoly, toPoly, dft, inverseDft, dftMult)+import qualified Data.Vector.Unboxed as U+import GHC.TypeNats (KnownNat, natVal, type (+), type (^))+import Test.Tasty+import Test.Tasty.QuickCheck hiding (scale, numTests)++import Dense ()++testSuite :: TestTree+testSuite = testGroup "DFT"+  [ testGroup "dft matches reference"+    [ dftMatchesRef (0 :: Mod (2 ^ 0 + 1))+    , dftMatchesRef (2 :: Mod (2 ^ 1 + 1))+    , dftMatchesRef (2 :: Mod (2 ^ 2 + 1))+    , dftMatchesRef (3 :: Mod (2 ^ 4 + 1))+    , dftMatchesRef (3 :: Mod (2 ^ 8 + 1))+    ]+  , testGroup "dft is invertible"+    [ dftIsInvertible (0 :: Mod (2 ^ 0 + 1))+    , dftIsInvertible (2 :: Mod (2 ^ 1 + 1))+    , dftIsInvertible (2 :: Mod (2 ^ 2 + 1))+    , dftIsInvertible (3 :: Mod (2 ^ 4 + 1))+    , dftIsInvertible (3 :: Mod (2 ^ 8 + 1))+    ]+  , testProperty "dftMult matches reference" dftMultMatchesRef+  ]++dftMatchesRef :: KnownNat n1 => Mod n1 -> TestTree+dftMatchesRef primRoot = testProperty (show n) $ do+  xs <- U.replicateM n arbitrary+  pure $ dft primRoot xs === dftRef primRoot xs+  where+    n = fromIntegral (natVal primRoot - 1)++dftRef :: (Num a, U.Unbox a) => a -> U.Vector a -> U.Vector a+dftRef primRoot xs = U.generate (U.length xs) $+  \k -> sum (map (\j -> xs U.! j * primRoot ^ (j * k)) [0..n-1])+  where+    n = U.length xs++dftIsInvertible :: KnownNat n1 => Mod n1 -> TestTree+dftIsInvertible primRoot = testProperty (show n) $ do+  xs <- U.replicateM n arbitrary+  let ys = dft primRoot xs+      zs = inverseDft primRoot ys+  pure $ xs === zs+  where+    n = fromIntegral (natVal primRoot - 1)++dftMultMatchesRef :: UPoly Int -> UPoly Int -> Property+dftMultMatchesRef xs ys = zs === dftZs+  where+    xs', ys', dftZs' :: UPoly (Complex Double)+    xs' = toPoly $ U.map fromIntegral $ unPoly xs+    ys' = toPoly $ U.map fromIntegral $ unPoly ys+    dftZs' = dftMult (\k -> cis (2 * pi / fromIntegral k)) xs' ys'++    zs, dftZs :: UPoly (Complex Int)+    zs  = toPoly $ U.map (:+ 0) $ unPoly $ xs * ys+    dftZs  = toPoly $ U.map (\(x :+ y) -> round x :+ round y) $ unPoly dftZs'

@@ -1,24 +1,22 @@ {-# LANGUAGE DataKinds                  #-} {-# LANGUAGE FlexibleContexts           #-} {-# LANGUAGE FlexibleInstances          #-}-{-# LANGUAGE GeneralizedNewtypeDeriving #-} {-# LANGUAGE ScopedTypeVariables        #-} -{-# OPTIONS_GHC -fno-warn-orphans #-}- module Dense   ( testSuite   , ShortPoly(..)   ) where -import Prelude hiding (gcd, quotRem, rem)+import Prelude hiding (gcd, quotRem, quot, rem)+import Control.Exception import Data.Euclidean (Euclidean(..), GcdDomain(..)) import Data.Int-import Data.Mod+import Data.Mod.Word import Data.Poly import qualified Data.Poly.Semiring as S import Data.Proxy-import Data.Semiring (Semiring)+import Data.Semiring (Semiring(..)) import qualified Data.Vector as V import qualified Data.Vector.Generic as G import qualified Data.Vector.Unboxed as U@@ -28,29 +26,17 @@ import Quaternion import TestUtils -instance (Eq a, Semiring a, Arbitrary a, G.Vector v a) => Arbitrary (Poly v a) where-  arbitrary = S.toPoly . G.fromList <$> arbitrary-  shrink = fmap (S.toPoly . G.fromList) . shrink . G.toList . unPoly--instance (Eq a, Semiring a, Arbitrary a, G.Vector v a) => Arbitrary (PolyOverField (Poly v a)) where-  arbitrary = PolyOverField . S.toPoly . G.fromList . (\xs -> take (length xs `mod` 10) xs) <$> arbitrary-  shrink = fmap (PolyOverField . S.toPoly . G.fromList) . shrink . G.toList . unPoly . unPolyOverField--newtype ShortPoly a = ShortPoly { unShortPoly :: a }-  deriving (Eq, Show, Semiring, GcdDomain, Euclidean)--instance (Eq a, Semiring a, Arbitrary a, G.Vector v a) => Arbitrary (ShortPoly (Poly v a)) where-  arbitrary = ShortPoly . S.toPoly . G.fromList . (\xs -> take (length xs `mod` 10) xs) <$> arbitrary-  shrink = fmap (ShortPoly . S.toPoly . G.fromList) . shrink . G.toList . unPoly . unShortPoly- testSuite :: TestTree testSuite = testGroup "Dense"-    [ arithmeticTests-    , otherTests-    , lawsTests-    , evalTests-    , derivTests-    ]+  [ arithmeticTests+  , otherTests+  , divideByZeroTests+  , lawsTests+  , evalTests+  , derivTests+  , patternTests+  , conversionTests+  ]  lawsTests :: TestTree lawsTests = testGroup "Laws"@@ -58,54 +44,54 @@  semiringTests :: [TestTree] semiringTests =-  [ mySemiringLaws (Proxy :: Proxy (Poly U.Vector ()))-  , mySemiringLaws (Proxy :: Proxy (Poly U.Vector Int8))-  , mySemiringLaws (Proxy :: Proxy (Poly V.Vector Integer))-  , mySemiringLaws (Proxy :: Proxy (Poly U.Vector (Quaternion Int)))+  [ mySemiringLaws (Proxy :: Proxy (UPoly ()))+  , mySemiringLaws (Proxy :: Proxy (UPoly Int8))+  , mySemiringLaws (Proxy :: Proxy (VPoly Integer))+  , mySemiringLaws (Proxy :: Proxy (UPoly (Quaternion Int)))   ]  ringTests :: [TestTree] ringTests =-  [ myRingLaws (Proxy :: Proxy (Poly U.Vector ()))-  , myRingLaws (Proxy :: Proxy (Poly U.Vector Int8))-  , myRingLaws (Proxy :: Proxy (Poly V.Vector Integer))-  , myRingLaws (Proxy :: Proxy (Poly U.Vector (Quaternion Int)))+  [ myRingLaws (Proxy :: Proxy (UPoly ()))+  , myRingLaws (Proxy :: Proxy (UPoly Int8))+  , myRingLaws (Proxy :: Proxy (VPoly Integer))+  , myRingLaws (Proxy :: Proxy (UPoly (Quaternion Int)))   ]  numTests :: [TestTree] numTests =-  [ myNumLaws (Proxy :: Proxy (Poly U.Vector Int8))-  , myNumLaws (Proxy :: Proxy (Poly V.Vector Integer))-  , myNumLaws (Proxy :: Proxy (Poly U.Vector (Quaternion Int)))+  [ myNumLaws (Proxy :: Proxy (UPoly Int8))+  , myNumLaws (Proxy :: Proxy (VPoly Integer))+  , myNumLaws (Proxy :: Proxy (UPoly (Quaternion Int)))   ]  gcdDomainTests :: [TestTree] gcdDomainTests =-  [ myGcdDomainLaws (Proxy :: Proxy (ShortPoly (Poly V.Vector Integer)))-  , myGcdDomainLaws (Proxy :: Proxy (PolyOverField (Poly V.Vector (Mod 3))))-  , myGcdDomainLaws (Proxy :: Proxy (PolyOverField (Poly V.Vector Rational)))+  [ myGcdDomainLaws (Proxy :: Proxy (ShortPoly (VPoly Integer)))+  , myGcdDomainLaws (Proxy :: Proxy (ShortPoly (UPoly (Mod 3))))+  , myGcdDomainLaws (Proxy :: Proxy (ShortPoly (VPoly Rational)))   ]  euclideanTests :: [TestTree] euclideanTests =-  [ myEuclideanLaws (Proxy :: Proxy (ShortPoly (Poly V.Vector (Mod 3))))-  , myEuclideanLaws (Proxy :: Proxy (ShortPoly (Poly V.Vector Rational)))+  [ myEuclideanLaws (Proxy :: Proxy (ShortPoly (UPoly (Mod 3))))+  , myEuclideanLaws (Proxy :: Proxy (ShortPoly (VPoly Rational)))   ]  isListTests :: [TestTree] isListTests =-  [ myIsListLaws (Proxy :: Proxy (Poly U.Vector ()))-  , myIsListLaws (Proxy :: Proxy (Poly U.Vector Int8))-  , myIsListLaws (Proxy :: Proxy (Poly V.Vector Integer))-  , myIsListLaws (Proxy :: Proxy (Poly U.Vector (Quaternion Int)))+  [ myIsListLaws (Proxy :: Proxy (UPoly ()))+  , myIsListLaws (Proxy :: Proxy (UPoly Int8))+  , myIsListLaws (Proxy :: Proxy (VPoly Integer))+  , myIsListLaws (Proxy :: Proxy (UPoly (Quaternion Int)))   ]  showTests :: [TestTree] showTests =-  [ myShowLaws (Proxy :: Proxy (Poly U.Vector ()))-  , myShowLaws (Proxy :: Proxy (Poly U.Vector Int8))-  , myShowLaws (Proxy :: Proxy (Poly V.Vector Integer))-  , myShowLaws (Proxy :: Proxy (Poly U.Vector (Quaternion Int)))+  [ myShowLaws (Proxy :: Proxy (UPoly ()))+  , myShowLaws (Proxy :: Proxy (UPoly Int8))+  , myShowLaws (Proxy :: Proxy (VPoly Integer))+  , myShowLaws (Proxy :: Proxy (UPoly (Quaternion Int)))   ]  arithmeticTests :: TestTree@@ -119,6 +105,9 @@   , testProperty "multiplication matches reference" $     \(xs :: [Int]) ys -> toPoly (V.fromList (mulRef xs ys)) ===       toPoly (V.fromList xs) * toPoly (V.fromList ys)+  , tenTimesLess $+    testProperty "quotRemFractional matches quotRem" $+    \(xs :: VPoly Rational) ys -> ys /= 0 ==> quotRemFractional xs ys === quotRem xs ys   ]  addRef :: Num a => [a] -> [a] -> [a]@@ -158,16 +147,30 @@   , tenTimesLess $     testProperty "scale matches multiplication by monomial" $     \p c (xs :: UPoly a) -> scale p c xs === monomial p c * xs+  , tenTimesLess $+    testProperty "scale' matches multiplication by monomial'" $+    \p c (xs :: UPoly a) -> S.scale p c xs === S.monomial p c * xs   ]  monomialRef :: Num a => Word -> a -> [a] monomialRef p c = replicate (fromIntegral p) 0 ++ [c] +divideByZeroTests :: TestTree+divideByZeroTests = testGroup "divideByZero"+  [ testProperty "quotRem" $ testProp ((uncurry (+) .) . quotRem)+  , testProperty "quot"    $ testProp quot+  , testProperty "rem"     $ testProp rem+  , testProperty "divide"  $ testProp divide+  , testProperty "degree"  $ once $ degree (0 :: VPoly Rational) === 0+  ]+  where+    testProp f xs = ioProperty ((== Left DivideByZero) <$> try (evaluate (xs `f` (0 :: VPoly Rational))))+ evalTests :: TestTree evalTests = testGroup "eval" $ concat-  [ evalTestGroup  (Proxy :: Proxy (Poly U.Vector Int8))-  , evalTestGroup  (Proxy :: Proxy (Poly V.Vector Integer))-  , substTestGroup (Proxy :: Proxy (Poly U.Vector Int8))+  [ evalTestGroup  (Proxy :: Proxy (UPoly Int8))+  , evalTestGroup  (Proxy :: Proxy (VPoly Integer))+  , substTestGroup (Proxy :: Proxy (UPoly Int8))   ]  evalTestGroup@@ -177,18 +180,18 @@   -> [TestTree] evalTestGroup _ =   [ testProperty "eval (p + q) r = eval p r + eval q r" $-    \p q r -> e (p + q) r === e p r + e q r+    \(ShortPoly p) (ShortPoly q) r -> e (p + q) r === e p r + e q r   , testProperty "eval (p * q) r = eval p r * eval q r" $-    \p q r -> e (p * q) r === e p r * e q r+    \(ShortPoly p) (ShortPoly q) r -> e (p * q) r === e p r * e q r   , testProperty "eval x p = p" $     \p -> e X p === p   , testProperty "eval (monomial 0 c) p = c" $     \c p -> e (monomial 0 c) p === c    , testProperty "eval' (p + q) r = eval' p r + eval' q r" $-    \p q r -> e' (p + q) r === e' p r + e' q r+    \(ShortPoly p) (ShortPoly q) r -> e' (p + q) r === e' p r + e' q r   , testProperty "eval' (p * q) r = eval' p r * eval' q r" $-    \p q r -> e' (p * q) r === e' p r * e' q r+    \(ShortPoly p) (ShortPoly q) r -> e' (p * q) r === e' p r * e' q r   , testProperty "eval' x p = p" $     \p -> e' S.X p === p   , testProperty "eval' (S.monomial 0 c) p = c" $@@ -209,14 +212,14 @@ substTestGroup _ =   [ tenTimesLess $ tenTimesLess $ tenTimesLess $     testProperty "subst (p + q) r = subst p r + subst q r" $-    \p q r -> e (p + q) r === e p r + e q r+    \p q (ShortPoly r) -> e (p + q) r === e p r + e q r   , testProperty "subst x p = p" $     \p -> e X p === p   , testProperty "subst (monomial 0 c) p = monomial 0 c" $     \c p -> e (monomial 0 c) p === monomial 0 c   , tenTimesLess $ tenTimesLess $ tenTimesLess $     testProperty "subst' (p + q) r = subst' p r + subst' q r" $-    \p q r -> e' (p + q) r === e' p r + e' q r+    \p q (ShortPoly r) -> e' (p + q) r === e' p r + e' q r   , testProperty "subst' x p = p" $     \p -> e' S.X p === p   , testProperty "subst' (S.monomial 0 c) p = S.monomial 0 c" $@@ -231,21 +234,49 @@ derivTests :: TestTree derivTests = testGroup "deriv"   [ testProperty "deriv = S.deriv" $-    \(p :: Poly V.Vector Integer) -> deriv p === S.deriv p+    \(p :: VPoly Integer) -> deriv p === S.deriv p   , testProperty "integral = S.integral" $-    \(p :: Poly V.Vector Rational) -> integral p === S.integral p+    \(p :: VPoly Rational) -> integral p === S.integral p   , testProperty "deriv . integral = id" $-    \(p :: Poly V.Vector Rational) -> deriv (integral p) === p+    \(p :: VPoly Rational) -> deriv (integral p) === p   , testProperty "deriv c = 0" $-    \c -> deriv (monomial 0 c :: Poly V.Vector Int) === 0+    \c -> deriv (monomial 0 c :: UPoly Int) === 0   , testProperty "deriv cX = c" $-    \c -> deriv (monomial 0 c * X :: Poly V.Vector Int) === monomial 0 c+    \c -> deriv (monomial 0 c * X :: UPoly Int) === monomial 0 c   , testProperty "deriv (p + q) = deriv p + deriv q" $-    \p q -> deriv (p + q) === (deriv p + deriv q :: Poly V.Vector Int)+    \p q -> deriv (p + q) === (deriv p + deriv q :: UPoly Int)   , testProperty "deriv (p * q) = p * deriv q + q * deriv p" $-    \p q -> deriv (p * q) === (p * deriv q + q * deriv p :: Poly V.Vector Int)+    \p q -> deriv (p * q) === (p * deriv q + q * deriv p :: UPoly Int)   , tenTimesLess $ tenTimesLess $ tenTimesLess $     testProperty "deriv (subst p q) = deriv q * subst (deriv p) q" $-    \(p :: Poly V.Vector Int) (q :: Poly U.Vector Int) ->+    \(ShortPoly (p :: UPoly Int)) (ShortPoly (q :: UPoly Int)) ->       deriv (subst p q) === deriv q * subst (deriv p) q+  ]++patternTests :: TestTree+patternTests = testGroup "pattern"+  [ testProperty "X  :: UPoly Int" $ once $+    case (monomial 1 1 :: UPoly Int) of X -> True; _ -> False+  , testProperty "X  :: UPoly Int" $ once $+    (X :: UPoly Int) === monomial 1 1+  , testProperty "X' :: UPoly Int" $ once $+    case (S.monomial 1 1 :: UPoly Int) of S.X -> True; _ -> False+  , testProperty "X' :: UPoly Int" $ once $+    (S.X :: UPoly Int) === S.monomial 1 1+  , testProperty "X' :: UPoly ()" $ once $+    case (zero :: UPoly ()) of S.X -> True; _ -> False+  , testProperty "X' :: UPoly ()" $ once $+    (S.X :: UPoly ()) === zero+  ]++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" $+    \(xs :: UPoly Int8) -> xs === toPoly (unPoly xs)+  , testProperty "S.toPoly . S.unPoly = id" $+    \(xs :: UPoly Int8) -> xs === S.toPoly (S.unPoly xs)   ]

@@ -1,52 +1,35 @@ {-# LANGUAGE FlexibleContexts           #-} {-# LANGUAGE FlexibleInstances          #-}-{-# LANGUAGE GeneralizedNewtypeDeriving #-} {-# LANGUAGE ScopedTypeVariables        #-} -{-# OPTIONS_GHC -fno-warn-orphans #-}- module DenseLaurent   ( testSuite   ) where -import Prelude hiding (gcd, quotRem, rem)-import Data.Euclidean (Euclidean(..), GcdDomain, Field)+import Prelude hiding (gcd, quotRem, quot, rem)+import Control.Exception+import Data.Euclidean (GcdDomain(..), Field) import Data.Int import qualified Data.Poly import Data.Poly.Laurent import Data.Proxy import Data.Semiring (Semiring(..))-import qualified Data.Vector as V import qualified Data.Vector.Generic as G import qualified Data.Vector.Unboxed as U import Test.Tasty import Test.Tasty.QuickCheck hiding (scale, numTests) -import Dense (ShortPoly(..)) import Quaternion import TestUtils -instance (Eq a, Semiring a, Arbitrary a, G.Vector v a) => Arbitrary (Laurent v a) where-  arbitrary = toLaurent <$> ((`rem` 10) <$> arbitrary) <*> arbitrary-  shrink = fmap (uncurry toLaurent) . shrink . unLaurent--instance (Eq a, Semiring a, Arbitrary a, G.Vector v a) => Arbitrary (LaurentOverField (Laurent v a)) where-  arbitrary = (LaurentOverField .) . toLaurent <$> ((`rem` 10) <$> arbitrary) <*> (Data.Poly.unPolyOverField <$> arbitrary)-  shrink = fmap (LaurentOverField . uncurry toLaurent . fmap Data.Poly.unPolyOverField) . shrink . fmap Data.Poly.PolyOverField . unLaurent . unLaurentOverField--newtype ShortLaurent a = ShortLaurent { unShortLaurent :: a }-  deriving (Eq, Show, Semiring, GcdDomain)--instance (Eq a, Semiring a, Arbitrary a, G.Vector v a) => Arbitrary (ShortLaurent (Laurent v a)) where-  arbitrary = (ShortLaurent .) . toLaurent <$> ((`rem` 10) <$> arbitrary) <*> (unShortPoly <$> arbitrary)-  shrink = fmap (ShortLaurent . uncurry toLaurent . fmap unShortPoly) . shrink . fmap ShortPoly . unLaurent . unShortLaurent- testSuite :: TestTree testSuite = testGroup "DenseLaurent"   [ otherTests+  , divideByZeroTests   , lawsTests   , evalTests   , derivTests+  , patternTests   ]  lawsTests :: TestTree@@ -55,39 +38,39 @@  semiringTests :: [TestTree] semiringTests =-  [ mySemiringLaws (Proxy :: Proxy (Laurent U.Vector ()))-  , mySemiringLaws (Proxy :: Proxy (Laurent U.Vector Int8))-  , mySemiringLaws (Proxy :: Proxy (Laurent V.Vector Integer))-  , mySemiringLaws (Proxy :: Proxy (Laurent U.Vector (Quaternion Int)))+  [ mySemiringLaws (Proxy :: Proxy (ULaurent ()))+  , mySemiringLaws (Proxy :: Proxy (ULaurent Int8))+  , mySemiringLaws (Proxy :: Proxy (VLaurent Integer))+  , mySemiringLaws (Proxy :: Proxy (ULaurent (Quaternion Int)))   ]  ringTests :: [TestTree] ringTests =-  [ myRingLaws (Proxy :: Proxy (Laurent U.Vector ()))-  , myRingLaws (Proxy :: Proxy (Laurent U.Vector Int8))-  , myRingLaws (Proxy :: Proxy (Laurent V.Vector Integer))-  , myRingLaws (Proxy :: Proxy (Laurent U.Vector (Quaternion Int)))+  [ myRingLaws (Proxy :: Proxy (ULaurent ()))+  , myRingLaws (Proxy :: Proxy (ULaurent Int8))+  , myRingLaws (Proxy :: Proxy (VLaurent Integer))+  , myRingLaws (Proxy :: Proxy (ULaurent (Quaternion Int)))   ]  numTests :: [TestTree] numTests =-  [ myNumLaws (Proxy :: Proxy (Laurent U.Vector Int8))-  , myNumLaws (Proxy :: Proxy (Laurent V.Vector Integer))-  , myNumLaws (Proxy :: Proxy (Laurent U.Vector (Quaternion Int)))+  [ myNumLaws (Proxy :: Proxy (ULaurent Int8))+  , myNumLaws (Proxy :: Proxy (VLaurent Integer))+  , myNumLaws (Proxy :: Proxy (ULaurent (Quaternion Int)))   ]  gcdDomainTests :: [TestTree] gcdDomainTests =-  [ myGcdDomainLaws (Proxy :: Proxy (ShortLaurent (Laurent V.Vector Integer)))-  , myGcdDomainLaws (Proxy :: Proxy (LaurentOverField (Laurent V.Vector Rational)))+  [ myGcdDomainLaws (Proxy :: Proxy (ShortPoly (VLaurent Integer)))+  , myGcdDomainLaws (Proxy :: Proxy (ShortPoly (VLaurent Rational)))   ]  showTests :: [TestTree] showTests =-  [ myShowLaws (Proxy :: Proxy (Laurent U.Vector ()))-  , myShowLaws (Proxy :: Proxy (Laurent U.Vector Int8))-  , myShowLaws (Proxy :: Proxy (Laurent V.Vector Integer))-  , myShowLaws (Proxy :: Proxy (Laurent U.Vector (Quaternion Int)))+  [ myShowLaws (Proxy :: Proxy (ULaurent ()))+  , myShowLaws (Proxy :: Proxy (ULaurent Int8))+  , myShowLaws (Proxy :: Proxy (VLaurent Integer))+  , myShowLaws (Proxy :: Proxy (ULaurent (Quaternion Int)))   ]  otherTests :: TestTree@@ -109,12 +92,22 @@   , tenTimesLess $     testProperty "scale matches multiplication by monomial" $     \p c (xs :: ULaurent a) -> scale p c xs === monomial p c * xs+  , tenTimesLess $+    testProperty "toLaurent . unLaurent" $+    \(xs :: ULaurent a) -> uncurry toLaurent (unLaurent xs) === xs   ] +divideByZeroTests :: TestTree+divideByZeroTests = testGroup "divideByZero"+  [ testProperty "divide"  $ testProp divide+  ]+  where+    testProp f xs = ioProperty ((== Left DivideByZero) <$> try (evaluate (xs `f` (0 :: VLaurent Rational))))+ evalTests :: TestTree evalTests = testGroup "eval" $ concat-  [ evalTestGroup  (Proxy :: Proxy (Laurent V.Vector Rational))-  , substTestGroup (Proxy :: Proxy (Laurent U.Vector Int8))+  [ evalTestGroup  (Proxy :: Proxy (VLaurent Rational))+  , substTestGroup (Proxy :: Proxy (ULaurent Int8))   ]  evalTestGroup@@ -124,9 +117,9 @@   -> [TestTree] evalTestGroup _ =   [ testProperty "eval (p + q) r = eval p r + eval q r" $-    \p q r -> e (p `plus` q) r === e p r `plus` e q r+    \(ShortPoly p) (ShortPoly q) r -> e (p `plus` q) r === e p r `plus` e q r   , testProperty "eval (p * q) r = eval p r * eval q r" $-    \p q r -> e (p `times` q) r === e p r `times` e q r+    \(ShortPoly p) (ShortPoly q) r -> e (p `times` q) r === e p r `times` e q r   , testProperty "eval x p = p" $     \p -> e X p === p   , testProperty "eval (monomial 0 c) p = c" $@@ -144,7 +137,7 @@ substTestGroup _ =   [ tenTimesLess $ tenTimesLess $ tenTimesLess $     testProperty "subst (p + q) r = subst p r + subst q r" $-    \p q r -> e (p + q) r === e p r + e q r+    \p q (ShortPoly r) -> e (p + q) r === e p r + e q r   , testProperty "subst x p = p" $     \p -> e Data.Poly.X p === p   , testProperty "subst (monomial 0 c) p = monomial 0 c" $@@ -157,15 +150,35 @@ derivTests :: TestTree derivTests = testGroup "deriv"   [ testProperty "deriv c = 0" $-    \c -> deriv (monomial 0 c :: Laurent V.Vector Int) === 0+    \c -> deriv (monomial 0 c :: ULaurent Int) === 0   , testProperty "deriv cX = c" $-    \c -> deriv (monomial 0 c * X :: Laurent V.Vector Int) === monomial 0 c+    \c -> deriv (monomial 0 c * X :: ULaurent Int) === monomial 0 c   , testProperty "deriv (p + q) = deriv p + deriv q" $-    \p q -> deriv (p + q) === (deriv p + deriv q :: Laurent V.Vector Int)+    \p q -> deriv (p + q) === (deriv p + deriv q :: ULaurent Int)   , testProperty "deriv (p * q) = p * deriv q + q * deriv p" $-    \p q -> deriv (p * q) === (p * deriv q + q * deriv p :: Laurent V.Vector Int)+    \p q -> deriv (p * q) === (p * deriv q + q * deriv p :: ULaurent Int)   , tenTimesLess $ tenTimesLess $ tenTimesLess $     testProperty "deriv (subst p q) = deriv q * subst (deriv p) q" $-    \(p :: Data.Poly.Poly V.Vector Int) (q :: Laurent U.Vector Int) ->+    \(ShortPoly (p :: Data.Poly.UPoly Int)) (ShortPoly (q :: ULaurent Int)) ->       deriv (subst p q) === deriv q * subst (Data.Poly.deriv p) q+  ]++patternTests :: TestTree+patternTests = testGroup "pattern"+  [ testProperty "X  :: ULaurent Int" $ once $+    case (monomial 1 1 :: ULaurent Int) of X -> True; _ -> False+  , testProperty "X  :: ULaurent Int" $ once $+    (X :: ULaurent Int) === monomial 1 1+  , testProperty "X :: ULaurent ()" $ once $+    case (zero :: ULaurent ()) of X -> True; _ -> False+  , testProperty "X :: ULaurent ()" $ once $+    (X :: ULaurent ()) === zero+  , testProperty "X^-k" $+    \(NonNegative j) k -> ((X^j)^-k :: ULaurent Int) === monomial (- j * k) 1+  , testProperty "^-" $+    \(p :: ULaurent Int) (NonNegative k) -> ioProperty $ do+      et <- try (evaluate (p^-k)) :: IO (Either PatternMatchFail (ULaurent Int))+      pure $ case et of+        Left{}  -> True+        Right t -> p^k * t == one   ]

@@ -2,17 +2,23 @@  import Test.Tasty -import qualified Dense as Dense-import qualified DenseLaurent as DenseLaurent-import qualified Orthogonal as Orthogonal-import qualified Sparse as Sparse-import qualified SparseLaurent as SparseLaurent+import qualified Dense+import qualified DenseLaurent+import qualified DFT+import qualified Multi+import qualified MultiLaurent+import qualified Orthogonal+import qualified Sparse+import qualified SparseLaurent  main :: IO () main = defaultMain $ testGroup "All"     [ Dense.testSuite     , DenseLaurent.testSuite+    , DFT.testSuite     , Sparse.testSuite     , SparseLaurent.testSuite+    , Multi.testSuite+    , MultiLaurent.testSuite     , Orthogonal.testSuite     ]

@@ -0,0 +1,307 @@+{-# LANGUAGE DataKinds                  #-}+{-# LANGUAGE FlexibleContexts           #-}+{-# LANGUAGE FlexibleInstances          #-}+{-# LANGUAGE ScopedTypeVariables        #-}+{-# LANGUAGE UndecidableInstances       #-}++module Multi+  ( testSuite+  ) where++import Prelude hiding (gcd, quotRem, rem)+import Control.Exception+import Data.Euclidean (GcdDomain(..))+import Data.Function+import Data.Int+import Data.List (groupBy, sortOn)+import Data.Mod.Word+import Data.Proxy+import Data.Semiring (Semiring(..))+import qualified Data.Vector as V+import qualified Data.Vector.Generic as G+import qualified Data.Vector.Generic.Sized as SG+import qualified Data.Vector.Sized as SV+import qualified Data.Vector.Unboxed as U+import qualified Data.Vector.Unboxed.Sized as SU+import Test.Tasty+import Test.Tasty.QuickCheck hiding (scale, numTests)++import Data.Poly.Multi+import qualified Data.Poly.Multi.Semiring as S++import Quaternion+import TestUtils++testSuite :: TestTree+testSuite = testGroup "Multi"+  [ arithmeticTests+  , otherTests+  , divideByZeroTests+  , lawsTests+  , evalTests+  , derivTests+  , patternTests+  , conversionTests+  ]++lawsTests :: TestTree+lawsTests = testGroup "Laws"+  $ semiringTests ++ ringTests ++ numTests ++ gcdDomainTests ++ isListTests ++ showTests++semiringTests :: [TestTree]+semiringTests =+  [ mySemiringLaws (Proxy :: Proxy (UMultiPoly 3 ()))+  , mySemiringLaws (Proxy :: Proxy (ShortPoly (UMultiPoly 2 Int8)))+  , mySemiringLaws (Proxy :: Proxy (ShortPoly (VMultiPoly 2 Integer)))+  , tenTimesLess+  $ mySemiringLaws (Proxy :: Proxy (ShortPoly (UMultiPoly 2 (Quaternion Int))))+  ]++ringTests :: [TestTree]+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)))+  ]++numTests :: [TestTree]+numTests =+  [ myNumLaws (Proxy :: Proxy (ShortPoly (UMultiPoly 2 Int8)))+  , myNumLaws (Proxy :: Proxy (ShortPoly (VMultiPoly 2 Integer)))+  , tenTimesLess+  $ myNumLaws (Proxy :: Proxy (ShortPoly (UMultiPoly 2 (Quaternion Int))))+  ]++gcdDomainTests :: [TestTree]+gcdDomainTests =+  [ myGcdDomainLaws (Proxy :: Proxy (ShortPoly (VMultiPoly 3 Integer)))+  , tenTimesLess+  $ myGcdDomainLaws (Proxy :: Proxy (ShortPoly (VMultiPoly 3 (Mod 3))))+  , tenTimesLess+  $ myGcdDomainLaws (Proxy :: Proxy (ShortPoly (VMultiPoly 3 Rational)))+  ]++isListTests :: [TestTree]+isListTests =+  [ myIsListLaws (Proxy :: Proxy (UMultiPoly 3 ()))+  , myIsListLaws (Proxy :: Proxy (UMultiPoly 3 Int8))+  , myIsListLaws (Proxy :: Proxy (VMultiPoly 3 Integer))+  , tenTimesLess+  $ myIsListLaws (Proxy :: Proxy (UMultiPoly 3 (Quaternion Int)))+  ]++showTests :: [TestTree]+showTests =+  [ myShowLaws (Proxy :: Proxy (UMultiPoly 4 ()))+  , myShowLaws (Proxy :: Proxy (UMultiPoly 4 Int8))+  , myShowLaws (Proxy :: Proxy (VMultiPoly 4 Integer))+  , tenTimesLess+  $ myShowLaws (Proxy :: Proxy (UMultiPoly 4 (Quaternion Int)))+  ]++arithmeticTests :: TestTree+arithmeticTests = testGroup "Arithmetic"+  [ testProperty "addition matches reference" $+    \(xs :: [(SU.Vector 3 Word, Int)]) ys -> toMultiPoly (V.fromList (addRef xs ys)) ===+      toMultiPoly (V.fromList xs) + toMultiPoly (V.fromList ys)+  , testProperty "subtraction matches reference" $+    \(xs :: [(SU.Vector 3 Word, Int)]) ys -> toMultiPoly (V.fromList (subRef xs ys)) ===+      toMultiPoly (V.fromList xs) - toMultiPoly (V.fromList ys)+  , tenTimesLess $+    testProperty "multiplication matches reference" $+    \(xs :: [(SU.Vector 3 Word, Int)]) ys -> toMultiPoly (V.fromList (mulRef xs ys)) ===+      toMultiPoly (V.fromList xs) * toMultiPoly (V.fromList ys)+  ]++addRef :: (Num a, Ord t) => [(t, a)] -> [(t, a)] -> [(t, a)]+addRef [] ys = ys+addRef xs [] = xs+addRef xs@((xp, xc) : xs') ys@((yp, yc) : ys') =+  case xp `compare` yp of+    LT -> (xp, xc) : addRef xs' ys+    EQ -> (xp, xc + yc) : addRef xs' ys'+    GT -> (yp, yc) : addRef xs ys'++subRef :: (Num a, Ord t) => [(t, a)] -> [(t, a)] -> [(t, a)]+subRef [] ys = map (fmap negate) ys+subRef xs [] = xs+subRef xs@((xp, xc) : xs') ys@((yp, yc) : ys') =+  case xp `compare` yp of+    LT -> (xp, xc) : subRef xs' ys+    EQ -> (xp, xc - yc) : subRef xs' ys'+    GT -> (yp, negate yc) : subRef xs ys'++mulRef :: (Num a, Ord t, Num t) => [(t, a)] -> [(t, a)] -> [(t, a)]+mulRef xs ys+  = map (\ws -> (fst (head ws), sum (map snd ws)))+  $ groupBy ((==) `on` fst)+  $ sortOn fst+  $ [ (xp + yp, xc * yc) | (xp, xc) <- xs, (yp, yc) <- ys ]++otherTests :: TestTree+otherTests = testGroup "other" $ concat+  [ otherTestGroup (Proxy :: Proxy Int8)+  , otherTestGroup (Proxy :: Proxy (Quaternion Int))+  ]++otherTestGroup+  :: forall a.+     (Eq a, Show a, Semiring a, Num a, Arbitrary a, U.Unbox a, G.Vector U.Vector a)+  => Proxy a+  -> [TestTree]+otherTestGroup _ =+  [ testProperty "monomial matches reference" $+    \(ps :: SU.Vector 3 Word) (c :: a) -> monomial ps c === toMultiPoly (V.fromList (monomialRef ps c))+  , tenTimesLess $+    testProperty "scale matches multiplication by monomial" $+    \ps c (xs :: UMultiPoly 3 a) -> scale ps c xs === monomial ps c * xs+  , tenTimesLess $+    testProperty "scale' matches multiplication by monomial" $+    \ps c (xs :: UMultiPoly 3 a) -> S.scale ps c xs === S.monomial ps c * xs+  ]++monomialRef :: Num a => t -> a -> [(t, a)]+monomialRef p c = [(p, c)]++divideByZeroTests :: TestTree+divideByZeroTests = testGroup "divideByZero"+  [ testProperty "divide"  $ testProp divide+  ]+  where+    testProp f xs = ioProperty ((== Left DivideByZero) <$> try (evaluate (xs `f` (0 :: VMultiPoly 3 Rational))))++evalTests :: TestTree+evalTests = testGroup "eval" $ concat+  [ evalTestGroup  (Proxy :: Proxy (UMultiPoly 3 Int8))+  , evalTestGroup  (Proxy :: Proxy (VMultiPoly 3 Integer))+  , substTestGroup (Proxy :: Proxy (UMultiPoly 3 Int8))+  ]++evalTestGroup+  :: forall v a.+     (Eq a, Num a, Semiring a, Arbitrary a, Show a, Eq (v (SU.Vector 3 Word, a)), Show (v (SU.Vector 3 Word, a)), G.Vector v (SU.Vector 3 Word, a))+  => Proxy (MultiPoly v 3 a)+  -> [TestTree]+evalTestGroup _ =+  [ testProperty "eval (p + q) rs = eval p rs + eval q rs" $+    \(ShortPoly p) (ShortPoly q) rs -> e (p + q) rs === e p rs + e q rs+  , testProperty "eval (p * q) rs = eval p rs * eval q rs" $+    \(ShortPoly p) (ShortPoly q) rs -> e (p * q) rs === e p rs * e q rs+  , testProperty "eval x p = p" $+    \p -> e X (SV.fromTuple (p, undefined, undefined)) === p+  , testProperty "eval (monomial 0 c) p = c" $+    \c ps -> e (monomial 0 c) ps === c++  , testProperty "eval' (p + q) rs = eval' p rs + eval' q rs" $+    \(ShortPoly p) (ShortPoly q) rs -> e' (p + q) rs === e' p rs + e' q rs+  , testProperty "eval' (p * q) rs = eval' p rs * eval' q rs" $+    \(ShortPoly p) (ShortPoly q) rs -> e' (p * q) rs === e' p rs * e' q rs+  , testProperty "eval' x p = p" $+    \p -> e' S.X (SV.fromTuple (p, undefined, undefined)) === p+  , testProperty "eval' (monomial 0 c) p = c" $+    \c ps -> e' (monomial 0 c) ps === c+  ]+  where+    e :: MultiPoly v 3 a -> SV.Vector 3 a -> a+    e = eval+    e' :: MultiPoly v 3 a -> SV.Vector 3 a -> a+    e' = S.eval++substTestGroup+  :: forall v a.+     (Eq a, Num a, Semiring a, Arbitrary a, Show a, Eq (v (SU.Vector 3 Word, a)), Show (v (SU.Vector 3 Word, a)), G.Vector v (SU.Vector 3 Word, a))+  => Proxy (MultiPoly v 3 a)+  -> [TestTree]+substTestGroup _ =+  [ testProperty "subst x p = p" $+    \p -> e X (SV.fromTuple (p, undefined, undefined)) === p+  , testProperty "subst (monomial 0 c) ps = monomial 0 c" $+    \c ps -> e (monomial 0 c) ps === monomial 0 c+  , testProperty "subst' x p = p" $+    \p -> e' S.X (SV.fromTuple (p, undefined, undefined)) === p+  , testProperty "subst' (S.monomial 0 c) ps = S.monomial 0 c" $+    \c ps -> e' (S.monomial 0 c) ps === S.monomial 0 c+  ]+  where+    e :: MultiPoly v 3 a -> SV.Vector 3 (MultiPoly v 3 a) -> MultiPoly v 3 a+    e = subst+    e' :: MultiPoly v 3 a -> SV.Vector 3 (MultiPoly v 3 a) -> MultiPoly v 3 a+    e' = S.subst++derivTests :: TestTree+derivTests = testGroup "deriv"+  [ testProperty "deriv = S.deriv" $+    \k (p :: VMultiPoly 3 Integer) -> deriv k p === S.deriv k p+  , testProperty "integral = S.integral" $+    \k (p :: VMultiPoly 3 Rational) -> integral k p === S.integral k p+  , testProperty "deriv . integral = id" $+    \k (p :: VMultiPoly 3 Rational) ->+      deriv k (integral k p) === p+  , testProperty "deriv c = 0" $+    \k c ->+      deriv k (monomial 0 c :: UMultiPoly 3 Int) === 0+  , testProperty "deriv cX = c" $+    \c ->+      deriv 0 (monomial 0 c * X :: UMultiPoly 3 Int) === monomial 0 c+  , testProperty "deriv (p + q) = deriv p + deriv q" $+    \k p q ->+      deriv k (p + q) === (deriv k p + deriv k q :: UMultiPoly 3 Int)+  , testProperty "deriv (p * q) = p * deriv q + q * deriv p" $+    \k p q ->+      deriv k (p * q) === (p * deriv k q + q * deriv k p :: UMultiPoly 3 Int)+  ]++patternTests :: TestTree+patternTests = testGroup "pattern"+  [ testProperty "X  :: UMultiPoly Int" $ once $+    case (monomial 1 1 :: UMultiPoly 1 Int) of X -> True; _ -> False+  , testProperty "X  :: UMultiPoly Int" $ once $+    (X :: UMultiPoly 1 Int) === monomial 1 1+  , testProperty "S.X  :: UMultiPoly Int8" $ once $+    case (S.monomial 1 1 :: UMultiPoly 1 Int8) of S.X -> True; _ -> False+  , testProperty "S.X  :: UMultiPoly Int8" $ once $+    (S.X :: UMultiPoly 1 Int8) === S.monomial 1 1+  , testProperty "X :: UMultiPoly ()" $ once $+    case (zero :: UMultiPoly 1 ()) of S.X -> True; _ -> False+  , testProperty "X :: UMultiPoly ()" $ once $+    (S.X :: UMultiPoly 1 ()) === zero++  , testProperty "Y  :: UMultiPoly Int" $ once $+    case (monomial (SG.fromTuple (0, 1)) 1 :: UMultiPoly 2 Int) of Y -> True; _ -> False+  , testProperty "Y  :: UMultiPoly Int" $ once $+    (Y :: UMultiPoly 2 Int) === monomial (SG.fromTuple (0, 1)) 1+  , testProperty "S.Y  :: UMultiPoly Int8" $ once $+    case (S.monomial (SG.fromTuple (0, 1)) 1 :: UMultiPoly 2 Int8) of S.Y -> True; _ -> False+  , testProperty "S.Y  :: UMultiPoly Int8" $ once $+    (S.Y :: UMultiPoly 2 Int8) === S.monomial (SG.fromTuple (0, 1)) 1+  , testProperty "Y :: UMultiPoly ()" $ once $+    case (zero :: UMultiPoly 2 ()) of S.Y -> True; _ -> False+  , testProperty "Y :: UMultiPoly ()" $ once $+    (S.Y :: UMultiPoly 2 ()) === zero++  , testProperty "Z  :: UMultiPoly Int" $ once $+    case (monomial (SG.fromTuple (0, 0, 1)) 1 :: UMultiPoly 3 Int) of Z -> True; _ -> False+  , testProperty "Z  :: UMultiPoly Int" $ once $+    (Z :: UMultiPoly 3 Int) === monomial (SG.fromTuple (0, 0, 1)) 1+  , testProperty "S.Z  :: UMultiPoly Int8" $ once $+    case (S.monomial (SG.fromTuple (0, 0, 1)) 1 :: UMultiPoly 3 Int) of S.Z -> True; _ -> False+  , testProperty "S.Z  :: UMultiPoly Int" $ once $+    (S.Z :: UMultiPoly 3 Int) === S.monomial (SG.fromTuple (0, 0, 1)) 1+  , testProperty "Z :: UMultiPoly ()" $ once $+    case (zero :: UMultiPoly 3 ()) of S.Z -> True; _ -> False+  , testProperty "Z :: UMultiPoly ()" $ once $+    (S.Z :: UMultiPoly 3 ()) === zero+  ]++conversionTests :: TestTree+conversionTests = testGroup "conversions"+  [ testProperty "unsegregate . segregate = id" $+    \(xs :: UMultiPoly 3 Int8) -> xs === unsegregate (segregate xs)+  , testProperty "segregate . unsegregate = id" $+    \xs -> xs === segregate (unsegregate xs :: UMultiPoly 3 Int8)+  , testProperty "toMultiPoly . unMultiPoly = id" $+    \(xs :: UMultiPoly 3 Int8) -> xs === toMultiPoly (unMultiPoly xs)+  , testProperty "S.toMultiPoly . S.unMultiPoly = id" $+    \(xs :: UMultiPoly 3 Int8) -> xs === S.toMultiPoly (S.unMultiPoly xs)+  ]

@@ -0,0 +1,222 @@+{-# LANGUAGE DataKinds                  #-}+{-# LANGUAGE FlexibleContexts           #-}+{-# LANGUAGE FlexibleInstances          #-}+{-# LANGUAGE ScopedTypeVariables        #-}+{-# LANGUAGE UndecidableInstances       #-}++module MultiLaurent+  ( testSuite+  ) where++import Prelude hiding (gcd, quotRem, quot, rem)+import Control.Exception+import Data.Euclidean (GcdDomain(..), Field)+import Data.Int+import qualified Data.Poly.Multi+import Data.Poly.Multi.Laurent+import Data.Proxy+import Data.Semiring (Semiring(..))+import qualified Data.Vector.Generic as G+import qualified Data.Vector.Generic.Sized as SG+import qualified Data.Vector.Sized as SV+import qualified Data.Vector.Unboxed as U+import qualified Data.Vector.Unboxed.Sized as SU+import Test.Tasty+import Test.Tasty.QuickCheck hiding (scale, numTests)++import Quaternion+import TestUtils++testSuite :: TestTree+testSuite = testGroup "MultiLaurent"+  [ otherTests+  , divideByZeroTests+  , lawsTests+  , evalTests+  , derivTests+  , patternTests+  , conversionTests+  ]++lawsTests :: TestTree+lawsTests = testGroup "Laws"+  $ semiringTests ++ ringTests ++ numTests ++ gcdDomainTests ++ isListTests ++ showTests++semiringTests :: [TestTree]+semiringTests =+  [ mySemiringLaws (Proxy :: Proxy (UMultiLaurent 3 ()))+  , mySemiringLaws (Proxy :: Proxy (ShortPoly (UMultiLaurent 2 Int8)))+  , mySemiringLaws (Proxy :: Proxy (ShortPoly (VMultiLaurent 2 Integer)))+  , tenTimesLess+  $ mySemiringLaws (Proxy :: Proxy (ShortPoly (UMultiLaurent 2 (Quaternion Int))))+  ]++ringTests :: [TestTree]+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)))+  ]++numTests :: [TestTree]+numTests =+  [ myNumLaws (Proxy :: Proxy (ShortPoly (UMultiLaurent 2 Int8)))+  , myNumLaws (Proxy :: Proxy (ShortPoly (VMultiLaurent 2 Integer)))+  , tenTimesLess+  $ myNumLaws (Proxy :: Proxy (ShortPoly (UMultiLaurent 2 (Quaternion Int))))+  ]++gcdDomainTests :: [TestTree]+gcdDomainTests =+  [ myGcdDomainLaws (Proxy :: Proxy (ShortPoly (VMultiLaurent 3 Integer)))+  , tenTimesLess+  $ myGcdDomainLaws (Proxy :: Proxy (ShortPoly (VMultiLaurent 3 Rational)))+  ]++isListTests :: [TestTree]+isListTests =+  [ myIsListLaws (Proxy :: Proxy (UMultiLaurent 3 ()))+  , myIsListLaws (Proxy :: Proxy (UMultiLaurent 3 Int8))+  , myIsListLaws (Proxy :: Proxy (VMultiLaurent 3 Integer))+  , tenTimesLess+  $ myIsListLaws (Proxy :: Proxy (UMultiLaurent 3 (Quaternion Int)))+  ]++showTests :: [TestTree]+showTests =+  [ myShowLaws (Proxy :: Proxy (UMultiLaurent 4 ()))+  , myShowLaws (Proxy :: Proxy (UMultiLaurent 4 Int8))+  , myShowLaws (Proxy :: Proxy (VMultiLaurent 4 Integer))+  , tenTimesLess+  $ myShowLaws (Proxy :: Proxy (UMultiLaurent 4 (Quaternion Int)))+  ]++otherTests :: TestTree+otherTests = testGroup "other" $ concat+  [ otherTestGroup (Proxy :: Proxy Int8)+  , otherTestGroup (Proxy :: Proxy (Quaternion Int))+  ]++otherTestGroup+  :: forall a.+     (Eq a, Show a, Semiring a, Num a, Arbitrary a, U.Unbox a, G.Vector U.Vector a)+  => Proxy a+  -> [TestTree]+otherTestGroup _ =+  [ testProperty "scale matches multiplication by monomial" $+    \p c (xs :: UMultiLaurent 3 a) -> scale p c xs === monomial p c * xs+  , tenTimesLess $+    testProperty "toMultiLaurent . unMultiLaurent" $+    \(xs :: UMultiLaurent 3 a) -> uncurry toMultiLaurent (unMultiLaurent xs) === xs+  ]++divideByZeroTests :: TestTree+divideByZeroTests = testGroup "divideByZero"+  [ testProperty "divide"  $ testProp divide+  ]+  where+    testProp f xs = ioProperty ((== Left DivideByZero) <$> try (evaluate (xs `f` (0 :: VMultiLaurent 3 Rational))))++evalTests :: TestTree+evalTests = testGroup "eval" $ concat+  [ evalTestGroup  (Proxy :: Proxy (VMultiLaurent 3 Rational))+  , substTestGroup (Proxy :: Proxy (UMultiLaurent 3 Int8))+  ]++evalTestGroup+  :: forall v a.+     (Eq a, Field a, Arbitrary a, Show a, Eq (v (Word, a)), Show (v (Word, a)), G.Vector v (Word, a), Eq (v (SU.Vector 3 Word, a)), Show (v (SU.Vector 3 Word, a)), G.Vector v (SU.Vector 3 Word, a))+  => Proxy (MultiLaurent v 3 a)+  -> [TestTree]+evalTestGroup _ =+  [ testProperty "eval (p + q) r = eval p r + eval q r" $+    \(ShortPoly p) (ShortPoly q) r -> e (p `plus` q) r === e p r `plus` e q r+  , testProperty "eval (p * q) r = eval p r * eval q r" $+    \(ShortPoly p) (ShortPoly q) r -> e (p `times` q) r === e p r `times` e q r+  , testProperty "eval x p = p" $+    \p -> e X (SV.fromTuple (p, undefined, undefined)) === p+  , testProperty "eval (monomial 0 c) p = c" $+    \c p -> e (monomial 0 c) p === c+  ]+  where+    e :: MultiLaurent v 3 a -> SV.Vector 3 a -> a+    e = eval++substTestGroup+  :: forall v a.+     (Eq a, Num a, Semiring a, Arbitrary a, Show a, Eq (v (SU.Vector 3 Word, a)), Show (v (Word, a)), G.Vector v (Word, a), G.Vector v (SU.Vector 3 Word, a))+  => Proxy (MultiLaurent v 3 a)+  -> [TestTree]+substTestGroup _ =+  [ testProperty "subst x p = p" $+    \p -> e Data.Poly.Multi.X (SV.fromTuple (p, undefined, undefined)) === p+  , testProperty "subst (monomial 0 c) p = monomial 0 c" $+    \c p -> e (Data.Poly.Multi.monomial 0 c) p === monomial 0 c+  ]+  where+    e :: Data.Poly.Multi.MultiPoly v 3 a -> SV.Vector 3 (MultiLaurent v 3 a) -> MultiLaurent v 3 a+    e = subst++derivTests :: TestTree+derivTests = testGroup "deriv"+  [ testProperty "deriv c = 0" $+    \k c -> deriv k (monomial 0 c :: UMultiLaurent 3 Int) === 0+  , testProperty "deriv cX = c" $+    \c -> deriv 0 (monomial 0 c * X :: UMultiLaurent 3 Int) === monomial 0 c+  , testProperty "deriv (p + q) = deriv p + deriv q" $+    \k p q -> deriv k (p + q) === (deriv k p + deriv k q :: UMultiLaurent 3 Int)+  , testProperty "deriv (p * q) = p * deriv q + q * deriv p" $+    \k p q -> deriv k (p * q) === (p * deriv k q + q * deriv k p :: UMultiLaurent 3 Int)+  ]++patternTests :: TestTree+patternTests = testGroup "pattern"+  [ testProperty "X  :: UMultiLaurent Int" $ once $+    case (monomial 1 1 :: UMultiLaurent 1 Int) of X -> True; _ -> False+  , testProperty "X  :: UMultiLaurent Int" $ once $+    (X :: UMultiLaurent 1 Int) === monomial 1 1+  , testProperty "X :: UMultiLaurent ()" $ once $+    case (zero :: UMultiLaurent 1 ()) of X -> True; _ -> False+  , testProperty "X :: UMultiLaurent ()" $ once $+    (X :: UMultiLaurent 1 ()) === zero++  , testProperty "Y  :: UMultiLaurent Int" $ once $+    case (monomial (SG.fromTuple (0, 1)) 1 :: UMultiLaurent 2 Int) of Y -> True; _ -> False+  , testProperty "Y  :: UMultiLaurent Int" $ once $+    (Y :: UMultiLaurent 2 Int) === monomial (SG.fromTuple (0, 1)) 1+  , testProperty "Y :: UMultiLaurent ()" $ once $+    case (zero :: UMultiLaurent 2 ()) of Y -> True; _ -> False+  , testProperty "Y :: UMultiLaurent ()" $ once $+    (Y :: UMultiLaurent 2 ()) === zero++  , testProperty "Z  :: UMultiLaurent Int" $ once $+    case (monomial (SG.fromTuple (0, 0, 1)) 1 :: UMultiLaurent 3 Int) of Z -> True; _ -> False+  , testProperty "Z  :: UMultiLaurent Int" $ once $+    (Z :: UMultiLaurent 3 Int) === monomial (SG.fromTuple (0, 0, 1)) 1+  , testProperty "Z :: UMultiLaurent ()" $ once $+    case (zero :: UMultiLaurent 3 ()) of Z -> True; _ -> False+  , testProperty "Z :: UMultiLaurent ()" $ once $+    (Z :: UMultiLaurent 3 ()) === zero++  , testProperty "X^-k" $+    \(NonNegative j) k -> ((X^j)^-k :: UMultiLaurent 1 Int) === monomial (SG.singleton (- j * k)) 1+  , testProperty "Y^-k" $+    \(NonNegative j) k -> ((Y^j)^-k :: UMultiLaurent 2 Int) === monomial (SG.fromTuple (0, - j * k)) 1+  , testProperty "Z^-k" $+    \(NonNegative j) k -> ((Z^j)^-k :: UMultiLaurent 3 Int) === monomial (SG.fromTuple (0, 0, - j * k)) 1+  , testProperty "^-" $+    \(p :: UMultiLaurent 3 Int) (NonNegative k) -> ioProperty $ do+      et <- try (evaluate (p^-k)) :: IO (Either PatternMatchFail (UMultiLaurent 3 Int))+      pure $ case et of+        Left{}  -> True+        Right t -> p^k * t == one+  ]++conversionTests :: TestTree+conversionTests = testGroup "conversions"+  [ testProperty "unsegregate . segregate = id" $+    \(xs :: UMultiLaurent 3 Int8) -> xs === unsegregate (segregate xs)+  , testProperty "segregate . unsegregate = id" $+    \xs -> xs === segregate (unsegregate xs :: UMultiLaurent 3 Int8)+  ]

@@ -130,7 +130,7 @@   [ integral11 (x * y) === 0 | (x : xs) <- tails polys, y <- xs ]   where     polys :: [VPoly Rational]-    polys = take limit $ legendre+    polys = take limit legendre  hermiteProbRef :: [VPoly Integer] hermiteProbRef = iterate (\he -> [0, 1] * he - deriv he) 1

@@ -18,16 +18,14 @@   ) where  import Prelude hiding (negate)-import Control.Monad import Data.Semiring (Semiring(..), Ring(..), minus) import GHC.Generics import Test.Tasty.QuickCheck hiding (scale) -import Data.Vector.Unboxed (Vector)+import Data.Vector.Unboxed (Vector, Unbox) import qualified Data.Vector.Generic as G import Data.Vector.Unboxed.Mutable (MVector) import qualified Data.Vector.Generic.Mutable as M-import Data.Vector.Unboxed (Unbox)  data Quaternion a = Quaternion !a !a !a !a   deriving (Eq, Ord, Show, Generic)@@ -63,9 +61,9 @@ newtype instance MVector s (Quaternion a) = MV_Quaternion (MVector s (a, a, a, a)) newtype instance Vector    (Quaternion a) = V_Quaternion  (Vector    (a, a, a, a)) -instance (Unbox a) => Unbox (Quaternion a)+instance Unbox a => Unbox (Quaternion a) -instance (Unbox a) => M.MVector MVector (Quaternion a) where+instance Unbox a => M.MVector MVector (Quaternion a) where   {-# INLINE basicLength #-}   {-# INLINE basicUnsafeSlice #-}   {-# INLINE basicOverlaps #-}@@ -81,30 +79,30 @@   basicLength (MV_Quaternion v) = M.basicLength v   basicUnsafeSlice i n (MV_Quaternion v) = MV_Quaternion $ M.basicUnsafeSlice i n v   basicOverlaps (MV_Quaternion v1) (MV_Quaternion v2) = M.basicOverlaps v1 v2-  basicUnsafeNew n = MV_Quaternion `liftM` M.basicUnsafeNew n+  basicUnsafeNew n = MV_Quaternion `fmap` M.basicUnsafeNew n   basicInitialize (MV_Quaternion v) = M.basicInitialize v-  basicUnsafeReplicate n (Quaternion a b c d) = MV_Quaternion `liftM` M.basicUnsafeReplicate n (a, b, c, d)-  basicUnsafeRead (MV_Quaternion v) i = (\(a, b, c, d) -> Quaternion a b c d) `liftM` M.basicUnsafeRead v i+  basicUnsafeReplicate n (Quaternion a b c d) = MV_Quaternion `fmap` M.basicUnsafeReplicate n (a, b, c, d)+  basicUnsafeRead (MV_Quaternion v) i = (\(a, b, c, d) -> Quaternion a b c d) `fmap` M.basicUnsafeRead v i   basicUnsafeWrite (MV_Quaternion v) i (Quaternion a b c d) = M.basicUnsafeWrite v i (a, b, c, d)   basicClear (MV_Quaternion v) = M.basicClear v   basicSet (MV_Quaternion v) (Quaternion a b c d) = M.basicSet v (a, b, c, d)   basicUnsafeCopy (MV_Quaternion v1) (MV_Quaternion v2) = M.basicUnsafeCopy v1 v2   basicUnsafeMove (MV_Quaternion v1) (MV_Quaternion v2) = M.basicUnsafeMove v1 v2-  basicUnsafeGrow (MV_Quaternion v) n = MV_Quaternion `liftM` M.basicUnsafeGrow v n+  basicUnsafeGrow (MV_Quaternion v) n = MV_Quaternion `fmap` M.basicUnsafeGrow v n -instance (Unbox a) => G.Vector Vector (Quaternion a) where+instance Unbox a => G.Vector Vector (Quaternion a) where   {-# INLINE basicUnsafeFreeze #-}   {-# INLINE basicUnsafeThaw #-}   {-# INLINE basicLength #-}   {-# INLINE basicUnsafeSlice #-}   {-# INLINE basicUnsafeIndexM #-}   {-# INLINE elemseq #-}-  basicUnsafeFreeze (MV_Quaternion v) = V_Quaternion `liftM` G.basicUnsafeFreeze v-  basicUnsafeThaw (V_Quaternion v) = MV_Quaternion `liftM` G.basicUnsafeThaw v+  basicUnsafeFreeze (MV_Quaternion v) = V_Quaternion `fmap` G.basicUnsafeFreeze v+  basicUnsafeThaw (V_Quaternion v) = MV_Quaternion `fmap` G.basicUnsafeThaw v   basicLength (V_Quaternion v) = G.basicLength v   basicUnsafeSlice i n (V_Quaternion v) = V_Quaternion $ G.basicUnsafeSlice i n v   basicUnsafeIndexM (V_Quaternion v) i-                = (\(a, b, c, d) -> Quaternion a b c d) `liftM` G.basicUnsafeIndexM v i+                = (\(a, b, c, d) -> Quaternion a b c d) `fmap` G.basicUnsafeIndexM v i   basicUnsafeCopy (MV_Quaternion mv) (V_Quaternion v)                 = G.basicUnsafeCopy mv v   elemseq _ (Quaternion a b c d) z = G.elemseq (undefined :: Vector a) a

@@ -1,55 +1,46 @@ {-# LANGUAGE DataKinds                  #-} {-# LANGUAGE FlexibleContexts           #-} {-# LANGUAGE FlexibleInstances          #-}-{-# LANGUAGE GeneralizedNewtypeDeriving #-} {-# LANGUAGE ScopedTypeVariables        #-} {-# LANGUAGE UndecidableInstances       #-} -{-# OPTIONS_GHC -fno-warn-orphans #-}- module Sparse   ( testSuite   , ShortPoly(..)   ) where -import Prelude hiding (gcd, quotRem, rem)+import Prelude hiding (gcd, quotRem, quot, rem)+import Control.Exception import Data.Euclidean (Euclidean(..), GcdDomain(..)) import Data.Function import Data.Int import Data.List (groupBy, sortOn)-import Data.Mod+import Data.Mod.Word import Data.Poly.Sparse import qualified Data.Poly.Sparse.Semiring as S import Data.Proxy-import Data.Semiring (Semiring)+import Data.Semiring (Semiring(..)) import qualified Data.Vector as V import qualified Data.Vector.Generic as G import qualified Data.Vector.Unboxed as U+import qualified Data.Vector.Unboxed.Sized as SU import Test.Tasty import Test.Tasty.QuickCheck hiding (scale, numTests)  import Quaternion import TestUtils -instance (Eq a, Semiring a, Arbitrary a, G.Vector v (Word, a)) => Arbitrary (Poly v a) where-  arbitrary = S.toPoly . G.fromList <$> arbitrary-  shrink = fmap (S.toPoly . G.fromList) . shrink . G.toList . unPoly--newtype ShortPoly a = ShortPoly { unShortPoly :: a }-  deriving (Eq, Show, Semiring, GcdDomain, Euclidean)--instance (Eq a, Semiring a, Arbitrary a, G.Vector v (Word, a)) => Arbitrary (ShortPoly (Poly v a)) where-  arbitrary = ShortPoly . S.toPoly . G.fromList . (\xs -> take (length xs `mod` 5) xs) <$> arbitrary-  shrink = fmap (ShortPoly . S.toPoly . G.fromList) . shrink . G.toList . unPoly . unShortPoly- testSuite :: TestTree testSuite = testGroup "Sparse"-    [ arithmeticTests-    , otherTests-    , lawsTests-    , evalTests-    , derivTests-    ]+  [ arithmeticTests+  , otherTests+  , divideByZeroTests+  , lawsTests+  , evalTests+  , derivTests+  , patternTests+  , conversionTests+  ]  lawsTests :: TestTree lawsTests = testGroup "Laws"@@ -57,60 +48,60 @@  semiringTests :: [TestTree] semiringTests =-  [ mySemiringLaws (Proxy :: Proxy (Poly U.Vector ()))-  , mySemiringLaws (Proxy :: Proxy (Poly U.Vector Int8))-  , mySemiringLaws (Proxy :: Proxy (Poly V.Vector Integer))+  [ mySemiringLaws (Proxy :: Proxy (UPoly ()))+  , mySemiringLaws (Proxy :: Proxy (UPoly Int8))+  , mySemiringLaws (Proxy :: Proxy (VPoly Integer))   , tenTimesLess-  $ mySemiringLaws (Proxy :: Proxy (Poly U.Vector (Quaternion Int)))+  $ mySemiringLaws (Proxy :: Proxy (UPoly (Quaternion Int)))   ]  ringTests :: [TestTree] ringTests =-  [ myRingLaws (Proxy :: Proxy (Poly U.Vector ()))-  , myRingLaws (Proxy :: Proxy (Poly U.Vector Int8))-  , myRingLaws (Proxy :: Proxy (Poly V.Vector Integer))-  , myRingLaws (Proxy :: Proxy (Poly U.Vector (Quaternion Int)))+  [ myRingLaws (Proxy :: Proxy (UPoly ()))+  , myRingLaws (Proxy :: Proxy (UPoly Int8))+  , myRingLaws (Proxy :: Proxy (VPoly Integer))+  , myRingLaws (Proxy :: Proxy (UPoly (Quaternion Int)))   ]  numTests :: [TestTree] numTests =-  [ myNumLaws (Proxy :: Proxy (Poly U.Vector Int8))-  , myNumLaws (Proxy :: Proxy (Poly V.Vector Integer))+  [ myNumLaws (Proxy :: Proxy (UPoly Int8))+  , myNumLaws (Proxy :: Proxy (VPoly Integer))   , tenTimesLess-  $ myNumLaws (Proxy :: Proxy (Poly U.Vector (Quaternion Int)))+  $ myNumLaws (Proxy :: Proxy (UPoly (Quaternion Int)))   ]  gcdDomainTests :: [TestTree] gcdDomainTests =-  [ myGcdDomainLaws (Proxy :: Proxy (ShortPoly (Poly V.Vector Integer)))+  [ myGcdDomainLaws (Proxy :: Proxy (ShortPoly (VPoly Integer)))   , tenTimesLess-  $ myGcdDomainLaws (Proxy :: Proxy (ShortPoly (Poly V.Vector (Mod 3))))+  $ myGcdDomainLaws (Proxy :: Proxy (ShortPoly (UPoly (Mod 3))))   , tenTimesLess-  $ myGcdDomainLaws (Proxy :: Proxy (ShortPoly (Poly V.Vector Rational)))+  $ myGcdDomainLaws (Proxy :: Proxy (ShortPoly (VPoly Rational)))   ]  euclideanTests :: [TestTree] euclideanTests =-  [ myEuclideanLaws (Proxy :: Proxy (ShortPoly (Poly V.Vector (Mod 3))))-  , myEuclideanLaws (Proxy :: Proxy (ShortPoly (Poly V.Vector Rational)))+  [ myEuclideanLaws (Proxy :: Proxy (ShortPoly (UPoly (Mod 3))))+  , myEuclideanLaws (Proxy :: Proxy (ShortPoly (VPoly Rational)))   ]  isListTests :: [TestTree] isListTests =-  [ myIsListLaws (Proxy :: Proxy (Poly U.Vector ()))-  , myIsListLaws (Proxy :: Proxy (Poly U.Vector Int8))-  , myIsListLaws (Proxy :: Proxy (Poly V.Vector Integer))+  [ myIsListLaws (Proxy :: Proxy (UPoly ()))+  , myIsListLaws (Proxy :: Proxy (UPoly Int8))+  , myIsListLaws (Proxy :: Proxy (VPoly Integer))   , tenTimesLess-  $ myIsListLaws (Proxy :: Proxy (Poly U.Vector (Quaternion Int)))+  $ myIsListLaws (Proxy :: Proxy (UPoly (Quaternion Int)))   ]  showTests :: [TestTree] showTests =-  [ myShowLaws (Proxy :: Proxy (Poly U.Vector ()))-  , myShowLaws (Proxy :: Proxy (Poly U.Vector Int8))-  , myShowLaws (Proxy :: Proxy (Poly V.Vector Integer))+  [ myShowLaws (Proxy :: Proxy (UPoly ()))+  , myShowLaws (Proxy :: Proxy (UPoly Int8))+  , myShowLaws (Proxy :: Proxy (VPoly Integer))   , tenTimesLess-  $ myShowLaws (Proxy :: Proxy (Poly U.Vector (Quaternion Int)))+  $ myShowLaws (Proxy :: Proxy (UPoly (Quaternion Int)))   ]  arithmeticTests :: TestTree@@ -125,6 +116,9 @@     testProperty "multiplication matches reference" $     \(xs :: [(Word, Int)]) ys -> toPoly (V.fromList (mulRef xs ys)) ===       toPoly (V.fromList xs) * toPoly (V.fromList ys)+  , tenTimesLess $+    testProperty "quotRemFractional matches quotRem" $+    \(xs :: VPoly Rational) ys -> ys /= 0 ==> quotRemFractional xs ys === quotRem xs ys   ]  addRef :: Num a => [(Word, a)] -> [(Word, a)] -> [(Word, a)]@@ -173,37 +167,51 @@   , tenTimesLess $     testProperty "scale matches multiplication by monomial" $     \p c (xs :: UPoly a) -> scale p c xs === monomial p c * xs+  , tenTimesLess $+    testProperty "scale' matches multiplication by monomial'" $+    \p c (xs :: UPoly a) -> S.scale p c xs === S.monomial p c * xs   ]  monomialRef :: Num a => Word -> a -> [(Word, a)] monomialRef p c = [(p, c)] +divideByZeroTests :: TestTree+divideByZeroTests = testGroup "divideByZero"+  [ testProperty "quotRem" $ testProp ((uncurry (+) .) . quotRem)+  , testProperty "quot"    $ testProp quot+  , testProperty "rem"     $ testProp rem+  , testProperty "divide"  $ testProp divide+  , testProperty "degree"  $ once $ degree (0 :: VPoly Rational) === 0+  ]+  where+    testProp f xs = ioProperty ((== Left DivideByZero) <$> try (evaluate (xs `f` (0 :: VPoly Rational))))+ evalTests :: TestTree evalTests = testGroup "eval" $ concat-  [ evalTestGroup  (Proxy :: Proxy (Poly U.Vector Int8))-  , evalTestGroup  (Proxy :: Proxy (Poly V.Vector Integer))-  , substTestGroup (Proxy :: Proxy (Poly U.Vector Int8))+  [ evalTestGroup  (Proxy :: Proxy (UPoly Int8))+  , evalTestGroup  (Proxy :: Proxy (VPoly Integer))+  , substTestGroup (Proxy :: Proxy (UPoly Int8))   ]  evalTestGroup   :: forall v a.-     (Eq a, Num a, Semiring a, Arbitrary a, Show a, Eq (v (Word, a)), Show (v (Word, a)), G.Vector v (Word, a))+     (Eq a, Num a, Semiring a, Arbitrary a, Show a, Eq (v (Word, a)), Show (v (Word, a)), G.Vector v (Word, a), G.Vector v (SU.Vector 1 Word, a))   => Proxy (Poly v a)   -> [TestTree] evalTestGroup _ =   [ testProperty "eval (p + q) r = eval p r + eval q r" $-    \p q r -> e (p + q) r === e p r + e q r+    \(ShortPoly p) (ShortPoly q) r -> e (p + q) r === e p r + e q r   , testProperty "eval (p * q) r = eval p r * eval q r" $-    \p q r -> e (p * q) r === e p r * e q r+    \(ShortPoly p) (ShortPoly q) r -> e (p * q) r === e p r * e q r   , testProperty "eval x p = p" $     \p -> e X p === p   , testProperty "eval (monomial 0 c) p = c" $     \c p -> e (monomial 0 c) p === c    , testProperty "eval' (p + q) r = eval' p r + eval' q r" $-    \p q r -> e' (p + q) r === e' p r + e' q r+    \(ShortPoly p) (ShortPoly q) r -> e' (p + q) r === e' p r + e' q r   , testProperty "eval' (p * q) r = eval' p r * eval' q r" $-    \p q r -> e' (p * q) r === e' p r * e' q r+    \(ShortPoly p) (ShortPoly q) r -> e' (p * q) r === e' p r * e' q r   , testProperty "eval' x p = p" $     \p -> e' S.X p === p   , testProperty "eval' (S.monomial 0 c) p = c" $@@ -218,7 +226,7 @@  substTestGroup   :: forall v a.-     (Eq a, Num a, Semiring a, Arbitrary a, Show a, Eq (v (Word, a)), Show (v (Word, a)), G.Vector v (Word, a))+     (Eq a, Num a, Semiring a, Arbitrary a, Show a, Eq (v (SU.Vector 1 Word, a)), Show (v (SU.Vector 1 Word, a)), G.Vector v (Word, a), G.Vector v (SU.Vector 1 Word, a))   => Proxy (Poly v a)   -> [TestTree] substTestGroup _ =@@ -240,20 +248,45 @@ derivTests :: TestTree derivTests = testGroup "deriv"   [ testProperty "deriv = S.deriv" $-    \(p :: Poly V.Vector Integer) -> deriv p === S.deriv p+    \(p :: VPoly Integer) -> deriv p === S.deriv p   , testProperty "integral = S.integral" $-    \(p :: Poly V.Vector Rational) -> integral p === S.integral p+    \(p :: VPoly Rational) -> integral p === S.integral p   , testProperty "deriv . integral = id" $-    \(p :: Poly V.Vector Rational) -> deriv (integral p) === p+    \(p :: VPoly Rational) -> deriv (integral p) === p   , testProperty "deriv c = 0" $-    \c -> deriv (monomial 0 c :: Poly V.Vector Int) === 0+    \c -> deriv (monomial 0 c :: UPoly Int) === 0   , testProperty "deriv cX = c" $-    \c -> deriv (monomial 0 c * X :: Poly V.Vector Int) === monomial 0 c+    \c -> deriv (monomial 0 c * X :: UPoly Int) === monomial 0 c   , testProperty "deriv (p + q) = deriv p + deriv q" $-    \p q -> deriv (p + q) === (deriv p + deriv q :: Poly V.Vector Int)+    \p q -> deriv (p + q) === (deriv p + deriv q :: UPoly Int)   , testProperty "deriv (p * q) = p * deriv q + q * deriv p" $-    \p q -> deriv (p * q) === (p * deriv q + q * deriv p :: Poly V.Vector Int)-  -- , testProperty "deriv (subst p q) = deriv q * subst (deriv p) q" $-  --   \(p :: Poly V.Vector Int) (q :: Poly U.Vector Int) ->-  --     deriv (subst p q) === deriv q * subst (deriv p) q+    \p q -> deriv (p * q) === (p * deriv q + q * deriv p :: UPoly Int)+  ]++patternTests :: TestTree+patternTests = testGroup "pattern"+  [ testProperty "X  :: UPoly Int" $ once $+    case (monomial 1 1 :: UPoly Int) of X -> True; _ -> False+  , testProperty "X  :: UPoly Int" $ once $+    (X :: UPoly Int) === monomial 1 1+  , testProperty "X' :: UPoly Int" $ once $+    case (S.monomial 1 1 :: UPoly Int) of S.X -> True; _ -> False+  , testProperty "X' :: UPoly Int" $ once $+    (S.X :: UPoly Int) === S.monomial 1 1+  , testProperty "X' :: UPoly ()" $ once $+    case (zero :: UPoly ()) of S.X -> True; _ -> False+  , testProperty "X' :: UPoly ()" $ once $+    (S.X :: UPoly ()) === zero+  ]++conversionTests :: TestTree+conversionTests = testGroup "conversions"+  [ testProperty "denseToSparse . sparseToDense = id" $+    \(xs :: UPoly Int8) -> xs === denseToSparse (sparseToDense xs)+  , testProperty "denseToSparse' . sparseToDense' = id" $+    \(xs :: UPoly Int8) -> xs === S.denseToSparse (S.sparseToDense xs)+  , 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)   ]

@@ -1,50 +1,39 @@+{-# LANGUAGE DataKinds                  #-} {-# LANGUAGE FlexibleContexts           #-} {-# LANGUAGE FlexibleInstances          #-}-{-# LANGUAGE GeneralizedNewtypeDeriving #-} {-# LANGUAGE ScopedTypeVariables        #-} {-# LANGUAGE UndecidableInstances       #-} -{-# OPTIONS_GHC -fno-warn-orphans #-}- module SparseLaurent   ( testSuite   ) where -import Prelude hiding (gcd, quotRem, rem)-import Data.Euclidean (Euclidean(..), GcdDomain(..), Field)+import Prelude hiding (gcd, quotRem, quot, rem)+import Control.Exception+import Data.Euclidean (GcdDomain(..), Field) import Data.Int import qualified Data.Poly.Sparse import Data.Poly.Sparse.Laurent import Data.Proxy import Data.Semiring (Semiring(..))-import qualified Data.Vector as V import qualified Data.Vector.Generic as G import qualified Data.Vector.Unboxed as U+import qualified Data.Vector.Unboxed.Sized as SU import Test.Tasty import Test.Tasty.QuickCheck hiding (scale, numTests)  import Quaternion-import Sparse (ShortPoly(..)) import TestUtils -instance (Eq a, Semiring a, Arbitrary a, G.Vector v (Word, a)) => Arbitrary (Laurent v a) where-  arbitrary = toLaurent <$> ((`rem` 10) <$> arbitrary) <*> arbitrary-  shrink = fmap (uncurry toLaurent) . shrink . unLaurent--newtype ShortLaurent a = ShortLaurent { unShortLaurent :: a }-  deriving (Eq, Show, Semiring, GcdDomain)--instance (Eq a, Semiring a, Arbitrary a, G.Vector v (Word, a)) => Arbitrary (ShortLaurent (Laurent v a)) where-  arbitrary = (ShortLaurent .) . toLaurent <$> ((`rem` 10) <$> arbitrary) <*> (unShortPoly <$> arbitrary)-  shrink = fmap (ShortLaurent . uncurry toLaurent . fmap unShortPoly) . shrink . fmap ShortPoly . unLaurent . unShortLaurent- testSuite :: TestTree testSuite = testGroup "SparseLaurent"-    [ otherTests-    , lawsTests-    , evalTests-    , derivTests-    ]+  [ otherTests+  , divideByZeroTests+  , lawsTests+  , evalTests+  , derivTests+  , patternTests+  ]  lawsTests :: TestTree lawsTests = testGroup "Laws"@@ -52,52 +41,52 @@  semiringTests :: [TestTree] semiringTests =-  [ mySemiringLaws (Proxy :: Proxy (Laurent U.Vector ()))-  , mySemiringLaws (Proxy :: Proxy (Laurent U.Vector Int8))-  , mySemiringLaws (Proxy :: Proxy (Laurent V.Vector Integer))+  [ mySemiringLaws (Proxy :: Proxy (ULaurent ()))+  , mySemiringLaws (Proxy :: Proxy (ULaurent Int8))+  , mySemiringLaws (Proxy :: Proxy (VLaurent Integer))   , tenTimesLess-  $ mySemiringLaws (Proxy :: Proxy (Laurent U.Vector (Quaternion Int)))+  $ mySemiringLaws (Proxy :: Proxy (ULaurent (Quaternion Int)))   ]  ringTests :: [TestTree] ringTests =-  [ myRingLaws (Proxy :: Proxy (Laurent U.Vector ()))-  , myRingLaws (Proxy :: Proxy (Laurent U.Vector Int8))-  , myRingLaws (Proxy :: Proxy (Laurent V.Vector Integer))-  , myRingLaws (Proxy :: Proxy (Laurent U.Vector (Quaternion Int)))+  [ myRingLaws (Proxy :: Proxy (ULaurent ()))+  , myRingLaws (Proxy :: Proxy (ULaurent Int8))+  , myRingLaws (Proxy :: Proxy (VLaurent Integer))+  , myRingLaws (Proxy :: Proxy (ULaurent (Quaternion Int)))   ]  numTests :: [TestTree] numTests =-  [ myNumLaws (Proxy :: Proxy (Laurent U.Vector Int8))-  , myNumLaws (Proxy :: Proxy (Laurent V.Vector Integer))+  [ myNumLaws (Proxy :: Proxy (ULaurent Int8))+  , myNumLaws (Proxy :: Proxy (VLaurent Integer))   , tenTimesLess-  $ myNumLaws (Proxy :: Proxy (Laurent U.Vector (Quaternion Int)))+  $ myNumLaws (Proxy :: Proxy (ULaurent (Quaternion Int)))   ]  gcdDomainTests :: [TestTree] gcdDomainTests =-  [ myGcdDomainLaws (Proxy :: Proxy (ShortLaurent (Laurent V.Vector Integer)))+  [ myGcdDomainLaws (Proxy :: Proxy (ShortPoly (VLaurent Integer)))   , tenTimesLess-  $ myGcdDomainLaws (Proxy :: Proxy (ShortLaurent (Laurent V.Vector Rational)))+  $ myGcdDomainLaws (Proxy :: Proxy (ShortPoly (VLaurent Rational)))   ]  isListTests :: [TestTree] isListTests =-  [ myIsListLaws (Proxy :: Proxy (Laurent U.Vector ()))-  , myIsListLaws (Proxy :: Proxy (Laurent U.Vector Int8))-  , myIsListLaws (Proxy :: Proxy (Laurent V.Vector Integer))+  [ myIsListLaws (Proxy :: Proxy (ULaurent ()))+  , myIsListLaws (Proxy :: Proxy (ULaurent Int8))+  , myIsListLaws (Proxy :: Proxy (VLaurent Integer))   , tenTimesLess-  $ myIsListLaws (Proxy :: Proxy (Laurent U.Vector (Quaternion Int)))+  $ myIsListLaws (Proxy :: Proxy (ULaurent (Quaternion Int)))   ]  showTests :: [TestTree] showTests =-  [ myShowLaws (Proxy :: Proxy (Laurent U.Vector ()))-  , myShowLaws (Proxy :: Proxy (Laurent U.Vector Int8))-  , myShowLaws (Proxy :: Proxy (Laurent V.Vector Integer))+  [ myShowLaws (Proxy :: Proxy (ULaurent ()))+  , myShowLaws (Proxy :: Proxy (ULaurent Int8))+  , myShowLaws (Proxy :: Proxy (VLaurent Integer))   , tenTimesLess-  $ myShowLaws (Proxy :: Proxy (Laurent U.Vector (Quaternion Int)))+  $ myShowLaws (Proxy :: Proxy (ULaurent (Quaternion Int)))   ]  otherTests :: TestTree@@ -119,24 +108,34 @@   , tenTimesLess $     testProperty "scale matches multiplication by monomial" $     \p c (xs :: ULaurent a) -> scale p c xs === monomial p c * xs+  , tenTimesLess $+    testProperty "toLaurent . unLaurent" $+    \(xs :: ULaurent a) -> uncurry toLaurent (unLaurent xs) === xs   ] +divideByZeroTests :: TestTree+divideByZeroTests = testGroup "divideByZero"+  [ testProperty "divide"  $ testProp divide+  ]+  where+    testProp f xs = ioProperty ((== Left DivideByZero) <$> try (evaluate (xs `f` (0 :: VLaurent Rational))))+ evalTests :: TestTree evalTests = testGroup "eval" $ concat-  [ evalTestGroup  (Proxy :: Proxy (Laurent V.Vector Rational))-  , substTestGroup (Proxy :: Proxy (Laurent U.Vector Int8))+  [ evalTestGroup  (Proxy :: Proxy (VLaurent Rational))+  , substTestGroup (Proxy :: Proxy (ULaurent Int8))   ]  evalTestGroup   :: forall v a.-     (Eq a, Field a, Arbitrary a, Show a, Eq (v (Word, a)), Show (v (Word, a)), G.Vector v (Word, a))+     (Eq a, Field a, Arbitrary a, Show a, Eq (v (Word, a)), Show (v (Word, a)), G.Vector v (Word, a), G.Vector v (SU.Vector 1 Word, a))   => Proxy (Laurent v a)   -> [TestTree] evalTestGroup _ =   [ testProperty "eval (p + q) r = eval p r + eval q r" $-    \p q r -> e (p `plus` q) r === e p r `plus` e q r+    \(ShortPoly p) (ShortPoly q) r -> e (p `plus` q) r === e p r `plus` e q r   , testProperty "eval (p * q) r = eval p r * eval q r" $-    \p q r -> e (p `times` q) r === e p r `times` e q r+    \(ShortPoly p) (ShortPoly q) r -> e (p `times` q) r === e p r `times` e q r   , testProperty "eval x p = p" $     \p -> e X p === p   , testProperty "eval (monomial 0 c) p = c" $@@ -148,7 +147,7 @@  substTestGroup   :: forall v a.-     (Eq a, Num a, Semiring a, Arbitrary a, Show a, Eq (v (Word, a)), Show (v (Word, a)), G.Vector v (Word, a))+     (Eq a, Num a, Semiring a, Arbitrary a, Show a, Eq (v (SU.Vector 1 Word, a)), Show (v (Word, a)), G.Vector v (Word, a), G.Vector v (SU.Vector 1 Word, a))   => Proxy (Laurent v a)   -> [TestTree] substTestGroup _ =@@ -164,14 +163,31 @@ derivTests :: TestTree derivTests = testGroup "deriv"   [ testProperty "deriv c = 0" $-    \c -> deriv (monomial 0 c :: Laurent V.Vector Int) === 0+    \c -> deriv (monomial 0 c :: ULaurent Int) === 0   , testProperty "deriv cX = c" $-    \c -> deriv (monomial 0 c * X :: Laurent V.Vector Int) === monomial 0 c+    \c -> deriv (monomial 0 c * X :: ULaurent Int) === monomial 0 c   , testProperty "deriv (p + q) = deriv p + deriv q" $-    \p q -> deriv (p + q) === (deriv p + deriv q :: Laurent V.Vector Int)+    \p q -> deriv (p + q) === (deriv p + deriv q :: ULaurent Int)   , testProperty "deriv (p * q) = p * deriv q + q * deriv p" $-    \p q -> deriv (p * q) === (p * deriv q + q * deriv p :: Laurent V.Vector Int)-  -- , testProperty "deriv (subst p q) = deriv q * subst (deriv p) q" $-  --   \(p :: Laurent V.Vector Int) (q :: Laurent U.Vector Int) ->-  --     deriv (subst p q) === deriv q * subst (deriv p) q+    \p q -> deriv (p * q) === (p * deriv q + q * deriv p :: ULaurent Int)+  ]++patternTests :: TestTree+patternTests = testGroup "pattern"+  [ testProperty "X  :: ULaurent Int" $ once $+    case (monomial 1 1 :: ULaurent Int) of X -> True; _ -> False+  , testProperty "X  :: ULaurent Int" $ once $+    (X :: ULaurent Int) === monomial 1 1+  , testProperty "X :: ULaurent ()" $ once $+    case (zero :: ULaurent ()) of X -> True; _ -> False+  , testProperty "X :: ULaurent ()" $ once $+    (X :: ULaurent ()) === zero+  , testProperty "X^-k" $+    \(NonNegative j) k -> ((X^j)^-k :: ULaurent Int) === monomial (- j * k) 1+  , testProperty "^-" $+    \(p :: ULaurent Int) (NonNegative k) -> ioProperty $ do+      et <- try (evaluate (p^-k)) :: IO (Either PatternMatchFail (ULaurent Int))+      pure $ case et of+        Left{}  -> True+        Right t -> p^k * t == one   ]

@@ -1,10 +1,15 @@-{-# LANGUAGE CPP              #-}-{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE DataKinds                  #-}+{-# LANGUAGE FlexibleContexts           #-}+{-# LANGUAGE FlexibleInstances          #-}+{-# LANGUAGE GeneralizedNewtypeDeriving #-}+{-# LANGUAGE ScopedTypeVariables        #-}+{-# LANGUAGE UndecidableInstances       #-}  {-# OPTIONS_GHC -fno-warn-orphans #-}  module TestUtils-  ( tenTimesLess+  ( ShortPoly(..)+  , tenTimesLess   , mySemiringLaws   , myRingLaws   , myNumLaws@@ -14,25 +19,75 @@   , myShowLaws   ) where +import Prelude hiding (lcm, rem)+import Control.Arrow import Data.Euclidean-import Data.Mod+import Data.Finite+import Data.Mod.Word import Data.Proxy-import Data.Semiring (Semiring, Ring)+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.Tasty import Test.Tasty.QuickCheck -#if MIN_VERSION_base(4,10,0)-import GHC.TypeNats (KnownNat)-#else-import GHC.TypeLits (KnownNat)-#endif+import qualified Data.Poly.Semiring as Dense+import qualified Data.Poly.Laurent as DenseLaurent+import Data.Poly.Multi.Semiring+import qualified Data.Poly.Multi.Laurent as MultiLaurent +newtype ShortPoly a = ShortPoly { unShortPoly :: a }+  deriving (Eq, Show, Semiring, GcdDomain, Euclidean, Num)+ instance KnownNat m => Arbitrary (Mod m) where   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++instance (Eq a, Semiring a, Arbitrary a, G.Vector v a) => Arbitrary (ShortPoly (Dense.Poly v a)) where+  arbitrary = ShortPoly . Dense.toPoly . G.fromList . (\xs -> take (length xs `mod` 10) xs) <$> arbitrary+  shrink = fmap (ShortPoly . Dense.toPoly . G.fromList) . shrink . G.toList . Dense.unPoly . unShortPoly++instance (Eq a, Semiring a, Arbitrary a, G.Vector v a) => Arbitrary (DenseLaurent.Laurent v a) where+  arbitrary = DenseLaurent.toLaurent <$> ((`rem` 10) <$> arbitrary) <*> arbitrary+  shrink = fmap (uncurry DenseLaurent.toLaurent) . shrink . DenseLaurent.unLaurent++instance (Eq a, Semiring a, Arbitrary a, G.Vector v a) => Arbitrary (ShortPoly (DenseLaurent.Laurent v a)) where+  arbitrary = (ShortPoly .) . DenseLaurent.toLaurent <$> ((`rem` 10) <$> arbitrary) <*> (unShortPoly <$> arbitrary)+  shrink = fmap (ShortPoly . uncurry DenseLaurent.toLaurent . fmap unShortPoly) . shrink . fmap ShortPoly . DenseLaurent.unLaurent . unShortPoly++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++instance (Eq a, Semiring a, Arbitrary a, KnownNat n, G.Vector v (SU.Vector n Word, a)) => Arbitrary (ShortPoly (MultiPoly v n a)) where+  arbitrary = ShortPoly . toMultiPoly . G.fromList . (\xs -> take (length xs `mod` 4) (map (first (SU.map (`mod` 3))) xs)) <$> arbitrary+  shrink = fmap (ShortPoly . toMultiPoly . G.fromList) . shrink . G.toList . unMultiPoly . unShortPoly++instance (Eq a, Semiring a, Arbitrary a, KnownNat n, G.Vector v (Word, a), G.Vector v (SU.Vector n Word, a)) => Arbitrary (MultiLaurent.MultiLaurent v n a) where+  arbitrary = MultiLaurent.toMultiLaurent <$> (SU.map (`rem` 10) <$> arbitrary) <*> arbitrary+  shrink = fmap (uncurry MultiLaurent.toMultiLaurent) . shrink . MultiLaurent.unMultiLaurent++instance (Eq a, Semiring a, Arbitrary a, KnownNat n, G.Vector v (Word, a), G.Vector v (SU.Vector n Word, a)) => Arbitrary (ShortPoly (MultiLaurent.MultiLaurent v n a)) where+  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++-------------------------------------------------------------------------------+ tenTimesLess :: TestTree -> TestTree tenTimesLess = adjustOption $   \(QuickCheckTests n) -> QuickCheckTests (max 100 (n `div` 10))@@ -76,8 +131,8 @@       where         test = uncurry testProperty pair -myGcdDomainLaws :: (Eq a, GcdDomain a, Arbitrary a, Show a) => Proxy a -> TestTree-myGcdDomainLaws proxy = testGroup tpclss $ map tune props+myGcdDomainLaws :: forall a. (Eq a, GcdDomain a, Arbitrary a, Show a) => Proxy a -> TestTree+myGcdDomainLaws proxy = testGroup tpclss $ map tune $ lcm0 : props   where     Laws tpclss props = gcdDomainLaws proxy @@ -90,6 +145,8 @@       _ -> test       where         test = uncurry testProperty pair++    lcm0 = ("lcm0", property $ \(x :: a) -> lcm x zero === zero .&&. lcm zero x === zero)  myEuclideanLaws :: (Eq a, Euclidean a, Arbitrary a, Show a) => Proxy a -> TestTree myEuclideanLaws proxy = testGroup tpclss $ map (uncurry testProperty) props

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