diff --git a/README.md b/README.md
--- a/README.md
+++ b/README.md
@@ -32,9 +32,9 @@
 1 * X^2 + (-3) * X + 2
 ```
 
-(Unfortunately, a type is often ambiguous and must be given explicitly.)
+(Unfortunately, types are often ambiguous and must be given explicitly.)
 
-While being convenient to read and write in REPL, `X` is relatively slow. The fastest approach is to use `toPoly`, providing it with a vector of coefficients (head is the constant term):
+While being convenient to read and write in REPL, `X` is relatively slow. The fastest approach is to use `toPoly`, providing it with a vector of coefficients (constant term first):
 
 ```haskell
 > toPoly (Data.Vector.Unboxed.fromList [2, -3, 1 :: Int])
@@ -99,7 +99,7 @@
 
 ## Deconstruction
 
-Use `unPoly` to deconstruct a polynomial to a vector of coefficients (head is the constant term):
+Use `unPoly` to deconstruct a polynomial to a vector of coefficients (constant term first):
 
 ```haskell
 > unPoly (X^2 - 3 * X + 2 :: UPoly Int)
diff --git a/bench/DenseBench.hs b/bench/DenseBench.hs
--- a/bench/DenseBench.hs
+++ b/bench/DenseBench.hs
@@ -1,4 +1,3 @@
-{-# LANGUAGE CPP                        #-}
 {-# LANGUAGE GeneralizedNewtypeDeriving #-}
 {-# LANGUAGE RankNTypes                 #-}
 {-# LANGUAGE TypeApplications           #-}
@@ -9,15 +8,13 @@
 
 import Prelude hiding (quotRem, gcd)
 import Gauge.Main
-import Data.Poly
-import qualified Data.Vector.Unboxed as U
-#if MIN_VERSION_semirings(0,5,2)
 import Data.Euclidean (Euclidean(..), GcdDomain(..), Field)
+import Data.Poly
 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
+import qualified Data.Vector.Unboxed as U
 
 benchSuite :: Benchmark
 benchSuite = bgroup "dense" $ concat
@@ -26,14 +23,10 @@
   , map benchEval     [100, 1000, 10000]
   , map benchDeriv    [100, 1000, 10000]
   , map benchIntegral [100, 1000, 10000]
-#if MIN_VERSION_semirings(0,5,2)
   , 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
   ]
 
 benchAdd :: Int -> Benchmark
@@ -51,28 +44,18 @@
 benchIntegral :: Int -> Benchmark
 benchIntegral k = bench ("integral/" ++ show k) $ nf doIntegral k
 
-#if MIN_VERSION_semirings(0,5,2)
-
 benchQuotRem :: Int -> Benchmark
 benchQuotRem k = bench ("quotRem/" ++ show k) $ nf doQuotRem k
 
 benchGcd :: Int -> Benchmark
 benchGcd k = bench ("gcd/" ++ show k) $ nf doGcd 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 (getMod2 . doGcdExt @Mod2) k
-
 benchGcdFracM :: Int -> Benchmark
 benchGcdFracM k = bench ("gcdFrac/Mod2/" ++ show k) $ nf (getMod2 . doGcdFrac @Mod2) k
 
-#endif
-
 doBinOp :: (forall a. Num a => a -> a -> a) -> Int -> Int
 doBinOp op n = U.sum zs
   where
@@ -98,8 +81,6 @@
     xs = toPoly $ U.generate n ((* 2) . fromIntegral)
     zs = unPoly $ integral xs
 
-#if MIN_VERSION_semirings(0,5,2)
-
 gen1 :: Ring a => Int -> a
 gen1 k = S.fromIntegral (truncate (pi * fromIntegral k :: Double) `mod` (k + 1))
 
@@ -120,18 +101,9 @@
     ys = toPoly $ V.generate n gen2
     gs = unPoly $ xs `gcd` ys
 
-doGcdExt :: (Eq a, Field a) => Int -> a
-doGcdExt n = V.foldl' plus zero gs
-  where
-    xs = S.toPoly $ V.generate n gen1
-    ys = S.toPoly $ V.generate n gen2
-    gs = unPoly $ fst $ xs `gcdExt` ys
-
 doGcdFrac :: (Eq a, Field a) => Int -> a
 doGcdFrac n = V.foldl' plus zero gs
   where
     xs = PolyOverField $ S.toPoly $ V.generate n gen1
     ys = PolyOverField $ S.toPoly $ V.generate n gen2
     gs = unPoly $ unPolyOverField $ xs `gcd` ys
-
-#endif
diff --git a/changelog.md b/changelog.md
--- a/changelog.md
+++ b/changelog.md
@@ -1,3 +1,10 @@
+# 0.4.0.0
+
+* Implement Laurent polynomials.
+* Implement orthogonal polynomials.
+* Decomission extended GCD, use `Data.Euclidean.gcdExt`.
+* Decomission `PolyOverFractional`, use `PolyOverField`.
+
 # 0.3.3.0
 
 * Add function `subst`.
diff --git a/poly.cabal b/poly.cabal
--- a/poly.cabal
+++ b/poly.cabal
@@ -1,5 +1,5 @@
 name: poly
-version: 0.3.3.0
+version: 0.4.0.0
 synopsis: Polynomials
 description:
   Polynomials backed by `Vector`.
@@ -8,14 +8,14 @@
 license-file: LICENSE
 author: Andrew Lelechenko
 maintainer: andrew.lelechenko@gmail.com
-copyright: 2019 Andrew Lelechenko
+copyright: 2019-2020 Andrew Lelechenko
 category: Math, Numerical
 build-type: Simple
-extra-source-files: README.md
 cabal-version: >=1.10
-tested-with: GHC ==8.0.2 GHC ==8.2.2 GHC ==8.4.4 GHC ==8.6.5 GHC ==8.8.1
+tested-with: GHC ==8.0.2 GHC ==8.2.2 GHC ==8.4.4 GHC ==8.6.5 GHC ==8.8.3 GHC ==8.10.1
 extra-source-files:
   changelog.md
+  README.md
 
 source-repository head
   type: git
@@ -25,8 +25,11 @@
   hs-source-dirs: src
   exposed-modules:
     Data.Poly
+    Data.Poly.Laurent
+    Data.Poly.Orthogonal
     Data.Poly.Semiring
     Data.Poly.Sparse
+    Data.Poly.Sparse.Laurent
     Data.Poly.Sparse.Semiring
   other-modules:
     Data.Poly.Internal.Dense
@@ -40,31 +43,36 @@
     base >= 4.9 && < 5,
     deepseq >= 1.1 && < 1.5,
     primitive >= 0.6,
-    semirings >= 0.2,
+    semirings >= 0.5.2,
     vector >= 0.12.0.2,
-    vector-algorithms >= 0.7
+    vector-algorithms >= 0.8.0.3
   default-language: Haskell2010
-  ghc-options: -Wall
+  ghc-options: -Wall -Wcompat
 
 test-suite poly-tests
   type: exitcode-stdio-1.0
   main-is: Main.hs
   other-modules:
     Dense
+    DenseLaurent
+    Orthogonal
     Quaternion
     Sparse
+    SparseLaurent
+    TestUtils
   build-depends:
     base >=4.9 && <5,
+    mod,
     poly,
     QuickCheck >=2.12,
-    quickcheck-classes >=0.5,
-    semirings >= 0.2,
+    quickcheck-classes >=0.6.3,
+    semirings >= 0.5.2,
     tasty >= 0.11,
     tasty-quickcheck >= 0.8,
     vector >= 0.12.0.2
   default-language: Haskell2010
   hs-source-dirs: test
-  ghc-options: -Wall
+  ghc-options: -Wall -Wcompat -threaded -rtsopts
 
 benchmark poly-gauge
   build-depends:
@@ -81,4 +89,4 @@
     SparseBench
   default-language: Haskell2010
   hs-source-dirs: bench
-  ghc-options: -Wall
+  ghc-options: -Wall -Wcompat
diff --git a/src/Data/Poly.hs b/src/Data/Poly.hs
--- a/src/Data/Poly.hs
+++ b/src/Data/Poly.hs
@@ -7,7 +7,6 @@
 -- Dense polynomials and a 'Num'-based interface.
 --
 
-{-# LANGUAGE CPP             #-}
 {-# LANGUAGE PatternSynonyms #-}
 
 module Data.Poly
@@ -16,7 +15,6 @@
   , UPoly
   , unPoly
   , leading
-  -- * Num interface
   , toPoly
   , monomial
   , scale
@@ -25,19 +23,10 @@
   , subst
   , deriv
   , integral
-#if MIN_VERSION_semirings(0,4,2)
-  -- * 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.Field (gcdExt)
+import Data.Poly.Internal.Dense.Field ()
 import Data.Poly.Internal.Dense.GcdDomain ()
 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
@@ -7,7 +7,6 @@
 -- Dense polynomials of one variable.
 --
 
-{-# LANGUAGE CPP                        #-}
 {-# LANGUAGE FlexibleInstances          #-}
 {-# LANGUAGE GeneralizedNewtypeDeriving #-}
 {-# LANGUAGE PatternSynonyms            #-}
@@ -37,11 +36,10 @@
   , pattern X'
   , eval'
   , subst'
+  , substitute'
   , deriv'
-#if MIN_VERSION_semirings(0,5,0)
   , unscale'
   , integral'
-#endif
   ) where
 
 import Prelude hiding (quotRem, quot, rem, gcd, lcm, (^))
@@ -50,6 +48,7 @@
 import Control.Monad.Primitive
 import Control.Monad.ST
 import Data.Bits
+import Data.Euclidean (Euclidean, Field, quot)
 import Data.List (foldl', intersperse)
 import Data.Semiring (Semiring(..), Ring())
 import qualified Data.Semiring as Semiring
@@ -58,13 +57,6 @@
 import qualified Data.Vector.Generic.Mutable as MG
 import qualified Data.Vector.Unboxed as U
 import GHC.Exts
-#if !MIN_VERSION_semirings(0,4,0)
-import Data.Semigroup
-import Numeric.Natural
-#endif
-#if MIN_VERSION_semirings(0,5,0)
-import Data.Euclidean (Euclidean, Field, quot)
-#endif
 
 -- | Polynomials of one variable with coefficients from @a@,
 -- backed by a 'G.Vector' @v@ (boxed, unboxed, storable, etc.).
@@ -171,13 +163,11 @@
   {-# INLINE plus #-}
   {-# INLINE times #-}
 
-#if MIN_VERSION_semirings(0,4,0)
   fromNatural n = if n' == zero then zero else Poly $ G.singleton n'
     where
       n' :: a
       n' = fromNatural n
   {-# INLINE fromNatural #-}
-#endif
 
 instance (Eq a, Ring a, G.Vector v a) => Ring (Poly v a) where
   negate (Poly xs) = Poly $ G.map Semiring.negate xs
@@ -365,7 +355,6 @@
 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
@@ -374,7 +363,6 @@
     MG.unsafeWrite zs k (G.unsafeIndex xs (k + fromIntegral yp) `quot` yc)
   G.unsafeFreeze zs
 {-# INLINABLE unscale' #-}
-#endif
 
 data StrictPair a b = !a :*: !b
 
@@ -433,17 +421,6 @@
   | otherwise = toPoly' $ G.imap (\i x -> fromNatural (fromIntegral (i + 1)) `times` x) $ G.tail xs
 {-# INLINE deriv' #-}
 
-#if !MIN_VERSION_semirings(0,4,0)
-fromNatural :: Semiring a => Natural -> a
-fromNatural 0 = zero
-fromNatural n = getAdd' (stimes n (Add' one))
-
-newtype Add' a = Add' { getAdd' :: a }
-
-instance Semiring a => Semigroup (Add' a) where
-  Add' a <> Add' b = Add' (a `plus` b)
-#endif
-
 -- | Compute an indefinite integral of a polynomial,
 -- setting constant term to zero.
 --
@@ -462,7 +439,6 @@
       lenXs = G.length xs
 {-# INLINABLE 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
@@ -475,7 +451,6 @@
     where
       lenXs = G.length xs
 {-# INLINABLE 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
--- a/src/Data/Poly/Internal/Dense/Field.hs
+++ b/src/Data/Poly/Internal/Dense/Field.hs
@@ -8,7 +8,6 @@
 --
 
 {-# LANGUAGE ConstraintKinds            #-}
-{-# LANGUAGE CPP                        #-}
 {-# LANGUAGE FlexibleInstances          #-}
 {-# LANGUAGE PatternSynonyms            #-}
 {-# LANGUAGE ScopedTypeVariables        #-}
@@ -16,24 +15,17 @@
 
 {-# 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 Prelude hiding (quotRem, quot, rem, gcd, recip)
 import Control.Exception
 import Control.Monad
 import Control.Monad.Primitive
 import Control.Monad.ST
-import Data.Euclidean (Euclidean(..))
-#if !MIN_VERSION_semirings(0,5,0)
-import Data.Semiring (Ring)
-#else
-import Data.Euclidean (Field)
-#endif
+import Data.Euclidean (Euclidean(..), Field)
+import Data.Field (recip)
 import Data.Semiring (times, minus, zero, one)
 import qualified Data.Vector.Generic as G
 import qualified Data.Vector.Generic.Mutable as MG
@@ -41,10 +33,6 @@
 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)
 
@@ -57,32 +45,40 @@
   {-# INLINE rem #-}
 
 quotientAndRemainder
-  :: (Field a, G.Vector v a)
+  :: (Eq a, 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)
+  | lenXs < lenYs = (G.empty, xs)
+  | lenYs == 0 = throw DivideByZero
+  | lenYs == 1 = let invY = recip (G.unsafeHead ys) in
+                 (G.map (`times` invY) xs, G.empty)
   | 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
+    let yLast = G.unsafeLast ys
+        invYLast = recip yLast
     forM_ [lenQs - 1, lenQs - 2 .. 0] $ \i -> do
       r <- MG.unsafeRead rs (lenYs - 1 + i)
-      let q = r `quot` G.unsafeLast ys
+      let q = if yLast == one then r else r `times` invYLast
       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)
+      MG.unsafeWrite rs (lenYs - 1 + i) zero
+      forM_ [0 .. lenYs - 2] $ \k -> do
+        let y = G.unsafeIndex ys k
+        when (y /= zero) $
+          MG.unsafeModify rs (\c -> c `minus` q `times` y) (i + k)
     let rs' = MG.unsafeSlice 0 lenYs rs
     (,) <$> G.unsafeFreeze qs <*> G.unsafeFreeze rs'
+  where
+    lenXs = G.length xs
+    lenYs = G.length ys
+    lenQs = lenXs - lenYs + 1
 {-# INLINABLE quotientAndRemainder #-}
 
 remainder
-  :: (Field a, G.Vector v a)
+  :: (Eq a, Field a, G.Vector v a)
   => v a
   -> v a
   -> v a
@@ -96,25 +92,29 @@
 {-# INLINABLE remainder #-}
 
 remainderM
-  :: (PrimMonad m, Field a, G.Vector v a)
+  :: (PrimMonad m, Eq a, 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 ()
+  | lenXs < lenYs = pure ()
+  | lenYs == 0 = throw DivideByZero
+  | lenYs == 1 = MG.set xs zero
   | otherwise = do
-    let lenXs = MG.length xs
-        lenYs = MG.length ys
-        lenQs = lenXs - lenYs + 1
     yLast <- MG.unsafeRead ys (lenYs - 1)
+    let invYLast = recip yLast
     forM_ [lenQs - 1, lenQs - 2 .. 0] $ \i -> do
       r <- MG.unsafeRead xs (lenYs - 1 + i)
-      forM_ [0 .. lenYs - 1] $ \k -> do
+      MG.unsafeWrite xs (lenYs - 1 + i) zero
+      let q = if yLast == one then r else r `times` invYLast
+      forM_ [0 .. lenYs - 2] $ \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)
+        when (y /= zero) $
+          MG.unsafeModify xs (\c -> c `minus` q `times` y) (i + k)
+  where
+    lenXs = MG.length xs
+    lenYs = MG.length ys
+    lenQs = lenXs - lenYs + 1
 {-# INLINABLE remainderM #-}
 
 fieldGcd
@@ -139,53 +139,3 @@
     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'
-{-# INLINABLE 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/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
@@ -7,7 +7,6 @@
 -- GcdDomain for GcdDomain underlying
 --
 
-{-# LANGUAGE CPP                        #-}
 {-# LANGUAGE FlexibleInstances          #-}
 {-# LANGUAGE PatternSynonyms            #-}
 {-# LANGUAGE ScopedTypeVariables        #-}
@@ -18,8 +17,6 @@
 module Data.Poly.Internal.Dense.GcdDomain
   () where
 
-#if MIN_VERSION_semirings(0,4,2)
-
 import Prelude hiding (gcd, lcm, (^))
 import Control.Exception
 import Control.Monad
@@ -170,5 +167,3 @@
 
     go (lenQs - 1)
 {-# INLINABLE quotient #-}
-
-#endif
diff --git a/src/Data/Poly/Internal/PolyOverField.hs b/src/Data/Poly/Internal/PolyOverField.hs
--- a/src/Data/Poly/Internal/PolyOverField.hs
+++ b/src/Data/Poly/Internal/PolyOverField.hs
@@ -7,19 +7,13 @@
 -- 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, (^))
@@ -36,23 +30,6 @@
 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 #-}
@@ -67,9 +44,3 @@
   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/Sparse.hs b/src/Data/Poly/Internal/Sparse.hs
--- a/src/Data/Poly/Internal/Sparse.hs
+++ b/src/Data/Poly/Internal/Sparse.hs
@@ -7,7 +7,6 @@
 -- Sparse polynomials of one variable.
 --
 
-{-# LANGUAGE CPP                        #-}
 {-# LANGUAGE FlexibleContexts           #-}
 {-# LANGUAGE GeneralizedNewtypeDeriving #-}
 {-# LANGUAGE PatternSynonyms            #-}
@@ -38,10 +37,9 @@
   , pattern X'
   , eval'
   , subst'
+  , substitute'
   , deriv'
-#if MIN_VERSION_semirings(0,5,0)
   , integral'
-#endif
   ) where
 
 import Prelude hiding (quot)
@@ -50,6 +48,7 @@
 import Control.Monad.Primitive
 import Control.Monad.ST
 import Data.Bits
+import Data.Euclidean (Field, quot)
 import Data.List (intersperse)
 import Data.Ord
 import Data.Semiring (Semiring(..), Ring())
@@ -60,13 +59,6 @@
 import qualified Data.Vector.Unboxed as U
 import qualified Data.Vector.Algorithms.Tim as Tim
 import GHC.Exts
-#if !MIN_VERSION_semirings(0,4,0)
-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.).
@@ -215,13 +207,11 @@
   {-# INLINE plus #-}
   {-# INLINE times #-}
 
-#if MIN_VERSION_semirings(0,4,0)
   fromNatural n = if n' == zero then zero else Poly $ G.singleton (0, n')
     where
       n' :: a
       n' = fromNatural n
   {-# INLINE fromNatural #-}
-#endif
 
 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
@@ -528,17 +518,6 @@
   xs
 {-# INLINE deriv' #-}
 
-#if !MIN_VERSION_semirings(0,4,0)
-fromNatural :: Semiring a => Natural -> a
-fromNatural 0 = zero
-fromNatural n = getAdd' (stimes n (Add' one))
-
-newtype Add' a = Add' { getAdd' :: a }
-
-instance Semiring a => Semigroup (Add' a) where
-  Add' a <> Add' b = Add' (a `plus` b)
-#endif
-
 derivPoly
   :: G.Vector v (Word, a)
   => (a -> Bool)
@@ -575,13 +554,11 @@
   $ 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
--- a/src/Data/Poly/Internal/Sparse/Field.hs
+++ b/src/Data/Poly/Internal/Sparse/Field.hs
@@ -8,7 +8,6 @@
 --
 
 {-# LANGUAGE ConstraintKinds            #-}
-{-# LANGUAGE CPP                        #-}
 {-# LANGUAGE FlexibleContexts           #-}
 {-# LANGUAGE FlexibleInstances          #-}
 {-# LANGUAGE PatternSynonyms            #-}
@@ -18,31 +17,18 @@
 
 {-# OPTIONS_GHC -fno-warn-orphans #-}
 
-#if MIN_VERSION_semirings(0,4,2)
-
-module Data.Poly.Internal.Sparse.Field
-  ( gcdExt
-  ) where
+module Data.Poly.Internal.Sparse.Field () where
 
 import Prelude hiding (quotRem, quot, rem, gcd)
 import Control.Arrow
 import Control.Exception
-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 Data.Euclidean (Euclidean(..), Field)
+import Data.Semiring (minus, plus, times, zero)
 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
@@ -68,53 +54,3 @@
           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'
-{-# INLINABLE 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/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
@@ -7,7 +7,6 @@
 -- GcdDomain for GcdDomain underlying
 --
 
-{-# LANGUAGE CPP                        #-}
 {-# LANGUAGE FlexibleContexts           #-}
 {-# LANGUAGE FlexibleInstances          #-}
 {-# LANGUAGE PatternSynonyms            #-}
@@ -20,8 +19,6 @@
 module Data.Poly.Internal.Sparse.GcdDomain
   () where
 
-#if MIN_VERSION_semirings(0,4,2)
-
 import Prelude hiding (gcd, lcm, (^))
 import Control.Exception
 import Data.Euclidean
@@ -75,5 +72,3 @@
         gx = fromMaybe err $ divide g xc
         gy = fromMaybe err $ divide g yc
         err = error "gcd: violated internal invariant"
-
-#endif
diff --git a/src/Data/Poly/Laurent.hs b/src/Data/Poly/Laurent.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Poly/Laurent.hs
@@ -0,0 +1,284 @@
+-- |
+-- Module:      Data.Poly.Laurent
+-- Copyright:   (c) 2020 Andrew Lelechenko
+-- Licence:     BSD3
+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>
+--
+-- <https://en.wikipedia.org/wiki/Laurent_polynomial Laurent polynomials>.
+--
+
+{-# LANGUAGE FlexibleInstances          #-}
+{-# LANGUAGE GeneralizedNewtypeDeriving #-}
+{-# LANGUAGE PatternSynonyms            #-}
+{-# LANGUAGE ScopedTypeVariables        #-}
+{-# LANGUAGE ViewPatterns               #-}
+
+module Data.Poly.Laurent
+  ( Laurent
+  , VLaurent
+  , ULaurent
+  , unLaurent
+  , toLaurent
+  , leading
+  , monomial
+  , scale
+  , pattern X
+  , (^-)
+  , eval
+  , subst
+  , deriv
+  , LaurentOverField(..)
+  ) where
+
+import Prelude hiding (quotRem, quot, rem, gcd)
+import Control.Arrow (first)
+import Control.DeepSeq (NFData(..))
+import Data.Euclidean (GcdDomain(..), Euclidean(..), Field)
+import Data.List (intersperse)
+import Data.Semiring (Semiring(..), Ring())
+import qualified Data.Semiring as Semiring
+import qualified Data.Vector as V
+import qualified Data.Vector.Generic as G
+import qualified Data.Vector.Unboxed as U
+
+import Data.Poly.Internal.Dense (Poly(..))
+import qualified Data.Poly.Internal.Dense as Dense
+import Data.Poly.Internal.Dense.Field ()
+import Data.Poly.Internal.Dense.GcdDomain ()
+import Data.Poly.Internal.PolyOverField
+
+-- | <https://en.wikipedia.org/wiki/Laurent_polynomial Laurent polynomials>
+-- of one variable with coefficients from @a@,
+-- backed by a 'G.Vector' @v@ (boxed, unboxed, storable, etc.).
+--
+-- Use pattern 'X' and operator '^-' for construction:
+--
+-- >>> (X + 1) + (X^-1 - 1) :: VLaurent Integer
+-- 1 * X + 0 + 1 * X^-1
+-- >>> (X + 1) * (1 - X^-1) :: ULaurent Int
+-- 1 * X + 0 + (-1) * X^-1
+--
+-- Polynomials are stored normalized, without leading
+-- and trailing
+-- zero coefficients, so 0 * X + 1 + 0 * X^-1 equals to 1.
+--
+-- 'Ord' instance does not make much sense mathematically,
+-- it is defined only for the sake of 'Data.Set.Set', 'Data.Map.Map', etc.
+--
+data Laurent v a = Laurent !Int !(Poly v a)
+  deriving (Eq, Ord)
+
+-- | Deconstruct a 'Laurent' polynomial into an offset (largest possible)
+-- and a regular polynomial.
+--
+-- >>> unLaurent (2 * X + 1 :: ULaurent Int)
+-- (0,2 * X + 1)
+-- >>> unLaurent (1 + 2 * X^-1 :: ULaurent Int)
+-- (-1,1 * X + 2)
+-- >>> unLaurent (2 * X^2 + X :: ULaurent Int)
+-- (1,2 * X + 1)
+-- >>> unLaurent (0 :: ULaurent Int)
+-- (0,0)
+unLaurent :: Laurent v a -> (Int, Poly v a)
+unLaurent (Laurent off poly) = (off, poly)
+
+-- | Construct 'Laurent' polynomial from an offset and a regular polynomial.
+-- One can imagine it as 'Data.Poly.scale'', but allowing negative offsets.
+--
+-- >>> toLaurent 2 (2 * Data.Poly.X + 1) :: ULaurent Int
+-- 2 * X^3 + 1 * X^2
+-- >>> toLaurent (-2) (2 * Data.Poly.X + 1) :: ULaurent Int
+-- 2 * X^-1 + 1 * X^-2
+toLaurent
+  :: (Eq a, Semiring a, G.Vector v a)
+  => Int
+  -> Poly v a
+  -> Laurent v a
+toLaurent off (Poly xs) = go 0
+  where
+    go k
+      | k >= G.length xs
+      = Laurent 0 zero
+      | G.unsafeIndex xs k == zero
+      = go (k + 1)
+      | otherwise
+      = Laurent (off + k) (Poly (G.unsafeDrop k xs))
+{-# INLINE toLaurent #-}
+
+toLaurentNum
+  :: (Eq a, Num a, G.Vector v a)
+  => Int
+  -> Poly v a
+  -> Laurent v a
+toLaurentNum off (Poly xs) = go 0
+  where
+    go k
+      | k >= G.length xs
+      = Laurent 0 0
+      | G.unsafeIndex xs k == 0
+      = go (k + 1)
+      | otherwise
+      = Laurent (off + k) (Poly (G.unsafeDrop k xs))
+{-# INLINE toLaurentNum #-}
+
+instance NFData (v a) => NFData (Laurent v a) where
+  rnf (Laurent off poly) = rnf off `seq` rnf poly
+
+instance (Show a, G.Vector v a) => Show (Laurent v a) where
+  showsPrec d (Laurent off poly)
+    | G.null (unPoly poly)
+      = showString "0"
+    | otherwise
+      = showParen (d > 0)
+      $ foldl (.) id
+      $ intersperse (showString " + ")
+      $ G.ifoldl (\acc i c -> showCoeff (i + off) c : acc) []
+      $ unPoly poly
+    where
+      showCoeff 0 c = showsPrec 7 c
+      showCoeff 1 c = showsPrec 7 c . showString " * X"
+      showCoeff i c = showsPrec 7 c . showString (" * X^" ++ show i)
+
+-- | Laurent polynomials backed by boxed vectors.
+type VLaurent = Laurent V.Vector
+
+-- | Laurent polynomials backed by unboxed vectors.
+type ULaurent = Laurent U.Vector
+
+-- | Return a leading power and coefficient of a non-zero polynomial.
+--
+-- >>> leading ((2 * X + 1) * (2 * X^2 - 1) :: ULaurent Int)
+-- Just (3,4)
+-- >>> leading (0 :: ULaurent Int)
+-- Nothing
+leading :: G.Vector v a => Laurent v a -> Maybe (Int, a)
+leading (Laurent off poly) = first ((+ off) . fromIntegral) <$> Dense.leading poly
+
+-- | Note that 'abs' = 'id' and 'signum' = 'const' 1.
+instance (Eq a, Num a, G.Vector v a) => Num (Laurent v a) where
+  Laurent off1 poly1 * Laurent off2 poly2 = toLaurentNum (off1 + off2) (poly1 * poly2)
+  Laurent off1 poly1 + Laurent off2 poly2 = case off1 `compare` off2 of
+    LT -> toLaurentNum off1 (poly1 + Dense.scale (fromIntegral $ off2 - off1) 1 poly2)
+    EQ -> toLaurentNum off1 (poly1 + poly2)
+    GT -> toLaurentNum off2 (Dense.scale (fromIntegral $ off1 - off2) 1 poly1 + poly2)
+  Laurent off1 poly1 - Laurent off2 poly2 = case off1 `compare` off2 of
+    LT -> toLaurentNum off1 (poly1 - Dense.scale (fromIntegral $ off2 - off1) 1 poly2)
+    EQ -> toLaurentNum off1 (poly1 - poly2)
+    GT -> toLaurentNum off2 (Dense.scale (fromIntegral $ off1 - off2) 1 poly1 - poly2)
+  negate (Laurent off poly) = Laurent off (negate poly)
+  abs = id
+  signum = const 1
+  fromInteger n = Laurent 0 (fromInteger n)
+  {-# INLINE (+) #-}
+  {-# INLINE (-) #-}
+  {-# INLINE negate #-}
+  {-# INLINE fromInteger #-}
+  {-# INLINE (*) #-}
+
+instance (Eq a, Semiring a, G.Vector v a) => Semiring (Laurent v a) where
+  zero = Laurent 0 zero
+  one  = Laurent 0 one
+  Laurent off1 poly1 `times` Laurent off2 poly2 =
+    toLaurent (off1 + off2) (poly1 `times` poly2)
+  Laurent off1 poly1 `plus` Laurent off2 poly2 = case off1 `compare` off2 of
+    LT -> toLaurent off1 (poly1 `plus` Dense.scale' (fromIntegral $ off2 - off1) one poly2)
+    EQ -> toLaurent off1 (poly1 `plus` poly2)
+    GT -> toLaurent off2 (Dense.scale' (fromIntegral $ off1 - off2) one poly1 `plus` poly2)
+  fromNatural n = Laurent 0 (fromNatural n)
+  {-# INLINE zero #-}
+  {-# INLINE one #-}
+  {-# INLINE plus #-}
+  {-# INLINE times #-}
+  {-# INLINE fromNatural #-}
+
+instance (Eq a, Ring a, G.Vector v a) => Ring (Laurent v a) where
+  negate (Laurent off poly) = Laurent off (Semiring.negate poly)
+
+-- | Create a monomial from a power and a coefficient.
+monomial :: (Eq a, Semiring a, G.Vector v a) => Int -> a -> Laurent v a
+monomial p c
+  | c == zero = Laurent 0 zero
+  | otherwise = Laurent p (Dense.monomial' 0 c)
+{-# INLINE monomial #-}
+
+-- | Multiply a polynomial by a monomial, expressed as a power and a coefficient.
+--
+-- >>> scale 2 3 (X^2 + 1) :: ULaurent Int
+-- 3 * X^4 + 0 * X^3 + 3 * X^2 + 0 * X + 0
+scale :: (Eq a, Semiring a, G.Vector v a) => Int -> a -> Laurent v a -> Laurent v a
+scale yp yc (Laurent off poly) = toLaurent (off + yp) (Dense.scale' 0 yc poly)
+
+-- | Evaluate at a given point.
+--
+-- >>> eval (X^2 + 1 :: ULaurent Int) 3
+-- 10
+eval :: (Field a, G.Vector v a) => Laurent v a -> a -> a
+eval (Laurent off poly) x = Dense.eval' poly x `times`
+  (if off >= 0 then x Semiring.^ off else quot one x Semiring.^ (- off))
+{-# INLINE eval #-}
+
+-- | Substitute another polynomial instead of 'Data.Poly.X'.
+--
+-- >>> subst (X^2 + 1 :: UPoly Int) (X + 1 :: ULaurent Int)
+-- 1 * X^2 + 2 * X + 2
+subst :: (Eq a, Semiring a, G.Vector v a, G.Vector w a) => Poly v a -> Laurent w a -> Laurent w a
+subst = Dense.substitute' (scale 0)
+{-# INLINE subst #-}
+
+-- | Take a derivative.
+--
+-- >>> deriv (X^3 + 3 * X) :: ULaurent Int
+-- 3 * X^2 + 0 * X + 3
+deriv :: (Eq a, Ring a, G.Vector v a) => Laurent v a -> Laurent v a
+deriv (Laurent off (Poly xs)) =
+  toLaurent (off - 1) $ Dense.toPoly' $ G.imap (times . Semiring.fromIntegral . (+ off)) xs
+{-# INLINE deriv #-}
+
+-- | Create an identity polynomial.
+pattern X :: (Eq a, Semiring a, G.Vector v a, Eq (v a)) => Laurent v a
+pattern X <- ((==) var -> True)
+  where X = var
+
+var :: forall a v. (Eq a, Semiring a, G.Vector v a, Eq (v a)) => Laurent v a
+var
+  | (one :: a) == zero = Laurent 0 zero
+  | otherwise          = Laurent 1 one
+{-# INLINE var #-}
+
+-- | This operator can be applied only to 'X',
+-- but is instrumental to express Laurent polynomials in mathematical fashion:
+--
+-- >>> X + 2 + 3 * X^-1 :: ULaurent Int
+-- 1 * X + 2 + 3 * X^(-1)
+(^-)
+  :: (Eq a, Semiring a, G.Vector v a, Eq (v a))
+  => Laurent v a
+  -> Int
+  -> Laurent v a
+X^-n = monomial (negate n) one
+_^-_ = error "(^-) can be applied only to X"
+
+-- | Consider using 'LaurentOverField' wrapper,
+-- which provides a much faster implementation of
+-- 'Data.Euclidean.gcd' for polynomials over 'Field'.
+instance (Eq a, Ring a, GcdDomain a, Eq (v a), G.Vector v a) => GcdDomain (Laurent v a) where
+  divide (Laurent off1 poly1) (Laurent off2 poly2) =
+    toLaurent (off1 - off2) <$> divide poly1 poly2
+  {-# INLINE divide #-}
+
+  gcd (Laurent _ poly1) (Laurent _ poly2) =
+    toLaurent 0 (gcd poly1 poly2)
+  {-# INLINE gcd #-}
+
+-- | Wrapper for Laurent polynomials over 'Field',
+-- providing a faster 'GcdDomain' instance.
+newtype LaurentOverField laurent = LaurentOverField { unLaurentOverField :: laurent }
+  deriving (Eq, NFData, Num, Ord, Ring, Semiring, Show)
+
+instance (Eq a, Eq (v a), Field a, G.Vector v a) => GcdDomain (LaurentOverField (Laurent v a)) where
+  divide (LaurentOverField (Laurent off1 poly1)) (LaurentOverField (Laurent off2 poly2)) =
+    LaurentOverField . toLaurent (off1 - off2) . unPolyOverField <$> divide (PolyOverField poly1) (PolyOverField poly2)
+
+  gcd (LaurentOverField (Laurent _ poly1)) (LaurentOverField (Laurent _ poly2)) =
+    LaurentOverField (toLaurent 0 (unPolyOverField (gcd (PolyOverField poly1) (PolyOverField poly2))))
+  {-# INLINE gcd #-}
diff --git a/src/Data/Poly/Orthogonal.hs b/src/Data/Poly/Orthogonal.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Poly/Orthogonal.hs
@@ -0,0 +1,128 @@
+-- |
+-- Module:      Data.Poly.Orthogonal
+-- Copyright:   (c) 2019 Andrew Lelechenko
+-- Licence:     BSD3
+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>
+--
+-- Classical orthogonal polynomials.
+--
+
+{-# LANGUAGE OverloadedLists     #-}
+{-# LANGUAGE RebindableSyntax    #-}
+
+module Data.Poly.Orthogonal
+  ( legendre
+  , legendreShifted
+  , gegenbauer
+  , jacobi
+  , chebyshev1
+  , chebyshev2
+  , hermiteProb
+  , hermitePhys
+  , laguerre
+  , laguerreGen
+  ) where
+
+import Prelude hiding (quot, Num(..), fromIntegral)
+import Data.Euclidean
+import Data.Semiring
+import Data.Poly.Semiring
+import Data.Poly.Internal.Dense (unscale')
+import Data.Vector.Generic (Vector, fromListN)
+
+-- | <https://en.wikipedia.org/wiki/Legendre_polynomials Legendre polynomials>.
+--
+-- >>> take 3 legendre :: [Data.Poly.VPoly Double]
+-- [1.0,1.0 * X + 0.0,1.5 * X^2 + 0.0 * X + (-0.5)]
+legendre :: (Eq a, Field a, Vector v a) => [Poly v a]
+legendre = map (flip subst' (toPoly [1 `quot` 2, 1 `quot` 2])) legendreShifted
+  where
+    subst' :: (Eq a, Semiring a, Vector v a) => Poly v a -> Poly v a -> Poly v a
+    subst' = subst
+
+-- | <https://en.wikipedia.org/wiki/Legendre_polynomials#Shifted_Legendre_polynomials Shifted Legendre polynomials>.
+--
+-- >>> take 3 legendreShifted :: [Data.Poly.VPoly Integer]
+-- [1,2 * X + (-1),6 * X^2 + (-6) * X + 1]
+legendreShifted :: (Eq a, Euclidean a, Ring a, Vector v a) => [Poly v a]
+legendreShifted = xs
+  where
+    xs = 1 : toPoly [-1, 2] : zipWith3 rec (iterate (+ 1) 1) xs (tail xs)
+    rec n pm1 p = unscale' 0 (n + 1) (toPoly [-1 - 2 * n, 2 + 4 * n] * p - scale 0 n pm1)
+
+-- | <https://en.wikipedia.org/wiki/Gegenbauer_polynomials Gegenbauer polynomials>.
+gegenbauer :: (Eq a, Field a, Vector v a) => a -> [Poly v a]
+gegenbauer g = jacobi a a
+  where
+    a = g - 1 `quot` 2
+
+-- | <https://en.wikipedia.org/wiki/Jacobi_polynomials Jacobi polynomials>.
+jacobi :: (Eq a, Field a, Vector v a) => a -> a -> [Poly v a]
+jacobi a b = xs
+  where
+    x0 = 1
+    x1 = toPoly [(a - b) `quot` 2, (a + b + 2) `quot` 2]
+    xs = x0 : x1 : zipWith3 rec (iterate (+ 1) 2) xs (tail xs)
+    rec n pm1 p = toPoly [d, c] * p - scale 0 cm1 pm1
+      where
+        cp1 = 2 * n * (n + a + b) * (2 * n + a + b - 2)
+        q   = (2 * n + a + b - 1) `quot` cp1
+        c   = q * ((2 * n + a + b) * (2 * n + a + b - 2))
+        d   = q * (a * a - b * b)
+        cm1 = 2 * (n + a - 1) * (n + b - 1) * (2 * n + a + b) `quot` cp1
+
+-- | <https://en.wikipedia.org/wiki/Chebyshev_polynomials Chebyshev polynomials>
+-- of the first kind.
+--
+-- >>> take 3 chebyshev1 :: [VPoly Integer]
+-- [1,1 * X + 0,2 * X^2 + 0 * X + (-1)]
+chebyshev1 :: (Eq a, Ring a, Vector v a) => [Poly v a]
+chebyshev1 = xs
+  where
+    xs = 1 : monomial 1 1 : zipWith (\pm1 p -> scale 1 2 p - pm1) xs (tail xs)
+
+-- | <https://en.wikipedia.org/wiki/Chebyshev_polynomials Chebyshev polynomials>
+-- of the second kind.
+--
+-- >>> take 3 chebyshev2 :: [VPoly Integer]
+-- [1,2 * X + 0,4 * X^2 + 0 * X + (-1)]
+chebyshev2 :: (Eq a, Ring a, Vector v a) => [Poly v a]
+chebyshev2 = xs
+  where
+    xs = 1 : monomial 1 2 : zipWith (\pm1 p -> scale 1 2 p - pm1) xs (tail xs)
+
+-- | Probabilists' <https://en.wikipedia.org/wiki/Hermite_polynomials Hermite polynomials>.
+--
+-- >>> take 3 hermiteProb :: [VPoly Integer]
+-- [1,1 * X + 0,1 * X^2 + 0 * X + (-1)]
+hermiteProb :: (Eq a, Ring a, Vector v a) => [Poly v a]
+hermiteProb = xs
+  where
+    xs = 1 : monomial 1 1 : zipWith3 rec (iterate (+ 1) 1) xs (tail xs)
+    rec n pm1 p = scale 1 1 p - scale 0 n pm1
+
+-- | Physicists' <https://en.wikipedia.org/wiki/Hermite_polynomials Hermite polynomials>.
+--
+-- >>> take 3 hermitePhys :: [VPoly Double]
+-- [1,2 * X + 0,4 * X^2 + 0 * X + (-2)]
+hermitePhys :: (Eq a, Ring a, Vector v a) => [Poly v a]
+hermitePhys = xs
+  where
+    xs = 1 : monomial 1 2 : zipWith3 rec (iterate (+ 1) 1) xs (tail xs)
+    rec n pm1 p = scale 1 2 p - scale 0 (2 * n) pm1
+
+-- | <https://en.wikipedia.org/wiki/Laguerre_polynomials Laguerre polynomials>.
+--
+-- >>> take 3 laguerre :: [VPoly Double]
+-- [1.0,(-1.0) * X + 1.0,0.5 * X^2 + (-2.0) * X + 1.0]
+laguerre :: (Eq a, Field a, Vector v a) => [Poly v a]
+laguerre = laguerreGen 0
+
+-- | <https://en.wikipedia.org/wiki/Laguerre_polynomials#Generalized_Laguerre_polynomials Generalized Laguerre polynomials>
+laguerreGen :: (Eq a, Field a, Vector v a) => a -> [Poly v a]
+laguerreGen a = xs
+  where
+    x0 = 1
+    x1 = toPoly [1 + a, -1]
+    xs = x0 : x1 : zipWith3 rec (iterate (+ 1) 1) xs (tail xs)
+    rec n pm1 p = toPoly [(2 * n + 1 + a) `quot` (n + 1), -1 `quot` (n + 1)] * p - scale 0 ((n + a) `quot` (n + 1)) pm1
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,7 +7,6 @@
 -- Dense polynomials and a 'Semiring'-based interface.
 --
 
-{-# LANGUAGE CPP             #-}
 {-# LANGUAGE PatternSynonyms #-}
 
 module Data.Poly.Semiring
@@ -16,7 +15,6 @@
   , UPoly
   , unPoly
   , leading
-  -- * Semiring interface
   , toPoly
   , monomial
   , scale
@@ -24,32 +22,19 @@
   , eval
   , subst
   , deriv
-#if MIN_VERSION_semirings(0,5,0)
   , integral
-#endif
-#if MIN_VERSION_semirings(0,4,2)
-  -- * Polynomials over 'Field'
   , PolyOverField(..)
-  , gcdExt
-  , PolyOverFractional
-  , pattern PolyOverFractional
-  , unPolyOverFractional
-#endif
   ) where
 
+import Data.Euclidean (Field)
 import Data.Semiring (Semiring)
 import qualified Data.Vector.Generic as G
 
 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.Field (gcdExt)
+import Data.Poly.Internal.Dense.Field ()
 import Data.Poly.Internal.Dense.GcdDomain ()
 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).
@@ -98,7 +83,6 @@
 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.
 --
@@ -106,4 +90,3 @@
 -- 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,6 @@
 -- Sparse polynomials with 'Num' instance.
 --
 
-{-# LANGUAGE CPP             #-}
 {-# LANGUAGE PatternSynonyms #-}
 
 module Data.Poly.Sparse
@@ -16,7 +15,6 @@
   , UPoly
   , unPoly
   , leading
-  -- * Num interface
   , toPoly
   , monomial
   , scale
@@ -25,14 +23,8 @@
   , subst
   , deriv
   , integral
-#if MIN_VERSION_semirings(0,4,2)
-  -- * Polynomials over 'Field'
-  , gcdExt
-#endif
   ) where
 
 import Data.Poly.Internal.Sparse
-#if MIN_VERSION_semirings(0,4,2)
-import Data.Poly.Internal.Sparse.Field (gcdExt)
+import Data.Poly.Internal.Sparse.Field ()
 import Data.Poly.Internal.Sparse.GcdDomain ()
-#endif
diff --git a/src/Data/Poly/Sparse/Laurent.hs b/src/Data/Poly/Sparse/Laurent.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Poly/Sparse/Laurent.hs
@@ -0,0 +1,283 @@
+-- |
+-- Module:      Data.Poly.Sparse.Laurent
+-- Copyright:   (c) 2020 Andrew Lelechenko
+-- Licence:     BSD3
+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>
+--
+-- Sparse <https://en.wikipedia.org/wiki/Laurent_polynomial Laurent polynomials>.
+--
+
+{-# LANGUAGE FlexibleContexts           #-}
+{-# LANGUAGE FlexibleInstances          #-}
+{-# LANGUAGE GeneralizedNewtypeDeriving #-}
+{-# LANGUAGE PatternSynonyms            #-}
+{-# LANGUAGE ScopedTypeVariables        #-}
+{-# LANGUAGE StandaloneDeriving         #-}
+{-# LANGUAGE TypeFamilies               #-}
+{-# LANGUAGE UndecidableInstances       #-}
+{-# LANGUAGE ViewPatterns               #-}
+
+module Data.Poly.Sparse.Laurent
+  ( Laurent
+  , VLaurent
+  , ULaurent
+  , unLaurent
+  , toLaurent
+  , leading
+  , monomial
+  , scale
+  , pattern X
+  , (^-)
+  , eval
+  , subst
+  , deriv
+  ) where
+
+import Prelude hiding (quotRem, quot, rem, gcd)
+import Control.Arrow (first)
+import Control.DeepSeq (NFData(..))
+import Data.Euclidean (GcdDomain(..), Euclidean(..), Field)
+import Data.List (intersperse)
+import Data.Ord
+import Data.Semiring (Semiring(..), Ring())
+import qualified Data.Semiring as Semiring
+import qualified Data.Vector as V
+import qualified Data.Vector.Generic as G
+import qualified Data.Vector.Unboxed as U
+import GHC.Exts
+
+import Data.Poly.Internal.Sparse (Poly(..))
+import qualified Data.Poly.Internal.Sparse as Sparse
+import Data.Poly.Internal.Sparse.Field ()
+import Data.Poly.Internal.Sparse.GcdDomain ()
+
+-- | <https://en.wikipedia.org/wiki/Laurent_polynomial Laurent polynomials>
+-- of one variable with coefficients from @a@,
+-- backed by a 'G.Vector' @v@ (boxed, unboxed, storable, etc.).
+--
+-- Use pattern 'X' and operator '^-' for construction:
+--
+-- >>> (X + 1) + (X^-1 - 1) :: VLaurent Integer
+-- 1 * X + 1 * X^-1
+-- >>> (X + 1) * (1 - X^-1) :: ULaurent Int
+-- 1 * X + (-1) * X^-1
+--
+-- Polynomials are stored normalized, without
+-- zero coefficients, so 0 * X + 1 + 0 * X^-1 equals to 1.
+--
+-- 'Ord' instance does not make much sense mathematically,
+-- it is defined only for the sake of 'Data.Set.Set', 'Data.Map.Map', etc.
+--
+data Laurent v a = Laurent !Int !(Poly v a)
+
+deriving instance Eq  (v (Word, a)) => Eq  (Laurent v a)
+deriving instance Ord (v (Word, a)) => Ord (Laurent v a)
+
+instance (Eq a, Semiring a, G.Vector v (Word, a)) => IsList (Laurent v a) where
+  type Item (Laurent v a) = (Int, a)
+
+  fromList xs = toLaurent minPow (fromList ys)
+    where
+      minPow = minimum $ maxBound : map fst xs
+      ys = map (first (fromIntegral . (subtract minPow))) xs
+
+  fromListN n xs = toLaurent minPow (fromListN n ys)
+    where
+      minPow = minimum $ maxBound : map fst xs
+      ys = map (first (fromIntegral . (subtract minPow))) xs
+
+  toList (Laurent off poly) =
+    map (first ((+ off) . fromIntegral)) $ G.toList $ unPoly poly
+
+-- | Deconstruct a 'Laurent' polynomial into an offset (largest possible)
+-- and a regular polynomial.
+--
+-- >>> unLaurent (2 * X + 1 :: ULaurent Int)
+-- (0,2 * X + 1)
+-- >>> unLaurent (1 + 2 * X^-1 :: ULaurent Int)
+-- (-1,1 * X + 2)
+-- >>> unLaurent (2 * X^2 + X :: ULaurent Int)
+-- (1,2 * X + 1)
+-- >>> unLaurent (0 :: ULaurent Int)
+-- (0,0)
+unLaurent :: Laurent v a -> (Int, Poly v a)
+unLaurent (Laurent off poly) = (off, poly)
+
+-- | Construct 'Laurent' polynomial from an offset and a regular polynomial.
+-- One can imagine it as 'Data.Poly.Sparse.scale'', but allowing negative offsets.
+--
+-- >>> toLaurent 2 (2 * Data.Poly.Sparse.X + 1) :: ULaurent Int
+-- 2 * X^3 + 1 * X^2
+-- >>> toLaurent (-2) (2 * Data.Poly.Sparse.X + 1) :: ULaurent Int
+-- 2 * X^-1 + 1 * X^-2
+toLaurent
+  :: (Eq a, Semiring a, G.Vector v (Word, a))
+  => Int
+  -> Poly v a
+  -> Laurent v a
+toLaurent off (Poly xs)
+  | G.null xs = Laurent 0 zero
+  | otherwise = Laurent (off + fromIntegral minPow) (Poly ys)
+    where
+      minPow = fst $ G.minimumBy (comparing fst) xs
+      ys = if minPow == 0 then xs else G.map (first (subtract minPow)) xs
+{-# INLINE toLaurent #-}
+
+toLaurentNum
+  :: (Eq a, Num a, G.Vector v (Word, a))
+  => Int
+  -> Poly v a
+  -> Laurent v a
+toLaurentNum off (Poly xs)
+  | G.null xs = Laurent 0 0
+  | otherwise = Laurent (off + fromIntegral minPow) (Poly ys)
+    where
+      minPow = fst $ G.minimumBy (comparing fst) xs
+      ys = if minPow == 0 then xs else G.map (first (subtract minPow)) xs
+{-# INLINE toLaurentNum #-}
+
+instance NFData (v (Word, a)) => NFData (Laurent v a) where
+  rnf (Laurent off poly) = rnf off `seq` rnf poly
+
+instance (Show a, G.Vector v (Word, a)) => Show (Laurent v a) where
+  showsPrec d (Laurent off poly)
+    | G.null (unPoly poly)
+      = showString "0"
+    | otherwise
+      = showParen (d > 0)
+      $ foldl (.) id
+      $ intersperse (showString " + ")
+      $ G.ifoldl (\acc i c -> showCoeff (i + off) c : acc) []
+      $ unPoly poly
+    where
+      showCoeff 0 c = showsPrec 7 c
+      showCoeff 1 c = showsPrec 7 c . showString " * X"
+      showCoeff i c = showsPrec 7 c . showString (" * X^" ++ show i)
+
+-- | Laurent polynomials backed by boxed vectors.
+type VLaurent = Laurent V.Vector
+
+-- | Laurent polynomials backed by unboxed vectors.
+type ULaurent = Laurent U.Vector
+
+-- | Return a leading power and coefficient of a non-zero polynomial.
+--
+-- >>> leading ((2 * X + 1) * (2 * X^2 - 1) :: ULaurent Int)
+-- Just (3,4)
+-- >>> leading (0 :: ULaurent Int)
+-- Nothing
+leading :: G.Vector v (Word, a) => Laurent v a -> Maybe (Int, a)
+leading (Laurent off poly) = first ((+ off) . fromIntegral) <$> Sparse.leading poly
+
+-- | Note that 'abs' = 'id' and 'signum' = 'const' 1.
+instance (Eq a, Num a, G.Vector v (Word, a)) => Num (Laurent v a) where
+  Laurent off1 poly1 * Laurent off2 poly2 = toLaurentNum (off1 + off2) (poly1 * poly2)
+  Laurent off1 poly1 + Laurent off2 poly2 = case off1 `compare` off2 of
+    LT -> toLaurentNum off1 (poly1 + Sparse.scale (fromIntegral $ off2 - off1) 1 poly2)
+    EQ -> toLaurentNum off1 (poly1 + poly2)
+    GT -> toLaurentNum off2 (Sparse.scale (fromIntegral $ off1 - off2) 1 poly1 + poly2)
+  Laurent off1 poly1 - Laurent off2 poly2 = case off1 `compare` off2 of
+    LT -> toLaurentNum off1 (poly1 - Sparse.scale (fromIntegral $ off2 - off1) 1 poly2)
+    EQ -> toLaurentNum off1 (poly1 - poly2)
+    GT -> toLaurentNum off2 (Sparse.scale (fromIntegral $ off1 - off2) 1 poly1 - poly2)
+  negate (Laurent off poly) = Laurent off (negate poly)
+  abs = id
+  signum = const 1
+  fromInteger n = Laurent 0 (fromInteger n)
+  {-# INLINE (+) #-}
+  {-# INLINE (-) #-}
+  {-# INLINE negate #-}
+  {-# INLINE fromInteger #-}
+  {-# INLINE (*) #-}
+
+instance (Eq a, Semiring a, G.Vector v (Word, a)) => Semiring (Laurent v a) where
+  zero = Laurent 0 zero
+  one  = Laurent 0 one
+  Laurent off1 poly1 `times` Laurent off2 poly2 =
+    toLaurent (off1 + off2) (poly1 `times` poly2)
+  Laurent off1 poly1 `plus` Laurent off2 poly2 = case off1 `compare` off2 of
+    LT -> toLaurent off1 (poly1 `plus` Sparse.scale' (fromIntegral $ off2 - off1) one poly2)
+    EQ -> toLaurent off1 (poly1 `plus` poly2)
+    GT -> toLaurent off2 (Sparse.scale' (fromIntegral $ off1 - off2) one poly1 `plus` poly2)
+  fromNatural n = Laurent 0 (fromNatural n)
+  {-# INLINE zero #-}
+  {-# INLINE one #-}
+  {-# INLINE plus #-}
+  {-# INLINE times #-}
+  {-# INLINE fromNatural #-}
+
+instance (Eq a, Ring a, G.Vector v (Word, a)) => Ring (Laurent v a) where
+  negate (Laurent off poly) = Laurent off (Semiring.negate poly)
+
+-- | Create a monomial from a power and a coefficient.
+monomial :: (Eq a, Semiring a, G.Vector v (Word, a)) => Int -> a -> Laurent v a
+monomial p c
+  | c == zero = Laurent 0 zero
+  | otherwise = Laurent p (Sparse.monomial' 0 c)
+{-# INLINE monomial #-}
+
+-- | Multiply a polynomial by a monomial, expressed as a power and a coefficient.
+--
+-- >>> scale 2 3 (X^2 + 1) :: ULaurent Int
+-- 3 * X^4 + 3 * X^2
+scale :: (Eq a, Semiring a, G.Vector v (Word, a)) => Int -> a -> Laurent v a -> Laurent v a
+scale yp yc (Laurent off poly) = toLaurent (off + yp) (Sparse.scale' 0 yc poly)
+
+-- | Evaluate at a given point.
+--
+-- >>> eval (X^2 + 1 :: ULaurent Int) 3
+-- 10
+eval :: (Field a, G.Vector v (Word, a)) => Laurent v a -> a -> a
+eval (Laurent off poly) x = Sparse.eval' poly x `times`
+  (if off >= 0 then x Semiring.^ off else quot one x Semiring.^ (- off))
+{-# INLINE eval #-}
+
+-- | Substitute another polynomial instead of 'Data.Poly.Sparse.X'.
+--
+-- >>> subst (X^2 + 1 :: UPoly Int) (X + 1 :: ULaurent Int)
+-- 1 * X^2 + 2 * X + 2
+subst :: (Eq a, Semiring a, G.Vector v (Word, a), G.Vector w (Word, a)) => Poly v a -> Laurent w a -> Laurent w a
+subst = Sparse.substitute' (scale 0)
+{-# INLINE subst #-}
+
+-- | Take a derivative.
+--
+-- >>> deriv (X^3 + 3 * X) :: ULaurent Int
+-- 3 * X^2 + 3
+deriv :: (Eq a, Ring a, G.Vector v (Word, a)) => Laurent v a -> Laurent v a
+deriv (Laurent off (Poly xs)) =
+  toLaurent (off - 1) $ Sparse.toPoly' $ G.map (\(i, x) -> (i, x `times` Semiring.fromIntegral (fromIntegral i + off))) xs
+{-# INLINE deriv #-}
+
+-- | Create an identity polynomial.
+pattern X :: (Eq a, Semiring a, G.Vector v (Word, a), Eq (v (Word, a))) => Laurent v a
+pattern X <- ((==) var -> True)
+  where X = var
+
+var :: forall a v. (Eq a, Semiring a, G.Vector v (Word, a), Eq (v (Word, a))) => Laurent v a
+var
+  | (one :: a) == zero = Laurent 0 zero
+  | otherwise          = Laurent 1 one
+{-# INLINE var #-}
+
+-- | This operator can be applied only to 'X',
+-- but is instrumental to express Laurent polynomials in mathematical fashion:
+--
+-- >>> X + 2 + 3 * X^-1 :: ULaurent Int
+-- 1 * X + 2 + 3 * X^(-1)
+(^-)
+  :: (Eq a, Semiring a, G.Vector v (Word, a), Eq (v (Word, a)))
+  => Laurent v a
+  -> Int
+  -> Laurent v a
+X^-n = monomial (negate n) one
+_^-_ = error "(^-) can be applied only to X"
+
+instance (Eq a, Ring a, GcdDomain a, Eq (v (Word, a)), G.Vector v (Word, a)) => GcdDomain (Laurent v a) where
+  divide (Laurent off1 poly1) (Laurent off2 poly2) =
+    toLaurent (off1 - off2) <$> divide poly1 poly2
+  {-# INLINE divide #-}
+
+  gcd (Laurent _ poly1) (Laurent _ poly2) =
+    toLaurent 0 (gcd poly1 poly2)
+  {-# INLINE gcd #-}
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,7 +7,6 @@
 -- Sparse polynomials with 'Semiring' instance.
 --
 
-{-# LANGUAGE CPP              #-}
 {-# LANGUAGE FlexibleContexts #-}
 {-# LANGUAGE PatternSynonyms  #-}
 
@@ -17,7 +16,6 @@
   , UPoly
   , unPoly
   , leading
-  -- * Semiring interface
   , toPoly
   , monomial
   , scale
@@ -25,27 +23,17 @@
   , eval
   , subst
   , 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.Euclidean (Field)
 import Data.Semiring (Semiring)
 import qualified Data.Vector.Generic as G
 
 import Data.Poly.Internal.Sparse (Poly(..), VPoly, UPoly, leading)
 import qualified Data.Poly.Internal.Sparse as Sparse
-#if MIN_VERSION_semirings(0,4,2)
-import Data.Poly.Internal.Sparse.Field (gcdExt)
+import Data.Poly.Internal.Sparse.Field ()
 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).
@@ -94,7 +82,6 @@
 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.
 --
@@ -102,4 +89,3 @@
 -- 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
@@ -1,4 +1,4 @@
-{-# LANGUAGE CPP                        #-}
+{-# LANGUAGE DataKinds                  #-}
 {-# LANGUAGE FlexibleContexts           #-}
 {-# LANGUAGE FlexibleInstances          #-}
 {-# LANGUAGE GeneralizedNewtypeDeriving #-}
@@ -8,14 +8,13 @@
 
 module Dense
   ( testSuite
+  , ShortPoly(..)
   ) where
 
 import Prelude hiding (gcd, quotRem, rem)
-#if MIN_VERSION_semirings(0,4,2)
 import Data.Euclidean (Euclidean(..), GcdDomain(..))
-#endif
 import Data.Int
-import Data.Maybe
+import Data.Mod
 import Data.Poly
 import qualified Data.Poly.Semiring as S
 import Data.Proxy
@@ -25,30 +24,20 @@
 import qualified Data.Vector.Unboxed as U
 import Test.Tasty
 import Test.Tasty.QuickCheck hiding (scale, numTests)
-import Test.QuickCheck.Classes
 
 import Quaternion
+import TestUtils
 
 instance (Eq a, Semiring a, Arbitrary a, G.Vector v a) => Arbitrary (Poly v a) where
   arbitrary = S.toPoly . G.fromList <$> arbitrary
   shrink = fmap (S.toPoly . G.fromList) . shrink . G.toList . unPoly
 
-#if MIN_VERSION_semirings(0,4,2)
 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 }
-  deriving
-    ( Eq
-    , Show
-    , Semiring
-#if MIN_VERSION_semirings(0,4,2)
-    , GcdDomain
-    , Euclidean
-#endif
-    )
+  deriving (Eq, Show, Semiring, GcdDomain, Euclidean)
 
 instance (Eq a, Semiring a, Arbitrary a, G.Vector v a) => Arbitrary (ShortPoly (Poly v a)) where
   arbitrary = ShortPoly . S.toPoly . G.fromList . (\xs -> take (length xs `mod` 10) xs) <$> arbitrary
@@ -61,95 +50,62 @@
     , lawsTests
     , evalTests
     , derivTests
-#if MIN_VERSION_semirings(0,4,2)
-    , gcdExtTests
-#endif
     ]
 
 lawsTests :: TestTree
 lawsTests = testGroup "Laws"
-  [ semiringTests
-  , ringTests
-  , numTests
-  , euclideanTests
-  , isListTests
-  , showTests
+  $ semiringTests ++ ringTests ++ numTests ++ euclideanTests ++ gcdDomainTests ++ isListTests ++ showTests
+
+semiringTests :: [TestTree]
+semiringTests =
+  [ mySemiringLaws (Proxy :: Proxy (Poly U.Vector ()))
+  , mySemiringLaws (Proxy :: Proxy (Poly U.Vector Int8))
+  , mySemiringLaws (Proxy :: Proxy (Poly V.Vector Integer))
+  , mySemiringLaws (Proxy :: Proxy (Poly U.Vector (Quaternion Int)))
   ]
 
-semiringTests :: TestTree
-semiringTests
-  = testGroup "Semiring"
-  $ map (uncurry testProperty)
-  $ concatMap lawsProperties
-  [ semiringLaws (Proxy :: Proxy (Poly U.Vector ()))
-  , 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 =
+  [ myRingLaws (Proxy :: Proxy (Poly U.Vector ()))
+  , myRingLaws (Proxy :: Proxy (Poly U.Vector Int8))
+  , myRingLaws (Proxy :: Proxy (Poly V.Vector Integer))
+  , myRingLaws (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 Int8))
-  , ringLaws (Proxy :: Proxy (Poly V.Vector Integer))
-  , ringLaws (Proxy :: Proxy (Poly U.Vector (Quaternion Int)))
-#endif
+numTests :: [TestTree]
+numTests =
+  [ myNumLaws (Proxy :: Proxy (Poly U.Vector Int8))
+  , myNumLaws (Proxy :: Proxy (Poly V.Vector Integer))
+  , myNumLaws (Proxy :: Proxy (Poly U.Vector (Quaternion Int)))
   ]
 
-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
+gcdDomainTests :: [TestTree]
+gcdDomainTests =
+  [ myGcdDomainLaws (Proxy :: Proxy (ShortPoly (Poly V.Vector Integer)))
+  , myGcdDomainLaws (Proxy :: Proxy (PolyOverField (Poly V.Vector (Mod 3))))
+  , myGcdDomainLaws (Proxy :: Proxy (PolyOverField (Poly V.Vector Rational)))
   ]
 
-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
+euclideanTests :: [TestTree]
+euclideanTests =
+  [ myEuclideanLaws (Proxy :: Proxy (ShortPoly (Poly V.Vector (Mod 3))))
+  , myEuclideanLaws (Proxy :: Proxy (ShortPoly (Poly V.Vector Rational)))
   ]
 
-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)))
+isListTests :: [TestTree]
+isListTests =
+  [ myIsListLaws (Proxy :: Proxy (Poly U.Vector ()))
+  , myIsListLaws (Proxy :: Proxy (Poly U.Vector Int8))
+  , myIsListLaws (Proxy :: Proxy (Poly V.Vector Integer))
+  , myIsListLaws (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
+showTests :: [TestTree]
+showTests =
+  [ myShowLaws (Proxy :: Proxy (Poly U.Vector ()))
+  , myShowLaws (Proxy :: Proxy (Poly U.Vector Int8))
+  , myShowLaws (Proxy :: Proxy (Poly V.Vector Integer))
+  , myShowLaws (Proxy :: Proxy (Poly U.Vector (Quaternion Int)))
   ]
 
 arithmeticTests :: TestTree
@@ -199,7 +155,8 @@
     \p c -> c /= 0 ==> leading (monomial p c :: UPoly a) === Just (p, c)
   , testProperty "monomial matches reference" $
     \p (c :: a) -> monomial p c === toPoly (V.fromList (monomialRef p c))
-  , testProperty "scale matches multiplication by monomial" $
+  , tenTimesLess $
+    testProperty "scale matches multiplication by monomial" $
     \p c (xs :: UPoly a) -> scale p c xs === monomial p c * xs
   ]
 
@@ -250,13 +207,15 @@
   => Proxy (Poly v a)
   -> [TestTree]
 substTestGroup _ =
-  [ testProperty "subst (p + q) r = subst p r + subst q r" $
+  [ tenTimesLess $ tenTimesLess $ tenTimesLess $
+    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" $
+  , tenTimesLess $ tenTimesLess $ tenTimesLess $
+    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
@@ -273,10 +232,8 @@
 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" $
@@ -287,26 +244,8 @@
     \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)
-  , testProperty "deriv (subst p q) = deriv q * subst (deriv p) q" $
+  , tenTimesLess $ tenTimesLess $ tenTimesLess $
+    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)
-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/DenseLaurent.hs b/test/DenseLaurent.hs
new file mode 100644
--- /dev/null
+++ b/test/DenseLaurent.hs
@@ -0,0 +1,171 @@
+{-# LANGUAGE FlexibleContexts           #-}
+{-# LANGUAGE FlexibleInstances          #-}
+{-# LANGUAGE GeneralizedNewtypeDeriving #-}
+{-# LANGUAGE ScopedTypeVariables        #-}
+
+{-# OPTIONS_GHC -fno-warn-orphans #-}
+
+module DenseLaurent
+  ( testSuite
+  ) where
+
+import Prelude hiding (gcd, quotRem, rem)
+import Data.Euclidean (Euclidean(..), GcdDomain, Field)
+import Data.Int
+import qualified Data.Poly
+import Data.Poly.Laurent
+import Data.Proxy
+import Data.Semiring (Semiring(..))
+import qualified Data.Vector as V
+import qualified Data.Vector.Generic as G
+import qualified Data.Vector.Unboxed as U
+import Test.Tasty
+import Test.Tasty.QuickCheck hiding (scale, numTests)
+
+import Dense (ShortPoly(..))
+import Quaternion
+import TestUtils
+
+instance (Eq a, Semiring a, Arbitrary a, G.Vector v a) => Arbitrary (Laurent v a) where
+  arbitrary = toLaurent <$> ((`rem` 10) <$> arbitrary) <*> arbitrary
+  shrink = fmap (uncurry toLaurent) . shrink . unLaurent
+
+instance (Eq a, Semiring a, Arbitrary a, G.Vector v a) => Arbitrary (LaurentOverField (Laurent v a)) where
+  arbitrary = (LaurentOverField .) . toLaurent <$> ((`rem` 10) <$> arbitrary) <*> (Data.Poly.unPolyOverField <$> arbitrary)
+  shrink = fmap (LaurentOverField . uncurry toLaurent . fmap Data.Poly.unPolyOverField) . shrink . fmap Data.Poly.PolyOverField . unLaurent . unLaurentOverField
+
+newtype ShortLaurent a = ShortLaurent { unShortLaurent :: a }
+  deriving (Eq, Show, Semiring, GcdDomain)
+
+instance (Eq a, Semiring a, Arbitrary a, G.Vector v a) => Arbitrary (ShortLaurent (Laurent v a)) where
+  arbitrary = (ShortLaurent .) . toLaurent <$> ((`rem` 10) <$> arbitrary) <*> (unShortPoly <$> arbitrary)
+  shrink = fmap (ShortLaurent . uncurry toLaurent . fmap unShortPoly) . shrink . fmap ShortPoly . unLaurent . unShortLaurent
+
+testSuite :: TestTree
+testSuite = testGroup "DenseLaurent"
+  [ otherTests
+  , lawsTests
+  , evalTests
+  , derivTests
+  ]
+
+lawsTests :: TestTree
+lawsTests = testGroup "Laws"
+  $ semiringTests ++ ringTests ++ numTests ++ gcdDomainTests ++ showTests
+
+semiringTests :: [TestTree]
+semiringTests =
+  [ mySemiringLaws (Proxy :: Proxy (Laurent U.Vector ()))
+  , mySemiringLaws (Proxy :: Proxy (Laurent U.Vector Int8))
+  , mySemiringLaws (Proxy :: Proxy (Laurent V.Vector Integer))
+  , mySemiringLaws (Proxy :: Proxy (Laurent U.Vector (Quaternion Int)))
+  ]
+
+ringTests :: [TestTree]
+ringTests =
+  [ myRingLaws (Proxy :: Proxy (Laurent U.Vector ()))
+  , myRingLaws (Proxy :: Proxy (Laurent U.Vector Int8))
+  , myRingLaws (Proxy :: Proxy (Laurent V.Vector Integer))
+  , myRingLaws (Proxy :: Proxy (Laurent U.Vector (Quaternion Int)))
+  ]
+
+numTests :: [TestTree]
+numTests =
+  [ myNumLaws (Proxy :: Proxy (Laurent U.Vector Int8))
+  , myNumLaws (Proxy :: Proxy (Laurent V.Vector Integer))
+  , myNumLaws (Proxy :: Proxy (Laurent U.Vector (Quaternion Int)))
+  ]
+
+gcdDomainTests :: [TestTree]
+gcdDomainTests =
+  [ myGcdDomainLaws (Proxy :: Proxy (ShortLaurent (Laurent V.Vector Integer)))
+  , myGcdDomainLaws (Proxy :: Proxy (LaurentOverField (Laurent V.Vector Rational)))
+  ]
+
+showTests :: [TestTree]
+showTests =
+  [ myShowLaws (Proxy :: Proxy (Laurent U.Vector ()))
+  , myShowLaws (Proxy :: Proxy (Laurent U.Vector Int8))
+  , myShowLaws (Proxy :: Proxy (Laurent V.Vector Integer))
+  , myShowLaws (Proxy :: Proxy (Laurent U.Vector (Quaternion Int)))
+  ]
+
+otherTests :: TestTree
+otherTests = testGroup "other" $ concat
+  [ otherTestGroup (Proxy :: Proxy Int8)
+  , otherTestGroup (Proxy :: Proxy (Quaternion Int))
+  ]
+
+otherTestGroup
+  :: forall a.
+     (Eq a, Show a, Semiring a, Num a, Arbitrary a, U.Unbox a, G.Vector U.Vector a)
+  => Proxy a
+  -> [TestTree]
+otherTestGroup _ =
+  [ testProperty "leading p 0 == Nothing" $
+    \p -> leading (monomial p 0 :: ULaurent a) === Nothing
+  , testProperty "leading . monomial = id" $
+    \p c -> c /= 0 ==> leading (monomial p c :: ULaurent a) === Just (p, c)
+  , tenTimesLess $
+    testProperty "scale matches multiplication by monomial" $
+    \p c (xs :: ULaurent a) -> scale p c xs === monomial p c * xs
+  ]
+
+evalTests :: TestTree
+evalTests = testGroup "eval" $ concat
+  [ evalTestGroup  (Proxy :: Proxy (Laurent V.Vector Rational))
+  , substTestGroup (Proxy :: Proxy (Laurent U.Vector Int8))
+  ]
+
+evalTestGroup
+  :: forall v a.
+     (Eq a, Field a, Arbitrary a, Show a, Eq (v a), Show (v a), G.Vector v a)
+  => Proxy (Laurent v a)
+  -> [TestTree]
+evalTestGroup _ =
+  [ testProperty "eval (p + q) r = eval p r + eval q r" $
+    \p q r -> e (p `plus` q) r === e p r `plus` e q r
+  , testProperty "eval (p * q) r = eval p r * eval q r" $
+    \p q r -> e (p `times` q) r === e p r `times` e q r
+  , testProperty "eval x p = p" $
+    \p -> e X p === p
+  , testProperty "eval (monomial 0 c) p = c" $
+    \c p -> e (monomial 0 c) p === c
+  ]
+  where
+    e :: Laurent v a -> a -> a
+    e = 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 (Laurent v a)
+  -> [TestTree]
+substTestGroup _ =
+  [ tenTimesLess $ tenTimesLess $ tenTimesLess $
+    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 Data.Poly.X p === p
+  , testProperty "subst (monomial 0 c) p = monomial 0 c" $
+    \c p -> e (Data.Poly.monomial 0 c) p === monomial 0 c
+  ]
+  where
+    e :: Data.Poly.Poly v a -> Laurent v a -> Laurent v a
+    e = subst
+
+derivTests :: TestTree
+derivTests = testGroup "deriv"
+  [ testProperty "deriv c = 0" $
+    \c -> deriv (monomial 0 c :: Laurent V.Vector Int) === 0
+  , testProperty "deriv cX = c" $
+    \c -> deriv (monomial 0 c * X :: Laurent V.Vector Int) === monomial 0 c
+  , testProperty "deriv (p + q) = deriv p + deriv q" $
+    \p q -> deriv (p + q) === (deriv p + deriv q :: Laurent V.Vector Int)
+  , testProperty "deriv (p * q) = p * deriv q + q * deriv p" $
+    \p q -> deriv (p * q) === (p * deriv q + q * deriv p :: Laurent V.Vector Int)
+  , tenTimesLess $ tenTimesLess $ tenTimesLess $
+    testProperty "deriv (subst p q) = deriv q * subst (deriv p) q" $
+    \(p :: Data.Poly.Poly V.Vector Int) (q :: Laurent U.Vector Int) ->
+      deriv (subst p q) === deriv q * subst (Data.Poly.deriv p) q
+  ]
diff --git a/test/Main.hs b/test/Main.hs
--- a/test/Main.hs
+++ b/test/Main.hs
@@ -3,10 +3,16 @@
 import Test.Tasty
 
 import qualified Dense as Dense
+import qualified DenseLaurent as DenseLaurent
+import qualified Orthogonal as Orthogonal
 import qualified Sparse as Sparse
+import qualified SparseLaurent as SparseLaurent
 
 main :: IO ()
 main = defaultMain $ testGroup "All"
     [ Dense.testSuite
+    , DenseLaurent.testSuite
     , Sparse.testSuite
+    , SparseLaurent.testSuite
+    , Orthogonal.testSuite
     ]
diff --git a/test/Orthogonal.hs b/test/Orthogonal.hs
new file mode 100644
--- /dev/null
+++ b/test/Orthogonal.hs
@@ -0,0 +1,155 @@
+{-# LANGUAGE OverloadedLists #-}
+
+module Orthogonal
+  ( testSuite
+  ) where
+
+import Test.Tasty
+
+import Data.List (foldl', tails)
+import Data.Poly (VPoly, deriv, eval, integral)
+import Data.Poly.Orthogonal
+import Test.Tasty.QuickCheck
+
+testSuite :: TestTree
+testSuite = testGroup "Orthogonal"
+  [ testGroup "differential equations"
+    [ testProperty "jacobi"      prop_jacobi_de
+    , testProperty "gegenbauer"  prop_gegenbauer_de
+    , testProperty "legendre"    prop_legendre_de
+    , testProperty "chebyshev1"  prop_chebyshev1_de
+    , testProperty "chebyshev2"  prop_chebyshev2_de
+    , testProperty "hermitePhys" prop_hermitePhys_de
+    , testProperty "laguerre"    prop_laguerre_de
+    , testProperty "laguerreGen" prop_laguerreGen_de
+    ]
+  , testGroup "normalization"
+    [ testProperty "jacobi"     prop_jacobi_norm
+    , testProperty "gegenbauer" prop_gegenbauer_norm
+    , testProperty "legendre"   prop_legendre_norm
+    , testProperty "chebyshev1" prop_chebyshev1_norm
+    , testProperty "chebyshev2" prop_chebyshev2_norm
+    ]
+  , testGroup "orthogonality"
+    [ testProperty "legendre"   prop_legendre_orth
+    ]
+  , testGroup "Hermite"
+    [ testProperty "hermiteProb" prop_hermiteProb
+    , testProperty "hermitePhys" prop_hermitePhys
+    ]
+  ]
+
+prop_jacobi_de :: Rational -> Rational -> Property
+prop_jacobi_de a b = foldl' (.&&.) (property True) $
+  zipWith (((=== 0) .) . de) [0..limit] (jacobi a b)
+  where
+    de :: Rational -> VPoly Rational -> VPoly Rational
+    de n y = [1, 0, -1] * deriv (deriv y)
+           + [b - a, - (a + b + 2)] * deriv y
+           + [n * (n + a + b + 1)] * y
+
+prop_gegenbauer_de :: Rational -> Property
+prop_gegenbauer_de g = foldl' (.&&.) (property True) $
+  zipWith (((=== 0) .) . de) [0..limit] (gegenbauer g)
+  where
+    de :: Rational -> VPoly Rational -> VPoly Rational
+    de n y = [1, 0, -1] * deriv (deriv y)
+           + [0, - (2 * g + 1)] * deriv y
+           + [n * (n + 2 * g)] * y
+
+prop_legendre_de :: Property
+prop_legendre_de = once $ foldl' (.&&.) (property True) $
+  zipWith (((=== 0) .) . de) [0..limit] legendre
+  where
+    de :: Rational -> VPoly Rational -> VPoly Rational
+    de n y = deriv ([1, 0, -1] * deriv y) + [n * (n + 1)] * y
+
+prop_chebyshev1_de :: Property
+prop_chebyshev1_de = once $ foldl' (.&&.) (property True) $
+  zipWith (((=== 0) .) . de) [0..limit] chebyshev1
+  where
+    de :: Integer -> VPoly Integer -> VPoly Integer
+    de n y = [1, 0, -1] * deriv (deriv y) + [0, -1] * deriv y + [n * n] * y
+
+prop_chebyshev2_de :: Property
+prop_chebyshev2_de = once $ foldl' (.&&.) (property True) $
+  zipWith (((=== 0) .) . de) [0..limit] chebyshev2
+  where
+    de :: Integer -> VPoly Integer -> VPoly Integer
+    de n y = [1, 0, -1] * deriv (deriv y) + [0, -3] * deriv y + [n * (n + 2)] * y
+
+prop_hermitePhys_de :: Property
+prop_hermitePhys_de = once $ foldl' (.&&.) (property True) $
+  zipWith (((=== 0) .) . de) [0..limit] hermitePhys
+  where
+    de :: Integer -> VPoly Integer -> VPoly Integer
+    de n y = deriv (deriv y) + [0, -2] * deriv y + [2 * n] * y
+
+prop_laguerre_de :: Property
+prop_laguerre_de = once $ foldl' (.&&.) (property True) $
+  zipWith (((=== 0) .) . de) [0..limit] laguerre
+  where
+    de :: Rational -> VPoly Rational -> VPoly Rational
+    de n y = [0, 1] * deriv (deriv y) + [1, -1] * deriv y + [n] * y
+
+prop_laguerreGen_de :: Rational -> Property
+prop_laguerreGen_de a  = foldl' (.&&.) (property True) $
+  zipWith (((=== 0) .) . de) [0..limit] (laguerreGen a)
+  where
+    de :: Rational -> VPoly Rational -> VPoly Rational
+    de n y = [0, 1] * deriv (deriv y) + [1 + a, -1] * deriv y + [n] * y
+
+prop_jacobi_norm :: Rational -> Rational -> Property
+prop_jacobi_norm a b = foldl' (.&&.) (property True) $
+  zipWith (\n y -> norm n === eval y 1) [0..limit] (jacobi a b :: [VPoly Rational])
+  where
+    prod n x = product $ take n $ iterate (subtract 1) (fromIntegral n + x)
+    norm n = prod n a / prod n 0
+
+prop_gegenbauer_norm :: Rational -> Property
+prop_gegenbauer_norm a = foldl' (.&&.) (property True) $
+  zipWith (\n y -> norm n === eval y 1) [0..limit] (gegenbauer a :: [VPoly Rational])
+  where
+    prod n x = product $ take n $ iterate (subtract 1) (fromIntegral n + x)
+    norm n = prod n (a - 1 / 2) / prod n 0
+
+prop_legendre_norm :: Property
+prop_legendre_norm = once $ foldl' (.&&.) (property True) $
+  map ((=== 1) . flip eval 1) (take limit legendre :: [VPoly Rational])
+
+prop_chebyshev1_norm :: Property
+prop_chebyshev1_norm = once $ foldl' (.&&.) (property True) $
+  map ((=== 1) . flip eval 1) (take limit chebyshev1 :: [VPoly Integer])
+
+prop_chebyshev2_norm :: Property
+prop_chebyshev2_norm = once $ foldl' (.&&.) (property True) $
+  zipWith (\n y -> n + 1 === eval y 1) [0..limit] (chebyshev2 :: [VPoly Integer])
+
+prop_legendre_orth :: Property
+prop_legendre_orth = once $ foldl' (.&&.) (property True) $
+  [ integral11 (x * y) === 0 | (x : xs) <- tails polys, y <- xs ]
+  where
+    polys :: [VPoly Rational]
+    polys = take limit $ legendre
+
+hermiteProbRef :: [VPoly Integer]
+hermiteProbRef = iterate (\he -> [0, 1] * he - deriv he) 1
+
+hermitePhysRef :: [VPoly Integer]
+hermitePhysRef = iterate (\h -> [0, 2] * h - deriv h) 1
+
+prop_hermiteProb :: Property
+prop_hermiteProb = once $ foldl' (.&&.) (property True) $
+  take limit $ zipWith (===) hermiteProb hermiteProbRef
+
+prop_hermitePhys :: Property
+prop_hermitePhys = once $ foldl' (.&&.) (property True) $
+  take limit $ zipWith (===) hermitePhys hermitePhysRef
+
+integral11 :: VPoly Rational -> Rational
+integral11 x = eval y 1 - eval y (-1)
+  where
+    y = integral x
+
+limit :: Num a => a
+limit = 10
diff --git a/test/Quaternion.hs b/test/Quaternion.hs
--- a/test/Quaternion.hs
+++ b/test/Quaternion.hs
@@ -29,7 +29,7 @@
 import qualified Data.Vector.Generic.Mutable as M
 import Data.Vector.Unboxed (Unbox)
 
-data Quaternion a = Quaternion a a a a
+data Quaternion a = Quaternion !a !a !a !a
   deriving (Eq, Ord, Show, Generic)
 
 instance Ring a => Semiring (Quaternion a) where
diff --git a/test/Sparse.hs b/test/Sparse.hs
--- a/test/Sparse.hs
+++ b/test/Sparse.hs
@@ -1,4 +1,4 @@
-{-# LANGUAGE CPP                        #-}
+{-# LANGUAGE DataKinds                  #-}
 {-# LANGUAGE FlexibleContexts           #-}
 {-# LANGUAGE FlexibleInstances          #-}
 {-# LANGUAGE GeneralizedNewtypeDeriving #-}
@@ -9,16 +9,15 @@
 
 module Sparse
   ( testSuite
+  , ShortPoly(..)
   ) where
 
 import Prelude hiding (gcd, quotRem, rem)
-#if MIN_VERSION_semirings(0,4,2)
 import Data.Euclidean (Euclidean(..), GcdDomain(..))
-#endif
 import Data.Function
 import Data.Int
-import Data.List
-import Data.Maybe
+import Data.List (groupBy, sortOn)
+import Data.Mod
 import Data.Poly.Sparse
 import qualified Data.Poly.Sparse.Semiring as S
 import Data.Proxy
@@ -28,24 +27,16 @@
 import qualified Data.Vector.Unboxed as U
 import Test.Tasty
 import Test.Tasty.QuickCheck hiding (scale, numTests)
-import Test.QuickCheck.Classes
 
 import Quaternion
+import TestUtils
 
 instance (Eq a, Semiring a, Arbitrary a, G.Vector v (Word, a)) => Arbitrary (Poly v a) where
   arbitrary = S.toPoly . G.fromList <$> arbitrary
   shrink = fmap (S.toPoly . G.fromList) . shrink . G.toList . unPoly
 
 newtype ShortPoly a = ShortPoly { unShortPoly :: a }
-  deriving
-    ( Eq
-    , Show
-    , Semiring
-#if MIN_VERSION_semirings(0,4,2)
-    , GcdDomain
-    , Euclidean
-#endif
-    )
+  deriving (Eq, Show, Semiring, GcdDomain, Euclidean)
 
 instance (Eq a, Semiring a, Arbitrary a, G.Vector v (Word, a)) => Arbitrary (ShortPoly (Poly v a)) where
   arbitrary = ShortPoly . S.toPoly . G.fromList . (\xs -> take (length xs `mod` 5) xs) <$> arbitrary
@@ -58,94 +49,68 @@
     , lawsTests
     , evalTests
     , derivTests
-#if MIN_VERSION_semirings(0,4,2)
-    , gcdExtTests
-#endif
     ]
 
 lawsTests :: TestTree
 lawsTests = testGroup "Laws"
-  [ semiringTests
-  , ringTests
-  , numTests
-  , euclideanTests
-  , isListTests
-  , showTests
+  $ semiringTests ++ ringTests ++ numTests ++ euclideanTests ++ gcdDomainTests ++ isListTests ++ showTests
+
+semiringTests :: [TestTree]
+semiringTests =
+  [ mySemiringLaws (Proxy :: Proxy (Poly U.Vector ()))
+  , mySemiringLaws (Proxy :: Proxy (Poly U.Vector Int8))
+  , mySemiringLaws (Proxy :: Proxy (Poly V.Vector Integer))
+  , tenTimesLess
+  $ mySemiringLaws (Proxy :: Proxy (Poly U.Vector (Quaternion Int)))
   ]
 
-semiringTests :: TestTree
-semiringTests
-  = testGroup "Semiring"
-  $ map (uncurry testProperty)
-  $ concatMap lawsProperties
-  [ semiringLaws (Proxy :: Proxy (Poly U.Vector ()))
-  , 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 =
+  [ myRingLaws (Proxy :: Proxy (Poly U.Vector ()))
+  , myRingLaws (Proxy :: Proxy (Poly U.Vector Int8))
+  , myRingLaws (Proxy :: Proxy (Poly V.Vector Integer))
+  , myRingLaws (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 Int8))
-  , ringLaws (Proxy :: Proxy (Poly V.Vector Integer))
-  , ringLaws (Proxy :: Proxy (Poly U.Vector (Quaternion Int)))
-#endif
+numTests :: [TestTree]
+numTests =
+  [ myNumLaws (Proxy :: Proxy (Poly U.Vector Int8))
+  , myNumLaws (Proxy :: Proxy (Poly V.Vector Integer))
+  , tenTimesLess
+  $ myNumLaws (Proxy :: Proxy (Poly U.Vector (Quaternion Int)))
   ]
 
-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
+gcdDomainTests :: [TestTree]
+gcdDomainTests =
+  [ myGcdDomainLaws (Proxy :: Proxy (ShortPoly (Poly V.Vector Integer)))
+  , tenTimesLess
+  $ myGcdDomainLaws (Proxy :: Proxy (ShortPoly (Poly V.Vector (Mod 3))))
+  , tenTimesLess
+  $ myGcdDomainLaws (Proxy :: Proxy (ShortPoly (Poly V.Vector Rational)))
   ]
 
-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
+euclideanTests :: [TestTree]
+euclideanTests =
+  [ myEuclideanLaws (Proxy :: Proxy (ShortPoly (Poly V.Vector (Mod 3))))
+  , myEuclideanLaws (Proxy :: Proxy (ShortPoly (Poly V.Vector Rational)))
   ]
 
-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)))
+isListTests :: [TestTree]
+isListTests =
+  [ myIsListLaws (Proxy :: Proxy (Poly U.Vector ()))
+  , myIsListLaws (Proxy :: Proxy (Poly U.Vector Int8))
+  , myIsListLaws (Proxy :: Proxy (Poly V.Vector Integer))
+  , tenTimesLess
+  $ myIsListLaws (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
+showTests :: [TestTree]
+showTests =
+  [ myShowLaws (Proxy :: Proxy (Poly U.Vector ()))
+  , myShowLaws (Proxy :: Proxy (Poly U.Vector Int8))
+  , myShowLaws (Proxy :: Proxy (Poly V.Vector Integer))
+  , tenTimesLess
+  $ myShowLaws (Proxy :: Proxy (Poly U.Vector (Quaternion Int)))
   ]
 
 arithmeticTests :: TestTree
@@ -156,7 +121,8 @@
   , testProperty "subtraction matches reference" $
     \(xs :: [(Word, Int)]) ys -> toPoly (V.fromList (subRef xs ys)) ===
       toPoly (V.fromList xs) - toPoly (V.fromList ys)
-  , testProperty "multiplication matches reference" $
+  , tenTimesLess $
+    testProperty "multiplication matches reference" $
     \(xs :: [(Word, Int)]) ys -> toPoly (V.fromList (mulRef xs ys)) ===
       toPoly (V.fromList xs) * toPoly (V.fromList ys)
   ]
@@ -204,7 +170,8 @@
     \p c -> c /= 0 ==> leading (monomial p c :: UPoly a) === Just (p, c)
   , testProperty "monomial matches reference" $
     \p (c :: a) -> monomial p c === toPoly (V.fromList (monomialRef p c))
-  , testProperty "scale matches multiplication by monomial" $
+  , tenTimesLess $
+    testProperty "scale matches multiplication by monomial" $
     \p c (xs :: UPoly a) -> scale p c xs === monomial p c * xs
   ]
 
@@ -274,12 +241,10 @@
 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" $
@@ -292,22 +257,3 @@
   --   \(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)
-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/SparseLaurent.hs b/test/SparseLaurent.hs
new file mode 100644
--- /dev/null
+++ b/test/SparseLaurent.hs
@@ -0,0 +1,177 @@
+{-# LANGUAGE FlexibleContexts           #-}
+{-# LANGUAGE FlexibleInstances          #-}
+{-# LANGUAGE GeneralizedNewtypeDeriving #-}
+{-# LANGUAGE ScopedTypeVariables        #-}
+{-# LANGUAGE UndecidableInstances       #-}
+
+{-# OPTIONS_GHC -fno-warn-orphans #-}
+
+module SparseLaurent
+  ( testSuite
+  ) where
+
+import Prelude hiding (gcd, quotRem, rem)
+import Data.Euclidean (Euclidean(..), GcdDomain(..), Field)
+import Data.Int
+import qualified Data.Poly.Sparse
+import Data.Poly.Sparse.Laurent
+import Data.Proxy
+import Data.Semiring (Semiring(..))
+import qualified Data.Vector as V
+import qualified Data.Vector.Generic as G
+import qualified Data.Vector.Unboxed as U
+import Test.Tasty
+import Test.Tasty.QuickCheck hiding (scale, numTests)
+
+import Quaternion
+import Sparse (ShortPoly(..))
+import TestUtils
+
+instance (Eq a, Semiring a, Arbitrary a, G.Vector v (Word, a)) => Arbitrary (Laurent v a) where
+  arbitrary = toLaurent <$> ((`rem` 10) <$> arbitrary) <*> arbitrary
+  shrink = fmap (uncurry toLaurent) . shrink . unLaurent
+
+newtype ShortLaurent a = ShortLaurent { unShortLaurent :: a }
+  deriving (Eq, Show, Semiring, GcdDomain)
+
+instance (Eq a, Semiring a, Arbitrary a, G.Vector v (Word, a)) => Arbitrary (ShortLaurent (Laurent v a)) where
+  arbitrary = (ShortLaurent .) . toLaurent <$> ((`rem` 10) <$> arbitrary) <*> (unShortPoly <$> arbitrary)
+  shrink = fmap (ShortLaurent . uncurry toLaurent . fmap unShortPoly) . shrink . fmap ShortPoly . unLaurent . unShortLaurent
+
+testSuite :: TestTree
+testSuite = testGroup "SparseLaurent"
+    [ otherTests
+    , lawsTests
+    , evalTests
+    , derivTests
+    ]
+
+lawsTests :: TestTree
+lawsTests = testGroup "Laws"
+  $ semiringTests ++ ringTests ++ numTests ++ gcdDomainTests ++ isListTests ++ showTests
+
+semiringTests :: [TestTree]
+semiringTests =
+  [ mySemiringLaws (Proxy :: Proxy (Laurent U.Vector ()))
+  , mySemiringLaws (Proxy :: Proxy (Laurent U.Vector Int8))
+  , mySemiringLaws (Proxy :: Proxy (Laurent V.Vector Integer))
+  , tenTimesLess
+  $ mySemiringLaws (Proxy :: Proxy (Laurent U.Vector (Quaternion Int)))
+  ]
+
+ringTests :: [TestTree]
+ringTests =
+  [ myRingLaws (Proxy :: Proxy (Laurent U.Vector ()))
+  , myRingLaws (Proxy :: Proxy (Laurent U.Vector Int8))
+  , myRingLaws (Proxy :: Proxy (Laurent V.Vector Integer))
+  , myRingLaws (Proxy :: Proxy (Laurent U.Vector (Quaternion Int)))
+  ]
+
+numTests :: [TestTree]
+numTests =
+  [ myNumLaws (Proxy :: Proxy (Laurent U.Vector Int8))
+  , myNumLaws (Proxy :: Proxy (Laurent V.Vector Integer))
+  , tenTimesLess
+  $ myNumLaws (Proxy :: Proxy (Laurent U.Vector (Quaternion Int)))
+  ]
+
+gcdDomainTests :: [TestTree]
+gcdDomainTests =
+  [ myGcdDomainLaws (Proxy :: Proxy (ShortLaurent (Laurent V.Vector Integer)))
+  , tenTimesLess
+  $ myGcdDomainLaws (Proxy :: Proxy (ShortLaurent (Laurent V.Vector Rational)))
+  ]
+
+isListTests :: [TestTree]
+isListTests =
+  [ myIsListLaws (Proxy :: Proxy (Laurent U.Vector ()))
+  , myIsListLaws (Proxy :: Proxy (Laurent U.Vector Int8))
+  , myIsListLaws (Proxy :: Proxy (Laurent V.Vector Integer))
+  , tenTimesLess
+  $ myIsListLaws (Proxy :: Proxy (Laurent U.Vector (Quaternion Int)))
+  ]
+
+showTests :: [TestTree]
+showTests =
+  [ myShowLaws (Proxy :: Proxy (Laurent U.Vector ()))
+  , myShowLaws (Proxy :: Proxy (Laurent U.Vector Int8))
+  , myShowLaws (Proxy :: Proxy (Laurent V.Vector Integer))
+  , tenTimesLess
+  $ myShowLaws (Proxy :: Proxy (Laurent U.Vector (Quaternion Int)))
+  ]
+
+otherTests :: TestTree
+otherTests = testGroup "other" $ concat
+  [ otherTestGroup (Proxy :: Proxy Int8)
+  , otherTestGroup (Proxy :: Proxy (Quaternion Int))
+  ]
+
+otherTestGroup
+  :: forall a.
+     (Eq a, Show a, Semiring a, Num a, Arbitrary a, U.Unbox a, G.Vector U.Vector a)
+  => Proxy a
+  -> [TestTree]
+otherTestGroup _ =
+  [ testProperty "leading p 0 == Nothing" $
+    \p -> leading (monomial p 0 :: ULaurent a) === Nothing
+  , testProperty "leading . monomial = id" $
+    \p c -> c /= 0 ==> leading (monomial p c :: ULaurent a) === Just (p, c)
+  , tenTimesLess $
+    testProperty "scale matches multiplication by monomial" $
+    \p c (xs :: ULaurent a) -> scale p c xs === monomial p c * xs
+  ]
+
+evalTests :: TestTree
+evalTests = testGroup "eval" $ concat
+  [ evalTestGroup  (Proxy :: Proxy (Laurent V.Vector Rational))
+  , substTestGroup (Proxy :: Proxy (Laurent U.Vector Int8))
+  ]
+
+evalTestGroup
+  :: forall v a.
+     (Eq a, Field a, Arbitrary a, Show a, Eq (v (Word, a)), Show (v (Word, a)), G.Vector v (Word, a))
+  => Proxy (Laurent v a)
+  -> [TestTree]
+evalTestGroup _ =
+  [ testProperty "eval (p + q) r = eval p r + eval q r" $
+    \p q r -> e (p `plus` q) r === e p r `plus` e q r
+  , testProperty "eval (p * q) r = eval p r * eval q r" $
+    \p q r -> e (p `times` q) r === e p r `times` e q r
+  , testProperty "eval x p = p" $
+    \p -> e X p === p
+  , testProperty "eval (monomial 0 c) p = c" $
+    \c p -> e (monomial 0 c) p === c
+  ]
+  where
+    e :: Laurent v a -> a -> a
+    e = 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 (Laurent v a)
+  -> [TestTree]
+substTestGroup _ =
+  [ testProperty "subst x p = p" $
+    \p -> e Data.Poly.Sparse.X p === p
+  , testProperty "subst (monomial 0 c) p = monomial 0 c" $
+    \c p -> e (Data.Poly.Sparse.monomial 0 c) p === monomial 0 c
+  ]
+  where
+    e :: Data.Poly.Sparse.Poly v a -> Laurent v a -> Laurent v a
+    e = subst
+
+derivTests :: TestTree
+derivTests = testGroup "deriv"
+  [ testProperty "deriv c = 0" $
+    \c -> deriv (monomial 0 c :: Laurent V.Vector Int) === 0
+  , testProperty "deriv cX = c" $
+    \c -> deriv (monomial 0 c * X :: Laurent V.Vector Int) === monomial 0 c
+  , testProperty "deriv (p + q) = deriv p + deriv q" $
+    \p q -> deriv (p + q) === (deriv p + deriv q :: Laurent V.Vector Int)
+  , testProperty "deriv (p * q) = p * deriv q + q * deriv p" $
+    \p q -> deriv (p * q) === (p * deriv q + q * deriv p :: Laurent V.Vector Int)
+  -- , testProperty "deriv (subst p q) = deriv q * subst (deriv p) q" $
+  --   \(p :: Laurent V.Vector Int) (q :: Laurent U.Vector Int) ->
+  --     deriv (subst p q) === deriv q * subst (deriv p) q
+  ]
diff --git a/test/TestUtils.hs b/test/TestUtils.hs
new file mode 100644
--- /dev/null
+++ b/test/TestUtils.hs
@@ -0,0 +1,113 @@
+{-# LANGUAGE CPP              #-}
+{-# LANGUAGE FlexibleContexts #-}
+
+{-# OPTIONS_GHC -fno-warn-orphans #-}
+
+module TestUtils
+  ( tenTimesLess
+  , mySemiringLaws
+  , myRingLaws
+  , myNumLaws
+  , myGcdDomainLaws
+  , myEuclideanLaws
+  , myIsListLaws
+  , myShowLaws
+  ) where
+
+import Data.Euclidean
+import Data.Mod
+import Data.Proxy
+import Data.Semiring (Semiring, Ring)
+import GHC.Exts
+import Test.QuickCheck.Classes
+import Test.Tasty
+import Test.Tasty.QuickCheck
+
+#if MIN_VERSION_base(4,10,0)
+import GHC.TypeNats (KnownNat)
+#else
+import GHC.TypeLits (KnownNat)
+#endif
+
+instance KnownNat m => Arbitrary (Mod m) where
+  arbitrary = oneof [arbitraryBoundedEnum, fromInteger <$> arbitrary]
+  shrink = map fromInteger . shrink . toInteger . unMod
+
+tenTimesLess :: TestTree -> TestTree
+tenTimesLess = adjustOption $
+  \(QuickCheckTests n) -> QuickCheckTests (max 100 (n `div` 10))
+
+mySemiringLaws :: (Eq a, Semiring a, Arbitrary a, Show a) => Proxy a -> TestTree
+mySemiringLaws proxy = testGroup tpclss $ map tune props
+  where
+    Laws tpclss props = semiringLaws proxy
+
+    tune pair = case fst pair of
+      "Multiplicative Associativity" ->
+        tenTimesLess test
+      "Multiplication Left Distributes Over Addition" ->
+        tenTimesLess test
+      "Multiplication Right Distributes Over Addition" ->
+        tenTimesLess test
+      _ -> test
+      where
+        test = uncurry testProperty pair
+
+myRingLaws :: (Eq a, Ring a, Arbitrary a, Show a) => Proxy a -> TestTree
+myRingLaws proxy = testGroup tpclss $ map (uncurry testProperty) props
+  where
+    Laws tpclss props = ringLaws proxy
+
+myNumLaws :: (Eq a, Num a, Arbitrary a, Show a) => Proxy a -> TestTree
+myNumLaws proxy = testGroup tpclss $ map tune props
+  where
+    Laws tpclss props = numLaws proxy
+
+    tune pair = case fst pair of
+      "Multiplicative Associativity" ->
+        tenTimesLess test
+      "Multiplication Left Distributes Over Addition" ->
+        tenTimesLess test
+      "Multiplication Right Distributes Over Addition" ->
+        tenTimesLess test
+      "Subtraction" ->
+        tenTimesLess test
+      _ -> test
+      where
+        test = uncurry testProperty pair
+
+myGcdDomainLaws :: (Eq a, GcdDomain a, Arbitrary a, Show a) => Proxy a -> TestTree
+myGcdDomainLaws proxy = testGroup tpclss $ map tune props
+  where
+    Laws tpclss props = gcdDomainLaws proxy
+
+    tune pair = case fst pair of
+      "gcd1"    -> tenTimesLess $ tenTimesLess test
+      "gcd2"    -> tenTimesLess $ tenTimesLess test
+      "lcm1"    -> tenTimesLess $ tenTimesLess $ tenTimesLess test
+      "lcm2"    -> tenTimesLess test
+      "coprime" -> tenTimesLess $ tenTimesLess test
+      _ -> test
+      where
+        test = uncurry testProperty pair
+
+myEuclideanLaws :: (Eq a, Euclidean a, Arbitrary a, Show a) => Proxy a -> TestTree
+myEuclideanLaws proxy = testGroup tpclss $ map (uncurry testProperty) props
+  where
+    Laws tpclss props = euclideanLaws proxy
+
+myIsListLaws :: (Eq a, IsList a, Arbitrary a, Show a, Show (Item a), Arbitrary (Item a)) => Proxy a -> TestTree
+myIsListLaws proxy = testGroup tpclss $ map (uncurry testProperty) props
+  where
+    Laws tpclss props = isListLaws proxy
+
+myShowLaws :: (Eq a, Arbitrary a, Show a) => Proxy a -> TestTree
+myShowLaws proxy = testGroup tpclss $ map tune props
+  where
+    Laws tpclss props = showLaws proxy
+
+    tune pair = case fst pair of
+      "Equivariance: showList" -> tenTimesLess $ tenTimesLess test
+      _ -> test
+      where
+        test = uncurry testProperty pair