diff --git a/README.md b/README.md
--- a/README.md
+++ b/README.md
@@ -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.
diff --git a/bench/Bench.hs b/bench/Bench.hs
new file mode 100644
--- /dev/null
+++ b/bench/Bench.hs
@@ -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
+  ]
diff --git a/bench/DenseBench.hs b/bench/DenseBench.hs
new file mode 100644
--- /dev/null
+++ b/bench/DenseBench.hs
@@ -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
diff --git a/bench/SparseBench.hs b/bench/SparseBench.hs
new file mode 100644
--- /dev/null
+++ b/bench/SparseBench.hs
@@ -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
+
diff --git a/changelog.md b/changelog.md
--- a/changelog.md
+++ b/changelog.md
@@ -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.(-)`.
diff --git a/poly.cabal b/poly.cabal
--- a/poly.cabal
+++ b/poly.cabal
@@ -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
diff --git a/src/Data/Poly.hs b/src/Data/Poly.hs
--- a/src/Data/Poly.hs
+++ b/src/Data/Poly.hs
@@ -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
diff --git a/src/Data/Poly/Internal/Dense.hs b/src/Data/Poly/Internal/Dense.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Poly/Internal/Dense.hs
@@ -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' #-}
diff --git a/src/Data/Poly/Internal/Dense/Fractional.hs b/src/Data/Poly/Internal/Dense/Fractional.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Poly/Internal/Dense/Fractional.hs
@@ -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 #-}
diff --git a/src/Data/Poly/Internal/Dense/GcdDomain.hs b/src/Data/Poly/Internal/Dense/GcdDomain.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Poly/Internal/Dense/GcdDomain.hs
@@ -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 #-}
diff --git a/src/Data/Poly/Internal/PolyOverFractional.hs b/src/Data/Poly/Internal/PolyOverFractional.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Poly/Internal/PolyOverFractional.hs
@@ -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 #-}
diff --git a/src/Data/Poly/Internal/Sparse.hs b/src/Data/Poly/Internal/Sparse.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Poly/Internal/Sparse.hs
@@ -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' #-}
diff --git a/src/Data/Poly/Internal/Sparse/Fractional.hs b/src/Data/Poly/Internal/Sparse/Fractional.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Poly/Internal/Sparse/Fractional.hs
@@ -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 #-}
diff --git a/src/Data/Poly/Internal/Sparse/GcdDomain.hs b/src/Data/Poly/Internal/Sparse/GcdDomain.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Poly/Internal/Sparse/GcdDomain.hs
@@ -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"
diff --git a/src/Data/Poly/Semiring.hs b/src/Data/Poly/Semiring.hs
--- a/src/Data/Poly/Semiring.hs
+++ b/src/Data/Poly/Semiring.hs
@@ -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
diff --git a/src/Data/Poly/Sparse.hs b/src/Data/Poly/Sparse.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Poly/Sparse.hs
@@ -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 ()
diff --git a/src/Data/Poly/Sparse/Semiring.hs b/src/Data/Poly/Sparse/Semiring.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Poly/Sparse/Semiring.hs
@@ -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'
diff --git a/src/Data/Poly/Uni/Dense.hs b/src/Data/Poly/Uni/Dense.hs
deleted file mode 100644
--- a/src/Data/Poly/Uni/Dense.hs
+++ /dev/null
@@ -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]
diff --git a/test/Dense.hs b/test/Dense.hs
new file mode 100644
--- /dev/null
+++ b/test/Dense.hs
@@ -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
+  ]
diff --git a/test/Main.hs b/test/Main.hs
--- a/test/Main.hs
+++ b/test/Main.hs
@@ -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
-  -- ]
diff --git a/test/Sparse.hs b/test/Sparse.hs
new file mode 100644
--- /dev/null
+++ b/test/Sparse.hs
@@ -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
+  ]