diff --git a/bench/DenseBench.hs b/bench/DenseBench.hs
--- a/bench/DenseBench.hs
+++ b/bench/DenseBench.hs
@@ -8,17 +8,14 @@
   ) where
 
 import Prelude hiding (quotRem, gcd)
-import Control.DeepSeq
 import Gauge.Main
 import Data.Poly
 import qualified Data.Vector.Unboxed as U
-#if MIN_VERSION_semirings(0,4,2)
-import Control.Exception
-import Data.Bits
-import Data.Coerce
-import Data.Euclidean
-import Data.Semiring (Semiring(..), Ring, isZero)
-import qualified Data.Semiring as S (negate)
+#if MIN_VERSION_semirings(0,5,2)
+import Data.Euclidean (Euclidean(..), GcdDomain(..), Field)
+import qualified Data.Poly.Semiring as S (toPoly)
+import Data.Semiring (Semiring(..), Ring, Mod2(..))
+import qualified Data.Semiring as S (fromIntegral)
 import qualified Data.Vector as V
 #endif
 
@@ -29,7 +26,7 @@
   , map benchEval     [100, 1000, 10000]
   , map benchDeriv    [100, 1000, 10000]
   , map benchIntegral [100, 1000, 10000]
-#if MIN_VERSION_semirings(0,4,2)
+#if MIN_VERSION_semirings(0,5,2)
   , map benchQuotRem    [10, 100]
   , map benchGcd        [10, 100]
   , map benchGcdExtRat  [10, 20, 40]
@@ -54,7 +51,7 @@
 benchIntegral :: Int -> Benchmark
 benchIntegral k = bench ("integral/" ++ show k) $ nf doIntegral k
 
-#if MIN_VERSION_semirings(0,4,2)
+#if MIN_VERSION_semirings(0,5,2)
 
 benchQuotRem :: Int -> Benchmark
 benchQuotRem k = bench ("quotRem/" ++ show k) $ nf doQuotRem k
@@ -69,10 +66,10 @@
 benchGcdFracRat k = bench ("gcdFrac/Rational/" ++ show k) $ nf (doGcdFrac @Rational) k
 
 benchGcdExtM :: Int -> Benchmark
-benchGcdExtM k = bench ("gcdExt/Mod2/" ++ show k) $ nf (doGcdExt @Mod2) k
+benchGcdExtM k = bench ("gcdExt/Mod2/" ++ show k) $ nf (getMod2 . doGcdExt @Mod2) k
 
 benchGcdFracM :: Int -> Benchmark
-benchGcdFracM k = bench ("gcdFrac/Mod2/" ++ show k) $ nf (doGcdFrac @Mod2) k
+benchGcdFracM k = bench ("gcdFrac/Mod2/" ++ show k) $ nf (getMod2 . doGcdFrac @Mod2) k
 
 #endif
 
@@ -101,13 +98,13 @@
     xs = toPoly $ U.generate n ((* 2) . fromIntegral)
     zs = unPoly $ integral xs
 
-#if MIN_VERSION_semirings(0,4,2)
+#if MIN_VERSION_semirings(0,5,2)
 
-gen1 :: Num a => Int -> a
-gen1 k = fromIntegral (truncate (pi * fromIntegral k :: Double) `mod` (k + 1))
+gen1 :: Ring a => Int -> a
+gen1 k = S.fromIntegral (truncate (pi * fromIntegral k :: Double) `mod` (k + 1))
 
-gen2 :: Num a => Int -> a
-gen2 k = fromIntegral (truncate (exp 1.0 * fromIntegral k :: Double) `mod` (k + 1))
+gen2 :: Ring a => Int -> a
+gen2 k = S.fromIntegral (truncate (exp 1.0 * fromIntegral k :: Double) `mod` (k + 1))
 
 doQuotRem :: Int -> Double
 doQuotRem n = U.sum (unPoly qs) + U.sum (unPoly rs)
@@ -123,48 +120,18 @@
     ys = toPoly $ V.generate n gen2
     gs = unPoly $ xs `gcd` ys
 
