poly 0.2.0.0 → 0.3.0.0

raw patch · 21 files changed

+2201/−494 lines, 21 filesdep +gaugedep +vector-algorithmsdep ~basedep ~semiringsdep ~vectorPVP ok

version bump matches the API change (PVP)

Dependencies added: gauge, vector-algorithms

Dependency ranges changed: base, semirings, vector

API changes (from Hackage documentation)

Files

• README.md +107/−3
• bench/Bench.hs +13/−0
• bench/DenseBench.hs +94/−0
• bench/SparseBench.hs +70/−0
• changelog.md +7/−0
• poly.cabal +33/−5
• src/Data/Poly.hs +12/−3
• src/Data/Poly/Internal/Dense.hs +364/−0
• src/Data/Poly/Internal/Dense/Fractional.hs +129/−0
• src/Data/Poly/Internal/Dense/GcdDomain.hs +173/−0
• src/Data/Poly/Internal/PolyOverFractional.hs +46/−0
• src/Data/Poly/Internal/Sparse.hs +532/−0
• src/Data/Poly/Internal/Sparse/Fractional.hs +69/−0
• src/Data/Poly/Internal/Sparse/GcdDomain.hs +78/−0
• src/Data/Poly/Semiring.hs +20/−6
• src/Data/Poly/Sparse.hs +30/−0
• src/Data/Poly/Sparse/Semiring.hs +76/−0
• src/Data/Poly/Uni/Dense.hs +0/−356
• test/Dense.hs +173/−0
• test/Main.hs +4/−121
• test/Sparse.hs +171/−0

@@ -1,6 +1,6 @@ # poly [![Build Status](https://travis-ci.org/Bodigrim/poly.svg)](https://travis-ci.org/Bodigrim/poly) [![Hackage](http://img.shields.io/hackage/v/poly.svg)](https://hackage.haskell.org/package/poly) -Polynomials with `Num` and `Semiring` instances, backed by `Vector`.+Univariate polynomials, backed by `Vector`.  ```haskell > (X + 1) + (X - 1) :: VPoly Integer@@ -8,13 +8,117 @@  > (X + 1) * (X - 1) :: UPoly Int 1 * X^2 + 0 * X + (-1)+``` +## Vectors++`Poly v a` is polymorphic over a container `v`, implementing `Vector` interface, and coefficients of type `a`. Usually `v` is either a boxed vector from `Data.Vector` or an unboxed vector from `Data.Vector.Unboxed`. Use unboxed vectors whenever possible, e. g., when coefficients are `Int` or `Double`.++There are handy type synonyms:++```haskell+type VPoly a = Poly Data.Vector.Vector         a+type UPoly a = Poly Data.Vector.Unboxed.Vector a+```++## Construction++The simplest way to construct a polynomial is using the pattern `X`:++```haskell+> X^2 - 3*X + 2 :: UPoly Int+1 * X^2 + (-3) * X + 2+```++(Unfortunately, a type is 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):++```haskell+> toPoly (Data.Vector.Unboxed.fromList [2, -3, 1 :: Int])+1 * X^2 + (-3) * X + 2+```++There is a shortcut to construct a monomial:++```haskell+> monomial 2 3 :: UPoly Int+3 * X^2 + 0 * X + 0+```++## Operations++Most operations are provided by means of instances, like `Eq` and `Num`. For example,++```haskell+> (X^2 + 1) * (X^2 - 1) :: UPoly Int+1 * X^4 + 0 * X^3 + 0 * X^2 + 0 * X + (-1)+```++One can also find convenient to `scale` by monomial (cf. `monomial` above):++```haskell+> scale 2 3 (X^2 + 1) :: UPoly Int+3 * X^4 + 0 * X^3 + 3 * X^2 + 0 * X + 0+```++While `Poly` cannot be made an instance of `Integral` (because there is no meaningful `toInteger`),+it is an instance of `GcdDomain` and `Euclidean` from `semirings` package. These type classes+cover main functionality of `Integral`, providing division with remainder and `gcd` / `lcm`:++```haskell+> Data.Euclidean.gcd (X^2 + 7 * X + 6) (X^2 - 5 * X - 6) :: Data.Poly.UPoly Int+1 * X + 1++> Data.Euclidean.quotRem (X^3 + 2) (X^2 - 1 :: Data.Poly.UPoly Double)+(1.0 * X + 0.0,1.0 * X + 2.0)+```++Miscellaneous utilities include `eval` for evaluation at a given value of indeterminate,+and reciprocals `deriv` / `integral`:++```haskell > eval (X^2 + 1 :: UPoly Int) 3 10  > eval (X^2 + 1 :: VPoly (UPoly Int)) (X + 1) 1 * X^2 + 2 * X + 2 -> deriv (X^3 + 3 * X) :: UPoly Int-3 * X^2 + 0 * X + 3+> deriv (X^3 + 3 * X) :: UPoly Double+3.0 * X^2 + 0.0 * X + 3.0++> integral (3 * X^2 + 3) :: UPoly Double+1.0 * X^3 + 0.0 * X^2 + 3.0 * X + 0.0 ```++## Deconstruction++Use `unPoly` to deconstruct a polynomial to a vector of coefficients (head is the constant term):++```haskell+> unPoly (X^2 - 3 * X + 2 :: UPoly Int)+[2,-3,1]+```++Further, `leading` is a shortcut to to obtain the leading term of a non-zero polynomial,+expressed as a power and a coefficient:++```haskell+> leading (X^2 - 3 * X + 2 :: UPoly Int)+Just (2,1)+```++## Flavours++The same API is exposed in four flavours:++* `Data.Poly` provides dense polynomials with `Num`-based interface.+  This is a default choice for most users.++* `Data.Poly.Semiring` provides dense polynomials with `Semiring`-based interface.++* `Data.Poly.Sparse` provides sparse polynomials with `Num`-based interface.+  Besides that, you may find it easier to use in REPL+  because of a more readable `Show` instance, skipping zero coefficients.++* `Data.Poly.Sparse.Semiring` provides sparse polynomials with `Semiring`-based interface.

