diff --git a/README.md b/README.md
--- a/README.md
+++ b/README.md
@@ -76,10 +76,10 @@
 cover main functionality of `Integral`, providing division with remainder and `gcd` / `lcm`:
 
 ```haskell
-> Data.Euclidean.gcd (X^2 + 7 * X + 6) (X^2 - 5 * X - 6) :: Data.Poly.UPoly Int
+> Data.Euclidean.gcd (X^2 + 7 * X + 6) (X^2 - 5 * X - 6) :: UPoly Int
 1 * X + 1
 
-> Data.Euclidean.quotRem (X^3 + 2) (X^2 - 1 :: Data.Poly.UPoly Double)
+> Data.Euclidean.quotRem (X^3 + 2) (X^2 - 1 :: UPoly Double)
 (1.0 * X + 0.0,1.0 * X + 2.0)
 ```
 
diff --git a/bench/DenseBench.hs b/bench/DenseBench.hs
--- a/bench/DenseBench.hs
+++ b/bench/DenseBench.hs
@@ -1,16 +1,24 @@
-{-# LANGUAGE CPP        #-}
-{-# LANGUAGE RankNTypes #-}
+{-# LANGUAGE CPP                        #-}
+{-# LANGUAGE GeneralizedNewtypeDeriving #-}
+{-# LANGUAGE RankNTypes                 #-}
+{-# LANGUAGE TypeApplications           #-}
 
 module DenseBench
   ( benchSuite
   ) where
 
 import Prelude hiding (quotRem, gcd)
+import Control.DeepSeq
 import Gauge.Main
 import Data.Poly
 import qualified Data.Vector.Unboxed as U
 #if MIN_VERSION_semirings(0,4,2)
+import Control.Exception
+import Data.Bits
+import Data.Coerce
 import Data.Euclidean
+import Data.Semiring (Semiring(..), Ring, isZero)
+import qualified Data.Semiring as S (negate)
 import qualified Data.Vector as V
 #endif
 
@@ -22,9 +30,12 @@
   , map benchDeriv    [100, 1000, 10000]
   , map benchIntegral [100, 1000, 10000]
 #if MIN_VERSION_semirings(0,4,2)
-  , map benchQuotRem  [10, 100]
-  , map benchGcdFrac  [10, 100]
-  , map benchGcd      [10, 100]
+  , map benchQuotRem    [10, 100]
+  , map benchGcd        [10, 100]
+  , map benchGcdExtRat  [10, 20, 40]
+  , map benchGcdFracRat [10, 20, 40]
+  , map benchGcdExtM    [10, 100, 1000]
+  , map benchGcdFracM   [10, 100, 1000]
 #endif
   ]
 
@@ -51,9 +62,18 @@
 benchGcd :: Int -> Benchmark
 benchGcd k = bench ("gcd/" ++ show k) $ nf doGcd k
 
-benchGcdFrac :: Int -> Benchmark
-benchGcdFrac k = bench ("gcdFrac/" ++ show k) $ nf doGcdFrac k
+benchGcdExtRat :: Int -> Benchmark
+benchGcdExtRat k = bench ("gcdExt/Rational/" ++ show k) $ nf (doGcdExt @Rational) k
 
+benchGcdFracRat :: Int -> Benchmark
+benchGcdFracRat k = bench ("gcdFrac/Rational/" ++ show k) $ nf (doGcdFrac @Rational) k
+
+benchGcdExtM :: Int -> Benchmark
+benchGcdExtM k = bench ("gcdExt/Mod2/" ++ show k) $ nf (doGcdExt @Mod2) k
+
+benchGcdFracM :: Int -> Benchmark
+benchGcdFracM k = bench ("gcdFrac/Mod2/" ++ show k) $ nf (doGcdFrac @Mod2) k
+
 #endif
 
 doBinOp :: (forall a. Num a => a -> a -> a) -> Int -> Int
@@ -83,25 +103,68 @@
 
 #if MIN_VERSION_semirings(0,4,2)
 
+gen1 :: Num a => Int -> a
+gen1 k = fromIntegral (truncate (pi * fromIntegral k :: Double) `mod` (k + 1))
+
+gen2 :: Num a => Int -> a
+gen2 k = fromIntegral (truncate (exp 1.0 * fromIntegral k :: Double) `mod` (k + 1))
+
 doQuotRem :: Int -> Double
 doQuotRem n = U.sum (unPoly qs) + U.sum (unPoly rs)
   where
-    xs = toPoly $ U.generate (2 * n) ((+ 1.0) . (* 2.0) . fromIntegral)
-    ys = toPoly $ U.generate n       ((+ 2.0) . (* 3.0) . fromIntegral)
+    xs = toPoly $ U.generate (2 * n) gen1
+    ys = toPoly $ U.generate n       gen2
     (qs, rs) = xs `quotRem` ys
 
 doGcd :: Int -> Integer
 doGcd n = V.sum gs
   where
-    xs = toPoly $ V.generate n ((+ 1) . (* 2) . fromIntegral)
-    ys = toPoly $ V.generate n ((+ 2) . (* 3) . fromIntegral)
+    xs = toPoly $ V.generate n gen1
+    ys = toPoly $ V.generate n gen2
     gs = unPoly $ xs `gcd` ys
 
-doGcdFrac :: Int -> Rational
+doGcdExt :: (Eq a, Num a, Field a) => Int -> a
+doGcdExt n = V.sum gs
+  where
+    xs = toPoly $ V.generate n gen1
+    ys = toPoly $ V.generate n gen2
+    gs = unPoly $ fst $ xs `gcdExt` ys
+
+doGcdFrac :: (Eq a, Num a, Field a) => Int -> a
 doGcdFrac n = V.sum gs
   where
-    xs = PolyOverFractional $ toPoly $ V.generate n ((+ 1) . (* 2) . fromIntegral)
-    ys = PolyOverFractional $ toPoly $ V.generate n ((+ 2) . (* 3) . fromIntegral)
-    gs = unPoly $ unPolyOverFractional $ xs `gcd` ys
+    xs = PolyOverField $ toPoly $ V.generate n gen1
+    ys = PolyOverField $ toPoly $ V.generate n gen2
+    gs = unPoly $ unPolyOverField $ xs `gcd` ys
+
+-- | Inspired by 'semirings'.
+newtype Mod2 = Mod2 { _getMod2 :: Bool }
+  deriving (Eq, NFData)
+
+instance Num Mod2 where
+  (+) = coerce (xor @Bool)
+  (*) = coerce (&&)
+  negate = id
+  abs    = id
+  signum = id
+  fromInteger = Mod2 . odd
+
+instance Semiring Mod2 where
+  plus  = coerce (xor @Bool)
+  times = coerce (&&)
+  fromNatural = Mod2 . odd
+
+instance Ring Mod2 where
+  negate = id
+
+instance GcdDomain Mod2 where
+
+instance Euclidean Mod2 where
+  degree = const 0
+  quotRem x y
+    | isZero y  = throw DivideByZero
+    | otherwise = (x, zero)
+
+instance Field Mod2
 
 #endif
diff --git a/bench/SparseBench.hs b/bench/SparseBench.hs
--- a/bench/SparseBench.hs
+++ b/bench/SparseBench.hs
@@ -67,4 +67,3 @@
 doIntegral xs = U.foldl' (\acc (_, x) -> acc + x) 0 zs
   where
     zs = unPoly $ integral xs
-
diff --git a/changelog.md b/changelog.md
--- a/changelog.md
+++ b/changelog.md
@@ -1,3 +1,10 @@
+# 0.3.2.0
+
+* Add `NFData` instance.
+* Implement extended GCD.
+* Rename `PolyOverFractional` to `PolyOverField`.
+* Add `integral` with `Semiring`-based interface.
+
 # 0.3.1.0
 
 * Implement Karatsuba multiplication.
diff --git a/poly.cabal b/poly.cabal
--- a/poly.cabal
+++ b/poly.cabal
@@ -1,5 +1,5 @@
 name: poly
-version: 0.3.1.0
+version: 0.3.2.0
 synopsis: Polynomials
 description:
   Polynomials backed by `Vector`.
@@ -30,14 +30,15 @@
     Data.Poly.Sparse.Semiring
   other-modules:
     Data.Poly.Internal.Dense
-    Data.Poly.Internal.Dense.Fractional
+    Data.Poly.Internal.Dense.Field
     Data.Poly.Internal.Dense.GcdDomain
-    Data.Poly.Internal.PolyOverFractional
+    Data.Poly.Internal.PolyOverField
     Data.Poly.Internal.Sparse
-    Data.Poly.Internal.Sparse.Fractional
+    Data.Poly.Internal.Sparse.Field
     Data.Poly.Internal.Sparse.GcdDomain
   build-depends:
     base >= 4.9 && < 5,
+    deepseq >= 1.1 && < 1.5,
     primitive >= 0.6,
     semirings >= 0.2,
     vector >= 0.12.0.2,
@@ -68,6 +69,7 @@
 benchmark poly-gauge
   build-depends:
     base >=4.9 && <5,
+    deepseq >= 1.1 && < 1.5,
     gauge >= 0.1,
     poly,
     semirings >= 0.2,
diff --git a/src/Data/Poly.hs b/src/Data/Poly.hs
--- a/src/Data/Poly.hs
+++ b/src/Data/Poly.hs
@@ -7,10 +7,8 @@
 -- Dense polynomials and a 'Num'-based interface.
 --
 
-{-# LANGUAGE CPP                        #-}
-{-# LANGUAGE FlexibleInstances          #-}
-{-# LANGUAGE GeneralizedNewtypeDeriving #-}
-{-# LANGUAGE PatternSynonyms            #-}
+{-# LANGUAGE CPP             #-}
+{-# LANGUAGE PatternSynonyms #-}
 
 module Data.Poly
   ( Poly
@@ -27,14 +25,18 @@
   , deriv
   , integral
 #if MIN_VERSION_semirings(0,4,2)
-  -- * Fractional coefficients
-  , PolyOverFractional(..)
+  -- * Polynomials over 'Field'
+  , PolyOverField(..)
+  , gcdExt
+  , PolyOverFractional
+  , pattern PolyOverFractional
+  , unPolyOverFractional
 #endif
   ) where
 
 import Data.Poly.Internal.Dense
 #if MIN_VERSION_semirings(0,4,2)
-import Data.Poly.Internal.Dense.Fractional ()
+import Data.Poly.Internal.Dense.Field (gcdExt)
 import Data.Poly.Internal.Dense.GcdDomain ()
-import Data.Poly.Internal.PolyOverFractional
+import Data.Poly.Internal.PolyOverField
 #endif
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
@@ -36,14 +36,19 @@
   , pattern X'
   , eval'
   , deriv'
+#if MIN_VERSION_semirings(0,5,0)
+  , integral'
+#endif
   ) where
 
-import Prelude hiding (quotRem, rem, gcd, lcm, (^))
+import Prelude hiding (quotRem, quot, rem, gcd, lcm, (^))
+import Control.DeepSeq (NFData)
 import Control.Monad
 import Control.Monad.Primitive
 import Control.Monad.ST
+import Data.Bits
 import Data.List (foldl', intersperse)
-import Data.Semiring (Semiring(..))
+import Data.Semiring (Semiring(..), Ring())
 import qualified Data.Semiring as Semiring
 import qualified Data.Vector as V
 import qualified Data.Vector.Generic as G
@@ -54,6 +59,9 @@
 import Data.Semigroup
 import Numeric.Natural
 #endif
+#if MIN_VERSION_semirings(0,5,0)
+import Data.Euclidean (Field, quot)
+#endif
 
 -- | Polynomials of one variable with coefficients from @a@,
 -- backed by a 'G.Vector' @v@ (boxed, unboxed, storable, etc.).
@@ -76,7 +84,7 @@
   -- ^ Convert 'Poly' to a vector of coefficients
   -- (first element corresponds to a constant term).
   }
-  deriving (Eq, Ord)
+  deriving (Eq, NFData, Ord)
 
 instance (Eq a, Semiring a, G.Vector v a) => IsList (Poly v a) where
   type Item (Poly v a) = a
@@ -139,7 +147,7 @@
   abs = id
   signum = const 1
   fromInteger n = case fromInteger n of
-    0 -> Poly $ G.empty
+    0 -> Poly G.empty
     m -> Poly $ G.singleton m
   Poly xs * Poly ys = toPoly $ karatsuba xs ys
   {-# INLINE (+) #-}
@@ -168,7 +176,7 @@
   {-# INLINE fromNatural #-}
 #endif
 
-instance (Eq a, Semiring.Ring a, G.Vector v a) => Semiring.Ring (Poly v a) where
+instance (Eq a, Ring a, G.Vector v a) => Ring (Poly v a) where
   negate (Poly xs) = Poly $ G.map Semiring.negate xs
 
 dropWhileEnd
@@ -176,7 +184,7 @@
   => (a -> Bool)
   -> v a
   -> v a
-dropWhileEnd p xs = G.basicUnsafeSlice 0 (go (G.basicLength xs)) xs
+dropWhileEnd p xs = G.unsafeSlice 0 (go (G.length xs)) xs
   where
     go 0 = 0
     go n = if p (G.unsafeIndex xs (n - 1)) then go (n - 1) else n
@@ -187,12 +195,12 @@
   => (a -> Bool)
   -> G.Mutable v (PrimState m) a
   -> m (G.Mutable v (PrimState m) a)
-dropWhileEndM p xs = go (MG.basicLength xs)
+dropWhileEndM p xs = go (MG.length xs)
   where
-    go 0 = pure $ MG.basicUnsafeSlice 0 0 xs
+    go 0 = pure $ MG.unsafeSlice 0 0 xs
     go n = do
       x <- MG.unsafeRead xs (n - 1)
-      if p x then go (n - 1) else pure (MG.basicUnsafeSlice 0 n xs)
+      if p x then go (n - 1) else pure (MG.unsafeSlice 0 n xs)
 {-# INLINE dropWhileEndM #-}
 
 plusPoly
@@ -202,17 +210,17 @@
   -> v a
   -> v a
 plusPoly add xs ys = runST $ do
-  let lenXs = G.basicLength xs
-      lenYs = G.basicLength ys
+  let lenXs = G.length xs
+      lenYs = G.length ys
       lenMn = lenXs `min` lenYs
       lenMx = lenXs `max` lenYs
 
-  zs <- MG.basicUnsafeNew lenMx
+  zs <- MG.unsafeNew lenMx
   forM_ [0 .. lenMn - 1] $ \i ->
     MG.unsafeWrite zs i (add (G.unsafeIndex xs i) (G.unsafeIndex ys i))
   G.unsafeCopy
-    (MG.basicUnsafeSlice lenMn (lenMx - lenMn) zs)
-    (G.basicUnsafeSlice  lenMn (lenMx - lenMn) (if lenXs <= lenYs then ys else xs))
+    (MG.unsafeSlice lenMn (lenMx - lenMn) zs)
+    (G.unsafeSlice  lenMn (lenMx - lenMn) (if lenXs <= lenYs then ys else xs))
 
   G.unsafeFreeze zs
 {-# INLINE plusPoly #-}
@@ -225,12 +233,12 @@
   -> v a
   -> v a
 minusPoly neg sub xs ys = runST $ do
-  let lenXs = G.basicLength xs
-      lenYs = G.basicLength ys
+  let lenXs = G.length xs
+      lenYs = G.length ys
       lenMn = lenXs `min` lenYs
       lenMx = lenXs `max` lenYs
 
-  zs <- MG.basicUnsafeNew lenMx
+  zs <- MG.unsafeNew lenMx
   forM_ [0 .. lenMn - 1] $ \i ->
     MG.unsafeWrite zs i (sub (G.unsafeIndex xs i) (G.unsafeIndex ys i))
 
@@ -238,8 +246,8 @@
     then forM_ [lenXs .. lenYs - 1] $ \i ->
       MG.unsafeWrite zs i (neg (G.unsafeIndex ys i))
     else G.unsafeCopy
-      (MG.basicUnsafeSlice lenYs (lenXs - lenYs) zs)
-      (G.basicUnsafeSlice  lenYs (lenXs - lenYs) xs)
+      (MG.unsafeSlice lenYs (lenXs - lenYs) zs)
+      (G.unsafeSlice  lenYs (lenXs - lenYs) xs)
 
   G.unsafeFreeze zs
 {-# INLINE minusPoly #-}
@@ -256,31 +264,31 @@
   | lenXs <= karatsubaThreshold || lenYs <= karatsubaThreshold
   = convolution 0 (+) (*) xs ys
   | otherwise = runST $ do
-    zs <- MG.basicUnsafeNew lenZs
+    zs <- MG.unsafeNew lenZs
     forM_ [0 .. lenZs - 1] $ \k -> do
-      let z0 = if k < G.basicLength zs0
+      let z0 = if k < G.length zs0
                then G.unsafeIndex zs0 k
                else 0
-          z11 = if k - m >= 0 && k - m < G.basicLength zs11
+          z11 = if k - m >= 0 && k - m < G.length zs11
                then G.unsafeIndex zs11 (k - m)
                else 0
-          z10 = if k - m >= 0 && k - m < G.basicLength zs0
+          z10 = if k - m >= 0 && k - m < G.length zs0
                then G.unsafeIndex zs0 (k - m)
                else 0
-          z12 = if k - m >= 0 && k - m < G.basicLength zs2
+          z12 = if k - m >= 0 && k - m < G.length zs2
                then G.unsafeIndex zs2 (k - m)
                else 0
-          z2 = if k - 2 * m >= 0 && k - 2 * m < G.basicLength zs2
+          z2 = if k - 2 * m >= 0 && k - 2 * m < G.length zs2
                then G.unsafeIndex zs2 (k - 2 * m)
                else 0
       MG.unsafeWrite zs k (z0 + (z11 - z10 - z12) + z2)
     G.unsafeFreeze zs
   where
-    lenXs = G.basicLength xs
-    lenYs = G.basicLength ys
+    lenXs = G.length xs
+    lenYs = G.length ys
     lenZs = lenXs + lenYs - 1
 
-    m    = ((lenXs `min` lenYs) + 1) `quot` 2
+    m    = ((lenXs `min` lenYs) + 1) `shiftR` 1
 
     xs0  = G.slice 0 m xs
     xs1  = G.slice m (lenXs - m) xs
@@ -303,18 +311,15 @@
   -> v a
   -> v a
 convolution zer add mul xs ys
-  | G.null xs || G.null ys = G.empty
-  | otherwise = runST $ do
-    let lenXs = G.basicLength xs
-        lenYs = G.basicLength ys
-        lenZs = lenXs + lenYs - 1
-    zs <- MG.basicUnsafeNew lenZs
-    forM_ [0 .. lenZs - 1] $ \k -> do
-      let is = [max (k - lenYs + 1) 0 .. min k (lenXs - 1)]
-          acc = foldl' add zer $ flip map is $ \i ->
-            mul (G.unsafeIndex xs i) (G.unsafeIndex ys (k - i))
-      MG.unsafeWrite zs k acc
-    G.unsafeFreeze zs
+  | lenXs == 0 || lenYs == 0 = G.empty
+  | otherwise = G.generate lenZs $ \k -> foldl'
+    (\acc i -> acc `add` mul (G.unsafeIndex xs i) (G.unsafeIndex ys (k - i)))
+    zer
+    [max (k - lenYs + 1) 0 .. min k (lenXs - 1)]
+  where
+    lenXs = G.length xs
+    lenYs = G.length ys
+    lenZs = lenXs + lenYs - 1
 {-# INLINE convolution #-}
 
 -- | Create a monomial from a power and a coefficient.
@@ -338,8 +343,8 @@
   -> v a
   -> v a
 scaleInternal zer mul yp yc xs = runST $ do
-  let lenXs = G.basicLength xs
-  zs <- MG.basicUnsafeNew (fromIntegral yp + lenXs)
+  let lenXs = G.length xs
+  zs <- MG.unsafeNew (fromIntegral yp + lenXs)
   forM_ [0 .. fromIntegral yp - 1] $ \k ->
     MG.unsafeWrite zs k zer
   forM_ [0 .. lenXs - 1] $ \k ->
@@ -372,12 +377,12 @@
 -- 1 * X^2 + 2 * X + 2
 eval :: (Num a, G.Vector v a) => Poly v a -> a -> a
 eval (Poly cs) x = fst' $
-  G.foldl' (\(acc :*: xn) cn -> (acc + cn * xn :*: x * xn)) (0 :*: 1) cs
+  G.foldl' (\(acc :*: xn) cn -> acc + cn * xn :*: x * xn) (0 :*: 1) cs
 {-# INLINE eval #-}
 
 eval' :: (Semiring a, G.Vector v a) => Poly v a -> a -> a
 eval' (Poly cs) x = fst' $
-  G.foldl' (\(acc :*: xn) cn -> (acc `plus` cn `times` xn :*: x `times` xn)) (zero :*: one) cs
+  G.foldl' (\(acc :*: xn) cn -> acc `plus` cn `times` xn :*: x `times` xn) (zero :*: one) cs
 {-# INLINE eval' #-}
 
 -- | Take a derivative.
@@ -416,14 +421,29 @@
 integral (Poly xs)
   | G.null xs = Poly G.empty
   | otherwise = toPoly $ runST $ do
-    zs <- MG.basicUnsafeNew (lenXs + 1)
+    zs <- MG.unsafeNew (lenXs + 1)
     MG.unsafeWrite zs 0 0
     forM_ [0 .. lenXs - 1] $ \i ->
       MG.unsafeWrite zs (i + 1) (G.unsafeIndex xs i * recip (fromIntegral i + 1))
     G.unsafeFreeze zs
     where
-      lenXs = G.basicLength xs
+      lenXs = G.length xs
 {-# INLINE integral #-}
+
+#if MIN_VERSION_semirings(0,5,0)
+integral' :: (Eq a, Field a, G.Vector v a) => Poly v a -> Poly v a
+integral' (Poly xs)
+  | G.null xs = Poly G.empty
+  | otherwise = toPoly' $ runST $ do
+    zs <- MG.unsafeNew (lenXs + 1)
+    MG.unsafeWrite zs zero zero
+    forM_ [0 .. lenXs - 1] $ \i ->
+      MG.unsafeWrite zs (i + 1) (G.unsafeIndex xs i `quot` Semiring.fromIntegral (i + 1))
+    G.unsafeFreeze zs
+    where
+      lenXs = G.length xs
+{-# INLINE integral' #-}
+#endif
 
 -- | Create an identity polynomial.
 pattern X :: (Eq a, Num a, G.Vector v a, Eq (v a)) => Poly v a
diff --git a/src/Data/Poly/Internal/Dense/Field.hs b/src/Data/Poly/Internal/Dense/Field.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Poly/Internal/Dense/Field.hs
@@ -0,0 +1,189 @@
+-- |
+-- Module:      Data.Poly.Internal.Dense.Field
+-- Copyright:   (c) 2019 Andrew Lelechenko
+-- Licence:     BSD3
+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>
+--
+-- GcdDomain for Field underlying
+--
+
+{-# LANGUAGE ConstraintKinds            #-}
+{-# LANGUAGE CPP                        #-}
+{-# LANGUAGE FlexibleInstances          #-}
+{-# LANGUAGE PatternSynonyms            #-}
+{-# LANGUAGE ScopedTypeVariables        #-}
+{-# LANGUAGE TypeFamilies               #-}
+
+{-# OPTIONS_GHC -fno-warn-orphans #-}
+
+#if MIN_VERSION_semirings(0,4,2)
+
+module Data.Poly.Internal.Dense.Field
+  ( fieldGcd
+  , gcdExt
+  ) where
+
+import Prelude hiding (quotRem, quot, rem, gcd)
+import Control.Exception
+import Control.Monad
+import Control.Monad.Primitive
+import Control.Monad.ST
+import Data.Euclidean
+#if !MIN_VERSION_semirings(0,5,0)
+import Data.Semiring (Ring)
+#endif
+import Data.Semiring (times, minus, zero, one)
+import qualified Data.Vector.Generic as G
+import qualified Data.Vector.Generic.Mutable as MG
+
+import Data.Poly.Internal.Dense
+import Data.Poly.Internal.Dense.GcdDomain ()
+
+#if !MIN_VERSION_semirings(0,5,0)
+type Field a = (Euclidean a, Ring a, Fractional a)
+#endif
+
+instance (Eq a, Eq (v a), Field a, G.Vector v a) => Euclidean (Poly v a) where
+  degree (Poly xs) = fromIntegral (G.length xs)
+
+  quotRem (Poly xs) (Poly ys) = (toPoly' qs, toPoly' rs)
+    where
+      (qs, rs) = quotientAndRemainder xs ys
+  {-# INLINE quotRem #-}
+
+  rem (Poly xs) (Poly ys) = toPoly' $ remainder xs ys
+  {-# INLINE rem #-}
+
+quotientAndRemainder
+  :: (Field a, G.Vector v a)
+  => v a
+  -> v a
+  -> (v a, v a)
+quotientAndRemainder xs ys
+  | G.null ys = throw DivideByZero
+  | G.length xs < G.length ys = (G.empty, xs)
+  | otherwise = runST $ do
+    let lenXs = G.length xs
+        lenYs = G.length ys
+        lenQs = lenXs - lenYs + 1
+    qs <- MG.unsafeNew lenQs
+    rs <- MG.unsafeNew lenXs
+    G.unsafeCopy rs xs
+    forM_ [lenQs - 1, lenQs - 2 .. 0] $ \i -> do
+      r <- MG.unsafeRead rs (lenYs - 1 + i)
+      let q = r `quot` G.unsafeLast ys
+      MG.unsafeWrite qs i q
+      forM_ [0 .. lenYs - 1] $ \k -> do
+        MG.unsafeModify rs (\c -> c `minus` q `times` G.unsafeIndex ys k) (i + k)
+    let rs' = MG.unsafeSlice 0 lenYs rs
+    (,) <$> G.unsafeFreeze qs <*> G.unsafeFreeze rs'
+{-# INLINE quotientAndRemainder #-}
+
+remainder
+  :: (Field a, G.Vector v a)
+  => v a
+  -> v a
+  -> v a
+remainder xs ys
+  | G.null ys = throw DivideByZero
+  | otherwise = runST $ do
+    rs <- G.thaw xs
+    ys' <- G.unsafeThaw ys
+    remainderM rs ys'
+    G.unsafeFreeze $ MG.unsafeSlice 0 (G.length xs `min` G.length ys) rs
+{-# INLINE remainder #-}
+
+remainderM
+  :: (PrimMonad m, Field a, G.Vector v a)
+  => G.Mutable v (PrimState m) a
+  -> G.Mutable v (PrimState m) a
+  -> m ()
+remainderM xs ys
+  | MG.null ys = throw DivideByZero
+  | MG.length xs < MG.length ys = pure ()
+  | otherwise = do
+    let lenXs = MG.length xs
+        lenYs = MG.length ys
+        lenQs = lenXs - lenYs + 1
+    yLast <- MG.unsafeRead ys (lenYs - 1)
+    forM_ [lenQs - 1, lenQs - 2 .. 0] $ \i -> do
+      r <- MG.unsafeRead xs (lenYs - 1 + i)
+      forM_ [0 .. lenYs - 1] $ \k -> do
+        y <- MG.unsafeRead ys k
+        -- do not move r / yLast outside the loop,
+        -- because of numerical instability
+        MG.unsafeModify xs (\c -> c `minus` r `times` y `quot` yLast) (i + k)
+{-# INLINE remainderM #-}
+
+fieldGcd
+  :: (Eq a, Field a, G.Vector v a)
+  => Poly v a
+  -> Poly v a
+  -> Poly v a
+fieldGcd (Poly xs) (Poly ys) = toPoly' $ runST $ do
+  xs' <- G.thaw xs
+  ys' <- G.thaw ys
+  gcdM xs' ys'
+{-# INLINE fieldGcd #-}
+
+gcdM
+  :: (PrimMonad m, Eq a, Field a, G.Vector v a)
+  => G.Mutable v (PrimState m) a
+  -> G.Mutable v (PrimState m) a
+  -> m (v a)
+gcdM xs ys = do
+  ys' <- dropWhileEndM (== zero) ys
+  if MG.null ys' then G.unsafeFreeze xs else do
+    remainderM xs ys'
+    gcdM ys' xs
+{-# INLINE gcdM #-}
+
+-- | Execute the extended Euclidean algorithm.
+-- For polynomials @a@ and @b@, compute their unique greatest common divisor @g@
+-- and the unique coefficient polynomial @s@ satisfying @as + bt = g@,
+-- such that either @g@ is monic, or @g = 0@ and @s@ is monic, or @g = s = 0@.
+--
+-- >>> gcdExt (X^2 + 1 :: UPoly Double) (X^3 + 3 * X :: UPoly Double)
+-- (1.0, 0.5 * X^2 + (-0.0) * X + 1.0)
+-- >>> gcdExt (X^3 + 3 * X :: UPoly Double) (3 * X^4 + 3 * X^2 :: UPoly Double)
+-- (1.0 * X + 0.0,(-0.16666666666666666) * X^2 + (-0.0) * X + 0.3333333333333333)
+gcdExt
+  :: (Eq a, Field a, G.Vector v a, Eq (v a))
+  => Poly v a
+  -> Poly v a
+  -> (Poly v a, Poly v a)
+gcdExt xs ys = case scaleMonic gs of
+  Just (c', gs') -> (gs', scale' zero c' ss)
+  Nothing -> case scaleMonic ss of
+    Just (_, ss') -> (zero, ss')
+    Nothing -> (zero, zero)
+  where
+    (gs, ss) = go ys xs zero one
+      where
+        go r' r s' s
+          | r' == zero = (r, s)
+          | otherwise  = case r `quotRem` r' of
+            (q, r'') -> go r'' r' (s `minus` q `times` s') s'
+{-# INLINE gcdExt #-}
+
+-- | Scale a non-zero polynomial such that its leading coefficient is one,
+-- returning the reciprocal of the leading coefficient in the scaling.
+--
+-- >>> scaleMonic (X^3 + 3 * X :: UPoly Double)
+-- Just (1.0, 1.0 * X^3 + 0.0 * X^2 + 3.0 * X + 0.0)
+-- >>> scaleMonic (3 * X^4 + 3 * X^2 :: UPoly Double)
+-- Just (0.3333333333333333, 1.0 * X^4 + 0.0 * X^3 + 1.0 * X^2 + 0.0 * X + 0.0)
+scaleMonic
+  :: (Eq a, Field a, G.Vector v a, Eq (v a))
+  => Poly v a
+  -> Maybe (a, Poly v a)
+scaleMonic xs = case leading xs of
+  Nothing -> Nothing
+  Just (_, c) -> let c' = one `quot` c in Just (c', scale' zero c' xs)
+{-# INLINE scaleMonic #-}
+
+#else
+
+module Data.Poly.Internal.Dense.Field () where
+
+#endif
diff --git a/src/Data/Poly/Internal/Dense/Fractional.hs b/src/Data/Poly/Internal/Dense/Fractional.hs
deleted file mode 100644
--- a/src/Data/Poly/Internal/Dense/Fractional.hs
+++ /dev/null
@@ -1,138 +0,0 @@
--- |
--- Module:      Data.Poly.Internal.Dense.Fractional
--- Copyright:   (c) 2019 Andrew Lelechenko
--- Licence:     BSD3
--- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>
---
--- GcdDomain for Fractional underlying
---
-
-{-# LANGUAGE CPP                        #-}
-{-# LANGUAGE FlexibleInstances          #-}
-{-# LANGUAGE GeneralizedNewtypeDeriving #-}
-{-# LANGUAGE PatternSynonyms            #-}
-{-# LANGUAGE ScopedTypeVariables        #-}
-{-# LANGUAGE TypeFamilies               #-}
-{-# LANGUAGE ViewPatterns               #-}
-
-{-# OPTIONS_GHC -fno-warn-orphans #-}
-
-#if MIN_VERSION_semirings(0,4,2)
-
-module Data.Poly.Internal.Dense.Fractional
-  ( fractionalGcd
-  ) where
-
-import Prelude hiding (rem, gcd)
-import Control.Exception
-import Control.Monad
-import Control.Monad.Primitive
-import Control.Monad.ST
-import Data.Euclidean
-import qualified Data.Semiring as Semiring
-import qualified Data.Vector.Generic as G
-import qualified Data.Vector.Generic.Mutable as MG
-
-import Data.Poly.Internal.Dense
-import Data.Poly.Internal.Dense.GcdDomain ()
-
-instance (Eq a, Eq (v a), Semiring.Ring a, GcdDomain a, Fractional a, G.Vector v a) => Euclidean (Poly v a) where
-  degree (Poly xs) = fromIntegral (G.basicLength xs)
-
-  quotRem (Poly xs) (Poly ys) = (toPoly qs, toPoly rs)
-    where
-      (qs, rs) = quotientAndRemainder xs ys
-  {-# INLINE quotRem #-}
-
-  rem (Poly xs) (Poly ys) = toPoly $ remainder xs ys
-  {-# INLINE rem #-}
-
-quotientAndRemainder
-  :: (Fractional a, G.Vector v a)
-  => v a
-  -> v a
-  -> (v a, v a)
-quotientAndRemainder xs ys
-  | G.null ys = throw DivideByZero
-  | G.basicLength xs < G.basicLength ys = (G.empty, xs)
-  | otherwise = runST $ do
-    let lenXs = G.basicLength xs
-        lenYs = G.basicLength ys
-        lenQs = lenXs - lenYs + 1
-    qs <- MG.basicUnsafeNew lenQs
-    rs <- MG.basicUnsafeNew lenXs
-    G.unsafeCopy rs xs
-    forM_ [lenQs - 1, lenQs - 2 .. 0] $ \i -> do
-      r <- MG.unsafeRead rs (lenYs - 1 + i)
-      let q = r / G.unsafeLast ys
-      MG.unsafeWrite qs i q
-      forM_ [0 .. lenYs - 1] $ \k -> do
-        MG.unsafeModify rs (\c -> c - q * G.unsafeIndex ys k) (i + k)
-    let rs' = MG.basicUnsafeSlice 0 lenYs rs
-    (,) <$> G.unsafeFreeze qs <*> G.unsafeFreeze rs'
-{-# INLINE quotientAndRemainder #-}
-
-remainder
-  :: (Fractional a, G.Vector v a)
-  => v a
-  -> v a
-  -> v a
-remainder xs ys
-  | G.null ys = throw DivideByZero
-  | otherwise = runST $ do
-    rs <- G.thaw xs
-    ys' <- G.unsafeThaw ys
-    remainderM rs ys'
-    G.unsafeFreeze $ MG.basicUnsafeSlice 0 (G.basicLength xs `min` G.basicLength ys) rs
-{-# INLINE remainder #-}
-
-remainderM
-  :: (PrimMonad m, Fractional a, G.Vector v a)
-  => G.Mutable v (PrimState m) a
-  -> G.Mutable v (PrimState m) a
-  -> m ()
-remainderM xs ys
-  | MG.null ys = throw DivideByZero
-  | MG.basicLength xs < MG.basicLength ys = pure ()
-  | otherwise = do
-    let lenXs = MG.basicLength xs
-        lenYs = MG.basicLength ys
-        lenQs = lenXs - lenYs + 1
-    yLast <- MG.unsafeRead ys (lenYs - 1)
-    forM_ [lenQs - 1, lenQs - 2 .. 0] $ \i -> do
-      r <- MG.unsafeRead xs (lenYs - 1 + i)
-      forM_ [0 .. lenYs - 1] $ \k -> do
-        y <- MG.unsafeRead ys k
-        -- do not move r / yLast outside the loop,
-        -- because of numerical instability
-        MG.unsafeModify xs (\c -> c - r * y / yLast) (i + k)
-{-# INLINE remainderM #-}
-
-fractionalGcd
-  :: (Eq a, Fractional a, G.Vector v a)
-  => Poly v a
-  -> Poly v a
-  -> Poly v a
-fractionalGcd (Poly xs) (Poly ys) = toPoly $ runST $ do
-  xs' <- G.thaw xs
-  ys' <- G.thaw ys
-  gcdM xs' ys'
-{-# INLINE fractionalGcd #-}
-
-gcdM
-  :: (PrimMonad m, Eq a, Fractional a, G.Vector v a)
-  => G.Mutable v (PrimState m) a
-  -> G.Mutable v (PrimState m) a
-  -> m (v a)
-gcdM xs ys = do
-  ys' <- dropWhileEndM (== 0) ys
-  if MG.null ys' then G.unsafeFreeze xs else do
-    remainderM xs ys'
-    gcdM ys' xs
-{-# INLINE gcdM #-}
-
-#else
-
-module Data.Poly.Internal.Dense.Fractional () where
-
-#endif
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
@@ -9,11 +9,9 @@
 
 {-# LANGUAGE CPP                        #-}
 {-# LANGUAGE FlexibleInstances          #-}
-{-# LANGUAGE GeneralizedNewtypeDeriving #-}
 {-# LANGUAGE PatternSynonyms            #-}
 {-# LANGUAGE ScopedTypeVariables        #-}
 {-# LANGUAGE TypeFamilies               #-}
-{-# LANGUAGE ViewPatterns               #-}
 
 {-# OPTIONS_GHC -fno-warn-orphans #-}
 
@@ -28,18 +26,16 @@
 import Control.Monad.Primitive
 import Control.Monad.ST
 import Data.Euclidean
-import Data.Semiring (Semiring(..), isZero)
-import qualified Data.Semiring as Semiring
+import Data.Semiring (Semiring(..), Ring(), isZero, minus)
 import qualified Data.Vector.Generic as G
 import qualified Data.Vector.Generic.Mutable as MG
 
 import Data.Poly.Internal.Dense
 
--- | Consider using 'Data.Poly.Semiring.PolyOverFractional' wrapper,
+-- | Consider using 'Data.Poly.Semiring.PolyOverField' wrapper,
 -- which provides a much faster implementation of
--- 'Data.Euclidean.gcd' for 'Fractional'
--- coefficients.
-instance (Eq a, Semiring.Ring a, GcdDomain a, Eq (v a), G.Vector v a) => GcdDomain (Poly v a) where
+-- 'Data.Euclidean.gcd' for polynomials over 'Field'.
+instance (Eq a, Ring a, GcdDomain a, Eq (v a), G.Vector v a) => GcdDomain (Poly v a) where
   divide (Poly xs) (Poly ys) =
     toPoly' <$> quotient xs ys
 
@@ -50,7 +46,7 @@
   {-# INLINE gcd #-}
 
 gcdNonEmpty
-  :: (Eq a, Semiring.Ring a, GcdDomain a, G.Vector v a)
+  :: (Eq a, Ring a, GcdDomain a, G.Vector v a)
   => v a
   -> v a
   -> v a
@@ -62,7 +58,7 @@
     ys' <- G.thaw ys
     zs' <- gcdM xs' ys'
 
-    let lenZs = MG.basicLength zs'
+    let lenZs = MG.length zs'
         go acc 0 = pure acc
         go acc n = do
           t <- MG.unsafeRead zs' (n - 1)
@@ -80,7 +76,7 @@
     G.unsafeFreeze zs'
 
 gcdM
-  :: (PrimMonad m, Eq a, Semiring.Ring a, GcdDomain a, G.Vector v a)
+  :: (PrimMonad m, Eq a, Ring a, GcdDomain a, G.Vector v a)
   => G.Mutable v (PrimState m) a
   -> G.Mutable v (PrimState m) a
   -> m (G.Mutable v (PrimState m) a)
@@ -88,8 +84,8 @@
   | MG.null xs = pure ys
   | MG.null ys = pure xs
   | otherwise = do
-  let lenXs = MG.basicLength xs
-      lenYs = MG.basicLength ys
+  let lenXs = MG.length xs
+      lenYs = MG.length ys
   xLast <- MG.unsafeRead xs (lenXs - 1)
   yLast <- MG.unsafeRead ys (lenYs - 1)
   let z = xLast `lcm` yLast
@@ -105,7 +101,7 @@
       x <- MG.unsafeRead xs i
       MG.unsafeModify
         ys
-        (\y -> (y `times` zy) `plus` Semiring.negate (x `times` zx))
+        (\y -> (y `times` zy) `minus` x `times` zx)
         (i + lenYs - lenXs)
     forM_ [0 .. lenYs - lenXs - 1] $
       MG.unsafeModify ys (`times` zy)
@@ -116,7 +112,7 @@
       y <- MG.unsafeRead ys i
       MG.unsafeModify
         xs
-        (\x -> (x `times` zx) `plus` Semiring.negate (y `times` zy))
+        (\x -> (x `times` zx) `minus` y `times` zy)
         (i + lenXs - lenYs)
     forM_ [0 .. lenXs - lenYs - 1] $
       MG.unsafeModify xs (`times` zx)
@@ -128,7 +124,7 @@
   :: (Eq a, Semiring a, PrimMonad m, G.Vector v a)
   => G.Mutable v (PrimState m) a
   -> m Bool
-isZeroM xs = go (MG.basicLength xs)
+isZeroM xs = go (MG.length xs)
   where
     go 0 = pure True
     go n = do
@@ -137,20 +133,20 @@
 {-# INLINE isZeroM #-}
 
 quotient
-  :: (Eq a, Eq (v a), Semiring.Ring a, GcdDomain a, G.Vector v a)
+  :: (Eq a, Eq (v a), Ring a, GcdDomain a, G.Vector v a)
   => v a
   -> v a
   -> Maybe (v a)
 quotient xs ys
   | G.null ys = throw DivideByZero
   | G.null xs = Just xs
-  | G.basicLength xs < G.basicLength ys = Nothing
+  | G.length xs < G.length ys = Nothing
   | otherwise = runST $ do
-    let lenXs = G.basicLength xs
-        lenYs = G.basicLength ys
+    let lenXs = G.length xs
+        lenYs = G.length ys
         lenQs = lenXs - lenYs + 1
-    qs <- MG.basicUnsafeNew lenQs
-    rs <- MG.basicUnsafeNew lenXs
+    qs <- MG.unsafeNew lenQs
+    rs <- MG.unsafeNew lenXs
     G.unsafeCopy rs xs
 
     let go i
@@ -168,7 +164,7 @@
                 forM_ [0 .. lenYs - 1] $ \k -> do
                   MG.unsafeModify
                     rs
-                    (\c -> c `plus` (Semiring.negate $ q `times` G.unsafeIndex ys k))
+                    (\c -> c `minus` q `times` G.unsafeIndex ys k)
                     (i + k)
                 go (i - 1)
 
diff --git a/src/Data/Poly/Internal/PolyOverField.hs b/src/Data/Poly/Internal/PolyOverField.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Poly/Internal/PolyOverField.hs
@@ -0,0 +1,75 @@
+-- |
+-- Module:      Data.Poly.Internal.PolyOverField
+-- Copyright:   (c) 2019 Andrew Lelechenko
+-- Licence:     BSD3
+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>
+--
+-- Wrapper with a more efficient 'Euclidean' instance.
+--
+
+{-# LANGUAGE CPP                        #-}
+{-# LANGUAGE ConstraintKinds            #-}
+{-# LANGUAGE FlexibleInstances          #-}
+{-# LANGUAGE GeneralizedNewtypeDeriving #-}
+{-# LANGUAGE PatternSynonyms            #-}
+
+#if MIN_VERSION_semirings(0,4,2)
+
+module Data.Poly.Internal.PolyOverField
+  ( PolyOverField(..)
+  , PolyOverFractional
+  , pattern PolyOverFractional
+  , unPolyOverFractional
+  ) where
+
+import Prelude hiding (quotRem, quot, rem, gcd, lcm, (^))
+import Control.DeepSeq (NFData)
+import Data.Euclidean
+import Data.Semiring
+import qualified Data.Vector.Generic as G
+
+import qualified Data.Poly.Internal.Dense as Dense
+import qualified Data.Poly.Internal.Dense.Field as Dense (fieldGcd)
+
+-- | Wrapper for polynomials over 'Field',
+-- providing a faster 'GcdDomain' instance.
+newtype PolyOverField poly = PolyOverField { unPolyOverField :: poly }
+  deriving (Eq, NFData, Num, Ord, Ring, Semiring, Show)
+
+-- |
+type PolyOverFractional = PolyOverField
+{-# DEPRECATED PolyOverFractional "Use 'PolyOverField'" #-}
+
+-- |
+pattern PolyOverFractional :: poly -> PolyOverField poly
+pattern PolyOverFractional poly = PolyOverField poly
+
+-- |
+unPolyOverFractional :: PolyOverField poly -> poly
+unPolyOverFractional = unPolyOverField
+{-# DEPRECATED unPolyOverFractional "Use 'unPolyOverField'" #-}
+
+#if !MIN_VERSION_semirings(0,5,0)
+type Field a = (Euclidean a, Ring a, Fractional a)
+#endif
+
+instance (Eq a, Eq (v a), Field a, G.Vector v a) => GcdDomain (PolyOverField (Dense.Poly v a)) where
+  gcd (PolyOverField x) (PolyOverField y) = PolyOverField (Dense.fieldGcd x y)
+  {-# INLINE gcd #-}
+
+instance (Eq a, Eq (v a), Field a, G.Vector v a) => Euclidean (PolyOverField (Dense.Poly v a)) where
+  degree (PolyOverField x) =
+    degree x
+  quotRem (PolyOverField x) (PolyOverField y) =
+    let (q, r) = quotRem x y in
+      (PolyOverField q, PolyOverField r)
+  {-# INLINE quotRem #-}
+  rem (PolyOverField x) (PolyOverField y) =
+    PolyOverField (rem x y)
+  {-# INLINE rem #-}
+
+#else
+
+module Data.Poly.Internal.PolyOverField () where
+
+#endif
diff --git a/src/Data/Poly/Internal/PolyOverFractional.hs b/src/Data/Poly/Internal/PolyOverFractional.hs
deleted file mode 100644
--- a/src/Data/Poly/Internal/PolyOverFractional.hs
+++ /dev/null
@@ -1,55 +0,0 @@
--- |
--- Module:      Data.Poly.Internal.PolyOverFractional
--- Copyright:   (c) 2019 Andrew Lelechenko
--- Licence:     BSD3
--- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>
---
--- Wrapper with a more efficient 'Euclidean' instance.
---
-
-{-# LANGUAGE CPP                        #-}
-{-# LANGUAGE FlexibleInstances          #-}
-{-# LANGUAGE GeneralizedNewtypeDeriving #-}
-{-# LANGUAGE UndecidableInstances       #-}
-
-#if MIN_VERSION_semirings(0,4,2)
-
-module Data.Poly.Internal.PolyOverFractional
-  ( PolyOverFractional(..)
-  ) where
-
-import Prelude hiding (quotRem, quot, rem, gcd, lcm, (^))
-import Data.Euclidean
-import Data.Semiring
-import qualified Data.Semiring as Semiring
-import qualified Data.Vector.Generic as G
-
-import qualified Data.Poly.Internal.Dense as Dense
-import qualified Data.Poly.Internal.Dense.Fractional as Dense (fractionalGcd)
-
--- | Wrapper over polynomials,
--- providing a faster 'GcdDomain' instance,
--- when coefficients are 'Fractional'.
-newtype PolyOverFractional poly = PolyOverFractional { unPolyOverFractional :: poly }
-  deriving (Eq, Ord, Show, Num, Semiring, Semiring.Ring)
-
-instance (Eq a, Eq (v a), Semiring.Ring a, GcdDomain a, Fractional a, G.Vector v a) => GcdDomain (PolyOverFractional (Dense.Poly v a)) where
-  gcd (PolyOverFractional x) (PolyOverFractional y) = PolyOverFractional (Dense.fractionalGcd x y)
-  {-# INLINE gcd #-}
-
-instance (Eq a, Eq (v a), Semiring.Ring a, GcdDomain a, Fractional a, G.Vector v a) => Euclidean (PolyOverFractional (Dense.Poly v a)) where
-  degree (PolyOverFractional x) =
-    degree x
-  quotRem (PolyOverFractional x) (PolyOverFractional y) =
-    let (q, r) = quotRem x y in
-      (PolyOverFractional q, PolyOverFractional r)
-  {-# INLINE quotRem #-}
-  rem (PolyOverFractional x) (PolyOverFractional y) =
-    PolyOverFractional (rem x y)
-  {-# INLINE rem #-}
-
-#else
-
-module Data.Poly.Internal.PolyOverFractional () where
-
-#endif
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
@@ -7,14 +7,15 @@
 -- Sparse polynomials of one variable.
 --
 
-{-# LANGUAGE CPP                  #-}
-{-# LANGUAGE FlexibleContexts     #-}
-{-# LANGUAGE PatternSynonyms      #-}
-{-# LANGUAGE ScopedTypeVariables  #-}
-{-# LANGUAGE StandaloneDeriving   #-}
-{-# LANGUAGE TypeFamilies         #-}
-{-# LANGUAGE UndecidableInstances #-}
-{-# LANGUAGE ViewPatterns         #-}
+{-# LANGUAGE CPP                        #-}
+{-# LANGUAGE FlexibleContexts           #-}
+{-# LANGUAGE GeneralizedNewtypeDeriving #-}
+{-# LANGUAGE PatternSynonyms            #-}
+{-# LANGUAGE ScopedTypeVariables        #-}
+{-# LANGUAGE StandaloneDeriving         #-}
+{-# LANGUAGE TypeFamilies               #-}
+{-# LANGUAGE UndecidableInstances       #-}
+{-# LANGUAGE ViewPatterns               #-}
 
 module Data.Poly.Internal.Sparse
   ( Poly(..)
@@ -36,14 +37,20 @@
   , pattern X'
   , eval'
   , deriv'
+#if MIN_VERSION_semirings(0,5,0)
+  , integral'
+#endif
   ) where
 
+import Prelude hiding (quot)
+import Control.DeepSeq (NFData)
 import Control.Monad
 import Control.Monad.Primitive
 import Control.Monad.ST
+import Data.Bits
 import Data.List (intersperse)
 import Data.Ord
-import Data.Semiring (Semiring(..))
+import Data.Semiring (Semiring(..), Ring())
 import qualified Data.Semiring as Semiring
 import qualified Data.Vector as V
 import qualified Data.Vector.Generic as G
@@ -55,6 +62,9 @@
 import Data.Semigroup
 import Numeric.Natural
 #endif
+#if MIN_VERSION_semirings(0,5,0)
+import Data.Euclidean (Field, quot)
+#endif
 
 -- | Polynomials of one variable with coefficients from @a@,
 -- backed by a 'G.Vector' @v@ (boxed, unboxed, storable, etc.).
@@ -78,8 +88,9 @@
   -- (first element corresponds to a constant term).
   }
 
-deriving instance Eq   (v (Word, a)) => Eq   (Poly v a)
-deriving instance Ord  (v (Word, a)) => Ord  (Poly v a)
+deriving instance Eq     (v (Word, a)) => Eq     (Poly v a)
+deriving instance Ord    (v (Word, a)) => Ord    (Poly v a)
+deriving instance NFData (v (Word, a)) => NFData (Poly v a)
 
 instance (Eq a, Semiring a, G.Vector v (Word, a)) => IsList (Poly v a) where
   type Item (Poly v a) = (Word, a)
@@ -143,7 +154,7 @@
   | otherwise = runST $ do
     ws <- G.thaw vs
     l' <- normalizeM p add ws
-    G.unsafeFreeze $ MG.basicUnsafeSlice 0 l' ws
+    G.unsafeFreeze $ MG.unsafeSlice 0 l' ws
 
 normalizeM
   :: (PrimMonad m, G.Vector v (Word, a))
@@ -152,9 +163,9 @@
   -> G.Mutable v (PrimState m) (Word, a)
   -> m Int
 normalizeM p add ws = do
-    let l = MG.basicLength ws
+    let l = MG.length ws
     let go i j acc@(accP, accC)
-          | j >= l = do
+          | j >= l =
             if p accC
               then do
                 MG.write ws i acc
@@ -181,7 +192,7 @@
   abs = id
   signum = const 1
   fromInteger n = case fromInteger n of
-    0 -> Poly $ G.empty
+    0 -> Poly G.empty
     m -> Poly $ G.singleton (0, m)
   Poly xs * Poly ys = Poly $ convolution (/= 0) (+) (*) xs ys
   {-# INLINE (+) #-}
@@ -210,7 +221,7 @@
   {-# INLINE fromNatural #-}
 #endif
 
-instance (Eq a, Semiring.Ring a, G.Vector v (Word, a)) => Semiring.Ring (Poly v a) where
+instance (Eq a, Ring a, G.Vector v (Word, a)) => Ring (Poly v a) where
   negate (Poly xs) = Poly $ G.map (fmap Semiring.negate) xs
 
 plusPoly
@@ -221,9 +232,9 @@
   -> v (Word, a)
   -> v (Word, a)
 plusPoly p add xs ys = runST $ do
-  zs <- MG.basicUnsafeNew (G.basicLength xs + G.basicLength ys)
+  zs <- MG.unsafeNew (G.length xs + G.length ys)
   lenZs <- plusPolyM p add xs ys zs
-  G.unsafeFreeze $ MG.basicUnsafeSlice 0 lenZs zs
+  G.unsafeFreeze $ MG.unsafeSlice 0 lenZs zs
 {-# INLINE plusPoly #-}
 
 plusPolyM
@@ -236,20 +247,20 @@
   -> m Int
 plusPolyM p add xs ys zs = go 0 0 0
   where
-    lenXs = G.basicLength xs
-    lenYs = G.basicLength ys
+    lenXs = G.length xs
+    lenYs = G.length ys
 
     go ix iy iz
       | ix == lenXs, iy == lenYs = pure iz
       | ix == lenXs = do
         G.unsafeCopy
-          (MG.basicUnsafeSlice iz (lenYs - iy) zs)
-          (G.basicUnsafeSlice iy (lenYs - iy) ys)
+          (MG.unsafeSlice iz (lenYs - iy) zs)
+          (G.unsafeSlice iy (lenYs - iy) ys)
         pure $ iz + lenYs - iy
       | iy == lenYs = do
         G.unsafeCopy
-          (MG.basicUnsafeSlice iz (lenXs - ix) zs)
-          (G.basicUnsafeSlice ix (lenXs - ix) xs)
+          (MG.unsafeSlice iz (lenXs - ix) zs)
+          (G.unsafeSlice ix (lenXs - ix) xs)
         pure $ iz + lenXs - ix
       | (xp, xc) <- G.unsafeIndex xs ix
       , (yp, yc) <- G.unsafeIndex ys iy
@@ -278,7 +289,7 @@
   -> v (Word, a)
   -> v (Word, a)
 minusPoly p neg sub xs ys = runST $ do
-  zs <- MG.basicUnsafeNew (lenXs + lenYs)
+  zs <- MG.unsafeNew (lenXs + lenYs)
   let go ix iy iz
         | ix == lenXs, iy == lenYs = pure iz
         | ix == lenXs = do
@@ -288,8 +299,8 @@
           pure $ iz + lenYs - iy
         | iy == lenYs = do
           G.unsafeCopy
-            (MG.basicUnsafeSlice iz (lenXs - ix) zs)
-            (G.basicUnsafeSlice ix (lenXs - ix) xs)
+            (MG.unsafeSlice iz (lenXs - ix) zs)
+            (G.unsafeSlice ix (lenXs - ix) xs)
           pure $ iz + lenXs - ix
         | (xp, xc) <- G.unsafeIndex xs ix
         , (yp, yc) <- G.unsafeIndex ys iy
@@ -308,10 +319,10 @@
             MG.unsafeWrite zs iz (yp, neg yc)
             go ix (iy + 1) (iz + 1)
   lenZs <- go 0 0 0
-  G.unsafeFreeze $ MG.basicUnsafeSlice 0 lenZs zs
+  G.unsafeFreeze $ MG.unsafeSlice 0 lenZs zs
   where
-    lenXs = G.basicLength xs
-    lenYs = G.basicLength ys
+    lenXs = G.length xs
+    lenYs = G.length ys
 {-# INLINE minusPoly #-}
 
 scaleM
@@ -324,7 +335,7 @@
   -> m Int
 scaleM p mul xs (yp, yc) zs = go 0 0
   where
-    lenXs = G.basicLength xs
+    lenXs = G.length xs
 
     go ix iz
       | ix == lenXs = pure iz
@@ -347,9 +358,9 @@
   -> Poly v a
   -> Poly v a
 scaleInternal p mul yp yc (Poly xs) = runST $ do
-  zs <- MG.basicUnsafeNew (G.basicLength xs)
+  zs <- MG.unsafeNew (G.length xs)
   len <- scaleM p (flip mul) xs (yp, yc) zs
-  fmap Poly $ G.unsafeFreeze $ MG.basicUnsafeSlice 0 len zs
+  fmap Poly $ G.unsafeFreeze $ MG.unsafeSlice 0 len zs
 {-# INLINE scaleInternal #-}
 
 -- | Multiply a polynomial by a monomial, expressed as a power and a coefficient.
@@ -372,29 +383,29 @@
   -> v (Word, a)
   -> v (Word, a)
 convolution p add mult xs ys
-  | G.basicLength xs >= G.basicLength ys
+  | G.length xs >= G.length ys
   = go mult xs ys
   | otherwise
   = go (flip mult) ys xs
   where
     go :: (a -> a -> a) -> v (Word, a) -> v (Word, a) -> v (Word, a)
     go mul long short = runST $ do
-      let lenLong   = G.basicLength long
-          lenShort  = G.basicLength short
+      let lenLong   = G.length long
+          lenShort  = G.length short
           lenBuffer = lenLong * lenShort
-      slices <- MG.basicUnsafeNew lenShort
-      buffer <- MG.basicUnsafeNew lenBuffer
+      slices <- MG.unsafeNew lenShort
+      buffer <- MG.unsafeNew lenBuffer
 
       forM_ [0 .. lenShort - 1] $ \iShort -> do
         let (pShort, cShort) = G.unsafeIndex short iShort
             from = iShort * lenLong
-            bufferSlice = MG.basicUnsafeSlice from lenLong buffer
+            bufferSlice = MG.unsafeSlice from lenLong buffer
         len <- scaleM p mul long (pShort, cShort) bufferSlice
         MG.unsafeWrite slices iShort (from, len)
 
       slices' <- G.unsafeFreeze slices
       buffer' <- G.unsafeFreeze buffer
-      bufferNew <- MG.basicUnsafeNew lenBuffer
+      bufferNew <- MG.unsafeNew lenBuffer
       gogo slices' buffer' bufferNew
 
     gogo
@@ -404,29 +415,29 @@
       -> G.Mutable v (PrimState m) (Word, a)
       -> m (v (Word, a))
     gogo slices buffer bufferNew
-      | G.basicLength slices == 0
+      | G.length slices == 0
       = pure G.empty
-      | G.basicLength slices == 1
+      | G.length slices == 1
       , (from, len) <- G.unsafeIndex slices 0
-      = pure $ G.basicUnsafeSlice from len buffer
+      = pure $ G.unsafeSlice from len buffer
       | otherwise = do
-        let nSlices = G.basicLength slices
-        slicesNew <- MG.basicUnsafeNew ((nSlices + 1) `quot` 2)
-        forM_ [0 .. (nSlices - 2) `quot` 2] $ \i -> do
+        let nSlices = G.length slices
+        slicesNew <- MG.unsafeNew ((nSlices + 1) `shiftR` 1)
+        forM_ [0 .. (nSlices - 2) `shiftR` 1] $ \i -> do
           let (from1, len1) = G.unsafeIndex slices (2 * i)
               (from2, len2) = G.unsafeIndex slices (2 * i + 1)
-              slice1 = G.basicUnsafeSlice from1 len1 buffer
-              slice2 = G.basicUnsafeSlice from2 len2 buffer
-              slice3 = MG.basicUnsafeSlice from1 (len1 + len2) bufferNew
+              slice1 = G.unsafeSlice from1 len1 buffer
+              slice2 = G.unsafeSlice from2 len2 buffer
+              slice3 = MG.unsafeSlice from1 (len1 + len2) bufferNew
           len3 <- plusPolyM p add slice1 slice2 slice3
           MG.unsafeWrite slicesNew i (from1, len3)
 
         when (odd nSlices) $ do
           let (from, len) = G.unsafeIndex slices (nSlices - 1)
-              slice1 = G.basicUnsafeSlice from len buffer
-              slice3 = MG.basicUnsafeSlice from len bufferNew
+              slice1 = G.unsafeSlice from len buffer
+              slice3 = MG.unsafeSlice from len bufferNew
           G.unsafeCopy slice3 slice1
-          MG.unsafeWrite slicesNew (nSlices `quot` 2) (from, len)
+          MG.unsafeWrite slicesNew (nSlices `shiftR` 1) (from, len)
 
         slicesNew' <- G.unsafeFreeze slicesNew
         buffer'    <- G.unsafeThaw   buffer
@@ -509,8 +520,8 @@
 derivPoly p mul xs
   | G.null xs = G.empty
   | otherwise = runST $ do
-    let lenXs = G.basicLength xs
-    zs <- MG.basicUnsafeNew lenXs
+    let lenXs = G.length xs
+    zs <- MG.unsafeNew lenXs
     let go ix iz
           | ix == lenXs = pure iz
           | (xp, xc) <- G.unsafeIndex xs ix
@@ -522,7 +533,7 @@
             else
               go (ix + 1) iz
     lenZs <- go 0 0
-    G.unsafeFreeze $ MG.basicUnsafeSlice 0 lenZs zs
+    G.unsafeFreeze $ MG.unsafeSlice 0 lenZs zs
 {-# INLINE derivPoly #-}
 
 -- | Compute an indefinite integral of a polynomial,
@@ -535,6 +546,14 @@
   = Poly
   $ G.map (\(p, c) -> (p + 1, c / (fromIntegral p + 1))) xs
 {-# INLINE integral #-}
+
+#if MIN_VERSION_semirings(0,5,0)
+integral' :: (Eq a, Field a, G.Vector v (Word, a)) => Poly v a -> Poly v a
+integral' (Poly xs)
+  = Poly
+  $ G.map (\(p, c) -> (p + 1, c `quot` Semiring.fromIntegral (p + 1))) xs
+{-# INLINE integral' #-}
+#endif
 
 -- | Create an identity polynomial.
 pattern X :: (Eq a, Num a, G.Vector v (Word, a), Eq (v (Word, a))) => Poly v a
diff --git a/src/Data/Poly/Internal/Sparse/Field.hs b/src/Data/Poly/Internal/Sparse/Field.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Poly/Internal/Sparse/Field.hs
@@ -0,0 +1,118 @@
+-- |
+-- Module:      Data.Poly.Internal.Sparse.Field
+-- Copyright:   (c) 2019 Andrew Lelechenko
+-- Licence:     BSD3
+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>
+--
+-- GcdDomain for Field underlying
+--
+
+{-# LANGUAGE ConstraintKinds            #-}
+{-# LANGUAGE CPP                        #-}
+{-# LANGUAGE FlexibleContexts           #-}
+{-# LANGUAGE FlexibleInstances          #-}
+{-# LANGUAGE PatternSynonyms            #-}
+{-# LANGUAGE ScopedTypeVariables        #-}
+{-# LANGUAGE TypeFamilies               #-}
+{-# LANGUAGE UndecidableInstances       #-}
+
+{-# OPTIONS_GHC -fno-warn-orphans #-}
+
+#if MIN_VERSION_semirings(0,4,2)
+
+module Data.Poly.Internal.Sparse.Field
+  ( gcdExt
+  ) where
+
+import Prelude hiding (quotRem, quot, rem, gcd)
+import Control.Arrow
+import Control.Exception
+import Data.Euclidean
+#if !MIN_VERSION_semirings(0,5,0)
+import Data.Semiring (Ring)
+#endif
+import Data.Semiring (minus, plus, times, zero, one)
+import qualified Data.Vector.Generic as G
+
+import Data.Poly.Internal.Sparse
+import Data.Poly.Internal.Sparse.GcdDomain ()
+
+#if !MIN_VERSION_semirings(0,5,0)
+type Field a = (Euclidean a, Ring a, Fractional a)
+#endif
+
+instance (Eq a, Eq (v (Word, a)), Field a, G.Vector v (Word, a)) => Euclidean (Poly v a) where
+  degree (Poly xs)
+    | G.null xs = 0
+    | otherwise = 1 + fromIntegral (fst (G.last xs))
+
+  quotRem = quotientRemainder
+
+quotientRemainder
+  :: (Eq a, Field a, G.Vector v (Word, a))
+  => Poly v a
+  -> Poly v a
+  -> (Poly v a, Poly v a)
+quotientRemainder ts ys = case leading ys of
+  Nothing -> throw DivideByZero
+  Just (yp, yc) -> go ts
+    where
+      go xs = case leading xs of
+        Nothing -> (zero, zero)
+        Just (xp, xc) -> case xp `compare` yp of
+          LT -> (zero, xs)
+          EQ -> (zs, xs')
+          GT -> first (`plus` zs) $ go xs'
+          where
+            zs = Poly $ G.singleton (xp `minus` yp, xc `quot` yc)
+            xs' = xs `minus` zs `times` ys
+
+-- | Execute the extended Euclidean algorithm.
+-- For polynomials @a@ and @b@, compute their unique greatest common divisor @g@
+-- and the unique coefficient polynomial @s@ satisfying @as + bt = g@,
+-- such that either @g@ is monic, or @g = 0@ and @s@ is monic, or @g = s = 0@.
+--
+-- >>> gcdExt (X^2 + 1 :: UPoly Double) (X^3 + 3 * X :: UPoly Double)
+-- (1.0, 0.5 * X^2 + (-0.0) * X + 1.0)
+-- >>> gcdExt (X^3 + 3 * X :: UPoly Double) (3 * X^4 + 3 * X^2 :: UPoly Double)
+-- (1.0 * X + 0.0,(-0.16666666666666666) * X^2 + (-0.0) * X + 0.3333333333333333)
+gcdExt
+  :: (Eq a, Field a, G.Vector v (Word, a), Eq (v (Word, a)))
+  => Poly v a
+  -> Poly v a
+  -> (Poly v a, Poly v a)
+gcdExt xs ys = case scaleMonic gs of
+  Just (c', gs') -> (gs', scale' zero c' ss)
+  Nothing -> case scaleMonic ss of
+    Just (_, ss') -> (zero, ss')
+    Nothing -> (zero, zero)
+  where
+    (gs, ss) = go ys xs zero one
+      where
+        go r' r s' s
+          | r' == zero = (r, s)
+          | otherwise  = case r `quotRem` r' of
+            (q, r'') -> go r'' r' (s `minus` q `times` s') s'
+{-# INLINE gcdExt #-}
+
+-- | Scale a non-zero polynomial such that its leading coefficient is one,
+-- returning the reciprocal of the leading coefficient in the scaling.
+--
+-- >>> scaleMonic (X^3 + 3 * X :: UPoly Double)
+-- Just (1.0, 1.0 * X^3 + 0.0 * X^2 + 3.0 * X + 0.0)
+-- >>> scaleMonic (3 * X^4 + 3 * X^2 :: UPoly Double)
+-- Just (0.3333333333333333, 1.0 * X^4 + 0.0 * X^3 + 1.0 * X^2 + 0.0 * X + 0.0)
+scaleMonic
+  :: (Eq a, Field a, G.Vector v (Word, a), Eq (v (Word, a)))
+  => Poly v a
+  -> Maybe (a, Poly v a)
+scaleMonic xs = case leading xs of
+  Nothing -> Nothing
+  Just (_, c) -> let c' = one `quot` c in Just (c', scale' zero c' xs)
+{-# INLINE scaleMonic #-}
+
+#else
+
+module Data.Poly.Internal.Sparse.Field () where
+
+#endif
diff --git a/src/Data/Poly/Internal/Sparse/Fractional.hs b/src/Data/Poly/Internal/Sparse/Fractional.hs
deleted file mode 100644
--- a/src/Data/Poly/Internal/Sparse/Fractional.hs
+++ /dev/null
@@ -1,78 +0,0 @@
--- |
--- Module:      Data.Poly.Internal.Sparse.Fractional
--- Copyright:   (c) 2019 Andrew Lelechenko
--- Licence:     BSD3
--- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>
---
--- GcdDomain for Fractional underlying
---
-
-{-# LANGUAGE CPP                        #-}
-{-# LANGUAGE FlexibleContexts           #-}
-{-# LANGUAGE FlexibleInstances          #-}
-{-# LANGUAGE GeneralizedNewtypeDeriving #-}
-{-# LANGUAGE PatternSynonyms            #-}
-{-# LANGUAGE ScopedTypeVariables        #-}
-{-# LANGUAGE TypeFamilies               #-}
-{-# LANGUAGE UndecidableInstances       #-}
-{-# LANGUAGE ViewPatterns               #-}
-
-{-# OPTIONS_GHC -fno-warn-orphans #-}
-
-#if MIN_VERSION_semirings(0,4,2)
-
-module Data.Poly.Internal.Sparse.Fractional
-  ( fractionalGcd
-  ) where
-
-import Prelude hiding (quotRem, rem, gcd)
-import Control.Arrow
-import Control.Exception
-import Data.Euclidean
-import qualified Data.Semiring as Semiring
-import qualified Data.Vector.Generic as G
-
-import Data.Poly.Internal.Sparse
-import Data.Poly.Internal.Sparse.GcdDomain ()
-
-instance (Eq a, Eq (v (Word, a)), Semiring.Ring a, GcdDomain a, Fractional a, G.Vector v (Word, a)) => Euclidean (Poly v a) where
-  degree (Poly xs)
-    | G.null xs = 0
-    | otherwise = 1 + fromIntegral (fst (G.last xs))
-
-  quotRem = quotientRemainder
-
-quotientRemainder
-  :: (Eq a, Fractional a, G.Vector v (Word, a))
-  => Poly v a
-  -> Poly v a
-  -> (Poly v a, Poly v a)
-quotientRemainder ts ys = case leading ys of
-  Nothing -> throw DivideByZero
-  Just (yp, yc) -> go ts
-    where
-      go xs = case leading xs of
-        Nothing -> (0, 0)
-        Just (xp, xc) -> case xp `compare` yp of
-          LT -> (0, xs)
-          EQ -> (zs, xs')
-          GT -> first (+ zs) $ go xs'
-          where
-            zs = Poly $ G.singleton (xp - yp, xc / yc)
-            xs' = xs - zs * ys
-
-fractionalGcd
-  :: (Eq a, Fractional a, G.Vector v (Word, a))
-  => Poly v a
-  -> Poly v a
-  -> Poly v a
-fractionalGcd xs ys
-  | G.null (unPoly ys) = xs
-  | otherwise = fractionalGcd ys $ snd $ quotientRemainder xs ys
-{-# INLINE fractionalGcd #-}
-
-#else
-
-module Data.Poly.Internal.Sparse.Fractional () where
-
-#endif
diff --git a/src/Data/Poly/Internal/Sparse/GcdDomain.hs b/src/Data/Poly/Internal/Sparse/GcdDomain.hs
--- a/src/Data/Poly/Internal/Sparse/GcdDomain.hs
+++ b/src/Data/Poly/Internal/Sparse/GcdDomain.hs
@@ -10,12 +10,10 @@
 {-# LANGUAGE CPP                        #-}
 {-# LANGUAGE FlexibleContexts           #-}
 {-# LANGUAGE FlexibleInstances          #-}
-{-# LANGUAGE GeneralizedNewtypeDeriving #-}
 {-# LANGUAGE PatternSynonyms            #-}
 {-# LANGUAGE ScopedTypeVariables        #-}
 {-# LANGUAGE TypeFamilies               #-}
 {-# LANGUAGE UndecidableInstances       #-}
-{-# LANGUAGE ViewPatterns               #-}
 
 {-# OPTIONS_GHC -fno-warn-orphans #-}
 
@@ -28,17 +26,15 @@
 import Control.Exception
 import Data.Euclidean
 import Data.Maybe
-import Data.Semiring (Semiring(..))
-import qualified Data.Semiring as Semiring
+import Data.Semiring (Semiring(..), Ring(), minus)
 import qualified Data.Vector.Generic as G
 
 import Data.Poly.Internal.Sparse
 
--- | Consider using 'Data.Poly.Sparse.Semiring.PolyOverFractional' wrapper,
+-- | Consider using 'Data.Poly.Sparse.Semiring.PolyOverField' wrapper,
 -- which provides a much faster implementation of
--- 'Data.Euclidean.gcd' for 'Fractional'
--- coefficients.
-instance (Eq a, Semiring.Ring a, GcdDomain a, Eq (v (Word, a)), G.Vector v (Word, a)) => GcdDomain (Poly v a) where
+-- 'Data.Euclidean.gcd' for polynomials over 'Field'.
+instance (Eq a, Ring a, GcdDomain a, Eq (v (Word, a)), G.Vector v (Word, a)) => GcdDomain (Poly v a) where
   divide xs ys = case leading ys of
     Nothing -> throw DivideByZero
     Just (yp, yc) -> case leading xs of
@@ -48,7 +44,7 @@
         | otherwise -> do
           zc <- divide xc yc
           let z = Poly $ G.singleton (xp - yp, zc)
-          rest <- divide (xs `plus` Semiring.negate z `times` ys) ys
+          rest <- divide (xs `minus` z `times` ys) ys
           pure $ rest `plus` z
 
   gcd xs ys
@@ -62,7 +58,7 @@
         xy = monomial' 0 (gcd (cont xs) (cont ys))
 
 gcdHelper
-  :: (Eq a, Semiring.Ring a, GcdDomain a, G.Vector v (Word, a))
+  :: (Eq a, Ring a, GcdDomain a, G.Vector v (Word, a))
   => Poly v a
   -> Poly v a
   -> Poly v a
@@ -71,9 +67,9 @@
   Just (xp, xc) -> case leading ys of
     Nothing -> xs
     Just (yp, yc) -> case xp `compare` yp of
-      LT -> gcdHelper xs (ys `times` monomial' 0 gy `plus` Semiring.negate (xs `times` monomial' (yp - xp) gx))
-      EQ -> gcdHelper xs (ys `times` monomial' 0 gy `plus` Semiring.negate (xs `times` monomial' 0 gx))
-      GT -> gcdHelper (xs `times` monomial' 0 gx `plus` Semiring.negate (ys `times` monomial' (xp - yp) gy)) ys
+      LT -> gcdHelper xs (ys `times` monomial' 0 gy `minus` xs `times` monomial' (yp - xp) gx)
+      EQ -> gcdHelper xs (ys `times` monomial' 0 gy `minus` xs `times` monomial' 0 gx)
+      GT -> gcdHelper (xs `times` monomial' 0 gx `minus` ys `times` monomial' (xp - yp) gy) ys
       where
         g = lcm xc yc
         gx = fromMaybe err $ divide g xc
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
@@ -7,8 +7,8 @@
 -- Dense polynomials and a 'Semiring'-based interface.
 --
 
-{-# LANGUAGE CPP                 #-}
-{-# LANGUAGE PatternSynonyms     #-}
+{-# LANGUAGE CPP             #-}
+{-# LANGUAGE PatternSynonyms #-}
 
 module Data.Poly.Semiring
   ( Poly
@@ -23,9 +23,16 @@
   , pattern X
   , eval
   , deriv
+#if MIN_VERSION_semirings(0,5,0)
+  , integral
+#endif
 #if MIN_VERSION_semirings(0,4,2)
-  -- * Fractional coefficients
-  , PolyOverFractional(..)
+  -- * Polynomials over 'Field'
+  , PolyOverField(..)
+  , gcdExt
+  , PolyOverFractional
+  , pattern PolyOverFractional
+  , unPolyOverFractional
 #endif
   ) where
 
@@ -35,10 +42,13 @@
 import Data.Poly.Internal.Dense (Poly(..), VPoly, UPoly, leading)
 import qualified Data.Poly.Internal.Dense as Dense
 #if MIN_VERSION_semirings(0,4,2)
-import Data.Poly.Internal.Dense.Fractional ()
+import Data.Poly.Internal.Dense.Field (gcdExt)
 import Data.Poly.Internal.Dense.GcdDomain ()
-import Data.Poly.Internal.PolyOverFractional
+import Data.Poly.Internal.PolyOverField
 #endif
+#if MIN_VERSION_semirings(0,5,0)
+import Data.Euclidean (Field)
+#endif
 
 -- | Make 'Poly' from a vector of coefficients
 -- (first element corresponds to a constant term).
@@ -81,3 +91,13 @@
 -- 3 * X^2 + 0 * X + 3
 deriv :: (Eq a, Semiring a, G.Vector v a) => Poly v a -> Poly v a
 deriv = Dense.deriv'
+
+#if MIN_VERSION_semirings(0,5,0)
+-- | Compute an indefinite integral of a polynomial,
+-- setting constant term to zero.
+--
+-- >>> integral (3 * X^2 + 3) :: UPoly Double
+-- 1.0 * X^3 + 0.0 * X^2 + 3.0 * X + 0.0
+integral :: (Eq a, Field a, G.Vector v a) => Poly v a -> Poly v a
+integral = Dense.integral'
+#endif
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
@@ -7,7 +7,8 @@
 -- Sparse polynomials with 'Num' instance.
 --
 
-{-# LANGUAGE PatternSynonyms     #-}
+{-# LANGUAGE CPP             #-}
+{-# LANGUAGE PatternSynonyms #-}
 
 module Data.Poly.Sparse
   ( Poly
@@ -23,8 +24,14 @@
   , eval
   , deriv
   , integral
+#if MIN_VERSION_semirings(0,4,2)
+  -- * Polynomials over 'Field'
+  , gcdExt
+#endif
   ) where
 
 import Data.Poly.Internal.Sparse
-import Data.Poly.Internal.Sparse.Fractional ()
+#if MIN_VERSION_semirings(0,4,2)
+import Data.Poly.Internal.Sparse.Field (gcdExt)
 import Data.Poly.Internal.Sparse.GcdDomain ()
+#endif
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
@@ -7,8 +7,9 @@
 -- Sparse polynomials with 'Semiring' instance.
 --
 
-{-# LANGUAGE FlexibleContexts    #-}
-{-# LANGUAGE PatternSynonyms     #-}
+{-# LANGUAGE CPP              #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE PatternSynonyms  #-}
 
 module Data.Poly.Sparse.Semiring
   ( Poly
@@ -23,6 +24,13 @@
   , pattern X
   , eval
   , deriv
+#if MIN_VERSION_semirings(0,5,0)
+  , integral
+#endif
+#if MIN_VERSION_semirings(0,4,2)
+  -- * Polynomials over 'Field'
+  , gcdExt
+#endif
   ) where
 
 import Data.Semiring (Semiring)
@@ -30,8 +38,13 @@
 
 import Data.Poly.Internal.Sparse (Poly(..), VPoly, UPoly, leading)
 import qualified Data.Poly.Internal.Sparse as Sparse
-import Data.Poly.Internal.Sparse.Fractional ()
+#if MIN_VERSION_semirings(0,4,2)
+import Data.Poly.Internal.Sparse.Field (gcdExt)
 import Data.Poly.Internal.Sparse.GcdDomain ()
+#endif
+#if MIN_VERSION_semirings(0,5,0)
+import Data.Euclidean (Field)
+#endif
 
 -- | Make 'Poly' from a list of (power, coefficient) pairs.
 -- (first element corresponds to a constant term).
@@ -74,3 +87,13 @@
 -- 3 * X^2 + 3
 deriv :: (Eq a, Semiring a, G.Vector v (Word, a)) => Poly v a -> Poly v a
 deriv = Sparse.deriv'
+
+#if MIN_VERSION_semirings(0,5,0)
+-- | Compute an indefinite integral of a polynomial,
+-- setting constant term to zero.
+--
+-- >>> integral (3 * X^2 + 3) :: UPoly Double
+-- 1.0 * X^3 + 3.0 * X
+integral :: (Eq a, Field a, G.Vector v (Word, a)) => Poly v a -> Poly v a
+integral = Sparse.integral'
+#endif
diff --git a/test/Dense.hs b/test/Dense.hs
--- a/test/Dense.hs
+++ b/test/Dense.hs
@@ -10,11 +10,12 @@
   ( testSuite
   ) where
 
-import Prelude hiding (quotRem)
+import Prelude hiding (gcd, quotRem, rem)
 #if MIN_VERSION_semirings(0,4,2)
 import Data.Euclidean
 #endif
 import Data.Int
+import Data.Maybe
 import Data.Poly
 import qualified Data.Poly.Semiring as S
 import Data.Proxy
@@ -23,7 +24,7 @@
 import qualified Data.Vector.Generic as G
 import qualified Data.Vector.Unboxed as U
 import Test.Tasty
-import Test.Tasty.QuickCheck hiding (scale)
+import Test.Tasty.QuickCheck hiding (scale, numTests)
 import Test.QuickCheck.Classes
 
 import Quaternion
@@ -33,9 +34,9 @@
   shrink = fmap (S.toPoly . G.fromList) . shrink . G.toList . unPoly
 
 #if MIN_VERSION_semirings(0,4,2)
-instance (Eq a, Semiring a, Arbitrary a, G.Vector v a) => Arbitrary (PolyOverFractional (Poly v a)) where
-  arbitrary = PolyOverFractional . S.toPoly . G.fromList . (\xs -> take (length xs `mod` 10) xs) <$> arbitrary
-  shrink = fmap (PolyOverFractional . S.toPoly . G.fromList) . shrink . G.toList . unPoly . unPolyOverFractional
+instance (Eq a, Semiring a, Arbitrary a, G.Vector v a) => Arbitrary (PolyOverField (Poly v a)) where
+  arbitrary = PolyOverField . S.toPoly . G.fromList . (\xs -> take (length xs `mod` 10) xs) <$> arbitrary
+  shrink = fmap (PolyOverField . S.toPoly . G.fromList) . shrink . G.toList . unPoly . unPolyOverField
 #endif
 
 newtype ShortPoly a = ShortPoly { unShortPoly :: a }
@@ -57,14 +58,24 @@
 testSuite = testGroup "Dense"
     [ arithmeticTests
     , otherTests
-    , semiringTests
+    , lawsTests
     , evalTests
     , derivTests
 #if MIN_VERSION_semirings(0,4,2)
-    -- , euclideanTests
+    , gcdExtTests
 #endif
     ]
 
+lawsTests :: TestTree
+lawsTests = testGroup "Laws"
+  [ semiringTests
+  , ringTests
+  , numTests
+  , euclideanTests
+  , isListTests
+  , showTests
+  ]
+
 semiringTests :: TestTree
 semiringTests
   = testGroup "Semiring"
@@ -74,26 +85,73 @@
   , semiringLaws (Proxy :: Proxy (Poly U.Vector Int8))
   , semiringLaws (Proxy :: Proxy (Poly V.Vector Integer))
   , semiringLaws (Proxy :: Proxy (Poly U.Vector (Quaternion Int)))
+  ]
+
+ringTests :: TestTree
+ringTests
+  = testGroup "Ring"
+  $ map (uncurry testProperty)
+  $ concatMap lawsProperties
+  [
 #if MIN_VERSION_quickcheck_classes(0,6,1)
-  , ringLaws (Proxy :: Proxy (Poly U.Vector ()))
+    ringLaws (Proxy :: Proxy (Poly U.Vector ()))
   , ringLaws (Proxy :: Proxy (Poly U.Vector Int8))
   , ringLaws (Proxy :: Proxy (Poly V.Vector Integer))
   , ringLaws (Proxy :: Proxy (Poly U.Vector (Quaternion Int)))
 #endif
   ]
 
-#if MIN_VERSION_semirings(0,4,2)
--- euclideanTests :: TestTree
--- euclideanTests
---   = testGroup "Euclidean"
---   $ map (uncurry testProperty)
---   $ concatMap lawsProperties
---   [ gcdDomainLaws (Proxy :: Proxy (ShortPoly (Poly V.Vector Integer)))
---   , gcdDomainLaws (Proxy :: Proxy (PolyOverFractional (Poly V.Vector Rational)))
---   , euclideanLaws (Proxy :: Proxy (ShortPoly (Poly V.Vector Rational)))
---   ]
+numTests :: TestTree
+numTests
+  = testGroup "Num"
+  $ map (uncurry testProperty)
+  $ concatMap lawsProperties
+  [
+#if MIN_VERSION_quickcheck_classes(0,6,3)
+    numLaws (Proxy :: Proxy (Poly U.Vector Int8))
+  , numLaws (Proxy :: Proxy (Poly V.Vector Integer))
+  , numLaws (Proxy :: Proxy (Poly U.Vector (Quaternion Int)))
 #endif
+  ]
 
+euclideanTests :: TestTree
+euclideanTests
+  = testGroup "Euclidean"
+  $ map (uncurry testProperty)
+  $ concatMap lawsProperties
+  [
+#if MIN_VERSION_semirings(0,4,2) && MIN_VERSION_quickcheck_classes(0,6,3)
+    gcdDomainLaws (Proxy :: Proxy (ShortPoly (Poly V.Vector Integer)))
+  , gcdDomainLaws (Proxy :: Proxy (PolyOverField (Poly V.Vector Rational)))
+  , euclideanLaws (Proxy :: Proxy (ShortPoly (Poly V.Vector Rational)))
+#endif
+  ]
+
+isListTests :: TestTree
+isListTests
+  = testGroup "IsList"
+  $ map (uncurry testProperty)
+  $ concatMap lawsProperties
+  [ isListLaws (Proxy :: Proxy (Poly U.Vector ()))
+  , isListLaws (Proxy :: Proxy (Poly U.Vector Int8))
+  , isListLaws (Proxy :: Proxy (Poly V.Vector Integer))
+  , isListLaws (Proxy :: Proxy (Poly U.Vector (Quaternion Int)))
+  ]
+
+showTests :: TestTree
+showTests
+  = testGroup "Show"
+  $ map (uncurry testProperty)
+  $ concatMap lawsProperties
+  [
+#if MIN_VERSION_quickcheck_classes(0,6,0)
+    showLaws (Proxy :: Proxy (Poly U.Vector ()))
+  , showLaws (Proxy :: Proxy (Poly U.Vector Int8))
+  , showLaws (Proxy :: Proxy (Poly V.Vector Integer))
+  , showLaws (Proxy :: Proxy (Poly U.Vector (Quaternion Int)))
+#endif
+  ]
+
 arithmeticTests :: TestTree
 arithmeticTests = testGroup "Arithmetic"
   [ testProperty "addition matches reference" $
@@ -189,6 +247,8 @@
 derivTests = testGroup "deriv"
   [ testProperty "deriv = S.deriv" $
     \(p :: Poly V.Vector Integer) -> deriv p === S.deriv p
+  , testProperty "integral = S.integral" $
+    \(p :: Poly V.Vector Rational) -> integral p === S.integral p
   , testProperty "deriv . integral = id" $
     \(p :: Poly V.Vector Rational) -> deriv (integral p) === p
   , testProperty "deriv c = 0" $
@@ -205,3 +265,22 @@
   --     deriv (eval (toPoly $ fmap (monomial 0) $ unPoly p) q) ===
   --       deriv q * eval (toPoly $ fmap (monomial 0) $ unPoly $ deriv p) q
   ]
+
+#if MIN_VERSION_semirings(0,4,2)
+gcdExtTests :: TestTree
+gcdExtTests = localOption (QuickCheckMaxSize 12) $ testGroup "gcdExt"
+  [ testProperty "gcdExt == S.gcdExt" $
+    \(a :: Poly V.Vector Rational) b ->
+      gcdExt a b === S.gcdExt a b
+  , testProperty "g == as (mod b) for gcdExt" $
+    \(a :: Poly V.Vector Rational) b -> b /= 0 ==>
+      uncurry ((. flip rem b) . (===) . flip rem b) ((* a) <$> gcdExt a b)
+  , testProperty "fst . gcdExt == gcd (mod units)" $
+    \(a :: Poly V.Vector Rational) b ->
+      fst (gcdExt a b) `sameUpToUnits` gcd a b
+  ]
+
+sameUpToUnits :: (Eq a, GcdDomain a) => a -> a -> Bool
+sameUpToUnits x y = x == y ||
+  isJust (x `divide` y) && isJust (y `divide` x)
+#endif
diff --git a/test/Sparse.hs b/test/Sparse.hs
--- a/test/Sparse.hs
+++ b/test/Sparse.hs
@@ -11,13 +11,14 @@
   ( testSuite
   ) where
 
-import Prelude hiding (quotRem)
+import Prelude hiding (gcd, quotRem, rem)
 #if MIN_VERSION_semirings(0,4,2)
 import Data.Euclidean
 #endif
 import Data.Function
 import Data.Int
 import Data.List
+import Data.Maybe
 import Data.Poly.Sparse
 import qualified Data.Poly.Sparse.Semiring as S
 import Data.Proxy
@@ -26,7 +27,7 @@
 import qualified Data.Vector.Generic as G
 import qualified Data.Vector.Unboxed as U
 import Test.Tasty
-import Test.Tasty.QuickCheck hiding (scale)
+import Test.Tasty.QuickCheck hiding (scale, numTests)
 import Test.QuickCheck.Classes
 
 import Quaternion
@@ -54,11 +55,24 @@
 testSuite = testGroup "Sparse"
     [ arithmeticTests
     , otherTests
-    , semiringTests
+    , lawsTests
     , evalTests
     , derivTests
+#if MIN_VERSION_semirings(0,4,2)
+    , gcdExtTests
+#endif
     ]
 
+lawsTests :: TestTree
+lawsTests = testGroup "Laws"
+  [ semiringTests
+  , ringTests
+  , numTests
+  , euclideanTests
+  , isListTests
+  , showTests
+  ]
+
 semiringTests :: TestTree
 semiringTests
   = testGroup "Semiring"
@@ -68,14 +82,72 @@
   , semiringLaws (Proxy :: Proxy (Poly U.Vector Int8))
   , semiringLaws (Proxy :: Proxy (Poly V.Vector Integer))
   , semiringLaws (Proxy :: Proxy (Poly U.Vector (Quaternion Int)))
+  ]
+
+ringTests :: TestTree
+ringTests
+  = testGroup "Ring"
+  $ map (uncurry testProperty)
+  $ concatMap lawsProperties
+  [
 #if MIN_VERSION_quickcheck_classes(0,6,1)
-  , ringLaws (Proxy :: Proxy (Poly U.Vector ()))
+    ringLaws (Proxy :: Proxy (Poly U.Vector ()))
   , ringLaws (Proxy :: Proxy (Poly U.Vector Int8))
   , ringLaws (Proxy :: Proxy (Poly V.Vector Integer))
   , ringLaws (Proxy :: Proxy (Poly U.Vector (Quaternion Int)))
 #endif
   ]
 
+numTests :: TestTree
+numTests
+  = testGroup "Num"
+  $ map (uncurry testProperty)
+  $ concatMap lawsProperties
+  [
+#if MIN_VERSION_quickcheck_classes(0,6,3)
+    numLaws (Proxy :: Proxy (Poly U.Vector Int8))
+  , numLaws (Proxy :: Proxy (Poly V.Vector Integer))
+  , numLaws (Proxy :: Proxy (Poly U.Vector (Quaternion Int)))
+#endif
+  ]
+
+euclideanTests :: TestTree
+euclideanTests
+  = testGroup "Euclidean"
+  $ map (uncurry testProperty)
+  $ concatMap lawsProperties
+  [
+#if MIN_VERSION_semirings(0,4,2) && MIN_VERSION_quickcheck_classes(0,6,3)
+    gcdDomainLaws (Proxy :: Proxy (ShortPoly (Poly V.Vector Integer)))
+  , euclideanLaws (Proxy :: Proxy (ShortPoly (Poly V.Vector Rational)))
+#endif
+  ]
+
+isListTests :: TestTree
+isListTests
+  = testGroup "IsList"
+  $ map (uncurry testProperty)
+  $ concatMap lawsProperties
+  [ isListLaws (Proxy :: Proxy (Poly U.Vector ()))
+  , isListLaws (Proxy :: Proxy (Poly U.Vector Int8))
+  , isListLaws (Proxy :: Proxy (Poly V.Vector Integer))
+  , isListLaws (Proxy :: Proxy (Poly U.Vector (Quaternion Int)))
+  ]
+
+showTests :: TestTree
+showTests
+  = testGroup "Show"
+  $ map (uncurry testProperty)
+  $ concatMap lawsProperties
+  [
+#if MIN_VERSION_quickcheck_classes(0,6,0)
+    showLaws (Proxy :: Proxy (Poly U.Vector ()))
+  , showLaws (Proxy :: Proxy (Poly U.Vector Int8))
+  , showLaws (Proxy :: Proxy (Poly V.Vector Integer))
+  , showLaws (Proxy :: Proxy (Poly U.Vector (Quaternion Int)))
+#endif
+  ]
+
 arithmeticTests :: TestTree
 arithmeticTests = testGroup "Arithmetic"
   [ testProperty "addition matches reference" $
@@ -180,6 +252,8 @@
 derivTests = testGroup "deriv"
   [ testProperty "deriv = S.deriv" $
     \(p :: Poly V.Vector Integer) -> deriv p === S.deriv p
+  , testProperty "integral = S.integral" $
+    \(p :: Poly V.Vector Rational) -> integral p === S.integral p
   , testProperty "deriv . integral = id" $
     \(p :: Poly V.Vector Rational) -> deriv (integral p) === p
   , testProperty "deriv c = 0" $
@@ -196,3 +270,22 @@
   --     deriv (eval (toPoly $ fmap (fmap $ monomial 0) $ unPoly p) q) ===
   --       deriv q * eval (toPoly $ fmap (fmap $ monomial 0) $ unPoly $ deriv p) q
   ]
+
+#if MIN_VERSION_semirings(0,4,2)
+gcdExtTests :: TestTree
+gcdExtTests = localOption (QuickCheckMaxSize 12) $ testGroup "gcdExt"
+  [ testProperty "gcdExt == S.gcdExt" $
+    \(a :: Poly V.Vector Rational) b ->
+      gcdExt a b === S.gcdExt a b
+  , testProperty "g == as (mod b) for gcdExt" $
+    \(a :: Poly V.Vector Rational) b -> b /= 0 ==>
+      uncurry ((. flip rem b) . (===) . flip rem b) ((* a) <$> gcdExt a b)
+  , testProperty "fst . gcdExt == gcd (mod units)" $
+    \(a :: Poly V.Vector Rational) b ->
+      fst (gcdExt a b) `sameUpToUnits` gcd a b
+  ]
+
+sameUpToUnits :: (Eq a, GcdDomain a) => a -> a -> Bool
+sameUpToUnits x y = x == y ||
+  isJust (x `divide` y) && isJust (y `divide` x)
+#endif