diff --git a/README.md b/README.md --- a/README.md +++ b/README.md @@ -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. diff --git a/bench/DenseBench.hs b/bench/DenseBench.hs --- a/bench/DenseBench.hs +++ b/bench/DenseBench.hs @@ -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 diff --git a/bench/SparseBench.hs b/bench/SparseBench.hs --- a/bench/SparseBench.hs +++ b/bench/SparseBench.hs @@ -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 diff --git a/changelog.md b/changelog.md --- a/changelog.md +++ b/changelog.md @@ -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. diff --git a/poly.cabal b/poly.cabal --- a/poly.cabal +++ b/poly.cabal @@ -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 diff --git a/src/Data/Poly.hs b/src/Data/Poly.hs --- a/src/Data/Poly.hs +++ b/src/Data/Poly.hs @@ -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 diff --git a/src/Data/Poly/Internal/Convert.hs b/src/Data/Poly/Internal/Convert.hs new file mode 100644 --- /dev/null +++ b/src/Data/Poly/Internal/Convert.hs @@ -0,0 +1,88 @@ +-- | +-- Module: Data.Poly.Internal.Convert +-- Copyright: (c) 2020 Andrew Lelechenko +-- Licence: BSD3 +-- Maintainer: Andrew Lelechenko +-- +-- 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 diff --git a/src/Data/Poly/Internal/Dense.hs b/src/Data/Poly/Internal/Dense.hs --- a/src/Data/Poly/Internal/Dense.hs +++ b/src/Data/Poly/Internal/Dense.hs @@ -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' #-} diff --git a/src/Data/Poly/Internal/Dense/DFT.hs b/src/Data/Poly/Internal/Dense/DFT.hs new file mode 100644 --- /dev/null +++ b/src/Data/Poly/Internal/Dense/DFT.hs @@ -0,0 +1,80 @@ +-- | +-- Module: Data.Poly.Internal.Dense.FFT +-- Copyright: (c) 2020 Andrew Lelechenko +-- Licence: BSD3 +-- Maintainer: Andrew Lelechenko +-- +-- 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 + +-- | +-- \( 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 +-- \( 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 diff --git a/src/Data/Poly/Internal/Dense/Field.hs b/src/Data/Poly/Internal/Dense/Field.hs --- a/src/Data/Poly/Internal/Dense/Field.hs +++ b/src/Data/Poly/Internal/Dense/Field.hs @@ -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 #-} diff --git a/src/Data/Poly/Internal/Dense/GcdDomain.hs b/src/Data/Poly/Internal/Dense/GcdDomain.hs --- a/src/Data/Poly/Internal/Dense/GcdDomain.hs +++ b/src/Data/Poly/Internal/Dense/GcdDomain.hs @@ -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) diff --git a/src/Data/Poly/Internal/Dense/Laurent.hs b/src/Data/Poly/Internal/Dense/Laurent.hs new file mode 100644 --- /dev/null +++ b/src/Data/Poly/Internal/Dense/Laurent.hs @@ -0,0 +1,288 @@ +-- | +-- Module: Data.Poly.Internal.Dense.Laurent +-- Copyright: (c) 2020 Andrew Lelechenko +-- Licence: BSD3 +-- Maintainer: Andrew Lelechenko +-- +-- . +-- + +{-# 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 () + +-- | +-- 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 #-} diff --git a/src/Data/Poly/Internal/Multi.hs b/src/Data/Poly/Internal/Multi.hs new file mode 100644 --- /dev/null +++ b/src/Data/Poly/Internal/Multi.hs @@ -0,0 +1,553 @@ +-- | +-- Module: Data.Poly.Internal.Multi +-- Copyright: (c) 2020 Andrew Lelechenko +-- Licence: BSD3 +-- Maintainer: Andrew Lelechenko +-- + +{-# 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 diff --git a/src/Data/Poly/Internal/Multi/Core.hs b/src/Data/Poly/Internal/Multi/Core.hs new file mode 100644 --- /dev/null +++ b/src/Data/Poly/Internal/Multi/Core.hs @@ -0,0 +1,311 @@ +-- | +-- Module: Data.Poly.Internal.Multi.Core +-- Copyright: (c) 2019 Andrew Lelechenko +-- Licence: BSD3 +-- Maintainer: Andrew Lelechenko +-- +-- 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 #-} diff --git a/src/Data/Poly/Internal/Multi/Field.hs b/src/Data/Poly/Internal/Multi/Field.hs new file mode 100644 --- /dev/null +++ b/src/Data/Poly/Internal/Multi/Field.hs @@ -0,0 +1,73 @@ +-- | +-- Module: Data.Poly.Internal.Multi.Field +-- Copyright: (c) 2019 Andrew Lelechenko +-- Licence: BSD3 +-- Maintainer: Andrew Lelechenko +-- +-- 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) diff --git a/src/Data/Poly/Internal/Multi/GcdDomain.hs b/src/Data/Poly/Internal/Multi/GcdDomain.hs new file mode 100644 --- /dev/null +++ b/src/Data/Poly/Internal/Multi/GcdDomain.hs @@ -0,0 +1,179 @@ +-- | +-- Module: Data.Poly.Internal.Multi.GcdDomain +-- Copyright: (c) 2019 Andrew Lelechenko +-- Licence: BSD3 +-- Maintainer: Andrew Lelechenko +-- +-- 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 diff --git a/src/Data/Poly/Internal/Multi/Laurent.hs b/src/Data/Poly/Internal/Multi/Laurent.hs new file mode 100644 --- /dev/null +++ b/src/Data/Poly/Internal/Multi/Laurent.hs @@ -0,0 +1,507 @@ +-- | +-- Module: Data.Poly.Internal.Multi.Laurent +-- Copyright: (c) 2020 Andrew Lelechenko +-- Licence: BSD3 +-- Maintainer: Andrew Lelechenko +-- +-- Sparse multivariate +-- . +-- + +{-# 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 +-- +-- 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 + +-- | +-- 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 diff --git a/src/Data/Poly/Internal/PolyOverField.hs b/src/Data/Poly/Internal/PolyOverField.hs deleted file mode 100644 --- a/src/Data/Poly/Internal/PolyOverField.hs +++ /dev/null @@ -1,46 +0,0 @@ --- | --- Module: Data.Poly.Internal.PolyOverField --- Copyright: (c) 2019 Andrew Lelechenko --- Licence: BSD3 --- Maintainer: Andrew Lelechenko --- --- 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 #-} diff --git a/src/Data/Poly/Internal/Sparse.hs b/src/Data/Poly/Internal/Sparse.hs deleted file mode 100644 --- a/src/Data/Poly/Internal/Sparse.hs +++ /dev/null @@ -1,583 +0,0 @@ --- | --- Module: Data.Poly.Internal.Sparse --- Copyright: (c) 2019 Andrew Lelechenko --- Licence: BSD3 --- Maintainer: Andrew Lelechenko --- --- 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' #-} diff --git a/src/Data/Poly/Internal/Sparse/Field.hs b/src/Data/Poly/Internal/Sparse/Field.hs deleted file mode 100644 --- a/src/Data/Poly/Internal/Sparse/Field.hs +++ /dev/null @@ -1,56 +0,0 @@ --- | --- Module: Data.Poly.Internal.Sparse.Field --- Copyright: (c) 2019 Andrew Lelechenko --- Licence: BSD3 --- Maintainer: Andrew Lelechenko --- --- 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 diff --git a/src/Data/Poly/Internal/Sparse/GcdDomain.hs b/src/Data/Poly/Internal/Sparse/GcdDomain.hs deleted file mode 100644 --- a/src/Data/Poly/Internal/Sparse/GcdDomain.hs +++ /dev/null @@ -1,74 +0,0 @@ --- | --- Module: Data.Poly.Internal.Sparse.GcdDomain --- Copyright: (c) 2019 Andrew Lelechenko --- Licence: BSD3 --- Maintainer: Andrew Lelechenko --- --- 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" diff --git a/src/Data/Poly/Laurent.hs b/src/Data/Poly/Laurent.hs --- a/src/Data/Poly/Laurent.hs +++ b/src/Data/Poly/Laurent.hs @@ -7,11 +7,7 @@ -- . -- -{-# 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 - --- | --- 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 diff --git a/src/Data/Poly/Multi.hs b/src/Data/Poly/Multi.hs new file mode 100644 --- /dev/null +++ b/src/Data/Poly/Multi.hs @@ -0,0 +1,35 @@ +-- | +-- Module: Data.Poly.Multi +-- Copyright: (c) 2020 Andrew Lelechenko +-- Licence: BSD3 +-- Maintainer: Andrew Lelechenko +-- +-- 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 () diff --git a/src/Data/Poly/Multi/Laurent.hs b/src/Data/Poly/Multi/Laurent.hs new file mode 100644 --- /dev/null +++ b/src/Data/Poly/Multi/Laurent.hs @@ -0,0 +1,32 @@ +-- | +-- Module: Data.Poly.Multi.Laurent +-- Copyright: (c) 2020 Andrew Lelechenko +-- Licence: BSD3 +-- Maintainer: Andrew Lelechenko +-- +-- Sparse multivariate +-- . +-- + +{-# 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 diff --git a/src/Data/Poly/Multi/Semiring.hs b/src/Data/Poly/Multi/Semiring.hs new file mode 100644 --- /dev/null +++ b/src/Data/Poly/Multi/Semiring.hs @@ -0,0 +1,157 @@ +-- | +-- Module: Data.Poly.Multi.Semiring +-- Copyright: (c) 2020 Andrew Lelechenko +-- Licence: BSD3 +-- Maintainer: Andrew Lelechenko +-- +-- 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' diff --git a/src/Data/Poly/Orthogonal.hs b/src/Data/Poly/Orthogonal.hs --- a/src/Data/Poly/Orthogonal.hs +++ b/src/Data/Poly/Orthogonal.hs @@ -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' . -- -- >>> 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 diff --git a/src/Data/Poly/Semiring.hs b/src/Data/Poly/Semiring.hs --- a/src/Data/Poly/Semiring.hs +++ b/src/Data/Poly/Semiring.hs @@ -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 +-- . +-- 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' diff --git a/src/Data/Poly/Sparse.hs b/src/Data/Poly/Sparse.hs --- a/src/Data/Poly/Sparse.hs +++ b/src/Data/Poly/Sparse.hs @@ -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 diff --git a/src/Data/Poly/Sparse/Laurent.hs b/src/Data/Poly/Sparse/Laurent.hs --- a/src/Data/Poly/Sparse/Laurent.hs +++ b/src/Data/Poly/Sparse/Laurent.hs @@ -4,18 +4,13 @@ -- Licence: BSD3 -- Maintainer: Andrew Lelechenko -- --- Sparse . +-- Sparse +-- . -- +{-# 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 () - --- | --- 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 diff --git a/src/Data/Poly/Sparse/Semiring.hs b/src/Data/Poly/Sparse/Semiring.hs --- a/src/Data/Poly/Sparse/Semiring.hs +++ b/src/Data/Poly/Sparse/Semiring.hs @@ -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' diff --git a/test/DFT.hs b/test/DFT.hs new file mode 100644 --- /dev/null +++ b/test/DFT.hs @@ -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' diff --git a/test/Dense.hs b/test/Dense.hs --- a/test/Dense.hs +++ b/test/Dense.hs @@ -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) ] diff --git a/test/DenseLaurent.hs b/test/DenseLaurent.hs --- a/test/DenseLaurent.hs +++ b/test/DenseLaurent.hs @@ -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 ] diff --git a/test/Main.hs b/test/Main.hs --- a/test/Main.hs +++ b/test/Main.hs @@ -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 ] diff --git a/test/Multi.hs b/test/Multi.hs new file mode 100644 --- /dev/null +++ b/test/Multi.hs @@ -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) + ] diff --git a/test/MultiLaurent.hs b/test/MultiLaurent.hs new file mode 100644 --- /dev/null +++ b/test/MultiLaurent.hs @@ -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) + ] diff --git a/test/Orthogonal.hs b/test/Orthogonal.hs --- a/test/Orthogonal.hs +++ b/test/Orthogonal.hs @@ -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 diff --git a/test/Quaternion.hs b/test/Quaternion.hs --- a/test/Quaternion.hs +++ b/test/Quaternion.hs @@ -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 diff --git a/test/Sparse.hs b/test/Sparse.hs --- a/test/Sparse.hs +++ b/test/Sparse.hs @@ -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) ] diff --git a/test/SparseLaurent.hs b/test/SparseLaurent.hs --- a/test/SparseLaurent.hs +++ b/test/SparseLaurent.hs @@ -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 ] diff --git a/test/TestUtils.hs b/test/TestUtils.hs --- a/test/TestUtils.hs +++ b/test/TestUtils.hs @@ -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 diff --git a/test/doctests.hs b/test/doctests.hs new file mode 100644 --- /dev/null +++ b/test/doctests.hs @@ -0,0 +1,4 @@ +import Test.DocTest (doctest) + +main :: IO () +main = doctest ["src"]