packages feed

poly 0.3.3.0 → 0.4.0.0

raw patch · 26 files changed

+1506/−568 lines, 26 filesdep +moddep ~quickcheck-classesdep ~semirings

Dependencies added: mod

Dependency ranges changed: quickcheck-classes, semirings

Files

README.md view
@@ -32,9 +32,9 @@ 1 * X^2 + (-3) * X + 2 ``` -(Unfortunately, a type is often ambiguous and must be given explicitly.)+(Unfortunately, types are often ambiguous and must be given explicitly.) -While being convenient to read and write in REPL, `X` is relatively slow. The fastest approach is to use `toPoly`, providing it with a vector of coefficients (head is the constant term):+While being convenient to read and write in REPL, `X` is relatively slow. The fastest approach is to use `toPoly`, providing it with a vector of coefficients (constant term first):  ```haskell > toPoly (Data.Vector.Unboxed.fromList [2, -3, 1 :: Int])@@ -99,7 +99,7 @@  ## Deconstruction -Use `unPoly` to deconstruct a polynomial to a vector of coefficients (head is the constant term):+Use `unPoly` to deconstruct a polynomial to a vector of coefficients (constant term first):  ```haskell > unPoly (X^2 - 3 * X + 2 :: UPoly Int)
bench/DenseBench.hs view
@@ -1,4 +1,3 @@-{-# LANGUAGE CPP                        #-} {-# LANGUAGE GeneralizedNewtypeDeriving #-} {-# LANGUAGE RankNTypes                 #-} {-# LANGUAGE TypeApplications           #-}@@ -9,15 +8,13 @@  import Prelude hiding (quotRem, gcd) import Gauge.Main-import Data.Poly-import qualified Data.Vector.Unboxed as U-#if MIN_VERSION_semirings(0,5,2) import Data.Euclidean (Euclidean(..), GcdDomain(..), Field)+import Data.Poly import qualified Data.Poly.Semiring as S (toPoly) import Data.Semiring (Semiring(..), Ring, Mod2(..)) import qualified Data.Semiring as S (fromIntegral) import qualified Data.Vector as V-#endif+import qualified Data.Vector.Unboxed as U  benchSuite :: Benchmark benchSuite = bgroup "dense" $ concat@@ -26,14 +23,10 @@   , map benchEval     [100, 1000, 10000]   , map benchDeriv    [100, 1000, 10000]   , map benchIntegral [100, 1000, 10000]-#if MIN_VERSION_semirings(0,5,2)   , map benchQuotRem    [10, 100]   , map benchGcd        [10, 100]-  , map benchGcdExtRat  [10, 20, 40]   , map benchGcdFracRat [10, 20, 40]-  , map benchGcdExtM    [10, 100, 1000]   , map benchGcdFracM   [10, 100, 1000]-#endif   ]  benchAdd :: Int -> Benchmark@@ -51,28 +44,18 @@ benchIntegral :: Int -> Benchmark benchIntegral k = bench ("integral/" ++ show k) $ nf doIntegral k -#if MIN_VERSION_semirings(0,5,2)- benchQuotRem :: Int -> Benchmark benchQuotRem k = bench ("quotRem/" ++ show k) $ nf doQuotRem k  benchGcd :: Int -> Benchmark benchGcd k = bench ("gcd/" ++ show k) $ nf doGcd k -benchGcdExtRat :: Int -> Benchmark-benchGcdExtRat k = bench ("gcdExt/Rational/" ++ show k) $ nf (doGcdExt @Rational) k- benchGcdFracRat :: Int -> Benchmark benchGcdFracRat k = bench ("gcdFrac/Rational/" ++ show k) $ nf (doGcdFrac @Rational) k -benchGcdExtM :: Int -> Benchmark-benchGcdExtM k = bench ("gcdExt/Mod2/" ++ show k) $ nf (getMod2 . doGcdExt @Mod2) k- benchGcdFracM :: Int -> Benchmark benchGcdFracM k = bench ("gcdFrac/Mod2/" ++ show k) $ nf (getMod2 . doGcdFrac @Mod2) k -#endif- doBinOp :: (forall a. Num a => a -> a -> a) -> Int -> Int doBinOp op n = U.sum zs   where@@ -98,8 +81,6 @@     xs = toPoly $ U.generate n ((* 2) . fromIntegral)     zs = unPoly $ integral xs -#if MIN_VERSION_semirings(0,5,2)- gen1 :: Ring a => Int -> a gen1 k = S.fromIntegral (truncate (pi * fromIntegral k :: Double) `mod` (k + 1)) @@ -120,18 +101,9 @@     ys = toPoly $ V.generate n gen2     gs = unPoly $ xs `gcd` ys -doGcdExt :: (Eq a, Field a) => Int -> a-doGcdExt n = V.foldl' plus zero gs-  where-    xs = S.toPoly $ V.generate n gen1-    ys = S.toPoly $ V.generate n gen2-    gs = unPoly $ fst $ xs `gcdExt` 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--#endif
changelog.md view
@@ -1,3 +1,10 @@+# 0.4.0.0++* Implement Laurent polynomials.+* Implement orthogonal polynomials.+* Decomission extended GCD, use `Data.Euclidean.gcdExt`.+* Decomission `PolyOverFractional`, use `PolyOverField`.+ # 0.3.3.0  * Add function `subst`.
poly.cabal view
@@ -1,5 +1,5 @@ name: poly-version: 0.3.3.0+version: 0.4.0.0 synopsis: Polynomials description:   Polynomials backed by `Vector`.@@ -8,14 +8,14 @@ license-file: LICENSE author: Andrew Lelechenko maintainer: andrew.lelechenko@gmail.com-copyright: 2019 Andrew Lelechenko+copyright: 2019-2020 Andrew Lelechenko category: Math, Numerical build-type: Simple-extra-source-files: README.md 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.1+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 extra-source-files:   changelog.md+  README.md  source-repository head   type: git@@ -25,8 +25,11 @@   hs-source-dirs: src   exposed-modules:     Data.Poly+    Data.Poly.Laurent+    Data.Poly.Orthogonal     Data.Poly.Semiring     Data.Poly.Sparse+    Data.Poly.Sparse.Laurent     Data.Poly.Sparse.Semiring   other-modules:     Data.Poly.Internal.Dense@@ -40,31 +43,36 @@     base >= 4.9 && < 5,     deepseq >= 1.1 && < 1.5,     primitive >= 0.6,-    semirings >= 0.2,+    semirings >= 0.5.2,     vector >= 0.12.0.2,-    vector-algorithms >= 0.7+    vector-algorithms >= 0.8.0.3   default-language: Haskell2010-  ghc-options: -Wall+  ghc-options: -Wall -Wcompat  test-suite poly-tests   type: exitcode-stdio-1.0   main-is: Main.hs   other-modules:     Dense+    DenseLaurent+    Orthogonal     Quaternion     Sparse+    SparseLaurent+    TestUtils   build-depends:     base >=4.9 && <5,+    mod,     poly,     QuickCheck >=2.12,-    quickcheck-classes >=0.5,-    semirings >= 0.2,+    quickcheck-classes >=0.6.3,+    semirings >= 0.5.2,     tasty >= 0.11,     tasty-quickcheck >= 0.8,     vector >= 0.12.0.2   default-language: Haskell2010   hs-source-dirs: test-  ghc-options: -Wall+  ghc-options: -Wall -Wcompat -threaded -rtsopts  benchmark poly-gauge   build-depends:@@ -81,4 +89,4 @@     SparseBench   default-language: Haskell2010   hs-source-dirs: bench-  ghc-options: -Wall+  ghc-options: -Wall -Wcompat
src/Data/Poly.hs view
@@ -7,7 +7,6 @@ -- Dense polynomials and a 'Num'-based interface. -- -{-# LANGUAGE CPP             #-} {-# LANGUAGE PatternSynonyms #-}  module Data.Poly@@ -16,7 +15,6 @@   , UPoly   , unPoly   , leading-  -- * Num interface   , toPoly   , monomial   , scale@@ -25,19 +23,10 @@   , subst   , deriv   , integral-#if MIN_VERSION_semirings(0,4,2)-  -- * Polynomials over 'Field'   , PolyOverField(..)-  , gcdExt-  , PolyOverFractional-  , pattern PolyOverFractional-  , unPolyOverFractional-#endif   ) where  import Data.Poly.Internal.Dense-#if MIN_VERSION_semirings(0,4,2)-import Data.Poly.Internal.Dense.Field (gcdExt)+import Data.Poly.Internal.Dense.Field () import Data.Poly.Internal.Dense.GcdDomain () import Data.Poly.Internal.PolyOverField-#endif
src/Data/Poly/Internal/Dense.hs view
@@ -7,7 +7,6 @@ -- Dense polynomials of one variable. -- -{-# LANGUAGE CPP                        #-} {-# LANGUAGE FlexibleInstances          #-} {-# LANGUAGE GeneralizedNewtypeDeriving #-} {-# LANGUAGE PatternSynonyms            #-}@@ -37,11 +36,10 @@   , pattern X'   , eval'   , subst'+  , substitute'   , deriv'-#if MIN_VERSION_semirings(0,5,0)   , unscale'   , integral'-#endif   ) where  import Prelude hiding (quotRem, quot, rem, gcd, lcm, (^))@@ -50,6 +48,7 @@ import Control.Monad.Primitive import Control.Monad.ST import Data.Bits+import Data.Euclidean (Euclidean, Field, quot) import Data.List (foldl', intersperse) import Data.Semiring (Semiring(..), Ring()) import qualified Data.Semiring as Semiring@@ -58,13 +57,6 @@ import qualified Data.Vector.Generic.Mutable as MG import qualified Data.Vector.Unboxed as U import GHC.Exts-#if !MIN_VERSION_semirings(0,4,0)-import Data.Semigroup-import Numeric.Natural-#endif-#if MIN_VERSION_semirings(0,5,0)-import Data.Euclidean (Euclidean, Field, quot)-#endif  -- | Polynomials of one variable with coefficients from @a@, -- backed by a 'G.Vector' @v@ (boxed, unboxed, storable, etc.).@@ -171,13 +163,11 @@   {-# INLINE plus #-}   {-# INLINE times #-} -#if MIN_VERSION_semirings(0,4,0)   fromNatural n = if n' == zero then zero else Poly $ G.singleton n'     where       n' :: a       n' = fromNatural n   {-# INLINE fromNatural #-}-#endif  instance (Eq a, Ring a, G.Vector v a) => Ring (Poly v a) where   negate (Poly xs) = Poly $ G.map Semiring.negate xs@@ -365,7 +355,6 @@ scale' :: (Eq a, Semiring a, G.Vector v a) => Word -> a -> Poly v a -> Poly v a scale' yp yc (Poly xs) = toPoly' $ scaleInternal zero times yp yc xs -#if MIN_VERSION_semirings(0,5,0) unscale' :: (Eq a, Euclidean a, G.Vector v a) => Word -> a -> Poly v a -> Poly v a unscale' yp yc (Poly xs) = toPoly' $ runST $ do   let lenZs = G.length xs - fromIntegral yp@@ -374,7 +363,6 @@     MG.unsafeWrite zs k (G.unsafeIndex xs (k + fromIntegral yp) `quot` yc)   G.unsafeFreeze zs {-# INLINABLE unscale' #-}-#endif  data StrictPair a b = !a :*: !b @@ -433,17 +421,6 @@   | otherwise = toPoly' $ G.imap (\i x -> fromNatural (fromIntegral (i + 1)) `times` x) $ G.tail xs {-# INLINE deriv' #-} -#if !MIN_VERSION_semirings(0,4,0)-fromNatural :: Semiring a => Natural -> a-fromNatural 0 = zero-fromNatural n = getAdd' (stimes n (Add' one))--newtype Add' a = Add' { getAdd' :: a }--instance Semiring a => Semigroup (Add' a) where-  Add' a <> Add' b = Add' (a `plus` b)-#endif- -- | Compute an indefinite integral of a polynomial, -- setting constant term to zero. --@@ -462,7 +439,6 @@       lenXs = G.length xs {-# INLINABLE integral #-} -#if MIN_VERSION_semirings(0,5,0) integral' :: (Eq a, Field a, G.Vector v a) => Poly v a -> Poly v a integral' (Poly xs)   | G.null xs = Poly G.empty@@ -475,7 +451,6 @@     where       lenXs = G.length xs {-# INLINABLE integral' #-}-#endif  -- | Create an identity polynomial. pattern X :: (Eq a, Num a, G.Vector v a, Eq (v a)) => Poly v a
src/Data/Poly/Internal/Dense/Field.hs view
@@ -8,7 +8,6 @@ --  {-# LANGUAGE ConstraintKinds            #-}-{-# LANGUAGE CPP                        #-} {-# LANGUAGE FlexibleInstances          #-} {-# LANGUAGE PatternSynonyms            #-} {-# LANGUAGE ScopedTypeVariables        #-}@@ -16,24 +15,17 @@  {-# OPTIONS_GHC -fno-warn-orphans #-} -#if MIN_VERSION_semirings(0,4,2)- module Data.Poly.Internal.Dense.Field   ( fieldGcd-  , gcdExt   ) where -import Prelude hiding (quotRem, quot, rem, gcd)+import Prelude hiding (quotRem, quot, rem, gcd, recip) import Control.Exception import Control.Monad import Control.Monad.Primitive import Control.Monad.ST-import Data.Euclidean (Euclidean(..))-#if !MIN_VERSION_semirings(0,5,0)-import Data.Semiring (Ring)-#else-import Data.Euclidean (Field)-#endif+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@@ -41,10 +33,6 @@ import Data.Poly.Internal.Dense import Data.Poly.Internal.Dense.GcdDomain () -#if !MIN_VERSION_semirings(0,5,0)-type Field a = (Euclidean a, Ring a, Fractional a)-#endif- instance (Eq a, Eq (v a), Field a, G.Vector v a) => Euclidean (Poly v a) where   degree (Poly xs) = fromIntegral (G.length xs) @@ -57,32 +45,40 @@   {-# INLINE rem #-}  quotientAndRemainder-  :: (Field a, G.Vector v a)+  :: (Eq a, Field a, G.Vector v a)   => v a   -> v a   -> (v a, v a) quotientAndRemainder xs ys-  | G.null ys = throw DivideByZero-  | G.length xs < G.length ys = (G.empty, xs)+  | 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)   | otherwise = runST $ do-    let lenXs = G.length xs-        lenYs = G.length ys-        lenQs = lenXs - lenYs + 1     qs <- MG.unsafeNew lenQs     rs <- MG.unsafeNew lenXs     G.unsafeCopy rs xs+    let yLast = G.unsafeLast ys+        invYLast = recip yLast     forM_ [lenQs - 1, lenQs - 2 .. 0] $ \i -> do       r <- MG.unsafeRead rs (lenYs - 1 + i)-      let q = r `quot` G.unsafeLast ys+      let q = if yLast == one then r else r `times` invYLast       MG.unsafeWrite qs i q-      forM_ [0 .. lenYs - 1] $ \k -> do-        MG.unsafeModify rs (\c -> c `minus` q `times` G.unsafeIndex ys k) (i + k)+      MG.unsafeWrite rs (lenYs - 1 + i) zero+      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)     let rs' = MG.unsafeSlice 0 lenYs rs     (,) <$> G.unsafeFreeze qs <*> G.unsafeFreeze rs'+  where+    lenXs = G.length xs+    lenYs = G.length ys+    lenQs = lenXs - lenYs + 1 {-# INLINABLE quotientAndRemainder #-}  remainder-  :: (Field a, G.Vector v a)+  :: (Eq a, Field a, G.Vector v a)   => v a   -> v a   -> v a@@ -96,25 +92,29 @@ {-# INLINABLE remainder #-}  remainderM-  :: (PrimMonad m, Field a, G.Vector v a)+  :: (PrimMonad m, Eq a, Field a, G.Vector v a)   => G.Mutable v (PrimState m) a   -> G.Mutable v (PrimState m) a   -> m () remainderM xs ys-  | MG.null ys = throw DivideByZero-  | MG.length xs < MG.length ys = pure ()+  | lenXs < lenYs = pure ()+  | lenYs == 0 = throw DivideByZero+  | lenYs == 1 = MG.set xs zero   | otherwise = do-    let lenXs = MG.length xs-        lenYs = MG.length ys-        lenQs = lenXs - lenYs + 1     yLast <- MG.unsafeRead ys (lenYs - 1)+    let invYLast = recip yLast     forM_ [lenQs - 1, lenQs - 2 .. 0] $ \i -> do       r <- MG.unsafeRead xs (lenYs - 1 + i)-      forM_ [0 .. lenYs - 1] $ \k -> do+      MG.unsafeWrite xs (lenYs - 1 + i) zero+      let q = if yLast == one then r else r `times` invYLast+      forM_ [0 .. lenYs - 2] $ \k -> do         y <- MG.unsafeRead ys k-        -- do not move r / yLast outside the loop,-        -- because of numerical instability-        MG.unsafeModify xs (\c -> c `minus` r `times` y `quot` yLast) (i + k)+        when (y /= zero) $+          MG.unsafeModify xs (\c -> c `minus` q `times` y) (i + k)+  where+    lenXs = MG.length xs+    lenYs = MG.length ys+    lenQs = lenXs - lenYs + 1 {-# INLINABLE remainderM #-}  fieldGcd@@ -139,53 +139,3 @@     remainderM xs ys'     gcdM ys' xs {-# INLINE gcdM #-}---- | Execute the extended Euclidean algorithm.--- For polynomials @a@ and @b@, compute their unique greatest common divisor @g@--- and the unique coefficient polynomial @s@ satisfying @as + bt = g@,--- such that either @g@ is monic, or @g = 0@ and @s@ is monic, or @g = s = 0@.------ >>> gcdExt (X^2 + 1 :: UPoly Double) (X^3 + 3 * X :: UPoly Double)--- (1.0, 0.5 * X^2 + (-0.0) * X + 1.0)--- >>> gcdExt (X^3 + 3 * X :: UPoly Double) (3 * X^4 + 3 * X^2 :: UPoly Double)--- (1.0 * X + 0.0,(-0.16666666666666666) * X^2 + (-0.0) * X + 0.3333333333333333)-gcdExt-  :: (Eq a, Field a, G.Vector v a, Eq (v a))-  => Poly v a-  -> Poly v a-  -> (Poly v a, Poly v a)-gcdExt xs ys = case scaleMonic gs of-  Just (c', gs') -> (gs', scale' zero c' ss)-  Nothing -> case scaleMonic ss of-    Just (_, ss') -> (zero, ss')-    Nothing -> (zero, zero)-  where-    (gs, ss) = go ys xs zero one-      where-        go r' r s' s-          | r' == zero = (r, s)-          | otherwise  = case r `quotRem` r' of-            (q, r'') -> go r'' r' (s `minus` q `times` s') s'-{-# INLINABLE gcdExt #-}---- | Scale a non-zero polynomial such that its leading coefficient is one,--- returning the reciprocal of the leading coefficient in the scaling.------ >>> scaleMonic (X^3 + 3 * X :: UPoly Double)--- Just (1.0, 1.0 * X^3 + 0.0 * X^2 + 3.0 * X + 0.0)--- >>> scaleMonic (3 * X^4 + 3 * X^2 :: UPoly Double)--- Just (0.3333333333333333, 1.0 * X^4 + 0.0 * X^3 + 1.0 * X^2 + 0.0 * X + 0.0)-scaleMonic-  :: (Eq a, Field a, G.Vector v a, Eq (v a))-  => Poly v a-  -> Maybe (a, Poly v a)-scaleMonic xs = case leading xs of-  Nothing -> Nothing-  Just (_, c) -> let c' = one `quot` c in Just (c', scale' zero c' xs)-{-# INLINE scaleMonic #-}--#else--module Data.Poly.Internal.Dense.Field () where--#endif
src/Data/Poly/Internal/Dense/GcdDomain.hs view
@@ -7,7 +7,6 @@ -- GcdDomain for GcdDomain underlying -- -{-# LANGUAGE CPP                        #-} {-# LANGUAGE FlexibleInstances          #-} {-# LANGUAGE PatternSynonyms            #-} {-# LANGUAGE ScopedTypeVariables        #-}@@ -18,8 +17,6 @@ module Data.Poly.Internal.Dense.GcdDomain   () where -#if MIN_VERSION_semirings(0,4,2)- import Prelude hiding (gcd, lcm, (^)) import Control.Exception import Control.Monad@@ -170,5 +167,3 @@      go (lenQs - 1) {-# INLINABLE quotient #-}--#endif
src/Data/Poly/Internal/PolyOverField.hs view
@@ -7,19 +7,13 @@ -- Wrapper with a more efficient 'Euclidean' instance. -- -{-# LANGUAGE CPP                        #-} {-# LANGUAGE ConstraintKinds            #-} {-# LANGUAGE FlexibleInstances          #-} {-# LANGUAGE GeneralizedNewtypeDeriving #-} {-# LANGUAGE PatternSynonyms            #-} -#if MIN_VERSION_semirings(0,4,2)- module Data.Poly.Internal.PolyOverField   ( PolyOverField(..)-  , PolyOverFractional-  , pattern PolyOverFractional-  , unPolyOverFractional   ) where  import Prelude hiding (quotRem, quot, rem, gcd, lcm, (^))@@ -36,23 +30,6 @@ newtype PolyOverField poly = PolyOverField { unPolyOverField :: poly }   deriving (Eq, NFData, Num, Ord, Ring, Semiring, Show) --- |-type PolyOverFractional = PolyOverField-{-# DEPRECATED PolyOverFractional "Use 'PolyOverField'" #-}---- |-pattern PolyOverFractional :: poly -> PolyOverField poly-pattern PolyOverFractional poly = PolyOverField poly---- |-unPolyOverFractional :: PolyOverField poly -> poly-unPolyOverFractional = unPolyOverField-{-# DEPRECATED unPolyOverFractional "Use 'unPolyOverField'" #-}--#if !MIN_VERSION_semirings(0,5,0)-type Field a = (Euclidean a, Ring a, Fractional a)-#endif- 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 #-}@@ -67,9 +44,3 @@   rem (PolyOverField x) (PolyOverField y) =     PolyOverField (rem x y)   {-# INLINE rem #-}--#else--module Data.Poly.Internal.PolyOverField () where--#endif
src/Data/Poly/Internal/Sparse.hs view
@@ -7,7 +7,6 @@ -- Sparse polynomials of one variable. -- -{-# LANGUAGE CPP                        #-} {-# LANGUAGE FlexibleContexts           #-} {-# LANGUAGE GeneralizedNewtypeDeriving #-} {-# LANGUAGE PatternSynonyms            #-}@@ -38,10 +37,9 @@   , pattern X'   , eval'   , subst'+  , substitute'   , deriv'-#if MIN_VERSION_semirings(0,5,0)   , integral'-#endif   ) where  import Prelude hiding (quot)@@ -50,6 +48,7 @@ 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())@@ -60,13 +59,6 @@ import qualified Data.Vector.Unboxed as U import qualified Data.Vector.Algorithms.Tim as Tim import GHC.Exts-#if !MIN_VERSION_semirings(0,4,0)-import Data.Semigroup-import Numeric.Natural-#endif-#if MIN_VERSION_semirings(0,5,0)-import Data.Euclidean (Field, quot)-#endif  -- | Polynomials of one variable with coefficients from @a@, -- backed by a 'G.Vector' @v@ (boxed, unboxed, storable, etc.).@@ -215,13 +207,11 @@   {-# INLINE plus #-}   {-# INLINE times #-} -#if MIN_VERSION_semirings(0,4,0)   fromNatural n = if n' == zero then zero else Poly $ G.singleton (0, n')     where       n' :: a       n' = fromNatural n   {-# INLINE fromNatural #-}-#endif  instance (Eq a, Ring a, G.Vector v (Word, a)) => Ring (Poly v a) where   negate (Poly xs) = Poly $ G.map (fmap Semiring.negate) xs@@ -528,17 +518,6 @@   xs {-# INLINE deriv' #-} -#if !MIN_VERSION_semirings(0,4,0)-fromNatural :: Semiring a => Natural -> a-fromNatural 0 = zero-fromNatural n = getAdd' (stimes n (Add' one))--newtype Add' a = Add' { getAdd' :: a }--instance Semiring a => Semigroup (Add' a) where-  Add' a <> Add' b = Add' (a `plus` b)-#endif- derivPoly   :: G.Vector v (Word, a)   => (a -> Bool)@@ -575,13 +554,11 @@   $ G.map (\(p, c) -> (p + 1, c / (fromIntegral p + 1))) xs {-# INLINE integral #-} -#if MIN_VERSION_semirings(0,5,0) 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' #-}-#endif  -- | Create an identity polynomial. pattern X :: (Eq a, Num a, G.Vector v (Word, a), Eq (v (Word, a))) => Poly v a
src/Data/Poly/Internal/Sparse/Field.hs view
@@ -8,7 +8,6 @@ --  {-# LANGUAGE ConstraintKinds            #-}-{-# LANGUAGE CPP                        #-} {-# LANGUAGE FlexibleContexts           #-} {-# LANGUAGE FlexibleInstances          #-} {-# LANGUAGE PatternSynonyms            #-}@@ -18,31 +17,18 @@  {-# OPTIONS_GHC -fno-warn-orphans #-} -#if MIN_VERSION_semirings(0,4,2)--module Data.Poly.Internal.Sparse.Field-  ( gcdExt-  ) where+module Data.Poly.Internal.Sparse.Field () where  import Prelude hiding (quotRem, quot, rem, gcd) import Control.Arrow import Control.Exception-import Data.Euclidean (Euclidean(..))-#if !MIN_VERSION_semirings(0,5,0)-import Data.Semiring (Ring)-#else-import Data.Euclidean (Field)-#endif-import Data.Semiring (minus, plus, times, zero, one)+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 () -#if !MIN_VERSION_semirings(0,5,0)-type Field a = (Euclidean a, Ring a, Fractional a)-#endif- 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@@ -68,53 +54,3 @@           where             zs = Poly $ G.singleton (xp `minus` yp, xc `quot` yc)             xs' = xs `minus` zs `times` ys---- | Execute the extended Euclidean algorithm.--- For polynomials @a@ and @b@, compute their unique greatest common divisor @g@--- and the unique coefficient polynomial @s@ satisfying @as + bt = g@,--- such that either @g@ is monic, or @g = 0@ and @s@ is monic, or @g = s = 0@.------ >>> gcdExt (X^2 + 1 :: UPoly Double) (X^3 + 3 * X :: UPoly Double)--- (1.0, 0.5 * X^2 + (-0.0) * X + 1.0)--- >>> gcdExt (X^3 + 3 * X :: UPoly Double) (3 * X^4 + 3 * X^2 :: UPoly Double)--- (1.0 * X + 0.0,(-0.16666666666666666) * X^2 + (-0.0) * X + 0.3333333333333333)-gcdExt-  :: (Eq a, Field a, G.Vector v (Word, a), Eq (v (Word, a)))-  => Poly v a-  -> Poly v a-  -> (Poly v a, Poly v a)-gcdExt xs ys = case scaleMonic gs of-  Just (c', gs') -> (gs', scale' zero c' ss)-  Nothing -> case scaleMonic ss of-    Just (_, ss') -> (zero, ss')-    Nothing -> (zero, zero)-  where-    (gs, ss) = go ys xs zero one-      where-        go r' r s' s-          | r' == zero = (r, s)-          | otherwise  = case r `quotRem` r' of-            (q, r'') -> go r'' r' (s `minus` q `times` s') s'-{-# INLINABLE gcdExt #-}---- | Scale a non-zero polynomial such that its leading coefficient is one,--- returning the reciprocal of the leading coefficient in the scaling.------ >>> scaleMonic (X^3 + 3 * X :: UPoly Double)--- Just (1.0, 1.0 * X^3 + 0.0 * X^2 + 3.0 * X + 0.0)--- >>> scaleMonic (3 * X^4 + 3 * X^2 :: UPoly Double)--- Just (0.3333333333333333, 1.0 * X^4 + 0.0 * X^3 + 1.0 * X^2 + 0.0 * X + 0.0)-scaleMonic-  :: (Eq a, Field a, G.Vector v (Word, a), Eq (v (Word, a)))-  => Poly v a-  -> Maybe (a, Poly v a)-scaleMonic xs = case leading xs of-  Nothing -> Nothing-  Just (_, c) -> let c' = one `quot` c in Just (c', scale' zero c' xs)-{-# INLINE scaleMonic #-}--#else--module Data.Poly.Internal.Sparse.Field () where--#endif
src/Data/Poly/Internal/Sparse/GcdDomain.hs view
@@ -7,7 +7,6 @@ -- GcdDomain for GcdDomain underlying -- -{-# LANGUAGE CPP                        #-} {-# LANGUAGE FlexibleContexts           #-} {-# LANGUAGE FlexibleInstances          #-} {-# LANGUAGE PatternSynonyms            #-}@@ -20,8 +19,6 @@ module Data.Poly.Internal.Sparse.GcdDomain   () where -#if MIN_VERSION_semirings(0,4,2)- import Prelude hiding (gcd, lcm, (^)) import Control.Exception import Data.Euclidean@@ -75,5 +72,3 @@         gx = fromMaybe err $ divide g xc         gy = fromMaybe err $ divide g yc         err = error "gcd: violated internal invariant"--#endif
+ src/Data/Poly/Laurent.hs view
@@ -0,0 +1,284 @@+-- |+-- Module:      Data.Poly.Laurent+-- Copyright:   (c) 2020 Andrew Lelechenko+-- Licence:     BSD3+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+--+-- <https://en.wikipedia.org/wiki/Laurent_polynomial Laurent polynomials>.+--++{-# LANGUAGE FlexibleInstances          #-}+{-# LANGUAGE GeneralizedNewtypeDeriving #-}+{-# LANGUAGE PatternSynonyms            #-}+{-# LANGUAGE ScopedTypeVariables        #-}+{-# LANGUAGE ViewPatterns               #-}++module Data.Poly.Laurent+  ( Laurent+  , VLaurent+  , ULaurent+  , unLaurent+  , toLaurent+  , leading+  , monomial+  , scale+  , pattern X+  , (^-)+  , eval+  , subst+  , deriv+  , LaurentOverField(..)+  ) where++import Prelude hiding (quotRem, quot, rem, gcd)+import Control.Arrow (first)+import Control.DeepSeq (NFData(..))+import Data.Euclidean (GcdDomain(..), Euclidean(..), Field)+import Data.List (intersperse)+import Data.Semiring (Semiring(..), Ring())+import qualified Data.Semiring as Semiring+import qualified Data.Vector as V+import qualified Data.Vector.Generic as G+import qualified Data.Vector.Unboxed as U++import Data.Poly.Internal.Dense (Poly(..))+import qualified Data.Poly.Internal.Dense as Dense+import Data.Poly.Internal.Dense.Field ()+import Data.Poly.Internal.Dense.GcdDomain ()+import Data.Poly.Internal.PolyOverField++-- | <https://en.wikipedia.org/wiki/Laurent_polynomial Laurent polynomials>+-- of one variable with coefficients from @a@,+-- backed by a 'G.Vector' @v@ (boxed, unboxed, storable, etc.).+--+-- Use pattern 'X' and operator '^-' for construction:+--+-- >>> (X + 1) + (X^-1 - 1) :: VLaurent Integer+-- 1 * X + 0 + 1 * X^-1+-- >>> (X + 1) * (1 - X^-1) :: ULaurent Int+-- 1 * X + 0 + (-1) * X^-1+--+-- Polynomials are stored normalized, without leading+-- and trailing+-- zero coefficients, so 0 * X + 1 + 0 * X^-1 equals to 1.+--+-- 'Ord' instance does not make much sense mathematically,+-- it is defined only for the sake of 'Data.Set.Set', 'Data.Map.Map', etc.+--+data Laurent v a = Laurent !Int !(Poly v a)+  deriving (Eq, Ord)++-- | Deconstruct a 'Laurent' polynomial into an offset (largest possible)+-- and a regular polynomial.+--+-- >>> unLaurent (2 * X + 1 :: ULaurent Int)+-- (0,2 * X + 1)+-- >>> unLaurent (1 + 2 * X^-1 :: ULaurent Int)+-- (-1,1 * X + 2)+-- >>> unLaurent (2 * X^2 + X :: ULaurent Int)+-- (1,2 * X + 1)+-- >>> unLaurent (0 :: ULaurent Int)+-- (0,0)+unLaurent :: Laurent v a -> (Int, Poly v a)+unLaurent (Laurent off poly) = (off, poly)++-- | Construct 'Laurent' polynomial from an offset and a regular polynomial.+-- One can imagine it as 'Data.Poly.scale'', but allowing negative offsets.+--+-- >>> toLaurent 2 (2 * Data.Poly.X + 1) :: ULaurent Int+-- 2 * X^3 + 1 * X^2+-- >>> toLaurent (-2) (2 * Data.Poly.X + 1) :: ULaurent Int+-- 2 * X^-1 + 1 * X^-2+toLaurent+  :: (Eq a, Semiring a, G.Vector v a)+  => Int+  -> Poly v a+  -> Laurent v a+toLaurent off (Poly xs) = go 0+  where+    go k+      | k >= G.length xs+      = Laurent 0 zero+      | G.unsafeIndex xs k == zero+      = go (k + 1)+      | otherwise+      = Laurent (off + k) (Poly (G.unsafeDrop k xs))+{-# INLINE toLaurent #-}++toLaurentNum+  :: (Eq a, Num a, G.Vector v a)+  => Int+  -> Poly v a+  -> Laurent v a+toLaurentNum off (Poly xs) = go 0+  where+    go k+      | k >= G.length xs+      = Laurent 0 0+      | G.unsafeIndex xs k == 0+      = go (k + 1)+      | otherwise+      = Laurent (off + k) (Poly (G.unsafeDrop k xs))+{-# INLINE toLaurentNum #-}++instance NFData (v a) => NFData (Laurent v a) where+  rnf (Laurent off poly) = rnf off `seq` rnf poly++instance (Show a, G.Vector v a) => Show (Laurent v a) where+  showsPrec d (Laurent off poly)+    | G.null (unPoly poly)+      = showString "0"+    | otherwise+      = showParen (d > 0)+      $ foldl (.) id+      $ intersperse (showString " + ")+      $ G.ifoldl (\acc i c -> showCoeff (i + off) c : acc) []+      $ unPoly poly+    where+      showCoeff 0 c = showsPrec 7 c+      showCoeff 1 c = showsPrec 7 c . showString " * X"+      showCoeff i c = showsPrec 7 c . showString (" * X^" ++ show i)++-- | Laurent polynomials backed by boxed vectors.+type VLaurent = Laurent V.Vector++-- | Laurent polynomials backed by unboxed vectors.+type ULaurent = Laurent U.Vector++-- | Return a leading power and coefficient of a non-zero polynomial.+--+-- >>> leading ((2 * X + 1) * (2 * X^2 - 1) :: ULaurent Int)+-- Just (3,4)+-- >>> leading (0 :: ULaurent Int)+-- Nothing+leading :: G.Vector v a => Laurent v a -> Maybe (Int, a)+leading (Laurent off poly) = first ((+ off) . fromIntegral) <$> Dense.leading poly++-- | Note that 'abs' = 'id' and 'signum' = 'const' 1.+instance (Eq a, Num a, G.Vector v a) => Num (Laurent v a) where+  Laurent off1 poly1 * Laurent off2 poly2 = toLaurentNum (off1 + off2) (poly1 * poly2)+  Laurent off1 poly1 + Laurent off2 poly2 = case off1 `compare` off2 of+    LT -> toLaurentNum off1 (poly1 + Dense.scale (fromIntegral $ off2 - off1) 1 poly2)+    EQ -> toLaurentNum off1 (poly1 + poly2)+    GT -> toLaurentNum off2 (Dense.scale (fromIntegral $ off1 - off2) 1 poly1 + poly2)+  Laurent off1 poly1 - Laurent off2 poly2 = case off1 `compare` off2 of+    LT -> toLaurentNum off1 (poly1 - Dense.scale (fromIntegral $ off2 - off1) 1 poly2)+    EQ -> toLaurentNum off1 (poly1 - poly2)+    GT -> toLaurentNum off2 (Dense.scale (fromIntegral $ off1 - off2) 1 poly1 - poly2)+  negate (Laurent off poly) = Laurent off (negate poly)+  abs = id+  signum = const 1+  fromInteger n = Laurent 0 (fromInteger n)+  {-# INLINE (+) #-}+  {-# INLINE (-) #-}+  {-# INLINE negate #-}+  {-# INLINE fromInteger #-}+  {-# INLINE (*) #-}++instance (Eq a, Semiring a, G.Vector v a) => Semiring (Laurent v a) where+  zero = Laurent 0 zero+  one  = Laurent 0 one+  Laurent off1 poly1 `times` Laurent off2 poly2 =+    toLaurent (off1 + off2) (poly1 `times` poly2)+  Laurent off1 poly1 `plus` Laurent off2 poly2 = case off1 `compare` off2 of+    LT -> toLaurent off1 (poly1 `plus` Dense.scale' (fromIntegral $ off2 - off1) one poly2)+    EQ -> toLaurent off1 (poly1 `plus` poly2)+    GT -> toLaurent off2 (Dense.scale' (fromIntegral $ off1 - off2) one poly1 `plus` poly2)+  fromNatural n = Laurent 0 (fromNatural n)+  {-# INLINE zero #-}+  {-# INLINE one #-}+  {-# INLINE plus #-}+  {-# INLINE times #-}+  {-# INLINE fromNatural #-}++instance (Eq a, Ring a, G.Vector v a) => Ring (Laurent v a) where+  negate (Laurent off poly) = Laurent off (Semiring.negate poly)++-- | Create a monomial from a power and a coefficient.+monomial :: (Eq a, Semiring a, G.Vector v a) => Int -> a -> Laurent v a+monomial p c+  | c == zero = Laurent 0 zero+  | otherwise = Laurent p (Dense.monomial' 0 c)+{-# INLINE monomial #-}++-- | Multiply a polynomial by a monomial, expressed as a power and a coefficient.+--+-- >>> scale 2 3 (X^2 + 1) :: ULaurent Int+-- 3 * X^4 + 0 * X^3 + 3 * X^2 + 0 * X + 0+scale :: (Eq a, Semiring a, G.Vector v a) => Int -> a -> Laurent v a -> Laurent v a+scale yp yc (Laurent off poly) = toLaurent (off + yp) (Dense.scale' 0 yc poly)++-- | Evaluate at a given point.+--+-- >>> eval (X^2 + 1 :: ULaurent Int) 3+-- 10+eval :: (Field a, G.Vector v a) => Laurent v a -> a -> a+eval (Laurent off poly) x = Dense.eval' poly x `times`+  (if off >= 0 then x Semiring.^ off else quot one x Semiring.^ (- off))+{-# INLINE eval #-}++-- | Substitute another polynomial instead of 'Data.Poly.X'.+--+-- >>> subst (X^2 + 1 :: UPoly Int) (X + 1 :: ULaurent Int)+-- 1 * X^2 + 2 * X + 2+subst :: (Eq a, Semiring a, G.Vector v a, G.Vector w a) => Poly v a -> Laurent w a -> Laurent w a+subst = Dense.substitute' (scale 0)+{-# INLINE subst #-}++-- | Take a derivative.+--+-- >>> deriv (X^3 + 3 * X) :: ULaurent Int+-- 3 * X^2 + 0 * X + 3+deriv :: (Eq a, Ring a, G.Vector v a) => Laurent v a -> Laurent v a+deriv (Laurent off (Poly xs)) =+  toLaurent (off - 1) $ Dense.toPoly' $ G.imap (times . Semiring.fromIntegral . (+ off)) xs+{-# INLINE deriv #-}++-- | Create an identity polynomial.+pattern X :: (Eq a, Semiring a, G.Vector v a, Eq (v a)) => Laurent v a+pattern X <- ((==) var -> True)+  where X = var++var :: forall a v. (Eq a, Semiring a, G.Vector v a, Eq (v a)) => Laurent v a+var+  | (one :: a) == zero = Laurent 0 zero+  | otherwise          = Laurent 1 one+{-# INLINE var #-}++-- | This operator can be applied only to 'X',+-- but is instrumental to express Laurent polynomials in mathematical fashion:+--+-- >>> X + 2 + 3 * X^-1 :: ULaurent Int+-- 1 * X + 2 + 3 * X^(-1)+(^-)+  :: (Eq a, Semiring a, G.Vector v a, Eq (v a))+  => Laurent v a+  -> Int+  -> Laurent v a+X^-n = monomial (negate n) one+_^-_ = error "(^-) can be applied only to X"++-- | Consider using 'LaurentOverField' wrapper,+-- which provides a much faster implementation of+-- 'Data.Euclidean.gcd' for polynomials over 'Field'.+instance (Eq a, Ring a, GcdDomain a, Eq (v a), G.Vector v a) => GcdDomain (Laurent v a) where+  divide (Laurent off1 poly1) (Laurent off2 poly2) =+    toLaurent (off1 - off2) <$> divide poly1 poly2+  {-# INLINE divide #-}++  gcd (Laurent _ poly1) (Laurent _ poly2) =+    toLaurent 0 (gcd poly1 poly2)+  {-# INLINE gcd #-}++-- | Wrapper for Laurent polynomials over 'Field',+-- providing a faster 'GcdDomain' instance.+newtype LaurentOverField laurent = LaurentOverField { unLaurentOverField :: laurent }+  deriving (Eq, NFData, Num, Ord, Ring, Semiring, Show)++instance (Eq a, Eq (v a), Field a, G.Vector v a) => GcdDomain (LaurentOverField (Laurent v a)) where+  divide (LaurentOverField (Laurent off1 poly1)) (LaurentOverField (Laurent off2 poly2)) =+    LaurentOverField . toLaurent (off1 - off2) . unPolyOverField <$> divide (PolyOverField poly1) (PolyOverField poly2)++  gcd (LaurentOverField (Laurent _ poly1)) (LaurentOverField (Laurent _ poly2)) =+    LaurentOverField (toLaurent 0 (unPolyOverField (gcd (PolyOverField poly1) (PolyOverField poly2))))+  {-# INLINE gcd #-}
+ src/Data/Poly/Orthogonal.hs view
@@ -0,0 +1,128 @@+-- |+-- Module:      Data.Poly.Orthogonal+-- Copyright:   (c) 2019 Andrew Lelechenko+-- Licence:     BSD3+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+--+-- Classical orthogonal polynomials.+--++{-# LANGUAGE OverloadedLists     #-}+{-# LANGUAGE RebindableSyntax    #-}++module Data.Poly.Orthogonal+  ( legendre+  , legendreShifted+  , gegenbauer+  , jacobi+  , chebyshev1+  , chebyshev2+  , hermiteProb+  , hermitePhys+  , laguerre+  , laguerreGen+  ) where++import Prelude hiding (quot, Num(..), fromIntegral)+import Data.Euclidean+import Data.Semiring+import Data.Poly.Semiring+import Data.Poly.Internal.Dense (unscale')+import Data.Vector.Generic (Vector, fromListN)++-- | <https://en.wikipedia.org/wiki/Legendre_polynomials Legendre polynomials>.+--+-- >>> 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+  where+    subst' :: (Eq a, Semiring a, Vector v a) => Poly v a -> Poly v a -> Poly v a+    subst' = subst++-- | <https://en.wikipedia.org/wiki/Legendre_polynomials#Shifted_Legendre_polynomials Shifted Legendre polynomials>.+--+-- >>> take 3 legendreShifted :: [Data.Poly.VPoly Integer]+-- [1,2 * X + (-1),6 * X^2 + (-6) * X + 1]+legendreShifted :: (Eq a, Euclidean a, Ring a, Vector v a) => [Poly v a]+legendreShifted = xs+  where+    xs = 1 : toPoly [-1, 2] : zipWith3 rec (iterate (+ 1) 1) xs (tail xs)+    rec n pm1 p = unscale' 0 (n + 1) (toPoly [-1 - 2 * n, 2 + 4 * n] * p - scale 0 n pm1)++-- | <https://en.wikipedia.org/wiki/Gegenbauer_polynomials Gegenbauer polynomials>.+gegenbauer :: (Eq a, Field a, Vector v a) => a -> [Poly v a]+gegenbauer g = jacobi a a+  where+    a = g - 1 `quot` 2++-- | <https://en.wikipedia.org/wiki/Jacobi_polynomials Jacobi polynomials>.+jacobi :: (Eq a, Field a, Vector v a) => a -> a -> [Poly v a]+jacobi a b = xs+  where+    x0 = 1+    x1 = toPoly [(a - b) `quot` 2, (a + b + 2) `quot` 2]+    xs = x0 : x1 : zipWith3 rec (iterate (+ 1) 2) xs (tail xs)+    rec n pm1 p = toPoly [d, c] * p - scale 0 cm1 pm1+      where+        cp1 = 2 * n * (n + a + b) * (2 * n + a + b - 2)+        q   = (2 * n + a + b - 1) `quot` cp1+        c   = q * ((2 * n + a + b) * (2 * n + a + b - 2))+        d   = q * (a * a - b * b)+        cm1 = 2 * (n + a - 1) * (n + b - 1) * (2 * n + a + b) `quot` cp1++-- | <https://en.wikipedia.org/wiki/Chebyshev_polynomials Chebyshev polynomials>+-- of the first kind.+--+-- >>> take 3 chebyshev1 :: [VPoly Integer]+-- [1,1 * X + 0,2 * X^2 + 0 * X + (-1)]+chebyshev1 :: (Eq a, Ring a, Vector v a) => [Poly v a]+chebyshev1 = xs+  where+    xs = 1 : monomial 1 1 : zipWith (\pm1 p -> scale 1 2 p - pm1) xs (tail xs)++-- | <https://en.wikipedia.org/wiki/Chebyshev_polynomials Chebyshev polynomials>+-- of the second kind.+--+-- >>> take 3 chebyshev2 :: [VPoly Integer]+-- [1,2 * X + 0,4 * X^2 + 0 * X + (-1)]+chebyshev2 :: (Eq a, Ring a, Vector v a) => [Poly v a]+chebyshev2 = xs+  where+    xs = 1 : monomial 1 2 : zipWith (\pm1 p -> scale 1 2 p - pm1) xs (tail xs)++-- | Probabilists' <https://en.wikipedia.org/wiki/Hermite_polynomials Hermite polynomials>.+--+-- >>> take 3 hermiteProb :: [VPoly Integer]+-- [1,1 * X + 0,1 * X^2 + 0 * X + (-1)]+hermiteProb :: (Eq a, Ring a, Vector v a) => [Poly v a]+hermiteProb = xs+  where+    xs = 1 : monomial 1 1 : zipWith3 rec (iterate (+ 1) 1) xs (tail xs)+    rec n pm1 p = scale 1 1 p - scale 0 n pm1++-- | Physicists' <https://en.wikipedia.org/wiki/Hermite_polynomials Hermite polynomials>.+--+-- >>> take 3 hermitePhys :: [VPoly Double]+-- [1,2 * X + 0,4 * X^2 + 0 * X + (-2)]+hermitePhys :: (Eq a, Ring a, Vector v a) => [Poly v a]+hermitePhys = xs+  where+    xs = 1 : monomial 1 2 : zipWith3 rec (iterate (+ 1) 1) xs (tail xs)+    rec n pm1 p = scale 1 2 p - scale 0 (2 * n) pm1++-- | <https://en.wikipedia.org/wiki/Laguerre_polynomials Laguerre polynomials>.+--+-- >>> take 3 laguerre :: [VPoly Double]+-- [1.0,(-1.0) * X + 1.0,0.5 * X^2 + (-2.0) * X + 1.0]+laguerre :: (Eq a, Field a, Vector v a) => [Poly v a]+laguerre = laguerreGen 0++-- | <https://en.wikipedia.org/wiki/Laguerre_polynomials#Generalized_Laguerre_polynomials Generalized Laguerre polynomials>+laguerreGen :: (Eq a, Field a, Vector v a) => a -> [Poly v a]+laguerreGen a = xs+  where+    x0 = 1+    x1 = toPoly [1 + a, -1]+    xs = x0 : x1 : zipWith3 rec (iterate (+ 1) 1) xs (tail xs)+    rec n pm1 p = toPoly [(2 * n + 1 + a) `quot` (n + 1), -1 `quot` (n + 1)] * p - scale 0 ((n + a) `quot` (n + 1)) pm1
src/Data/Poly/Semiring.hs view
@@ -7,7 +7,6 @@ -- Dense polynomials and a 'Semiring'-based interface. -- -{-# LANGUAGE CPP             #-} {-# LANGUAGE PatternSynonyms #-}  module Data.Poly.Semiring@@ -16,7 +15,6 @@   , UPoly   , unPoly   , leading-  -- * Semiring interface   , toPoly   , monomial   , scale@@ -24,32 +22,19 @@   , eval   , subst   , deriv-#if MIN_VERSION_semirings(0,5,0)   , integral-#endif-#if MIN_VERSION_semirings(0,4,2)-  -- * Polynomials over 'Field'   , PolyOverField(..)-  , gcdExt-  , PolyOverFractional-  , pattern PolyOverFractional-  , unPolyOverFractional-#endif   ) where +import Data.Euclidean (Field) import Data.Semiring (Semiring) import qualified Data.Vector.Generic as G  import Data.Poly.Internal.Dense (Poly(..), VPoly, UPoly, leading) import qualified Data.Poly.Internal.Dense as Dense-#if MIN_VERSION_semirings(0,4,2)-import Data.Poly.Internal.Dense.Field (gcdExt)+import Data.Poly.Internal.Dense.Field () import Data.Poly.Internal.Dense.GcdDomain () import Data.Poly.Internal.PolyOverField-#endif-#if MIN_VERSION_semirings(0,5,0)-import Data.Euclidean (Field)-#endif  -- | Make 'Poly' from a vector of coefficients -- (first element corresponds to a constant term).@@ -98,7 +83,6 @@ deriv :: (Eq a, Semiring a, G.Vector v a) => Poly v a -> Poly v a deriv = Dense.deriv' -#if MIN_VERSION_semirings(0,5,0) -- | Compute an indefinite integral of a polynomial, -- setting constant term to zero. --@@ -106,4 +90,3 @@ -- 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'-#endif
src/Data/Poly/Sparse.hs view
@@ -7,7 +7,6 @@ -- Sparse polynomials with 'Num' instance. -- -{-# LANGUAGE CPP             #-} {-# LANGUAGE PatternSynonyms #-}  module Data.Poly.Sparse@@ -16,7 +15,6 @@   , UPoly   , unPoly   , leading-  -- * Num interface   , toPoly   , monomial   , scale@@ -25,14 +23,8 @@   , subst   , deriv   , integral-#if MIN_VERSION_semirings(0,4,2)-  -- * Polynomials over 'Field'-  , gcdExt-#endif   ) where  import Data.Poly.Internal.Sparse-#if MIN_VERSION_semirings(0,4,2)-import Data.Poly.Internal.Sparse.Field (gcdExt)+import Data.Poly.Internal.Sparse.Field () import Data.Poly.Internal.Sparse.GcdDomain ()-#endif
+ src/Data/Poly/Sparse/Laurent.hs view
@@ -0,0 +1,283 @@+-- |+-- Module:      Data.Poly.Sparse.Laurent+-- Copyright:   (c) 2020 Andrew Lelechenko+-- Licence:     BSD3+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+--+-- Sparse <https://en.wikipedia.org/wiki/Laurent_polynomial Laurent polynomials>.+--++{-# LANGUAGE FlexibleContexts           #-}+{-# LANGUAGE FlexibleInstances          #-}+{-# LANGUAGE GeneralizedNewtypeDeriving #-}+{-# LANGUAGE PatternSynonyms            #-}+{-# LANGUAGE ScopedTypeVariables        #-}+{-# LANGUAGE StandaloneDeriving         #-}+{-# LANGUAGE TypeFamilies               #-}+{-# LANGUAGE UndecidableInstances       #-}+{-# LANGUAGE ViewPatterns               #-}++module Data.Poly.Sparse.Laurent+  ( Laurent+  , VLaurent+  , ULaurent+  , unLaurent+  , toLaurent+  , leading+  , monomial+  , scale+  , pattern X+  , (^-)+  , eval+  , subst+  , 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 qualified Data.Vector.Generic as G+import qualified Data.Vector.Unboxed as U+import GHC.Exts++import Data.Poly.Internal.Sparse (Poly(..))+import qualified Data.Poly.Internal.Sparse as Sparse+import Data.Poly.Internal.Sparse.Field ()+import Data.Poly.Internal.Sparse.GcdDomain ()++-- | <https://en.wikipedia.org/wiki/Laurent_polynomial Laurent polynomials>+-- of one variable with coefficients from @a@,+-- backed by a 'G.Vector' @v@ (boxed, unboxed, storable, etc.).+--+-- Use pattern 'X' and operator '^-' for construction:+--+-- >>> (X + 1) + (X^-1 - 1) :: VLaurent Integer+-- 1 * X + 1 * X^-1+-- >>> (X + 1) * (1 - X^-1) :: ULaurent Int+-- 1 * X + (-1) * X^-1+--+-- Polynomials are stored normalized, without+-- zero coefficients, so 0 * X + 1 + 0 * X^-1 equals to 1.+--+-- 'Ord' instance does not make much sense mathematically,+-- it is defined only for the sake of 'Data.Set.Set', 'Data.Map.Map', etc.+--+data Laurent v a = Laurent !Int !(Poly v a)++deriving instance Eq  (v (Word, a)) => Eq  (Laurent v a)+deriving instance Ord (v (Word, a)) => Ord (Laurent v a)++instance (Eq a, Semiring a, G.Vector v (Word, a)) => IsList (Laurent v a) where+  type Item (Laurent v a) = (Int, a)++  fromList xs = toLaurent minPow (fromList ys)+    where+      minPow = minimum $ maxBound : map fst xs+      ys = map (first (fromIntegral . (subtract minPow))) xs++  fromListN n xs = toLaurent minPow (fromListN n ys)+    where+      minPow = minimum $ maxBound : map fst xs+      ys = map (first (fromIntegral . (subtract minPow))) xs++  toList (Laurent off poly) =+    map (first ((+ off) . fromIntegral)) $ G.toList $ unPoly poly++-- | 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.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))+  => Int+  -> Poly v 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 #-}++toLaurentNum+  :: (Eq a, Num a, G.Vector v (Word, a))+  => Int+  -> Poly v 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)++-- | 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 #-}++-- | Multiply a polynomial by a monomial, expressed as a power and a coefficient.+--+-- >>> 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)++-- | 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 #-}++-- | Substitute another polynomial instead of 'Data.Poly.Sparse.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 #-}++-- | 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)))+  => 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 #-}
src/Data/Poly/Sparse/Semiring.hs view
@@ -7,7 +7,6 @@ -- Sparse polynomials with 'Semiring' instance. -- -{-# LANGUAGE CPP              #-} {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE PatternSynonyms  #-} @@ -17,7 +16,6 @@   , UPoly   , unPoly   , leading-  -- * Semiring interface   , toPoly   , monomial   , scale@@ -25,27 +23,17 @@   , eval   , subst   , deriv-#if MIN_VERSION_semirings(0,5,0)   , integral-#endif-#if MIN_VERSION_semirings(0,4,2)-  -- * Polynomials over 'Field'-  , gcdExt-#endif   ) where +import Data.Euclidean (Field) import Data.Semiring (Semiring) import qualified Data.Vector.Generic as G  import Data.Poly.Internal.Sparse (Poly(..), VPoly, UPoly, leading) import qualified Data.Poly.Internal.Sparse as Sparse-#if MIN_VERSION_semirings(0,4,2)-import Data.Poly.Internal.Sparse.Field (gcdExt)+import Data.Poly.Internal.Sparse.Field () import Data.Poly.Internal.Sparse.GcdDomain ()-#endif-#if MIN_VERSION_semirings(0,5,0)-import Data.Euclidean (Field)-#endif  -- | Make 'Poly' from a list of (power, coefficient) pairs. -- (first element corresponds to a constant term).@@ -94,7 +82,6 @@ deriv :: (Eq a, Semiring a, G.Vector v (Word, a)) => Poly v a -> Poly v a deriv = Sparse.deriv' -#if MIN_VERSION_semirings(0,5,0) -- | Compute an indefinite integral of a polynomial, -- setting constant term to zero. --@@ -102,4 +89,3 @@ -- 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'-#endif
test/Dense.hs view
@@ -1,4 +1,4 @@-{-# LANGUAGE CPP                        #-}+{-# LANGUAGE DataKinds                  #-} {-# LANGUAGE FlexibleContexts           #-} {-# LANGUAGE FlexibleInstances          #-} {-# LANGUAGE GeneralizedNewtypeDeriving #-}@@ -8,14 +8,13 @@  module Dense   ( testSuite+  , ShortPoly(..)   ) where  import Prelude hiding (gcd, quotRem, rem)-#if MIN_VERSION_semirings(0,4,2) import Data.Euclidean (Euclidean(..), GcdDomain(..))-#endif import Data.Int-import Data.Maybe+import Data.Mod import Data.Poly import qualified Data.Poly.Semiring as S import Data.Proxy@@ -25,30 +24,20 @@ import qualified Data.Vector.Unboxed as U import Test.Tasty import Test.Tasty.QuickCheck hiding (scale, numTests)-import Test.QuickCheck.Classes  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 -#if MIN_VERSION_semirings(0,4,2) 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-#endif  newtype ShortPoly a = ShortPoly { unShortPoly :: a }-  deriving-    ( Eq-    , Show-    , Semiring-#if MIN_VERSION_semirings(0,4,2)-    , GcdDomain-    , Euclidean-#endif-    )+  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@@ -61,95 +50,62 @@     , lawsTests     , evalTests     , derivTests-#if MIN_VERSION_semirings(0,4,2)-    , gcdExtTests-#endif     ]  lawsTests :: TestTree lawsTests = testGroup "Laws"-  [ semiringTests-  , ringTests-  , numTests-  , euclideanTests-  , isListTests-  , showTests+  $ semiringTests ++ ringTests ++ numTests ++ euclideanTests ++ gcdDomainTests ++ isListTests ++ showTests++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)))   ] -semiringTests :: TestTree-semiringTests-  = testGroup "Semiring"-  $ map (uncurry testProperty)-  $ concatMap lawsProperties-  [ semiringLaws (Proxy :: Proxy (Poly U.Vector ()))-  , semiringLaws (Proxy :: Proxy (Poly U.Vector Int8))-  , semiringLaws (Proxy :: Proxy (Poly V.Vector Integer))-  , semiringLaws (Proxy :: Proxy (Poly U.Vector (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)))   ] -ringTests :: TestTree-ringTests-  = testGroup "Ring"-  $ map (uncurry testProperty)-  $ concatMap lawsProperties-  [-#if MIN_VERSION_quickcheck_classes(0,6,1)-    ringLaws (Proxy :: Proxy (Poly U.Vector ()))-  , ringLaws (Proxy :: Proxy (Poly U.Vector Int8))-  , ringLaws (Proxy :: Proxy (Poly V.Vector Integer))-  , ringLaws (Proxy :: Proxy (Poly U.Vector (Quaternion Int)))-#endif+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)))   ] -numTests :: TestTree-numTests-  = testGroup "Num"-  $ map (uncurry testProperty)-  $ concatMap lawsProperties-  [-#if MIN_VERSION_quickcheck_classes(0,6,3)-    numLaws (Proxy :: Proxy (Poly U.Vector Int8))-  , numLaws (Proxy :: Proxy (Poly V.Vector Integer))-  , numLaws (Proxy :: Proxy (Poly U.Vector (Quaternion Int)))-#endif+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)))   ] -euclideanTests :: TestTree-euclideanTests-  = testGroup "Euclidean"-  $ map (uncurry testProperty)-  $ concatMap lawsProperties-  [-#if MIN_VERSION_semirings(0,4,2) && MIN_VERSION_quickcheck_classes(0,6,3)-    gcdDomainLaws (Proxy :: Proxy (ShortPoly (Poly V.Vector Integer)))-  , gcdDomainLaws (Proxy :: Proxy (PolyOverField (Poly V.Vector Rational)))-  , euclideanLaws (Proxy :: Proxy (ShortPoly (Poly V.Vector Rational)))-#endif+euclideanTests :: [TestTree]+euclideanTests =+  [ myEuclideanLaws (Proxy :: Proxy (ShortPoly (Poly V.Vector (Mod 3))))+  , myEuclideanLaws (Proxy :: Proxy (ShortPoly (Poly V.Vector Rational)))   ] -isListTests :: TestTree-isListTests-  = testGroup "IsList"-  $ map (uncurry testProperty)-  $ concatMap lawsProperties-  [ isListLaws (Proxy :: Proxy (Poly U.Vector ()))-  , isListLaws (Proxy :: Proxy (Poly U.Vector Int8))-  , isListLaws (Proxy :: Proxy (Poly V.Vector Integer))-  , isListLaws (Proxy :: Proxy (Poly U.Vector (Quaternion Int)))+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)))   ] -showTests :: TestTree-showTests-  = testGroup "Show"-  $ map (uncurry testProperty)-  $ concatMap lawsProperties-  [-#if MIN_VERSION_quickcheck_classes(0,6,0)-    showLaws (Proxy :: Proxy (Poly U.Vector ()))-  , showLaws (Proxy :: Proxy (Poly U.Vector Int8))-  , showLaws (Proxy :: Proxy (Poly V.Vector Integer))-  , showLaws (Proxy :: Proxy (Poly U.Vector (Quaternion Int)))-#endif+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)))   ]  arithmeticTests :: TestTree@@ -199,7 +155,8 @@     \p c -> c /= 0 ==> leading (monomial p c :: UPoly a) === Just (p, c)   , testProperty "monomial matches reference" $     \p (c :: a) -> monomial p c === toPoly (V.fromList (monomialRef p c))-  , testProperty "scale matches multiplication by monomial" $+  , tenTimesLess $+    testProperty "scale matches multiplication by monomial" $     \p c (xs :: UPoly a) -> scale p c xs === monomial p c * xs   ] @@ -250,13 +207,15 @@   => Proxy (Poly v a)   -> [TestTree] substTestGroup _ =-  [ testProperty "subst (p + q) r = subst p r + subst q r" $+  [ 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   , 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-  , testProperty "subst' (p + q) r = subst' p r + subst' q r" $+  , 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   , testProperty "subst' x p = p" $     \p -> e' S.X p === p@@ -273,10 +232,8 @@ derivTests = testGroup "deriv"   [ testProperty "deriv = S.deriv" $     \(p :: Poly V.Vector Integer) -> deriv p === S.deriv p-#if MIN_VERSION_semirings(0,5,0)   , testProperty "integral = S.integral" $     \(p :: Poly V.Vector Rational) -> integral p === S.integral p-#endif   , testProperty "deriv . integral = id" $     \(p :: Poly V.Vector Rational) -> deriv (integral p) === p   , testProperty "deriv c = 0" $@@ -287,26 +244,8 @@     \p q -> deriv (p + q) === (deriv p + deriv q :: Poly V.Vector 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" $+  , tenTimesLess $ tenTimesLess $ tenTimesLess $+    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   ]--#if MIN_VERSION_semirings(0,4,2)-gcdExtTests :: TestTree-gcdExtTests = localOption (QuickCheckMaxSize 12) $ testGroup "gcdExt"-  [ testProperty "gcdExt == S.gcdExt" $-    \(a :: Poly V.Vector Rational) b ->-      gcdExt a b === S.gcdExt a b-  , testProperty "g == as (mod b) for gcdExt" $-    \(a :: Poly V.Vector Rational) b -> b /= 0 ==>-      uncurry ((. flip rem b) . (===) . flip rem b) ((* a) <$> gcdExt a b)-  , testProperty "fst . gcdExt == gcd (mod units)" $-    \(a :: Poly V.Vector Rational) b ->-      fst (gcdExt a b) `sameUpToUnits` gcd a b-  ]--sameUpToUnits :: (Eq a, GcdDomain a) => a -> a -> Bool-sameUpToUnits x y = x == y ||-  isJust (x `divide` y) && isJust (y `divide` x)-#endif
+ test/DenseLaurent.hs view
@@ -0,0 +1,171 @@+{-# 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 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+  , lawsTests+  , evalTests+  , derivTests+  ]++lawsTests :: TestTree+lawsTests = testGroup "Laws"+  $ semiringTests ++ ringTests ++ numTests ++ gcdDomainTests ++ showTests++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)))+  ]++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)))+  ]++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)))+  ]++gcdDomainTests :: [TestTree]+gcdDomainTests =+  [ myGcdDomainLaws (Proxy :: Proxy (ShortLaurent (Laurent V.Vector Integer)))+  , myGcdDomainLaws (Proxy :: Proxy (LaurentOverField (Laurent V.Vector 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)))+  ]++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 "leading p 0 == Nothing" $+    \p -> leading (monomial p 0 :: ULaurent a) === Nothing+  , testProperty "leading . monomial = id" $+    \p c -> c /= 0 ==> leading (monomial p c :: ULaurent a) === Just (p, c)+  , tenTimesLess $+    testProperty "scale matches multiplication by monomial" $+    \p c (xs :: ULaurent a) -> scale p c xs === monomial p c * xs+  ]++evalTests :: TestTree+evalTests = testGroup "eval" $ concat+  [ evalTestGroup  (Proxy :: Proxy (Laurent V.Vector Rational))+  , substTestGroup (Proxy :: Proxy (Laurent U.Vector Int8))+  ]++evalTestGroup+  :: forall v a.+     (Eq a, Field a, Arbitrary a, Show a, Eq (v a), Show (v a), G.Vector v 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+  , 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+  , 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+  ]+  where+    e :: Laurent v a -> a -> a+    e = eval++substTestGroup+  :: forall v a.+     (Eq a, Num a, Semiring a, Arbitrary a, Show a, Eq (v a), Show (v a), G.Vector v a)+  => Proxy (Laurent v a)+  -> [TestTree]+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+  , testProperty "subst x p = p" $+    \p -> e Data.Poly.X p === p+  , testProperty "subst (monomial 0 c) p = monomial 0 c" $+    \c p -> e (Data.Poly.monomial 0 c) p === monomial 0 c+  ]+  where+    e :: Data.Poly.Poly v a -> Laurent v a -> Laurent v a+    e = subst++derivTests :: TestTree+derivTests = testGroup "deriv"+  [ testProperty "deriv c = 0" $+    \c -> deriv (monomial 0 c :: Laurent V.Vector Int) === 0+  , testProperty "deriv cX = c" $+    \c -> deriv (monomial 0 c * X :: Laurent V.Vector 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)+  , 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)+  , 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) ->+      deriv (subst p q) === deriv q * subst (Data.Poly.deriv p) q+  ]
test/Main.hs view
@@ -3,10 +3,16 @@ 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  main :: IO () main = defaultMain $ testGroup "All"     [ Dense.testSuite+    , DenseLaurent.testSuite     , Sparse.testSuite+    , SparseLaurent.testSuite+    , Orthogonal.testSuite     ]
+ test/Orthogonal.hs view
@@ -0,0 +1,155 @@+{-# LANGUAGE OverloadedLists #-}++module Orthogonal+  ( testSuite+  ) where++import Test.Tasty++import Data.List (foldl', tails)+import Data.Poly (VPoly, deriv, eval, integral)+import Data.Poly.Orthogonal+import Test.Tasty.QuickCheck++testSuite :: TestTree+testSuite = testGroup "Orthogonal"+  [ testGroup "differential equations"+    [ testProperty "jacobi"      prop_jacobi_de+    , testProperty "gegenbauer"  prop_gegenbauer_de+    , testProperty "legendre"    prop_legendre_de+    , testProperty "chebyshev1"  prop_chebyshev1_de+    , testProperty "chebyshev2"  prop_chebyshev2_de+    , testProperty "hermitePhys" prop_hermitePhys_de+    , testProperty "laguerre"    prop_laguerre_de+    , testProperty "laguerreGen" prop_laguerreGen_de+    ]+  , testGroup "normalization"+    [ testProperty "jacobi"     prop_jacobi_norm+    , testProperty "gegenbauer" prop_gegenbauer_norm+    , testProperty "legendre"   prop_legendre_norm+    , testProperty "chebyshev1" prop_chebyshev1_norm+    , testProperty "chebyshev2" prop_chebyshev2_norm+    ]+  , testGroup "orthogonality"+    [ testProperty "legendre"   prop_legendre_orth+    ]+  , testGroup "Hermite"+    [ testProperty "hermiteProb" prop_hermiteProb+    , testProperty "hermitePhys" prop_hermitePhys+    ]+  ]++prop_jacobi_de :: Rational -> Rational -> Property+prop_jacobi_de a b = foldl' (.&&.) (property True) $+  zipWith (((=== 0) .) . de) [0..limit] (jacobi a b)+  where+    de :: Rational -> VPoly Rational -> VPoly Rational+    de n y = [1, 0, -1] * deriv (deriv y)+           + [b - a, - (a + b + 2)] * deriv y+           + [n * (n + a + b + 1)] * y++prop_gegenbauer_de :: Rational -> Property+prop_gegenbauer_de g = foldl' (.&&.) (property True) $+  zipWith (((=== 0) .) . de) [0..limit] (gegenbauer g)+  where+    de :: Rational -> VPoly Rational -> VPoly Rational+    de n y = [1, 0, -1] * deriv (deriv y)+           + [0, - (2 * g + 1)] * deriv y+           + [n * (n + 2 * g)] * y++prop_legendre_de :: Property+prop_legendre_de = once $ foldl' (.&&.) (property True) $+  zipWith (((=== 0) .) . de) [0..limit] legendre+  where+    de :: Rational -> VPoly Rational -> VPoly Rational+    de n y = deriv ([1, 0, -1] * deriv y) + [n * (n + 1)] * y++prop_chebyshev1_de :: Property+prop_chebyshev1_de = once $ foldl' (.&&.) (property True) $+  zipWith (((=== 0) .) . de) [0..limit] chebyshev1+  where+    de :: Integer -> VPoly Integer -> VPoly Integer+    de n y = [1, 0, -1] * deriv (deriv y) + [0, -1] * deriv y + [n * n] * y++prop_chebyshev2_de :: Property+prop_chebyshev2_de = once $ foldl' (.&&.) (property True) $+  zipWith (((=== 0) .) . de) [0..limit] chebyshev2+  where+    de :: Integer -> VPoly Integer -> VPoly Integer+    de n y = [1, 0, -1] * deriv (deriv y) + [0, -3] * deriv y + [n * (n + 2)] * y++prop_hermitePhys_de :: Property+prop_hermitePhys_de = once $ foldl' (.&&.) (property True) $+  zipWith (((=== 0) .) . de) [0..limit] hermitePhys+  where+    de :: Integer -> VPoly Integer -> VPoly Integer+    de n y = deriv (deriv y) + [0, -2] * deriv y + [2 * n] * y++prop_laguerre_de :: Property+prop_laguerre_de = once $ foldl' (.&&.) (property True) $+  zipWith (((=== 0) .) . de) [0..limit] laguerre+  where+    de :: Rational -> VPoly Rational -> VPoly Rational+    de n y = [0, 1] * deriv (deriv y) + [1, -1] * deriv y + [n] * y++prop_laguerreGen_de :: Rational -> Property+prop_laguerreGen_de a  = foldl' (.&&.) (property True) $+  zipWith (((=== 0) .) . de) [0..limit] (laguerreGen a)+  where+    de :: Rational -> VPoly Rational -> VPoly Rational+    de n y = [0, 1] * deriv (deriv y) + [1 + a, -1] * deriv y + [n] * y++prop_jacobi_norm :: Rational -> Rational -> Property+prop_jacobi_norm a b = foldl' (.&&.) (property True) $+  zipWith (\n y -> norm n === eval y 1) [0..limit] (jacobi a b :: [VPoly Rational])+  where+    prod n x = product $ take n $ iterate (subtract 1) (fromIntegral n + x)+    norm n = prod n a / prod n 0++prop_gegenbauer_norm :: Rational -> Property+prop_gegenbauer_norm a = foldl' (.&&.) (property True) $+  zipWith (\n y -> norm n === eval y 1) [0..limit] (gegenbauer a :: [VPoly Rational])+  where+    prod n x = product $ take n $ iterate (subtract 1) (fromIntegral n + x)+    norm n = prod n (a - 1 / 2) / prod n 0++prop_legendre_norm :: Property+prop_legendre_norm = once $ foldl' (.&&.) (property True) $+  map ((=== 1) . flip eval 1) (take limit legendre :: [VPoly Rational])++prop_chebyshev1_norm :: Property+prop_chebyshev1_norm = once $ foldl' (.&&.) (property True) $+  map ((=== 1) . flip eval 1) (take limit chebyshev1 :: [VPoly Integer])++prop_chebyshev2_norm :: Property+prop_chebyshev2_norm = once $ foldl' (.&&.) (property True) $+  zipWith (\n y -> n + 1 === eval y 1) [0..limit] (chebyshev2 :: [VPoly Integer])++prop_legendre_orth :: Property+prop_legendre_orth = once $ foldl' (.&&.) (property True) $+  [ integral11 (x * y) === 0 | (x : xs) <- tails polys, y <- xs ]+  where+    polys :: [VPoly Rational]+    polys = take limit $ legendre++hermiteProbRef :: [VPoly Integer]+hermiteProbRef = iterate (\he -> [0, 1] * he - deriv he) 1++hermitePhysRef :: [VPoly Integer]+hermitePhysRef = iterate (\h -> [0, 2] * h - deriv h) 1++prop_hermiteProb :: Property+prop_hermiteProb = once $ foldl' (.&&.) (property True) $+  take limit $ zipWith (===) hermiteProb hermiteProbRef++prop_hermitePhys :: Property+prop_hermitePhys = once $ foldl' (.&&.) (property True) $+  take limit $ zipWith (===) hermitePhys hermitePhysRef++integral11 :: VPoly Rational -> Rational+integral11 x = eval y 1 - eval y (-1)+  where+    y = integral x++limit :: Num a => a+limit = 10
test/Quaternion.hs view
@@ -29,7 +29,7 @@ import qualified Data.Vector.Generic.Mutable as M import Data.Vector.Unboxed (Unbox) -data Quaternion a = Quaternion a a a a+data Quaternion a = Quaternion !a !a !a !a   deriving (Eq, Ord, Show, Generic)  instance Ring a => Semiring (Quaternion a) where
test/Sparse.hs view
@@ -1,4 +1,4 @@-{-# LANGUAGE CPP                        #-}+{-# LANGUAGE DataKinds                  #-} {-# LANGUAGE FlexibleContexts           #-} {-# LANGUAGE FlexibleInstances          #-} {-# LANGUAGE GeneralizedNewtypeDeriving #-}@@ -9,16 +9,15 @@  module Sparse   ( testSuite+  , ShortPoly(..)   ) where  import Prelude hiding (gcd, quotRem, rem)-#if MIN_VERSION_semirings(0,4,2) import Data.Euclidean (Euclidean(..), GcdDomain(..))-#endif import Data.Function import Data.Int-import Data.List-import Data.Maybe+import Data.List (groupBy, sortOn)+import Data.Mod import Data.Poly.Sparse import qualified Data.Poly.Sparse.Semiring as S import Data.Proxy@@ -28,24 +27,16 @@ import qualified Data.Vector.Unboxed as U import Test.Tasty import Test.Tasty.QuickCheck hiding (scale, numTests)-import Test.QuickCheck.Classes  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-#if MIN_VERSION_semirings(0,4,2)-    , GcdDomain-    , Euclidean-#endif-    )+  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@@ -58,94 +49,68 @@     , lawsTests     , evalTests     , derivTests-#if MIN_VERSION_semirings(0,4,2)-    , gcdExtTests-#endif     ]  lawsTests :: TestTree lawsTests = testGroup "Laws"-  [ semiringTests-  , ringTests-  , numTests-  , euclideanTests-  , isListTests-  , showTests+  $ semiringTests ++ ringTests ++ numTests ++ euclideanTests ++ gcdDomainTests ++ isListTests ++ showTests++semiringTests :: [TestTree]+semiringTests =+  [ mySemiringLaws (Proxy :: Proxy (Poly U.Vector ()))+  , mySemiringLaws (Proxy :: Proxy (Poly U.Vector Int8))+  , mySemiringLaws (Proxy :: Proxy (Poly V.Vector Integer))+  , tenTimesLess+  $ mySemiringLaws (Proxy :: Proxy (Poly U.Vector (Quaternion Int)))   ] -semiringTests :: TestTree-semiringTests-  = testGroup "Semiring"-  $ map (uncurry testProperty)-  $ concatMap lawsProperties-  [ semiringLaws (Proxy :: Proxy (Poly U.Vector ()))-  , semiringLaws (Proxy :: Proxy (Poly U.Vector Int8))-  , semiringLaws (Proxy :: Proxy (Poly V.Vector Integer))-  , semiringLaws (Proxy :: Proxy (Poly U.Vector (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)))   ] -ringTests :: TestTree-ringTests-  = testGroup "Ring"-  $ map (uncurry testProperty)-  $ concatMap lawsProperties-  [-#if MIN_VERSION_quickcheck_classes(0,6,1)-    ringLaws (Proxy :: Proxy (Poly U.Vector ()))-  , ringLaws (Proxy :: Proxy (Poly U.Vector Int8))-  , ringLaws (Proxy :: Proxy (Poly V.Vector Integer))-  , ringLaws (Proxy :: Proxy (Poly U.Vector (Quaternion Int)))-#endif+numTests :: [TestTree]+numTests =+  [ myNumLaws (Proxy :: Proxy (Poly U.Vector Int8))+  , myNumLaws (Proxy :: Proxy (Poly V.Vector Integer))+  , tenTimesLess+  $ myNumLaws (Proxy :: Proxy (Poly U.Vector (Quaternion Int)))   ] -numTests :: TestTree-numTests-  = testGroup "Num"-  $ map (uncurry testProperty)-  $ concatMap lawsProperties-  [-#if MIN_VERSION_quickcheck_classes(0,6,3)-    numLaws (Proxy :: Proxy (Poly U.Vector Int8))-  , numLaws (Proxy :: Proxy (Poly V.Vector Integer))-  , numLaws (Proxy :: Proxy (Poly U.Vector (Quaternion Int)))-#endif+gcdDomainTests :: [TestTree]+gcdDomainTests =+  [ myGcdDomainLaws (Proxy :: Proxy (ShortPoly (Poly V.Vector Integer)))+  , tenTimesLess+  $ myGcdDomainLaws (Proxy :: Proxy (ShortPoly (Poly V.Vector (Mod 3))))+  , tenTimesLess+  $ myGcdDomainLaws (Proxy :: Proxy (ShortPoly (Poly V.Vector Rational)))   ] -euclideanTests :: TestTree-euclideanTests-  = testGroup "Euclidean"-  $ map (uncurry testProperty)-  $ concatMap lawsProperties-  [-#if MIN_VERSION_semirings(0,4,2) && MIN_VERSION_quickcheck_classes(0,6,3)-    gcdDomainLaws (Proxy :: Proxy (ShortPoly (Poly V.Vector Integer)))-  , euclideanLaws (Proxy :: Proxy (ShortPoly (Poly V.Vector Rational)))-#endif+euclideanTests :: [TestTree]+euclideanTests =+  [ myEuclideanLaws (Proxy :: Proxy (ShortPoly (Poly V.Vector (Mod 3))))+  , myEuclideanLaws (Proxy :: Proxy (ShortPoly (Poly V.Vector Rational)))   ] -isListTests :: TestTree-isListTests-  = testGroup "IsList"-  $ map (uncurry testProperty)-  $ concatMap lawsProperties-  [ isListLaws (Proxy :: Proxy (Poly U.Vector ()))-  , isListLaws (Proxy :: Proxy (Poly U.Vector Int8))-  , isListLaws (Proxy :: Proxy (Poly V.Vector Integer))-  , isListLaws (Proxy :: Proxy (Poly U.Vector (Quaternion Int)))+isListTests :: [TestTree]+isListTests =+  [ myIsListLaws (Proxy :: Proxy (Poly U.Vector ()))+  , myIsListLaws (Proxy :: Proxy (Poly U.Vector Int8))+  , myIsListLaws (Proxy :: Proxy (Poly V.Vector Integer))+  , tenTimesLess+  $ myIsListLaws (Proxy :: Proxy (Poly U.Vector (Quaternion Int)))   ] -showTests :: TestTree-showTests-  = testGroup "Show"-  $ map (uncurry testProperty)-  $ concatMap lawsProperties-  [-#if MIN_VERSION_quickcheck_classes(0,6,0)-    showLaws (Proxy :: Proxy (Poly U.Vector ()))-  , showLaws (Proxy :: Proxy (Poly U.Vector Int8))-  , showLaws (Proxy :: Proxy (Poly V.Vector Integer))-  , showLaws (Proxy :: Proxy (Poly U.Vector (Quaternion Int)))-#endif+showTests :: [TestTree]+showTests =+  [ myShowLaws (Proxy :: Proxy (Poly U.Vector ()))+  , myShowLaws (Proxy :: Proxy (Poly U.Vector Int8))+  , myShowLaws (Proxy :: Proxy (Poly V.Vector Integer))+  , tenTimesLess+  $ myShowLaws (Proxy :: Proxy (Poly U.Vector (Quaternion Int)))   ]  arithmeticTests :: TestTree@@ -156,7 +121,8 @@   , testProperty "subtraction matches reference" $     \(xs :: [(Word, Int)]) ys -> toPoly (V.fromList (subRef xs ys)) ===       toPoly (V.fromList xs) - toPoly (V.fromList ys)-  , testProperty "multiplication matches reference" $+  , tenTimesLess $+    testProperty "multiplication matches reference" $     \(xs :: [(Word, Int)]) ys -> toPoly (V.fromList (mulRef xs ys)) ===       toPoly (V.fromList xs) * toPoly (V.fromList ys)   ]@@ -204,7 +170,8 @@     \p c -> c /= 0 ==> leading (monomial p c :: UPoly a) === Just (p, c)   , testProperty "monomial matches reference" $     \p (c :: a) -> monomial p c === toPoly (V.fromList (monomialRef p c))-  , testProperty "scale matches multiplication by monomial" $+  , tenTimesLess $+    testProperty "scale matches multiplication by monomial" $     \p c (xs :: UPoly a) -> scale p c xs === monomial p c * xs   ] @@ -274,12 +241,10 @@ derivTests = testGroup "deriv"   [ testProperty "deriv = S.deriv" $     \(p :: Poly V.Vector Integer) -> deriv p === S.deriv p-#if MIN_VERSION_semirings(0,5,0)   , testProperty "integral = S.integral" $     \(p :: Poly V.Vector Rational) -> integral p === S.integral p   , testProperty "deriv . integral = id" $     \(p :: Poly V.Vector Rational) -> deriv (integral p) === p-#endif   , testProperty "deriv c = 0" $     \c -> deriv (monomial 0 c :: Poly V.Vector Int) === 0   , testProperty "deriv cX = c" $@@ -292,22 +257,3 @@   --   \(p :: Poly V.Vector Int) (q :: Poly U.Vector Int) ->   --     deriv (subst p q) === deriv q * subst (deriv p) q   ]--#if MIN_VERSION_semirings(0,4,2)-gcdExtTests :: TestTree-gcdExtTests = localOption (QuickCheckMaxSize 12) $ testGroup "gcdExt"-  [ testProperty "gcdExt == S.gcdExt" $-    \(a :: Poly V.Vector Rational) b ->-      gcdExt a b === S.gcdExt a b-  , testProperty "g == as (mod b) for gcdExt" $-    \(a :: Poly V.Vector Rational) b -> b /= 0 ==>-      uncurry ((. flip rem b) . (===) . flip rem b) ((* a) <$> gcdExt a b)-  , testProperty "fst . gcdExt == gcd (mod units)" $-    \(a :: Poly V.Vector Rational) b ->-      fst (gcdExt a b) `sameUpToUnits` gcd a b-  ]--sameUpToUnits :: (Eq a, GcdDomain a) => a -> a -> Bool-sameUpToUnits x y = x == y ||-  isJust (x `divide` y) && isJust (y `divide` x)-#endif
+ test/SparseLaurent.hs view
@@ -0,0 +1,177 @@+{-# 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 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 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+    ]++lawsTests :: TestTree+lawsTests = testGroup "Laws"+  $ semiringTests ++ ringTests ++ numTests ++ gcdDomainTests ++ isListTests ++ showTests++semiringTests :: [TestTree]+semiringTests =+  [ mySemiringLaws (Proxy :: Proxy (Laurent U.Vector ()))+  , mySemiringLaws (Proxy :: Proxy (Laurent U.Vector Int8))+  , mySemiringLaws (Proxy :: Proxy (Laurent V.Vector Integer))+  , tenTimesLess+  $ mySemiringLaws (Proxy :: Proxy (Laurent U.Vector (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)))+  ]++numTests :: [TestTree]+numTests =+  [ myNumLaws (Proxy :: Proxy (Laurent U.Vector Int8))+  , myNumLaws (Proxy :: Proxy (Laurent V.Vector Integer))+  , tenTimesLess+  $ myNumLaws (Proxy :: Proxy (Laurent U.Vector (Quaternion Int)))+  ]++gcdDomainTests :: [TestTree]+gcdDomainTests =+  [ myGcdDomainLaws (Proxy :: Proxy (ShortLaurent (Laurent V.Vector Integer)))+  , tenTimesLess+  $ myGcdDomainLaws (Proxy :: Proxy (ShortLaurent (Laurent V.Vector Rational)))+  ]++isListTests :: [TestTree]+isListTests =+  [ myIsListLaws (Proxy :: Proxy (Laurent U.Vector ()))+  , myIsListLaws (Proxy :: Proxy (Laurent U.Vector Int8))+  , myIsListLaws (Proxy :: Proxy (Laurent V.Vector Integer))+  , tenTimesLess+  $ myIsListLaws (Proxy :: Proxy (Laurent U.Vector (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))+  , tenTimesLess+  $ myShowLaws (Proxy :: Proxy (Laurent U.Vector (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 "leading p 0 == Nothing" $+    \p -> leading (monomial p 0 :: ULaurent a) === Nothing+  , testProperty "leading . monomial = id" $+    \p c -> c /= 0 ==> leading (monomial p c :: ULaurent a) === Just (p, c)+  , tenTimesLess $+    testProperty "scale matches multiplication by monomial" $+    \p c (xs :: ULaurent a) -> scale p c xs === monomial p c * xs+  ]++evalTests :: TestTree+evalTests = testGroup "eval" $ concat+  [ evalTestGroup  (Proxy :: Proxy (Laurent V.Vector Rational))+  , substTestGroup (Proxy :: Proxy (Laurent U.Vector 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))+  => 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+  , 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+  , 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+  ]+  where+    e :: Laurent v a -> a -> a+    e = eval++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))+  => Proxy (Laurent v a)+  -> [TestTree]+substTestGroup _ =+  [ testProperty "subst x p = p" $+    \p -> e Data.Poly.Sparse.X p === p+  , testProperty "subst (monomial 0 c) p = monomial 0 c" $+    \c p -> e (Data.Poly.Sparse.monomial 0 c) p === monomial 0 c+  ]+  where+    e :: Data.Poly.Sparse.Poly v a -> Laurent v a -> Laurent v a+    e = subst++derivTests :: TestTree+derivTests = testGroup "deriv"+  [ testProperty "deriv c = 0" $+    \c -> deriv (monomial 0 c :: Laurent V.Vector Int) === 0+  , testProperty "deriv cX = c" $+    \c -> deriv (monomial 0 c * X :: Laurent V.Vector 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)+  , 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+  ]
+ test/TestUtils.hs view
@@ -0,0 +1,113 @@+{-# LANGUAGE CPP              #-}+{-# LANGUAGE FlexibleContexts #-}++{-# OPTIONS_GHC -fno-warn-orphans #-}++module TestUtils+  ( tenTimesLess+  , mySemiringLaws+  , myRingLaws+  , myNumLaws+  , myGcdDomainLaws+  , myEuclideanLaws+  , myIsListLaws+  , myShowLaws+  ) where++import Data.Euclidean+import Data.Mod+import Data.Proxy+import Data.Semiring (Semiring, Ring)+import GHC.Exts+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++instance KnownNat m => Arbitrary (Mod m) where+  arbitrary = oneof [arbitraryBoundedEnum, fromInteger <$> arbitrary]+  shrink = map fromInteger . shrink . toInteger . unMod++tenTimesLess :: TestTree -> TestTree+tenTimesLess = adjustOption $+  \(QuickCheckTests n) -> QuickCheckTests (max 100 (n `div` 10))++mySemiringLaws :: (Eq a, Semiring a, Arbitrary a, Show a) => Proxy a -> TestTree+mySemiringLaws proxy = testGroup tpclss $ map tune props+  where+    Laws tpclss props = semiringLaws proxy++    tune pair = case fst pair of+      "Multiplicative Associativity" ->+        tenTimesLess test+      "Multiplication Left Distributes Over Addition" ->+        tenTimesLess test+      "Multiplication Right Distributes Over Addition" ->+        tenTimesLess test+      _ -> test+      where+        test = uncurry testProperty pair++myRingLaws :: (Eq a, Ring a, Arbitrary a, Show a) => Proxy a -> TestTree+myRingLaws proxy = testGroup tpclss $ map (uncurry testProperty) props+  where+    Laws tpclss props = ringLaws proxy++myNumLaws :: (Eq a, Num a, Arbitrary a, Show a) => Proxy a -> TestTree+myNumLaws proxy = testGroup tpclss $ map tune props+  where+    Laws tpclss props = numLaws proxy++    tune pair = case fst pair of+      "Multiplicative Associativity" ->+        tenTimesLess test+      "Multiplication Left Distributes Over Addition" ->+        tenTimesLess test+      "Multiplication Right Distributes Over Addition" ->+        tenTimesLess test+      "Subtraction" ->+        tenTimesLess test+      _ -> test+      where+        test = uncurry testProperty pair++myGcdDomainLaws :: (Eq a, GcdDomain a, Arbitrary a, Show a) => Proxy a -> TestTree+myGcdDomainLaws proxy = testGroup tpclss $ map tune props+  where+    Laws tpclss props = gcdDomainLaws proxy++    tune pair = case fst pair of+      "gcd1"    -> tenTimesLess $ tenTimesLess test+      "gcd2"    -> tenTimesLess $ tenTimesLess test+      "lcm1"    -> tenTimesLess $ tenTimesLess $ tenTimesLess test+      "lcm2"    -> tenTimesLess test+      "coprime" -> tenTimesLess $ tenTimesLess test+      _ -> test+      where+        test = uncurry testProperty pair++myEuclideanLaws :: (Eq a, Euclidean a, Arbitrary a, Show a) => Proxy a -> TestTree+myEuclideanLaws proxy = testGroup tpclss $ map (uncurry testProperty) props+  where+    Laws tpclss props = euclideanLaws proxy++myIsListLaws :: (Eq a, IsList a, Arbitrary a, Show a, Show (Item a), Arbitrary (Item a)) => Proxy a -> TestTree+myIsListLaws proxy = testGroup tpclss $ map (uncurry testProperty) props+  where+    Laws tpclss props = isListLaws proxy++myShowLaws :: (Eq a, Arbitrary a, Show a) => Proxy a -> TestTree+myShowLaws proxy = testGroup tpclss $ map tune props+  where+    Laws tpclss props = showLaws proxy++    tune pair = case fst pair of+      "Equivariance: showList" -> tenTimesLess $ tenTimesLess test+      _ -> test+      where+        test = uncurry testProperty pair