diff --git a/doctest/DoctestDriver.hs b/doctest/DoctestDriver.hs
new file mode 100644
--- /dev/null
+++ b/doctest/DoctestDriver.hs
@@ -0,0 +1,12 @@
+{-# LANGUAGE CPP #-}
+
+#if MIN_VERSION_GLASGOW_HASKELL(8,4,4,0)
+{-# OPTIONS_GHC -F -pgmF doctest-discover #-}
+#else
+module Main where
+
+import qualified System.IO as IO
+
+main :: IO ()
+main = IO.putStrLn "WARNING: doctest will not run on GHC versions earlier than 8.4.4"
+#endif
diff --git a/relation.cabal b/relation.cabal
--- a/relation.cabal
+++ b/relation.cabal
@@ -1,7 +1,7 @@
 cabal-version:      2.2
 
 name:               relation
-version:            0.5
+version:            0.5.1.0
 synopsis:           A data structure representing Relations on Sets.
 description:        A library to model relationships between two objects that are subclasses of Ord.
 
@@ -18,7 +18,7 @@
 category:           Data Structures
 stability:          Experimental
 build-type:         Simple
-tested-with:        GHC==7.4
+tested-with:        GHC == 8.2.2, GHC == 8.4.4, GHC == 8.6.5, GHC == 8.8.3
 extra-source-files: LICENSE
                     README.md
 
@@ -29,17 +29,22 @@
 common base                 { build-depends: base                 >= 4          && < 5      }
 
 common containers           { build-depends: containers           >= 0.5        && < 0.7    }
-common hedgehog             { build-depends: hedgehog             >= 0.5        && < 0.7    }
+common doctest              { build-depends: doctest              >= 0.16.2     && < 0.17   }
+common doctest-discover     { build-depends: doctest-discover     >= 0.2        && < 0.3    }
+common hedgehog             { build-depends: hedgehog             >= 0.5        && < 1.1    }
 common hspec                { build-depends: hspec                >= 2.4        && < 3      }
 common hw-hspec-hedgehog    { build-depends: hw-hspec-hedgehog    >= 0.1.0.4    && < 0.2    }
 
-common common
+common relation
+  build-depends: relation
+
+common config
   default-language:   Haskell2010
   ghc-options:        -Wall -O2
 
 library
-  import:   base, common
-          , containers
+  import:             base, config
+                    , containers
   hs-source-dirs:     src
   exposed-modules:    Data.Relation
                     , Data.Relation.Ops
@@ -47,18 +52,41 @@
                     , Data.Relation.Internal.Set
 
 test-suite relation-test
-  import:   base, common
-          , containers
-          , hedgehog
-          , hspec
-          , hw-hspec-hedgehog
+  import:             base, config
+                    , containers
+                    , hedgehog
+                    , hspec
+                    , hw-hspec-hedgehog
+                    , relation
   build-depends:      relation
   type:               exitcode-stdio-1.0
   main-is:            Spec.hs
-  build-depends:      relation
-  other-modules:      Data.RelationSpec
+  other-modules:      Data.Relation.Gen
+                    , Data.RelationSpec
                     , Paths_relation
   autogen-modules:    Paths_relation
   hs-source-dirs:     test
   ghc-options:        -threaded -rtsopts -with-rtsopts=-N
   build-tool-depends: hspec-discover:hspec-discover
+
+test-suite relation-doctest
+  import:               base, config
+                      , doctest
+                      , doctest-discover
+                      , relation
+  type:                 exitcode-stdio-1.0
+  ghc-options:          -threaded
+  main-is:              DoctestDriver.hs
+  HS-Source-Dirs:       doctest
+  build-tool-depends:   doctest-discover:doctest-discover
+
+test-suite doctest
+  import:               base, config
+                      , doctest
+                      , doctest-discover
+                      , relation
+  type:                 exitcode-stdio-1.0
+  ghc-options:          -threaded
+  main-is:              DoctestDriver.hs
+  HS-Source-Dirs:       doctest
+  build-tool-depends:   doctest-discover:doctest-discover
diff --git a/src/Data/Relation.hs b/src/Data/Relation.hs
--- a/src/Data/Relation.hs
+++ b/src/Data/Relation.hs
@@ -55,6 +55,8 @@
   , restrictRan     -- Restrict the range to that of the provided set
   , withoutDom      -- Restrict the domain to exclude elements of the provided set
   , withoutRan      -- Restrict the range to exclude elements of the provided set
+  , (<-<)           -- Compose two relations
+  , (>->)
 
     -- ** Conversion
   , toList          -- Construct a list from a relation
@@ -64,6 +66,7 @@
   ) where
 
 import Control.Monad          (MonadPlus, guard)
+import Data.Foldable          (fold)
 import Data.Functor           (Functor ((<$)))
 import Data.Map               (Map)
 import Data.Maybe             (fromMaybe)
@@ -81,15 +84,27 @@
 
 -- The size is calculated using the domain.
 -- |  @size r@ returns the number of tuples in the relation.
+-- >>> size (fromList []) == 0
+-- True
+-- >>> size (fromList [('a', 1)]) == 1
+-- True
 size :: Relation a b -> Int
 size r = M.foldr ((+) . S.size) 0 (R.domain r)
 
 -- | Construct a relation with no elements.
+-- >>> toList (fromList []) == []
+-- True
 empty :: Relation a b
 empty = Relation M.empty M.empty
 
 -- |
 -- The list must be formatted like: [(k1, v1), (k2, v2),..,(kn, vn)].
+-- >>> toList (fromList [('a', 1)]) == [('a', 1)]
+-- True
+-- >>> fromList [('a', 1), ('a', 1)] == fromList [('a', 1), ('a', 1)]
+-- True
+-- >>> fromList [('a', 1), ('b', 1)] == fromList [('a', 1), ('b', 1)]
+-- True
 fromList :: (Ord a, Ord b) => [(a, b)] -> Relation a b
 fromList xs = Relation
   { R.domain  = M.fromListWith S.union $ snd2Set    xs
@@ -98,13 +113,22 @@
   where snd2Set    = map (\(x, y) -> (x, S.singleton y))
         flipAndSet = map (\(x, y) -> (y, S.singleton x))
 
--- |
--- Builds a List from a Relation.
+-- | Builds a List from a Relation.
+-- >>> toList (fromList [('a', 1)]) == [('a', 1)]
+-- True
+-- >>> toList (fromList [('a', 1), ('b', 2)]) == [('a', 1), ('b', 2)]
+-- True
+-- >>> toList (fromList [('a', 1), ('a', 2)]) == [('a', 1), ('a', 2)]
+-- True
+-- >>> toList (fromList [('a', 1), ('a', 1)]) == [('a', 1)]
+-- True
 toList :: Relation a b -> [(a, b)]
 toList r = concatMap (\(x, y) -> zip (repeat x) (S.toList y)) (M.toList . R.domain $ r)
 
 -- |
 -- Builds a 'Relation' consiting of an association between: @x@ and @y@.
+-- >>> singleton 'a' 1 == fromList [('a', 1)]
+-- True
 singleton :: a -> b -> Relation a b
 singleton x y  = Relation
   { R.domain  = M.singleton x (S.singleton y)
@@ -112,6 +136,12 @@
   }
 
 -- | The 'Relation' that results from the union of two relations: @r@ and @s@.
+-- >>> fromList [('a', 1)] `union` fromList [('a', 1)] == fromList [('a', 1)]
+-- True
+-- >>> fromList [('a', 2)] `union` fromList [('a', 2)] == fromList [('a', 2)]
+-- True
+-- >>> fromList [('a', 1)] `union` fromList [('b', 2)] == fromList [('a', 1), ('b', 2)]
+-- True
 union :: (Ord a, Ord b) => Relation a b -> Relation a b -> Relation a b
 union r s = Relation
   { R.domain  = M.unionWith S.union (R.domain r) (R.domain s)
@@ -119,11 +149,27 @@
   }
 
 -- | Union a list of relations using the 'empty' relation.
+-- >>> unions [] == fromList []
+-- True
+-- >>> unions [fromList [('a', 1)]] == fromList [('a', 1)]
+-- True
+-- >>> unions [fromList [('a', 1)], fromList [('a', 1)]] == fromList [('a', 1)]
+-- True
+-- >>> unions [fromList [('a', 2)], fromList [('a', 2)]] == fromList [('a', 2)]
+-- True
+-- >>> unions [fromList [('a', 1)], fromList [('b', 2)]] == fromList [('a', 1), ('b', 2)]
+-- True
 unions :: (Ord a, Ord b) => [Relation a b] -> Relation a b
 unions = F.foldl' union empty
 
 -- | Intersection of two relations: @a@ and @b@ are related by @intersection r
 -- s@ exactly when @a@ and @b@ are related by @r@ and @s@.
+-- >>> fromList [('a', 1)] `intersection` fromList [('a', 1)] == fromList [('a', 1)]
+-- True
+-- >>> fromList [('a', 2)] `intersection` fromList [('a', 2)] == fromList [('a', 2)]
+-- True
+-- >>> fromList [('a', 1)] `intersection` fromList [('b', 2)] == fromList []
+-- True
 intersection :: (Ord a, Ord b) => Relation a b -> Relation a b -> Relation a b
 intersection r s = Relation
   { R.domain = doubleIntersect (R.domain r) (R.domain s)
@@ -230,3 +276,22 @@
   { R.domain  = M.mapMaybe (S.justUnlessEmpty . flip S.difference s) (R.domain r)
   , R.range   = M.withoutKeys (R.range r) s
   }
+
+-- | Compose two relations: right to left version.
+infixr 9 <-<
+(<-<) :: (Ord a, Ord b, Ord c) => Relation b c -> Relation a b -> Relation a c
+a <-< b = Relation
+  (compose (R.domain a) (R.domain b))
+  (compose (R.range b) (R.range a))
+ where
+  compose a' = M.mapMaybe
+    (S.justUnlessEmpty
+    . fold
+    . M.intersection a'
+    . M.fromSet (const ())
+    )
+
+-- | Compose two relations: left to right version.
+infixl 9 >->
+(>->) :: (Ord a, Ord b, Ord c) => Relation a b -> Relation b c -> Relation a c
+(>->) = flip (<-<)
diff --git a/test/Data/Relation/Gen.hs b/test/Data/Relation/Gen.hs
new file mode 100644
--- /dev/null
+++ b/test/Data/Relation/Gen.hs
@@ -0,0 +1,13 @@
+module Data.Relation.Gen
+  ( relation
+  ) where
+
+import Data.Relation
+import Hedgehog
+import Hedgehog.Range (Range)
+
+import qualified Hedgehog.Gen as G
+
+-- | Generate a relation given generators for the domain and range.
+relation :: (MonadGen m, Ord a, Ord b) => Range Int -> m a -> m b -> m (Relation a b)
+relation r ga gb = fmap fromList (G.list r ((,) <$> ga <*> gb))
diff --git a/test/Data/RelationSpec.hs b/test/Data/RelationSpec.hs
--- a/test/Data/RelationSpec.hs
+++ b/test/Data/RelationSpec.hs
@@ -7,10 +7,12 @@
 import Hedgehog
 import Test.Hspec
 
+import           Control.Monad (replicateM)
 import qualified Data.List              as L
 import qualified Data.Map               as M
 import qualified Data.Relation          as DR
 import qualified Data.Relation.Internal as DR
+import qualified Data.Relation.Gen      as GR
 import qualified Data.Set               as S
 import qualified Hedgehog.Gen           as G
 import qualified Hedgehog.Range         as R
@@ -193,3 +195,9 @@
         <*> G.alpha
       let r = DR.fromList as
       DR.withoutRan (DR.ran r) r === DR.empty
+    it "Compose associatively" $ require $ property $ do
+      ~[a,b,c] <- forAll $ replicateM 3 $ GR.relation
+        (R.linear 10 40)
+        (G.integral (R.linear (1 :: Integer) 10))
+        (G.integral (R.linear (1 :: Integer) 10))
+      ((a DR.<-< b) DR.<-< c) === (a DR.<-< (b DR.<-< c))