diff --git a/poly.cabal b/poly.cabal
index 80f49f2..a0b061c 100644
--- a/poly.cabal
+++ b/poly.cabal
@@ -26,23 +26,30 @@ library
   exposed-modules:
     Data.Poly
     Data.Poly.Laurent
-    Data.Poly.Orthogonal
     Data.Poly.Semiring
+    Data.Poly.Orthogonal
+
     Data.Poly.Sparse
     Data.Poly.Sparse.Laurent
-    Data.Poly.Sparse.Multi
     Data.Poly.Sparse.Semiring
+
+    Data.Poly.Multi
+    Data.Poly.Multi.Laurent
+    Data.Poly.Multi.Semiring
   other-modules:
     Data.Poly.Internal.Convert
+
     Data.Poly.Internal.Dense
     Data.Poly.Internal.Dense.Field
     Data.Poly.Internal.Dense.DFT
     Data.Poly.Internal.Dense.GcdDomain
     Data.Poly.Internal.Dense.Laurent
-    Data.Poly.Internal.Sparse
-    Data.Poly.Internal.Sparse.Field
-    Data.Poly.Internal.Sparse.GcdDomain
-    Data.Poly.Internal.Sparse.Laurent
+
+    Data.Poly.Internal.Multi
+    Data.Poly.Internal.Multi.Core
+    Data.Poly.Internal.Multi.Field
+    Data.Poly.Internal.Multi.GcdDomain
+    Data.Poly.Internal.Multi.Laurent
   build-depends:
     base >= 4.9 && < 5,
     deepseq >= 1.1 && < 1.5,
diff --git a/src/Data/Poly/Internal/Convert.hs b/src/Data/Poly/Internal/Convert.hs
index 5064a73..260f47b 100644
--- a/src/Data/Poly/Internal/Convert.hs
+++ b/src/Data/Poly/Internal/Convert.hs
@@ -7,6 +7,7 @@
 -- Conversions between polynomials.
 --
 