@@ -0,0 +1,13 @@+{-# LANGUAGE RankNTypes #-}++module Main where++import Gauge.Main+import qualified DenseBench as Dense+import qualified SparseBench as Sparse++main :: IO ()+main = defaultMain+  [ Dense.benchSuite+  , Sparse.benchSuite+  ]

@@ -0,0 +1,94 @@+{-# LANGUAGE RankNTypes #-}++module DenseBench+  ( benchSuite+  ) where++import Prelude hiding (quotRem, gcd)+import Gauge.Main+import Data.Euclidean+import Data.Poly+import qualified Data.Vector as V+import qualified Data.Vector.Unboxed as U++benchSuite :: Benchmark+benchSuite = bgroup "dense" $ concat+  [ map benchAdd      [100, 1000, 10000]+  , map benchMul      [10, 100]+  , map benchEval     [100, 1000, 10000]+  , map benchDeriv    [100, 1000, 10000]+  , map benchIntegral [100, 1000, 10000]+  , map benchQuotRem  [10, 100]+  , map benchGcdFrac  [10, 100]+  , map benchGcd      [10, 100]+  ]++benchAdd :: Int -> Benchmark+benchAdd k = bench ("add/" ++ show k) $ nf (doBinOp (+)) k++benchMul :: Int -> Benchmark+benchMul k = bench ("mul/" ++ show k) $ nf (doBinOp (*)) k++benchEval :: Int -> Benchmark+benchEval k = bench ("eval/" ++ show k) $ nf doEval k++benchDeriv :: Int -> Benchmark+benchDeriv k = bench ("deriv/" ++ show k) $ nf doDeriv k++benchIntegral :: Int -> Benchmark+benchIntegral k = bench ("integral/" ++ show k) $ nf doIntegral k++benchQuotRem :: Int -> Benchmark+benchQuotRem k = bench ("quotRem/" ++ show k) $ nf doQuotRem k++benchGcd :: Int -> Benchmark+benchGcd k = bench ("gcd/" ++ show k) $ nf doGcd k++benchGcdFrac :: Int -> Benchmark+benchGcdFrac k = bench ("gcdFrac/" ++ show k) $ nf doGcdFrac k++doBinOp :: (forall a. Num a => a -> a -> a) -> Int -> Int+doBinOp op n = U.sum zs+  where+    xs = toPoly $ U.generate n (* 2)+    ys = toPoly $ U.generate n (* 3)+    zs = unPoly $ xs `op` ys+{-# INLINE doBinOp #-}++doEval :: Int -> Int+doEval n = eval xs n+  where+    xs = toPoly $ U.generate n (* 2)++doDeriv :: Int -> Int+doDeriv n = U.sum zs+  where+    xs = toPoly $ U.generate n (* 2)+    zs = unPoly $ deriv xs++doIntegral :: Int -> Double+doIntegral n = U.sum zs+  where+    xs = toPoly $ U.generate n ((* 2) . fromIntegral)+    zs = unPoly $ integral xs++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)+    (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)+    gs = unPoly $ xs `gcd` ys++doGcdFrac :: Int -> Rational+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

@@ -0,0 +1,70 @@+{-# LANGUAGE RankNTypes #-}++{-# OPTIONS_GHC -fno-warn-type-defaults #-}++module SparseBench+  ( benchSuite+  ) where++import Gauge.Main+import Data.Poly.Sparse+import qualified Data.Vector.Unboxed as U++benchSuite :: Benchmark+benchSuite = bgroup "sparse" $ concat+  [ map benchAdd      $ zip3 tabs vecs2 vecs3+  , map benchMul      $ take 2 $ zip3 tabs vecs2 vecs3+  , map benchEval     $ zip tabs vecs2+  , map benchDeriv    $ zip tabs vecs2+  , map benchIntegral $ zip tabs vecs2'+  ]++tabs :: [Int]+tabs = [10, 100, 1000, 10000]++vecs2 :: [UPoly Int]+vecs2 = flip map tabs $+  \n -> toPoly $ U.generate n (\i -> (fromIntegral i ^ 2, i * 2))++vecs2' :: [UPoly Double]+vecs2' = flip map tabs $+  \n -> toPoly $ U.generate n (\i -> (fromIntegral i ^ 2, fromIntegral i * 2))++vecs3 :: [UPoly Int]+vecs3 = flip map tabs $+  \n -> toPoly $ U.generate n (\i -> (fromIntegral i ^ 3, i * 3))++benchAdd :: (Int, UPoly Int, UPoly Int) -> Benchmark+benchAdd (k, xs, ys) = bench ("add/" ++ show k) $ nf (doBinOp (+) xs) ys++benchMul :: (Int, UPoly Int, UPoly Int) -> Benchmark+benchMul (k, xs, ys) = bench ("mul/" ++ show k) $ nf (doBinOp (*) xs) ys++benchEval :: (Int, UPoly Int) -> Benchmark+benchEval (k, xs) = bench ("eval/" ++ show k) $ nf doEval xs++benchDeriv :: (Int, UPoly Int) -> Benchmark+benchDeriv (k, xs) = bench ("deriv/" ++ show k) $ nf doDeriv xs++benchIntegral :: (Int, UPoly Double) -> Benchmark+benchIntegral (k, xs) = bench ("integral/" ++ show k) $ nf doIntegral xs++doBinOp :: (forall a. Num a => a -> a -> a) -> UPoly Int -> UPoly Int -> Int+doBinOp op xs ys = U.foldl' (\acc (_, x) -> acc + x) 0 zs+  where+    zs = unPoly $ xs `op` ys+{-# INLINE doBinOp #-}++doEval :: UPoly Int -> Int+doEval xs = eval xs (U.length (unPoly xs))++doDeriv :: UPoly Int -> Int+doDeriv xs = U.foldl' (\acc (_, x) -> acc + x) 0 zs+  where+    zs = unPoly $ deriv xs++doIntegral :: UPoly Double -> Double+doIntegral xs = U.foldl' (\acc (_, x) -> acc + x) 0 zs+  where+    zs = unPoly $ integral xs+

@@ -1,3 +1,10 @@+# 0.3.0.0++* Implement sparse polynomials.+* Add `GcdDomain` and `Euclidean` instances.+* Add functions `leading`, `monomial`, `scale`.+* Remove function `constant`.+ # 0.2.0.0  * Fix a bug in `Num.(-)`.

@@ -1,8 +1,8 @@ name: poly-version: 0.2.0.0+version: 0.3.0.0 synopsis: Polynomials description:-  Polynomials with `Num` and `Semiring` instances, backed by `Vector`.+  Polynomials backed by `Vector`. homepage: https://github.com/Bodigrim/poly#readme license: BSD3 license-file: LICENSE@@ -26,19 +26,31 @@   exposed-modules:     Data.Poly     Data.Poly.Semiring+    Data.Poly.Sparse+    Data.Poly.Sparse.Semiring   other-modules:-    Data.Poly.Uni.Dense+    Data.Poly.Internal.Dense+    Data.Poly.Internal.Dense.Fractional+    Data.Poly.Internal.Dense.GcdDomain+    Data.Poly.Internal.PolyOverFractional+    Data.Poly.Internal.Sparse+    Data.Poly.Internal.Sparse.Fractional+    Data.Poly.Internal.Sparse.GcdDomain   build-depends:     base >= 4.9 && < 5,     primitive,-    semirings,-    vector+    semirings >= 0.4,+    vector,+    vector-algorithms   default-language: Haskell2010   ghc-options: -Wall  test-suite poly-tests   type: exitcode-stdio-1.0   main-is: Main.hs+  other-modules:+    Dense+    Sparse   build-depends:     base >=4.9 && <5,     poly,@@ -50,4 +62,20 @@     vector   default-language: Haskell2010   hs-source-dirs: test+  ghc-options: -Wall++benchmark poly-gauge+  build-depends:+    base,+    gauge,+    poly,+    semirings,+    vector+  type: exitcode-stdio-1.0+  main-is: Bench.hs+  other-modules:+    DenseBench+    SparseBench+  default-language: Haskell2010+  hs-source-dirs: bench   ghc-options: -Wall

@@ -7,20 +7,29 @@ -- Dense polynomials and a 'Num'-based interface. -- -{-# LANGUAGE PatternSynonyms     #-}+{-# LANGUAGE FlexibleInstances          #-}+{-# LANGUAGE GeneralizedNewtypeDeriving #-}+{-# LANGUAGE PatternSynonyms            #-}  module Data.Poly   ( Poly   , VPoly   , UPoly   , unPoly+  , leading   -- * Num interface   , toPoly-  , constant+  , monomial+  , scale   , pattern X   , eval   , deriv   , integral+  -- * Fractional coefficients+  , PolyOverFractional(..)   ) where -import Data.Poly.Uni.Dense hiding (quotRem)+import Data.Poly.Internal.Dense+import Data.Poly.Internal.Dense.Fractional ()+import Data.Poly.Internal.Dense.GcdDomain ()+import Data.Poly.Internal.PolyOverFractional

@@ -0,0 +1,364 @@+-- |+-- Module:      Data.Poly.Internal.Dense+-- Copyright:   (c) 2019 Andrew Lelechenko+-- Licence:     BSD3+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+--+-- Dense polynomials of one variable.+--++{-# LANGUAGE FlexibleInstances          #-}+{-# LANGUAGE GeneralizedNewtypeDeriving #-}+{-# LANGUAGE PatternSynonyms            #-}+{-# LANGUAGE ScopedTypeVariables        #-}+{-# LANGUAGE TypeFamilies               #-}+{-# LANGUAGE ViewPatterns               #-}++module Data.Poly.Internal.Dense+  ( Poly(..)+  , VPoly+  , UPoly+  , leading+  , dropWhileEndM+  -- * Num interface+  , toPoly+  , monomial+  , scale+  , pattern X+  , eval+  , deriv+  , integral+  -- * Semiring interface+  , toPoly'+  , monomial'+  , scale'+  , pattern X'+  , eval'+  , deriv'+  ) where++import Prelude hiding (quotRem, quot, rem, gcd, lcm, (^))+import Control.Monad+import Control.Monad.Primitive+import Control.Monad.ST+import Data.List (foldl', intersperse)+import Data.Semiring (Semiring(..))+import qualified Data.Semiring as Semiring+import qualified Data.Vector as V+import qualified Data.Vector.Generic as G+import qualified Data.Vector.Generic.Mutable as MG+import qualified Data.Vector.Unboxed as U++-- | Polynomials of one variable with coefficients from @a@,+-- backed by a 'G.Vector' @v@ (boxed, unboxed, storable, etc.).+--+-- Use pattern 'X' for construction:+--+-- >>> (X + 1) + (X - 1) :: VPoly Integer+-- 2 * X + 0+-- >>> (X + 1) * (X - 1) :: UPoly Int+-- 1 * X^2 + 0 * X + (-1)+--+-- Polynomials are stored normalized, without leading+-- zero coefficients, so 0 * 'X' + 1 equals to 1.+--+-- 'Ord' instance does not make much sense mathematically,+-- it is defined only for the sake of 'Data.Set.Set', 'Data.Map.Map', etc.+--+newtype Poly v a = Poly+  { unPoly :: v a+  -- ^ Convert 'Poly' to a vector of coefficients+  -- (first element corresponds to a constant term).+  }+  deriving (Eq, Ord)++instance (Show a, G.Vector v a) => Show (Poly v a) where+  showsPrec d (Poly xs)+    | G.null xs+      = showString "0"+    | G.length xs == 1+      = showsPrec d (G.head xs)+    | otherwise+      = showParen (d > 0)+      $ foldl (.) id+      $ intersperse (showString " + ")+      $ G.ifoldl (\acc i c -> showCoeff i c : acc) [] xs+    where+      showCoeff 0 c = showsPrec 7 c+      showCoeff 1 c = showsPrec 7 c . showString " * X"+      showCoeff i c = showsPrec 7 c . showString " * X^" . showsPrec 7 i++-- | Polynomials backed by boxed vectors.+type VPoly = Poly V.Vector++-- | Polynomials backed by unboxed vectors.+type UPoly = Poly U.Vector++-- | Make 'Poly' from a list of coefficients+-- (first element corresponds to a constant term).+--+-- >>> :set -XOverloadedLists+-- >>> toPoly [1,2,3] :: VPoly Integer+-- 3 * X^2 + 2 * X + 1+-- >>> toPoly [0,0,0] :: UPoly Int+-- 0+toPoly :: (Eq a, Num a, G.Vector v a) => v a -> Poly v a+toPoly = Poly . dropWhileEnd (== 0)++toPoly' :: (Eq a, Semiring a, G.Vector v a) => v a -> Poly v a+toPoly' = Poly . dropWhileEnd (== zero)++-- | Return a leading power and coefficient of a non-zero polynomial.+--+-- >>> leading ((2 * X + 1) * (2 * X^2 - 1) :: UPoly Int)+-- Just (3,4)+-- >>> leading (0 :: UPoly Int)+-- Nothing+leading :: G.Vector v a => Poly v a -> Maybe (Word, a)+leading (Poly v)+  | G.null v  = Nothing+  | otherwise = Just (fromIntegral (G.length v - 1), G.last v)++-- | Note that 'abs' = 'id' and 'signum' = 'const' 1.+instance (Eq a, Num a, G.Vector v a) => Num (Poly v a) where+  Poly xs + Poly ys = toPoly $ plusPoly (+) xs ys+  Poly xs - Poly ys = toPoly $ minusPoly negate (-) xs ys+  negate (Poly xs) = Poly $ G.map negate xs+  abs = id+  signum = const 1+  fromInteger n = case fromInteger n of+    0 -> Poly $ G.empty+    m -> Poly $ G.singleton m+  Poly xs * Poly ys = toPoly $ convolution 0 (+) (*) xs ys+  {-# INLINE (+) #-}+  {-# INLINE (-) #-}+  {-# INLINE negate #-}+  {-# INLINE fromInteger #-}+  {-# INLINE (*) #-}++instance (Eq a, Semiring a, G.Vector v a) => Semiring (Poly v a) where+  zero = Poly G.empty+  one+    | (one :: a) == zero = zero+    | otherwise = Poly $ G.singleton one+  plus (Poly xs) (Poly ys) = toPoly' $ plusPoly plus xs ys+  times (Poly xs) (Poly ys) = toPoly' $ convolution zero plus times xs ys+  {-# INLINE zero #-}+  {-# INLINE one #-}+  {-# INLINE plus #-}+  {-# INLINE times #-}++instance (Eq a, Semiring.Ring a, G.Vector v a) => Semiring.Ring (Poly v a) where+  negate (Poly xs) = Poly $ G.map Semiring.negate xs++dropWhileEnd+  :: G.Vector v a+  => (a -> Bool)+  -> v a+  -> v a+dropWhileEnd p xs = G.basicUnsafeSlice 0 (go (G.basicLength xs)) xs+  where+    go 0 = 0+    go n = if p (G.unsafeIndex xs (n - 1)) then go (n - 1) else n+{-# INLINE dropWhileEnd #-}++dropWhileEndM+  :: (PrimMonad m, G.Vector v a)+  => (a -> Bool)+  -> G.Mutable v (PrimState m) a+  -> m (G.Mutable v (PrimState m) a)+dropWhileEndM p xs = go (MG.basicLength xs)+  where+    go 0 = pure $ MG.basicUnsafeSlice 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)+{-# INLINE dropWhileEndM #-}++plusPoly+  :: G.Vector v a+  => (a -> a -> a)+  -> v a+  -> v a+  -> v a+plusPoly add xs ys = runST $ do+  let lenXs = G.basicLength xs+      lenYs = G.basicLength ys+      lenMn = lenXs `min` lenYs+      lenMx = lenXs `max` lenYs++  zs <- MG.basicUnsafeNew 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))++  G.unsafeFreeze zs+{-# INLINE plusPoly #-}++minusPoly+  :: G.Vector v a+  => (a -> a)+  -> (a -> a -> a)+  -> v a+  -> v a+  -> v a+minusPoly neg sub xs ys = runST $ do+  let lenXs = G.basicLength xs+      lenYs = G.basicLength ys+      lenMn = lenXs `min` lenYs+      lenMx = lenXs `max` lenYs++  zs <- MG.basicUnsafeNew lenMx+  forM_ [0 .. lenMn - 1] $ \i ->+    MG.unsafeWrite zs i (sub (G.unsafeIndex xs i) (G.unsafeIndex ys i))++  if lenXs < lenYs+    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)++  G.unsafeFreeze zs+{-# INLINE minusPoly #-}++convolution+  :: G.Vector v a+  => a+  -> (a -> a -> a)+  -> (a -> a -> a)+  -> v a+  -> 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+{-# INLINE convolution #-}++-- | Create a monomial from a power and a coefficient.+monomial :: (Eq a, Num a, G.Vector v a) => Word -> a -> Poly v a+monomial _ 0 = Poly G.empty+monomial p c = Poly $ G.generate (fromIntegral p + 1) (\i -> if i == fromIntegral p then c else 0)++monomial' :: (Eq a, Semiring a, G.Vector v a) => Word -> a -> Poly v a+monomial' p c+  | c == zero = Poly G.empty+  | otherwise = Poly $ G.generate (fromIntegral p + 1) (\i -> if i == fromIntegral p then c else zero)++scaleInternal+  :: (Eq a, G.Vector v a)+  => a+  -> (a -> a -> a)+  -> Word+  -> a+  -> Poly v a+  -> v a+scaleInternal zer mul yp yc (Poly xs) = runST $ do+  let lenXs = G.basicLength xs+  zs <- MG.basicUnsafeNew (fromIntegral yp + lenXs)+  forM_ [0 .. fromIntegral yp - 1] $ \k ->+    MG.unsafeWrite zs k zer+  forM_ [0 .. lenXs - 1] $ \k ->+    MG.unsafeWrite zs (fromIntegral yp + k) (mul yc $ G.unsafeIndex xs k)+  G.unsafeFreeze zs++-- | Multiply a polynomial by a monomial, expressed as a power and a coefficient.+--+-- >>> scale 2 3 (X^2 + 1) :: UPoly Int+-- 3 * X^4 + 0 * X^3 + 3 * X^2 + 0 * X + 0+scale :: (Eq a, Num a, G.Vector v a) => Word -> a -> Poly v a -> Poly v a+scale yp yc xs = toPoly $ scaleInternal 0 (*) yp yc xs++scale' :: (Eq a, Semiring a, G.Vector v a) => Word -> a -> Poly v a -> Poly v a+scale' yp yc xs = toPoly' $ scaleInternal zero times yp yc xs++data StrictPair a b = !a :*: !b++infixr 1 :*:++fst' :: StrictPair a b -> a+fst' (a :*: _) = a++-- | Evaluate at a given point.+--+-- >>> eval (X^2 + 1 :: UPoly Int) 3+-- 10+-- >>> eval (X^2 + 1 :: VPoly (UPoly Int)) (X + 1)+-- 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+{-# 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+{-# INLINE eval' #-}++-- | Take a derivative.+--+-- >>> deriv (X^3 + 3 * X) :: UPoly Int+-- 3 * X^2 + 0 * X + 3+deriv :: (Eq a, Num a, G.Vector v a) => Poly v a -> Poly v a+deriv (Poly xs)+  | G.null xs = Poly G.empty+  | otherwise = toPoly $ G.imap (\i x -> fromIntegral (i + 1) * x) $ G.tail xs+{-# INLINE deriv #-}++deriv' :: (Eq a, Semiring a, G.Vector v a) => Poly v a -> Poly v a+deriv' (Poly xs)+  | G.null xs = Poly G.empty+  | otherwise = toPoly' $ G.imap (\i x -> fromNatural (fromIntegral (i + 1)) `times` x) $ G.tail xs+{-# INLINE deriv' #-}++-- | Compute an indefinite integral of a polynomial,+-- setting constant term to zero.+--+-- >>> integral (3 * X^2 + 3) :: UPoly Double+-- 1.0 * X^3 + 0.0 * X^2 + 3.0 * X + 0.0+integral :: (Eq a, Fractional a, G.Vector v a) => Poly v a -> Poly v a+integral (Poly xs)+  | G.null xs = Poly G.empty+  | otherwise = toPoly $ runST $ do+    zs <- MG.basicUnsafeNew (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+{-# INLINE integral #-}++-- | Create an identity polynomial.+pattern X :: (Eq a, Num a, G.Vector v a, Eq (v a)) => Poly v a+pattern X <- ((==) var -> True)+  where X = var++var :: forall a v. (Eq a, Num a, G.Vector v a, Eq (v a)) => Poly v a+var+  | (1 :: a) == 0 = Poly G.empty+  | otherwise     = Poly $ G.fromList [0, 1]+{-# INLINE var #-}++-- | Create an identity polynomial.+pattern X' :: (Eq a, Semiring a, G.Vector v a, Eq (v a)) => Poly v a+pattern X' <- ((==) var' -> True)+  where X' = var'++var' :: forall a v. (Eq a, Semiring a, G.Vector v a, Eq (v a)) => Poly v a+var'+  | (one :: a) == zero = Poly G.empty+  | otherwise          = Poly $ G.fromList [zero, one]+{-# INLINE var' #-}

@@ -0,0 +1,129 @@+-- |+-- 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 FlexibleInstances          #-}+{-# LANGUAGE GeneralizedNewtypeDeriving #-}+{-# LANGUAGE PatternSynonyms            #-}+{-# LANGUAGE ScopedTypeVariables        #-}+{-# LANGUAGE TypeFamilies               #-}+{-# LANGUAGE ViewPatterns               #-}++{-# OPTIONS_GHC -fno-warn-orphans #-}++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 #-}

@@ -0,0 +1,173 @@+-- |+-- Module:      Data.Poly.Internal.Dense.GcdDomain+-- Copyright:   (c) 2019 Andrew Lelechenko+-- Licence:     BSD3+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+--+-- GcdDomain for GcdDomain underlying+--++{-# LANGUAGE FlexibleInstances          #-}+{-# LANGUAGE GeneralizedNewtypeDeriving #-}+{-# LANGUAGE PatternSynonyms            #-}+{-# LANGUAGE ScopedTypeVariables        #-}+{-# LANGUAGE TypeFamilies               #-}+{-# LANGUAGE ViewPatterns               #-}++{-# OPTIONS_GHC -fno-warn-orphans #-}++module Data.Poly.Internal.Dense.GcdDomain+  () where++import Prelude hiding (gcd, lcm, (^))+import Control.Exception+import Control.Monad+import Control.Monad.Primitive+import Control.Monad.ST+import Data.Euclidean+import Data.Semiring (Semiring(..), isZero)+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++-- | Consider using 'Data.Poly.Semiring.PolyOverFractional' 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+  divide (Poly xs) (Poly ys) =+    toPoly' <$> quotient xs ys++  gcd (Poly xs) (Poly ys)+    | G.null xs = Poly ys+    | G.null ys = Poly xs+    | otherwise = toPoly' $ gcdNonEmpty xs ys+  {-# INLINE gcd #-}++gcdNonEmpty+  :: (Eq a, Semiring.Ring a, GcdDomain a, G.Vector v a)+  => v a+  -> v a+  -> v a+gcdNonEmpty xs ys = runST $ do+    let x = G.foldl1' gcd xs+        y = G.foldl1' gcd ys+        xy = x `gcd` y+    xs' <- G.thaw xs+    ys' <- G.thaw ys+    zs' <- gcdM xs' ys'++    let lenZs = MG.basicLength zs'+        go acc 0 = pure acc+        go acc n = do+          t <- MG.unsafeRead zs' (n - 1)+          go (acc `gcd` t) (n - 1)+    a <- MG.unsafeRead zs' (lenZs - 1)+    z <- go a (lenZs - 1)++    let err = error "gcdNonEmpty: violated internal invariant"+    forM_ [0 .. lenZs - 1] $ \i ->+      MG.unsafeModify+        zs'+        (\c -> maybe err (`times` xy) (c `divide` z))+        i++    G.unsafeFreeze zs'++gcdM+  :: (PrimMonad m, Eq a, Semiring.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)+gcdM xs ys+  | MG.null xs = pure ys+  | MG.null ys = pure xs+  | otherwise = do+  let lenXs = MG.basicLength xs+      lenYs = MG.basicLength ys+  xLast <- MG.unsafeRead xs (lenXs - 1)+  yLast <- MG.unsafeRead ys (lenYs - 1)+  let z = xLast `lcm` yLast+      zx = case z `divide` xLast of+        Nothing -> error "gcdM: highest coefficient is 0"+        Just t  -> t+      zy = case z `divide` yLast of+        Nothing -> error "gcdM: highest coefficient is 0"+        Just t  -> t++  if lenXs <= lenYs then do+    forM_ [0 .. lenXs - 1] $ \i -> do+      x <- MG.unsafeRead xs i+      MG.unsafeModify+        ys+        (\y -> (y `times` zy) `plus` Semiring.negate (x `times` zx))+        (i + lenYs - lenXs)+    forM_ [0 .. lenYs - lenXs - 1] $+      MG.unsafeModify ys (`times` zy)+    ys' <- dropWhileEndM isZero ys+    gcdM xs ys'+  else do+    forM_ [0 .. lenYs - 1] $ \i -> do+      y <- MG.unsafeRead ys i+      MG.unsafeModify+        xs+        (\x -> (x `times` zx) `plus` Semiring.negate (y `times` zy))+        (i + lenXs - lenYs)+    forM_ [0 .. lenXs - lenYs - 1] $+      MG.unsafeModify xs (`times` zx)+    xs' <- dropWhileEndM isZero xs+    gcdM xs' ys+{-# INLINE gcdM #-}++isZeroM+  :: (Eq a, Semiring a, PrimMonad m, G.Vector v a)+  => G.Mutable v (PrimState m) a+  -> m Bool+isZeroM xs = go (MG.basicLength xs)+  where+    go 0 = pure True+    go n = do+      x <- MG.unsafeRead xs (n - 1)+      if x == zero then go (n - 1) else pure False+{-# INLINE isZeroM #-}++quotient+  :: (Eq a, Eq (v a), Semiring.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+  | 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++    let go i+          | i < 0 = do+            b <- isZeroM rs+            if b+              then Just <$> G.unsafeFreeze qs+              else pure Nothing+          | otherwise = do+            r <- MG.unsafeRead rs (lenYs - 1 + i)+            case r `divide` G.unsafeLast ys of+              Nothing -> pure Nothing+              Just q -> do+                MG.unsafeWrite qs i q+                forM_ [0 .. lenYs - 1] $ \k -> do+                  MG.unsafeModify+                    rs+                    (\c -> c `plus` (Semiring.negate $ q `times` G.unsafeIndex ys k))+                    (i + k)+                go (i - 1)++    go (lenQs - 1)+{-# INLINE quotient #-}

@@ -0,0 +1,46 @@+-- |+-- 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 FlexibleInstances          #-}+{-# LANGUAGE GeneralizedNewtypeDeriving #-}+{-# LANGUAGE UndecidableInstances       #-}++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 #-}

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

@@ -0,0 +1,69 @@+-- |+-- 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 FlexibleContexts           #-}+{-# LANGUAGE FlexibleInstances          #-}+{-# LANGUAGE GeneralizedNewtypeDeriving #-}+{-# LANGUAGE PatternSynonyms            #-}+{-# LANGUAGE ScopedTypeVariables        #-}+{-# LANGUAGE TypeFamilies               #-}+{-# LANGUAGE UndecidableInstances       #-}+{-# LANGUAGE ViewPatterns               #-}++{-# OPTIONS_GHC -fno-warn-orphans #-}++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 #-}

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

@@ -14,19 +14,26 @@   , VPoly   , UPoly   , unPoly+  , leading   -- * Semiring interface   , toPoly-  , constant+  , monomial+  , scale   , pattern X   , eval   , deriv+  -- * Fractional coefficients+  , PolyOverFractional(..)   ) where  import Data.Semiring (Semiring) import qualified Data.Vector.Generic as G -import Data.Poly.Uni.Dense (Poly(..), VPoly, UPoly)-import qualified Data.Poly.Uni.Dense as Dense+import Data.Poly.Internal.Dense (Poly(..), VPoly, UPoly, leading)+import qualified Data.Poly.Internal.Dense as Dense+import Data.Poly.Internal.Dense.Fractional ()+import Data.Poly.Internal.Dense.GcdDomain ()+import Data.Poly.Internal.PolyOverFractional  -- | Make 'Poly' from a vector of coefficients -- (first element corresponds to a constant term).@@ -39,9 +46,16 @@ toPoly :: (Eq a, Semiring a, G.Vector v a) => v a -> Poly v a toPoly = Dense.toPoly' --- | Create a polynomial from a constant term.-constant :: (Eq a, Semiring a, G.Vector v a) => a -> Poly v a-constant = Dense.constant'+-- | Create a monomial from a power and a coefficient.+monomial :: (Eq a, Semiring a, G.Vector v a) => Word -> a -> Poly v a+monomial = Dense.monomial'++-- | Multiply a polynomial by a monomial, expressed as a power and a coefficient.+--+-- >>> scale 2 3 (X^2 + 1) :: UPoly Int+-- 3 * X^4 + 0 * X^3 + 3 * X^2 + 0 * X + 0+scale :: (Eq a, Semiring a, G.Vector v a) => Word -> a -> Poly v a -> Poly v a+scale = Dense.scale'  -- | Create an identity polynomial. pattern X :: (Eq a, Semiring a, G.Vector v a, Eq (v a)) => Poly v a

@@ -0,0 +1,30 @@+-- |+-- Module:      Data.Poly.Sparse+-- Copyright:   (c) 2019 Andrew Lelechenko+-- Licence:     BSD3+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+--+-- Sparse polynomials with 'Num' instance.+--++{-# LANGUAGE PatternSynonyms     #-}++module Data.Poly.Sparse+  ( Poly+  , VPoly+  , UPoly+  , unPoly+  , leading+  -- * Num interface+  , toPoly+  , monomial+  , scale+  , pattern X+  , eval+  , deriv+  , integral+  ) where++import Data.Poly.Internal.Sparse+import Data.Poly.Internal.Sparse.Fractional ()+import Data.Poly.Internal.Sparse.GcdDomain ()

@@ -0,0 +1,76 @@+-- |+-- Module:      Data.Poly.Sparse.Semiring+-- Copyright:   (c) 2019 Andrew Lelechenko+-- Licence:     BSD3+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+--+-- Sparse polynomials with 'Semiring' instance.+--++{-# LANGUAGE FlexibleContexts    #-}+{-# LANGUAGE PatternSynonyms     #-}++module Data.Poly.Sparse.Semiring+  ( Poly+  , VPoly+  , UPoly+  , unPoly+  , leading+  -- * Semiring interface+  , toPoly+  , monomial+  , scale+  , pattern X+  , eval+  , deriv+  ) where++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+import Data.Poly.Internal.Sparse.Fractional ()+import Data.Poly.Internal.Sparse.GcdDomain ()++-- | Make 'Poly' from a list of (power, coefficient) pairs.+-- (first element corresponds to a constant term).+--+-- >>> :set -XOverloadedLists+-- >>> toPoly [(0,1),(1,2),(2,3)] :: VPoly Integer+-- 3 * X^2 + 2 * X + 1+-- >>> S.toPoly [(0,0),(1,0),(2,0)] :: UPoly Int+-- 0+toPoly :: (Eq a, Semiring a, G.Vector v (Word, a)) => v (Word, a) -> Poly v a+toPoly = Sparse.toPoly'++-- | Create a monomial from a power and a coefficient.+monomial :: (Eq a, Semiring a, G.Vector v (Word, a)) => Word -> a -> Poly v a+monomial = Sparse.monomial'++-- | Multiply a polynomial by a monomial, expressed as a power and a coefficient.+--+-- >>> scale 2 3 (X^2 + 1) :: UPoly Int+-- 3 * X^4 + 3 * X^2+scale :: (Eq a, Semiring a, G.Vector v (Word, a)) => Word -> a -> Poly v a -> Poly v a+scale = Sparse.scale'++-- | Create an identity polynomial.+pattern X :: (Eq a, Semiring a, G.Vector v (Word, a), Eq (v (Word, a))) => Poly v a+pattern X = Sparse.X'++-- | Evaluate at a given point.+--+-- >>> eval (X^2 + 1 :: UPoly Int) 3+-- 10+-- >>> eval (X^2 + 1 :: VPoly (UPoly Int)) (X + 1)+-- 1 * X^2 + 2 * X + 2+eval :: (Semiring a, G.Vector v (Word, a)) => Poly v a -> a -> a+eval = Sparse.eval'++-- | Take a derivative.+--+-- >>> deriv (X^3 + 3 * X) :: UPoly Int+-- 3 * X^2 + 3+deriv :: (Eq a, Semiring a, G.Vector v (Word, a)) => Poly v a -> Poly v a+deriv = Sparse.deriv'

@@ -1,356 +0,0 @@--- |--- Module:      Data.Poly.Uni.Dense--- Copyright:   (c) 2019 Andrew Lelechenko--- Licence:     BSD3--- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>------ Dense polynomials of one variable.-----{-# LANGUAGE PatternSynonyms     #-}-{-# LANGUAGE ScopedTypeVariables #-}-{-# LANGUAGE ViewPatterns        #-}--module Data.Poly.Uni.Dense-  ( Poly-  , VPoly-  , UPoly-  , unPoly-  -- * Num interface-  , toPoly-  , constant-  , pattern X-  , eval-  , deriv-  , integral-  , quotRem-  -- * Semiring interface-  , toPoly'-  , constant'-  , pattern X'-  , eval'-  , deriv'-  ) where--import Prelude hiding (quotRem)-import Control.Exception-import Control.Monad-import Control.Monad.Primitive-import Control.Monad.ST-import Data.List (foldl', intersperse)-import Data.Semigroup (stimes)-import Data.Semiring (Semiring(..), Add(..))-import qualified Data.Semiring as Semiring-import qualified Data.Vector as V-import qualified Data.Vector.Generic as G-import qualified Data.Vector.Generic.Mutable as MG-import qualified Data.Vector.Unboxed as U---- | Polynomials of one variable with coefficients from @a@,--- backed by a 'G.Vector' @v@ (boxed, unboxed, storable, etc.).------ Use pattern 'X' for construction:------ >>> (X + 1) + (X - 1) :: VPoly Integer--- 2 * X + 0--- >>> (X + 1) * (X - 1) :: UPoly Int--- 1 * X^2 + 0 * X + (-1)------ Polynomials are stored normalized, without leading--- zero coefficients, so 0 * 'X' + 1 equals to 1.------ 'Ord' instance does not make much sense mathematically,--- it is defined only for the sake of 'Data.Set.Set', 'Data.Map.Map', etc.----newtype Poly v a = Poly-  { unPoly :: v a-  -- ^ Convert 'Poly' to a vector of coefficients-  -- (first element corresponds to a constant term).-  }-  deriving (Eq, Ord)--instance (Show a, G.Vector v a) => Show (Poly v a) where-  showsPrec d (Poly xs)-    | G.null xs-      = showString "0"-    | G.length xs == 1-      = showsPrec d (G.head xs)-    | otherwise-      = showParen (d > 0)-      $ foldl (.) id-      $ intersperse (showString " + ")-      $ G.ifoldl (\acc i c -> showCoeff i c : acc) [] xs-    where-      showCoeff 0 c = showsPrec 7 c-      showCoeff 1 c = showsPrec 7 c . showString " * X"-      showCoeff i c = showsPrec 7 c . showString " * X^" . showsPrec 7 i---- | Polynomials backed by boxed vectors.-type VPoly = Poly V.Vector---- | Polynomials backed by unboxed vectors.-type UPoly = Poly U.Vector---- | Make 'Poly' from a list of coefficients--- (first element corresponds to a constant term).------ >>> :set -XOverloadedLists--- >>> toPoly [1,2,3] :: VPoly Integer--- 3 * X^2 + 2 * X + 1--- >>> toPoly [0,0,0] :: UPoly Int--- 0-toPoly :: (Eq a, Num a, G.Vector v a) => v a -> Poly v a-toPoly = Poly . dropWhileEnd (== 0)--toPoly' :: (Eq a, Semiring a, G.Vector v a) => v a -> Poly v a-toPoly' = Poly . dropWhileEnd (== zero)--instance (Eq a, Num a, G.Vector v a) => Num (Poly v a) where-  Poly xs + Poly ys = toPoly $ plusPoly (+) xs ys-  Poly xs - Poly ys = toPoly $ minusPoly negate (-) xs ys-  negate (Poly xs) = Poly $ G.map negate xs-  abs = id-  signum = const 1-  fromInteger n = case fromInteger n of-    0 -> Poly $ G.empty-    m -> Poly $ G.singleton m-  Poly xs * Poly ys = toPoly $ convolution 0 (+) (*) xs ys--instance (Eq a, Semiring a, G.Vector v a) => Semiring (Poly v a) where-  zero = Poly G.empty-  one-    | (one :: a) == zero = zero-    | otherwise = Poly $ G.singleton one-  plus (Poly xs) (Poly ys) = toPoly' $ plusPoly plus xs ys-  times (Poly xs) (Poly ys) = toPoly' $ convolution zero plus times xs ys--instance (Eq a, Semiring.Ring a, G.Vector v a) => Semiring.Ring (Poly v a) where-  negate (Poly xs) = Poly $ G.map Semiring.negate xs--dropWhileEnd-  :: G.Vector v a-  => (a -> Bool)-  -> v a-  -> v a-dropWhileEnd p xs = G.basicUnsafeSlice 0 (go (G.basicLength xs)) xs-  where-    go 0 = 0-    go n = if p (G.unsafeIndex xs (n - 1)) then go (n - 1) else n--plusPoly-  :: G.Vector v a-  => (a -> a -> a)-  -> v a-  -> v a-  -> v a-plusPoly add xs ys = runST $ do-  zs <- MG.new (G.basicLength xs `max` G.basicLength ys)-  plusPolyM add xs ys zs-  G.unsafeFreeze zs--plusPolyM-  :: (PrimMonad m, G.Vector v a)-  => (a -> a -> a)-  -> v a-  -> v a-  -> G.Mutable v (PrimState m) a-  -> m ()-plusPolyM add xs ys zs = do-  let lenXs = G.basicLength xs-      lenYs = G.basicLength ys-  case lenXs `compare` lenYs of-    LT -> do-      forM_ [0 .. lenXs - 1] $ \i ->-        MG.unsafeWrite zs i (add (G.unsafeIndex xs i) (G.unsafeIndex ys i))-      G.unsafeCopy-        (MG.basicUnsafeSlice lenXs (lenYs - lenXs) zs)-        (G.basicUnsafeSlice  lenXs (lenYs - lenXs) ys)-    EQ -> do-      forM_ [0 .. lenXs - 1] $ \i ->-        MG.unsafeWrite zs i (add (G.unsafeIndex xs i) (G.unsafeIndex ys i))-    GT -> do-      forM_ [0 .. lenYs - 1] $ \i ->-        MG.unsafeWrite zs i (add (G.unsafeIndex xs i) (G.unsafeIndex ys i))-      G.unsafeCopy-        (MG.basicUnsafeSlice lenYs (lenXs - lenYs) zs)-        (G.basicUnsafeSlice  lenYs (lenXs - lenYs) xs)--minusPoly-  :: G.Vector v a-  => (a -> a)-  -> (a -> a -> a)-  -> v a-  -> v a-  -> v a-minusPoly neg sub xs ys = runST $ do-  zs <- MG.new (G.basicLength xs `max` G.basicLength ys)-  minusPolyM neg sub xs ys zs-  G.unsafeFreeze zs--minusPolyM-  :: (PrimMonad m, G.Vector v a)-  => (a -> a)-  -> (a -> a -> a)-  -> v a-  -> v a-  -> G.Mutable v (PrimState m) a-  -> m ()-minusPolyM neg sub xs ys zs = do-  let lenXs = G.basicLength xs-      lenYs = G.basicLength ys-  case lenXs `compare` lenYs of-    LT -> do-      forM_ [0 .. lenXs - 1] $ \i ->-        MG.unsafeWrite zs i (sub (G.unsafeIndex xs i) (G.unsafeIndex ys i))-      forM_ [lenXs .. lenYs - 1] $ \i ->-        MG.unsafeWrite zs i (neg (G.unsafeIndex ys i))-    EQ -> do-      forM_ [0 .. lenXs - 1] $ \i ->-        MG.unsafeWrite zs i (sub (G.unsafeIndex xs i) (G.unsafeIndex ys i))-    GT -> do-      forM_ [0 .. lenYs - 1] $ \i ->-        MG.unsafeWrite zs i (sub (G.unsafeIndex xs i) (G.unsafeIndex ys i))-      G.unsafeCopy-        (MG.basicUnsafeSlice lenYs (lenXs - lenYs) zs)-        (G.basicUnsafeSlice  lenYs (lenXs - lenYs) xs)--convolution-  :: G.Vector v a-  => a-  -> (a -> a -> a)-  -> (a -> a -> a)-  -> v a-  -> v a-  -> v a-convolution zer add mul xs ys-  | G.null xs || G.null ys = G.empty-  | otherwise = runST $ do-    zs <- MG.new 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-  where-    lenXs = G.basicLength xs-    lenYs = G.basicLength ys-    lenZs = lenXs + lenYs - 1---- | This is just a proof of concept,--- which should be replaced by a proper 'Euclidean' interface.-quotRem-  :: (Integral a, G.Vector v a)-  => Poly v a-  -> Poly v a-  -> (Poly v a, Poly v a)-quotRem (Poly xs) (Poly ys) = (toPoly qs, toPoly rs)-  where-    (qs, rs) = quotRem' xs ys--quotRem'-  :: (Integral a, G.Vector v a)-  => v a-  -> v a-  -> (v a, v a)-quotRem' 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.new lenQs-    rs <- MG.new lenXs-    G.unsafeCopy rs xs-    forM_ [lenQs - 1, lenQs - 2 .. 0] $ \i -> do-      let j = lenXs - 1 + i - (lenQs - 1)-      r <- MG.unsafeRead rs j-      let q = r `quot` G.unsafeLast ys-      MG.unsafeWrite qs i q-      forM_ [0 .. lenYs - 1] $ \k -> do-        MG.unsafeModify rs (\c -> c - q * G.unsafeIndex ys k) (j + k - lenYs + 1)-    (,) <$> G.unsafeFreeze qs <*> G.unsafeFreeze rs----- | Create a polynomial from a constant term.-constant :: (Eq a, Num a, G.Vector v a) => a -> Poly v a-constant 0 = Poly G.empty-constant c = Poly $ G.singleton c--constant' :: (Eq a, Semiring a, G.Vector v a) => a -> Poly v a-constant' c-  | c == zero = Poly G.empty-  | otherwise = Poly $ G.singleton c--data StrictPair a b = !a :*: !b--infixr 1 :*:--fst' :: StrictPair a b -> a-fst' (a :*: _) = a---- | Evaluate at a given point.------ >>> eval (X^2 + 1 :: UPoly Int) 3--- 10--- >>> eval (X^2 + 1 :: VPoly (UPoly Int)) (X + 1)--- 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--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---- | Take a derivative.------ >>> deriv (X^3 + 3 * X) :: UPoly Int--- 3 * X^2 + 0 * X + 3-deriv :: (Eq a, Num a, G.Vector v a) => Poly v a -> Poly v a-deriv (Poly xs)-  | G.null xs = Poly G.empty-  | otherwise = toPoly $ G.imap (\i x -> fromIntegral (i + 1) * x) $ G.tail xs--deriv' :: (Eq a, Semiring a, G.Vector v a) => Poly v a -> Poly v a-deriv' (Poly xs)-  | G.null xs = Poly G.empty-  | otherwise = toPoly' $ G.imap (\i x -> getAdd (stimes (i + 1) (Add x))) $ G.tail xs---- | Compute an indefinite integral of a polynomial,--- setting constant term to zero.------ >>> integral (constant 3.0 * X^2 + constant 3.0) :: UPoly Double--- 1.0 * X^3 + 0.0 * X^2 + 3.0 * X + 0.0-integral :: (Eq a, Fractional a, G.Vector v a) => Poly v a -> Poly v a-integral (Poly xs)-  | G.null xs = Poly G.empty-  | otherwise = toPoly $ runST $ do-    zs <- MG.new (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---- | Create an identity polynomial.-pattern X :: (Eq a, Num a, G.Vector v a, Eq (v a)) => Poly v a-pattern X <- ((==) var -> True)-  where X = var--var :: forall a v. (Eq a, Num a, G.Vector v a, Eq (v a)) => Poly v a-var-  | (1 :: a) == 0 = Poly G.empty-  | otherwise     = Poly $ G.fromList [0, 1]---- | Create an identity polynomial.-pattern X' :: (Eq a, Semiring a, G.Vector v a, Eq (v a)) => Poly v a-pattern X' <- ((==) var' -> True)-  where X' = var'--var' :: forall a v. (Eq a, Semiring a, G.Vector v a, Eq (v a)) => Poly v a-var'-  | (one :: a) == zero = Poly G.empty-  | otherwise          = Poly $ G.fromList [zero, one]

@@ -0,0 +1,173 @@+{-# LANGUAGE FlexibleInstances          #-}+{-# LANGUAGE GeneralizedNewtypeDeriving #-}+{-# LANGUAGE ScopedTypeVariables        #-}++{-# OPTIONS_GHC -fno-warn-orphans #-}++module Dense+  ( testSuite+  ) where++import Prelude hiding (quotRem)+import Data.Euclidean+import Data.Int+import Data.Poly+import qualified Data.Poly.Semiring as S+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)+import Test.QuickCheck.Classes++instance (Eq a, Semiring a, Arbitrary a, G.Vector v a) => Arbitrary (Poly v a) where+  arbitrary = S.toPoly . G.fromList <$> arbitrary+  shrink = fmap (S.toPoly . G.fromList) . shrink . G.toList . unPoly++instance (Eq a, Semiring a, Arbitrary a, G.Vector v a) => Arbitrary (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++newtype ShortPoly a = ShortPoly { unShortPoly :: a }+  deriving (Eq, Show, Semiring, GcdDomain, Euclidean)++instance (Eq a, Semiring a, Arbitrary a, G.Vector v a) => Arbitrary (ShortPoly (Poly v a)) where+  arbitrary = ShortPoly . S.toPoly . G.fromList . (\xs -> take (length xs `mod` 10) xs) <$> arbitrary+  shrink = fmap (ShortPoly . S.toPoly . G.fromList) . shrink . G.toList . unPoly . unShortPoly++testSuite :: TestTree+testSuite = testGroup "Dense"+    [ arithmeticTests+    , otherTests+    , semiringTests+    , evalTests+    , derivTests+    -- , euclideanTests+    ]++semiringTests :: TestTree+semiringTests+  = testGroup "Semiring"+  $ map (uncurry testProperty)+  $ concatMap lawsProperties+  [ semiringLaws (Proxy :: Proxy (Poly U.Vector ()))+  ,     ringLaws (Proxy :: Proxy (Poly U.Vector ()))+  , semiringLaws (Proxy :: Proxy (Poly U.Vector Int8))+  ,     ringLaws (Proxy :: Proxy (Poly U.Vector Int8))+  , semiringLaws (Proxy :: Proxy (Poly V.Vector Integer))+  ,     ringLaws (Proxy :: Proxy (Poly V.Vector Integer))+  ]++-- 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)))+--   ]++arithmeticTests :: TestTree+arithmeticTests = testGroup "Arithmetic"+  [ testProperty "addition matches reference" $+    \(xs :: [Int]) ys -> toPoly (V.fromList (addRef xs ys)) ===+      toPoly (V.fromList xs) + toPoly (V.fromList ys)+  , testProperty "subtraction matches reference" $+    \(xs :: [Int]) ys -> toPoly (V.fromList (subRef xs ys)) ===+      toPoly (V.fromList xs) - toPoly (V.fromList ys)+  , testProperty "multiplication matches reference" $+    \(xs :: [Int]) ys -> toPoly (V.fromList (mulRef xs ys)) ===+      toPoly (V.fromList xs) * toPoly (V.fromList ys)+  ]++addRef :: Num a => [a] -> [a] -> [a]+addRef [] ys = ys+addRef xs [] = xs+addRef (x : xs) (y : ys) = (x + y) : addRef xs ys++subRef :: Num a => [a] -> [a] -> [a]+subRef [] ys = map negate ys+subRef xs [] = xs+subRef (x : xs) (y : ys) = (x - y) : subRef xs ys++mulRef :: Num a => [a] -> [a] -> [a]+mulRef xs ys+  = foldl addRef []+  $ zipWith (\x zs -> map (* x) zs) xs+  $ iterate (0 :) ys++otherTests :: TestTree+otherTests = testGroup "Other"+  [ testProperty "leading p 0 == Nothing" $+    \p -> leading (monomial p 0 :: UPoly Int) === Nothing+  , testProperty "leading . monomial = id" $+    \p c -> c /= 0 ==> leading (monomial p c :: UPoly Int) === Just (p, c)+  , testProperty "monomial matches reference" $+    \p (c :: Int) -> monomial p c === toPoly (V.fromList (monomialRef p c))+  , testProperty "scale matches multiplication by monomial" $+    \p c (xs :: UPoly Int) -> scale p c xs === monomial p c * xs+  ]++monomialRef :: Num a => Word -> a -> [a]+monomialRef p c = replicate (fromIntegral p) 0 ++ [c]++evalTests :: TestTree+evalTests = testGroup "eval" $ concat+  [ evalTestGroup (Proxy :: Proxy (Poly U.Vector Int8))+  , evalTestGroup (Proxy :: Proxy (Poly V.Vector Integer))+  ]++evalTestGroup+  :: forall v a.+     (Eq a, Num a, Semiring a, Arbitrary a, Show a, Eq (v a), Show (v a), G.Vector v a)+  => Proxy (Poly v a)+  -> [TestTree]+evalTestGroup _ =+  [ testProperty "eval (p + q) r = eval p r + eval q r" $+    \p q r -> e (p + q) r === e p r + e q r+  , testProperty "eval (p * q) r = eval p r * eval q r" $+    \p q r -> e (p * q) r === e p r * e q r+  , testProperty "eval x p = p" $+    \p -> e X p === p+  , testProperty "eval (monomial 0 c) p = c" $+    \c p -> e (monomial 0 c) p === c++  , testProperty "eval' (p + q) r = eval' p r + eval' q r" $+    \p q r -> e' (p + q) r === e' p r + e' q r+  , testProperty "eval' (p * q) r = eval' p r * eval' q r" $+    \p q r -> e' (p * q) r === e' p r * e' q r+  , testProperty "eval' x p = p" $+    \p -> e' S.X p === p+  , testProperty "eval' (S.monomial 0 c) p = c" $+    \c p -> e' (S.monomial 0 c) p === c+  ]++  where+    e :: Poly v a -> a -> a+    e = eval+    e' :: Poly v a -> a -> a+    e' = S.eval++derivTests :: TestTree+derivTests = testGroup "deriv"+  [ testProperty "deriv = S.deriv" $+    \(p :: Poly V.Vector Integer) -> deriv p === S.deriv p+  , testProperty "deriv . integral = id" $+    \(p :: Poly V.Vector Rational) -> deriv (integral p) === p+  , testProperty "deriv c = 0" $+    \c -> deriv (monomial 0 c :: Poly V.Vector Int) === 0+  , testProperty "deriv cX = c" $+    \c -> deriv (monomial 0 c * X :: Poly V.Vector Int) === monomial 0 c+  , testProperty "deriv (p + q) = deriv p + deriv q" $+    \p q -> deriv (p + q) === (deriv p + deriv q :: Poly V.Vector Int)+  , 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)+  -- The following property takes too long for a regular test-suite+  -- , testProperty "deriv (eval p q) = deriv q * eval (deriv p) q" $+  --   \(p :: Poly V.Vector Int) (q :: Poly U.Vector Int) ->+  --     deriv (eval (toPoly $ fmap (monomial 0) $ unPoly p) q) ===+  --       deriv q * eval (toPoly $ fmap (monomial 0) $ unPoly $ deriv p) q+  ]

@@ -1,129 +1,12 @@-{-# LANGUAGE ScopedTypeVariables #-}--{-# OPTIONS_GHC -fno-warn-orphans #-}- module Main where -import Prelude hiding (quotRem)-import Data.Int-import Data.Poly-import qualified Data.Poly.Semiring as S-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-import Test.QuickCheck.Classes (lawsProperties, semiringLaws, ringLaws) -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+import qualified Dense as Dense+import qualified Sparse as Sparse  main :: IO () main = defaultMain $ testGroup "All"-    [ arithmeticTests-    , semiringTests-    , evalTests-    , derivTests-    , quotRemTests+    [ Dense.testSuite+    , Sparse.testSuite     ]--semiringTests :: TestTree-semiringTests-  = testGroup "Semiring"-  $ map (uncurry testProperty)-  $ concatMap lawsProperties-  [ semiringLaws (Proxy :: Proxy (Poly U.Vector ()))-  ,     ringLaws (Proxy :: Proxy (Poly U.Vector ()))-  , semiringLaws (Proxy :: Proxy (Poly U.Vector Int8))-  ,     ringLaws (Proxy :: Proxy (Poly U.Vector Int8))-  , semiringLaws (Proxy :: Proxy (Poly V.Vector Integer))-  ,     ringLaws (Proxy :: Proxy (Poly V.Vector Integer))-  ]--arithmeticTests :: TestTree-arithmeticTests = testGroup "Arithmetic"-  [ testProperty "addition matches reference" $-    \(xs :: [Int]) ys -> toPoly (V.fromList (addRef xs ys)) ===-      toPoly (V.fromList xs) + toPoly (V.fromList ys)-  , testProperty "subtraction matches reference" $-    \(xs :: [Int]) ys -> toPoly (V.fromList (subRef xs ys)) ===-      toPoly (V.fromList xs) - toPoly (V.fromList ys)-  ]--addRef :: Num a => [a] -> [a] -> [a]-addRef [] ys = ys-addRef xs [] = xs-addRef (x : xs) (y : ys) = (x + y) : addRef xs ys--subRef :: Num a => [a] -> [a] -> [a]-subRef [] ys = map negate ys-subRef xs [] = xs-subRef (x : xs) (y : ys) = (x - y) : subRef xs ys--evalTests :: TestTree-evalTests = testGroup "eval" $ concat-  [ evalTestGroup (Proxy :: Proxy (Poly U.Vector Int8))-  , evalTestGroup (Proxy :: Proxy (Poly V.Vector Integer))-  ]--evalTestGroup-  :: forall v a.-     (Eq a, Num a, Semiring a, Arbitrary a, Show a, Eq (v a), Show (v a), G.Vector v a)-  => Proxy (Poly v a)-  -> [TestTree]-evalTestGroup _ =-  [ testProperty "eval (p + q) r = eval p r + eval q r" $-    \p q r -> e (p + q) r === e p r + e q r-  , testProperty "eval (p * q) r = eval p r * eval q r" $-    \p q r -> e (p * q) r === e p r * e q r-  , testProperty "eval x p = p" $-    \p -> e X p === p-  , testProperty "eval (constant c) p = c" $-    \c p -> e (constant c) p === c--  , testProperty "eval' (p + q) r = eval' p r + eval' q r" $-    \p q r -> e' (p + q) r === e' p r + e' q r-  , testProperty "eval' (p * q) r = eval' p r * eval' q r" $-    \p q r -> e' (p * q) r === e' p r * e' q r-  , testProperty "eval' x p = p" $-    \p -> e' S.X p === p-  , testProperty "eval' (S.constant c) p = c" $-    \c p -> e' (S.constant c) p === c-  ]--  where-    e :: Poly v a -> a -> a-    e = eval-    e' :: Poly v a -> a -> a-    e' = S.eval--derivTests :: TestTree-derivTests = testGroup "deriv"-  [ testProperty "deriv = S.deriv" $-    \(p :: Poly V.Vector Integer) -> deriv p === S.deriv p-  , testProperty "deriv . integral = id" $-    \(p :: Poly V.Vector Rational) -> deriv (integral p) === p-  , testProperty "deriv c = 0" $-    \c -> deriv (constant c :: Poly V.Vector Int) === 0-  , testProperty "deriv cX = c" $-    \c -> deriv (constant c * X :: Poly V.Vector Int) === constant c-  , testProperty "deriv (p + q) = deriv p + deriv q" $-    \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)-  -- The following property takes too long for a regular test-suite-  -- , testProperty "deriv (eval p q) = deriv q * eval (deriv p) q" $-  --   \(p :: Poly V.Vector Int) (q :: Poly U.Vector Int) ->-  --     deriv (eval (toPoly $ fmap constant $ unPoly p) q) ===-  --       deriv q * eval (toPoly $ fmap constant $ unPoly $ deriv p) q-  ]--quotRemTests :: TestTree-quotRemTests = testGroup "quotRem" []-  -- [ testProperty "(q, r) = x `quotRem` y ==> q * y + r == x" $-  --   \(x :: Poly U.Vector Int) y -> let (q, r) = x `quotRem` y in-  --     y === 0 .||. q * y + r === x-  -- ]

@@ -0,0 +1,171 @@+{-# LANGUAGE FlexibleContexts           #-}+{-# LANGUAGE FlexibleInstances          #-}+{-# LANGUAGE GeneralizedNewtypeDeriving #-}+{-# LANGUAGE ScopedTypeVariables        #-}+{-# LANGUAGE UndecidableInstances       #-}++{-# OPTIONS_GHC -fno-warn-orphans #-}++module Sparse+  ( testSuite+  ) where++import Prelude hiding (quotRem)+import Data.Euclidean+import Data.Function+import Data.Int+import Data.List+import Data.Poly.Sparse+import qualified Data.Poly.Sparse.Semiring as S+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)+import Test.QuickCheck.Classes++instance (Eq a, Semiring a, Arbitrary a, G.Vector v (Word, a)) => Arbitrary (Poly v a) where+  arbitrary = S.toPoly . G.fromList <$> arbitrary+  shrink = fmap (S.toPoly . G.fromList) . shrink . G.toList . unPoly++newtype ShortPoly a = ShortPoly { unShortPoly :: a }+  deriving (Eq, Show, Semiring, GcdDomain, Euclidean)++instance (Eq a, Semiring a, Arbitrary a, G.Vector v (Word, a)) => Arbitrary (ShortPoly (Poly v a)) where+  arbitrary = ShortPoly . S.toPoly . G.fromList . (\xs -> take (length xs `mod` 5) xs) <$> arbitrary+  shrink = fmap (ShortPoly . S.toPoly . G.fromList) . shrink . G.toList . unPoly . unShortPoly++testSuite :: TestTree+testSuite = testGroup "Sparse"+    [ arithmeticTests+    , otherTests+    , semiringTests+    , evalTests+    , derivTests+    ]++semiringTests :: TestTree+semiringTests+  = testGroup "Semiring"+  $ map (uncurry testProperty)+  $ concatMap lawsProperties+  [ semiringLaws (Proxy :: Proxy (Poly U.Vector ()))+  ,     ringLaws (Proxy :: Proxy (Poly U.Vector ()))+  , semiringLaws (Proxy :: Proxy (Poly U.Vector Int8))+  ,     ringLaws (Proxy :: Proxy (Poly U.Vector Int8))+  , semiringLaws (Proxy :: Proxy (Poly V.Vector Integer))+  ,     ringLaws (Proxy :: Proxy (Poly V.Vector Integer))+  ]++arithmeticTests :: TestTree+arithmeticTests = testGroup "Arithmetic"+  [ testProperty "addition matches reference" $+    \(xs :: [(Word, Int)]) ys -> toPoly (V.fromList (addRef xs ys)) ===+      toPoly (V.fromList xs) + toPoly (V.fromList ys)+  , 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" $+    \(xs :: [(Word, Int)]) ys -> toPoly (V.fromList (mulRef xs ys)) ===+      toPoly (V.fromList xs) * toPoly (V.fromList ys)+  ]++addRef :: Num a => [(Word, a)] -> [(Word, a)] -> [(Word, a)]+addRef [] ys = ys+addRef xs [] = xs+addRef xs@((xp, xc) : xs') ys@((yp, yc) : ys') =+  case xp `compare` yp of+    LT -> (xp, xc) : addRef xs' ys+    EQ -> (xp, xc + yc) : addRef xs' ys'+    GT -> (yp, yc) : addRef xs ys'++subRef :: Num a => [(Word, a)] -> [(Word, a)] -> [(Word, a)]+subRef [] ys = map (fmap negate) ys+subRef xs [] = xs+subRef xs@((xp, xc) : xs') ys@((yp, yc) : ys') =+  case xp `compare` yp of+    LT -> (xp, xc) : subRef xs' ys+    EQ -> (xp, xc - yc) : subRef xs' ys'+    GT -> (yp, negate yc) : subRef xs ys'++mulRef :: Num a => [(Word, a)] -> [(Word, a)] -> [(Word, a)]+mulRef xs ys+  = map (\ws -> (fst (head ws), sum (map snd ws)))+  $ groupBy ((==) `on` fst)+  $ sortOn fst+  $ [ (xp + yp, xc * yc) | (xp, xc) <- xs, (yp, yc) <- ys ]++otherTests :: TestTree+otherTests = testGroup "Other"+  [ testProperty "leading p 0 == Nothing" $+    \p -> leading (monomial p 0 :: UPoly Int) === Nothing+  , testProperty "leading . monomial = id" $+    \p c -> c /= 0 ==> leading (monomial p c :: UPoly Int) === Just (p, c)+  , testProperty "monomial matches reference" $+    \p (c :: Int) -> monomial p c === toPoly (V.fromList (monomialRef p c))+  , testProperty "scale matches multiplication by monomial" $+    \p c (xs :: UPoly Int) -> scale p c xs === monomial p c * xs+  ]++monomialRef :: Num a => Word -> a -> [(Word, a)]+monomialRef p c = [(p, c)]++evalTests :: TestTree+evalTests = testGroup "eval" $ concat+  [ evalTestGroup (Proxy :: Proxy (Poly U.Vector Int8))+  , evalTestGroup (Proxy :: Proxy (Poly V.Vector Integer))+  ]++evalTestGroup+  :: forall v a.+     (Eq a, Num a, Semiring a, Arbitrary a, Show a, Eq (v (Word, a)), Show (v (Word, a)), G.Vector v (Word, a))+  => Proxy (Poly v a)+  -> [TestTree]+evalTestGroup _ =+  [ testProperty "eval (p + q) r = eval p r + eval q r" $+    \p q r -> e (p + q) r === e p r + e q r+  , testProperty "eval (p * q) r = eval p r * eval q r" $+    \p q r -> e (p * q) r === e p r * e q r+  , testProperty "eval x p = p" $+    \p -> e X p === p+  , testProperty "eval (monomial 0 c) p = c" $+    \c p -> e (monomial 0 c) p === c++  , testProperty "eval' (p + q) r = eval' p r + eval' q r" $+    \p q r -> e' (p + q) r === e' p r + e' q r+  , testProperty "eval' (p * q) r = eval' p r * eval' q r" $+    \p q r -> e' (p * q) r === e' p r * e' q r+  , testProperty "eval' x p = p" $+    \p -> e' S.X p === p+  , testProperty "eval' (S.monomial 0 c) p = c" $+    \c p -> e' (S.monomial 0 c) p === c+  ]++  where+    e :: Poly v a -> a -> a+    e = eval+    e' :: Poly v a -> a -> a+    e' = S.eval++derivTests :: TestTree+derivTests = testGroup "deriv"+  [ testProperty "deriv = S.deriv" $+    \(p :: Poly V.Vector Integer) -> deriv p === S.deriv p+  , testProperty "deriv . integral = id" $+    \(p :: Poly V.Vector Rational) -> deriv (integral p) === p+  , testProperty "deriv c = 0" $+    \c -> deriv (monomial 0 c :: Poly V.Vector Int) === 0+  , testProperty "deriv cX = c" $+    \c -> deriv (monomial 0 c * X :: Poly V.Vector Int) === monomial 0 c+  , testProperty "deriv (p + q) = deriv p + deriv q" $+    \p q -> deriv (p + q) === (deriv p + deriv q :: Poly V.Vector Int)+  , 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)+  -- The following property takes too long for a regular test-suite+  -- , testProperty "deriv (eval p q) = deriv q * eval (deriv p) q" $+  --   \(p :: Poly V.Vector Int) (q :: Poly U.Vector Int) ->+  --     deriv (eval (toPoly $ fmap (fmap $ monomial 0) $ unPoly p) q) ===+  --       deriv q * eval (toPoly $ fmap (fmap $ monomial 0) $ unPoly $ deriv p) q+  ]