-doGcdExt :: (Eq a, Num a, Field a) => Int -> a
-doGcdExt n = V.sum gs
+doGcdExt :: (Eq a, Field a) => Int -> a
+doGcdExt n = V.foldl' plus zero gs
   where
-    xs = toPoly $ V.generate n gen1
-    ys = toPoly $ V.generate n gen2
+    xs = S.toPoly $ V.generate n gen1
+    ys = S.toPoly $ V.generate n gen2
     gs = unPoly $ fst $ xs `gcdExt` ys
 
-doGcdFrac :: (Eq a, Num a, Field a) => Int -> a
-doGcdFrac n = V.sum gs
+doGcdFrac :: (Eq a, Field a) => Int -> a
+doGcdFrac n = V.foldl' plus zero gs
   where
-    xs = PolyOverField $ toPoly $ V.generate n gen1
-    ys = PolyOverField $ toPoly $ V.generate n gen2
+    xs = PolyOverField $ S.toPoly $ V.generate n gen1
+    ys = PolyOverField $ S.toPoly $ V.generate n gen2
     gs = unPoly $ unPolyOverField $ xs `gcd` ys
-
--- | Inspired by 'semirings'.
-newtype Mod2 = Mod2 { _getMod2 :: Bool }
-  deriving (Eq, NFData)
-
-instance Num Mod2 where
-  (+) = coerce (xor @Bool)
-  (*) = coerce (&&)
-  negate = id
-  abs    = id
-  signum = id
-  fromInteger = Mod2 . odd
-
-instance Semiring Mod2 where
-  plus  = coerce (xor @Bool)
-  times = coerce (&&)
-  fromNatural = Mod2 . odd
-
-instance Ring Mod2 where
-  negate = id
-
-instance GcdDomain Mod2 where
-
-instance Euclidean Mod2 where
-  degree = const 0
-  quotRem x y
-    | isZero y  = throw DivideByZero
-    | otherwise = (x, zero)
-
-instance Field Mod2
 
 #endif
diff --git a/changelog.md b/changelog.md
--- a/changelog.md
+++ b/changelog.md
@@ -1,3 +1,8 @@
+# 0.3.3.0
+
+* Add function `subst`.
+* Fix compatibility issues.
+
 # 0.3.2.0
 
 * Add `NFData` instance.
diff --git a/poly.cabal b/poly.cabal
--- a/poly.cabal
+++ b/poly.cabal
@@ -1,5 +1,5 @@
 name: poly
-version: 0.3.2.0
+version: 0.3.3.0
 synopsis: Polynomials
 description:
   Polynomials backed by `Vector`.
diff --git a/src/Data/Poly.hs b/src/Data/Poly.hs
--- a/src/Data/Poly.hs
+++ b/src/Data/Poly.hs
@@ -22,6 +22,7 @@
   , scale
   , pattern X
   , eval
+  , subst
   , deriv
   , integral
 #if MIN_VERSION_semirings(0,4,2)
diff --git a/src/Data/Poly/Internal/Dense.hs b/src/Data/Poly/Internal/Dense.hs
--- a/src/Data/Poly/Internal/Dense.hs
+++ b/src/Data/Poly/Internal/Dense.hs
@@ -27,6 +27,7 @@
   , scale
   , pattern X
   , eval
+  , subst
   , deriv
   , integral
   -- * Semiring interface
@@ -35,8 +36,10 @@
   , scale'
   , pattern X'
   , eval'
+  , subst'
   , deriv'
 #if MIN_VERSION_semirings(0,5,0)
+  , unscale'
   , integral'
 #endif
   ) where
@@ -60,7 +63,7 @@
 import Numeric.Natural
 #endif
 #if MIN_VERSION_semirings(0,5,0)
-import Data.Euclidean (Field, quot)
+import Data.Euclidean (Euclidean, Field, quot)
 #endif
 
 -- | Polynomials of one variable with coefficients from @a@,