+{-# LANGUAGE DataKinds        #-}
 {-# LANGUAGE FlexibleContexts #-}
 
 module Data.Poly.Internal.Convert
@@ -20,46 +21,67 @@ import Control.Monad.ST
 import Data.Semiring (Semiring(..))
 import qualified Data.Vector.Generic as G
 import qualified Data.Vector.Generic.Mutable as MG
+import qualified Data.Vector.Unboxed.Sized as SU
 
 import qualified Data.Poly.Internal.Dense as Dense
-import qualified Data.Poly.Internal.Sparse as Sparse
+import qualified Data.Poly.Internal.Multi as Sparse
 
 -- | Convert from dense to sparse polynomials.
 --
 -- >>> :set -XFlexibleContexts
 -- >>> denseToSparse (1 + Data.Poly.X^2) :: Data.Poly.Sparse.UPoly Int
 -- 1 * X^2 + 1
-denseToSparse :: (Eq a, Num a, G.Vector v a, G.Vector v (Word, a)) => Dense.Poly v a -> Sparse.Poly v a
+denseToSparse
+  :: (Eq a, Num a, G.Vector v a, G.Vector v (SU.Vector 1 Word, a))
+  => Dense.Poly v a
+  -> Sparse.Poly v a
 denseToSparse = denseToSparseInternal 0
 
-denseToSparse' :: (Eq a, Semiring a, G.Vector v a, G.Vector v (Word, a)) => Dense.Poly v a -> Sparse.Poly v a
+denseToSparse'
+  :: (Eq a, Semiring a, G.Vector v a, G.Vector v (SU.Vector 1 Word, a))
+  => Dense.Poly v a
+  -> Sparse.Poly v a
 denseToSparse' = denseToSparseInternal zero
 
-denseToSparseInternal :: (Eq a, G.Vector v a, G.Vector v (Word, a)) => a -> Dense.Poly v a -> Sparse.Poly v a
-denseToSparseInternal z = Sparse.Poly . G.imapMaybe (\i c -> if c == z then Nothing else Just (fromIntegral i, c)) . Dense.unPoly
+denseToSparseInternal
+  :: (Eq a, G.Vector v a, G.Vector v (SU.Vector 1 Word, a))
+  => a
+  -> Dense.Poly v a
+  -> Sparse.Poly v a
+denseToSparseInternal z = Sparse.MultiPoly . G.imapMaybe (\i c -> if c == z then Nothing else Just (fromIntegral i, c)) . Dense.unPoly
 
 -- | Convert from sparse to dense polynomials.
 --
 -- >>> :set -XFlexibleContexts
 -- >>> sparseToDense (1 + Data.Poly.Sparse.X^2) :: Data.Poly.UPoly Int
 -- 1 * X^2 + 0 * X + 1
-sparseToDense :: (Num a, G.Vector v a, G.Vector v (Word, a)) => Sparse.Poly v a -> Dense.Poly v a
+sparseToDense
+  :: (Num a, G.Vector v a, G.Vector v (SU.Vector 1 Word, a))
+  => Sparse.Poly v a
+  -> Dense.Poly v a
 sparseToDense = sparseToDenseInternal 0
 
-sparseToDense' :: (Semiring a, G.Vector v a, G.Vector v (Word, a)) => Sparse.Poly v a -> Dense.Poly v a
+sparseToDense'
+  :: (Semiring a, G.Vector v a, G.Vector v (SU.Vector 1 Word, a))
+  => Sparse.Poly v a
+  -> Dense.Poly v a
 sparseToDense' = sparseToDenseInternal zero
 
-sparseToDenseInternal :: (G.Vector v a, G.Vector v (Word, a)) => a -> Sparse.Poly v a -> Dense.Poly v a
-sparseToDenseInternal z (Sparse.Poly xs)
+sparseToDenseInternal
+  :: (G.Vector v a, G.Vector v (SU.Vector 1 Word, a))
+  => a
+  -> Sparse.Poly v a
+  -> Dense.Poly v a
+sparseToDenseInternal z (Sparse.MultiPoly xs)
   | G.null xs = Dense.Poly G.empty
   | otherwise = runST $ do
-  let len = fromIntegral (fst (G.unsafeLast xs) + 1)
+  let len = fromIntegral (SU.head (fst (G.unsafeLast xs)) + 1)
   ys <- MG.unsafeNew len
   MG.set ys z
   let go xi yi
         | xi >= G.length xs = pure ()
         | (yi', c) <- G.unsafeIndex xs xi
-        , yi == fromIntegral yi'
+        , yi == fromIntegral (SU.head yi')
         = MG.unsafeWrite ys yi c >> go (xi + 1) (yi + 1)
         | otherwise = go xi (yi + 1)
   go 0 0
diff --git a/src/Data/Poly/Sparse/Multi.hs b/src/Data/Poly/Internal/Multi.hs
similarity index 68%
rename from src/Data/Poly/Sparse/Multi.hs
rename to src/Data/Poly/Internal/Multi.hs
index 425721d..8bc283e 100644
--- a/src/Data/Poly/Sparse/Multi.hs
+++ b/src/Data/Poly/Internal/Multi.hs
@@ -1,5 +1,5 @@
 -- |
--- Module:      Data.Poly.Sparse.Multi
+-- Module:      Data.Poly.Internal.Multi
 -- Copyright:   (c) 2020 Andrew Lelechenko
 -- Licence:     BSD3
 -- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>
@@ -25,24 +25,41 @@
 {-# LANGUAGE QuantifiedConstraints      #-}
 #endif
 
-module Data.Poly.Sparse.Multi
-  ( MultiPoly
+{-# OPTIONS_GHC -fno-warn-orphans #-}
+
+module Data.Poly.Internal.Multi
+  ( MultiPoly(..)
   , VMultiPoly
   , UMultiPoly
-  , unMultiPoly
   , toMultiPoly
+  , toMultiPoly'
+  , leading
   , monomial
+  , monomial'
   , scale
+  , scale'
   , pattern X
   , pattern Y
   , pattern Z
+  , pattern X'
+  , pattern Y'
+  , pattern Z'
   , eval
+  , eval'
   , subst
+  , subst'
+  , substitute
+  , substitute'
   , deriv
+  , deriv'
   , integral
+  , integral'
+  -- * Univariate polynomials
+  , Poly
+  , VPoly
+  , UPoly
+  , unPoly
   -- * Conversions
-  , polyToMultiPoly
-  , multiPolyToPoly
   , segregate
   , unsegregate
   ) where
@@ -50,42 +67,34 @@ module Data.Poly.Sparse.Multi
 import Prelude hiding (quot, gcd)
 import Control.Arrow
 import Control.DeepSeq
-import Control.Exception
-import Data.Euclidean (GcdDomain(..), Field, quot)
+import Data.Euclidean (Field, quot)
 import Data.Finite
 import Data.Kind
 import Data.List (intersperse)
-import Data.Proxy
 import Data.Semiring (Semiring(..), Ring())
 import qualified Data.Semiring as Semiring
-import Data.Type.Equality
-import qualified Data.Vector.Generic as G
-import qualified Data.Vector.Unboxed as U
 import qualified Data.Vector as V
+import qualified Data.Vector.Generic as G
 import qualified Data.Vector.Generic.Sized as SG
+import qualified Data.Vector.Unboxed as U
 import qualified Data.Vector.Unboxed.Sized as SU
 import qualified Data.Vector.Sized as SV
 import GHC.Exts (IsList(..))
-import Unsafe.Coerce
 
-import Data.Poly.Internal.Sparse (Poly(..), VPoly, normalize, plusPoly, minusPoly, convolution, scaleInternal, derivPoly)
-import Data.Poly.Internal.Sparse.GcdDomain ()
+import Data.Poly.Internal.Multi.Core
 
 #if MIN_VERSION_base(4,10,0)
-import GHC.TypeNats (Nat, KnownNat, type (+), type (<=), SomeNat(..), natVal, someNatVal)
+import GHC.TypeNats (KnownNat, Nat, type (+), type (<=))
 #else
-import GHC.TypeLits (Nat, KnownNat, type (+), type (<=), SomeNat(..), natVal)
-import qualified GHC.TypeLits as TL
-import Data.Maybe
-
-someNatVal :: Integer -> SomeNat
-someNatVal = fromJust . TL.someNatVal
+import GHC.TypeLits (KnownNat, Nat, type (+), type (<=))
 #endif
 
 -- | Polynomials of @n@ variables with coefficients from @a@,
 -- backed by a 'G.Vector' @v@ (boxed, unboxed, storable, etc.).
 --
--- Use patterns 'X', 'Y' and 'Z' for construction:
+-- Use patterns 'Data.Poly.Multi.X',
+-- 'Data.Poly.Multi.Y' and
+-- 'Data.Poly.Multi.Z' for construction:
 --
 -- >>> :set -XDataKinds
 -- >>> (X + 1) + (Y - 1) + Z :: VMultiPoly 3 Integer
@@ -94,7 +103,7 @@ someNatVal = fromJust . TL.someNatVal
 -- 1 * X * Y + (-1) * X + 1 * Y + (-1)
 --
 -- Polynomials are stored normalized, without
--- zero coefficients, so 0 * 'X' + 1 equals to 1.
+-- zero coefficients, so 0 * 'Data.Poly.Multi.X' + 1 equals to 1.
 --
 -- 'Ord' instance does not make much sense mathematically,
 -- it is defined only for the sake of 'Data.Set.Set', 'Data.Map.Map', etc.
@@ -108,10 +117,47 @@ deriving instance Eq     (v (SU.Vector n Word, a)) => Eq     (MultiPoly v n a)
 deriving instance Ord    (v (SU.Vector n Word, a)) => Ord    (MultiPoly v n a)
 deriving instance NFData (v (SU.Vector n Word, a)) => NFData (MultiPoly v n a)
 
+-- | Polynomials backed by boxed vectors.
+type VMultiPoly n = MultiPoly V.Vector n
+
+-- | Polynomials backed by unboxed vectors.
+type UMultiPoly n = MultiPoly U.Vector n
+
+-- | Polynomials of one variable with coefficients from @a@,
+-- backed by a 'G.Vector' @v@ (boxed, unboxed, storable, etc.).
+--
+-- Use pattern 'Data.Poly.Multi.X' for construction:
+--
+-- >>> (X + 1) + (X - 1) :: VPoly Integer
+-- 2 * X
+-- >>> (X + 1) * (X - 1) :: UPoly Int
+-- 1 * X^2 + (-1)
+--
+-- Polynomials are stored normalized, without
+-- zero coefficients, so 0 * 'Data.Poly.Multi.X' + 1 equals to 1.
+--
+-- 'Ord' instance does not make much sense mathematically,
+-- it is defined only for the sake of 'Data.Set.Set', 'Data.Map.Map', etc.
+--
+type Poly v a = MultiPoly v 1 a
+
+-- | Polynomials backed by boxed vectors.
+type VPoly a = VMultiPoly 1 a
+
+-- | Polynomials backed by unboxed vectors.
+type UPoly a = UMultiPoly 1 a
+
+-- | Convert 'Poly' to a vector of coefficients.
+unPoly
+  :: (G.Vector v (Word, a), G.Vector v (SU.Vector 1 Word, a))
+  => Poly v a
+  -> v (Word, a)
+unPoly = G.map (first SU.head) . unMultiPoly
+
 instance (Eq a, Semiring a, G.Vector v (SU.Vector n Word, a)) => IsList (MultiPoly v n a) where
   type Item (MultiPoly v n a) = (SU.Vector n Word, a)
-  fromList = toMultiPoly . G.fromList
-  fromListN = (toMultiPoly .) . G.fromListN
+  fromList = toMultiPoly' . G.fromList
+  fromListN = (toMultiPoly' .) . G.fromListN
   toList = G.toList . unMultiPoly
 
 instance (Show a, KnownNat n, G.Vector v (SU.Vector n Word, a)) => Show (MultiPoly v n a) where
@@ -141,12 +187,6 @@ instance (Show a, KnownNat n, G.Vector v (SU.Vector n Word, a)) => Show (MultiPo
         2 -> "Z"
         k -> "X" ++ show k
 
--- | Polynomials backed by boxed vectors.
-type VMultiPoly n = MultiPoly V.Vector n
-
--- | Polynomials backed by unboxed vectors.
-type UMultiPoly n = MultiPoly U.Vector n
-
 -- | Make 'MultiPoly' from a list of (powers, coefficient) pairs.
 --
 -- >>> :set -XOverloadedLists -XDataKinds
@@ -156,10 +196,16 @@ type UMultiPoly n = MultiPoly U.Vector n
 -- >>> toMultiPoly [(fromTuple (0,0),0),(fromTuple (0,1),0),(fromTuple (1,0),0)] :: UMultiPoly 2 Int
 -- 0
 toMultiPoly
+  :: (Eq a, Num a, G.Vector v (SU.Vector n Word, a))
+  => v (SU.Vector n Word, a)
+  -> MultiPoly v n a
+toMultiPoly = MultiPoly . normalize (/= 0) (+)
+
+toMultiPoly'
   :: (Eq a, Semiring a, G.Vector v (SU.Vector n Word, a))
   => v (SU.Vector n Word, a)
   -> MultiPoly v n a
-toMultiPoly = MultiPoly . normalize (/= zero) plus
+toMultiPoly' = MultiPoly . normalize (/= zero) plus
 
 -- | Note that 'abs' = 'id' and 'signum' = 'const' 1.
 instance (Eq a, Num a, KnownNat n, G.Vector v (SU.Vector n Word, a)) => Num (MultiPoly v n a) where
@@ -199,6 +245,18 @@ instance (Eq a, Semiring a, KnownNat n, G.Vector v (SU.Vector n Word, a)) => Sem
 instance (Eq a, Ring a, KnownNat n, G.Vector v (SU.Vector n Word, a)) => Ring (MultiPoly v n a) where
   negate (MultiPoly xs) = MultiPoly $ G.map (fmap Semiring.negate) xs
 
+-- | Return a leading power and coefficient of a non-zero polynomial.
+--
+-- >>> import Data.Poly.Sparse (UPoly)
+-- >>> leading ((2 * X + 1) * (2 * X^2 - 1) :: UPoly Int)
+-- Just (3,4)
+-- >>> leading (0 :: UPoly Int)
+-- Nothing
+leading :: G.Vector v (SU.Vector 1 Word, a) => Poly v a -> Maybe (Word, a)
+leading (MultiPoly v)
+  | G.null v  = Nothing
+  | otherwise = Just $ first SU.head $ G.last v
+
 -- | Multiply a polynomial by a monomial, expressed as powers and a coefficient.
 --
 -- >>> :set -XDataKinds
@@ -206,20 +264,37 @@ instance (Eq a, Ring a, KnownNat n, G.Vector v (SU.Vector n Word, a)) => Ring (M
 -- >>> scale (fromTuple (1, 1)) 3 (X + Y) :: UMultiPoly 2 Int
 -- 3 * X^2 * Y + 3 * X * Y^2
 scale
+  :: (Eq a, Num a, KnownNat n, G.Vector v (SU.Vector n Word, a))
+  => SU.Vector n Word
+  -> a
+  -> MultiPoly v n a
+  -> MultiPoly v n a
+scale yp yc = MultiPoly . scaleInternal (/= 0) (*) yp yc . unMultiPoly
+
+scale'
   :: (Eq a, Semiring a, KnownNat n, G.Vector v (SU.Vector n Word, a))
   => SU.Vector n Word
   -> a
   -> MultiPoly v n a
   -> MultiPoly v n a
-scale yp yc = MultiPoly . scaleInternal (/= zero) times yp yc . unMultiPoly
+scale' yp yc = MultiPoly . scaleInternal (/= zero) times yp yc . unMultiPoly
 
 -- | Create a monomial from powers and a coefficient.
 monomial
-  :: (Eq a, Semiring a, G.Vector v (SU.Vector n Word, a))
+  :: (Eq a, Num a, G.Vector v (SU.Vector n Word, a))
   => SU.Vector n Word
   -> a
   -> MultiPoly v n a
 monomial p c
+  | c == 0    = MultiPoly G.empty
+  | otherwise = MultiPoly $ G.singleton (p, c)
+
+monomial'
+  :: (Eq a, Semiring a, G.Vector v (SU.Vector n Word, a))
+  => SU.Vector n Word
+  -> a
+  -> MultiPoly v n a
+monomial' p c
   | c == zero = MultiPoly G.empty
   | otherwise = MultiPoly $ G.singleton (p, c)
 
@@ -230,13 +305,21 @@ monomial p c
 -- >>> eval (X^2 + Y^2 :: UMultiPoly 2 Int) (fromTuple (3, 4) :: Data.Vector.Sized.Vector 2 Int)
 -- 25
 eval
-  :: (Semiring a, G.Vector v (SU.Vector n Word, a), G.Vector u a)
+  :: (Num a, G.Vector v (SU.Vector n Word, a), G.Vector u a)
   => MultiPoly v n a
   -> SG.Vector u n a
   -> a
-eval = substitute times
+eval = substitute (*)
 {-# INLINE eval #-}
 
+eval'
+  :: (Semiring a, G.Vector v (SU.Vector n Word, a), G.Vector u a)
+  => MultiPoly v n a
+  -> SG.Vector u n a
+  -> a
+eval' = substitute' times
+{-# INLINE eval' #-}
+
 -- | Substitute another polynomials instead of variables.
 --
 -- >>> :set -XDataKinds
@@ -244,28 +327,52 @@ eval = substitute times
 -- >>> subst (X^2 + Y^2 + Z^2 :: UMultiPoly 3 Int) (fromTuple (X + 1, Y + 1, X + Y :: UMultiPoly 2 Int))
 -- 2 * X^2 + 2 * X * Y + 2 * X + 2 * Y^2 + 2 * Y + 2
 subst
-  :: (Eq a, Semiring a, KnownNat m, G.Vector v (SU.Vector n Word, a), G.Vector w (SU.Vector m Word, a))
+  :: (Eq a, Num a, KnownNat m, G.Vector v (SU.Vector n Word, a), G.Vector w (SU.Vector m Word, a))
   => MultiPoly v n a
   -> SV.Vector n (MultiPoly w m a)
   -> MultiPoly w m a
 subst = substitute (scale 0)
 {-# INLINE subst #-}
 
+subst'
+  :: (Eq a, Semiring a, KnownNat m, G.Vector v (SU.Vector n Word, a), G.Vector w (SU.Vector m Word, a))
+  => MultiPoly v n a
+  -> SV.Vector n (MultiPoly w m a)
+  -> MultiPoly w m a
+subst' = substitute' (scale' 0)
+{-# INLINE subst' #-}
+
 substitute
+  :: forall v u n a b.
+     (G.Vector v (SU.Vector n Word, a), G.Vector u b, Num b)
+  => (a -> b -> b)
+  -> MultiPoly v n a
+  -> SG.Vector u n b
+  -> b
+substitute f (MultiPoly cs) xs = G.foldl' go 0 cs
+  where
+    go :: b -> (SU.Vector n Word, a) -> b
+    go acc (ps, c) = acc + f c (doMonom ps)
+
+    doMonom :: SU.Vector n Word -> b
+    doMonom = SU.ifoldl' (\acc i p -> acc * ((xs `SG.index` i) ^ p)) 1
+{-# INLINE substitute #-}
+
+substitute'
   :: forall v u n a b.
      (G.Vector v (SU.Vector n Word, a), G.Vector u b, Semiring b)
   => (a -> b -> b)
   -> MultiPoly v n a
   -> SG.Vector u n b
   -> b
-substitute f (MultiPoly cs) xs = G.foldl' go zero cs
+substitute' f (MultiPoly cs) xs = G.foldl' go zero cs
   where
     go :: b -> (SU.Vector n Word, a) -> b
     go acc (ps, c) = acc `plus` f c (doMonom ps)
 
     doMonom :: SU.Vector n Word -> b
     doMonom = SU.ifoldl' (\acc i p -> acc `times` ((xs `SG.index` i) Semiring.^ p)) one
-{-# INLINE substitute #-}
+{-# INLINE substitute' #-}
 
 -- | Take a derivative with respect to the /i/-th variable.
 --
@@ -275,11 +382,22 @@ substitute f (MultiPoly cs) xs = G.foldl' go zero cs
 -- >>> deriv 1 (X^3 + 3 * Y) :: UMultiPoly 2 Int
 -- 3
 deriv
-  :: (Eq a, Semiring a, G.Vector v (SU.Vector n Word, a))
+  :: (Eq a, Num a, G.Vector v (SU.Vector n Word, a))
   => Finite n
   -> MultiPoly v n a
   -> MultiPoly v n a
 deriv i (MultiPoly xs) = MultiPoly $ derivPoly
+  (/= 0)
+  (\ps -> ps SU.// [(i, ps `SU.index` i - 1)])
+  (\ps c -> fromIntegral (ps `SU.index` i) * c)
+  xs
+
+deriv'
+  :: (Eq a, Semiring a, G.Vector v (SU.Vector n Word, a))
+  => Finite n
+  -> MultiPoly v n a
+  -> MultiPoly v n a
+deriv' i (MultiPoly xs) = MultiPoly $ derivPoly
   (/= zero)
   (\ps -> ps SU.// [(i, ps `SU.index` i - 1)])
   (\ps c -> fromNatural (fromIntegral (ps `SU.index` i)) `times` c)
@@ -295,71 +413,115 @@ deriv i (MultiPoly xs) = MultiPoly $ derivPoly
 -- >>> integral 1 (3 * X^2 + 2 * Y) :: UMultiPoly 2 Double
 -- 3.0 * X^2 * Y + 1.0 * Y^2
 integral
-  :: (Field a, G.Vector v (SU.Vector n Word, a))
+  :: (Fractional a, G.Vector v (SU.Vector n Word, a))
   => Finite n
   -> MultiPoly v n a
   -> MultiPoly v n a
 integral i (MultiPoly xs)
+  = MultiPoly
+  $ G.map (\(ps, c) -> let p = ps `SU.index` i in
+    (ps SU.// [(i, p + 1)], c / fromIntegral (p + 1))) xs
+
+integral'
+  :: (Field a, G.Vector v (SU.Vector n Word, a))
+  => Finite n
+  -> MultiPoly v n a
+  -> MultiPoly v n a
+integral' i (MultiPoly xs)
   = MultiPoly
   $ G.map (\(ps, c) -> let p = ps `SU.index` i in
     (ps SU.// [(i, p + 1)], c `quot` Semiring.fromIntegral (p + 1))) xs
 
 -- | Create a polynomial equal to the first variable.
 pattern X
-  :: (Eq a, Semiring a, KnownNat n, 1 <= n, G.Vector v (SU.Vector n Word, a))
+  :: (Eq a, Num a, KnownNat n, 1 <= n, G.Vector v (SU.Vector n Word, a))
   => MultiPoly v n a
 pattern X <- (isVar 0 -> True)
   where X = var 0
 
+pattern X'
+  :: (Eq a, Semiring a, KnownNat n, 1 <= n, G.Vector v (SU.Vector n Word, a))
+  => MultiPoly v n a
+pattern X' <- (isVar' 0 -> True)
+  where X' = var' 0
+
 -- | Create a polynomial equal to the second variable.
 pattern Y
-  :: (Eq a, Semiring a, KnownNat n, 2 <= n, G.Vector v (SU.Vector n Word, a))
+  :: (Eq a, Num a, KnownNat n, 2 <= n, G.Vector v (SU.Vector n Word, a))
   => MultiPoly v n a
 pattern Y <- (isVar 1 -> True)
   where Y = var 1
 
+pattern Y'
+  :: (Eq a, Semiring a, KnownNat n, 2 <= n, G.Vector v (SU.Vector n Word, a))
+  => MultiPoly v n a
+pattern Y' <- (isVar' 1 -> True)
+  where Y' = var' 1
+
 -- | Create a polynomial equal to the third variable.
 pattern Z
-  :: (Eq a, Semiring a, KnownNat n, 3 <= n, G.Vector v (SU.Vector n Word, a))
+  :: (Eq a, Num a, KnownNat n, 3 <= n, G.Vector v (SU.Vector n Word, a))
   => MultiPoly v n a
 pattern Z <- (isVar 2 -> True)
   where Z = var 2
 
+pattern Z'
+  :: (Eq a, Semiring a, KnownNat n, 3 <= n, G.Vector v (SU.Vector n Word, a))
+  => MultiPoly v n a
+pattern Z' <- (isVar' 2 -> True)
+  where Z' = var' 2
+
 var
   :: forall v n a.
-     (Eq a, Semiring a, KnownNat n, G.Vector v (SU.Vector n Word, a))
+     (Eq a, Num a, KnownNat n, G.Vector v (SU.Vector n Word, a))
   => Finite n
   -> MultiPoly v n a
 var i
+  | (1 :: a) == 0 = MultiPoly G.empty
+  | otherwise     = MultiPoly $ G.singleton
+    (SU.generate (\j -> if i == j then 1 else 0), 1)
+{-# INLINE var #-}
+
+var'
+  :: forall v n a.
+     (Eq a, Semiring a, KnownNat n, G.Vector v (SU.Vector n Word, a))
+  => Finite n
+  -> MultiPoly v n a
+var' i
   | (one :: a) == zero = MultiPoly G.empty
   | otherwise          = MultiPoly $ G.singleton
     (SU.generate (\j -> if i == j then 1 else 0), one)
-{-# INLINE var #-}
+{-# INLINE var' #-}
 
 isVar
   :: forall v n a.
-     (Eq a, Semiring a, KnownNat n, G.Vector v (SU.Vector n Word, a))
+     (Eq a, Num a, KnownNat n, G.Vector v (SU.Vector n Word, a))
   => Finite n
   -> MultiPoly v n a
   -> Bool
 isVar i (MultiPoly xs)
+  | (1 :: a) == 0
+  = G.null xs
+  | otherwise
+  = G.length xs == 1 && G.unsafeHead xs == (SU.generate (\j -> if i == j then 1 else 0), 1)
+{-# INLINE isVar #-}
+
+isVar'
+  :: forall v n a.
+     (Eq a, Semiring a, KnownNat n, G.Vector v (SU.Vector n Word, a))
+  => Finite n
+  -> MultiPoly v n a
+  -> Bool
+isVar' i (MultiPoly xs)
   | (one :: a) == zero
   = G.null xs
   | otherwise
   = G.length xs == 1 && G.unsafeHead xs == (SU.generate (\j -> if i == j then 1 else 0), one)
-{-# INLINE isVar #-}
+{-# INLINE isVar' #-}
 
 -------------------------------------------------------------------------------
 -- GcdDomain
 
-data IsSucc n where
-  IsSucc :: KnownNat m => n :~: 1 + m -> IsSucc n
-
--- | This is unsafe when n ~ 0.
-isSucc :: forall n. KnownNat n => IsSucc n
-isSucc = case someNatVal (natVal (Proxy :: Proxy n) - 1) of
-  SomeNat (_ :: Proxy m) -> IsSucc (unsafeCoerce Refl :: n :~: 1 + m)
-
 groupOn :: (G.Vector v a, Eq b) => (a -> b) -> v a -> [v a]
 groupOn f = go
   where
@@ -372,22 +534,6 @@ groupOn f = go
           fy = f (G.unsafeHead xs)
           mk = G.findIndex ((/= fy) . f) (G.unsafeTail xs)
 
--- | Convert a sparse univariate polynomial
--- to a multivariate polynomial of one variable.
-polyToMultiPoly
-  :: (G.Vector v (Word, a), G.Vector v (SU.Vector 1 Word, a))
-  => Poly v a
-  -> MultiPoly v 1 a
-polyToMultiPoly = MultiPoly . G.map (first SU.singleton) . unPoly
-
--- | Convert a multivariate polynomial of one variable
--- to a sparse univariate polynomial.
-multiPolyToPoly
-  :: (G.Vector v (Word, a), G.Vector v (SU.Vector 1 Word, a))
-  => MultiPoly v 1 a
-  -> Poly v a
-multiPolyToPoly = Poly . G.map (first SU.head) . unMultiPoly
-
 -- | Interpret a multivariate polynomial over 1+/m/ variables
 -- as a univariate polynomial, whose coefficients are
 -- multivariate polynomials over the last /m/ variables.
@@ -396,9 +542,9 @@ segregate
   => MultiPoly v (1 + m) a
   -> VPoly (MultiPoly v m a)
 segregate
-  = Poly
+  = MultiPoly
   . G.fromList
-  . map (\vs -> (SU.head (fst (G.unsafeHead vs)), MultiPoly $ G.map (first SU.tail) vs))
+  . map (\vs -> (SU.take (fst (G.unsafeHead vs)), MultiPoly $ G.map (first SU.tail) vs))
   . groupOn (SU.head . fst)
   . unMultiPoly
 
@@ -413,49 +559,6 @@ unsegregate
   = MultiPoly
   . G.concat
   . G.toList
-  . G.map (\(v, MultiPoly vs) -> G.map (first (SU.cons v)) vs)
-  . unPoly
+  . G.map (\(v, MultiPoly vs) -> G.map (first (v SU.++)) vs)
+  . unMultiPoly
 
-#if __GLASGOW_HASKELL__ >= 806
-instance (Eq a, Ring a, GcdDomain a, KnownNat n, forall m. KnownNat m => G.Vector v (SU.Vector m Word, a), forall m. KnownNat m => Eq (v (SU.Vector m Word, a))) => GcdDomain (MultiPoly v n a) where
-#else
-instance (Eq a, Ring a, GcdDomain a, KnownNat n) => GcdDomain (VMultiPoly n a) where
-#endif
-  divide xs ys
-    | G.null (unMultiPoly ys) = throw DivideByZero
-    | G.length (unMultiPoly ys) == 1 = divideSingleton xs (G.unsafeHead (unMultiPoly ys))
-    -- Polynomials of zero variables are necessarily constants,
-    -- so they have been dealt with above.
-    | otherwise = case isSucc :: IsSucc n of
-      IsSucc Refl -> unsegregate <$> segregate xs `divide` segregate ys
-  gcd xs ys
-    | G.null (unMultiPoly xs) = ys
-    | G.null (unMultiPoly ys) = xs
-    | G.length (unMultiPoly xs) == 1 = gcdSingleton (G.unsafeHead (unMultiPoly xs)) ys
-    | G.length (unMultiPoly ys) == 1 = gcdSingleton (G.unsafeHead (unMultiPoly ys)) xs
-    -- Polynomials of zero variables are necessarily constants,
-    -- so they have been dealt with above.
-    | otherwise = case isSucc :: IsSucc n of
-      IsSucc Refl -> unsegregate $ segregate xs `gcd` segregate ys
-
-divideSingleton
-  :: (GcdDomain a, G.Vector v (SU.Vector n Word, a))
-  => MultiPoly v n a
-  -> (SU.Vector n Word, a)
-  -> Maybe (MultiPoly v n a)
-divideSingleton (MultiPoly pcs) (p, c) = MultiPoly <$> G.mapM divideMonomial pcs
-  where
-    divideMonomial (p', c')
-      | SU.and (SU.zipWith (>=) p' p)
-      , Just c'' <- c' `divide` c
-      = Just (SU.zipWith (-) p' p, c'')
-      | otherwise
-      = Nothing
-
-gcdSingleton
-  :: (Eq a, GcdDomain a, G.Vector v (SU.Vector n Word, a))
-  => (SU.Vector n Word, a)
-  -> MultiPoly v n a
-  -> MultiPoly v n a
-gcdSingleton pc (MultiPoly pcs) = uncurry monomial $
-  G.foldl' (\(accP, accC) (p, c) -> (SU.zipWith min accP p, gcd accC c)) pc pcs
diff --git a/src/Data/Poly/Internal/Multi/Core.hs b/src/Data/Poly/Internal/Multi/Core.hs
new file mode 100644
index 0000000..8cd61d0
--- /dev/null
+++ b/src/Data/Poly/Internal/Multi/Core.hs
@@ -0,0 +1,315 @@
+-- |
+-- Module:      Data.Poly.Internal.Multi.Core
+-- Copyright:   (c) 2019 Andrew Lelechenko
+-- Licence:     BSD3
+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>
+--
+-- Sparse polynomials of one variable.
+--
+
+{-# LANGUAGE DataKinds                  #-}
+{-# LANGUAGE FlexibleContexts           #-}
+{-# LANGUAGE GeneralizedNewtypeDeriving #-}
+{-# LANGUAGE PatternSynonyms            #-}
+{-# LANGUAGE ScopedTypeVariables        #-}
+{-# LANGUAGE StandaloneDeriving         #-}
+{-# LANGUAGE TypeFamilies               #-}
+{-# LANGUAGE UndecidableInstances       #-}
+{-# LANGUAGE ViewPatterns               #-}
+
+module Data.Poly.Internal.Multi.Core
+  ( normalize
+  , plusPoly
+  , minusPoly
+  , convolution
+  , scaleInternal
+  , derivPoly
+  ) where
+
+import Control.Monad
+import Control.Monad.Primitive
+import Control.Monad.ST
+import Data.Bits
+import Data.Ord
+import qualified Data.Vector.Algorithms.Tim as Tim
+import qualified Data.Vector.Generic as G
+import qualified Data.Vector.Generic.Mutable as MG
+import qualified Data.Vector.Unboxed as U
+
+normalize
+  :: (G.Vector v (t, a), Ord t)
+  => (a -> Bool)
+  -> (a -> a -> a)
+  -> v (t, a)
+  -> v (t, a)
+normalize p add vs
+  | G.null vs = vs
+  | otherwise = runST $ do
+    ws <- G.thaw vs
+    l' <- normalizeM p add ws
+    G.unsafeFreeze $ MG.unsafeSlice 0 l' ws
+
+normalizeM
+  :: (PrimMonad m, G.Vector v (t, a), Ord t)
+  => (a -> Bool)
+  -> (a -> a -> a)
+  -> G.Mutable v (PrimState m) (t, a)
+  -> m Int
+normalizeM p add ws = do
+    let l = MG.length ws
+    let go i j acc@(accP, accC)
+          | j >= l =
+            if p accC
+              then do
+                MG.write ws i acc
+                pure $ i + 1
+              else pure i
+          | otherwise = do
+            v@(vp, vc) <- MG.unsafeRead ws j
+            if vp == accP
+              then go i (j + 1) (accP, accC `add` vc)
+              else if p accC
+                then do
+                  MG.write ws i acc
+                  go (i + 1) (j + 1) v
+                else go i (j + 1) v
+    Tim.sortBy (comparing fst) ws
+    wsHead <- MG.unsafeRead ws 0
+    go 0 1 wsHead
+
+plusPoly
+  :: (G.Vector v (t, a), Ord t)
+  => (a -> Bool)
+  -> (a -> a -> a)
+  -> v (t, a)
+  -> v (t, a)
+  -> v (t, a)
+plusPoly p add xs ys = runST $ do
+  zs <- MG.unsafeNew (G.length xs + G.length ys)
+  lenZs <- plusPolyM p add xs ys zs
+  G.unsafeFreeze $ MG.unsafeSlice 0 lenZs zs
+{-# INLINABLE plusPoly #-}
+
+plusPolyM
+  :: (PrimMonad m, G.Vector v (t, a), Ord t)
+  => (a -> Bool)
+  -> (a -> a -> a)
+  -> v (t, a)
+  -> v (t, a)
+  -> G.Mutable v (PrimState m) (t, a)
+  -> m Int
+plusPolyM p add xs ys zs = go 0 0 0
+  where
+    lenXs = G.length xs
+    lenYs = G.length ys
+
+    go ix iy iz
+      | ix == lenXs, iy == lenYs = pure iz
+      | ix == lenXs = do
+        G.unsafeCopy
+          (MG.unsafeSlice iz (lenYs - iy) zs)
+          (G.unsafeSlice iy (lenYs - iy) ys)
+        pure $ iz + lenYs - iy
+      | iy == lenYs = do
+        G.unsafeCopy
+          (MG.unsafeSlice iz (lenXs - ix) zs)
+          (G.unsafeSlice ix (lenXs - ix) xs)
+        pure $ iz + lenXs - ix
+      | (xp, xc) <- G.unsafeIndex xs ix
+      , (yp, yc) <- G.unsafeIndex ys iy
+      = case xp `compare` yp of
+        LT -> do
+          MG.unsafeWrite zs iz (xp, xc)
+          go (ix + 1) iy (iz + 1)
+        EQ -> do
+          let zc = xc `add` yc
+          if p zc then do
+            MG.unsafeWrite zs iz (xp, zc)
+            go (ix + 1) (iy + 1) (iz + 1)
+          else
+            go (ix + 1) (iy + 1) iz
+        GT -> do
+          MG.unsafeWrite zs iz (yp, yc)
+          go ix (iy + 1) (iz + 1)
+{-# INLINABLE plusPolyM #-}
+
+minusPoly
+  :: (G.Vector v (t, a), Ord t)
+  => (a -> Bool)
+  -> (a -> a)
+  -> (a -> a -> a)
+  -> v (t, a)
+  -> v (t, a)
+  -> v (t, a)
+minusPoly p neg sub xs ys = runST $ do
+  zs <- MG.unsafeNew (lenXs + lenYs)
+  let go ix iy iz
+        | ix == lenXs, iy == lenYs = pure iz
+        | ix == lenXs = do
+          forM_ [iy .. lenYs - 1] $ \i ->
+            MG.unsafeWrite zs (iz + i - iy)
+              (fmap neg (G.unsafeIndex ys i))
+          pure $ iz + lenYs - iy
+        | iy == lenYs = do
+          G.unsafeCopy
+            (MG.unsafeSlice iz (lenXs - ix) zs)
+            (G.unsafeSlice ix (lenXs - ix) xs)
+          pure $ iz + lenXs - ix
+        | (xp, xc) <- G.unsafeIndex xs ix
+        , (yp, yc) <- G.unsafeIndex ys iy
+        = case xp `compare` yp of
+          LT -> do
+            MG.unsafeWrite zs iz (xp, xc)
+            go (ix + 1) iy (iz + 1)
+          EQ -> do
+            let zc = xc `sub` yc
+            if p zc then do
+              MG.unsafeWrite zs iz (xp, zc)
+              go (ix + 1) (iy + 1) (iz + 1)
+            else
+              go (ix + 1) (iy + 1) iz
+          GT -> do
+            MG.unsafeWrite zs iz (yp, neg yc)
+            go ix (iy + 1) (iz + 1)
+  lenZs <- go 0 0 0
+  G.unsafeFreeze $ MG.unsafeSlice 0 lenZs zs
+  where
+    lenXs = G.length xs
+    lenYs = G.length ys
+{-# INLINABLE minusPoly #-}
+
+scaleM
+  :: (PrimMonad m, G.Vector v (t, a), Num t)
+  => (a -> Bool)
+  -> (a -> a -> a)
+  -> v (t, a)
+  -> (t, a)
+  -> G.Mutable v (PrimState m) (t, a)
+  -> m Int
+scaleM p mul xs (yp, yc) zs = go 0 0
+  where
+    lenXs = G.length xs
+
+    go ix iz
+      | ix == lenXs = pure iz
+      | (xp, xc) <- G.unsafeIndex xs ix
+      = do
+        let zc = xc `mul` yc
+        if p zc then do
+          MG.unsafeWrite zs iz (xp + yp, zc)
+          go (ix + 1) (iz + 1)
+        else
+          go (ix + 1) iz
+{-# INLINABLE scaleM #-}
+
+scaleInternal
+  :: (G.Vector v (t, a), Num t)
+  => (a -> Bool)
+  -> (a -> a -> a)
+  -> t
+  -> a
+  -> v (t, a)
+  -> v (t, a)
+scaleInternal p mul yp yc xs = runST $ do
+  zs <- MG.unsafeNew (G.length xs)
+  len <- scaleM p (flip mul) xs (yp, yc) zs
+  G.unsafeFreeze $ MG.unsafeSlice 0 len zs
+{-# INLINABLE scaleInternal #-}
+
+convolution
+  :: forall v t a.
+     (G.Vector v (t, a), Ord t, Num t)
+  => (a -> Bool)
+  -> (a -> a -> a)
+  -> (a -> a -> a)
+  -> v (t, a)
+  -> v (t, a)
+  -> v (t, a)
+convolution p add mult xs ys
+  | G.length xs >= G.length ys
+  = go mult xs ys
+  | otherwise
+  = go (flip mult) ys xs
+  where
+    go :: (a -> a -> a) -> v (t, a) -> v (t, a) -> v (t, a)
+    go mul long short = runST $ do
+      let lenLong   = G.length long
+          lenShort  = G.length short
+          lenBuffer = lenLong * lenShort
+      slices <- MG.unsafeNew lenShort
+      buffer <- MG.unsafeNew lenBuffer
+
+      forM_ [0 .. lenShort - 1] $ \iShort -> do
+        let (pShort, cShort) = G.unsafeIndex short iShort
+            from = iShort * lenLong
+            bufferSlice = MG.unsafeSlice from lenLong buffer
+        len <- scaleM p mul long (pShort, cShort) bufferSlice
+        MG.unsafeWrite slices iShort (from, len)
+
+      slices' <- G.unsafeFreeze slices
+      buffer' <- G.unsafeFreeze buffer
+      bufferNew <- MG.unsafeNew lenBuffer
+      gogo slices' buffer' bufferNew
+
+    gogo
+      :: PrimMonad m
+      => U.Vector (Int, Int)
+      -> v (t, a)
+      -> G.Mutable v (PrimState m) (t, a)
+      -> m (v (t, a))
+    gogo slices buffer bufferNew
+      | G.length slices == 0
+      = pure G.empty
+      | G.length slices == 1
+      , (from, len) <- G.unsafeIndex slices 0
+      = pure $ G.unsafeSlice from len buffer
+      | otherwise = do
+        let nSlices = G.length slices
+        slicesNew <- MG.unsafeNew ((nSlices + 1) `shiftR` 1)
+        forM_ [0 .. (nSlices - 2) `shiftR` 1] $ \i -> do
+          let (from1, len1) = G.unsafeIndex slices (2 * i)
+              (from2, len2) = G.unsafeIndex slices (2 * i + 1)
+              slice1 = G.unsafeSlice from1 len1 buffer
+              slice2 = G.unsafeSlice from2 len2 buffer
+              slice3 = MG.unsafeSlice from1 (len1 + len2) bufferNew
+          len3 <- plusPolyM p add slice1 slice2 slice3
+          MG.unsafeWrite slicesNew i (from1, len3)
+
+        when (odd nSlices) $ do
+          let (from, len) = G.unsafeIndex slices (nSlices - 1)
+              slice1 = G.unsafeSlice from len buffer
+              slice3 = MG.unsafeSlice from len bufferNew
+          G.unsafeCopy slice3 slice1
+          MG.unsafeWrite slicesNew (nSlices `shiftR` 1) (from, len)
+
+        slicesNew' <- G.unsafeFreeze slicesNew
+        buffer'    <- G.unsafeThaw   buffer
+        bufferNew' <- G.unsafeFreeze bufferNew
+        gogo slicesNew' bufferNew' buffer'
+{-# INLINABLE convolution #-}
+
+derivPoly
+  :: (G.Vector v (t, a))
+  => (a -> Bool)
+  -> (t -> t)
+  -> (t -> a -> a)
+  -> v (t, a)
+  -> v (t, a)
+derivPoly p dec mul xs
+  | G.null xs = G.empty
+  | otherwise = runST $ do
+    let lenXs = G.length xs
+    zs <- MG.unsafeNew lenXs
+    let go ix iz
+          | ix == lenXs = pure iz
+          | (xp, xc) <- G.unsafeIndex xs ix
+          = do
+            let zc = xp `mul` xc
+            if p zc then do
+              MG.unsafeWrite zs iz (dec xp, zc)
+              go (ix + 1) (iz + 1)
+            else
+              go (ix + 1) iz
+    lenZs <- go 0 0
+    G.unsafeFreeze $ MG.unsafeSlice 0 lenZs zs
+{-# INLINABLE derivPoly #-}
diff --git a/src/Data/Poly/Internal/Sparse/Field.hs b/src/Data/Poly/Internal/Multi/Field.hs
similarity index 63%
rename from src/Data/Poly/Internal/Sparse/Field.hs
rename to src/Data/Poly/Internal/Multi/Field.hs
index 2843432..bbcd02f 100644
--- a/src/Data/Poly/Internal/Sparse/Field.hs
+++ b/src/Data/Poly/Internal/Multi/Field.hs
@@ -1,13 +1,14 @@
 -- |
--- Module:      Data.Poly.Internal.Sparse.Field
+-- Module:      Data.Poly.Internal.Multi.Field
 -- Copyright:   (c) 2019 Andrew Lelechenko
 -- Licence:     BSD3
 -- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>
 --
--- GcdDomain for Field underlying
+-- Euclidean for Field underlying
 --
 
 {-# LANGUAGE ConstraintKinds            #-}
+{-# LANGUAGE DataKinds                  #-}
 {-# LANGUAGE FlexibleContexts           #-}
 {-# LANGUAGE FlexibleInstances          #-}
 {-# LANGUAGE ScopedTypeVariables        #-}
@@ -16,28 +17,29 @@
 
 {-# OPTIONS_GHC -fno-warn-orphans #-}
 
-module Data.Poly.Internal.Sparse.Field () where
+module Data.Poly.Internal.Multi.Field () where
 
 import Prelude hiding (quotRem, quot, rem, gcd)
 import Control.Arrow
 import Control.Exception
 import Data.Euclidean (Euclidean(..), Field)
-import Data.Semiring (minus, plus, times, zero)
+import Data.Semiring (Semiring(..), minus)
 import qualified Data.Vector.Generic as G
+import qualified Data.Vector.Unboxed.Sized as SU
 
-import Data.Poly.Internal.Sparse
-import Data.Poly.Internal.Sparse.GcdDomain ()
+import Data.Poly.Internal.Multi
+import Data.Poly.Internal.Multi.GcdDomain ()
 
 -- | Note that 'degree' 0 = 0.
-instance (Eq a, Field a, G.Vector v (Word, a)) => Euclidean (Poly v a) where
-  degree (Poly xs)
+instance (Eq a, Field a, G.Vector v (SU.Vector 1 Word, a)) => Euclidean (Poly v a) where
+  degree (MultiPoly xs)
     | G.null xs = 0
-    | otherwise = fromIntegral (fst (G.unsafeLast xs))
+    | otherwise = fromIntegral (SU.head (fst (G.unsafeLast xs)))
 
   quotRem = quotientRemainder
 
 quotientRemainder
-  :: (Eq a, Field a, G.Vector v (Word, a))
+  :: (Eq a, Field a, G.Vector v (SU.Vector 1 Word, a))
   => Poly v a
   -> Poly v a
   -> (Poly v a, Poly v a)
@@ -52,5 +54,5 @@ quotientRemainder ts ys = case leading ys of
           EQ -> (zs, xs')
           GT -> first (`plus` zs) $ go xs'
           where
-            zs = Poly $ G.singleton (xp `minus` yp, xc `quot` yc)
+            zs = MultiPoly $ G.singleton (SU.singleton (xp `minus` yp), xc `quot` yc)
             xs' = xs `minus` zs `times` ys
diff --git a/src/Data/Poly/Internal/Multi/GcdDomain.hs b/src/Data/Poly/Internal/Multi/GcdDomain.hs
new file mode 100644
index 0000000..0aa6ac0
--- /dev/null
+++ b/src/Data/Poly/Internal/Multi/GcdDomain.hs
@@ -0,0 +1,185 @@
+-- |
+-- Module:      Data.Poly.Internal.Multi.GcdDomain
+-- Copyright:   (c) 2019 Andrew Lelechenko
+-- Licence:     BSD3
+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>
+--
+-- GcdDomain for GcdDomain underlying
+--
+
+{-# LANGUAGE CPP                        #-}
+{-# LANGUAGE DataKinds                  #-}
+{-# LANGUAGE FlexibleContexts           #-}
+{-# LANGUAGE FlexibleInstances          #-}
+{-# LANGUAGE GADTs                      #-}
+{-# LANGUAGE ScopedTypeVariables        #-}
+{-# LANGUAGE TypeOperators              #-}
+{-# LANGUAGE TypeFamilies               #-}
+{-# LANGUAGE UndecidableInstances       #-}
+
+#if __GLASGOW_HASKELL__ >= 806
+{-# LANGUAGE QuantifiedConstraints      #-}
+#endif
+
+{-# OPTIONS_GHC -fno-warn-orphans #-}
+
+module Data.Poly.Internal.Multi.GcdDomain
+  () where
+
+import Prelude hiding (gcd, lcm, (^))
+import Control.Exception
+import Data.Euclidean
+import Data.Maybe
+import Data.Proxy
+import Data.Semiring (Semiring(..), Ring(), minus)
+import Data.Type.Equality
+import qualified Data.Vector.Generic as G
+import qualified Data.Vector.Unboxed.Sized as SU
+import Unsafe.Coerce
+
+import Data.Poly.Internal.Multi
+
+#if MIN_VERSION_base(4,10,0)
+import GHC.TypeNats (KnownNat, type (+), SomeNat(..), natVal, sameNat, someNatVal)
+#else
+import GHC.TypeLits (KnownNat, type (+), SomeNat(..), natVal, sameNat)
+import qualified GHC.TypeLits as TL
+import Data.Maybe
+
+someNatVal :: Integer -> SomeNat
+someNatVal = fromJust . TL.someNatVal
+#endif
+
+instance {-# OVERLAPPING #-} (Eq a, Ring a, GcdDomain a, G.Vector v (SU.Vector 1 Word, a)) => GcdDomain (Poly v a) where
+  divide xs ys
+    | G.null (unMultiPoly ys) = throw DivideByZero
+    | G.length (unMultiPoly ys) == 1 = divideSingleton xs (G.unsafeHead (unMultiPoly ys))
+    | otherwise = divide1 xs ys
+
+  gcd xs ys
+    | G.null (unMultiPoly xs) = ys
+    | G.null (unMultiPoly ys) = xs
+    | G.length (unMultiPoly xs) == 1 = gcdSingleton (G.unsafeHead (unMultiPoly xs)) ys
+    | G.length (unMultiPoly ys) == 1 = gcdSingleton (G.unsafeHead (unMultiPoly ys)) xs
+    | otherwise = gcd1 xs ys
+
+  lcm xs ys
+    | G.null (unMultiPoly xs) || G.null (unMultiPoly ys) = zero
+    | otherwise = (xs `divide'` gcd xs ys) `times` ys
+
+  coprime x y = isJust (one `divide` gcd x y)
+
+data IsSucc n where
+  IsSucc :: KnownNat m => n :~: 1 + m -> IsSucc n
+
+-- | This is unsafe when n ~ 0.
+isSucc :: forall n. KnownNat n => IsSucc n
+isSucc = case someNatVal (natVal (Proxy :: Proxy n) - 1) of
+  SomeNat (_ :: Proxy m) -> IsSucc (unsafeCoerce Refl :: n :~: 1 + m)
+
+#if __GLASGOW_HASKELL__ >= 806
+instance (Eq a, Ring a, GcdDomain a, KnownNat n, forall m. KnownNat m => G.Vector v (SU.Vector m Word, a), forall m. KnownNat m => Eq (v (SU.Vector m Word, a))) => GcdDomain (MultiPoly v n a) where
+#else
+instance (Eq a, Ring a, GcdDomain a, KnownNat n) => GcdDomain (VMultiPoly n a) where
+#endif
+  divide xs ys
+    | G.null (unMultiPoly ys) = throw DivideByZero
+    | G.length (unMultiPoly ys) == 1 = divideSingleton xs (G.unsafeHead (unMultiPoly ys))
+    -- Polynomials of zero variables are necessarily constants,
+    -- so they have been dealt with above.
+    | Just Refl <- sameNat (Proxy :: Proxy n) (Proxy :: Proxy 1)
+    = divide1 xs ys
+    | otherwise = case isSucc :: IsSucc n of
+      IsSucc Refl -> unsegregate <$> segregate xs `divide` segregate ys
+  gcd xs ys
+    | G.null (unMultiPoly xs) = ys
+    | G.null (unMultiPoly ys) = xs
+    | G.length (unMultiPoly xs) == 1 = gcdSingleton (G.unsafeHead (unMultiPoly xs)) ys
+    | G.length (unMultiPoly ys) == 1 = gcdSingleton (G.unsafeHead (unMultiPoly ys)) xs
+    -- Polynomials of zero variables are necessarily constants,
+    -- so they have been dealt with above.
+    | Just Refl <- sameNat (Proxy :: Proxy n) (Proxy :: Proxy 1)
+    = gcd1 xs ys
+    | otherwise = case isSucc :: IsSucc n of
+      IsSucc Refl -> unsegregate $ segregate xs `gcd` segregate ys
+
+divideSingleton
+  :: (GcdDomain a, G.Vector v (SU.Vector n Word, a))
+  => MultiPoly v n a
+  -> (SU.Vector n Word, a)
+  -> Maybe (MultiPoly v n a)
+divideSingleton (MultiPoly pcs) (p, c) = MultiPoly <$> G.mapM divideMonomial pcs
+  where
+    divideMonomial (p', c')
+      | SU.and (SU.zipWith (>=) p' p)
+      , Just c'' <- c' `divide` c
+      = Just (SU.zipWith (-) p' p, c'')
+      | otherwise
+      = Nothing
+
+gcdSingleton
+  :: (Eq a, GcdDomain a, G.Vector v (SU.Vector n Word, a))
+  => (SU.Vector n Word, a)
+  -> MultiPoly v n a
+  -> MultiPoly v n a
+gcdSingleton pc (MultiPoly pcs) = uncurry monomial' $
+  G.foldl' (\(accP, accC) (p, c) -> (SU.zipWith min accP p, gcd accC c)) pc pcs
+
+divide1
+  :: (Eq a, GcdDomain a, Ring a, G.Vector v (SU.Vector 1 Word, a))
+  => Poly v a
+  -> Poly v a
+  -> Maybe (Poly v a)
+divide1 xs ys = case leading ys of
+  Nothing -> throw DivideByZero
+  Just (yp, yc) -> case leading xs of
+    Nothing -> Just xs
+    Just (xp, xc)
+      | xp < yp -> Nothing
+      | otherwise -> do
+        zc <- divide xc yc
+        let z = MultiPoly $ G.singleton (SU.singleton (xp - yp), zc)
+        rest <- divide1 (xs `minus` z `times` ys) ys
+        pure $ rest `plus` z
+
+gcd1
+  :: (Eq a, GcdDomain a, Ring a, G.Vector v (SU.Vector 1 Word, a))
+  => Poly v a
+  -> Poly v a
+  -> Poly v a
+gcd1 x@(MultiPoly xs) y@(MultiPoly ys) =
+  times xy (divide1' z (monomial' 0 (content zs)))
+  where
+    z@(MultiPoly zs) = gcdHelper x y
+    xy = monomial' 0 (gcd (content xs) (content ys))
+    divide1' = (fromMaybe (error "gcd: violated internal invariant") .) . divide1
+
+content :: (GcdDomain a, G.Vector v (t, a)) => v (t, a) -> a
+content = G.foldl' (\acc (_, t) -> gcd acc t) zero
+
+gcdHelper
+  :: (Eq a, Ring a, GcdDomain a, G.Vector v (SU.Vector 1 Word, a))
+  => Poly v a
+  -> Poly v a
+  -> Poly v a
+gcdHelper xs ys = case (leading xs, leading ys) of
+  (Nothing, _) -> ys
+  (_, Nothing) -> xs
+  (Just (xp, xc), Just (yp, yc))
+    | yp <= xp
+    , Just xy <- xc `divide` yc
+    -> gcdHelper ys (xs `minus` ys `times` monomial' (SU.singleton (xp - yp)) xy)
+    | xp <= yp
+    , Just yx <- yc `divide` xc
+    -> gcdHelper xs (ys `minus` xs `times` monomial' (SU.singleton (yp - xp)) yx)
+    | yp <= xp
+    -> gcdHelper ys (xs `times` monomial' 0 gx `minus` ys `times` monomial' (SU.singleton (xp - yp)) gy)
+    | otherwise
+    -> gcdHelper xs (ys `times` monomial' 0 gy `minus` xs `times` monomial' (SU.singleton (yp - xp)) gx)
+    where
+      g = lcm xc yc
+      gx = divide' g xc
+      gy = divide' g yc
+
+divide' :: GcdDomain a => a -> a -> a
+divide' = (fromMaybe (error "gcd: violated internal invariant") .) . divide
diff --git a/src/Data/Poly/Internal/Sparse/Laurent.hs b/src/Data/Poly/Internal/Multi/Laurent.hs
similarity index 63%
rename from src/Data/Poly/Internal/Sparse/Laurent.hs
rename to src/Data/Poly/Internal/Multi/Laurent.hs
index 46e2faf..ca88ce1 100644
--- a/src/Data/Poly/Internal/Sparse/Laurent.hs
+++ b/src/Data/Poly/Internal/Multi/Laurent.hs
@@ -1,12 +1,14 @@
 -- |
--- Module:      Data.Poly.Internal.Sparse.Laurent
+-- Module:      Data.Poly.Internal.Multi.Laurent
 -- Copyright:   (c) 2020 Andrew Lelechenko
 -- Licence:     BSD3
 -- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>
 --
--- Sparse <https://en.wikipedia.org/wiki/Laurent_polynomial Laurent polynomials>.
+-- Multivariate sparse
+-- <https://en.wikipedia.org/wiki/Laurent_polynomial Laurent polynomials>.
 --
 
+{-# LANGUAGE DataKinds                  #-}
 {-# LANGUAGE FlexibleContexts           #-}
 {-# LANGUAGE FlexibleInstances          #-}
 {-# LANGUAGE PatternSynonyms            #-}
@@ -16,7 +18,7 @@
 {-# LANGUAGE UndecidableInstances       #-}
 {-# LANGUAGE ViewPatterns               #-}
 
-module Data.Poly.Internal.Sparse.Laurent
+module Data.Poly.Internal.Multi.Laurent
   ( Laurent
   , VLaurent
   , ULaurent
@@ -43,13 +45,15 @@ import Data.Semiring (Semiring(..), Ring())
 import qualified Data.Semiring as Semiring
 import qualified Data.Vector as V
 import qualified Data.Vector.Generic as G
+import qualified Data.Vector.Sized as SV
 import qualified Data.Vector.Unboxed as U
+import qualified Data.Vector.Unboxed.Sized as SU
 import GHC.Exts
 
-import Data.Poly.Internal.Sparse (Poly(..))
-import qualified Data.Poly.Internal.Sparse as Sparse
-import Data.Poly.Internal.Sparse.Field ()
-import Data.Poly.Internal.Sparse.GcdDomain ()
+import Data.Poly.Internal.Multi (Poly, MultiPoly(..))
+import qualified Data.Poly.Internal.Multi as Multi
+import Data.Poly.Internal.Multi.Field ()
+import Data.Poly.Internal.Multi.GcdDomain ()
 
 -- | <https://en.wikipedia.org/wiki/Laurent_polynomial Laurent polynomials>
 -- of one variable with coefficients from @a@,
@@ -70,10 +74,10 @@ import Data.Poly.Internal.Sparse.GcdDomain ()
 --
 data Laurent v a = Laurent !Int !(Poly v a)
 
-deriving instance Eq  (v (Word, a)) => Eq  (Laurent v a)
-deriving instance Ord (v (Word, a)) => Ord (Laurent v a)
+deriving instance Eq  (v (SU.Vector 1 Word, a)) => Eq  (Laurent v a)
+deriving instance Ord (v (SU.Vector 1 Word, a)) => Ord (Laurent v a)
 
-instance (Eq a, Semiring a, G.Vector v (Word, a)) => IsList (Laurent v a) where
+instance (Eq a, Semiring a, G.Vector v (SU.Vector 1 Word, a)) => IsList (Laurent v a) where
   type Item (Laurent v a) = (Int, a)
 
   fromList xs = toLaurent minPow (fromList ys)
@@ -86,8 +90,8 @@ instance (Eq a, Semiring a, G.Vector v (Word, a)) => IsList (Laurent v a) where
       minPow = minimum $ maxBound : map fst xs
       ys = map (first (fromIntegral . subtract minPow)) xs
 
-  toList (Laurent off poly) =
-    map (first ((+ off) . fromIntegral)) $ G.toList $ unPoly poly
+  toList (Laurent off (MultiPoly poly)) =
+    map (first ((+ off) . fromIntegral . SU.head)) $ G.toList poly
 
 -- | Deconstruct a 'Laurent' polynomial into an offset (largest possible)
 -- and a regular polynomial.
@@ -104,51 +108,50 @@ unLaurent :: Laurent v a -> (Int, Poly v a)
 unLaurent (Laurent off poly) = (off, poly)
 
 -- | Construct 'Laurent' polynomial from an offset and a regular polynomial.
--- One can imagine it as 'Data.Poly.Sparse.Semiring.scale', but allowing negative offsets.
+-- One can imagine it as 'Data.Poly.Multi.Semiring.scale', but allowing negative offsets.
 --
--- >>> toLaurent 2 (2 * Data.Poly.Sparse.X + 1) :: ULaurent Int
+-- >>> toLaurent 2 (2 * Data.Poly.Multi.X + 1) :: ULaurent Int
 -- 2 * X^3 + 1 * X^2
--- >>> toLaurent (-2) (2 * Data.Poly.Sparse.X + 1) :: ULaurent Int
+-- >>> toLaurent (-2) (2 * Data.Poly.Multi.X + 1) :: ULaurent Int
 -- 2 * X^-1 + 1 * X^-2
 toLaurent
-  :: (Eq a, Semiring a, G.Vector v (Word, a))
+  :: (Eq a, Semiring a, G.Vector v (SU.Vector 1 Word, a))
   => Int
   -> Poly v a
   -> Laurent v a
-toLaurent off (Poly xs)
+toLaurent off (MultiPoly xs)
   | G.null xs = Laurent 0 zero
-  | otherwise = Laurent (off + fromIntegral minPow) (Poly ys)
+  | otherwise = Laurent (off + fromIntegral (SU.head minPow)) (MultiPoly ys)
     where
       minPow = fst $ G.minimumBy (comparing fst) xs
       ys = if minPow == 0 then xs else G.map (first (subtract minPow)) xs
 {-# INLINE toLaurent #-}
 
 toLaurentNum
-  :: (Eq a, Num a, G.Vector v (Word, a))
+  :: (Eq a, Num a, G.Vector v (SU.Vector 1 Word, a))
   => Int
   -> Poly v a
   -> Laurent v a
-toLaurentNum off (Poly xs)
+toLaurentNum off (MultiPoly xs)
   | G.null xs = Laurent 0 0
-  | otherwise = Laurent (off + fromIntegral minPow) (Poly ys)
+  | otherwise = Laurent (off + fromIntegral (SU.head minPow)) (MultiPoly ys)
     where
       minPow = fst $ G.minimumBy (comparing fst) xs
       ys = if minPow == 0 then xs else G.map (first (subtract minPow)) xs
 {-# INLINE toLaurentNum #-}
 
-instance NFData (v (Word, a)) => NFData (Laurent v a) where
+instance NFData (v (SU.Vector 1 Word, a)) => NFData (Laurent v a) where
   rnf (Laurent off poly) = rnf off `seq` rnf poly
 
-instance (Show a, G.Vector v (Word, a)) => Show (Laurent v a) where
-  showsPrec d (Laurent off poly)
-    | G.null (unPoly poly)
+instance (Show a, G.Vector v (SU.Vector 1 Word, a)) => Show (Laurent v a) where
+  showsPrec d (Laurent off (MultiPoly poly))
+    | G.null poly
       = showString "0"
     | otherwise
       = showParen (d > 0)
       $ foldl (.) id
       $ intersperse (showString " + ")
-      $ G.foldl (\acc (i, c) -> showCoeff (fromIntegral i + off) c : acc) []
-      $ unPoly poly
+      $ G.foldl (\acc (i, c) -> showCoeff (fromIntegral (SU.head i) + off) c : acc) [] poly
     where
       showCoeff 0 c = showsPrec 7 c
       showCoeff 1 c = showsPrec 7 c . showString " * X"
@@ -166,20 +169,20 @@ type ULaurent = Laurent U.Vector
 -- Just (3,4)
 -- >>> leading (0 :: ULaurent Int)
 -- Nothing
-leading :: G.Vector v (Word, a) => Laurent v a -> Maybe (Int, a)
-leading (Laurent off poly) = first ((+ off) . fromIntegral) <$> Sparse.leading poly
+leading :: G.Vector v (SU.Vector 1 Word, a) => Laurent v a -> Maybe (Int, a)
+leading (Laurent off poly) = first ((+ off) . fromIntegral) <$> Multi.leading poly
 
 -- | Note that 'abs' = 'id' and 'signum' = 'const' 1.
-instance (Eq a, Num a, G.Vector v (Word, a)) => Num (Laurent v a) where
+instance (Eq a, Num a, G.Vector v (SU.Vector 1 Word, a)) => Num (Laurent v a) where
   Laurent off1 poly1 * Laurent off2 poly2 = toLaurentNum (off1 + off2) (poly1 * poly2)
   Laurent off1 poly1 + Laurent off2 poly2 = case off1 `compare` off2 of
-    LT -> toLaurentNum off1 (poly1 + Sparse.scale (fromIntegral $ off2 - off1) 1 poly2)
+    LT -> toLaurentNum off1 (poly1 + Multi.scale (fromIntegral $ off2 - off1) 1 poly2)
     EQ -> toLaurentNum off1 (poly1 + poly2)
-    GT -> toLaurentNum off2 (Sparse.scale (fromIntegral $ off1 - off2) 1 poly1 + poly2)
+    GT -> toLaurentNum off2 (Multi.scale (fromIntegral $ off1 - off2) 1 poly1 + poly2)
   Laurent off1 poly1 - Laurent off2 poly2 = case off1 `compare` off2 of
-    LT -> toLaurentNum off1 (poly1 - Sparse.scale (fromIntegral $ off2 - off1) 1 poly2)
+    LT -> toLaurentNum off1 (poly1 - Multi.scale (fromIntegral $ off2 - off1) 1 poly2)
     EQ -> toLaurentNum off1 (poly1 - poly2)
-    GT -> toLaurentNum off2 (Sparse.scale (fromIntegral $ off1 - off2) 1 poly1 - poly2)
+    GT -> toLaurentNum off2 (Multi.scale (fromIntegral $ off1 - off2) 1 poly1 - poly2)
   negate (Laurent off poly) = Laurent off (negate poly)
   abs = id
   signum = const 1
@@ -190,15 +193,15 @@ instance (Eq a, Num a, G.Vector v (Word, a)) => Num (Laurent v a) where
   {-# INLINE fromInteger #-}
   {-# INLINE (*) #-}
 
-instance (Eq a, Semiring a, G.Vector v (Word, a)) => Semiring (Laurent v a) where
+instance (Eq a, Semiring a, G.Vector v (SU.Vector 1 Word, a)) => Semiring (Laurent v a) where
   zero = Laurent 0 zero
   one  = Laurent 0 one
   Laurent off1 poly1 `times` Laurent off2 poly2 =
     toLaurent (off1 + off2) (poly1 `times` poly2)
   Laurent off1 poly1 `plus` Laurent off2 poly2 = case off1 `compare` off2 of
-    LT -> toLaurent off1 (poly1 `plus` Sparse.scale' (fromIntegral $ off2 - off1) one poly2)
+    LT -> toLaurent off1 (poly1 `plus` Multi.scale' (fromIntegral $ off2 - off1) one poly2)
     EQ -> toLaurent off1 (poly1 `plus` poly2)
-    GT -> toLaurent off2 (Sparse.scale' (fromIntegral $ off1 - off2) one poly1 `plus` poly2)
+    GT -> toLaurent off2 (Multi.scale' (fromIntegral $ off1 - off2) one poly1 `plus` poly2)
   fromNatural n = Laurent 0 (fromNatural n)
   {-# INLINE zero #-}
   {-# INLINE one #-}
@@ -206,63 +209,63 @@ instance (Eq a, Semiring a, G.Vector v (Word, a)) => Semiring (Laurent v a) wher
   {-# INLINE times #-}
   {-# INLINE fromNatural #-}
 
-instance (Eq a, Ring a, G.Vector v (Word, a)) => Ring (Laurent v a) where
+instance (Eq a, Ring a, G.Vector v (SU.Vector 1 Word, a)) => Ring (Laurent v a) where
   negate (Laurent off poly) = Laurent off (Semiring.negate poly)
 
 -- | Create a monomial from a power and a coefficient.
-monomial :: (Eq a, Semiring a, G.Vector v (Word, a)) => Int -> a -> Laurent v a
+monomial :: (Eq a, Semiring a, G.Vector v (SU.Vector 1 Word, a)) => Int -> a -> Laurent v a
 monomial p c
   | c == zero = Laurent 0 zero
-  | otherwise = Laurent p (Sparse.monomial' 0 c)
+  | otherwise = Laurent p (Multi.monomial' 0 c)
 {-# INLINE monomial #-}
 
 -- | Multiply a polynomial by a monomial, expressed as a power and a coefficient.
 --
 -- >>> scale 2 3 (X^2 + 1) :: ULaurent Double
 -- 3.0 * X^4 + 3.0 * X^2
-scale :: (Eq a, Semiring a, G.Vector v (Word, a)) => Int -> a -> Laurent v a -> Laurent v a
-scale yp yc (Laurent off poly) = toLaurent (off + yp) (Sparse.scale' 0 yc poly)
+scale :: (Eq a, Semiring a, G.Vector v (SU.Vector 1 Word, a)) => Int -> a -> Laurent v a -> Laurent v a
+scale yp yc (Laurent off poly) = toLaurent (off + yp) (Multi.scale' 0 yc poly)
 
 -- | Evaluate at a given point.
 --
 -- >>> eval (X^2 + 1 :: ULaurent Double) 3
 -- 10.0
-eval :: (Field a, G.Vector v (Word, a)) => Laurent v a -> a -> a
-eval (Laurent off poly) x = Sparse.eval' poly x `times`
+eval :: (Field a, G.Vector v (SU.Vector 1 Word, a)) => Laurent v a -> a -> a
+eval (Laurent off poly) x = Multi.eval' poly (SV.singleton x) `times`
   (if off >= 0 then x Semiring.^ off else quot one x Semiring.^ (- off))
 {-# INLINE eval #-}
 
--- | Substitute another polynomial instead of 'Data.Poly.Sparse.X'.
+-- | Substitute another polynomial instead of 'Data.Poly.Multi.X'.
 --
 -- >>> import Data.Poly.Sparse (UPoly)
 -- >>> subst (Data.Poly.Sparse.X^2 + 1 :: UPoly Int) (X + 1 :: ULaurent Int)
 -- 1 * X^2 + 2 * X + 2
-subst :: (Eq a, Semiring a, G.Vector v (Word, a), G.Vector w (Word, a)) => Poly v a -> Laurent w a -> Laurent w a
-subst = Sparse.substitute' (scale 0)
+subst :: (Eq a, Semiring a, G.Vector v (SU.Vector 1 Word, a), G.Vector w (SU.Vector 1 Word, a)) => Poly v a -> Laurent w a -> Laurent w a
+subst xs = Multi.substitute' (scale 0) xs . SV.singleton
 {-# INLINE subst #-}
 
 -- | Take a derivative.
 --
 -- >>> deriv (X^3 + 3 * X) :: ULaurent Int
 -- 3 * X^2 + 3
-deriv :: (Eq a, Ring a, G.Vector v (Word, a)) => Laurent v a -> Laurent v a
-deriv (Laurent off (Poly xs)) =
-  toLaurent (off - 1) $ Sparse.toPoly' $ G.map (\(i, x) -> (i, x `times` Semiring.fromIntegral (fromIntegral i + off))) xs
+deriv :: (Eq a, Ring a, G.Vector v (SU.Vector 1 Word, a)) => Laurent v a -> Laurent v a
+deriv (Laurent off (MultiPoly xs)) =
+  toLaurent (off - 1) $ Multi.toMultiPoly' $ G.map (\(i, x) -> (i, x `times` Semiring.fromIntegral (fromIntegral (SU.head i) + off))) xs
 {-# INLINE deriv #-}
 
 -- | Create an identity polynomial.
-pattern X :: (Eq a, Semiring a, G.Vector v (Word, a)) => Laurent v a
+pattern X :: (Eq a, Semiring a, G.Vector v (SU.Vector 1 Word, a)) => Laurent v a
 pattern X <- (isVar -> True)
   where X = var
 
-var :: forall a v. (Eq a, Semiring a, G.Vector v (Word, a)) => Laurent v a
+var :: forall a v. (Eq a, Semiring a, G.Vector v (SU.Vector 1 Word, a)) => Laurent v a
 var
   | (one :: a) == zero = Laurent 0 zero
   | otherwise          = Laurent 1 one
 {-# INLINE var #-}
 
-isVar :: forall v a. (Eq a, Semiring a, G.Vector v (Word, a)) => Laurent v a -> Bool
-isVar (Laurent off (Poly xs))
+isVar :: forall v a. (Eq a, Semiring a, G.Vector v (SU.Vector 1 Word, a)) => Laurent v a -> Bool
+isVar (Laurent off (MultiPoly xs))
   | (one :: a) == zero = off == 0 && G.null xs
   | otherwise          = off == 1 && G.length xs == 1 && G.unsafeHead xs == (0, one)
 {-# INLINE isVar #-}
@@ -273,14 +276,14 @@ isVar (Laurent off (Poly xs))
 -- >>> X + 2 + 3 * X^-1 :: ULaurent Int
 -- 1 * X + 2 + 3 * X^-1
 (^-)
-  :: (Eq a, Semiring a, G.Vector v (Word, a))
+  :: (Eq a, Semiring a, G.Vector v (SU.Vector 1 Word, a))
   => Laurent v a
   -> Int
   -> Laurent v a
 X^-n = monomial (negate n) one
 _^-_ = throw $ PatternMatchFail "(^-) can be applied only to X"
 
-instance (Eq a, Ring a, GcdDomain a, G.Vector v (Word, a)) => GcdDomain (Laurent v a) where
+instance (Eq a, Ring a, GcdDomain a, G.Vector v (SU.Vector 1 Word, a)) => GcdDomain (Laurent v a) where
   divide (Laurent off1 poly1) (Laurent off2 poly2) =
     toLaurent (off1 - off2) <$> divide poly1 poly2
   {-# INLINE divide #-}
diff --git a/src/Data/Poly/Internal/Sparse.hs b/src/Data/Poly/Internal/Sparse.hs
deleted file mode 100644
index b34448a..0000000
--- a/src/Data/Poly/Internal/Sparse.hs
+++ /dev/null
@@ -1,603 +0,0 @@
--- |
--- 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 GeneralizedNewtypeDeriving #-}
-{-# LANGUAGE PatternSynonyms            #-}
-{-# LANGUAGE ScopedTypeVariables        #-}
-{-# LANGUAGE StandaloneDeriving         #-}
-{-# LANGUAGE TypeFamilies               #-}
-{-# LANGUAGE UndecidableInstances       #-}
-{-# LANGUAGE ViewPatterns               #-}
-
-module Data.Poly.Internal.Sparse
-  ( Poly(..)
-  , VPoly
-  , UPoly
-  , leading
-  -- * Num interface
-  , toPoly
-  , monomial
-  , scale
-  , pattern X
-  , eval
-  , subst
-  , deriv
-  , integral
-  -- * Semiring interface
-  , toPoly'
-  , monomial'
-  , scale'
-  , pattern X'
-  , eval'
-  , subst'
-  , substitute'
-  , deriv'
-  , integral'
-  -- * Internals
-  , normalize
-  , plusPoly
-  , minusPoly
-  , convolution
-  , scaleInternal
-  , derivPoly
-  ) where
-
-import Prelude hiding (quot)
-import Control.DeepSeq (NFData)
-import Control.Monad
-import Control.Monad.Primitive
-import Control.Monad.ST
-import Data.Bits
-import Data.Euclidean (Field, quot)
-import Data.List (intersperse)
-import Data.Ord
-import Data.Semiring (Semiring(..), Ring())
-import qualified Data.Semiring as Semiring
-import qualified Data.Vector as V
-import qualified Data.Vector.Generic as G
-import qualified Data.Vector.Generic.Mutable as MG
-import qualified Data.Vector.Unboxed as U
-import qualified Data.Vector.Algorithms.Tim as Tim
-import GHC.Exts
-
--- | Polynomials of one variable with coefficients from @a@,
--- backed by a 'G.Vector' @v@ (boxed, unboxed, storable, etc.).
---
--- Use pattern 'X' for construction:
---
--- >>> (X + 1) + (X - 1) :: VPoly Integer
--- 2 * X
--- >>> (X + 1) * (X - 1) :: UPoly Int
--- 1 * X^2 + (-1)
---
--- Polynomials are stored normalized, without
--- zero coefficients, so 0 * 'X' + 1 equals to 1.
---
--- 'Ord' instance does not make much sense mathematically,
--- it is defined only for the sake of 'Data.Set.Set', 'Data.Map.Map', etc.
---
-newtype Poly v a = Poly
-  { unPoly :: v (Word, a)
-  -- ^ Convert 'Poly' to a vector of coefficients.
-  }
-
-deriving instance Eq     (v (Word, a)) => Eq     (Poly v a)
-deriving instance Ord    (v (Word, a)) => Ord    (Poly v a)
-deriving instance NFData (v (Word, a)) => NFData (Poly v a)
-
-instance (Eq a, Semiring a, G.Vector v (Word, a)) => IsList (Poly v a) where
-  type Item (Poly v a) = (Word, a)
-  fromList = toPoly' . G.fromList
-  fromListN = (toPoly' .) . G.fromListN
-  toList = G.toList . unPoly
-
-instance (Show a, G.Vector v (Word, a)) => Show (Poly v a) where
-  showsPrec d (Poly xs)
-    | G.null xs
-      = showString "0"
-    | otherwise
-      = showParen (d > 0)
-      $ foldl (.) id
-      $ intersperse (showString " + ")
-      $ G.foldl (\acc (i, c) -> showCoeff i c : acc) [] xs
-    where
-      showCoeff 0 c = showsPrec 7 c
-      showCoeff 1 c = showsPrec 7 c . showString " * X"
-      showCoeff i c = showsPrec 7 c . showString " * X^" . showsPrec 7 i
-
--- | Polynomials backed by boxed vectors.
-type VPoly = Poly V.Vector
-
--- | Polynomials backed by unboxed vectors.
-type UPoly = Poly U.Vector
-
--- | Make 'Poly' from a list of (power, coefficient) pairs.
---
--- >>> :set -XOverloadedLists
--- >>> toPoly [(0,1),(1,2),(2,3)] :: VPoly Integer
--- 3 * X^2 + 2 * X + 1
--- >>> toPoly [(0,0),(1,0),(2,0)] :: UPoly Int
--- 0
-toPoly :: (Eq a, Num a, G.Vector v (Word, a)) => 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 (t, a), Ord t)
-  => (a -> Bool)
-  -> (a -> a -> a)
-  -> v (t, a)
-  -> v (t, a)
-normalize p add vs
-  | G.null vs = vs
-  | otherwise = runST $ do
-    ws <- G.thaw vs
-    l' <- normalizeM p add ws
-    G.unsafeFreeze $ MG.unsafeSlice 0 l' ws
-
-normalizeM
-  :: (PrimMonad m, G.Vector v (t, a), Ord t)
-  => (a -> Bool)
-  -> (a -> a -> a)
-  -> G.Mutable v (PrimState m) (t, a)
-  -> m Int
-normalizeM p add ws = do
-    let l = MG.length ws
-    let go i j acc@(accP, accC)
-          | j >= l =
-            if p accC
-              then do
-                MG.write ws i acc
-                pure $ i + 1
-              else pure i
-          | otherwise = do
-            v@(vp, vc) <- MG.unsafeRead ws j
-            if vp == accP
-              then go i (j + 1) (accP, accC `add` vc)
-              else if p accC
-                then do
-                  MG.write ws i acc
-                  go (i + 1) (j + 1) v
-                else go i (j + 1) v
-    Tim.sortBy (comparing fst) ws
-    wsHead <- MG.unsafeRead ws 0
-    go 0 1 wsHead
-
--- | Note that 'abs' = 'id' and 'signum' = 'const' 1.
-instance (Eq a, Num a, G.Vector v (Word, a)) => Num (Poly v a) where
-  Poly xs + Poly ys = Poly $ plusPoly (/= 0) (+) xs ys
-  Poly xs - Poly ys = Poly $ minusPoly (/= 0) negate (-) xs ys
-  negate (Poly xs) = Poly $ G.map (fmap negate) xs
-  abs = id
-  signum = const 1
-  fromInteger n = case fromInteger n of
-    0 -> Poly G.empty
-    m -> Poly $ G.singleton (0, m)
-  Poly xs * Poly ys = Poly $ convolution (/= 0) (+) (*) xs ys
-  {-# INLINE (+) #-}
-  {-# INLINE (-) #-}
-  {-# INLINE negate #-}
-  {-# INLINE fromInteger #-}
-  {-# INLINE (*) #-}
-
-instance (Eq a, Semiring a, G.Vector v (Word, a)) => Semiring (Poly v a) where
-  zero = Poly G.empty
-  one
-    | (one :: a) == zero = zero
-    | otherwise = Poly $ G.singleton (0, one)
-  plus (Poly xs) (Poly ys) = Poly $ plusPoly (/= zero) plus xs ys
-  times (Poly xs) (Poly ys) = Poly $ convolution (/= zero) plus times xs ys
-  {-# INLINE zero #-}
-  {-# INLINE one #-}
-  {-# INLINE plus #-}
-  {-# INLINE times #-}
-
-  fromNatural n = if n' == zero then zero else Poly $ G.singleton (0, n')
-    where
-      n' :: a
-      n' = fromNatural n
-  {-# INLINE fromNatural #-}
-
-instance (Eq a, Ring a, G.Vector v (Word, a)) => Ring (Poly v a) where
-  negate (Poly xs) = Poly $ G.map (fmap Semiring.negate) xs
-
-plusPoly
-  :: (G.Vector v (t, a), Ord t)
-  => (a -> Bool)
-  -> (a -> a -> a)
-  -> v (t, a)
-  -> v (t, a)
-  -> v (t, a)
-plusPoly p add xs ys = runST $ do
-  zs <- MG.unsafeNew (G.length xs + G.length ys)
-  lenZs <- plusPolyM p add xs ys zs
-  G.unsafeFreeze $ MG.unsafeSlice 0 lenZs zs
-{-# INLINABLE plusPoly #-}
-
-plusPolyM
-  :: (PrimMonad m, G.Vector v (t, a), Ord t)
-  => (a -> Bool)
-  -> (a -> a -> a)
-  -> v (t, a)
-  -> v (t, a)
-  -> G.Mutable v (PrimState m) (t, a)
-  -> m Int
-plusPolyM p add xs ys zs = go 0 0 0
-  where
-    lenXs = G.length xs
-    lenYs = G.length ys
-
-    go ix iy iz
-      | ix == lenXs, iy == lenYs = pure iz
-      | ix == lenXs = do
-        G.unsafeCopy
-          (MG.unsafeSlice iz (lenYs - iy) zs)
-          (G.unsafeSlice iy (lenYs - iy) ys)
-        pure $ iz + lenYs - iy
-      | iy == lenYs = do
-        G.unsafeCopy
-          (MG.unsafeSlice iz (lenXs - ix) zs)
-          (G.unsafeSlice ix (lenXs - ix) xs)
-        pure $ iz + lenXs - ix
-      | (xp, xc) <- G.unsafeIndex xs ix
-      , (yp, yc) <- G.unsafeIndex ys iy
-      = case xp `compare` yp of
-        LT -> do
-          MG.unsafeWrite zs iz (xp, xc)
-          go (ix + 1) iy (iz + 1)
-        EQ -> do
-          let zc = xc `add` yc
-          if p zc then do
-            MG.unsafeWrite zs iz (xp, zc)
-            go (ix + 1) (iy + 1) (iz + 1)
-          else
-            go (ix + 1) (iy + 1) iz
-        GT -> do
-          MG.unsafeWrite zs iz (yp, yc)
-          go ix (iy + 1) (iz + 1)
-{-# INLINABLE plusPolyM #-}
-
-minusPoly
-  :: (G.Vector v (t, a), Ord t)
-  => (a -> Bool)
-  -> (a -> a)
-  -> (a -> a -> a)
-  -> v (t, a)
-  -> v (t, a)
-  -> v (t, a)
-minusPoly p neg sub xs ys = runST $ do
-  zs <- MG.unsafeNew (lenXs + lenYs)
-  let go ix iy iz
-        | ix == lenXs, iy == lenYs = pure iz
-        | ix == lenXs = do
-          forM_ [iy .. lenYs - 1] $ \i ->
-            MG.unsafeWrite zs (iz + i - iy)
-              (fmap neg (G.unsafeIndex ys i))
-          pure $ iz + lenYs - iy
-        | iy == lenYs = do
-          G.unsafeCopy
-            (MG.unsafeSlice iz (lenXs - ix) zs)
-            (G.unsafeSlice ix (lenXs - ix) xs)
-          pure $ iz + lenXs - ix
-        | (xp, xc) <- G.unsafeIndex xs ix
-        , (yp, yc) <- G.unsafeIndex ys iy
-        = case xp `compare` yp of
-          LT -> do
-            MG.unsafeWrite zs iz (xp, xc)
-            go (ix + 1) iy (iz + 1)
-          EQ -> do
-            let zc = xc `sub` yc
-            if p zc then do
-              MG.unsafeWrite zs iz (xp, zc)
-              go (ix + 1) (iy + 1) (iz + 1)
-            else
-              go (ix + 1) (iy + 1) iz
-          GT -> do
-            MG.unsafeWrite zs iz (yp, neg yc)
-            go ix (iy + 1) (iz + 1)
-  lenZs <- go 0 0 0
-  G.unsafeFreeze $ MG.unsafeSlice 0 lenZs zs
-  where
-    lenXs = G.length xs
-    lenYs = G.length ys
-{-# INLINABLE minusPoly #-}
-
-scaleM
-  :: (PrimMonad m, G.Vector v (t, a), Num t)
-  => (a -> Bool)
-  -> (a -> a -> a)
-  -> v (t, a)
-  -> (t, a)
-  -> G.Mutable v (PrimState m) (t, a)
-  -> m Int
-scaleM p mul xs (yp, yc) zs = go 0 0
-  where
-    lenXs = G.length xs
-
-    go ix iz
-      | ix == lenXs = pure iz
-      | (xp, xc) <- G.unsafeIndex xs ix
-      = do
-        let zc = xc `mul` yc
-        if p zc then do
-          MG.unsafeWrite zs iz (xp + yp, zc)
-          go (ix + 1) (iz + 1)
-        else
-          go (ix + 1) iz
-{-# INLINABLE scaleM #-}
-
-scaleInternal
-  :: (G.Vector v (t, a), Num t)
-  => (a -> Bool)
-  -> (a -> a -> a)
-  -> t
-  -> a
-  -> v (t, a)
-  -> v (t, a)
-scaleInternal p mul yp yc xs = runST $ do
-  zs <- MG.unsafeNew (G.length xs)
-  len <- scaleM p (flip mul) xs (yp, yc) zs
-  G.unsafeFreeze $ MG.unsafeSlice 0 len zs
-{-# INLINABLE scaleInternal #-}
-
--- | 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 yp yc = Poly . scaleInternal (/= 0) (*) yp yc . unPoly
-
-scale' :: (Eq a, Semiring a, G.Vector v (Word, a)) => Word -> a -> Poly v a -> Poly v a
-scale' yp yc = Poly . scaleInternal (/= zero) times yp yc . unPoly
-
-convolution
-  :: forall v t a.
-     (G.Vector v (t, a), Ord t, Num t)
-  => (a -> Bool)
-  -> (a -> a -> a)
-  -> (a -> a -> a)
-  -> v (t, a)
-  -> v (t, a)
-  -> v (t, a)
-convolution p add mult xs ys
-  | G.length xs >= G.length ys
-  = go mult xs ys
-  | otherwise
-  = go (flip mult) ys xs
-  where
-    go :: (a -> a -> a) -> v (t, a) -> v (t, a) -> v (t, a)
-    go mul long short = runST $ do
-      let lenLong   = G.length long
-          lenShort  = G.length short
-          lenBuffer = lenLong * lenShort
-      slices <- MG.unsafeNew lenShort
-      buffer <- MG.unsafeNew lenBuffer
-
-      forM_ [0 .. lenShort - 1] $ \iShort -> do
-        let (pShort, cShort) = G.unsafeIndex short iShort
-            from = iShort * lenLong
-            bufferSlice = MG.unsafeSlice from lenLong buffer
-        len <- scaleM p mul long (pShort, cShort) bufferSlice
-        MG.unsafeWrite slices iShort (from, len)
-
-      slices' <- G.unsafeFreeze slices
-      buffer' <- G.unsafeFreeze buffer
-      bufferNew <- MG.unsafeNew lenBuffer
-      gogo slices' buffer' bufferNew
-
-    gogo
-      :: PrimMonad m
-      => U.Vector (Int, Int)
-      -> v (t, a)
-      -> G.Mutable v (PrimState m) (t, a)
-      -> m (v (t, a))
-    gogo slices buffer bufferNew
-      | G.length slices == 0
-      = pure G.empty
-      | G.length slices == 1
-      , (from, len) <- G.unsafeIndex slices 0
-      = pure $ G.unsafeSlice from len buffer
-      | otherwise = do
-        let nSlices = G.length slices
-        slicesNew <- MG.unsafeNew ((nSlices + 1) `shiftR` 1)
-        forM_ [0 .. (nSlices - 2) `shiftR` 1] $ \i -> do
-          let (from1, len1) = G.unsafeIndex slices (2 * i)
-              (from2, len2) = G.unsafeIndex slices (2 * i + 1)
-              slice1 = G.unsafeSlice from1 len1 buffer
-              slice2 = G.unsafeSlice from2 len2 buffer
-              slice3 = MG.unsafeSlice from1 (len1 + len2) bufferNew
-          len3 <- plusPolyM p add slice1 slice2 slice3
-          MG.unsafeWrite slicesNew i (from1, len3)
-
-        when (odd nSlices) $ do
-          let (from, len) = G.unsafeIndex slices (nSlices - 1)
-              slice1 = G.unsafeSlice from len buffer
-              slice3 = MG.unsafeSlice from len bufferNew
-          G.unsafeCopy slice3 slice1
-          MG.unsafeWrite slicesNew (nSlices `shiftR` 1) (from, len)
-
-        slicesNew' <- G.unsafeFreeze slicesNew
-        buffer'    <- G.unsafeThaw   buffer
-        bufferNew' <- G.unsafeFreeze bufferNew
-        gogo slicesNew' bufferNew' buffer'
-{-# INLINABLE convolution #-}
-
--- | Create a monomial from a power and a coefficient.
-monomial :: (Eq a, Num a, G.Vector v (Word, a)) => Word -> a -> Poly v a
-monomial _ 0 = Poly G.empty
-monomial p c = Poly $ G.singleton (p, c)
-
-monomial' :: (Eq a, Semiring a, G.Vector v (Word, a)) => Word -> a -> Poly v a
-monomial' p c
-  | c == zero = Poly G.empty
-  | otherwise = Poly $ G.singleton (p, c)
-
-data Strict3 a b c = Strict3 !a !b !c
-
-fst3 :: Strict3 a b c -> a
-fst3 (Strict3 a _ _) = a
-
--- | Evaluate at a given point.
---
--- >>> eval (X^2 + 1 :: UPoly Int) 3
--- 10
-eval :: (Num a, G.Vector v (Word, a)) => Poly v a -> a -> a
-eval = substitute (*)
-{-# INLINE eval #-}
-
-eval' :: (Semiring a, G.Vector v (Word, a)) => Poly v a -> a -> a
-eval' = substitute' times
-{-# INLINE eval' #-}
-
--- | Substitute another polynomial instead of 'X'.
---
--- >>> subst (X^2 + 1 :: UPoly Int) (X + 1 :: UPoly Int)
--- 1 * X^2 + 2 * X + 2
-subst
-  :: (Eq a, Num a, G.Vector v (Word, a), G.Vector w (Word, a))
-  => Poly v a
-  -> Poly w a
-  -> Poly w a
-subst = substitute (scale 0)
-{-# INLINE subst #-}
-
-subst'
-  :: (Eq a, Semiring a, G.Vector v (Word, a), G.Vector w (Word, a))
-  => Poly v a
-  -> Poly w a
-  -> Poly w a
-subst' = substitute' (scale' 0)
-{-# INLINE subst' #-}
-
-substitute :: (G.Vector v (Word, a), Num b) => (a -> b -> b) -> Poly v a -> b -> b
-substitute f (Poly cs) x = fst3 $ G.foldl' go (Strict3 0 0 1) cs
-  where
-    go (Strict3 acc q xq) (p, c) =
-      let xp = xq * x ^ (p - q) in
-        Strict3 (acc + f c xp) p xp
-{-# INLINE substitute #-}
-
-substitute' :: (G.Vector v (Word, a), Semiring b) => (a -> b -> b) -> Poly v a -> b -> b
-substitute' f (Poly cs) x = fst3 $ G.foldl' go (Strict3 zero 0 one) cs
-  where
-    go (Strict3 acc q xq) (p, c) =
-      let xp = xq `times` (if p == q then one else x Semiring.^ (p - q)) in
-        Strict3 (acc `plus` f c xp) p xp
-{-# INLINE substitute' #-}
-
--- | Take a derivative.
---
--- >>> deriv (X^3 + 3 * X) :: UPoly Int
--- 3 * X^2 + 3
-deriv :: (Eq a, Num a, G.Vector v (Word, a)) => Poly v a -> Poly v a
-deriv (Poly xs) = Poly $ derivPoly
-  (/= 0)
-  pred
-  (\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)
-  pred
-  (\p c -> fromNatural (fromIntegral p) `times` c)
-  xs
-{-# INLINE deriv' #-}
-
-derivPoly
-  :: (G.Vector v (t, a))
-  => (a -> Bool)
-  -> (t -> t)
-  -> (t -> a -> a)
-  -> v (t, a)
-  -> v (t, a)
-derivPoly p dec mul xs
-  | G.null xs = G.empty
-  | otherwise = runST $ do
-    let lenXs = G.length xs
-    zs <- MG.unsafeNew lenXs
-    let go ix iz
-          | ix == lenXs = pure iz
-          | (xp, xc) <- G.unsafeIndex xs ix
-          = do
-            let zc = xp `mul` xc
-            if p zc then do
-              MG.unsafeWrite zs iz (dec xp, zc)
-              go (ix + 1) (iz + 1)
-            else
-              go (ix + 1) iz
-    lenZs <- go 0 0
-    G.unsafeFreeze $ MG.unsafeSlice 0 lenZs zs
-{-# INLINABLE derivPoly #-}
-
--- | Compute an indefinite integral of a polynomial,
--- setting constant term to zero.
---
--- >>> integral (3 * X^2 + 3) :: UPoly Double
--- 1.0 * X^3 + 3.0 * X
-integral :: (Fractional a, G.Vector v (Word, a)) => Poly v a -> Poly v a
-integral (Poly xs)
-  = Poly
-  $ G.map (\(p, c) -> (p + 1, c / (fromIntegral p + 1))) xs
-{-# INLINE integral #-}
-
-integral' :: (Field a, G.Vector v (Word, a)) => Poly v a -> Poly v a
-integral' (Poly xs)
-  = Poly
-  $ G.map (\(p, c) -> (p + 1, c `quot` Semiring.fromIntegral (p + 1))) xs
-{-# INLINE integral' #-}
-
--- | Create an identity polynomial.
-pattern X :: (Eq a, Num a, G.Vector v (Word, a)) => Poly v a
-pattern X <- (isVar -> True)
-  where X = var
-
-var :: forall a v. (Eq a, Num a, G.Vector v (Word, a)) => Poly v a
-var
-  | (1 :: a) == 0 = Poly G.empty
-  | otherwise     = Poly $ G.singleton (1, 1)
-{-# INLINE var #-}
-
-isVar :: forall v a. (Eq a, Num a, G.Vector v (Word, a)) => Poly v a -> Bool
-isVar (Poly xs)
-  | (1 :: a) == 0 = G.null xs
-  | otherwise     = G.length xs == 1 && G.unsafeHead xs == (1, 1)
-{-# INLINE isVar #-}
-
--- | Create an identity polynomial.
-pattern X' :: (Eq a, Semiring a, G.Vector v (Word, a)) => Poly v a
-pattern X' <- (isVar' -> True)
-  where X' = var'
-
-var' :: forall a v. (Eq a, Semiring a, G.Vector v (Word, a)) => Poly v a
-var'
-  | (one :: a) == zero = Poly G.empty
-  | otherwise          = Poly $ G.singleton (1, one)
-{-# INLINE var' #-}
-
-isVar' :: forall v a. (Eq a, Semiring a, G.Vector v (Word, a)) => Poly v a -> Bool
-isVar' (Poly xs)
-  | (one :: a) == zero = G.null xs
-  | otherwise          = G.length xs == 1 && G.unsafeHead xs == (1, one)
-{-# INLINE isVar' #-}
diff --git a/src/Data/Poly/Internal/Sparse/GcdDomain.hs b/src/Data/Poly/Internal/Sparse/GcdDomain.hs
deleted file mode 100644
index 1bc8ca9..0000000
--- a/src/Data/Poly/Internal/Sparse/GcdDomain.hs
+++ /dev/null
@@ -1,96 +0,0 @@
--- |
--- 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 ScopedTypeVariables        #-}
-{-# LANGUAGE TypeFamilies               #-}
-{-# LANGUAGE UndecidableInstances       #-}
-
-{-# OPTIONS_GHC -fno-warn-orphans #-}
-
-module Data.Poly.Internal.Sparse.GcdDomain
-  () where
-
-import Prelude hiding (gcd, lcm, (^))
-import Control.Exception
-import Data.Euclidean
-import Data.Maybe
-import Data.Semiring (Semiring(..), Ring(), minus)
-import qualified Data.Vector.Generic as G
-
-import Data.Poly.Internal.Sparse
-
-instance (Eq a, Ring a, GcdDomain a, G.Vector v (Word, a)) => GcdDomain (Poly v a) where
-  divide xs ys = case leading ys of
-    Nothing -> throw DivideByZero
-    Just (yp, yc) -> case leading xs of
-      Nothing -> Just xs
-      Just (xp, xc)
-        | xp < yp -> Nothing
-        | otherwise -> do
-          zc <- divide xc yc
-          let z = Poly $ G.singleton (xp - yp, zc)
-          rest <- divide (xs `minus` z `times` ys) ys
-          pure $ rest `plus` z
-
-  gcd xs ys
-    | G.null (unPoly xs) = ys
-    | G.null (unPoly ys) = xs
-    | G.length (unPoly xs) == 1 = gcdSingleton (G.unsafeHead (unPoly xs)) ys
-    | G.length (unPoly ys) == 1 = gcdSingleton (G.unsafeHead (unPoly ys)) xs
-    | otherwise = times xy (divide' zs (monomial' 0 (content zs)))
-      where
-        zs = gcdHelper xs ys
-        xy = monomial' 0 (gcd (content xs) (content ys))
-
-  lcm x@(Poly xs) y@(Poly ys)
-    | G.null xs || G.null ys = zero
-    | otherwise = (x `divide'` gcd x y) `times` y
-
-  coprime x y = isJust (one `divide` gcd x y)
-
-content :: (GcdDomain a, G.Vector v (Word, a)) => Poly v a -> a
-content (Poly ts) = G.foldl' (\acc (_, t) -> gcd acc t) zero ts
-
-gcdSingleton
-  :: (Eq a, GcdDomain a, G.Vector v (Word, a))
-  => (Word, a)
-  -> Poly v a
-  -> Poly v a
-gcdSingleton (p, c) pcs = monomial'
-  (min p (fst (G.unsafeHead (unPoly pcs))))
-  (gcd c (content pcs))
-
-gcdHelper
-  :: (Eq a, Ring a, GcdDomain a, G.Vector v (Word, a))
-  => Poly v a
-  -> Poly v a
-  -> Poly v a
-gcdHelper xs ys = case (leading xs, leading ys) of
-  (Nothing, _) -> ys
-  (_, Nothing) -> xs
-  (Just (xp, xc), Just (yp, yc))
-    | yp <= xp
-    , Just xy <- xc `divide` yc
-    -> gcdHelper ys (xs `minus` ys `times` monomial' (xp - yp) xy)
-    | xp <= yp
-    , Just yx <- yc `divide` xc
-    -> gcdHelper xs (ys `minus` xs `times` monomial' (yp - xp) yx)
-    | yp <= xp
-    -> gcdHelper ys (xs `times` monomial' 0 gx `minus` ys `times` monomial' (xp - yp) gy)
-    | otherwise
-    -> gcdHelper xs (ys `times` monomial' 0 gy `minus` xs `times` monomial' (yp - xp) gx)
-    where
-      g = lcm xc yc
-      gx = divide' g xc
-      gy = divide' g yc
-
-divide' :: GcdDomain a => a -> a -> a
-divide' = (fromMaybe (error "gcd: violated internal invariant") .) . divide
diff --git a/src/Data/Poly/Multi.hs b/src/Data/Poly/Multi.hs
new file mode 100644
index 0000000..29d1521
--- /dev/null
+++ b/src/Data/Poly/Multi.hs
@@ -0,0 +1,36 @@
+-- |
+-- Module:      Data.Poly.Multi
+-- Copyright:   (c) 2020 Andrew Lelechenko
+-- Licence:     BSD3
+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>
+--
+-- Multivariate sparse polynomials with 'Num' instance.
+--
+
+{-# LANGUAGE DataKinds        #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE PatternSynonyms  #-}
+
+module Data.Poly.Multi
+  ( MultiPoly
+  , VMultiPoly
+  , UMultiPoly
+  , unMultiPoly
+  , toMultiPoly
+  , monomial
+  , scale
+  , pattern X
+  , pattern Y
+  , pattern Z
+  , eval
+  , subst
+  , deriv
+  , integral
+  -- * Conversions
+  , segregate
+  , unsegregate
+  ) where
+
+import Data.Poly.Internal.Multi
+import Data.Poly.Internal.Multi.Field ()
+import Data.Poly.Internal.Multi.GcdDomain ()
diff --git a/src/Data/Poly/Multi/Laurent.hs b/src/Data/Poly/Multi/Laurent.hs
new file mode 100644
index 0000000..cc6a52f
--- /dev/null
+++ b/src/Data/Poly/Multi/Laurent.hs
@@ -0,0 +1,13 @@
+-- |
+-- Module:      Data.Poly.Multi.Laurent
+-- Copyright:   (c) 2020 Andrew Lelechenko
+-- Licence:     BSD3
+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>
+--
+-- Multivariate sparse
+-- <https://en.wikipedia.org/wiki/Laurent_polynomial Laurent polynomials>.
+--
+
+module Data.Poly.Multi.Laurent
+  (
+  ) where
diff --git a/src/Data/Poly/Multi/Semiring.hs b/src/Data/Poly/Multi/Semiring.hs
new file mode 100644
index 0000000..a41053f
--- /dev/null
+++ b/src/Data/Poly/Multi/Semiring.hs
@@ -0,0 +1,164 @@
+-- |
+-- Module:      Data.Poly.Multi.Semiring
+-- Copyright:   (c) 2020 Andrew Lelechenko
+-- Licence:     BSD3
+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>
+--
+-- Multivariate sparse polynomials with 'Semiring' instance.
+--
+
+{-# LANGUAGE CPP              #-}
+{-# LANGUAGE DataKinds        #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE PatternSynonyms  #-}
+{-# LANGUAGE TypeFamilies     #-}
+{-# LANGUAGE TypeOperators    #-}
+
+module Data.Poly.Multi.Semiring
+  ( MultiPoly
+  , VMultiPoly
+  , UMultiPoly
+  , unMultiPoly
+  , toMultiPoly
+  , monomial
+  , scale
+  , pattern X
+  , pattern Y
+  , pattern Z
+  , eval
+  , subst
+  , deriv
+  , integral
+  -- * Conversions
+  , segregate
+  , unsegregate
+  ) where
+
+import Data.Finite
+import Data.Euclidean (Field)
+import Data.Semiring (Semiring(..))
+import qualified Data.Vector.Generic as G
+import qualified Data.Vector.Generic.Sized as SG
+import qualified Data.Vector.Sized as SV
+import qualified Data.Vector.Unboxed.Sized as SU
+
+import Data.Poly.Internal.Multi (MultiPoly, VMultiPoly, UMultiPoly, unMultiPoly, segregate, unsegregate)
+import qualified Data.Poly.Internal.Multi as Multi
+import Data.Poly.Internal.Multi.Field ()
+import Data.Poly.Internal.Multi.GcdDomain ()
+
+#if MIN_VERSION_base(4,10,0)
+import GHC.TypeNats (KnownNat, type (<=))
+#else
+import GHC.TypeLits (KnownNat, type (<=))
+#endif
+
+-- | Make 'MultiPoly' from a list of (powers, coefficient) pairs.
+--
+-- >>> :set -XOverloadedLists -XDataKinds
+-- >>> import Data.Vector.Generic.Sized (fromTuple)
+-- >>> toMultiPoly [(fromTuple (0,0),1),(fromTuple (0,1),2),(fromTuple (1,0),3)] :: VMultiPoly 2 Integer
+-- 3 * X + 2 * Y + 1
+-- >>> toMultiPoly [(fromTuple (0,0),0),(fromTuple (0,1),0),(fromTuple (1,0),0)] :: UMultiPoly 2 Int
+-- 0
+toMultiPoly
+  :: (Eq a, Semiring a, G.Vector v (SU.Vector n Word, a))
+  => v (SU.Vector n Word, a)
+  -> MultiPoly v n a
+toMultiPoly = Multi.toMultiPoly'
+
+-- | Create a monomial from powers and a coefficient.
+monomial
+  :: (Eq a, Semiring a, G.Vector v (SU.Vector n Word, a))
+  => SU.Vector n Word
+  -> a
+  -> MultiPoly v n a
+monomial = Multi.monomial'
+
+-- | Multiply a polynomial by a monomial, expressed as powers and a coefficient.
+--
+-- >>> :set -XDataKinds
+-- >>> import Data.Vector.Generic.Sized (fromTuple)
+-- >>> scale (fromTuple (1, 1)) 3 (X + Y) :: UMultiPoly 2 Int
+-- 3 * X^2 * Y + 3 * X * Y^2
+scale
+  :: (Eq a, Semiring a, KnownNat n, G.Vector v (SU.Vector n Word, a))
+  => SU.Vector n Word
+  -> a
+  -> MultiPoly v n a
+  -> MultiPoly v n a
+scale = Multi.scale'
+
+-- | Create a polynomial equal to the first variable.
+pattern X
+  :: (Eq a, Semiring a, KnownNat n, 1 <= n, G.Vector v (SU.Vector n Word, a))
+  => MultiPoly v n a
+pattern X = Multi.X'
+
+-- | Create a polynomial equal to the second variable.
+pattern Y
+  :: (Eq a, Semiring a, KnownNat n, 2 <= n, G.Vector v (SU.Vector n Word, a))
+  => MultiPoly v n a
+pattern Y = Multi.Y'
+
+-- | Create a polynomial equal to the third variable.
+pattern Z
+  :: (Eq a, Semiring a, KnownNat n, 3 <= n, G.Vector v (SU.Vector n Word, a))
+  => MultiPoly v n a
+pattern Z = Multi.Z'
+
+-- | Evaluate at a given point.
+--
+-- >>> :set -XDataKinds
+-- >>> import Data.Vector.Generic.Sized (fromTuple)
+-- >>> eval (X^2 + Y^2 :: UMultiPoly 2 Int) (fromTuple (3, 4) :: Data.Vector.Sized.Vector 2 Int)
+-- 25
+eval
+  :: (Semiring a, G.Vector v (SU.Vector n Word, a), G.Vector u a)
+  => MultiPoly v n a
+  -> SG.Vector u n a
+  -> a
+eval = Multi.eval'
+
+-- | Substitute another polynomials instead of variables.
+--
+-- >>> :set -XDataKinds
+-- >>> import Data.Vector.Generic.Sized (fromTuple)
+-- >>> subst (X^2 + Y^2 + Z^2 :: UMultiPoly 3 Int) (fromTuple (X + 1, Y + 1, X + Y :: UMultiPoly 2 Int))
+-- 2 * X^2 + 2 * X * Y + 2 * X + 2 * Y^2 + 2 * Y + 2
+subst
+  :: (Eq a, Semiring a, KnownNat m, G.Vector v (SU.Vector n Word, a), G.Vector w (SU.Vector m Word, a))
+  => MultiPoly v n a
+  -> SV.Vector n (MultiPoly w m a)
+  -> MultiPoly w m a
+subst = Multi.subst'
+
+-- | Take a derivative with respect to the /i/-th variable.
+--
+-- >>> :set -XDataKinds
+-- >>> deriv 0 (X^3 + 3 * Y) :: UMultiPoly 2 Int
+-- 3 * X^2
+-- >>> deriv 1 (X^3 + 3 * Y) :: UMultiPoly 2 Int
+-- 3
+deriv
+  :: (Eq a, Semiring a, G.Vector v (SU.Vector n Word, a))
+  => Finite n
+  -> MultiPoly v n a
+  -> MultiPoly v n a
+deriv = Multi.deriv'
+
+-- | Compute an indefinite integral of a polynomial
+-- by the /i/-th variable,
+-- setting constant term to zero.
+--
+-- >>> :set -XDataKinds
+-- >>> integral 0 (3 * X^2 + 2 * Y) :: UMultiPoly 2 Double
+-- 1.0 * X^3 + 2.0 * X * Y
+-- >>> integral 1 (3 * X^2 + 2 * Y) :: UMultiPoly 2 Double
+-- 3.0 * X^2 * Y + 1.0 * Y^2
+integral
+  :: (Field a, G.Vector v (SU.Vector n Word, a))
+  => Finite n
+  -> MultiPoly v n a
+  -> MultiPoly v n a
+integral = Multi.integral'
diff --git a/src/Data/Poly/Semiring.hs b/src/Data/Poly/Semiring.hs
index 97e9f8f..7a1452f 100644
--- a/src/Data/Poly/Semiring.hs
+++ b/src/Data/Poly/Semiring.hs
@@ -7,6 +7,7 @@
 -- Dense polynomials and a 'Semiring'-based interface.
 --
 
+{-# LANGUAGE DataKinds        #-}
 {-# LANGUAGE FlexibleContexts #-}
 {-# LANGUAGE PatternSynonyms  #-}
 
@@ -37,6 +38,7 @@ import Data.Bits
 import Data.Euclidean (Field)
 import Data.Semiring (Semiring(..))
 import qualified Data.Vector.Generic as G
+import qualified Data.Vector.Unboxed.Sized as SU
 
 import qualified Data.Poly.Internal.Convert as Convert
 import Data.Poly.Internal.Dense (Poly(..), VPoly, UPoly, leading)
@@ -44,7 +46,7 @@ import qualified Data.Poly.Internal.Dense as Dense
 import Data.Poly.Internal.Dense.Field ()
 import Data.Poly.Internal.Dense.DFT
 import Data.Poly.Internal.Dense.GcdDomain ()
-import qualified Data.Poly.Internal.Sparse as Sparse
+import qualified Data.Poly.Internal.Multi as Sparse
 
 -- | Make 'Poly' from a vector of coefficients
 -- (first element corresponds to a constant term).
@@ -130,7 +132,7 @@ dftMult getPrimRoot (Poly xs) (Poly ys) =
 -- >>> :set -XFlexibleContexts
 -- >>> denseToSparse (1 `plus` Data.Poly.X^2) :: Data.Poly.Sparse.UPoly Int
 -- 1 * X^2 + 1
-denseToSparse :: (Eq a, Semiring a, G.Vector v a, G.Vector v (Word, a)) => Dense.Poly v a -> Sparse.Poly v a
+denseToSparse :: (Eq a, Semiring a, G.Vector v a, G.Vector v (SU.Vector 1 Word, a)) => Dense.Poly v a -> Sparse.Poly v a
 denseToSparse = Convert.denseToSparse'
 
 -- | Convert from sparse to dense polynomials.
@@ -138,5 +140,5 @@ denseToSparse = Convert.denseToSparse'
 -- >>> :set -XFlexibleContexts
 -- >>> sparseToDense (1 `plus` Data.Poly.Sparse.X^2) :: Data.Poly.UPoly Int
 -- 1 * X^2 + 0 * X + 1
-sparseToDense :: (Semiring a, G.Vector v a, G.Vector v (Word, a)) => Sparse.Poly v a -> Dense.Poly v a
+sparseToDense :: (Semiring a, G.Vector v a, G.Vector v (SU.Vector 1 Word, a)) => Sparse.Poly v a -> Dense.Poly v a
 sparseToDense = Convert.sparseToDense'
diff --git a/src/Data/Poly/Sparse.hs b/src/Data/Poly/Sparse.hs
index 283ef2d..ccdb67b 100644
--- a/src/Data/Poly/Sparse.hs
+++ b/src/Data/Poly/Sparse.hs
@@ -7,7 +7,9 @@
 -- Sparse polynomials with 'Num' instance.
 --
 
-{-# LANGUAGE PatternSynonyms #-}
+{-# LANGUAGE DataKinds        #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE PatternSynonyms  #-}
 
 module Data.Poly.Sparse
   ( Poly
@@ -28,7 +30,89 @@ module Data.Poly.Sparse
   , sparseToDense
   ) where
 
+import Control.Arrow
+import qualified Data.Vector.Generic as G
+import qualified Data.Vector.Unboxed.Sized as SU
+import qualified Data.Vector.Sized as SV
+
 import Data.Poly.Internal.Convert
-import Data.Poly.Internal.Sparse
-import Data.Poly.Internal.Sparse.Field ()
-import Data.Poly.Internal.Sparse.GcdDomain ()
+import Data.Poly.Internal.Multi (Poly, VPoly, UPoly, unPoly, leading, pattern X)
+import qualified Data.Poly.Internal.Multi as Multi
+import Data.Poly.Internal.Multi.Field ()
+import Data.Poly.Internal.Multi.GcdDomain ()
+
+-- | Make 'Poly' from a list of (power, coefficient) pairs.
+--
+-- >>> :set -XOverloadedLists
+-- >>> toPoly [(0,1),(1,2),(2,3)] :: VPoly Integer
+-- 3 * X^2 + 2 * X + 1
+-- >>> toPoly [(0,0),(1,0),(2,0)] :: UPoly Int
+-- 0
+toPoly
+  :: (Eq a, Num a, G.Vector v (Word, a), G.Vector v (SU.Vector 1 Word, a))
+  => v (Word, a)
+  -> Poly v a
+toPoly = Multi.toMultiPoly . G.map (first SU.singleton)
+
+-- | Create a monomial from a power and a coefficient.
+monomial
+  :: (Eq a, Num a, G.Vector v (SU.Vector 1 Word, a))
+  => Word
+  -> a
+  -> Poly v a
+monomial = Multi.monomial . SU.singleton
+
+-- | Multiply a polynomial by a monomial, expressed as a power and a coefficient.
+--
+-- >>> scale 2 3 (X^2 + 1) :: UPoly Int
+-- 3 * X^4 + 3 * X^2
+scale
+  :: (Eq a, Num a, G.Vector v (SU.Vector 1 Word, a))
+  => Word
+  -> a
+  -> Poly v a
+  -> Poly v a
+scale = Multi.scale . SU.singleton
+
+-- | Evaluate at a given point.
+--
+-- >>> eval (X^2 + 1 :: UPoly Int) 3
+-- 10
+eval
+  :: (Num a, G.Vector v (SU.Vector 1 Word, a))
+  => Poly v a
+  -> a
+  -> a
+eval p = Multi.eval p . SV.singleton
+
+-- | Substitute another polynomial instead of 'X'.
+--
+-- >>> subst (X^2 + 1 :: UPoly Int) (X + 1 :: UPoly Int)
+-- 1 * X^2 + 2 * X + 2
+subst
+  :: (Eq a, Num a, G.Vector v (SU.Vector 1 Word, a), G.Vector w (SU.Vector 1 Word, a))
+  => Poly v a
+  -> Poly w a
+  -> Poly w a
+subst p = Multi.subst p . SV.singleton
+
+-- | Take a derivative.
+--
+-- >>> deriv (X^3 + 3 * X) :: UPoly Int
+-- 3 * X^2 + 3
+deriv
+  :: (Eq a, Num a, G.Vector v (SU.Vector 1 Word, a))
+  => Poly v a
+  -> Poly v a
+deriv = Multi.deriv 0
+
+-- | Compute an indefinite integral of a polynomial,
+-- setting constant term to zero.
+--
+-- >>> integral (3 * X^2 + 3) :: UPoly Double
+-- 1.0 * X^3 + 3.0 * X
+integral
+  :: (Fractional a, G.Vector v (SU.Vector 1 Word, a))
+  => Poly v a
+  -> Poly v a
+integral = Multi.integral 0
diff --git a/src/Data/Poly/Sparse/Laurent.hs b/src/Data/Poly/Sparse/Laurent.hs
index 1236438..9e80a1a 100644
--- a/src/Data/Poly/Sparse/Laurent.hs
+++ b/src/Data/Poly/Sparse/Laurent.hs
@@ -4,10 +4,11 @@
 -- Licence:     BSD3
 -- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>
 --
--- Sparse <https://en.wikipedia.org/wiki/Laurent_polynomial Laurent polynomials>.
+-- Sparse
+-- <https://en.wikipedia.org/wiki/Laurent_polynomial Laurent polynomials>.
 --
 
-{-# LANGUAGE PatternSynonyms #-}
+{-# LANGUAGE PatternSynonyms            #-}
 
 module Data.Poly.Sparse.Laurent
   ( Laurent
@@ -25,4 +26,4 @@ module Data.Poly.Sparse.Laurent
   , deriv
   ) where
 
-import Data.Poly.Internal.Sparse.Laurent
+import Data.Poly.Internal.Multi.Laurent
diff --git a/src/Data/Poly/Sparse/Semiring.hs b/src/Data/Poly/Sparse/Semiring.hs
index c2ffd06..69d0bd4 100644
--- a/src/Data/Poly/Sparse/Semiring.hs
+++ b/src/Data/Poly/Sparse/Semiring.hs
@@ -7,6 +7,7 @@
 -- Sparse polynomials with 'Semiring' instance.
 --
 
+{-# LANGUAGE DataKinds        #-}
 {-# LANGUAGE FlexibleContexts #-}
 {-# LANGUAGE PatternSynonyms  #-}
 
@@ -29,16 +30,19 @@ module Data.Poly.Sparse.Semiring
   , sparseToDense
   ) where
 
+import Control.Arrow
 import Data.Euclidean (Field)
 import Data.Semiring (Semiring(..))
 import qualified Data.Vector.Generic as G
+import qualified Data.Vector.Unboxed.Sized as SU
+import qualified Data.Vector.Sized as SV
 
 import qualified Data.Poly.Internal.Convert as Convert
 import qualified Data.Poly.Internal.Dense as Dense
-import Data.Poly.Internal.Sparse (Poly(..), VPoly, UPoly, leading)
-import qualified Data.Poly.Internal.Sparse as Sparse
-import Data.Poly.Internal.Sparse.Field ()
-import Data.Poly.Internal.Sparse.GcdDomain ()
+import Data.Poly.Internal.Multi (Poly, VPoly, UPoly, unPoly, leading)
+import qualified Data.Poly.Internal.Multi as Multi
+import Data.Poly.Internal.Multi.Field ()
+import Data.Poly.Internal.Multi.GcdDomain ()
 
 -- | Make 'Poly' from a list of (power, coefficient) pairs.
 --
@@ -47,59 +51,90 @@ import Data.Poly.Internal.Sparse.GcdDomain ()
 -- 3 * X^2 + 2 * X + 1
 -- >>> toPoly [(0,0),(1,0),(2,0)] :: UPoly Int
 -- 0
-toPoly :: (Eq a, Semiring a, G.Vector v (Word, a)) => v (Word, a) -> Poly v a
-toPoly = Sparse.toPoly'
+toPoly
+  :: (Eq a, Semiring a, G.Vector v (Word, a), G.Vector v (SU.Vector 1 Word, a))
+  => v (Word, a)
+  -> Poly v a
+toPoly = Multi.toMultiPoly' . G.map (first SU.singleton)
 
 -- | Create a monomial from a power and a coefficient.
-monomial :: (Eq a, Semiring a, G.Vector v (Word, a)) => Word -> a -> Poly v a
-monomial = Sparse.monomial'
+monomial
+  :: (Eq a, Semiring a, G.Vector v (SU.Vector 1 Word, a))
+  => Word
+  -> a
+  -> Poly v a
+monomial = Multi.monomial' . SU.singleton
 
 -- | Multiply a polynomial by a monomial, expressed as a power and a coefficient.
 --
 -- >>> scale 2 3 (X^2 + 1) :: UPoly Int
 -- 3 * X^4 + 3 * X^2
-scale :: (Eq a, Semiring a, G.Vector v (Word, a)) => Word -> a -> Poly v a -> Poly v a
-scale = Sparse.scale'
+scale
+  :: (Eq a, Semiring a, G.Vector v (SU.Vector 1 Word, a))
+  => Word
+  -> a
+  -> Poly v a
+  -> Poly v a
+scale = Multi.scale' . SU.singleton
 
 -- | Create an identity polynomial.
-pattern X :: (Eq a, Semiring a, G.Vector v (Word, a)) => Poly v a
-pattern X = Sparse.X'
+pattern X
+  :: (Eq a, Semiring a, G.Vector v (SU.Vector 1 Word, a))
+  => Poly v a
+pattern X = Multi.X'
 
 -- | Evaluate at a given point.
 --
 -- >>> eval (X^2 + 1 :: UPoly Int) 3
 -- 10
-eval :: (Semiring a, G.Vector v (Word, a)) => Poly v a -> a -> a
-eval = Sparse.eval'
+eval
+  :: (Semiring a, G.Vector v (SU.Vector 1 Word, a))
+  => Poly v a
+  -> a
+  -> a
+eval p = Multi.eval' p . SV.singleton
 
 -- | Substitute another polynomial instead of 'X'.
 --
 -- >>> subst (X^2 + 1 :: UPoly Int) (X + 1 :: UPoly Int)
 -- 1 * X^2 + 2 * X + 2
-subst :: (Eq a, Semiring a, G.Vector v (Word, a), G.Vector w (Word, a)) => Poly v a -> Poly w a -> Poly w a
-subst = Sparse.subst'
+subst
+  :: (Eq a, Semiring a, G.Vector v (SU.Vector 1 Word, a), G.Vector w (SU.Vector 1 Word, a))
+  => Poly v a
+  -> Poly w a
+  -> Poly w a
+subst p = Multi.subst' p . SV.singleton
 
 -- | Take a derivative.
 --
 -- >>> deriv (X^3 + 3 * X) :: UPoly Int
 -- 3 * X^2 + 3
-deriv :: (Eq a, Semiring a, G.Vector v (Word, a)) => Poly v a -> Poly v a
-deriv = Sparse.deriv'
+deriv
+  :: (Eq a, Semiring a, G.Vector v (SU.Vector 1 Word, a))
+  => Poly v a
+  -> Poly v a
+deriv = Multi.deriv' 0
 
 -- | Compute an indefinite integral of a polynomial,
 -- setting constant term to zero.
 --
 -- >>> integral (3 * X^2 + 3) :: UPoly Double
 -- 1.0 * X^3 + 3.0 * X
-integral :: (Field a, G.Vector v (Word, a)) => Poly v a -> Poly v a
-integral = Sparse.integral'
+integral
+  :: (Field a, G.Vector v (SU.Vector 1 Word, a))
+  => Poly v a
+  -> Poly v a
+integral = Multi.integral' 0
 
 -- | Convert from dense to sparse polynomials.
 --
 -- >>> :set -XFlexibleContexts
 -- >>> denseToSparse (1 `plus` Data.Poly.X^2) :: Data.Poly.Sparse.UPoly Int
 -- 1 * X^2 + 1
-denseToSparse :: (Eq a, Semiring a, G.Vector v a, G.Vector v (Word, a)) => Dense.Poly v a -> Sparse.Poly v a
+denseToSparse
+  :: (Eq a, Semiring a, G.Vector v a, G.Vector v (SU.Vector 1 Word, a))
+  => Dense.Poly v a
+  -> Multi.Poly v a
 denseToSparse = Convert.denseToSparse'
 
 -- | Convert from sparse to dense polynomials.
@@ -107,5 +142,8 @@ denseToSparse = Convert.denseToSparse'
 -- >>> :set -XFlexibleContexts
 -- >>> sparseToDense (1 `plus` Data.Poly.Sparse.X^2) :: Data.Poly.UPoly Int
 -- 1 * X^2 + 0 * X + 1
-sparseToDense :: (Semiring a, G.Vector v a, G.Vector v (Word, a)) => Sparse.Poly v a -> Dense.Poly v a
+sparseToDense
+  :: (Semiring a, G.Vector v a, G.Vector v (SU.Vector 1 Word, a))
+  => Multi.Poly v a
+  -> Dense.Poly v a
 sparseToDense = Convert.sparseToDense'
diff --git a/test/DFT.hs b/test/DFT.hs
index 51804ea..7c5c4d5 100644
--- a/test/DFT.hs
+++ b/test/DFT.hs
@@ -2,8 +2,6 @@
 {-# LANGUAGE DataKinds                  #-}
 {-# LANGUAGE TypeOperators              #-}
 
-{-# OPTIONS_GHC -fno-warn-orphans #-}
-
 module DFT
   ( testSuite
   ) where
diff --git a/test/Dense.hs b/test/Dense.hs
index 749a956..18091ba 100644
--- a/test/Dense.hs
+++ b/test/Dense.hs
@@ -4,8 +4,6 @@
 {-# LANGUAGE GeneralizedNewtypeDeriving #-}
 {-# LANGUAGE ScopedTypeVariables        #-}
 
-{-# OPTIONS_GHC -fno-warn-orphans #-}
-
 module Dense
   ( testSuite
   , ShortPoly(..)
@@ -29,17 +27,6 @@ import Test.Tasty.QuickCheck hiding (scale, numTests)
 import Quaternion
 import TestUtils
 
-instance (Eq a, Semiring a, Arbitrary a, G.Vector v a) => Arbitrary (Poly v a) where
-  arbitrary = S.toPoly . G.fromList <$> arbitrary
-  shrink = fmap (S.toPoly . G.fromList) . shrink . G.toList . unPoly
-
-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
diff --git a/test/DenseLaurent.hs b/test/DenseLaurent.hs
index 52a56bd..aeb302a 100644
--- a/test/DenseLaurent.hs
+++ b/test/DenseLaurent.hs
@@ -3,8 +3,6 @@
 {-# LANGUAGE GeneralizedNewtypeDeriving #-}
 {-# LANGUAGE ScopedTypeVariables        #-}
 
-{-# OPTIONS_GHC -fno-warn-orphans #-}
-
 module DenseLaurent
   ( testSuite
   ) where
@@ -12,7 +10,7 @@ module DenseLaurent
 import Prelude hiding (gcd, quotRem, quot, rem)
 import Control.Exception
 import Data.Either
-import Data.Euclidean (Euclidean(..), GcdDomain(..), Field)
+import Data.Euclidean (GcdDomain(..), Field)
 import Data.Int
 import qualified Data.Poly
 import Data.Poly.Laurent
@@ -24,21 +22,9 @@ import qualified Data.Vector.Unboxed as U
 import Test.Tasty
 import Test.Tasty.QuickCheck hiding (scale, numTests)
 
-import Dense (ShortPoly(..))
 import Quaternion
 import TestUtils
 
-instance (Eq a, Semiring a, Arbitrary a, G.Vector v a) => Arbitrary (Laurent v a) where
-  arbitrary = toLaurent <$> ((`rem` 10) <$> arbitrary) <*> arbitrary
-  shrink = fmap (uncurry toLaurent) . shrink . unLaurent
-
-newtype ShortLaurent a = ShortLaurent { unShortLaurent :: a }
-  deriving (Eq, Show, Semiring, GcdDomain)
-
-instance (Eq a, Semiring a, Arbitrary a, G.Vector v a) => Arbitrary (ShortLaurent (Laurent v a)) where
-  arbitrary = (ShortLaurent .) . toLaurent <$> ((`rem` 10) <$> arbitrary) <*> (unShortPoly <$> arbitrary)
-  shrink = fmap (ShortLaurent . uncurry toLaurent . fmap unShortPoly) . shrink . fmap ShortPoly . unLaurent . unShortLaurent
-
 testSuite :: TestTree
 testSuite = testGroup "DenseLaurent"
   [ otherTests
@@ -78,8 +64,8 @@ numTests =
 
 gcdDomainTests :: [TestTree]
 gcdDomainTests =
-  [ myGcdDomainLaws (Proxy :: Proxy (ShortLaurent (VLaurent Integer)))
-  , myGcdDomainLaws (Proxy :: Proxy (ShortLaurent (VLaurent Rational)))
+  [ myGcdDomainLaws (Proxy :: Proxy (ShortPoly (VLaurent Integer)))
+  , myGcdDomainLaws (Proxy :: Proxy (ShortPoly (VLaurent Rational)))
   ]
 
 showTests :: [TestTree]
diff --git a/test/Multi.hs b/test/Multi.hs
index 4026147..eb46f14 100644
--- a/test/Multi.hs
+++ b/test/Multi.hs
@@ -6,22 +6,17 @@
 {-# LANGUAGE ScopedTypeVariables        #-}
 {-# LANGUAGE UndecidableInstances       #-}
 
-{-# OPTIONS_GHC -fno-warn-orphans #-}
-
 module Multi
   ( testSuite
   ) where
 
 import Prelude hiding (gcd, quotRem, rem)
-import Control.Arrow
 import Control.Exception
 import Data.Euclidean (GcdDomain(..))
-import Data.Finite
 import Data.Function
 import Data.Int
 import Data.List (groupBy, sortOn)
 import Data.Mod
-import Data.Poly.Sparse.Multi
 import Data.Proxy
 import Data.Semiring (Semiring(..))
 import qualified Data.Vector as V
@@ -33,34 +28,11 @@ import qualified Data.Vector.Unboxed.Sized as SU
 import Test.Tasty
 import Test.Tasty.QuickCheck hiding (scale, numTests)
 
-#if MIN_VERSION_base(4,10,0)
-import GHC.TypeNats (KnownNat)
-#else
-import GHC.TypeLits (KnownNat)
-#endif
+import Data.Poly.Multi.Semiring
 
 import Quaternion
-import Sparse ()
 import TestUtils
 
-instance KnownNat n => Arbitrary (Finite n) where
-  arbitrary = elements finites
-
-instance (Arbitrary a, KnownNat n, G.Vector v a) => Arbitrary (SG.Vector v n a) where
-  arbitrary = SG.replicateM arbitrary
-  shrink vs = [ vs SG.// [(i, x)] | i <- finites, x <- shrink (SG.index vs i) ]
-
-instance (Eq a, Semiring a, Arbitrary a, KnownNat n, G.Vector v (SU.Vector n Word, a)) => Arbitrary (MultiPoly v n a) where
-  arbitrary = toMultiPoly . G.fromList <$> arbitrary
-  shrink = fmap (toMultiPoly . G.fromList) . shrink . G.toList . unMultiPoly
-
-newtype ShortMultiPoly a = ShortMultiPoly { unShortMultiPoly :: a }
-  deriving (Eq, Show, Semiring, GcdDomain)
-
-instance (Eq a, Semiring a, Arbitrary a, KnownNat n, G.Vector v (SU.Vector n Word, a)) => Arbitrary (ShortMultiPoly (MultiPoly v n a)) where
-  arbitrary = ShortMultiPoly . toMultiPoly . G.fromList . (\xs -> take (length xs `mod` 4) (map (first (SU.map (`mod` 3))) xs)) <$> arbitrary
-  shrink = fmap (ShortMultiPoly . toMultiPoly . G.fromList) . shrink . G.toList . unMultiPoly . unShortMultiPoly
-
 testSuite :: TestTree
 testSuite = testGroup "Multi"
   [ arithmeticTests
@@ -104,11 +76,11 @@ numTests =
 
 gcdDomainTests :: [TestTree]
 gcdDomainTests =
-  [ myGcdDomainLaws (Proxy :: Proxy (ShortMultiPoly (VMultiPoly 3 Integer)))
+  [ myGcdDomainLaws (Proxy :: Proxy (ShortPoly (VMultiPoly 3 Integer)))
   , tenTimesLess
-  $ myGcdDomainLaws (Proxy :: Proxy (ShortMultiPoly (VMultiPoly 3 (Mod 3))))
+  $ myGcdDomainLaws (Proxy :: Proxy (ShortPoly (VMultiPoly 3 (Mod 3))))
   , tenTimesLess
-  $ myGcdDomainLaws (Proxy :: Proxy (ShortMultiPoly (VMultiPoly 3 Rational)))
+  $ myGcdDomainLaws (Proxy :: Proxy (ShortPoly (VMultiPoly 3 Rational)))
   ]
 
 isListTests :: [TestTree]
@@ -296,8 +268,4 @@ conversionTests = testGroup "conversions"
     \(xs :: UMultiPoly 3 Int8) -> xs === unsegregate (segregate xs)
   , testProperty "segregate . unsegregate = id" $
     \xs -> xs === segregate (unsegregate xs :: UMultiPoly 3 Int8)
-  , testProperty "multiPolyToPoly . polyToMultiPoly" $
-    \xs -> xs === multiPolyToPoly (polyToMultiPoly xs :: UMultiPoly 1 Int8)
-  , testProperty "polyToMultiPoly . multiPolyToPoly" $
-    \(xs :: UMultiPoly 1 Int8) -> xs === polyToMultiPoly (multiPolyToPoly xs)
   ]
diff --git a/test/Sparse.hs b/test/Sparse.hs
index f89a097..98f7742 100644
--- a/test/Sparse.hs
+++ b/test/Sparse.hs
@@ -5,8 +5,6 @@
 {-# LANGUAGE ScopedTypeVariables        #-}
 {-# LANGUAGE UndecidableInstances       #-}
 
-{-# OPTIONS_GHC -fno-warn-orphans #-}
-
 module Sparse
   ( testSuite
   , ShortPoly(..)
@@ -26,23 +24,13 @@ import Data.Semiring (Semiring(..))
 import qualified Data.Vector as V
 import qualified Data.Vector.Generic as G
 import qualified Data.Vector.Unboxed as U
+import qualified Data.Vector.Unboxed.Sized as SU
 import Test.Tasty
 import Test.Tasty.QuickCheck hiding (scale, numTests)
 
 import Quaternion
 import TestUtils
 
-instance (Eq a, Semiring a, Arbitrary a, G.Vector v (Word, a)) => Arbitrary (Poly v a) where
-  arbitrary = S.toPoly . G.fromList <$> arbitrary
-  shrink = fmap (S.toPoly . G.fromList) . shrink . G.toList . unPoly
-
-newtype ShortPoly a = ShortPoly { unShortPoly :: a }
-  deriving (Eq, Show, Semiring, GcdDomain, Euclidean)
-
-instance (Eq a, Semiring a, Arbitrary a, G.Vector v (Word, a)) => Arbitrary (ShortPoly (Poly v a)) where
-  arbitrary = ShortPoly . S.toPoly . G.fromList . (\xs -> take (length xs `mod` 5) xs) <$> arbitrary
-  shrink = fmap (ShortPoly . S.toPoly . G.fromList) . shrink . G.toList . unPoly . unShortPoly
-
 testSuite :: TestTree
 testSuite = testGroup "Sparse"
   [ arithmeticTests
@@ -205,7 +193,7 @@ evalTests = testGroup "eval" $ concat
 
 evalTestGroup
   :: forall v a.
-     (Eq a, Num a, Semiring a, Arbitrary a, Show a, Eq (v (Word, a)), Show (v (Word, a)), G.Vector v (Word, a))
+     (Eq a, Num a, Semiring a, Arbitrary a, Show a, Eq (v (Word, a)), Show (v (Word, a)), G.Vector v (Word, a), G.Vector v (SU.Vector 1 Word, a))
   => Proxy (Poly v a)
   -> [TestTree]
 evalTestGroup _ =
@@ -236,7 +224,7 @@ evalTestGroup _ =
 
 substTestGroup
   :: forall v a.
-     (Eq a, Num a, Semiring a, Arbitrary a, Show a, Eq (v (Word, a)), Show (v (Word, a)), G.Vector v (Word, a))
+     (Eq a, Num a, Semiring a, Arbitrary a, Show a, Eq (v (SU.Vector 1 Word, a)), Show (v (SU.Vector 1 Word, a)), G.Vector v (Word, a), G.Vector v (SU.Vector 1 Word, a))
   => Proxy (Poly v a)
   -> [TestTree]
 substTestGroup _ =
diff --git a/test/SparseLaurent.hs b/test/SparseLaurent.hs
index e681f94..06ed08b 100644
--- a/test/SparseLaurent.hs
+++ b/test/SparseLaurent.hs
@@ -1,11 +1,10 @@
+{-# LANGUAGE DataKinds                  #-}
 {-# LANGUAGE FlexibleContexts           #-}
 {-# LANGUAGE FlexibleInstances          #-}
 {-# LANGUAGE GeneralizedNewtypeDeriving #-}
 {-# LANGUAGE ScopedTypeVariables        #-}
 {-# LANGUAGE UndecidableInstances       #-}
 
-{-# OPTIONS_GHC -fno-warn-orphans #-}
-
 module SparseLaurent
   ( testSuite
   ) where
@@ -13,7 +12,7 @@ module SparseLaurent
 import Prelude hiding (gcd, quotRem, quot, rem)
 import Control.Exception
 import Data.Either
-import Data.Euclidean (Euclidean(..), GcdDomain(..), Field)
+import Data.Euclidean (GcdDomain(..), Field)
 import Data.Int
 import qualified Data.Poly.Sparse
 import Data.Poly.Sparse.Laurent
@@ -22,24 +21,13 @@ import Data.Semiring (Semiring(..))
 import qualified Data.Vector as V
 import qualified Data.Vector.Generic as G
 import qualified Data.Vector.Unboxed as U
+import qualified Data.Vector.Unboxed.Sized as SU
 import Test.Tasty
 import Test.Tasty.QuickCheck hiding (scale, numTests)
 
 import Quaternion
-import Sparse (ShortPoly(..))
 import TestUtils
 
-instance (Eq a, Semiring a, Arbitrary a, G.Vector v (Word, a)) => Arbitrary (Laurent v a) where
-  arbitrary = toLaurent <$> ((`rem` 10) <$> arbitrary) <*> arbitrary
-  shrink = fmap (uncurry toLaurent) . shrink . unLaurent
-
-newtype ShortLaurent a = ShortLaurent { unShortLaurent :: a }
-  deriving (Eq, Show, Semiring, GcdDomain)
-
-instance (Eq a, Semiring a, Arbitrary a, G.Vector v (Word, a)) => Arbitrary (ShortLaurent (Laurent v a)) where
-  arbitrary = (ShortLaurent .) . toLaurent <$> ((`rem` 10) <$> arbitrary) <*> (unShortPoly <$> arbitrary)
-  shrink = fmap (ShortLaurent . uncurry toLaurent . fmap unShortPoly) . shrink . fmap ShortPoly . unLaurent . unShortLaurent
-
 testSuite :: TestTree
 testSuite = testGroup "SparseLaurent"
   [ otherTests
@@ -81,9 +69,9 @@ numTests =
 
 gcdDomainTests :: [TestTree]
 gcdDomainTests =
-  [ myGcdDomainLaws (Proxy :: Proxy (ShortLaurent (Laurent V.Vector Integer)))
+  [ myGcdDomainLaws (Proxy :: Proxy (ShortPoly (Laurent V.Vector Integer)))
   , tenTimesLess
-  $ myGcdDomainLaws (Proxy :: Proxy (ShortLaurent (Laurent V.Vector Rational)))
+  $ myGcdDomainLaws (Proxy :: Proxy (ShortPoly (Laurent V.Vector Rational)))
   ]
 
 isListTests :: [TestTree]
@@ -143,7 +131,7 @@ evalTests = testGroup "eval" $ concat
 
 evalTestGroup
   :: forall v a.
-     (Eq a, Field a, Arbitrary a, Show a, Eq (v (Word, a)), Show (v (Word, a)), G.Vector v (Word, a))
+     (Eq a, Field a, Arbitrary a, Show a, Eq (v (Word, a)), Show (v (Word, a)), G.Vector v (Word, a), G.Vector v (SU.Vector 1 Word, a))
   => Proxy (Laurent v a)
   -> [TestTree]
 evalTestGroup _ =
@@ -162,7 +150,7 @@ evalTestGroup _ =
 
 substTestGroup
   :: forall v a.
-     (Eq a, Num a, Semiring a, Arbitrary a, Show a, Eq (v (Word, a)), Show (v (Word, a)), G.Vector v (Word, a))
+     (Eq a, Num a, Semiring a, Arbitrary a, Show a, Eq (v (SU.Vector 1 Word, a)), Show (v (Word, a)), G.Vector v (Word, a), G.Vector v (SU.Vector 1 Word, a))
   => Proxy (Laurent v a)
   -> [TestTree]
 substTestGroup _ =
diff --git a/test/TestUtils.hs b/test/TestUtils.hs
index 67a9aca..14378af 100644
--- a/test/TestUtils.hs
+++ b/test/TestUtils.hs
@@ -1,10 +1,15 @@
-{-# LANGUAGE CPP              #-}
-{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE CPP                        #-}
+{-# LANGUAGE DataKinds                  #-}
+{-# LANGUAGE FlexibleContexts           #-}
+{-# LANGUAGE FlexibleInstances          #-}
+{-# LANGUAGE GeneralizedNewtypeDeriving #-}
+{-# LANGUAGE UndecidableInstances       #-}
 
 {-# OPTIONS_GHC -fno-warn-orphans #-}
 
 module TestUtils
-  ( tenTimesLess
+  ( ShortPoly(..)
+  , tenTimesLess
   , mySemiringLaws
   , myRingLaws
   , myNumLaws
@@ -14,10 +19,15 @@ module TestUtils
   , myShowLaws
   ) where
 
-import Data.Euclidean
+import Control.Arrow
+import Data.Euclidean hiding (rem)
+import Data.Finite
 import Data.Mod
 import Data.Proxy
 import Data.Semiring (Semiring, Ring)
+import qualified Data.Vector.Generic as G
+import qualified Data.Vector.Generic.Sized as SG
+import qualified Data.Vector.Unboxed.Sized as SU
 import GHC.Exts
 import Test.QuickCheck.Classes
 import Test.Tasty
@@ -29,10 +39,59 @@ import GHC.TypeNats (KnownNat)
 import GHC.TypeLits (KnownNat)
 #endif
 
+import qualified Data.Poly.Semiring as Dense
+import qualified Data.Poly.Laurent as DenseLaurent
+import Data.Poly.Multi.Semiring
+import qualified Data.Poly.Sparse.Laurent as SparseLaurent
+
+newtype ShortPoly a = ShortPoly { unShortPoly :: a }
+  deriving (Eq, Show, Semiring, GcdDomain, Euclidean)
+
 instance KnownNat m => Arbitrary (Mod m) where
   arbitrary = oneof [arbitraryBoundedEnum, fromInteger <$> arbitrary]
   shrink = map fromInteger . shrink . toInteger . unMod
 
+instance KnownNat n => Arbitrary (Finite n) where
+  arbitrary = elements finites
+
+instance (Arbitrary a, KnownNat n, G.Vector v a) => Arbitrary (SG.Vector v n a) where
+  arbitrary = SG.replicateM arbitrary
+  shrink vs = [ vs SG.// [(i, x)] | i <- finites, x <- shrink (SG.index vs i) ]
+
+instance (Eq a, Semiring a, Arbitrary a, G.Vector v a) => Arbitrary (Dense.Poly v a) where
+  arbitrary = Dense.toPoly . G.fromList <$> arbitrary
+  shrink = fmap (Dense.toPoly . G.fromList) . shrink . G.toList . Dense.unPoly
+
+instance (Eq a, Semiring a, Arbitrary a, G.Vector v a) => Arbitrary (ShortPoly (Dense.Poly v a)) where
+  arbitrary = ShortPoly . Dense.toPoly . G.fromList . (\xs -> take (length xs `mod` 10) xs) <$> arbitrary
+  shrink = fmap (ShortPoly . Dense.toPoly . G.fromList) . shrink . G.toList . Dense.unPoly . unShortPoly
+
+instance (Eq a, Semiring a, Arbitrary a, G.Vector v a) => Arbitrary (DenseLaurent.Laurent v a) where
+  arbitrary = DenseLaurent.toLaurent <$> ((`rem` 10) <$> arbitrary) <*> arbitrary
+  shrink = fmap (uncurry DenseLaurent.toLaurent) . shrink . DenseLaurent.unLaurent
+
+instance (Eq a, Semiring a, Arbitrary a, G.Vector v a) => Arbitrary (ShortPoly (DenseLaurent.Laurent v a)) where
+  arbitrary = (ShortPoly .) . DenseLaurent.toLaurent <$> ((`rem` 10) <$> arbitrary) <*> (unShortPoly <$> arbitrary)
+  shrink = fmap (ShortPoly . uncurry DenseLaurent.toLaurent . fmap unShortPoly) . shrink . fmap ShortPoly . DenseLaurent.unLaurent . unShortPoly
+
+instance (Eq a, Semiring a, Arbitrary a, KnownNat n, G.Vector v (SU.Vector n Word, a)) => Arbitrary (MultiPoly v n a) where
+  arbitrary = toMultiPoly . G.fromList <$> arbitrary
+  shrink = fmap (toMultiPoly . G.fromList) . shrink . G.toList . unMultiPoly
+
+instance (Eq a, Semiring a, Arbitrary a, KnownNat n, G.Vector v (SU.Vector n Word, a)) => Arbitrary (ShortPoly (MultiPoly v n a)) where
+  arbitrary = ShortPoly . toMultiPoly . G.fromList . (\xs -> take (length xs `mod` 4) (map (first (SU.map (`mod` 3))) xs)) <$> arbitrary
+  shrink = fmap (ShortPoly . toMultiPoly . G.fromList) . shrink . G.toList . unMultiPoly . unShortPoly
+
+instance (Eq a, Semiring a, Arbitrary a, G.Vector v (Word, a), G.Vector v (SU.Vector 1 Word, a)) => Arbitrary (SparseLaurent.Laurent v a) where
+  arbitrary = SparseLaurent.toLaurent <$> ((`rem` 10) <$> arbitrary) <*> arbitrary
+  shrink = fmap (uncurry SparseLaurent.toLaurent) . shrink . SparseLaurent.unLaurent
+
+instance (Eq a, Semiring a, Arbitrary a, G.Vector v (Word, a), G.Vector v (SU.Vector 1 Word, a)) => Arbitrary (ShortPoly (SparseLaurent.Laurent v a)) where
+  arbitrary = (ShortPoly .) . SparseLaurent.toLaurent <$> ((`rem` 10) <$> arbitrary) <*> (unShortPoly <$> arbitrary)
+  shrink = fmap (ShortPoly . uncurry SparseLaurent.toLaurent . fmap unShortPoly) . shrink . fmap ShortPoly . SparseLaurent.unLaurent . unShortPoly
+
+
+-------------------------------------------------------------------------------
 tenTimesLess :: TestTree -> TestTree
 tenTimesLess = adjustOption $
   \(QuickCheckTests n) -> QuickCheckTests (max 100 (n `div` 10))