poly 0.3.1.0 → 0.3.2.0

raw patch · 21 files changed

+908/−471 lines, 21 filesdep +deepseqdep ~semiringsdep ~vectorPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

Dependencies added: deepseq

Dependency ranges changed: semirings, vector

API changes (from Hackage documentation)

Files

• README.md +2/−2
• bench/DenseBench.hs +78/−15
• bench/SparseBench.hs +0/−1
• changelog.md +7/−0
• poly.cabal +6/−4
• src/Data/Poly.hs +10/−8
• src/Data/Poly/Internal/Dense.hs +66/−46
• src/Data/Poly/Internal/Dense/Field.hs +189/−0
• src/Data/Poly/Internal/Dense/Fractional.hs +0/−138
• src/Data/Poly/Internal/Dense/GcdDomain.hs +19/−23
• src/Data/Poly/Internal/PolyOverField.hs +75/−0
• src/Data/Poly/Internal/PolyOverFractional.hs +0/−55
• src/Data/Poly/Internal/Sparse.hs +74/−55
• src/Data/Poly/Internal/Sparse/Field.hs +118/−0
• src/Data/Poly/Internal/Sparse/Fractional.hs +0/−78
• src/Data/Poly/Internal/Sparse/GcdDomain.hs +9/−13
• src/Data/Poly/Semiring.hs +26/−6
• src/Data/Poly/Sparse.hs +9/−2
• src/Data/Poly/Sparse/Semiring.hs +26/−3
• test/Dense.hs +97/−18
• test/Sparse.hs +97/−4

@@ -76,10 +76,10 @@ cover main functionality of `Integral`, providing division with remainder and `gcd` / `lcm`:  ```haskell-> Data.Euclidean.gcd (X^2 + 7 * X + 6) (X^2 - 5 * X - 6) :: Data.Poly.UPoly Int+> Data.Euclidean.gcd (X^2 + 7 * X + 6) (X^2 - 5 * X - 6) :: UPoly Int 1 * X + 1 -> Data.Euclidean.quotRem (X^3 + 2) (X^2 - 1 :: Data.Poly.UPoly Double)+> Data.Euclidean.quotRem (X^3 + 2) (X^2 - 1 :: UPoly Double) (1.0 * X + 0.0,1.0 * X + 2.0) ``` 