@@ -223,7 +226,7 @@
     (G.unsafeSlice  lenMn (lenMx - lenMn) (if lenXs <= lenYs then ys else xs))
 
   G.unsafeFreeze zs
-{-# INLINE plusPoly #-}
+{-# INLINABLE plusPoly #-}
 
 minusPoly
   :: G.Vector v a
@@ -250,7 +253,7 @@
       (G.unsafeSlice  lenYs (lenXs - lenYs) xs)
 
   G.unsafeFreeze zs
-{-# INLINE minusPoly #-}
+{-# INLINABLE minusPoly #-}
 
 karatsubaThreshold :: Int
 karatsubaThreshold = 32
@@ -300,7 +303,7 @@
     zs0  = karatsuba xs0 ys0
     zs2  = karatsuba xs1 ys1
     zs11 = karatsuba xs01 ys01
-{-# INLINE karatsuba #-}
+{-# INLINABLE karatsuba #-}
 
 convolution
   :: G.Vector v a
@@ -320,7 +323,7 @@
     lenXs = G.length xs
     lenYs = G.length ys
     lenZs = lenXs + lenYs - 1
-{-# INLINE convolution #-}
+{-# INLINABLE convolution #-}
 
 -- | Create a monomial from a power and a coefficient.
 monomial :: (Eq a, Num a, G.Vector v a) => Word -> a -> Poly v a
@@ -350,7 +353,7 @@
   forM_ [0 .. lenXs - 1] $ \k ->
     MG.unsafeWrite zs (fromIntegral yp + k) (mul yc $ G.unsafeIndex xs k)
   G.unsafeFreeze zs
-{-# INLINE scaleInternal #-}
+{-# INLINABLE scaleInternal #-}
 
 -- | Multiply a polynomial by a monomial, expressed as a power and a coefficient.
 --
@@ -362,6 +365,17 @@
 scale' :: (Eq a, Semiring a, G.Vector v a) => Word -> a -> Poly v a -> Poly v a
 scale' yp yc (Poly xs) = toPoly' $ scaleInternal zero times yp yc xs
 
+#if MIN_VERSION_semirings(0,5,0)
+unscale' :: (Eq a, Euclidean a, G.Vector v a) => Word -> a -> Poly v a -> Poly v a
+unscale' yp yc (Poly xs) = toPoly' $ runST $ do
+  let lenZs = G.length xs - fromIntegral yp
+  zs <- MG.unsafeNew lenZs
+  forM_ [0 .. lenZs - 1] $ \k ->
+    MG.unsafeWrite zs k (G.unsafeIndex xs (k + fromIntegral yp) `quot` yc)
+  G.unsafeFreeze zs
+{-# INLINABLE unscale' #-}
+#endif
+
 data StrictPair a b = !a :*: !b
 
 infixr 1 :*:
@@ -373,18 +387,36 @@
 --
 -- >>> eval (X^2 + 1 :: UPoly Int) 3
 -- 10
--- >>> eval (X^2 + 1 :: VPoly (UPoly Int)) (X + 1)
--- 1 * X^2 + 2 * X + 2
 eval :: (Num a, G.Vector v a) => Poly v a -> a -> a
-eval (Poly cs) x = fst' $
-  G.foldl' (\(acc :*: xn) cn -> acc + cn * xn :*: x * xn) (0 :*: 1) cs
+eval = substitute (*)
 {-# INLINE eval #-}
 
 eval' :: (Semiring a, G.Vector v a) => Poly v a -> a -> a
-eval' (Poly cs) x = fst' $
-  G.foldl' (\(acc :*: xn) cn -> acc `plus` cn `times` xn :*: x `times` xn) (zero :*: one) cs
+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 a, G.Vector w a) => Poly v a -> Poly w a -> Poly w a
+subst = substitute (scale 0)
+{-# INLINE subst #-}
+
+subst' :: (Eq a, Semiring a, G.Vector v a, G.Vector w a) => Poly v a -> Poly w a -> Poly w a
+subst' = substitute' (scale' 0)
+{-# INLINE subst' #-}
+
+substitute :: (G.Vector v a, Num b) => (a -> b -> b) -> Poly v a -> b -> b
+substitute f (Poly cs) x = fst' $
+  G.foldl' (\(acc :*: xn) cn -> acc + f cn xn :*: x * xn) (0 :*: 1) cs
+{-# INLINE substitute #-}
+
+substitute' :: (G.Vector v a, Semiring b) => (a -> b -> b) -> Poly v a -> b -> b
+substitute' f (Poly cs) x = fst' $
+  G.foldl' (\(acc :*: xn) cn -> acc `plus` f cn xn :*: x `times` xn) (zero :*: one) cs
+{-# INLINE substitute' #-}
+
 -- | Take a derivative.
 --
 -- >>> deriv (X^3 + 3 * X) :: UPoly Int
@@ -428,7 +460,7 @@
     G.unsafeFreeze zs
     where
       lenXs = G.length xs
-{-# INLINE integral #-}
+{-# INLINABLE integral #-}
 
 #if MIN_VERSION_semirings(0,5,0)
 integral' :: (Eq a, Field a, G.Vector v a) => Poly v a -> Poly v a
@@ -442,7 +474,7 @@
     G.unsafeFreeze zs
     where
       lenXs = G.length xs
-{-# INLINE integral' #-}
+{-# INLINABLE integral' #-}
 #endif
 
 -- | Create an identity polynomial.
diff --git a/src/Data/Poly/Internal/Dense/Field.hs b/src/Data/Poly/Internal/Dense/Field.hs
--- a/src/Data/Poly/Internal/Dense/Field.hs
+++ b/src/Data/Poly/Internal/Dense/Field.hs
@@ -28,9 +28,11 @@
 import Control.Monad
 import Control.Monad.Primitive
 import Control.Monad.ST
-import Data.Euclidean
+import Data.Euclidean (Euclidean(..))
 #if !MIN_VERSION_semirings(0,5,0)
 import Data.Semiring (Ring)
+#else
+import Data.Euclidean (Field)
 #endif
 import Data.Semiring (times, minus, zero, one)
 import qualified Data.Vector.Generic as G
@@ -77,7 +79,7 @@
         MG.unsafeModify rs (\c -> c `minus` q `times` G.unsafeIndex ys k) (i + k)
     let rs' = MG.unsafeSlice 0 lenYs rs
     (,) <$> G.unsafeFreeze qs <*> G.unsafeFreeze rs'
-{-# INLINE quotientAndRemainder #-}
+{-# INLINABLE quotientAndRemainder #-}
 
 remainder
   :: (Field a, G.Vector v a)
@@ -91,7 +93,7 @@
     ys' <- G.unsafeThaw ys
     remainderM rs ys'
     G.unsafeFreeze $ MG.unsafeSlice 0 (G.length xs `min` G.length ys) rs
-{-# INLINE remainder #-}
+{-# INLINABLE remainder #-}
 
 remainderM
   :: (PrimMonad m, Field a, G.Vector v a)
@@ -113,7 +115,7 @@
         -- do not move r / yLast outside the loop,
         -- because of numerical instability
         MG.unsafeModify xs (\c -> c `minus` r `times` y `quot` yLast) (i + k)
-{-# INLINE remainderM #-}
+{-# INLINABLE remainderM #-}
 
 fieldGcd
   :: (Eq a, Field a, G.Vector v a)
@@ -164,7 +166,7 @@
           | r' == zero = (r, s)
           | otherwise  = case r `quotRem` r' of
             (q, r'') -> go r'' r' (s `minus` q `times` s') s'
-{-# INLINE gcdExt #-}
+{-# INLINABLE gcdExt #-}
 
 -- | Scale a non-zero polynomial such that its leading coefficient is one,
 -- returning the reciprocal of the leading coefficient in the scaling.
diff --git a/src/Data/Poly/Internal/Dense/GcdDomain.hs b/src/Data/Poly/Internal/Dense/GcdDomain.hs
--- a/src/Data/Poly/Internal/Dense/GcdDomain.hs
+++ b/src/Data/Poly/Internal/Dense/GcdDomain.hs
@@ -118,7 +118,7 @@
       MG.unsafeModify xs (`times` zx)
     xs' <- dropWhileEndM isZero xs
     gcdM xs' ys
-{-# INLINE gcdM #-}
+{-# INLINABLE gcdM #-}
 
 isZeroM
   :: (Eq a, Semiring a, PrimMonad m, G.Vector v a)
@@ -169,6 +169,6 @@
                 go (i - 1)
 
     go (lenQs - 1)
-{-# INLINE quotient #-}
+{-# INLINABLE quotient #-}
 
 #endif
diff --git a/src/Data/Poly/Internal/Sparse.hs b/src/Data/Poly/Internal/Sparse.hs
--- a/src/Data/Poly/Internal/Sparse.hs
+++ b/src/Data/Poly/Internal/Sparse.hs
@@ -28,6 +28,7 @@
   , scale
   , pattern X
   , eval
+  , subst
   , deriv
   , integral
   -- * Semiring interface
@@ -36,6 +37,7 @@
   , scale'
   , pattern X'
   , eval'
+  , subst'
   , deriv'
 #if MIN_VERSION_semirings(0,5,0)
   , integral'
@@ -235,7 +237,7 @@
   zs <- MG.unsafeNew (G.length xs + G.length ys)
   lenZs <- plusPolyM p add xs ys zs
   G.unsafeFreeze $ MG.unsafeSlice 0 lenZs zs
-{-# INLINE plusPoly #-}
+{-# INLINABLE plusPoly #-}
 
 plusPolyM
   :: (PrimMonad m, G.Vector v (Word, a))
@@ -278,7 +280,7 @@
         GT -> do
           MG.unsafeWrite zs iz (yp, yc)
           go ix (iy + 1) (iz + 1)
-{-# INLINE plusPolyM #-}
+{-# INLINABLE plusPolyM #-}
 
 minusPoly
   :: G.Vector v (Word, a)
@@ -323,7 +325,7 @@
   where
     lenXs = G.length xs
     lenYs = G.length ys
-{-# INLINE minusPoly #-}
+{-# INLINABLE minusPoly #-}
 
 scaleM
   :: (PrimMonad m, G.Vector v (Word, a))
@@ -347,7 +349,7 @@
           go (ix + 1) (iz + 1)
         else
           go (ix + 1) iz
-{-# INLINE scaleM #-}
+{-# INLINABLE scaleM #-}
 
 scaleInternal
   :: G.Vector v (Word, a)
@@ -361,7 +363,7 @@
   zs <- MG.unsafeNew (G.length xs)
   len <- scaleM p (flip mul) xs (yp, yc) zs
   fmap Poly $ G.unsafeFreeze $ MG.unsafeSlice 0 len zs
-{-# INLINE scaleInternal #-}
+{-# INLINABLE scaleInternal #-}
 
 -- | Multiply a polynomial by a monomial, expressed as a power and a coefficient.
 --
@@ -443,7 +445,7 @@
         buffer'    <- G.unsafeThaw   buffer
         bufferNew' <- G.unsafeFreeze bufferNew
         gogo slicesNew' bufferNew' buffer'
-{-# INLINE convolution #-}
+{-# 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
@@ -464,23 +466,49 @@
 --
 -- >>> eval (X^2 + 1 :: UPoly Int) 3
 -- 10
--- >>> eval (X^2 + 1 :: VPoly (UPoly Int)) (X + 1)
--- 1 * X^2 + 2 * X + 2
 eval :: (Num a, G.Vector v (Word, a)) => Poly v a -> a -> a
-eval (Poly cs) x = fst3 $ G.foldl' go (Strict3 0 0 1) cs
+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 + c * xp) p xp
-{-# INLINE eval #-}
+        Strict3 (acc + f c xp) p xp
+{-# INLINE substitute #-}
 
-eval' :: (Semiring a, G.Vector v (Word, a)) => Poly v a -> a -> a
-eval' (Poly cs) x = fst3 $ G.foldl' go (Strict3 zero 0 one) cs
+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` c `times` xp) p xp
-{-# INLINE eval' #-}
+        Strict3 (acc `plus` f c xp) p xp
+{-# INLINE substitute' #-}
 
 -- | Take a derivative.
 --
@@ -534,7 +562,7 @@
               go (ix + 1) iz
     lenZs <- go 0 0
     G.unsafeFreeze $ MG.unsafeSlice 0 lenZs zs
-{-# INLINE derivPoly #-}
+{-# INLINABLE derivPoly #-}
 
 -- | Compute an indefinite integral of a polynomial,
 -- setting constant term to zero.
diff --git a/src/Data/Poly/Internal/Sparse/Field.hs b/src/Data/Poly/Internal/Sparse/Field.hs
--- a/src/Data/Poly/Internal/Sparse/Field.hs
+++ b/src/Data/Poly/Internal/Sparse/Field.hs
@@ -27,9 +27,11 @@
 import Prelude hiding (quotRem, quot, rem, gcd)
 import Control.Arrow
 import Control.Exception
-import Data.Euclidean
+import Data.Euclidean (Euclidean(..))
 #if !MIN_VERSION_semirings(0,5,0)
 import Data.Semiring (Ring)
+#else
+import Data.Euclidean (Field)
 #endif
 import Data.Semiring (minus, plus, times, zero, one)
 import qualified Data.Vector.Generic as G
@@ -93,7 +95,7 @@
           | r' == zero = (r, s)
           | otherwise  = case r `quotRem` r' of
             (q, r'') -> go r'' r' (s `minus` q `times` s') s'
-{-# INLINE gcdExt #-}
+{-# INLINABLE gcdExt #-}
 
 -- | Scale a non-zero polynomial such that its leading coefficient is one,
 -- returning the reciprocal of the leading coefficient in the scaling.
diff --git a/src/Data/Poly/Semiring.hs b/src/Data/Poly/Semiring.hs
--- a/src/Data/Poly/Semiring.hs
+++ b/src/Data/Poly/Semiring.hs
@@ -22,6 +22,7 @@
   , scale
   , pattern X
   , eval
+  , subst
   , deriv
 #if MIN_VERSION_semirings(0,5,0)
   , integral
@@ -80,10 +81,15 @@
 --
 -- >>> eval (X^2 + 1 :: UPoly Int) 3
 -- 10
--- >>> eval (X^2 + 1 :: VPoly (UPoly Int)) (X + 1)
--- 1 * X^2 + 2 * X + 2
 eval :: (Semiring a, G.Vector v a) => Poly v a -> a -> a
 eval = Dense.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, Semiring a, G.Vector v a, G.Vector w a) => Poly v a -> Poly w a -> Poly w a
+subst = Dense.subst'
 
 -- | Take a derivative.
 --
diff --git a/src/Data/Poly/Sparse.hs b/src/Data/Poly/Sparse.hs
--- a/src/Data/Poly/Sparse.hs
+++ b/src/Data/Poly/Sparse.hs
@@ -22,6 +22,7 @@
   , scale
   , pattern X
   , eval
+  , subst
   , deriv
   , integral
 #if MIN_VERSION_semirings(0,4,2)
diff --git a/src/Data/Poly/Sparse/Semiring.hs b/src/Data/Poly/Sparse/Semiring.hs
--- a/src/Data/Poly/Sparse/Semiring.hs
+++ b/src/Data/Poly/Sparse/Semiring.hs
@@ -23,6 +23,7 @@
   , scale
   , pattern X
   , eval
+  , subst
   , deriv
 #if MIN_VERSION_semirings(0,5,0)
   , integral
@@ -76,10 +77,15 @@
 --
 -- >>> eval (X^2 + 1 :: UPoly Int) 3
 -- 10
--- >>> eval (X^2 + 1 :: VPoly (UPoly Int)) (X + 1)
--- 1 * X^2 + 2 * X + 2
 eval :: (Semiring a, G.Vector v (Word, a)) => Poly v a -> a -> a
 eval = Sparse.eval'
+
+-- | 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'
 
 -- | Take a derivative.
 --
diff --git a/test/Dense.hs b/test/Dense.hs
--- a/test/Dense.hs
+++ b/test/Dense.hs
@@ -12,7 +12,7 @@
 
 import Prelude hiding (gcd, quotRem, rem)
 #if MIN_VERSION_semirings(0,4,2)
-import Data.Euclidean
+import Data.Euclidean (Euclidean(..), GcdDomain(..))
 #endif
 import Data.Int
 import Data.Maybe
@@ -208,8 +208,9 @@
 
 evalTests :: TestTree
 evalTests = testGroup "eval" $ concat
-  [ evalTestGroup (Proxy :: Proxy (Poly U.Vector Int8))
-  , evalTestGroup (Proxy :: Proxy (Poly V.Vector Integer))
+  [ evalTestGroup  (Proxy :: Proxy (Poly U.Vector Int8))
+  , evalTestGroup  (Proxy :: Proxy (Poly V.Vector Integer))
+  , substTestGroup (Proxy :: Proxy (Poly U.Vector Int8))
   ]
 
 evalTestGroup
@@ -243,12 +244,39 @@
     e' :: Poly v a -> a -> a
     e' = S.eval
 
+substTestGroup
+  :: forall v a.
+     (Eq a, Num a, Semiring a, Arbitrary a, Show a, Eq (v a), Show (v a), G.Vector v a)
+  => Proxy (Poly v a)
+  -> [TestTree]
+substTestGroup _ =
+  [ testProperty "subst (p + q) r = subst p r + subst q r" $
+    \p q r -> e (p + q) r === e p r + e q r
+  , testProperty "subst x p = p" $
+    \p -> e X p === p
+  , testProperty "subst (monomial 0 c) p = monomial 0 c" $
+    \c p -> e (monomial 0 c) p === monomial 0 c
+  , testProperty "subst' (p + q) r = subst' p r + subst' q r" $
+    \p q r -> e' (p + q) r === e' p r + e' q r
+  , testProperty "subst' x p = p" $
+    \p -> e' S.X p === p
+  , testProperty "subst' (S.monomial 0 c) p = S.monomial 0 c" $
+    \c p -> e' (S.monomial 0 c) p === S.monomial 0 c
+  ]
+  where
+    e :: Poly v a -> Poly v a -> Poly v a
+    e = subst
+    e' :: Poly v a -> Poly v a -> Poly v a
+    e' = S.subst
+
 derivTests :: TestTree
 derivTests = testGroup "deriv"
   [ testProperty "deriv = S.deriv" $
     \(p :: Poly V.Vector Integer) -> deriv p === S.deriv p
+#if MIN_VERSION_semirings(0,5,0)
   , testProperty "integral = S.integral" $
     \(p :: Poly V.Vector Rational) -> integral p === S.integral p
+#endif
   , testProperty "deriv . integral = id" $
     \(p :: Poly V.Vector Rational) -> deriv (integral p) === p
   , testProperty "deriv c = 0" $
@@ -259,11 +287,9 @@
     \p q -> deriv (p + q) === (deriv p + deriv q :: Poly V.Vector Int)
   , testProperty "deriv (p * q) = p * deriv q + q * deriv p" $
     \p q -> deriv (p * q) === (p * deriv q + q * deriv p :: Poly V.Vector Int)
-  -- The following property takes too long for a regular test-suite
-  -- , testProperty "deriv (eval p q) = deriv q * eval (deriv p) q" $
-  --   \(p :: Poly V.Vector Int) (q :: Poly U.Vector Int) ->
-  --     deriv (eval (toPoly $ fmap (monomial 0) $ unPoly p) q) ===
-  --       deriv q * eval (toPoly $ fmap (monomial 0) $ unPoly $ deriv p) q
+  , testProperty "deriv (subst p q) = deriv q * subst (deriv p) q" $
+    \(p :: Poly V.Vector Int) (q :: Poly U.Vector Int) ->
+      deriv (subst p q) === deriv q * subst (deriv p) q
   ]
 
 #if MIN_VERSION_semirings(0,4,2)
diff --git a/test/Sparse.hs b/test/Sparse.hs
--- a/test/Sparse.hs
+++ b/test/Sparse.hs
@@ -13,7 +13,7 @@
 
 import Prelude hiding (gcd, quotRem, rem)
 #if MIN_VERSION_semirings(0,4,2)
-import Data.Euclidean
+import Data.Euclidean (Euclidean(..), GcdDomain(..))
 #endif
 import Data.Function
 import Data.Int
@@ -213,8 +213,9 @@
 
 evalTests :: TestTree
 evalTests = testGroup "eval" $ concat
-  [ evalTestGroup (Proxy :: Proxy (Poly U.Vector Int8))
-  , evalTestGroup (Proxy :: Proxy (Poly V.Vector Integer))
+  [ evalTestGroup  (Proxy :: Proxy (Poly U.Vector Int8))
+  , evalTestGroup  (Proxy :: Proxy (Poly V.Vector Integer))
+  , substTestGroup (Proxy :: Proxy (Poly U.Vector Int8))
   ]
 
 evalTestGroup
@@ -248,14 +249,37 @@
     e' :: Poly v a -> a -> a
     e' = S.eval
 
+substTestGroup
+  :: forall v a.
+     (Eq a, Num a, Semiring a, Arbitrary a, Show a, Eq (v (Word, a)), Show (v (Word, a)), G.Vector v (Word, a))
+  => Proxy (Poly v a)
+  -> [TestTree]
+substTestGroup _ =
+  [ testProperty "subst x p = p" $
+    \p -> e X p === p
+  , testProperty "subst (monomial 0 c) p = monomial 0 c" $
+    \c p -> e (monomial 0 c) p === monomial 0 c
+  , testProperty "subst' x p = p" $
+    \p -> e' S.X p === p
+  , testProperty "subst' (S.monomial 0 c) p = S.monomial 0 c" $
+    \c p -> e' (S.monomial 0 c) p === S.monomial 0 c
+  ]
+  where
+    e :: Poly v a -> Poly v a -> Poly v a
+    e = subst
+    e' :: Poly v a -> Poly v a -> Poly v a
+    e' = S.subst
+
 derivTests :: TestTree
 derivTests = testGroup "deriv"
   [ testProperty "deriv = S.deriv" $
     \(p :: Poly V.Vector Integer) -> deriv p === S.deriv p
+#if MIN_VERSION_semirings(0,5,0)
   , testProperty "integral = S.integral" $
     \(p :: Poly V.Vector Rational) -> integral p === S.integral p
   , testProperty "deriv . integral = id" $
     \(p :: Poly V.Vector Rational) -> deriv (integral p) === p
+#endif
   , testProperty "deriv c = 0" $
     \c -> deriv (monomial 0 c :: Poly V.Vector Int) === 0
   , testProperty "deriv cX = c" $
@@ -264,11 +288,9 @@
     \p q -> deriv (p + q) === (deriv p + deriv q :: Poly V.Vector Int)
   , testProperty "deriv (p * q) = p * deriv q + q * deriv p" $
     \p q -> deriv (p * q) === (p * deriv q + q * deriv p :: Poly V.Vector Int)
-  -- The following property takes too long for a regular test-suite
-  -- , testProperty "deriv (eval p q) = deriv q * eval (deriv p) q" $
+  -- , testProperty "deriv (subst p q) = deriv q * subst (deriv p) q" $
   --   \(p :: Poly V.Vector Int) (q :: Poly U.Vector Int) ->
-  --     deriv (eval (toPoly $ fmap (fmap $ monomial 0) $ unPoly p) q) ===
-  --       deriv q * eval (toPoly $ fmap (fmap $ monomial 0) $ unPoly $ deriv p) q
+  --     deriv (subst p q) === deriv q * subst (deriv p) q
   ]
 
 #if MIN_VERSION_semirings(0,4,2)