@@ -1,16 +1,24 @@-{-# LANGUAGE CPP        #-}-{-# LANGUAGE RankNTypes #-}+{-# LANGUAGE CPP                        #-}+{-# LANGUAGE GeneralizedNewtypeDeriving #-}+{-# LANGUAGE RankNTypes                 #-}+{-# LANGUAGE TypeApplications           #-}  module DenseBench   ( benchSuite   ) where  import Prelude hiding (quotRem, gcd)+import Control.DeepSeq import Gauge.Main import Data.Poly import qualified Data.Vector.Unboxed as U #if MIN_VERSION_semirings(0,4,2)+import Control.Exception+import Data.Bits+import Data.Coerce import Data.Euclidean+import Data.Semiring (Semiring(..), Ring, isZero)+import qualified Data.Semiring as S (negate) import qualified Data.Vector as V #endif @@ -22,9 +30,12 @@   , map benchDeriv    [100, 1000, 10000]   , map benchIntegral [100, 1000, 10000] #if MIN_VERSION_semirings(0,4,2)-  , map benchQuotRem  [10, 100]-  , map benchGcdFrac  [10, 100]-  , map benchGcd      [10, 100]+  , 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   ] @@ -51,9 +62,18 @@ benchGcd :: Int -> Benchmark benchGcd k = bench ("gcd/" ++ show k) $ nf doGcd k -benchGcdFrac :: Int -> Benchmark-benchGcdFrac k = bench ("gcdFrac/" ++ show k) $ nf doGcdFrac 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 (doGcdExt @Mod2) k++benchGcdFracM :: Int -> Benchmark+benchGcdFracM k = bench ("gcdFrac/Mod2/" ++ show k) $ nf (doGcdFrac @Mod2) k+ #endif  doBinOp :: (forall a. Num a => a -> a -> a) -> Int -> Int@@ -83,25 +103,68 @@  #if MIN_VERSION_semirings(0,4,2) +gen1 :: Num a => Int -> a+gen1 k = fromIntegral (truncate (pi * fromIntegral k :: Double) `mod` (k + 1))++gen2 :: Num a => Int -> a+gen2 k = fromIntegral (truncate (exp 1.0 * fromIntegral k :: Double) `mod` (k + 1))+ doQuotRem :: Int -> Double doQuotRem n = U.sum (unPoly qs) + U.sum (unPoly rs)   where-    xs = toPoly $ U.generate (2 * n) ((+ 1.0) . (* 2.0) . fromIntegral)-    ys = toPoly $ U.generate n       ((+ 2.0) . (* 3.0) . fromIntegral)+    xs = toPoly $ U.generate (2 * n) gen1+    ys = toPoly $ U.generate n       gen2     (qs, rs) = xs `quotRem` ys  doGcd :: Int -> Integer doGcd n = V.sum gs   where-    xs = toPoly $ V.generate n ((+ 1) . (* 2) . fromIntegral)-    ys = toPoly $ V.generate n ((+ 2) . (* 3) . fromIntegral)+    xs = toPoly $ V.generate n gen1+    ys = toPoly $ V.generate n gen2     gs = unPoly $ xs `gcd` ys -doGcdFrac :: Int -> Rational+doGcdExt :: (Eq a, Num a, Field a) => Int -> a+doGcdExt n = V.sum gs+  where+    xs = toPoly $ V.generate n gen1+    ys = toPoly $ V.generate n gen2+    gs = unPoly $ fst $ xs `gcdExt` ys++doGcdFrac :: (Eq a, Num a, Field a) => Int -> a doGcdFrac n = V.sum gs   where-    xs = PolyOverFractional $ toPoly $ V.generate n ((+ 1) . (* 2) . fromIntegral)-    ys = PolyOverFractional $ toPoly $ V.generate n ((+ 2) . (* 3) . fromIntegral)-    gs = unPoly $ unPolyOverFractional $ xs `gcd` ys+    xs = PolyOverField $ toPoly $ V.generate n gen1+    ys = PolyOverField $ toPoly $ V.generate n gen2+    gs = unPoly $ unPolyOverField $ xs `gcd` ys++-- | Inspired by 'semirings'.+newtype Mod2 = Mod2 { _getMod2 :: Bool }+  deriving (Eq, NFData)++instance Num Mod2 where+  (+) = coerce (xor @Bool)+  (*) = coerce (&&)+  negate = id+  abs    = id+  signum = id+  fromInteger = Mod2 . odd++instance Semiring Mod2 where+  plus  = coerce (xor @Bool)+  times = coerce (&&)+  fromNatural = Mod2 . odd++instance Ring Mod2 where+  negate = id++instance GcdDomain Mod2 where++instance Euclidean Mod2 where+  degree = const 0+  quotRem x y+    | isZero y  = throw DivideByZero+    | otherwise = (x, zero)++instance Field Mod2  #endif

@@ -67,4 +67,3 @@ doIntegral xs = U.foldl' (\acc (_, x) -> acc + x) 0 zs   where     zs = unPoly $ integral xs-

@@ -1,3 +1,10 @@+# 0.3.2.0++* Add `NFData` instance.+* Implement extended GCD.+* Rename `PolyOverFractional` to `PolyOverField`.+* Add `integral` with `Semiring`-based interface.+ # 0.3.1.0  * Implement Karatsuba multiplication.

@@ -1,5 +1,5 @@ name: poly-version: 0.3.1.0+version: 0.3.2.0 synopsis: Polynomials description:   Polynomials backed by `Vector`.@@ -30,14 +30,15 @@     Data.Poly.Sparse.Semiring   other-modules:     Data.Poly.Internal.Dense-    Data.Poly.Internal.Dense.Fractional+    Data.Poly.Internal.Dense.Field     Data.Poly.Internal.Dense.GcdDomain-    Data.Poly.Internal.PolyOverFractional+    Data.Poly.Internal.PolyOverField     Data.Poly.Internal.Sparse-    Data.Poly.Internal.Sparse.Fractional+    Data.Poly.Internal.Sparse.Field     Data.Poly.Internal.Sparse.GcdDomain   build-depends:     base >= 4.9 && < 5,+    deepseq >= 1.1 && < 1.5,     primitive >= 0.6,     semirings >= 0.2,     vector >= 0.12.0.2,@@ -68,6 +69,7 @@ benchmark poly-gauge   build-depends:     base >=4.9 && <5,+    deepseq >= 1.1 && < 1.5,     gauge >= 0.1,     poly,     semirings >= 0.2,

@@ -7,10 +7,8 @@ -- Dense polynomials and a 'Num'-based interface. -- -{-# LANGUAGE CPP                        #-}-{-# LANGUAGE FlexibleInstances          #-}-{-# LANGUAGE GeneralizedNewtypeDeriving #-}-{-# LANGUAGE PatternSynonyms            #-}+{-# LANGUAGE CPP             #-}+{-# LANGUAGE PatternSynonyms #-}  module Data.Poly   ( Poly@@ -27,14 +25,18 @@   , deriv   , integral #if MIN_VERSION_semirings(0,4,2)-  -- * Fractional coefficients-  , PolyOverFractional(..)+  -- * 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.Fractional ()+import Data.Poly.Internal.Dense.Field (gcdExt) import Data.Poly.Internal.Dense.GcdDomain ()-import Data.Poly.Internal.PolyOverFractional+import Data.Poly.Internal.PolyOverField #endif

@@ -36,14 +36,19 @@   , pattern X'   , eval'   , deriv'+#if MIN_VERSION_semirings(0,5,0)+  , integral'+#endif   ) where -import Prelude hiding (quotRem, rem, gcd, lcm, (^))+import Prelude hiding (quotRem, quot, rem, gcd, lcm, (^))+import Control.DeepSeq (NFData) import Control.Monad import Control.Monad.Primitive import Control.Monad.ST+import Data.Bits import Data.List (foldl', intersperse)-import Data.Semiring (Semiring(..))+import Data.Semiring (Semiring(..), Ring()) import qualified Data.Semiring as Semiring import qualified Data.Vector as V import qualified Data.Vector.Generic as G@@ -54,6 +59,9 @@ 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.).@@ -76,7 +84,7 @@   -- ^ Convert 'Poly' to a vector of coefficients   -- (first element corresponds to a constant term).   }-  deriving (Eq, Ord)+  deriving (Eq, NFData, Ord)  instance (Eq a, Semiring a, G.Vector v a) => IsList (Poly v a) where   type Item (Poly v a) = a@@ -139,7 +147,7 @@   abs = id   signum = const 1   fromInteger n = case fromInteger n of-    0 -> Poly $ G.empty+    0 -> Poly G.empty     m -> Poly $ G.singleton m   Poly xs * Poly ys = toPoly $ karatsuba xs ys   {-# INLINE (+) #-}@@ -168,7 +176,7 @@   {-# INLINE fromNatural #-} #endif -instance (Eq a, Semiring.Ring a, G.Vector v a) => Semiring.Ring (Poly v a) where+instance (Eq a, Ring a, G.Vector v a) => Ring (Poly v a) where   negate (Poly xs) = Poly $ G.map Semiring.negate xs  dropWhileEnd@@ -176,7 +184,7 @@   => (a -> Bool)   -> v a   -> v a-dropWhileEnd p xs = G.basicUnsafeSlice 0 (go (G.basicLength xs)) xs+dropWhileEnd p xs = G.unsafeSlice 0 (go (G.length xs)) xs   where     go 0 = 0     go n = if p (G.unsafeIndex xs (n - 1)) then go (n - 1) else n@@ -187,12 +195,12 @@   => (a -> Bool)   -> G.Mutable v (PrimState m) a   -> m (G.Mutable v (PrimState m) a)-dropWhileEndM p xs = go (MG.basicLength xs)+dropWhileEndM p xs = go (MG.length xs)   where-    go 0 = pure $ MG.basicUnsafeSlice 0 0 xs+    go 0 = pure $ MG.unsafeSlice 0 0 xs     go n = do       x <- MG.unsafeRead xs (n - 1)-      if p x then go (n - 1) else pure (MG.basicUnsafeSlice 0 n xs)+      if p x then go (n - 1) else pure (MG.unsafeSlice 0 n xs) {-# INLINE dropWhileEndM #-}  plusPoly@@ -202,17 +210,17 @@   -> v a   -> v a plusPoly add xs ys = runST $ do-  let lenXs = G.basicLength xs-      lenYs = G.basicLength ys+  let lenXs = G.length xs+      lenYs = G.length ys       lenMn = lenXs `min` lenYs       lenMx = lenXs `max` lenYs -  zs <- MG.basicUnsafeNew lenMx+  zs <- MG.unsafeNew lenMx   forM_ [0 .. lenMn - 1] $ \i ->     MG.unsafeWrite zs i (add (G.unsafeIndex xs i) (G.unsafeIndex ys i))   G.unsafeCopy-    (MG.basicUnsafeSlice lenMn (lenMx - lenMn) zs)-    (G.basicUnsafeSlice  lenMn (lenMx - lenMn) (if lenXs <= lenYs then ys else xs))+    (MG.unsafeSlice lenMn (lenMx - lenMn) zs)+    (G.unsafeSlice  lenMn (lenMx - lenMn) (if lenXs <= lenYs then ys else xs))    G.unsafeFreeze zs {-# INLINE plusPoly #-}@@ -225,12 +233,12 @@   -> v a   -> v a minusPoly neg sub xs ys = runST $ do-  let lenXs = G.basicLength xs-      lenYs = G.basicLength ys+  let lenXs = G.length xs+      lenYs = G.length ys       lenMn = lenXs `min` lenYs       lenMx = lenXs `max` lenYs -  zs <- MG.basicUnsafeNew lenMx+  zs <- MG.unsafeNew lenMx   forM_ [0 .. lenMn - 1] $ \i ->     MG.unsafeWrite zs i (sub (G.unsafeIndex xs i) (G.unsafeIndex ys i)) @@ -238,8 +246,8 @@     then forM_ [lenXs .. lenYs - 1] $ \i ->       MG.unsafeWrite zs i (neg (G.unsafeIndex ys i))     else G.unsafeCopy-      (MG.basicUnsafeSlice lenYs (lenXs - lenYs) zs)-      (G.basicUnsafeSlice  lenYs (lenXs - lenYs) xs)+      (MG.unsafeSlice lenYs (lenXs - lenYs) zs)+      (G.unsafeSlice  lenYs (lenXs - lenYs) xs)    G.unsafeFreeze zs {-# INLINE minusPoly #-}@@ -256,31 +264,31 @@   | lenXs <= karatsubaThreshold || lenYs <= karatsubaThreshold   = convolution 0 (+) (*) xs ys   | otherwise = runST $ do-    zs <- MG.basicUnsafeNew lenZs+    zs <- MG.unsafeNew lenZs     forM_ [0 .. lenZs - 1] $ \k -> do-      let z0 = if k < G.basicLength zs0+      let z0 = if k < G.length zs0                then G.unsafeIndex zs0 k                else 0-          z11 = if k - m >= 0 && k - m < G.basicLength zs11+          z11 = if k - m >= 0 && k - m < G.length zs11                then G.unsafeIndex zs11 (k - m)                else 0-          z10 = if k - m >= 0 && k - m < G.basicLength zs0+          z10 = if k - m >= 0 && k - m < G.length zs0                then G.unsafeIndex zs0 (k - m)                else 0-          z12 = if k - m >= 0 && k - m < G.basicLength zs2+          z12 = if k - m >= 0 && k - m < G.length zs2                then G.unsafeIndex zs2 (k - m)                else 0-          z2 = if k - 2 * m >= 0 && k - 2 * m < G.basicLength zs2+          z2 = if k - 2 * m >= 0 && k - 2 * m < G.length zs2                then G.unsafeIndex zs2 (k - 2 * m)                else 0       MG.unsafeWrite zs k (z0 + (z11 - z10 - z12) + z2)     G.unsafeFreeze zs   where-    lenXs = G.basicLength xs-    lenYs = G.basicLength ys+    lenXs = G.length xs+    lenYs = G.length ys     lenZs = lenXs + lenYs - 1 -    m    = ((lenXs `min` lenYs) + 1) `quot` 2+    m    = ((lenXs `min` lenYs) + 1) `shiftR` 1      xs0  = G.slice 0 m xs     xs1  = G.slice m (lenXs - m) xs@@ -303,18 +311,15 @@   -> v a   -> v a convolution zer add mul xs ys-  | G.null xs || G.null ys = G.empty-  | otherwise = runST $ do-    let lenXs = G.basicLength xs-        lenYs = G.basicLength ys-        lenZs = lenXs + lenYs - 1-    zs <- MG.basicUnsafeNew lenZs-    forM_ [0 .. lenZs - 1] $ \k -> do-      let is = [max (k - lenYs + 1) 0 .. min k (lenXs - 1)]-          acc = foldl' add zer $ flip map is $ \i ->-            mul (G.unsafeIndex xs i) (G.unsafeIndex ys (k - i))-      MG.unsafeWrite zs k acc-    G.unsafeFreeze zs+  | lenXs == 0 || lenYs == 0 = G.empty+  | otherwise = G.generate lenZs $ \k -> foldl'+    (\acc i -> acc `add` mul (G.unsafeIndex xs i) (G.unsafeIndex ys (k - i)))+    zer+    [max (k - lenYs + 1) 0 .. min k (lenXs - 1)]+  where+    lenXs = G.length xs+    lenYs = G.length ys+    lenZs = lenXs + lenYs - 1 {-# INLINE convolution #-}  -- | Create a monomial from a power and a coefficient.@@ -338,8 +343,8 @@   -> v a   -> v a scaleInternal zer mul yp yc xs = runST $ do-  let lenXs = G.basicLength xs-  zs <- MG.basicUnsafeNew (fromIntegral yp + lenXs)+  let lenXs = G.length xs+  zs <- MG.unsafeNew (fromIntegral yp + lenXs)   forM_ [0 .. fromIntegral yp - 1] $ \k ->     MG.unsafeWrite zs k zer   forM_ [0 .. lenXs - 1] $ \k ->@@ -372,12 +377,12 @@ -- 1 * X^2 + 2 * X + 2 eval :: (Num a, G.Vector v a) => Poly v a -> a -> a eval (Poly cs) x = fst' $-  G.foldl' (\(acc :*: xn) cn -> (acc + cn * xn :*: x * xn)) (0 :*: 1) cs+  G.foldl' (\(acc :*: xn) cn -> acc + cn * xn :*: x * xn) (0 :*: 1) cs {-# INLINE eval #-}  eval' :: (Semiring a, G.Vector v a) => Poly v a -> a -> a eval' (Poly cs) x = fst' $-  G.foldl' (\(acc :*: xn) cn -> (acc `plus` cn `times` xn :*: x `times` xn)) (zero :*: one) cs+  G.foldl' (\(acc :*: xn) cn -> acc `plus` cn `times` xn :*: x `times` xn) (zero :*: one) cs {-# INLINE eval' #-}  -- | Take a derivative.@@ -416,14 +421,29 @@ integral (Poly xs)   | G.null xs = Poly G.empty   | otherwise = toPoly $ runST $ do-    zs <- MG.basicUnsafeNew (lenXs + 1)+    zs <- MG.unsafeNew (lenXs + 1)     MG.unsafeWrite zs 0 0     forM_ [0 .. lenXs - 1] $ \i ->       MG.unsafeWrite zs (i + 1) (G.unsafeIndex xs i * recip (fromIntegral i + 1))     G.unsafeFreeze zs     where-      lenXs = G.basicLength xs+      lenXs = G.length xs {-# INLINE 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+  | otherwise = toPoly' $ runST $ do+    zs <- MG.unsafeNew (lenXs + 1)+    MG.unsafeWrite zs zero zero+    forM_ [0 .. lenXs - 1] $ \i ->+      MG.unsafeWrite zs (i + 1) (G.unsafeIndex xs i `quot` Semiring.fromIntegral (i + 1))+    G.unsafeFreeze zs+    where+      lenXs = G.length xs+{-# INLINE integral' #-}+#endif  -- | Create an identity polynomial. pattern X :: (Eq a, Num a, G.Vector v a, Eq (v a)) => Poly v a

@@ -0,0 +1,189 @@+-- |+-- Module:      Data.Poly.Internal.Dense.Field+-- Copyright:   (c) 2019 Andrew Lelechenko+-- Licence:     BSD3+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+--+-- GcdDomain for Field underlying+--++{-# LANGUAGE ConstraintKinds            #-}+{-# LANGUAGE CPP                        #-}+{-# LANGUAGE FlexibleInstances          #-}+{-# LANGUAGE PatternSynonyms            #-}+{-# LANGUAGE ScopedTypeVariables        #-}+{-# LANGUAGE TypeFamilies               #-}++{-# 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 Control.Exception+import Control.Monad+import Control.Monad.Primitive+import Control.Monad.ST+import Data.Euclidean+#if !MIN_VERSION_semirings(0,5,0)+import Data.Semiring (Ring)+#endif+import Data.Semiring (times, minus, zero, one)+import qualified Data.Vector.Generic as G+import qualified Data.Vector.Generic.Mutable as MG++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)++  quotRem (Poly xs) (Poly ys) = (toPoly' qs, toPoly' rs)+    where+      (qs, rs) = quotientAndRemainder xs ys+  {-# INLINE quotRem #-}++  rem (Poly xs) (Poly ys) = toPoly' $ remainder xs ys+  {-# INLINE rem #-}++quotientAndRemainder+  :: (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)+  | 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+    forM_ [lenQs - 1, lenQs - 2 .. 0] $ \i -> do+      r <- MG.unsafeRead rs (lenYs - 1 + i)+      let q = r `quot` G.unsafeLast ys+      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)+    let rs' = MG.unsafeSlice 0 lenYs rs+    (,) <$> G.unsafeFreeze qs <*> G.unsafeFreeze rs'+{-# INLINE quotientAndRemainder #-}++remainder+  :: (Field a, G.Vector v a)+  => v a+  -> v a+  -> v a+remainder xs ys+  | G.null ys = throw DivideByZero+  | otherwise = runST $ do+    rs <- G.thaw xs+    ys' <- G.unsafeThaw ys+    remainderM rs ys'+    G.unsafeFreeze $ MG.unsafeSlice 0 (G.length xs `min` G.length ys) rs+{-# INLINE remainder #-}++remainderM+  :: (PrimMonad m, 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 ()+  | otherwise = do+    let lenXs = MG.length xs+        lenYs = MG.length ys+        lenQs = lenXs - lenYs + 1+    yLast <- MG.unsafeRead ys (lenYs - 1)+    forM_ [lenQs - 1, lenQs - 2 .. 0] $ \i -> do+      r <- MG.unsafeRead xs (lenYs - 1 + i)+      forM_ [0 .. lenYs - 1] $ \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)+{-# INLINE remainderM #-}++fieldGcd+  :: (Eq a, Field a, G.Vector v a)+  => Poly v a+  -> Poly v a+  -> Poly v a+fieldGcd (Poly xs) (Poly ys) = toPoly' $ runST $ do+  xs' <- G.thaw xs+  ys' <- G.thaw ys+  gcdM xs' ys'+{-# INLINE fieldGcd #-}++gcdM+  :: (PrimMonad m, Eq a, Field a, G.Vector v a)+  => G.Mutable v (PrimState m) a+  -> G.Mutable v (PrimState m) a+  -> m (v a)+gcdM xs ys = do+  ys' <- dropWhileEndM (== zero) ys+  if MG.null ys' then G.unsafeFreeze xs else do+    remainderM xs ys'+    gcdM ys' xs+{-# INLINE gcdM #-}++-- | 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'+{-# INLINE 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

@@ -1,138 +0,0 @@--- |--- Module:      Data.Poly.Internal.Dense.Fractional--- Copyright:   (c) 2019 Andrew Lelechenko--- Licence:     BSD3--- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>------ GcdDomain for Fractional underlying-----{-# LANGUAGE CPP                        #-}-{-# LANGUAGE FlexibleInstances          #-}-{-# LANGUAGE GeneralizedNewtypeDeriving #-}-{-# LANGUAGE PatternSynonyms            #-}-{-# LANGUAGE ScopedTypeVariables        #-}-{-# LANGUAGE TypeFamilies               #-}-{-# LANGUAGE ViewPatterns               #-}--{-# OPTIONS_GHC -fno-warn-orphans #-}--#if MIN_VERSION_semirings(0,4,2)--module Data.Poly.Internal.Dense.Fractional-  ( fractionalGcd-  ) where--import Prelude hiding (rem, gcd)-import Control.Exception-import Control.Monad-import Control.Monad.Primitive-import Control.Monad.ST-import Data.Euclidean-import qualified Data.Semiring as Semiring-import qualified Data.Vector.Generic as G-import qualified Data.Vector.Generic.Mutable as MG--import Data.Poly.Internal.Dense-import Data.Poly.Internal.Dense.GcdDomain ()--instance (Eq a, Eq (v a), Semiring.Ring a, GcdDomain a, Fractional a, G.Vector v a) => Euclidean (Poly v a) where-  degree (Poly xs) = fromIntegral (G.basicLength xs)--  quotRem (Poly xs) (Poly ys) = (toPoly qs, toPoly rs)-    where-      (qs, rs) = quotientAndRemainder xs ys-  {-# INLINE quotRem #-}--  rem (Poly xs) (Poly ys) = toPoly $ remainder xs ys-  {-# INLINE rem #-}--quotientAndRemainder-  :: (Fractional a, G.Vector v a)-  => v a-  -> v a-  -> (v a, v a)-quotientAndRemainder xs ys-  | G.null ys = throw DivideByZero-  | G.basicLength xs < G.basicLength ys = (G.empty, xs)-  | otherwise = runST $ do-    let lenXs = G.basicLength xs-        lenYs = G.basicLength ys-        lenQs = lenXs - lenYs + 1-    qs <- MG.basicUnsafeNew lenQs-    rs <- MG.basicUnsafeNew lenXs-    G.unsafeCopy rs xs-    forM_ [lenQs - 1, lenQs - 2 .. 0] $ \i -> do-      r <- MG.unsafeRead rs (lenYs - 1 + i)-      let q = r / G.unsafeLast ys-      MG.unsafeWrite qs i q-      forM_ [0 .. lenYs - 1] $ \k -> do-        MG.unsafeModify rs (\c -> c - q * G.unsafeIndex ys k) (i + k)-    let rs' = MG.basicUnsafeSlice 0 lenYs rs-    (,) <$> G.unsafeFreeze qs <*> G.unsafeFreeze rs'-{-# INLINE quotientAndRemainder #-}--remainder-  :: (Fractional a, G.Vector v a)-  => v a-  -> v a-  -> v a-remainder xs ys-  | G.null ys = throw DivideByZero-  | otherwise = runST $ do-    rs <- G.thaw xs-    ys' <- G.unsafeThaw ys-    remainderM rs ys'-    G.unsafeFreeze $ MG.basicUnsafeSlice 0 (G.basicLength xs `min` G.basicLength ys) rs-{-# INLINE remainder #-}--remainderM-  :: (PrimMonad m, Fractional 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.basicLength xs < MG.basicLength ys = pure ()-  | otherwise = do-    let lenXs = MG.basicLength xs-        lenYs = MG.basicLength ys-        lenQs = lenXs - lenYs + 1-    yLast <- MG.unsafeRead ys (lenYs - 1)-    forM_ [lenQs - 1, lenQs - 2 .. 0] $ \i -> do-      r <- MG.unsafeRead xs (lenYs - 1 + i)-      forM_ [0 .. lenYs - 1] $ \k -> do-        y <- MG.unsafeRead ys k-        -- do not move r / yLast outside the loop,-        -- because of numerical instability-        MG.unsafeModify xs (\c -> c - r * y / yLast) (i + k)-{-# INLINE remainderM #-}--fractionalGcd-  :: (Eq a, Fractional a, G.Vector v a)-  => Poly v a-  -> Poly v a-  -> Poly v a-fractionalGcd (Poly xs) (Poly ys) = toPoly $ runST $ do-  xs' <- G.thaw xs-  ys' <- G.thaw ys-  gcdM xs' ys'-{-# INLINE fractionalGcd #-}--gcdM-  :: (PrimMonad m, Eq a, Fractional a, G.Vector v a)-  => G.Mutable v (PrimState m) a-  -> G.Mutable v (PrimState m) a-  -> m (v a)-gcdM xs ys = do-  ys' <- dropWhileEndM (== 0) ys-  if MG.null ys' then G.unsafeFreeze xs else do-    remainderM xs ys'-    gcdM ys' xs-{-# INLINE gcdM #-}--#else--module Data.Poly.Internal.Dense.Fractional () where--#endif

@@ -9,11 +9,9 @@  {-# LANGUAGE CPP                        #-} {-# LANGUAGE FlexibleInstances          #-}-{-# LANGUAGE GeneralizedNewtypeDeriving #-} {-# LANGUAGE PatternSynonyms            #-} {-# LANGUAGE ScopedTypeVariables        #-} {-# LANGUAGE TypeFamilies               #-}-{-# LANGUAGE ViewPatterns               #-}  {-# OPTIONS_GHC -fno-warn-orphans #-} @@ -28,18 +26,16 @@ import Control.Monad.Primitive import Control.Monad.ST import Data.Euclidean-import Data.Semiring (Semiring(..), isZero)-import qualified Data.Semiring as Semiring+import Data.Semiring (Semiring(..), Ring(), isZero, minus) import qualified Data.Vector.Generic as G import qualified Data.Vector.Generic.Mutable as MG  import Data.Poly.Internal.Dense --- | Consider using 'Data.Poly.Semiring.PolyOverFractional' wrapper,+-- | Consider using 'Data.Poly.Semiring.PolyOverField' wrapper, -- which provides a much faster implementation of--- 'Data.Euclidean.gcd' for 'Fractional'--- coefficients.-instance (Eq a, Semiring.Ring a, GcdDomain a, Eq (v a), G.Vector v a) => GcdDomain (Poly v a) where+-- 'Data.Euclidean.gcd' for polynomials over 'Field'.+instance (Eq a, Ring a, GcdDomain a, Eq (v a), G.Vector v a) => GcdDomain (Poly v a) where   divide (Poly xs) (Poly ys) =     toPoly' <$> quotient xs ys @@ -50,7 +46,7 @@   {-# INLINE gcd #-}  gcdNonEmpty-  :: (Eq a, Semiring.Ring a, GcdDomain a, G.Vector v a)+  :: (Eq a, Ring a, GcdDomain a, G.Vector v a)   => v a   -> v a   -> v a@@ -62,7 +58,7 @@     ys' <- G.thaw ys     zs' <- gcdM xs' ys' -    let lenZs = MG.basicLength zs'+    let lenZs = MG.length zs'         go acc 0 = pure acc         go acc n = do           t <- MG.unsafeRead zs' (n - 1)@@ -80,7 +76,7 @@     G.unsafeFreeze zs'  gcdM-  :: (PrimMonad m, Eq a, Semiring.Ring a, GcdDomain a, G.Vector v a)+  :: (PrimMonad m, Eq a, Ring a, GcdDomain a, G.Vector v a)   => G.Mutable v (PrimState m) a   -> G.Mutable v (PrimState m) a   -> m (G.Mutable v (PrimState m) a)@@ -88,8 +84,8 @@   | MG.null xs = pure ys   | MG.null ys = pure xs   | otherwise = do-  let lenXs = MG.basicLength xs-      lenYs = MG.basicLength ys+  let lenXs = MG.length xs+      lenYs = MG.length ys   xLast <- MG.unsafeRead xs (lenXs - 1)   yLast <- MG.unsafeRead ys (lenYs - 1)   let z = xLast `lcm` yLast@@ -105,7 +101,7 @@       x <- MG.unsafeRead xs i       MG.unsafeModify         ys-        (\y -> (y `times` zy) `plus` Semiring.negate (x `times` zx))+        (\y -> (y `times` zy) `minus` x `times` zx)         (i + lenYs - lenXs)     forM_ [0 .. lenYs - lenXs - 1] $       MG.unsafeModify ys (`times` zy)@@ -116,7 +112,7 @@       y <- MG.unsafeRead ys i       MG.unsafeModify         xs-        (\x -> (x `times` zx) `plus` Semiring.negate (y `times` zy))+        (\x -> (x `times` zx) `minus` y `times` zy)         (i + lenXs - lenYs)     forM_ [0 .. lenXs - lenYs - 1] $       MG.unsafeModify xs (`times` zx)@@ -128,7 +124,7 @@   :: (Eq a, Semiring a, PrimMonad m, G.Vector v a)   => G.Mutable v (PrimState m) a   -> m Bool-isZeroM xs = go (MG.basicLength xs)+isZeroM xs = go (MG.length xs)   where     go 0 = pure True     go n = do@@ -137,20 +133,20 @@ {-# INLINE isZeroM #-}  quotient-  :: (Eq a, Eq (v a), Semiring.Ring a, GcdDomain a, G.Vector v a)+  :: (Eq a, Eq (v a), Ring a, GcdDomain a, G.Vector v a)   => v a   -> v a   -> Maybe (v a) quotient xs ys   | G.null ys = throw DivideByZero   | G.null xs = Just xs-  | G.basicLength xs < G.basicLength ys = Nothing+  | G.length xs < G.length ys = Nothing   | otherwise = runST $ do-    let lenXs = G.basicLength xs-        lenYs = G.basicLength ys+    let lenXs = G.length xs+        lenYs = G.length ys         lenQs = lenXs - lenYs + 1-    qs <- MG.basicUnsafeNew lenQs-    rs <- MG.basicUnsafeNew lenXs+    qs <- MG.unsafeNew lenQs+    rs <- MG.unsafeNew lenXs     G.unsafeCopy rs xs      let go i@@ -168,7 +164,7 @@                 forM_ [0 .. lenYs - 1] $ \k -> do                   MG.unsafeModify                     rs-                    (\c -> c `plus` (Semiring.negate $ q `times` G.unsafeIndex ys k))+                    (\c -> c `minus` q `times` G.unsafeIndex ys k)                     (i + k)                 go (i - 1) 

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

@@ -1,55 +0,0 @@--- |--- Module:      Data.Poly.Internal.PolyOverFractional--- Copyright:   (c) 2019 Andrew Lelechenko--- Licence:     BSD3--- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>------ Wrapper with a more efficient 'Euclidean' instance.-----{-# LANGUAGE CPP                        #-}-{-# LANGUAGE FlexibleInstances          #-}-{-# LANGUAGE GeneralizedNewtypeDeriving #-}-{-# LANGUAGE UndecidableInstances       #-}--#if MIN_VERSION_semirings(0,4,2)--module Data.Poly.Internal.PolyOverFractional-  ( PolyOverFractional(..)-  ) where--import Prelude hiding (quotRem, quot, rem, gcd, lcm, (^))-import Data.Euclidean-import Data.Semiring-import qualified Data.Semiring as Semiring-import qualified Data.Vector.Generic as G--import qualified Data.Poly.Internal.Dense as Dense-import qualified Data.Poly.Internal.Dense.Fractional as Dense (fractionalGcd)---- | Wrapper over polynomials,--- providing a faster 'GcdDomain' instance,--- when coefficients are 'Fractional'.-newtype PolyOverFractional poly = PolyOverFractional { unPolyOverFractional :: poly }-  deriving (Eq, Ord, Show, Num, Semiring, Semiring.Ring)--instance (Eq a, Eq (v a), Semiring.Ring a, GcdDomain a, Fractional a, G.Vector v a) => GcdDomain (PolyOverFractional (Dense.Poly v a)) where-  gcd (PolyOverFractional x) (PolyOverFractional y) = PolyOverFractional (Dense.fractionalGcd x y)-  {-# INLINE gcd #-}--instance (Eq a, Eq (v a), Semiring.Ring a, GcdDomain a, Fractional a, G.Vector v a) => Euclidean (PolyOverFractional (Dense.Poly v a)) where-  degree (PolyOverFractional x) =-    degree x-  quotRem (PolyOverFractional x) (PolyOverFractional y) =-    let (q, r) = quotRem x y in-      (PolyOverFractional q, PolyOverFractional r)-  {-# INLINE quotRem #-}-  rem (PolyOverFractional x) (PolyOverFractional y) =-    PolyOverFractional (rem x y)-  {-# INLINE rem #-}--#else--module Data.Poly.Internal.PolyOverFractional () where--#endif

@@ -7,14 +7,15 @@ -- Sparse polynomials of one variable. -- -{-# LANGUAGE CPP                  #-}-{-# LANGUAGE FlexibleContexts     #-}-{-# LANGUAGE PatternSynonyms      #-}-{-# LANGUAGE ScopedTypeVariables  #-}-{-# LANGUAGE StandaloneDeriving   #-}-{-# LANGUAGE TypeFamilies         #-}-{-# LANGUAGE UndecidableInstances #-}-{-# LANGUAGE ViewPatterns         #-}+{-# LANGUAGE CPP                        #-}+{-# LANGUAGE FlexibleContexts           #-}+{-# LANGUAGE GeneralizedNewtypeDeriving #-}+{-# LANGUAGE PatternSynonyms            #-}+{-# LANGUAGE ScopedTypeVariables        #-}+{-# LANGUAGE StandaloneDeriving         #-}+{-# LANGUAGE TypeFamilies               #-}+{-# LANGUAGE UndecidableInstances       #-}+{-# LANGUAGE ViewPatterns               #-}  module Data.Poly.Internal.Sparse   ( Poly(..)@@ -36,14 +37,20 @@   , pattern X'   , eval'   , deriv'+#if MIN_VERSION_semirings(0,5,0)+  , integral'+#endif   ) where +import Prelude hiding (quot)+import Control.DeepSeq (NFData) import Control.Monad import Control.Monad.Primitive import Control.Monad.ST+import Data.Bits import Data.List (intersperse) import Data.Ord-import Data.Semiring (Semiring(..))+import Data.Semiring (Semiring(..), Ring()) import qualified Data.Semiring as Semiring import qualified Data.Vector as V import qualified Data.Vector.Generic as G@@ -55,6 +62,9 @@ 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.).@@ -78,8 +88,9 @@   -- (first element corresponds to a constant term).   } -deriving instance Eq   (v (Word, a)) => Eq   (Poly v a)-deriving instance Ord  (v (Word, a)) => Ord  (Poly v a)+deriving instance Eq     (v (Word, a)) => Eq     (Poly v a)+deriving instance Ord    (v (Word, a)) => Ord    (Poly v a)+deriving instance NFData (v (Word, a)) => NFData (Poly v a)  instance (Eq a, Semiring a, G.Vector v (Word, a)) => IsList (Poly v a) where   type Item (Poly v a) = (Word, a)@@ -143,7 +154,7 @@   | otherwise = runST $ do     ws <- G.thaw vs     l' <- normalizeM p add ws-    G.unsafeFreeze $ MG.basicUnsafeSlice 0 l' ws+    G.unsafeFreeze $ MG.unsafeSlice 0 l' ws  normalizeM   :: (PrimMonad m, G.Vector v (Word, a))@@ -152,9 +163,9 @@   -> G.Mutable v (PrimState m) (Word, a)   -> m Int normalizeM p add ws = do-    let l = MG.basicLength ws+    let l = MG.length ws     let go i j acc@(accP, accC)-          | j >= l = do+          | j >= l =             if p accC               then do                 MG.write ws i acc@@ -181,7 +192,7 @@   abs = id   signum = const 1   fromInteger n = case fromInteger n of-    0 -> Poly $ G.empty+    0 -> Poly G.empty     m -> Poly $ G.singleton (0, m)   Poly xs * Poly ys = Poly $ convolution (/= 0) (+) (*) xs ys   {-# INLINE (+) #-}@@ -210,7 +221,7 @@   {-# INLINE fromNatural #-} #endif -instance (Eq a, Semiring.Ring a, G.Vector v (Word, a)) => Semiring.Ring (Poly v a) where+instance (Eq a, Ring a, G.Vector v (Word, a)) => Ring (Poly v a) where   negate (Poly xs) = Poly $ G.map (fmap Semiring.negate) xs  plusPoly@@ -221,9 +232,9 @@   -> v (Word, a)   -> v (Word, a) plusPoly p add xs ys = runST $ do-  zs <- MG.basicUnsafeNew (G.basicLength xs + G.basicLength ys)+  zs <- MG.unsafeNew (G.length xs + G.length ys)   lenZs <- plusPolyM p add xs ys zs-  G.unsafeFreeze $ MG.basicUnsafeSlice 0 lenZs zs+  G.unsafeFreeze $ MG.unsafeSlice 0 lenZs zs {-# INLINE plusPoly #-}  plusPolyM@@ -236,20 +247,20 @@   -> m Int plusPolyM p add xs ys zs = go 0 0 0   where-    lenXs = G.basicLength xs-    lenYs = G.basicLength ys+    lenXs = G.length xs+    lenYs = G.length ys      go ix iy iz       | ix == lenXs, iy == lenYs = pure iz       | ix == lenXs = do         G.unsafeCopy-          (MG.basicUnsafeSlice iz (lenYs - iy) zs)-          (G.basicUnsafeSlice iy (lenYs - iy) ys)+          (MG.unsafeSlice iz (lenYs - iy) zs)+          (G.unsafeSlice iy (lenYs - iy) ys)         pure $ iz + lenYs - iy       | iy == lenYs = do         G.unsafeCopy-          (MG.basicUnsafeSlice iz (lenXs - ix) zs)-          (G.basicUnsafeSlice ix (lenXs - ix) xs)+          (MG.unsafeSlice iz (lenXs - ix) zs)+          (G.unsafeSlice ix (lenXs - ix) xs)         pure $ iz + lenXs - ix       | (xp, xc) <- G.unsafeIndex xs ix       , (yp, yc) <- G.unsafeIndex ys iy@@ -278,7 +289,7 @@   -> v (Word, a)   -> v (Word, a) minusPoly p neg sub xs ys = runST $ do-  zs <- MG.basicUnsafeNew (lenXs + lenYs)+  zs <- MG.unsafeNew (lenXs + lenYs)   let go ix iy iz         | ix == lenXs, iy == lenYs = pure iz         | ix == lenXs = do@@ -288,8 +299,8 @@           pure $ iz + lenYs - iy         | iy == lenYs = do           G.unsafeCopy-            (MG.basicUnsafeSlice iz (lenXs - ix) zs)-            (G.basicUnsafeSlice ix (lenXs - ix) xs)+            (MG.unsafeSlice iz (lenXs - ix) zs)+            (G.unsafeSlice ix (lenXs - ix) xs)           pure $ iz + lenXs - ix         | (xp, xc) <- G.unsafeIndex xs ix         , (yp, yc) <- G.unsafeIndex ys iy@@ -308,10 +319,10 @@             MG.unsafeWrite zs iz (yp, neg yc)             go ix (iy + 1) (iz + 1)   lenZs <- go 0 0 0-  G.unsafeFreeze $ MG.basicUnsafeSlice 0 lenZs zs+  G.unsafeFreeze $ MG.unsafeSlice 0 lenZs zs   where-    lenXs = G.basicLength xs-    lenYs = G.basicLength ys+    lenXs = G.length xs+    lenYs = G.length ys {-# INLINE minusPoly #-}  scaleM@@ -324,7 +335,7 @@   -> m Int scaleM p mul xs (yp, yc) zs = go 0 0   where-    lenXs = G.basicLength xs+    lenXs = G.length xs      go ix iz       | ix == lenXs = pure iz@@ -347,9 +358,9 @@   -> Poly v a   -> Poly v a scaleInternal p mul yp yc (Poly xs) = runST $ do-  zs <- MG.basicUnsafeNew (G.basicLength xs)+  zs <- MG.unsafeNew (G.length xs)   len <- scaleM p (flip mul) xs (yp, yc) zs-  fmap Poly $ G.unsafeFreeze $ MG.basicUnsafeSlice 0 len zs+  fmap Poly $ G.unsafeFreeze $ MG.unsafeSlice 0 len zs {-# INLINE scaleInternal #-}  -- | Multiply a polynomial by a monomial, expressed as a power and a coefficient.@@ -372,29 +383,29 @@   -> v (Word, a)   -> v (Word, a) convolution p add mult xs ys-  | G.basicLength xs >= G.basicLength ys+  | G.length xs >= G.length ys   = go mult xs ys   | otherwise   = go (flip mult) ys xs   where     go :: (a -> a -> a) -> v (Word, a) -> v (Word, a) -> v (Word, a)     go mul long short = runST $ do-      let lenLong   = G.basicLength long-          lenShort  = G.basicLength short+      let lenLong   = G.length long+          lenShort  = G.length short           lenBuffer = lenLong * lenShort-      slices <- MG.basicUnsafeNew lenShort-      buffer <- MG.basicUnsafeNew lenBuffer+      slices <- MG.unsafeNew lenShort+      buffer <- MG.unsafeNew lenBuffer        forM_ [0 .. lenShort - 1] $ \iShort -> do         let (pShort, cShort) = G.unsafeIndex short iShort             from = iShort * lenLong-            bufferSlice = MG.basicUnsafeSlice from lenLong buffer+            bufferSlice = MG.unsafeSlice from lenLong buffer         len <- scaleM p mul long (pShort, cShort) bufferSlice         MG.unsafeWrite slices iShort (from, len)        slices' <- G.unsafeFreeze slices       buffer' <- G.unsafeFreeze buffer-      bufferNew <- MG.basicUnsafeNew lenBuffer+      bufferNew <- MG.unsafeNew lenBuffer       gogo slices' buffer' bufferNew      gogo@@ -404,29 +415,29 @@       -> G.Mutable v (PrimState m) (Word, a)       -> m (v (Word, a))     gogo slices buffer bufferNew-      | G.basicLength slices == 0+      | G.length slices == 0       = pure G.empty-      | G.basicLength slices == 1+      | G.length slices == 1       , (from, len) <- G.unsafeIndex slices 0-      = pure $ G.basicUnsafeSlice from len buffer+      = pure $ G.unsafeSlice from len buffer       | otherwise = do-        let nSlices = G.basicLength slices-        slicesNew <- MG.basicUnsafeNew ((nSlices + 1) `quot` 2)-        forM_ [0 .. (nSlices - 2) `quot` 2] $ \i -> do+        let nSlices = G.length slices+        slicesNew <- MG.unsafeNew ((nSlices + 1) `shiftR` 1)+        forM_ [0 .. (nSlices - 2) `shiftR` 1] $ \i -> do           let (from1, len1) = G.unsafeIndex slices (2 * i)               (from2, len2) = G.unsafeIndex slices (2 * i + 1)-              slice1 = G.basicUnsafeSlice from1 len1 buffer-              slice2 = G.basicUnsafeSlice from2 len2 buffer-              slice3 = MG.basicUnsafeSlice from1 (len1 + len2) bufferNew+              slice1 = G.unsafeSlice from1 len1 buffer+              slice2 = G.unsafeSlice from2 len2 buffer+              slice3 = MG.unsafeSlice from1 (len1 + len2) bufferNew           len3 <- plusPolyM p add slice1 slice2 slice3           MG.unsafeWrite slicesNew i (from1, len3)          when (odd nSlices) $ do           let (from, len) = G.unsafeIndex slices (nSlices - 1)-              slice1 = G.basicUnsafeSlice from len buffer-              slice3 = MG.basicUnsafeSlice from len bufferNew+              slice1 = G.unsafeSlice from len buffer+              slice3 = MG.unsafeSlice from len bufferNew           G.unsafeCopy slice3 slice1-          MG.unsafeWrite slicesNew (nSlices `quot` 2) (from, len)+          MG.unsafeWrite slicesNew (nSlices `shiftR` 1) (from, len)          slicesNew' <- G.unsafeFreeze slicesNew         buffer'    <- G.unsafeThaw   buffer@@ -509,8 +520,8 @@ derivPoly p mul xs   | G.null xs = G.empty   | otherwise = runST $ do-    let lenXs = G.basicLength xs-    zs <- MG.basicUnsafeNew lenXs+    let lenXs = G.length xs+    zs <- MG.unsafeNew lenXs     let go ix iz           | ix == lenXs = pure iz           | (xp, xc) <- G.unsafeIndex xs ix@@ -522,7 +533,7 @@             else               go (ix + 1) iz     lenZs <- go 0 0-    G.unsafeFreeze $ MG.basicUnsafeSlice 0 lenZs zs+    G.unsafeFreeze $ MG.unsafeSlice 0 lenZs zs {-# INLINE derivPoly #-}  -- | Compute an indefinite integral of a polynomial,@@ -535,6 +546,14 @@   = Poly   $ 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

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

@@ -1,78 +0,0 @@--- |--- Module:      Data.Poly.Internal.Sparse.Fractional--- Copyright:   (c) 2019 Andrew Lelechenko--- Licence:     BSD3--- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>------ GcdDomain for Fractional underlying-----{-# LANGUAGE CPP                        #-}-{-# LANGUAGE FlexibleContexts           #-}-{-# LANGUAGE FlexibleInstances          #-}-{-# LANGUAGE GeneralizedNewtypeDeriving #-}-{-# LANGUAGE PatternSynonyms            #-}-{-# LANGUAGE ScopedTypeVariables        #-}-{-# LANGUAGE TypeFamilies               #-}-{-# LANGUAGE UndecidableInstances       #-}-{-# LANGUAGE ViewPatterns               #-}--{-# OPTIONS_GHC -fno-warn-orphans #-}--#if MIN_VERSION_semirings(0,4,2)--module Data.Poly.Internal.Sparse.Fractional-  ( fractionalGcd-  ) where--import Prelude hiding (quotRem, rem, gcd)-import Control.Arrow-import Control.Exception-import Data.Euclidean-import qualified Data.Semiring as Semiring-import qualified Data.Vector.Generic as G--import Data.Poly.Internal.Sparse-import Data.Poly.Internal.Sparse.GcdDomain ()--instance (Eq a, Eq (v (Word, a)), Semiring.Ring a, GcdDomain a, Fractional a, G.Vector v (Word, a)) => Euclidean (Poly v a) where-  degree (Poly xs)-    | G.null xs = 0-    | otherwise = 1 + fromIntegral (fst (G.last xs))--  quotRem = quotientRemainder--quotientRemainder-  :: (Eq a, Fractional a, G.Vector v (Word, a))-  => Poly v a-  -> Poly v a-  -> (Poly v a, Poly v a)-quotientRemainder ts ys = case leading ys of-  Nothing -> throw DivideByZero-  Just (yp, yc) -> go ts-    where-      go xs = case leading xs of-        Nothing -> (0, 0)-        Just (xp, xc) -> case xp `compare` yp of-          LT -> (0, xs)-          EQ -> (zs, xs')-          GT -> first (+ zs) $ go xs'-          where-            zs = Poly $ G.singleton (xp - yp, xc / yc)-            xs' = xs - zs * ys--fractionalGcd-  :: (Eq a, Fractional a, G.Vector v (Word, a))-  => Poly v a-  -> Poly v a-  -> Poly v a-fractionalGcd xs ys-  | G.null (unPoly ys) = xs-  | otherwise = fractionalGcd ys $ snd $ quotientRemainder xs ys-{-# INLINE fractionalGcd #-}--#else--module Data.Poly.Internal.Sparse.Fractional () where--#endif

@@ -10,12 +10,10 @@ {-# LANGUAGE CPP                        #-} {-# LANGUAGE FlexibleContexts           #-} {-# LANGUAGE FlexibleInstances          #-}-{-# LANGUAGE GeneralizedNewtypeDeriving #-} {-# LANGUAGE PatternSynonyms            #-} {-# LANGUAGE ScopedTypeVariables        #-} {-# LANGUAGE TypeFamilies               #-} {-# LANGUAGE UndecidableInstances       #-}-{-# LANGUAGE ViewPatterns               #-}  {-# OPTIONS_GHC -fno-warn-orphans #-} @@ -28,17 +26,15 @@ import Control.Exception import Data.Euclidean import Data.Maybe-import Data.Semiring (Semiring(..))-import qualified Data.Semiring as Semiring+import Data.Semiring (Semiring(..), Ring(), minus) import qualified Data.Vector.Generic as G  import Data.Poly.Internal.Sparse --- | Consider using 'Data.Poly.Sparse.Semiring.PolyOverFractional' wrapper,+-- | Consider using 'Data.Poly.Sparse.Semiring.PolyOverField' wrapper, -- which provides a much faster implementation of--- 'Data.Euclidean.gcd' for 'Fractional'--- coefficients.-instance (Eq a, Semiring.Ring a, GcdDomain a, Eq (v (Word, a)), G.Vector v (Word, a)) => GcdDomain (Poly v a) where+-- 'Data.Euclidean.gcd' for polynomials over 'Field'.+instance (Eq a, Ring a, GcdDomain a, Eq (v (Word, a)), G.Vector v (Word, a)) => GcdDomain (Poly v a) where   divide xs ys = case leading ys of     Nothing -> throw DivideByZero     Just (yp, yc) -> case leading xs of@@ -48,7 +44,7 @@         | otherwise -> do           zc <- divide xc yc           let z = Poly $ G.singleton (xp - yp, zc)-          rest <- divide (xs `plus` Semiring.negate z `times` ys) ys+          rest <- divide (xs `minus` z `times` ys) ys           pure $ rest `plus` z    gcd xs ys@@ -62,7 +58,7 @@         xy = monomial' 0 (gcd (cont xs) (cont ys))  gcdHelper-  :: (Eq a, Semiring.Ring a, GcdDomain a, G.Vector v (Word, a))+  :: (Eq a, Ring a, GcdDomain a, G.Vector v (Word, a))   => Poly v a   -> Poly v a   -> Poly v a@@ -71,9 +67,9 @@   Just (xp, xc) -> case leading ys of     Nothing -> xs     Just (yp, yc) -> case xp `compare` yp of-      LT -> gcdHelper xs (ys `times` monomial' 0 gy `plus` Semiring.negate (xs `times` monomial' (yp - xp) gx))-      EQ -> gcdHelper xs (ys `times` monomial' 0 gy `plus` Semiring.negate (xs `times` monomial' 0 gx))-      GT -> gcdHelper (xs `times` monomial' 0 gx `plus` Semiring.negate (ys `times` monomial' (xp - yp) gy)) ys+      LT -> gcdHelper xs (ys `times` monomial' 0 gy `minus` xs `times` monomial' (yp - xp) gx)+      EQ -> gcdHelper xs (ys `times` monomial' 0 gy `minus` xs `times` monomial' 0 gx)+      GT -> gcdHelper (xs `times` monomial' 0 gx `minus` ys `times` monomial' (xp - yp) gy) ys       where         g = lcm xc yc         gx = fromMaybe err $ divide g xc

@@ -7,8 +7,8 @@ -- Dense polynomials and a 'Semiring'-based interface. -- -{-# LANGUAGE CPP                 #-}-{-# LANGUAGE PatternSynonyms     #-}+{-# LANGUAGE CPP             #-}+{-# LANGUAGE PatternSynonyms #-}  module Data.Poly.Semiring   ( Poly@@ -23,9 +23,16 @@   , pattern X   , eval   , deriv+#if MIN_VERSION_semirings(0,5,0)+  , integral+#endif #if MIN_VERSION_semirings(0,4,2)-  -- * Fractional coefficients-  , PolyOverFractional(..)+  -- * Polynomials over 'Field'+  , PolyOverField(..)+  , gcdExt+  , PolyOverFractional+  , pattern PolyOverFractional+  , unPolyOverFractional #endif   ) where @@ -35,10 +42,13 @@ 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.Fractional ()+import Data.Poly.Internal.Dense.Field (gcdExt) import Data.Poly.Internal.Dense.GcdDomain ()-import Data.Poly.Internal.PolyOverFractional+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).@@ -81,3 +91,13 @@ -- 3 * X^2 + 0 * X + 3 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.+--+-- >>> integral (3 * X^2 + 3) :: UPoly Double+-- 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

@@ -7,7 +7,8 @@ -- Sparse polynomials with 'Num' instance. -- -{-# LANGUAGE PatternSynonyms     #-}+{-# LANGUAGE CPP             #-}+{-# LANGUAGE PatternSynonyms #-}  module Data.Poly.Sparse   ( Poly@@ -23,8 +24,14 @@   , eval   , deriv   , integral+#if MIN_VERSION_semirings(0,4,2)+  -- * Polynomials over 'Field'+  , gcdExt+#endif   ) where  import Data.Poly.Internal.Sparse-import Data.Poly.Internal.Sparse.Fractional ()+#if MIN_VERSION_semirings(0,4,2)+import Data.Poly.Internal.Sparse.Field (gcdExt) import Data.Poly.Internal.Sparse.GcdDomain ()+#endif

@@ -7,8 +7,9 @@ -- Sparse polynomials with 'Semiring' instance. -- -{-# LANGUAGE FlexibleContexts    #-}-{-# LANGUAGE PatternSynonyms     #-}+{-# LANGUAGE CPP              #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE PatternSynonyms  #-}  module Data.Poly.Sparse.Semiring   ( Poly@@ -23,6 +24,13 @@   , pattern X   , eval   , 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.Semiring (Semiring)@@ -30,8 +38,13 @@  import Data.Poly.Internal.Sparse (Poly(..), VPoly, UPoly, leading) import qualified Data.Poly.Internal.Sparse as Sparse-import Data.Poly.Internal.Sparse.Fractional ()+#if MIN_VERSION_semirings(0,4,2)+import Data.Poly.Internal.Sparse.Field (gcdExt) 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).@@ -74,3 +87,13 @@ -- 3 * X^2 + 3 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.+--+-- >>> integral (3 * X^2 + 3) :: UPoly Double+-- 1.0 * X^3 + 3.0 * X+integral :: (Eq a, Field a, G.Vector v (Word, a)) => Poly v a -> Poly v a+integral = Sparse.integral'+#endif

@@ -10,11 +10,12 @@   ( testSuite   ) where -import Prelude hiding (quotRem)+import Prelude hiding (gcd, quotRem, rem) #if MIN_VERSION_semirings(0,4,2) import Data.Euclidean #endif import Data.Int+import Data.Maybe import Data.Poly import qualified Data.Poly.Semiring as S import Data.Proxy@@ -23,7 +24,7 @@ import qualified Data.Vector.Generic as G import qualified Data.Vector.Unboxed as U import Test.Tasty-import Test.Tasty.QuickCheck hiding (scale)+import Test.Tasty.QuickCheck hiding (scale, numTests) import Test.QuickCheck.Classes  import Quaternion@@ -33,9 +34,9 @@   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 (PolyOverFractional (Poly v a)) where-  arbitrary = PolyOverFractional . S.toPoly . G.fromList . (\xs -> take (length xs `mod` 10) xs) <$> arbitrary-  shrink = fmap (PolyOverFractional . S.toPoly . G.fromList) . shrink . G.toList . unPoly . unPolyOverFractional+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 }@@ -57,14 +58,24 @@ testSuite = testGroup "Dense"     [ arithmeticTests     , otherTests-    , semiringTests+    , lawsTests     , evalTests     , derivTests #if MIN_VERSION_semirings(0,4,2)-    -- , euclideanTests+    , gcdExtTests #endif     ] +lawsTests :: TestTree+lawsTests = testGroup "Laws"+  [ semiringTests+  , ringTests+  , numTests+  , euclideanTests+  , isListTests+  , showTests+  ]+ semiringTests :: TestTree semiringTests   = testGroup "Semiring"@@ -74,26 +85,73 @@   , 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+  = 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 ()))   , ringLaws (Proxy :: Proxy (Poly U.Vector Int8))   , ringLaws (Proxy :: Proxy (Poly V.Vector Integer))   , ringLaws (Proxy :: Proxy (Poly U.Vector (Quaternion Int))) #endif   ] -#if MIN_VERSION_semirings(0,4,2)--- euclideanTests :: TestTree--- euclideanTests---   = testGroup "Euclidean"---   $ map (uncurry testProperty)---   $ concatMap lawsProperties---   [ gcdDomainLaws (Proxy :: Proxy (ShortPoly (Poly V.Vector Integer)))---   , gcdDomainLaws (Proxy :: Proxy (PolyOverFractional (Poly V.Vector Rational)))---   , euclideanLaws (Proxy :: Proxy (ShortPoly (Poly V.Vector Rational)))---   ]+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+  ] +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+  ]++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)))+  ]++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+  ]+ arithmeticTests :: TestTree arithmeticTests = testGroup "Arithmetic"   [ testProperty "addition matches reference" $@@ -189,6 +247,8 @@ derivTests = testGroup "deriv"   [ testProperty "deriv = S.deriv" $     \(p :: Poly V.Vector Integer) -> deriv p === S.deriv p+  , 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   , testProperty "deriv c = 0" $@@ -205,3 +265,22 @@   --     deriv (eval (toPoly $ fmap (monomial 0) $ unPoly p) q) ===   --       deriv q * eval (toPoly $ fmap (monomial 0) $ unPoly $ 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

@@ -11,13 +11,14 @@   ( testSuite   ) where -import Prelude hiding (quotRem)+import Prelude hiding (gcd, quotRem, rem) #if MIN_VERSION_semirings(0,4,2) import Data.Euclidean #endif import Data.Function import Data.Int import Data.List+import Data.Maybe import Data.Poly.Sparse import qualified Data.Poly.Sparse.Semiring as S import Data.Proxy@@ -26,7 +27,7 @@ import qualified Data.Vector.Generic as G import qualified Data.Vector.Unboxed as U import Test.Tasty-import Test.Tasty.QuickCheck hiding (scale)+import Test.Tasty.QuickCheck hiding (scale, numTests) import Test.QuickCheck.Classes  import Quaternion@@ -54,11 +55,24 @@ testSuite = testGroup "Sparse"     [ arithmeticTests     , otherTests-    , semiringTests+    , lawsTests     , evalTests     , derivTests+#if MIN_VERSION_semirings(0,4,2)+    , gcdExtTests+#endif     ] +lawsTests :: TestTree+lawsTests = testGroup "Laws"+  [ semiringTests+  , ringTests+  , numTests+  , euclideanTests+  , isListTests+  , showTests+  ]+ semiringTests :: TestTree semiringTests   = testGroup "Semiring"@@ -68,14 +82,72 @@   , 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+  = 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 ()))   , 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+  = 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+  ]++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+  ]++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)))+  ]++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+  ]+ arithmeticTests :: TestTree arithmeticTests = testGroup "Arithmetic"   [ testProperty "addition matches reference" $@@ -180,6 +252,8 @@ derivTests = testGroup "deriv"   [ testProperty "deriv = S.deriv" $     \(p :: Poly V.Vector Integer) -> deriv p === S.deriv p+  , 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   , testProperty "deriv c = 0" $@@ -196,3 +270,22 @@   --     deriv (eval (toPoly $ fmap (fmap $ monomial 0) $ unPoly p) q) ===   --       deriv q * eval (toPoly $ fmap (fmap $ monomial 0) $ unPoly $ 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