diff --git a/Data/ChurchList.hs b/Data/ChurchList.hs
--- a/Data/ChurchList.hs
+++ b/Data/ChurchList.hs
@@ -97,3 +97,7 @@
 fromMaybe :: Maybe a -> ChurchList a
 fromMaybe Nothing = nil
 fromMaybe (Just x) = unit x
+
+{-# INLINE null #-}
+null :: ChurchList a -> Bool
+null = foldr (\_ _ -> False) True
diff --git a/Data/DynamicArray.hs b/Data/DynamicArray.hs
--- a/Data/DynamicArray.hs
+++ b/Data/DynamicArray.hs
@@ -41,24 +41,33 @@
     "}"
 
 -- | Create an empty array.
-newArray :: Default a => Array a
+newArray :: Array a
 newArray = runST $ do
-  marr <- P.newSmallArray 0 def
+  marr <- P.newSmallArray 0 undefined
   arr  <- P.unsafeFreezeSmallArray marr
   return (Array 0 arr)
 
 -- | Index into an array. O(1) time.
 {-# INLINE (!) #-}
 (!) :: Default a => Array a -> Int -> a
-arr ! n
+arr ! n = getWithDefault def n arr
+
+-- | Index into an array. O(1) time.
+{-# INLINE getWithDefault #-}
+getWithDefault :: a -> Int -> Array a -> a
+getWithDefault def n arr
   | 0 <= n && n < arraySize arr =
     P.indexSmallArray (arrayContents arr) n
   | otherwise = def
 
 -- | Update the array. O(n) time.
-{-# INLINEABLE update #-}
+{-# INLINE update #-}
 update :: Default a => Int -> a -> Array a -> Array a
-update n x arr = runST $ do
+update n x arr = updateWithDefault def n x arr
+
+{-# INLINEABLE updateWithDefault #-}
+updateWithDefault :: a -> Int -> a -> Array a -> Array a
+updateWithDefault def n x arr = runST $ do
   let size = arraySize arr `max` (n+1)
   marr <- P.newSmallArray size def
   P.copySmallArray marr 0 (arrayContents arr) 0 (arraySize arr)
diff --git a/Data/Heap.hs b/Data/Heap.hs
--- a/Data/Heap.hs
+++ b/Data/Heap.hs
@@ -2,44 +2,39 @@
 
 {-# LANGUAGE BangPatterns, ScopedTypeVariables #-}
 module Data.Heap(
-  Heap, empty, singleton, insert, removeMin, union, mapMaybe, size) where
+  Heap, empty, singleton, insert, removeMin, union, mapMaybe, size, toList) where
 
 -- | A heap.
 
--- Representation: the size of the heap, and the heap itself.
-data Heap a = Heap {-# UNPACK #-} !Int !(Heap1 a) deriving Show
 -- N.B.: arguments are not strict so code has to take care
 -- to force stuff appropriately.
-data Heap1 a = Nil | Node a (Heap1 a) (Heap1 a) deriving Show
+-- The Int field is the size of the heap.
+data Heap a = Nil | Node {-# UNPACK #-} !Int a (Heap a) (Heap a) deriving Show
 
 -- | Take the union of two heaps.
 {-# INLINEABLE union #-}
-union :: Ord a => Heap a -> Heap a -> Heap a
-union (Heap n1 h1) (Heap n2 h2) = Heap (n1+n2) (union1 h1 h2)
-
-{-# INLINEABLE union1 #-}
-union1 :: forall a. Ord a => Heap1 a -> Heap1 a -> Heap1 a
-union1 = u1
+union :: forall a. Ord a => Heap a -> Heap a -> Heap a
+union = u
   where
     -- The generated code is better when we do everything
-    -- through this u1 function instead of union1...
-    -- This is because u1 has no Ord constraint in its type.
-    u1 :: Heap1 a -> Heap1 a -> Heap1 a
-    u1 Nil h = h
-    u1 h Nil = h
-    u1 h1@(Node x1 l1 r1) h2@(Node x2 l2 r2)
-      | x1 <= x2 = (Node x1 $! u1 r1 h2) l1
-      | otherwise = (Node x2 $! u1 r2 h1) l2
+    -- through this u function instead of union...
+    -- This is because u has no Ord constraint in its type.
+    u :: Heap a -> Heap a -> Heap a
+    u Nil h = h
+    u h Nil = h
+    u h1@(Node s1 x1 l1 r1) h2@(Node s2 x2 l2 r2)
+      | x1 <= x2 = (Node (s1+s2) x1 $! u r1 h2) l1
+      | otherwise = (Node (s1+s2) x2 $! u r2 h1) l2
 
 -- | A singleton heap.
 {-# INLINE singleton #-}
 singleton :: a -> Heap a
-singleton !x = Heap 1 (Node x Nil Nil)
+singleton !x = Node 1 x Nil Nil
 
 -- | The empty heap.
 {-# INLINE empty #-}
 empty :: Heap a
-empty = Heap 0 Nil
+empty = Nil
 
 -- | Insert an element.
 {-# INLINEABLE insert #-}
@@ -49,60 +44,48 @@
 -- | Find and remove the minimum element.
 {-# INLINEABLE removeMin #-}
 removeMin :: Ord a => Heap a -> Maybe (a, Heap a)
-removeMin (Heap _ Nil) = Nothing
-removeMin (Heap n (Node x l r)) = Just (x, Heap (n-1) (union1 l r))
+removeMin Nil = Nothing
+removeMin (Node _ x l r) = Just (x, union l r)
 
+-- | Get the elements of a heap as a list, in unspecified order.
+toList :: Heap a -> [a]
+toList h = tl h []
+  where
+    tl Nil = id
+    tl (Node _ x l r) = (x:) . tl l . tl r
+
 -- | Map a function over a heap, removing all values which
 -- map to 'Nothing'. May be more efficient when the function
 -- being mapped is mostly monotonic.
 {-# INLINEABLE mapMaybe #-}
-mapMaybe :: Ord b => (a -> Maybe b) -> Heap a -> Heap b
-mapMaybe f (Heap _ h) = Heap (sz 0 h') h'
+mapMaybe :: forall a b. Ord b => (a -> Maybe b) -> Heap a -> Heap b
+mapMaybe f h = mm h
   where
-    -- Compute the size fairly efficiently.
-    sz !n Nil = n
-    sz !n (Node _ l r) = sz (sz (n+1) l) r
-
-    h' = mm h
-
+    mm :: Heap a -> Heap b
     mm Nil = Nil
-    mm (Node x l r) =
+    mm (Node _ x l r) =
       case f x of
         -- If the value maps to Nothing, get rid of it.
-        Nothing -> union1 l' r'
-        -- Otherwise, check if the heap invariant still holds
-        -- and sift downwards to restore it.
-        Just !y -> down y l' r'
+        Nothing -> union l' r'
+        -- If y is still the smallest in its subheap,
+        -- the calls to insert and union here will work without making
+        -- any recursive subcalls!
+        Just !y -> insert y l' `union` r'
       where
         !l' = mm l
         !r' = mm r
 
-    down x l@(Node y ll lr) r@(Node z rl rr)
-      -- Put the smallest of x, y and z at the root.
-      | y < x && y <= z =
-        (Node y $! down x ll lr) r
-      | z < x && z <= y =
-        Node z l $! down x rl rr
-    down x Nil (Node y l r)
-      -- Put the smallest of x and y at the root.
-      | y < x =
-        Node y Nil $! down x l r
-    down x (Node y l r) Nil
-      -- Put the smallest of x and y at the root.
-      | y < x =
-        (Node y $! down x l r) Nil
-    down x l r = Node x l r
-
 -- | Return the number of elements in the heap.
 {-# INLINE size #-}
 size :: Heap a -> Int
-size (Heap n _) = n
+size Nil = 0
+size (Node n _ _ _) = n
 
 -- Testing code:
 -- import Test.QuickCheck
 -- import qualified Data.List as List
 -- import qualified Data.Maybe as Maybe
-
+-- 
 -- instance (Arbitrary a, Ord a) => Arbitrary (Heap a) where
 --   arbitrary = sized arb
 --     where
@@ -113,42 +96,45 @@
 --            (n-1, union <$> arb' <*> arb')]
 --         where
 --           arb' = arb (n `div` 2)
-
--- toList :: Ord a => Heap a -> [a]
--- toList = List.unfoldr removeMin
-
+-- 
+-- toSortedList :: Ord a => Heap a -> [a]
+-- toSortedList = List.unfoldr removeMin
+-- 
 -- invariant :: Ord a => Heap a -> Bool
--- invariant h@(Heap n h1) =
---   n == length (toList h) && ord h1
+-- invariant h = ord h && sizeOK h
 --   where
 --     ord Nil = True
---     ord (Node x l r) = ord1 x l && ord1 x r
-
+--     ord (Node _ x l r) = ord1 x l && ord1 x r
+-- 
 --     ord1 _ Nil = True
---     ord1 x h@(Node y _ _) = x <= y && ord h
-
--- prop_1 h = withMaxSuccess 10000 $ invariant h
--- prop_2 x h = withMaxSuccess 10000 $ invariant (insert x h)
+--     ord1 x h@(Node _ y _ _) = x <= y && ord h
+-- 
+--     sizeOK Nil = size Nil == 0
+--     sizeOK (Node s _ l r) =
+--       s == size l + size r + 1
+-- 
+-- prop_1 h = withMaxSuccess 100000 $ invariant h
+-- prop_2 x h = withMaxSuccess 100000 $ invariant (insert x h)
 -- prop_3 h =
---   withMaxSuccess 1000 $
+--   withMaxSuccess 100000 $
 --   case removeMin h of
 --     Nothing -> discard
 --     Just (_, h) -> invariant h
--- prop_4 h = withMaxSuccess 10000 $ List.sort (toList h) == toList h
--- prop_5 x h = withMaxSuccess 10000 $ toList (insert x h) == List.insert x (toList h)
+-- prop_4 h = withMaxSuccess 100000 $ List.sort (toSortedList h) == toSortedList h
+-- prop_5 x h = withMaxSuccess 100000 $ toSortedList (insert x h) == List.insert x (toSortedList h)
 -- prop_6 x h =
---   withMaxSuccess 1000 $
+--   withMaxSuccess 100000 $
 --   case removeMin h of
 --     Nothing -> discard
---     Just (x, h') -> toList h == List.insert x (toList h')
--- prop_7 h1 h2 = withMaxSuccess 10000 $
+--     Just (x, h') -> toSortedList h == List.insert x (toSortedList h')
+-- prop_7 h1 h2 = withMaxSuccess 100000 $
 --   invariant (union h1 h2)
--- prop_8 h1 h2 = withMaxSuccess 10000 $
---   toList (union h1 h2) == List.sort (toList h1 ++ toList h2)
--- prop_9 (Blind f) h = withMaxSuccess 10000 $
+-- prop_8 h1 h2 = withMaxSuccess 100000 $
+--   toSortedList (union h1 h2) == List.sort (toSortedList h1 ++ toSortedList h2)
+-- prop_9 (Blind f) h = withMaxSuccess 100000 $
 --   invariant (mapMaybe f h)
 -- prop_10 (Blind f) h = withMaxSuccess 1000000 $
---   toList (mapMaybe f h) == List.sort (Maybe.mapMaybe f (toList h))
-
+--   toSortedList (mapMaybe f h) == List.sort (Maybe.mapMaybe f (toSortedList h))
+-- 
 -- return []
 -- main = $quickCheckAll
diff --git a/Twee.hs b/Twee.hs
--- a/Twee.hs
+++ b/Twee.hs
@@ -7,7 +7,7 @@
 import qualified Twee.Rule as Rule
 import Twee.Equation
 import qualified Twee.Proof as Proof
-import Twee.Proof(Axiom(..), Proof(..), ProvedGoal(..), provedGoal, certify, derivation, symm)
+import Twee.Proof(Axiom(..), Proof(..), Derivation, ProvedGoal(..), provedGoal, certify, derivation)
 import Twee.CP hiding (Config)
 import qualified Twee.CP as CP
 import Twee.Join hiding (Config, defaultConfig)
@@ -26,14 +26,16 @@
 import Data.Maybe
 import Data.List
 import Data.Function
-import qualified Data.Set as Set
-import Data.Set(Set)
+import qualified Data.Map.Strict as Map
+import Data.Map(Map)
 import Data.Int
 import Data.Ord
 import Control.Monad
 import Control.Monad.IO.Class
 import Control.Monad.Trans.Class
 import qualified Control.Monad.Trans.State.Strict as StateM
+import qualified Data.IntSet as IntSet
+import Data.IntSet(IntSet)
 
 ----------------------------------------------------------------------
 -- * Configuration and prover state.
@@ -47,9 +49,14 @@
     cfg_max_cp_depth           :: Int,
     cfg_simplify               :: Bool,
     cfg_renormalise_percent    :: Int,
+    cfg_cp_sample_size         :: Int,
+    cfg_renormalise_threshold  :: Int,
+    cfg_set_join_goals         :: Bool,
+    cfg_always_simplify        :: Bool,
+    cfg_complete_subsets       :: Bool,
     cfg_critical_pairs         :: CP.Config,
     cfg_join                   :: Join.Config,
-    cfg_proof_presentation     :: Proof.Config }
+    cfg_proof_presentation     :: Proof.Config f }
 
 -- | The prover state.
 data State f =
@@ -63,6 +70,12 @@
     st_next_active    :: {-# UNPACK #-} !Id,
     st_next_rule      :: {-# UNPACK #-} !RuleId,
     st_considered     :: {-# UNPACK #-} !Int64,
+    st_simplified_at  :: {-# UNPACK #-} !Id,
+    st_cp_sample      :: ![Maybe (Overlap f)],
+    st_cp_next_sample :: ![(Integer, Int)],
+    st_num_cps        :: !Integer,
+    st_not_complete   :: !IntSet,
+    st_complete       :: !(Index f (Rule f)),
     st_messages_rev   :: ![Message f] }
 
 -- | The default prover configuration.
@@ -74,6 +87,11 @@
     cfg_max_cp_depth = maxBound,
     cfg_simplify = True,
     cfg_renormalise_percent = 5,
+    cfg_renormalise_threshold = 20,
+    cfg_cp_sample_size = 100,
+    cfg_set_join_goals = True,
+    cfg_always_simplify = False,
+    cfg_complete_subsets = False,
     cfg_critical_pairs = CP.defaultConfig,
     cfg_join = Join.defaultConfig,
     cfg_proof_presentation = Proof.defaultConfig }
@@ -86,8 +104,8 @@
   cfg_max_cp_depth == maxBound
 
 -- | The initial state.
-initialState :: State f
-initialState =
+initialState :: Config f -> State f
+initialState Config{..} =
   State {
     st_rules = RuleIndex.empty,
     st_active_ids = IntMap.empty,
@@ -98,6 +116,12 @@
     st_next_active = 1,
     st_next_rule = 0,
     st_considered = 0,
+    st_simplified_at = 1,
+    st_cp_sample = [],
+    st_cp_next_sample = reservoir cfg_cp_sample_size,
+    st_num_cps = 0,
+    st_not_complete = IntSet.empty,
+    st_complete = Index.empty,
     st_messages_rev = [] }
 
 ----------------------------------------------------------------------
@@ -114,8 +138,12 @@
   | DeleteActive !(Active f)
     -- | The CP queue was simplified.
   | SimplifyQueue
+    -- | All except these axioms are complete (with a suitable-chosen subset of the rules).
+  | NotComplete !IntSet
     -- | The rules were reduced wrt each other.
   | Interreduce
+    -- | Status update: how many queued critical pairs there are.
+  | Status !Int
 
 instance Function f => Pretty (Message f) where
   pPrint (NewActive rule) = pPrint rule
@@ -125,8 +153,16 @@
     text "  (delete rule " <#> pPrint (active_id rule) <#> text ")"
   pPrint SimplifyQueue =
     text "  (simplifying queued critical pairs...)"
+  pPrint (NotComplete ax) =
+    case IntSet.toList ax of
+      [n] ->
+        text "  (axiom" <+> pPrint n <+> "is not completed yet)"
+      xs ->
+        text "  (axioms" <+> text (show xs) <+> "are not completed yet)"
   pPrint Interreduce =
     text "  (simplifying rules with respect to one another...)"
+  pPrint (Status n) =
+    text "  (" <#> pPrint n <+> text "queued critical pairs)"
 
 -- | Emit a message.
 message :: PrettyTerm f => Message f -> State f -> State f
@@ -159,9 +195,9 @@
 
 -- | Compute all critical pairs from a rule.
 {-# INLINEABLE makePassives #-}
+{-# SCC makePassives #-}
 makePassives :: Function f => Config f -> State f -> ActiveRule f -> [Passive Params]
 makePassives Config{..} State{..} rule =
-  {-# SCC makePassive #-}
   [ Passive (fromIntegral (score cfg_critical_pairs o)) (rule_rid rule1) (rule_rid rule2) (fromIntegral (overlap_pos o))
   | (rule1, rule2, o) <- overlaps (Depth cfg_max_cp_depth) (index_oriented st_rules) rules rule ]
   where
@@ -170,8 +206,9 @@
 -- | Turn a Passive back into an overlap.
 -- Doesn't try to simplify it.
 {-# INLINEABLE findPassive #-}
-findPassive :: forall f. Function f => Config f -> State f -> Passive Params -> Maybe (ActiveRule f, ActiveRule f, Overlap f)
-findPassive Config{..} State{..} Passive{..} = {-# SCC findPassive #-} do
+{-# SCC findPassive #-}
+findPassive :: forall f. Function f => State f -> Passive Params -> Maybe (ActiveRule f, ActiveRule f, Overlap f)
+findPassive State{..} Passive{..} = do
   rule1 <- IntMap.lookup (fromIntegral passive_rule1) st_rule_ids
   rule2 <- IntMap.lookup (fromIntegral passive_rule2) st_rule_ids
   let !depth = 1 + max (the rule1) (the rule2)
@@ -182,30 +219,37 @@
 
 -- | Renormalise a queued Passive.
 {-# INLINEABLE simplifyPassive #-}
+{-# SCC simplifyPassive #-}
 simplifyPassive :: Function f => Config f -> State f -> Passive Params -> Maybe (Passive Params)
-simplifyPassive config@Config{..} state@State{..} passive = {-# SCC simplifyPassive #-} do
-  (_, _, overlap) <- findPassive config state passive
+simplifyPassive Config{..} state@State{..} passive = do
+  (_, _, overlap) <- findPassive state passive
   overlap <- simplifyOverlap (index_oriented st_rules) overlap
   return passive {
     passive_score = fromIntegral $
       fromIntegral (passive_score passive) `intMin`
       score cfg_critical_pairs overlap }
 
+-- | Check if we should renormalise the queue.
+{-# INLINEABLE shouldSimplifyQueue #-}
+shouldSimplifyQueue :: Function f => Config f -> State f -> Bool
+shouldSimplifyQueue Config{..} State{..} =
+  length (filter isNothing st_cp_sample) * 100 >= cfg_renormalise_threshold * cfg_cp_sample_size
+
 -- | Renormalise the entire queue.
 {-# INLINEABLE simplifyQueue #-}
+{-# SCC simplifyQueue #-}
 simplifyQueue :: Function f => Config f -> State f -> State f
 simplifyQueue config state =
-  {-# SCC simplifyQueue #-}
-  state { st_queue = simp (st_queue state) }
+  resetSample config state { st_queue = simp (st_queue state) }
   where
     simp =
       Queue.mapMaybe (simplifyPassive config state)
 
 -- | Enqueue a set of critical pairs.
 {-# INLINEABLE enqueue #-}
+{-# SCC enqueue #-}
 enqueue :: Function f => State f -> RuleId -> [Passive Params] -> State f
 enqueue state rule passives =
-  {-# SCC enqueue #-}
   state { st_queue = Queue.insert rule passives (st_queue state) }
 
 -- | Dequeue a critical pair.
@@ -215,9 +259,9 @@
 --   * removing any orphans from the head of the queue
 --   * ignoring CPs that are too big
 {-# INLINEABLE dequeue #-}
+{-# SCC dequeue #-}
 dequeue :: Function f => Config f -> State f -> (Maybe (CriticalPair f, ActiveRule f, ActiveRule f), State f)
-dequeue config@Config{..} state@State{..} =
-  {-# SCC dequeue #-}
+dequeue Config{..} state@State{..} =
   case deq 0 st_queue of
     -- Explicitly make the queue empty, in case it e.g. contained a
     -- lot of orphans
@@ -228,14 +272,11 @@
   where
     deq !n queue = do
       (passive, queue) <- Queue.removeMin queue
-      case findPassive config state passive of
-        Just (rule1, rule2, overlap)
-          | passive_score passive >= 0,
-            Just Overlap{overlap_eqn = t :=: u} <-
-              simplifyOverlap (index_oriented st_rules) overlap,
-            fromMaybe True (cfg_accept_term <*> pure t),
+      case findPassive state passive of
+        Just (rule1, rule2, overlap@Overlap{overlap_eqn = t :=: u})
+          | fromMaybe True (cfg_accept_term <*> pure t),
             fromMaybe True (cfg_accept_term <*> pure u),
-            Just cp <- makeCriticalPair rule1 rule2 overlap ->
+            cp <- makeCriticalPair rule1 rule2 overlap ->
               return ((cp, rule1, rule2), n+1, queue)
         _ -> deq (n+1) queue
 
@@ -250,6 +291,7 @@
     active_rule  :: {-# UNPACK #-} !(Rule f),
     active_top   :: !(Maybe (Term f)),
     active_proof :: {-# UNPACK #-} !(Proof f),
+    active_max   :: !Max,
     -- A model in which the rule is false (used when reorienting)
     active_model :: !(Model f),
     active_rules :: ![ActiveRule f] }
@@ -259,6 +301,7 @@
   CriticalPair {
     cp_eqn = unorient active_rule,
     cp_depth = active_depth,
+    cp_max = active_max,
     cp_top = active_top,
     cp_proof = derivation active_proof }
 
@@ -268,19 +311,17 @@
     rule_active    :: {-# UNPACK #-} !Id,
     rule_rid       :: {-# UNPACK #-} !RuleId,
     rule_depth     :: {-# UNPACK #-} !Depth,
+    rule_max       :: !Max,
     rule_rule      :: {-# UNPACK #-} !(Rule f),
-    rule_proof     :: {-# UNPACK #-} !(Proof f),
     rule_positions :: !(Positions f) }
 
 instance PrettyTerm f => Symbolic (ActiveRule f) where
   type ConstantOf (ActiveRule f) = f
   termsDL ActiveRule{..} =
-    termsDL rule_rule `mplus`
-    termsDL (derivation rule_proof)
+    termsDL rule_rule
   subst_ sub r@ActiveRule{..} =
     r {
       rule_rule = rule',
-      rule_proof = certify (subst_ sub (derivation rule_proof)),
       rule_positions = positions (lhs rule') }
     where
       rule' = subst_ sub rule_rule
@@ -298,30 +339,79 @@
 instance Has (ActiveRule f) Id where the = rule_active
 instance Has (ActiveRule f) RuleId where the = rule_rid
 instance Has (ActiveRule f) Depth where the = rule_depth
+instance Has (ActiveRule f) Max where the = rule_max
 instance f ~ g => Has (ActiveRule f) (Rule g) where the = rule_rule
-instance f ~ g => Has (ActiveRule f) (Proof g) where the = rule_proof
 instance f ~ g => Has (ActiveRule f) (Positions g) where the = rule_positions
 
 newtype RuleId = RuleId Id deriving (Eq, Ord, Show, Num, Real, Integral, Enum)
 
 -- Add a new active.
 {-# INLINEABLE addActive #-}
+{-# SCC addActive #-}
 addActive :: Function f => Config f -> State f -> (Id -> RuleId -> RuleId -> Active f) -> State f
 addActive config state@State{..} active0 =
-  {-# SCC addActive #-}
   let
     active@Active{..} = active0 st_next_active st_next_rule (succ st_next_rule)
     state' =
       message (NewActive active) $
       addActiveOnly state{st_next_active = st_next_active+1, st_next_rule = st_next_rule+2} active
-  in if subsumed st_joinable st_rules (unorient active_rule) then
+  in if subsumed (st_joinable, st_complete) st_rules (unorient active_rule) then
     state
   else
-    normaliseGoals $
-    foldl' (uncurry . enqueue) state'
-      [ (the rule, makePassives config state' rule)
-      | rule <- active_rules ]
+    normaliseGoals config $
+    foldl' enqueueRule state' active_rules
+  where
+    enqueueRule state rule =
+      sample config (length passives) passives $
+      enqueue state (the rule) passives
+      where
+        passives = makePassives config state rule
 
+-- Update the list of sampled critical pairs.
+{-# INLINEABLE sample #-}
+sample :: Function f => Config f -> Int -> [Passive Params] -> State f -> State f
+sample cfg m passives state@State{st_cp_next_sample = ((n, pos):rest), ..}
+  | idx < fromIntegral m =
+    sample cfg m passives state {
+      st_cp_next_sample = rest,
+      st_cp_sample =
+        take pos st_cp_sample ++
+        [find (passives !! fromIntegral idx)] ++
+        drop (pos+1) st_cp_sample }
+  | otherwise = state{st_num_cps = st_num_cps + fromIntegral m}
+  where
+    idx = n - st_num_cps
+    find passive = do
+      (_, _, overlap) <- findPassive state passive
+      simplifyOverlap (index_oriented st_rules) overlap
+
+-- Reset the list of sampled critical pairs.
+{-# INLINEABLE resetSample #-}
+resetSample :: Function f => Config f -> State f -> State f
+resetSample cfg@Config{..} state@State{..} =
+  foldl' sample1 state' (Queue.toList st_queue)
+  where
+    state' =
+      state {
+        st_num_cps = 0,
+        st_cp_next_sample = reservoir cfg_cp_sample_size,
+        st_cp_sample = [] }
+
+    sample1 state (n, passives) = sample cfg n passives state
+
+-- Simplify the sampled critical pairs.
+-- (A sampled critical pair is replaced with Nothing if it can be
+-- simplified.)
+{-# INLINEABLE simplifySample #-}
+simplifySample :: Function f => State f -> State f
+simplifySample state@State{..} =
+  state{st_cp_sample = map (>>= simp) st_cp_sample}
+  where
+    simp overlap = do
+      overlap' <- simplifyOverlap (index_oriented st_rules) overlap
+      guard (overlap_eqn overlap == overlap_eqn overlap')
+      return overlap
+
 -- Add an active without generating critical pairs. Used in interreduction.
 {-# INLINEABLE addActiveOnly #-}
 addActiveOnly :: Function f => State f -> Active f -> State f
@@ -359,15 +449,15 @@
 -- Try to join a critical pair, but using a different set of critical
 -- pairs for normalisation.
 {-# INLINEABLE considerUsing #-}
+{-# SCC considerUsing #-}
 considerUsing ::
   Function f =>
   RuleIndex f (ActiveRule f) -> Config f -> State f -> CriticalPair f -> State f
 considerUsing rules config@Config{..} state@State{..} cp0 =
-  {-# SCC consider #-}
   -- Important to canonicalise the rule so that we don't get
   -- bigger and bigger variable indices over time
   let cp = canonicalise cp0 in
-  case joinCriticalPair cfg_join st_joinable rules Nothing cp of
+  case joinCriticalPair cfg_join (st_joinable, st_complete) rules Nothing cp of
     Right (mcp, cps) ->
       let
         state' = foldl' (considerUsing rules config) state cps
@@ -381,31 +471,32 @@
 {-# INLINEABLE addCP #-}
 addCP :: Function f => Config f -> Model f -> State f -> CriticalPair f -> State f
 addCP config model state@State{..} CriticalPair{..} =
-  addActive config state $ \n k1 k2 ->
   let
     pf = certify cp_proof
-    rule = orient cp_eqn
+    rule = orient cp_eqn pf
 
-    makeRule k r p =
+    makeRule n k r =
       ActiveRule {
         rule_active = n,
         rule_rid = k,
         rule_depth = cp_depth,
+        rule_max = cp_max,
         rule_rule = r rule,
-        rule_proof = p pf,
         rule_positions = positions (lhs (r rule)) }
   in
+  addActive config state $ \n k1 k2 ->
   Active {
     active_id = n,
     active_depth = cp_depth,
     active_rule = rule,
     active_model = model,
     active_top = cp_top,
+    active_max = cp_max,
     active_proof = pf,
     active_rules =
       usortBy (comparing (canonicalise . rule_rule)) $
-        makeRule k1 id id:
-        [ makeRule k2 backwards (certify . symm . derivation)
+        makeRule n k1 id:
+        [ makeRule n k2 backwards
         | not (oriented (orientation rule)) ] }
 
 -- Add a new equation.
@@ -416,6 +507,7 @@
     CriticalPair {
       cp_eqn = axiom_eqn axiom,
       cp_depth = 0,
+      cp_max = Max $ IntSet.fromList [axiom_number axiom | cfg_complete_subsets config],
       cp_top = Nothing,
       cp_proof = Proof.axiom axiom }
 
@@ -429,55 +521,122 @@
       Index.insert t (t :=: u) $
       Index.insert u (u :=: t) (st_joinable state) }
 
+-- Find a confluent subset of the rules.
+{-# INLINEABLE checkCompleteness #-}
+checkCompleteness :: Function f => Config f -> State f -> State f
+checkCompleteness _ state@State{..} | st_simplified_at == st_next_active = state
+checkCompleteness _config state =
+  state { st_not_complete = excluded,
+          st_complete = Index.fromListWith lhs rules }
+  where
+    maxSet s = if IntSet.null s then minBound else IntSet.findMax s
+    maxN = maximum [maxSet (unMax (active_max r)) | r <- IntMap.elems (st_active_ids state)]
+    excluded = go IntSet.empty
+    go excl
+      | m > maxN = excl
+      | otherwise = go (IntSet.insert m excl)
+      where
+        m = bound excl
+
+    bound excl = minimum . map (passiveMax excl) . concatMap snd . Queue.toList $ st_queue state
+
+    passiveMax excl p = fromMaybe maxBound $ do
+      (r1, r2, _) <- findPassive state p
+      let s = unMax (rule_max r1) `IntSet.union` unMax (rule_max r2)
+      guard (s `IntSet.disjoint` excl)
+      (n, _) <- IntSet.maxView s
+      return n
+    rules = map rule_rule (filter ok (IntMap.elems (st_rule_ids state)))
+    ok r = unMax (rule_max r) `IntSet.disjoint` excluded
+
 -- For goal terms we store the set of all their normal forms.
 -- Name and number are for information only.
 data Goal f =
   Goal {
-    goal_name     :: String,
-    goal_number   :: Int,
-    goal_eqn      :: Equation f,
-    goal_lhs      :: Set (Resulting f),
-    goal_rhs      :: Set (Resulting f) }
+    goal_name         :: String,
+    goal_number       :: Int,
+    goal_eqn          :: Equation f,
+    goal_expanded_lhs :: Map (Term f) (Derivation f),
+    goal_expanded_rhs :: Map (Term f) (Derivation f),
+    goal_lhs          :: Map (Term f) (Term f, Reduction f),
+    goal_rhs          :: Map (Term f) (Term f, Reduction f) }
+  deriving Show
 
 -- Add a new goal.
 {-# INLINEABLE addGoal #-}
 addGoal :: Function f => Config f -> State f -> Goal f -> State f
-addGoal _config state@State{..} goal =
-  normaliseGoals state { st_goals = goal:st_goals }
+addGoal config state@State{..} goal =
+  normaliseGoals config state { st_goals = goal:st_goals }
 
 -- Normalise all goals.
 {-# INLINEABLE normaliseGoals #-}
-normaliseGoals :: Function f => State f -> State f
-normaliseGoals state@State{..} =
+normaliseGoals :: Function f => Config f -> State f -> State f
+normaliseGoals Config{..} state@State{..} =
   state {
     st_goals =
-      map (goalMap (Rule.normalForms (rewrite reduces (index_all st_rules)))) st_goals }
+      map (goalMap (nf (rewrite reduces (index_all st_rules)))) st_goals }
   where
     goalMap f goal@Goal{..} =
       goal { goal_lhs = f goal_lhs, goal_rhs = f goal_rhs }
+    nf reduce goals
+      | cfg_set_join_goals =
+        let pair (t, red) = (fst (goals Map.! t), red) in
+        Map.map pair $ Rule.normalForms reduce (Map.map snd goals)
+      | otherwise =
+        Map.fromList $
+          [ (result t q, (u, r `trans` q))
+          | (t, (u, r)) <- Map.toList goals,
+            let q = Rule.normaliseWith (const True) reduce t ]
 
 -- Recompute all normal forms of all goals. Starts from the original goal term.
 -- Different from normalising all goals, because there may be an intermediate
 -- term on one of the reduction paths which we can now rewrite in a different
 -- way.
 {-# INLINEABLE recomputeGoals #-}
-recomputeGoals :: Function f => State f -> State f
-recomputeGoals state =
+recomputeGoals :: Function f => Config f -> State f -> State f
+recomputeGoals config state =
   -- Make this strict so that newTask can time it correctly
   forceList (map goal_lhs (st_goals state')) `seq`
   forceList (map goal_rhs (st_goals state')) `seq`
   state'
   where
     state' =
-      normaliseGoals (state { st_goals = map reset (st_goals state) })
-
-    reset goal@Goal{goal_eqn = t :=: u, ..} =
-      goal { goal_lhs = Set.singleton (reduce (Refl t)),
-             goal_rhs = Set.singleton (reduce (Refl u)) }
+      normaliseGoals config (state { st_goals = map resetGoal (st_goals state) })
 
     forceList [] = ()
     forceList (x:xs) = x `seq` forceList xs
 
+resetGoal :: Goal f -> Goal f
+resetGoal goal@Goal{..} =
+  goal { goal_lhs = expansions goal_expanded_lhs,
+         goal_rhs = expansions goal_expanded_rhs }
+  where
+    expansions m =
+      Map.mapWithKey (\t _ -> (t, [])) m
+
+-- Rewrite goal terms backwards using rewrite rules.
+{-# INLINEABLE rewriteGoalsBackwards #-}
+rewriteGoalsBackwards :: Function f => State f -> State f
+rewriteGoalsBackwards state =
+  state { st_goals = map backwardsGoal (st_goals state) }
+  where
+    backwardsGoal goal@Goal{..} =
+      resetGoal goal {
+        goal_expanded_lhs = backwardsMap goal_expanded_lhs,
+        goal_expanded_rhs = backwardsMap goal_expanded_rhs }
+    backwardsMap m =
+      Map.fromList $
+        Map.toList m ++
+        [ (ruleResult t r, p `Proof.trans` q)
+        | (t, p) <- Map.toList m,
+          r <- backwardsTerm t,
+          let q = ruleProof t r ]
+    backwardsTerm t = do
+      rule <- map the (Index.elems (RuleIndex.index_all (st_rules state)))
+      guard (usort (vars (lhs rule)) == usort (vars (rhs rule)))
+      [r] <- anywhere (tryRule (\_ _ -> True) (backwards rule)) t
+      return r
+
 -- Create a goal.
 {-# INLINE goal #-}
 goal :: Int -> String -> Equation f -> Goal f
@@ -486,8 +645,10 @@
     goal_name = name,
     goal_number = n,
     goal_eqn = t :=: u,
-    goal_lhs = Set.singleton (reduce (Refl t)),
-    goal_rhs = Set.singleton (reduce (Refl u)) }
+    goal_expanded_lhs = Map.singleton t (Proof.Refl t),
+    goal_expanded_rhs = Map.singleton u (Proof.Refl u),
+    goal_lhs = Map.singleton t (t, []),
+    goal_rhs = Map.singleton u (u, []) }
 
 ----------------------------------------------------------------------
 -- Interreduction.
@@ -495,17 +656,18 @@
 
 -- Simplify all rules.
 {-# INLINEABLE interreduce #-}
+{-# SCC interreduce #-}
 interreduce :: Function f => Config f -> State f -> State f
+interreduce _ state@State{..} | st_simplified_at == st_next_active = state
 interreduce config@Config{..} state =
-  {-# SCC interreduce #-}
   let
     state' =
       foldl' (interreduce1 config)
         -- Clear out st_joinable, since we don't know which
         -- equations have made use of each active.
-        state { st_joinable = Index.empty }
+        state { st_joinable = Index.empty, st_complete = Index.empty }
         (IntMap.elems (st_active_ids state))
-    in state' { st_joinable = st_joinable state }
+    in state' { st_joinable = st_joinable state, st_complete = st_complete state, st_simplified_at = st_next_active state' }
 
 {-# INLINEABLE interreduce1 #-}
 interreduce1 :: Function f => Config f -> State f -> Active f -> State f
@@ -514,7 +676,7 @@
   -- joinability, otherwise it will be trivially joinable.
   case
     joinCriticalPair cfg_join
-      (st_joinable state)
+      (Index.empty, Index.empty) -- (st_joinable state)
       (st_rules (deleteActive state active))
       (Just (active_model active)) (active_cp active)
   of
@@ -523,8 +685,8 @@
       message (DeleteActive active) $
       deleteActive state active
     Left (cp, model)
-      | not (cp_eqn cp `isInstanceOf` cp_eqn (active_cp active)) ->
-        flip (foldl' (addCP config model)) (split cp) $
+      | cp_eqn cp `simplerThan` cp_eqn (active_cp active) ->
+        flip (foldl' (consider config)) (split cp) $
         message (DeleteActive active) $
         deleteActive state active
       | model /= active_model active ->
@@ -532,10 +694,6 @@
         deleteActive state active
       | otherwise ->
         state
-  where
-    (t :=: u) `isInstanceOf` (t' :=: u') = isJust $ do
-      sub <- match t' t
-      matchIn sub u' u
 
 ----------------------------------------------------------------------
 -- The main loop.
@@ -550,22 +708,37 @@
 complete Output{..} config@Config{..} state =
   flip StateM.execStateT state $ do
     tasks <- sequence
-      [newTask 1 (fromIntegral cfg_renormalise_percent / 100) $ do
-         lift $ output_message SimplifyQueue
+      [newTask 10 (fromIntegral cfg_renormalise_percent / 100) $ do
          state <- StateM.get
-         StateM.put $! simplifyQueue config state,
+         when (shouldSimplifyQueue config state) $ do
+           lift $ output_message SimplifyQueue
+           StateM.put $! simplifyQueue config state,
+       newTask 1 0.02 $ do
+         when cfg_complete_subsets $ do
+           state <- StateM.get
+           let !state' = checkCompleteness config state
+           lift $ output_message (NotComplete (st_not_complete state'))
+           StateM.put $! state',
        newTask 1 0.05 $ do
          when cfg_simplify $ do
            lift $ output_message Interreduce
            state <- StateM.get
-           StateM.put $! interreduce config state,
+           StateM.put $! simplifySample $! interreduce config state,
        newTask 1 0.02 $ do
           state <- StateM.get
-          StateM.put $! recomputeGoals state ]
+          StateM.put $! recomputeGoals config state,
+       newTask 60 0.01 $ do
+          State{..} <- StateM.get
+          let !n = Queue.queueSize st_queue
+          lift $ output_message (Status n)]
 
     let
       loop = do
         progress <- StateM.state (complete1 config)
+        when cfg_always_simplify $ do
+          lift $ output_message Interreduce
+          state <- StateM.get
+          StateM.put $! simplifySample $! interreduce config state
         state <- StateM.get
         lift $ mapM_ output_message (messages state)
         StateM.put (clearMessages state)
@@ -592,18 +765,24 @@
 
 -- Return whatever goals we have proved and their proofs.
 {-# INLINEABLE solutions #-}
+{-# SCC solutions #-}
 solutions :: Function f => State f -> [ProvedGoal f]
-solutions State{..} = {-# SCC solutions #-} do
+solutions State{..} = do
   Goal{goal_lhs = ts, goal_rhs = us, ..} <- st_goals
-  guard (not (null (Set.intersection ts us)))
-  let t:_ = filter (`Set.member` us) (Set.toList ts)
-      u:_ = filter (== t) (Set.toList us)
+  let sols = Map.keys (Map.intersection ts us)
+  guard (not (null sols))
+  let sol:_ = sols
+  let t = ts Map.! sol
+      u = us Map.! sol
       -- Strict so that we check the proof before returning a solution
       !p =
         Proof.certify $
-          reductionProof (reduction t) `Proof.trans`
-          Proof.symm (reductionProof (reduction u))
+          expandedProof goal_expanded_lhs t `Proof.trans`
+          Proof.symm (expandedProof goal_expanded_rhs u)
   return (provedGoal goal_number goal_name p)
+  where
+    expandedProof m (t, red) =
+      m Map.! t `Proof.trans` reductionProof t red
 
 -- Return all current rewrite rules.
 {-# INLINEABLE rules #-}
@@ -623,14 +802,15 @@
     (progress, state') = complete1 cfg state
 
 {-# INLINEABLE normaliseTerm #-}
-normaliseTerm :: Function f => State f -> Term f -> Resulting f
+normaliseTerm :: Function f => State f -> Term f -> Reduction f
 normaliseTerm State{..} t =
   normaliseWith (const True) (rewrite reduces (index_all st_rules)) t
 
 {-# INLINEABLE normalForms #-}
-normalForms :: Function f => State f -> Term f -> Set (Resulting f)
+normalForms :: Function f => State f -> Term f -> Map (Term f) (Reduction f)
 normalForms State{..} t =
-  Rule.normalForms (rewrite reduces (index_all st_rules)) (Set.singleton (reduce (Refl t)))
+  Map.map snd $
+  Rule.normalForms (rewrite reduces (index_all st_rules)) (Map.singleton t [])
 
 {-# INLINEABLE simplifyTerm #-}
 simplifyTerm :: Function f => State f -> Term f -> Term f
diff --git a/Twee/Base.hs b/Twee/Base.hs
--- a/Twee/Base.hs
+++ b/Twee/Base.hs
@@ -8,23 +8,24 @@
   -- * The 'Symbolic' typeclass
   Symbolic(..), subst, terms,
   TermOf, TermListOf, SubstOf, TriangleSubstOf, BuilderOf, FunOf,
-  vars, isGround, funs, occ, occVar, canonicalise, renameAvoiding,
+  vars, isGround, funs, occ, occVar, canonicalise, renameAvoiding, renameManyAvoiding, freshVar,
   -- * General-purpose functionality
   Id(..), Has(..),
   -- * Typeclasses
-  Minimal(..), minimalTerm, isMinimal, erase,
-  Skolem(..), Arity(..), Sized(..), Ordered(..), lessThan, orientTerms, EqualsBonus(..), Strictness(..), Function, Extended(..)) where
+  Minimal(..), minimalTerm, isMinimal, erase, eraseExcept, ground,
+  Arity(..), Ordered(..), lessThan, orientTerms, EqualsBonus(..), Strictness(..), Function) where
 
 import Prelude hiding (lookup)
 import Control.Monad
 import qualified Data.DList as DList
 import Twee.Term hiding (subst, canonicalise)
 import qualified Twee.Term as Term
+import Twee.Utils
 import Twee.Pretty
 import Twee.Constraints hiding (funs)
 import Data.DList(DList)
-import Data.Typeable
 import Data.Int
+import Data.List
 import Data.Maybe
 import qualified Data.IntMap.Strict as IntMap
 
@@ -172,6 +173,22 @@
     (V x1, V x2) = boundLists (terms x)
     (V y1, V y2) = boundLists (terms y)
 
+-- | Return an x such that no variable >= x occurs in the argument.
+freshVar :: Symbolic a => a -> Int
+freshVar x
+  | x1 > x2 = 0 -- x is ground
+  | otherwise = x2+1
+  where
+    (V x1, V x2) = boundLists (terms x)
+
+{-# INLINEABLE renameManyAvoiding #-}
+renameManyAvoiding :: Symbolic a => [a] -> [a]
+renameManyAvoiding [] = []
+renameManyAvoiding (t:ts) = u:us
+  where
+    u = renameAvoiding us t
+    us = renameManyAvoiding ts
+  
 -- | Check if a term is the minimal constant.
 isMinimal :: Minimal f => Term f -> Bool
 isMinimal (App f Empty) | f == minimal = True
@@ -189,41 +206,28 @@
   where
     sub = fromMaybe undefined $ listToSubst [(x, minimalTerm) | x <- xs]
 
--- | Construction of Skolem constants.
-class Skolem f where
-  -- | Turn a variable into a Skolem constant.
-  skolem  :: Var -> Fun f
-  getSkolem :: Fun f -> Maybe Var
+-- | Erase all except a given set of variables from the argument, replacing them
+-- with the minimal constant.
+eraseExcept :: (Symbolic a, ConstantOf a ~ f, Minimal f) => [Var] -> a -> a
+eraseExcept xs t =
+  erase (usort (vars t) \\ xs) t
 
+-- | Replace all variables in the argument with the minimal constant.
+ground :: (Symbolic a, ConstantOf a ~ f, Minimal f) => a -> a
+ground t = erase (vars t) t
+
 -- | For types which have a notion of arity.
 class Arity f where
   -- | Measure the arity.
   arity :: f -> Int
 
-instance Arity f => Arity (Fun f) where
+instance (Labelled f, Arity f) => Arity (Fun f) where
   arity = arity . fun_value
 
 -- | For types which have a notion of size.
-class Sized a where
-  -- | Compute the size.
-  size  :: a -> Int
-
-instance Sized f => Sized (Fun f) where
-  size = size . fun_value
-
-instance Sized f => Sized (TermList f) where
-  size = aux 0
-    where
-      aux n Empty = n
-      aux n (ConsSym (App f _) t) = aux (n+size f) t
-      aux n (Cons (Var _) t) = aux (n+1) t
-
-instance Sized f => Sized (Term f) where
-  size = size . singleton
-
 -- | The collection of constraints which the type of function symbols must
 -- satisfy in order to be used by twee.
-type Function f = (Ordered f, Arity f, Sized f, Minimal f, Skolem f, PrettyTerm f, EqualsBonus f)
+type Function f = (Ordered f, Arity f, Minimal f, PrettyTerm f, EqualsBonus f, Labelled f)
 
 -- | A hack for encoding Horn clauses. See 'Twee.CP.Score'.
 -- The default implementation of 'hasEqualsBonus' should work OK.
@@ -235,54 +239,8 @@
   isTrue _ = False
   isFalse _ = False
 
-instance EqualsBonus f => EqualsBonus (Fun f) where
+instance (Labelled f, EqualsBonus f) => EqualsBonus (Fun f) where
   hasEqualsBonus = hasEqualsBonus . fun_value
   isEquals = isEquals . fun_value
   isTrue = isTrue . fun_value
   isFalse = isFalse . fun_value
-
--- | A function symbol extended with a minimal constant and Skolem functions.
--- Comes equipped with 'Minimal' and 'Skolem' instances.
-data Extended f =
-    -- | The minimal constant.
-    Minimal
-    -- | A Skolem function.
-  | Skolem Var
-    -- | An ordinary function symbol.
-  | Function f
-  deriving (Eq, Ord, Show, Functor)
-
-instance Pretty f => Pretty (Extended f) where
-  pPrintPrec _ _ Minimal = text "?"
-  pPrintPrec _ _ (Skolem (V n)) = text "sk" <#> pPrint n
-  pPrintPrec l p (Function f) = pPrintPrec l p f
-
-instance PrettyTerm f => PrettyTerm (Extended f) where
-  termStyle (Function f) = termStyle f
-  termStyle _ = uncurried
-
-instance Sized f => Sized (Extended f) where
-  size (Function f) = size f
-  size _ = 1
-
-instance Arity f => Arity (Extended f) where
-  arity (Function f) = arity f
-  arity _ = 0
-
-instance (Typeable f, Ord f) => Minimal (Extended f) where
-  minimal = fun Minimal
-
-instance (Typeable f, Ord f) => Skolem (Extended f) where
-  skolem x = fun (Skolem x)
-  getSkolem (F (Skolem x)) = Just x
-  getSkolem _ = Nothing
-
-instance EqualsBonus f => EqualsBonus (Extended f) where
-  hasEqualsBonus (Function f) = hasEqualsBonus f
-  hasEqualsBonus _ = False
-  isEquals (Function f) = isEquals f
-  isEquals _ = False
-  isTrue (Function f) = isTrue f
-  isTrue _ = False
-  isFalse (Function f) = isFalse f
-  isFalse _ = False
diff --git a/Twee/CP.hs b/Twee/CP.hs
--- a/Twee/CP.hs
+++ b/Twee/CP.hs
@@ -14,8 +14,13 @@
 import Twee.Utils
 import Twee.Equation
 import qualified Twee.Proof as Proof
-import Twee.Proof(Derivation, Proof, congPath)
+import Twee.Proof(Derivation, congPath)
+import Data.IntSet(IntSet)
+import qualified Data.IntSet as IntSet
 
+newtype Max = Max { unMax :: IntSet }
+  deriving (Eq, Ord, Show)
+
 -- | The set of positions at which a term can have critical overlaps.
 data Positions f = NilP | ConsP {-# UNPACK #-} !Int !(Positions f)
 type PositionsOf a = Positions (ConstantOf a)
@@ -30,7 +35,7 @@
     -- Consider only general superpositions.
     aux !_ !_ Empty = NilP
     aux n m (Cons (Var _) t) = aux (n+1) m t
-    aux n m (ConsSym t@App{} u)
+    aux n m ConsSym{hd = t@App{}, rest = u}
       | t `Set.member` m = aux (n+1) m u
       | otherwise = ConsP n (aux (n+1) (Set.insert t m) u)
 
@@ -66,7 +71,7 @@
 -- | Compute all overlaps of a rule with a set of rules.
 {-# INLINEABLE overlaps #-}
 overlaps ::
-  (Function f, Has a (Rule f), Has a (Positions f), Has a Depth) =>
+  forall a f. (Function f, Has a Id, Has a (Rule f), Has a (Positions f), Has a Depth) =>
   Depth -> Index f a -> [a] -> a -> [(a, a, Overlap f)]
 overlaps max_depth idx rules r =
   ChurchList.toList (overlapsChurch max_depth idx rules r)
@@ -98,17 +103,17 @@
 -- | Create an overlap at a particular position in a term.
 -- Doesn't simplify the overlap.
 {-# INLINE overlapAt #-}
+{-# SCC overlapAt #-}
 overlapAt :: Int -> Depth -> Rule f -> Rule f -> Maybe (Overlap f)
-overlapAt !n !depth (Rule _ !outer !outer') (Rule _ !inner !inner') = do
+overlapAt !n !depth (Rule _ _ !outer !outer') (Rule _ _ !inner !inner') = do
   let t = at n (singleton outer)
   sub <- unifyTri inner t
   let
-    top = {-# SCC overlap_top #-} termSubst sub outer
-    innerTerm = {-# SCC overlap_inner #-} termSubst sub inner
+    top = termSubst sub outer
+    innerTerm = termSubst sub inner
     -- Make sure to keep in sync with overlapProof
-    lhs = {-# SCC overlap_eqn_1 #-} termSubst sub outer'
-    rhs = {-# SCC overlap_eqn_2 #-}
-      buildReplacePositionSub sub n (singleton inner') (singleton outer)
+    lhs = termSubst sub outer'
+    rhs = buildReplacePositionSub sub n (singleton inner') (singleton outer)
 
   guard (lhs /= rhs)
   return Overlap {
@@ -122,7 +127,7 @@
 {-# INLINE simplifyOverlap #-}
 simplifyOverlap :: (Function f, Has a (Rule f)) => Index f a -> Overlap f -> Maybe (Overlap f)
 simplifyOverlap idx overlap@Overlap{overlap_eqn = lhs :=: rhs, ..}
-  | lhs == rhs'  = Nothing
+  | lhs == rhs   = Nothing
   | lhs' == rhs' = Nothing
   | otherwise = Just overlap{overlap_eqn = lhs' :=: rhs'}
   where
@@ -152,7 +157,7 @@
 defaultConfig :: Config
 defaultConfig =
   Config {
-    cfg_lhsweight = 3,
+    cfg_lhsweight = 4,
     cfg_rhsweight = 1,
     cfg_funweight = 7,
     cfg_varweight = 6,
@@ -180,19 +185,19 @@
     size' !n Empty = n
     size' n (Cons t ts)
       | len t > 1, t `isSubtermOfList` ts =
-        size' (n+cfg_dupcost+cfg_dupfactor*size t) ts
+        size' (n+cfg_dupcost+cfg_dupfactor*len t) ts
     size' n ts
       | Cons (App f ws@(Cons a (Cons b us))) vs <- ts,
-        hasEqualsBonus (fun_value f),
         not (isVar a),
         not (isVar b),
+        hasEqualsBonus (fun_value f),
         Just sub <- unify a b =
-        size' (n+cfg_funweight*size f) ws `min`
+        size' (n+cfg_funweight) ws `min`
         size' (size' (n+1) (subst sub us)) (subst sub vs)
     size' n (Cons (Var _) ts) =
       size' (n+cfg_varweight) ts
-    size' n (ConsSym (App f _) ts) =
-      size' (n+cfg_funweight*size f) ts
+    size' n ConsSym{hd = App{}, rest = ts} =
+      size' (n+cfg_funweight) ts
 
 ----------------------------------------------------------------------
 -- * Higher-level handling of critical pairs.
@@ -205,6 +210,7 @@
     cp_eqn   :: {-# UNPACK #-} !(Equation f),
     -- | The depth of the critical pair.
     cp_depth :: {-# UNPACK #-} !Depth,
+    cp_max :: !Max,
     -- | The critical term, if there is one.
     -- (Axioms do not have a critical term.)
     cp_top   :: !(Maybe (Term f)),
@@ -219,6 +225,7 @@
     CriticalPair {
       cp_eqn = subst_ sub cp_eqn,
       cp_depth = cp_depth,
+      cp_max = cp_max,
       cp_top = subst_ sub cp_top,
       cp_proof = subst_ sub cp_proof }
 
@@ -258,18 +265,21 @@
     [ CriticalPair {
         cp_eqn   = l :=: r',
         cp_depth = cp_depth,
+        cp_max   = cp_max,
         cp_top   = eraseExcept (vars l) cp_top,
         cp_proof = eraseExcept (vars l) cp_proof }
     | ord == Just GT ] ++
     [ CriticalPair {
         cp_eqn   = r :=: l',
         cp_depth = cp_depth,
+        cp_max   = cp_max,
         cp_top   = eraseExcept (vars r) cp_top,
         cp_proof = Proof.symm (eraseExcept (vars r) cp_proof) }
     | ord == Just LT ] ++
     [ CriticalPair {
         cp_eqn   = l' :=: r',
         cp_depth = cp_depth,
+        cp_max   = cp_max,
         cp_top   = eraseExcept (vars l) $ eraseExcept (vars r) cp_top,
         cp_proof = eraseExcept (vars l) $ eraseExcept (vars r) cp_proof }
     | ord == Nothing ] ++
@@ -278,12 +288,14 @@
     [ CriticalPair {
         cp_eqn   = l :=: l',
         cp_depth = cp_depth + 1,
+        cp_max   = cp_max,
         cp_top   = Nothing,
         cp_proof = cp_proof `Proof.trans` Proof.symm (erase ls cp_proof) }
     | not (null ls), ord /= Just GT ] ++
     [ CriticalPair {
         cp_eqn   = r :=: r',
         cp_depth = cp_depth + 1,
+        cp_max   = cp_max,
         cp_top   = Nothing,
         cp_proof = Proof.symm cp_proof `Proof.trans` erase rs cp_proof }
     | not (null rs), ord /= Just LT ]
@@ -300,30 +312,25 @@
 -- | Make a critical pair from two rules and an overlap.
 {-# INLINEABLE makeCriticalPair #-}
 makeCriticalPair ::
-  (Has a (Rule f), Has a (Proof f), Has a Id, Function f) =>
-  a -> a -> Overlap f -> Maybe (CriticalPair f)
-makeCriticalPair r1 r2 overlap@Overlap{..}
-  | lessEq overlap_top t = Nothing
-  | lessEq overlap_top u = Nothing
-  | otherwise =
-    Just $
-      CriticalPair overlap_eqn
-        overlap_depth
-        (Just overlap_top)
-        (overlapProof r1 r2 overlap)
-  where
-    t :=: u = overlap_eqn
+  forall f a. (Has a (Rule f), Has a Id, Has a Max, Function f) =>
+  a -> a -> Overlap f -> CriticalPair f
+makeCriticalPair r1 r2 overlap@Overlap{..} =
+  CriticalPair overlap_eqn
+    overlap_depth
+    (Max (unMax (the r1) `IntSet.union` unMax (the r2)))
+    (Just overlap_top)
+    (overlapProof r1 r2 overlap)
 
 -- | Return a proof for a critical pair.
 {-# INLINEABLE overlapProof #-}
 overlapProof ::
   forall a f.
-  (Has a (Rule f), Has a (Proof f), Has a Id) =>
+  (Has a (Rule f), Has a Id) =>
   a -> a -> Overlap f -> Derivation f
 overlapProof left right Overlap{..} =
-  Proof.symm (reductionProof (step left leftSub))
+  Proof.symm (ruleDerivation (subst leftSub (the left)))
   `Proof.trans`
-  congPath path overlap_top (reductionProof (step right rightSub))
+  congPath path overlap_top (ruleDerivation (subst rightSub (the right)))
   where
     Just leftSub = match (lhs (the left)) overlap_top
     Just rightSub = match (lhs (the right)) overlap_inner
diff --git a/Twee/Constraints.hs b/Twee/Constraints.hs
--- a/Twee/Constraints.hs
+++ b/Twee/Constraints.hs
@@ -24,7 +24,7 @@
     aux Empty = []
     aux (Cons (App f Empty) t) = Constant f:aux t
     aux (Cons (Var x) t) = Variable x:aux t
-    aux (ConsSym _ t) = aux t
+    aux ConsSym{rest = t} = aux t
 
 toTerm :: Atom f -> Term f
 toTerm (Constant f) = build (con f)
@@ -115,7 +115,7 @@
 norm :: Eq f => Branch f -> Atom f -> Atom f
 norm Branch{..} x = fromMaybe x (lookup x equals)
 
-contradictory :: (Minimal f, Ord f) => Branch f -> Bool
+contradictory :: (Minimal f, Ord f, Labelled f) => Branch f -> Bool
 contradictory Branch{..} =
   or [f == minimal | (_, Constant f) <- less] ||
   or [f /= g | (Constant f, Constant g) <- equals] ||
@@ -125,7 +125,7 @@
     cyclic (AcyclicSCC _) = False
     cyclic (CyclicSCC _) = True
 
-formAnd :: (Minimal f, Ordered f) => Formula f -> [Branch f] -> [Branch f]
+formAnd :: (Minimal f, Ordered f, Labelled f) => Formula f -> [Branch f] -> [Branch f]
 formAnd f bs = usort (bs >>= add f)
   where
     add (Less t u) b = addLess t u b
@@ -134,7 +134,7 @@
     add (And (f:fs)) b = add f b >>= add (And fs)
     add (Or fs) b = usort (concat [ add f b | f <- fs ])
 
-branches :: (Minimal f, Ordered f) => Formula f -> [Branch f]
+branches :: (Minimal f, Ordered f, Labelled f) => Formula f -> [Branch f]
 branches x = aux [x]
   where
     aux [] = [Branch [] [] []]
@@ -146,7 +146,7 @@
       concatMap (addLess t u) (aux xs) ++
       concatMap (addEquals u t) (aux xs)
 
-addLess :: (Minimal f, Ordered f) => Atom f -> Atom f -> Branch f -> [Branch f]
+addLess :: (Minimal f, Ordered f, Labelled f) => Atom f -> Atom f -> Branch f -> [Branch f]
 addLess _ (Constant min) _ | min == minimal = []
 addLess (Constant min) _ b | min == minimal = [b]
 addLess t0 u0 b@Branch{..} =
@@ -156,7 +156,7 @@
     t = norm b t0
     u = norm b u0
 
-addEquals :: (Minimal f, Ordered f) => Atom f -> Atom f -> Branch f -> [Branch f]
+addEquals :: (Minimal f, Ordered f, Labelled f) => Atom f -> Atom f -> Branch f -> [Branch f]
 addEquals t0 u0 b@Branch{..}
   | t == u || (t, u) `elem` equals = [b]
   | otherwise =
@@ -172,7 +172,7 @@
       | x == t = u
       | otherwise = x
 
-addTerm :: (Minimal f, Ordered f) => Atom f -> Branch f -> Branch f
+addTerm :: (Minimal f, Ordered f, Labelled f) => Atom f -> Branch f -> Branch f
 addTerm (Constant f) b
   | f `notElem` funs b =
     b {
@@ -251,7 +251,7 @@
   minimal :: Fun f
 
 {-# INLINE lessEqInModel #-}
-lessEqInModel :: (Minimal f, Ordered f) => Model f -> Atom f -> Atom f -> Maybe Strictness
+lessEqInModel :: (Minimal f, Ordered f, Labelled f) => Model f -> Atom f -> Atom f -> Maybe Strictness
 lessEqInModel (Model m) x y
   | Just (a, _) <- Map.lookup x m,
     Just (b, _) <- Map.lookup y m,
@@ -264,7 +264,7 @@
   | Constant a <- x, a == minimal = Just Nonstrict
   | otherwise = Nothing
 
-solve :: (Minimal f, Ordered f, PrettyTerm f) => [Atom f] -> Branch f -> Either (Model f) (Subst f)
+solve :: (Minimal f, Ordered f, PrettyTerm f, Labelled f) => [Atom f] -> Branch f -> Either (Model f) (Subst f)
 solve xs branch@Branch{..}
   | null equals && not (all true less) =
     error $ "Model " ++ prettyShow model ++ " is not a model of " ++ prettyShow branch ++ " (edges = " ++ prettyShow edges ++ ", vs = " ++ prettyShow vs ++ ")"
@@ -288,6 +288,7 @@
   -- | Check if the first term is less than or equal to the second in the given model,
   -- and decide whether the inequality is strict or nonstrict.
   lessIn :: Model f -> Term f -> Term f -> Maybe Strictness
+  lessEqSkolem :: Term f -> Term f -> Bool
 
 -- | Describes whether an inequality is strict or nonstrict.
 data Strictness =
diff --git a/Twee/Equation.hs b/Twee/Equation.hs
--- a/Twee/Equation.hs
+++ b/Twee/Equation.hs
@@ -3,7 +3,6 @@
 module Twee.Equation where
 
 import Twee.Base
-import Data.Maybe
 import Control.Monad
 
 --------------------------------------------------------------------------------
@@ -25,19 +24,13 @@
 instance PrettyTerm f => Pretty (Equation f) where
   pPrint (x :=: y) = pPrint x <+> text "=" <+> pPrint y
 
-instance Sized f => Sized (Equation f) where
-  size (x :=: y) = size x + size y
-
 -- | Order an equation roughly left-to-right.
 -- However, there is no guarantee that the result is oriented.
 order :: Function f => Equation f -> Equation f
 order (l :=: r)
   | l == r = l :=: r
-  | otherwise =
-    case compare (size l) (size r) of
-      LT -> r :=: l
-      GT -> l :=: r
-      EQ -> if lessEq l r then r :=: l else l :=: r
+  | lessEqSkolem l r = r :=: l
+  | otherwise = l :=: r
 
 -- | Apply a function to both sides of an equation.
 bothSides :: (Term f -> Term f') -> Equation f -> Equation f'
@@ -47,12 +40,17 @@
 trivial :: Eq f => Equation f -> Bool
 trivial (t :=: u) = t == u
 
+-- | A total order on equations. Equations with lesser terms are smaller.
 simplerThan :: Function f => Equation f -> Equation f -> Bool
 eq1 `simplerThan` eq2 =
-  t1 `lessEq` t2 &&
-  (isNothing (unify t1 t2) || (u1 `lessEq` u2))
+  --traceShow (hang (pPrint eq1) 2 (text "`simplerThan`" <+> pPrint eq2 <+> text "=" <+> pPrint res)) res
+  t1 `lessEqSkolem` t2 && (t1 /= t2 || ((u1 `lessEqSkolem` u2 && u1 /= u2)))
   where
-    t1 :=: u1 = skolemise eq1
-    t2 :=: u2 = skolemise eq2
+    t1 :=: u1 = canonicalise (order eq1)
+    t2 :=: u2 = canonicalise (order eq2)
 
-    skolemise = subst (con . skolem)
+-- | Match one equation against another.
+matchEquation :: Equation f -> Equation f -> Maybe (Subst f)
+matchEquation (pat1 :=: pat2) (t1 :=: t2) = do
+  sub <- match pat1 t1
+  matchIn sub pat2 t2
diff --git a/Twee/Index.hs b/Twee/Index.hs
--- a/Twee/Index.hs
+++ b/Twee/Index.hs
@@ -22,12 +22,12 @@
   lookup,
   matches,
   approxMatches,
-  elems) where
+  elems,
+  fromListWith) where
 
-import qualified Prelude
 import Prelude hiding (null, lookup)
 import Data.Maybe
-import Twee.Base hiding (var, fun, empty, size, singleton, prefix, funs, lookupList, lookup)
+import Twee.Base hiding (var, fun, empty, singleton, prefix, funs, lookupList, lookup)
 import qualified Twee.Term as Term
 import Data.DynamicArray
 import qualified Data.List as List
@@ -83,7 +83,7 @@
 instance Default (Index f a) where def = Nil
 
 -- To get predictable performance, the lookup function uses an explicit stack
--- instead of recursion to control backtracking.
+-- instead of a lazy list to control backtracking.
 data Stack f a =
   -- A normal stack frame: records the current index node and term.
   Frame {
@@ -98,13 +98,15 @@
   | Stop
 
 -- Turn a stack into a list of results.
+{-# SCC run #-}
 run :: Stack f a -> [a]
 run Stop = []
-run Frame{..} = run ({-# SCC run_inner #-} step frame_term frame_index frame_rest)
-run Yield{..} = {-# SCC run_found #-} yield_found ++ run yield_rest
+run Frame{..} = run (step frame_term frame_index frame_rest)
+run Yield{..} = yield_found ++ run yield_rest
 
 -- Execute a single stack frame.
 {-# INLINE step #-}
+{-# SCC step #-}
 step :: TermList f -> Index f a -> Stack f a -> Stack f a
 step !_ _ _ | False = undefined
 step t idx rest =
@@ -118,25 +120,14 @@
 
 -- The main work of 'step' goes on here.
 -- It is carefully tweaked to generate nice code,
--- including using UnsafeCons and only casing on each
--- term list exactly once.
+-- in particular casing on each term list exactly once.
 pref :: TermList f -> TermList f -> [a] -> Array (Index f a) -> Index f a -> Stack f a -> Stack f a
 pref !_ !_ _ !_ !_ _ | False = undefined
 pref search prefix here fun var rest =
   case search of
-    Empty ->
-      case prefix of
-        Empty ->
-          -- The search term matches this node.
-          case here of
-            [] -> rest
-            _ -> Yield here rest
-        _ ->
-          -- We've run out of search term - it doesn't match this node.
-          rest
-    UnsafeCons t ts ->
+    ConsSym{hd = t, tl = ts, rest = ts1} ->
       case prefix of
-        Cons u us ->
+        ConsSym{hd = u, tl = us, rest = us1} ->
           -- Check the search term against the prefix.
           case (t, u) of
             (_, Var _) ->
@@ -145,36 +136,38 @@
               pref ts us here fun var rest
             (App f _, App g _) | f == g ->
               -- Term and prefix start with same symbol, chop them off.
-               let
-                 UnsafeConsSym _ ts' = search
-                 UnsafeConsSym _ us' = prefix
-               in pref ts' us' here fun var rest
+               pref ts1 us1 here fun var rest
             _ ->
               -- Term and prefix don't match.
               rest
         _ ->
           -- We've exhausted the prefix, so let's descend into the tree.
           -- Seems to work better to explore the function node first.
-          let
-            tryVar =
+          case t of
+            App f _ ->
+              case (fun ! fun_id f, var) of
+                (Nil, Nil) ->
+                  rest
+                (Nil, Index{}) ->
+                  step ts var rest
+                (idx, Nil) ->
+                  step ts1 idx rest
+                (idx, Index{}) ->
+                  step ts1 idx (Frame ts var rest)
+            _ ->
               case var of
                 Nil -> rest
-                Index{} -> Frame ts var rest
-              where
-                UnsafeCons _ ts = search
-
-            tryFun =
-              case t of
-                App f _ ->
-                  case fun ! fun_id f of
-                    Nil -> tryVar
-                    idx -> Frame ts idx $! tryVar
-                _ ->
-                  tryVar
-              where
-                UnsafeConsSym t ts = search
-          in
-            tryFun
+                _ -> step ts var rest
+    Empty ->
+      case prefix of
+        Empty ->
+          -- The search term matches this node.
+          case here of
+            [] -> rest
+            _ -> Yield here rest
+        _ ->
+          -- We've run out of search term - it doesn't match this node.
+          rest
 
 -- | An empty index.
 empty :: Index f a
@@ -195,32 +188,44 @@
 singletonList t x = Index 0 t [x] newArray Nil
 
 -- | Insert an entry into the index.
+{-# SCC insert #-}
 insert :: Term f -> a -> Index f a -> Index f a
-insert !t x !idx = {-# SCC insert #-} aux (Term.singleton t) idx
+insert !t x !idx = aux (Term.singleton t) idx
   where
     aux t Nil = singletonList t x
-    aux (Cons t ts) idx@Index{prefix = Cons u us} | t == u =
-      withPrefix (Term.singleton t) (aux ts idx{prefix = us})
+    aux (Cons t ts) idx@Index{prefix = Cons u us}
+      | skeleton t == skeleton u =
+        withPrefix t (aux ts idx{prefix = us})
+    aux (ConsSym{hd = t, rest = ts}) idx@Index{prefix = ConsSym{hd = u, rest = us}}
+      | t `sameSymbolAs` u =
+        withPrefix (build (atom t)) (aux ts idx{prefix = us})
     aux t idx@Index{prefix = Cons{}} = aux t (expand idx)
 
     aux Empty idx =
       idx { size = 0, here = x:here idx }
-    aux t@(ConsSym (App f _) u) idx =
+    aux t@ConsSym{hd = App f _, rest = u} idx =
       idx {
         size = lenList t `min` size idx,
         fun  = update (fun_id f) idx' (fun idx) }
       where
         idx' = aux u (fun idx ! fun_id f)
-    aux t@(ConsSym (Var _) u) idx =
+    aux t@ConsSym{hd = Var _, rest = u} idx =
       idx {
         size = lenList t `min` size idx,
         var  = aux u (var idx) }
 
+    Var _ `sameSymbolAs` Var _ = True
+    App f _ `sameSymbolAs` App g _ = f == g
+    _ `sameSymbolAs` _ = False
+
+    skeleton t = build (subst (const (Term.var (V 0))) t)
+
+    atom (Var x) = Term.var x
+    atom (App f _) = con f
+
 -- Add a prefix to an index.
 -- Does not update the size field.
-{-# INLINE withPrefix #-}
-withPrefix :: TermList f -> Index f a -> Index f a
-withPrefix Empty idx = idx
+withPrefix :: Term f -> Index f a -> Index f a
 withPrefix _ Nil = Nil
 withPrefix t idx@Index{..} =
   idx{prefix = buildList (builder t `mappend` builder prefix)}
@@ -229,7 +234,7 @@
 -- giving an index which doesn't start with a prefix.
 {-# INLINE expand #-}
 expand :: Index f a -> Index f a
-expand idx@Index{size = size, prefix = ConsSym t ts} =
+expand idx@Index{size = size, prefix = ConsSym{hd = t, rest = ts}} =
   case t of
     Var _ ->
       Index {
@@ -248,12 +253,17 @@
 
 -- | Delete an entry from the index.
 {-# INLINEABLE delete #-}
+{-# SCC delete #-}
 delete :: Eq a => Term f -> a -> Index f a -> Index f a
-delete !t x !idx = {-# SCC delete #-} aux (Term.singleton t) idx
+delete !t x !idx = aux (Term.singleton t) idx
   where
     aux _ Nil = Nil
-    aux (Cons t ts) idx@Index{prefix = Cons u us} | t == u =
-      withPrefix (Term.singleton t) (aux ts idx{prefix = us})
+    aux (ConsSym{rest = ts}) idx@Index{prefix = u@ConsSym{rest = us}} =
+      -- The prefix must match, since the term ought to be in the index
+      -- (which is checked in the Empty case below).
+      case aux ts idx{prefix = us} of
+        Nil -> Nil
+        idx -> idx{prefix = u}
     aux _ idx@Index{prefix = Cons{}} = idx
 
     aux Empty idx
@@ -261,9 +271,9 @@
         idx { here = List.delete x (here idx) }
       | otherwise =
         error "deleted term not found in index"
-    aux (ConsSym (App f _) t) idx =
+    aux ConsSym{hd = App f _, rest = t} idx =
       idx { fun = update (fun_id f) (aux t (fun idx ! fun_id f)) (fun idx) }
-    aux (ConsSym (Var _) t) idx =
+    aux ConsSym{hd = Var _, rest = t} idx =
       idx { var = aux t (var idx) }
 
 -- | Look up a term in the index. Finds all key-value such that the search term
@@ -299,9 +309,9 @@
 approxMatches :: Term f -> Index f a -> [a]
 approxMatches t idx = approxMatchesList (Term.singleton t) idx
 
+{-# SCC approxMatchesList #-}
 approxMatchesList :: TermList f -> Index f a -> [a]
 approxMatchesList t idx =
-  {-# SCC approxMatchesList #-}
   run (Frame t idx Stop)
 
 -- | Return all elements of the index.
@@ -309,5 +319,9 @@
 elems Nil = []
 elems idx =
   here idx ++
-  concatMap elems (Prelude.map snd (toList (fun idx))) ++
+  concatMap elems (map snd (toList (fun idx))) ++
   elems (var idx)
+
+-- | Create an index from a list of items
+fromListWith :: (a -> Term f) -> [a] -> Index f a
+fromListWith f xs = foldr (\x -> insert (f x) x) empty xs
diff --git a/Twee/Join.hs b/Twee/Join.hs
--- a/Twee/Join.hs
+++ b/Twee/Join.hs
@@ -1,14 +1,13 @@
 -- | Tactics for joining critical pairs.
-{-# LANGUAGE FlexibleContexts, BangPatterns, RecordWildCards, TypeFamilies #-}
+{-# LANGUAGE FlexibleContexts, BangPatterns, RecordWildCards, TypeFamilies, ScopedTypeVariables #-}
 module Twee.Join where
 
 import Twee.Base
 import Twee.Rule
 import Twee.Equation
-import Twee.Proof(Proof)
 import qualified Twee.Proof as Proof
 import Twee.CP hiding (Config)
-import Twee.Constraints
+import Twee.Constraints hiding (funs)
 import qualified Twee.Index as Index
 import Twee.Index(Index)
 import Twee.Rule.Index(RuleIndex(..))
@@ -16,26 +15,29 @@
 import Data.Maybe
 import Data.Either
 import Data.Ord
-import qualified Data.Set as Set
+import qualified Data.Map.Strict as Map
 
 data Config =
   Config {
     cfg_ground_join :: !Bool,
-    cfg_use_connectedness :: !Bool,
+    cfg_use_connectedness_standalone :: !Bool,
+    cfg_use_connectedness_in_ground_joining :: !Bool,
     cfg_set_join :: !Bool }
 
 defaultConfig :: Config
 defaultConfig =
   Config {
     cfg_ground_join = True,
-    cfg_use_connectedness = True,
+    cfg_use_connectedness_standalone = True,
+    cfg_use_connectedness_in_ground_joining = False,
     cfg_set_join = False }
 
 {-# INLINEABLE joinCriticalPair #-}
+{-# SCC joinCriticalPair #-}
 joinCriticalPair ::
-  (Function f, Has a (Rule f), Has a (Proof f)) =>
+  (Function f, Has a (Rule f)) =>
   Config ->
-  Index f (Equation f) -> RuleIndex f a ->
+  (Index f (Equation f), Index f (Rule f)) -> RuleIndex f a ->
   Maybe (Model f) -> -- A model to try before checking ground joinability
   CriticalPair f ->
   Either
@@ -48,35 +50,49 @@
     -- after successfully joining all instances.
     (Maybe (CriticalPair f), [CriticalPair f])
 joinCriticalPair config eqns idx mmodel cp@CriticalPair{cp_eqn = t :=: u} =
-  {-# SCC joinCriticalPair #-}
   case allSteps config eqns idx cp of
     Nothing ->
       Right (Nothing, [])
     _ | cfg_set_join config &&
-        not (null $ Set.intersection
-          (normalForms (rewrite reduces (index_all idx)) (Set.singleton (reduce (Refl t))))
-          (normalForms (rewrite reduces (index_all idx)) (Set.singleton (reduce (Refl u))))) ->
+        not (null $ Map.intersection
+          (normalForms (rewrite reduces (index_all idx)) (Map.singleton t []))
+          (normalForms (rewrite reduces (index_all idx)) (Map.singleton u []))) ->
       Right (Just cp, [])
     Just cp ->
       case groundJoinFromMaybe config eqns idx mmodel (branches (And [])) cp of
         Left model -> Left (cp, model)
-        Right cps -> Right (Just cp, cps)
+        Right (mcp, cps) -> Right (mcp, cps)
 
 {-# INLINEABLE step1 #-}
 {-# INLINEABLE step2 #-}
 {-# INLINEABLE step3 #-}
 {-# INLINEABLE allSteps #-}
 step1, step2, step3, allSteps ::
-  (Function f, Has a (Rule f), Has a (Proof f)) =>
-  Config -> Index f (Equation f) -> RuleIndex f a -> CriticalPair f -> Maybe (CriticalPair f)
+  (Function f, Has a (Rule f)) =>
+  Config -> (Index f (Equation f), Index f (Rule f)) -> RuleIndex f a -> CriticalPair f -> Maybe (CriticalPair f)
+checkOrder :: Function f => CriticalPair f -> Maybe (CriticalPair f)
 allSteps config eqns idx cp =
   step1 config eqns idx cp >>=
   step2 config eqns idx >>=
+  checkOrder >>=
   step3 config eqns idx
-step1 _ eqns idx = joinWith eqns idx (\t _ -> normaliseWith (const True) (rewrite reducesOriented (index_oriented idx)) t)
-step2 _ eqns idx = joinWith eqns idx (\t _ -> normaliseWith (const True) (rewrite reduces (index_all idx)) t)
-step3 Config{..} eqns idx cp
-  | not cfg_use_connectedness = Just cp
+checkOrder cp
+  | tooBig cp = Nothing
+  | otherwise = Just cp
+  where
+    tooBig CriticalPair{cp_top = Just top, cp_eqn = t :=: u} =
+      lessEq top t || lessEq top u
+    tooBig _ = False
+step1 cfg eqns idx = joinWith cfg eqns idx (\t u -> normaliseWith (const True) (rewrite (ok t u) (index_oriented idx)) t)
+  where
+    --ok _ _ = reducesOriented
+   ok t u rule sub = reducesOriented rule sub && unorient rule `simplerThan` (t :=: u)
+step2 cfg eqns idx = joinWith cfg eqns idx (\t u -> normaliseWith (const True) (rewrite (ok t u) (index_all idx)) t)
+  where
+    --ok _ _ = reduces
+    ok t u rule sub = reduces rule sub && unorient rule `simplerThan` (t :=: u)
+step3 cfg@Config{..} eqns idx cp
+  | not cfg_use_connectedness_standalone = Just cp
   | otherwise =
     case cp_top cp of
       Just top ->
@@ -89,7 +105,7 @@
       _ -> Just cp
   where
     join (cp, top) =
-      joinWith eqns idx (\t u -> normaliseWith (`lessThan` top) (rewrite (ok t u) (index_all idx)) t) cp
+      joinWith cfg eqns idx (\t u -> normaliseWith (`lessThan` top) (rewrite (ok t u) (index_all idx)) t) cp
 
     ok t u rule sub =
       unorient rule `simplerThan` (t :=: u) &&
@@ -104,28 +120,35 @@
 
 {-# INLINEABLE joinWith #-}
 joinWith ::
-  (Has a (Rule f), Has a (Proof f)) =>
-  Index f (Equation f) -> RuleIndex f a -> (Term f -> Term f -> Resulting f) -> CriticalPair f -> Maybe (CriticalPair f)
-joinWith eqns idx reduce cp@CriticalPair{cp_eqn = lhs :=: rhs, ..}
+  (Function f, Has a (Rule f)) =>
+  Config -> (Index f (Equation f), Index f (Rule f)) -> RuleIndex f a -> (Term f -> Term f -> Reduction f) -> CriticalPair f -> Maybe (CriticalPair f)
+joinWith Config{..} eqns idx reduce cp@CriticalPair{cp_eqn = lhs :=: rhs, ..}
   | subsumed eqns idx eqn = Nothing
   | otherwise =
     Just cp {
       cp_eqn = eqn,
       cp_proof =
-        Proof.symm (reductionProof (reduction lred)) `Proof.trans`
+        Proof.symm (reductionProof lhs lred) `Proof.trans`
         cp_proof `Proof.trans`
-        reductionProof (reduction rred) }
+        reductionProof rhs rred }
   where
     lred = reduce lhs rhs
     rred = reduce rhs lhs
-    eqn = result lred :=: result rred
+    eqn = result lhs lred :=: result rhs rred
 
 {-# INLINEABLE subsumed #-}
 subsumed ::
-  (Has a (Rule f), Has a (Proof f)) =>
-  Index f (Equation f) -> RuleIndex f a -> Equation f -> Bool
-subsumed eqns idx (t :=: u)
+  (Has a (Rule f), Function f) =>
+  (Index f (Equation f), Index f (Rule f)) -> RuleIndex f a -> Equation f -> Bool
+subsumed (eqns, complete) idx (t :=: u)
   | t == u = True
+  | otherwise = subsumed1 eqns idx (norm t :=: norm u)
+  where
+    norm t
+      | Index.null complete = t
+      | otherwise = result t $ normaliseWith (const True) (rewrite reducesSkolem complete) t
+subsumed1 eqns idx (t :=: u)
+  | t == u = True
   | or [ rhs rule == u | rule <- Index.lookup t (index_all idx) ] = True
   | or [ rhs rule == t | rule <- Index.lookup u (index_all idx) ] = True
     -- No need to do this symmetrically because addJoinable adds
@@ -133,83 +156,94 @@
   | or [ u == subst sub u'
        | t' :=: u' <- Index.approxMatches t eqns,
          sub <- maybeToList (match t' t) ] = True
-subsumed eqns idx (App f ts :=: App g us)
+subsumed1 eqns idx (App f ts :=: App g us)
   | f == g =
     let
       sub Empty Empty = True
       sub (Cons t ts) (Cons u us) =
-        subsumed eqns idx (t :=: u) && sub ts us
+        subsumed1 eqns idx (t :=: u) && sub ts us
       sub _ _ = error "Function used with multiple arities"
     in
       sub ts us
-subsumed _ _ _ = False
+subsumed1 _ _ _ = False
 
 {-# INLINEABLE groundJoin #-}
 groundJoin ::
-  (Function f, Has a (Rule f), Has a (Proof f)) =>
-  Config -> Index f (Equation f) -> RuleIndex f a -> [Branch f] -> CriticalPair f -> Either (Model f) [CriticalPair f]
+  (Function f, Has a (Rule f)) =>
+  Config -> (Index f (Equation f), Index f (Rule f)) -> RuleIndex f a -> [Branch f] -> CriticalPair f -> Either (Model f) (Maybe (CriticalPair f), [CriticalPair f])
 groundJoin config eqns idx ctx cp@CriticalPair{cp_eqn = t :=: u, ..} =
   case partitionEithers (map (solve (usort (atoms t ++ atoms u))) ctx) of
     ([], instances) ->
       let cps = [ subst sub cp | sub <- instances ] in
-      Right (usortBy (comparing (canonicalise . order . cp_eqn)) cps)
+      Right (Just cp, usortBy (comparing (canonicalise . order . cp_eqn)) cps)
     (model:_, _) ->
       groundJoinFrom config eqns idx model ctx cp
 
 {-# INLINEABLE groundJoinFrom #-}
 groundJoinFrom ::
-  (Function f, Has a (Rule f), Has a (Proof f)) =>
-  Config -> Index f (Equation f) -> RuleIndex f a -> Model f -> [Branch f] -> CriticalPair f -> Either (Model f) [CriticalPair f]
+  (Function f, Has a (Rule f)) =>
+  Config -> (Index f (Equation f), Index f (Rule f)) -> RuleIndex f a -> Model f -> [Branch f] -> CriticalPair f -> Either (Model f) (Maybe (CriticalPair f), [CriticalPair f])
 groundJoinFrom config@Config{..} eqns idx model ctx cp@CriticalPair{cp_eqn = t :=: u, ..}
-  | not cfg_ground_join ||
-    (modelOK model && isJust (allSteps config eqns idx cp { cp_eqn = t' :=: u' })) = Left model
+  | not cfg_ground_join = Left model
+  | modelOK model && isJust (allSteps config' eqns idx cp { cp_eqn = t' :=: u' }) = Left model
   | otherwise =
-      let model1 = optimise model weakenModel (\m -> not (modelOK m) || (valid m (reduction nt) && valid m (reduction nu)))
-          model2 = optimise model1 weakenModel (\m -> not (modelOK m) || isNothing (allSteps config eqns idx cp { cp_eqn = result (normaliseIn m t u) :=: result (normaliseIn m u t) }))
+      let
+        model'
+          | modelOK model =
+            optimise weakenModel (\m -> modelOK m && isNothing (allSteps config' eqns idx cp { cp_eqn = result t (normaliseIn m t u) :=: result u (normaliseIn m u t) })) $
+            optimise weakenModel (\m -> modelOK m && (valid m nt && valid m nu)) model
+          | otherwise =
+            optimise weakenModel (not . modelOK) model
 
-          diag [] = Or []
-          diag (r:rs) = negateFormula r ||| (weaken r &&& diag rs)
-          weaken (LessEq t u) = Less t u
-          weaken x = x
-          ctx' = formAnd (diag (modelToLiterals model2)) ctx in
+        diag [] = Or []
+        diag (r:rs) = negateFormula r ||| (weaken r &&& diag rs)
+        weaken (LessEq t u) = Less t u
+        weaken x = x
+        ctx' = formAnd (diag (modelToLiterals model')) ctx in
 
-      groundJoin config eqns idx ctx' cp
+      case groundJoin config eqns idx ctx' cp of
+        Right (_, cps) | not (modelOK model) ->
+          Right (Nothing, cps)
+        res -> res
   where
-    normaliseIn m t u = normaliseWith (const True) (rewrite (ok t u m) (index_all idx)) t
+    config' = config{cfg_use_connectedness_standalone = False}
+    normaliseIn m t u =
+      case cp_top of
+        Just top | cfg_use_connectedness_in_ground_joining ->
+          normaliseWith (connectedIn m top) (rewrite (ok t u model) (index_all idx)) t
+        _ -> normaliseWith (const True) (rewrite (ok t u m) (index_all idx)) t
     ok t u m rule sub =
       reducesInModel m rule sub &&
       unorient rule `simplerThan` (t :=: u)
+    connectedIn m top t =
+      lessIn m t top == Just Strict
 
     nt = normaliseIn model t u
     nu = normaliseIn model u t
-    t' = result nt
-    u' = result nu
+    t' = result t nt
+    u' = result u nu
 
-    -- XXX not safe to exploit the top term if we then add the equation to
-    -- the joinable set. (It might then be used to join a CP with an entirely
-    -- different top term.)
-    modelOK _ = True
-{-    modelOK m =
+    modelOK m =
       case cp_top of
         Nothing -> True
         Just top ->
-          isNothing (lessIn m top t) && isNothing (lessIn m top u)-}
+          isNothing (lessIn m top t) && isNothing (lessIn m top u)
 
 {-# INLINEABLE groundJoinFromMaybe #-}
 groundJoinFromMaybe ::
-  (Function f, Has a (Rule f), Has a (Proof f)) =>
-  Config -> Index f (Equation f) -> RuleIndex f a -> Maybe (Model f) -> [Branch f] -> CriticalPair f -> Either (Model f) [CriticalPair f]
+  (Function f, Has a (Rule f)) =>
+  Config -> (Index f (Equation f), Index f (Rule f)) -> RuleIndex f a -> Maybe (Model f) -> [Branch f] -> CriticalPair f -> Either (Model f) (Maybe (CriticalPair f), [CriticalPair f])
 groundJoinFromMaybe config eqns idx Nothing = groundJoin config eqns idx
 groundJoinFromMaybe config eqns idx (Just model) = groundJoinFrom config eqns idx model
 
 {-# INLINEABLE valid #-}
 valid :: Function f => Model f -> Reduction f -> Bool
 valid model red =
-  and [ reducesInModel model rule sub
-      | Step _ rule sub <- steps red ]
+  and [ reducesInModel model rule emptySubst
+      | rule <- red ]
 
-optimise :: a -> (a -> [a]) -> (a -> Bool) -> a
-optimise x f p =
+optimise :: (a -> [a]) -> (a -> Bool) -> a -> a
+optimise f p x =
   case filter p (f x) of
-    y:_ -> optimise y f p
+    y:_ -> optimise f p y
     _   -> x
diff --git a/Twee/KBO.hs b/Twee/KBO.hs
--- a/Twee/KBO.hs
+++ b/Twee/KBO.hs
@@ -1,18 +1,45 @@
 -- | An implementation of Knuth-Bendix ordering.
 
-{-# LANGUAGE PatternGuards #-}
-module Twee.KBO(lessEq, lessIn) where
+{-# LANGUAGE PatternGuards, BangPatterns #-}
+module Twee.KBO(lessEq, lessIn, lessEqSkolem, Sized(..), Weighted(..)) where
 
-import Twee.Base hiding (lessEq, lessIn)
-import Data.List
-import Twee.Constraints hiding (lessEq, lessIn)
+import Twee.Base hiding (lessEq, lessIn, lessEqSkolem)
+import Twee.Equation
+import Twee.Constraints hiding (lessEq, lessIn, lessEqSkolem)
 import qualified Data.Map.Strict as Map
 import Data.Map.Strict(Map)
 import Data.Maybe
 import Control.Monad
+import Twee.Utils
 
+lessEqSkolem :: (Function f, Sized f, Weighted f) => Term f -> Term f -> Bool
+lessEqSkolem !t !u
+  | m < n = True
+  | m > n = False
+  where
+    m = size t
+    n = size u
+lessEqSkolem (App x Empty) _
+  | x == minimal = True
+lessEqSkolem _ (App x Empty)
+  | x == minimal = False
+lessEqSkolem (Var x) (Var y) = x <= y
+lessEqSkolem (Var _) _ = True
+lessEqSkolem _ (Var _) = False
+lessEqSkolem (App (F _ f) ts) (App (F _ g) us) =
+  case compare f g of
+    LT -> True
+    GT -> False
+    EQ ->
+      let loop Empty Empty = True
+          loop (Cons t ts) (Cons u us)
+            | t == u = loop ts us
+            | otherwise = lessEqSkolem t u
+      in loop ts us
+
 -- | Check if one term is less than another in KBO.
-lessEq :: Function f => Term f -> Term f -> Bool
+{-# SCC lessEq #-}
+lessEq :: (Function f, Sized f, Weighted f) => Term f -> Term f -> Bool
 lessEq (App f Empty) _ | f == minimal = True
 lessEq (Var x) (Var y) | x == y = True
 lessEq _ (Var _) = False
@@ -21,7 +48,7 @@
   (st < su ||
    (st == su && f << g) ||
    (st == su && f == g && lexLess ts us)) &&
-  xs `isSubsequenceOf` ys
+  xs `lessVars` ys
   where
     lexLess Empty Empty = True
     lexLess (Cons t ts) (Cons u us)
@@ -34,23 +61,30 @@
             | not (allSubst (\_ (Cons t Empty) -> isMinimal t) sub) -> error "weird term inequality"
             | otherwise -> lexLess (subst sub ts) (subst sub us)
     lexLess _ _ = error "incorrect function arity"
-    xs = sort (vars t)
-    ys = sort (vars u)
+    xs = weightedVars t
+    ys = weightedVars u
     st = size t
     su = size u
 
+    [] `lessVars` _ = True
+    ((x,k1):xs) `lessVars` ((y,k2):ys)
+      | x == y = k1 <= k2 && xs `lessVars` ys
+      | x > y  = ((x,k1):xs) `lessVars` ys
+    _ `lessVars` _ = False
+
 -- | Check if one term is less than another in a given model.
 
 -- See "notes/kbo under assumptions" for how this works.
 
-lessIn :: Function f => Model f -> Term f -> Term f -> Maybe Strictness
+{-# SCC lessIn #-}
+lessIn :: (Function f, Sized f, Weighted f) => Model f -> Term f -> Term f -> Maybe Strictness
 lessIn model t u =
   case sizeLessIn model t u of
     Nothing -> Nothing
     Just Strict -> Just Strict
     Just Nonstrict -> lexLessIn model t u
 
-sizeLessIn :: Function f => Model f -> Term f -> Term f -> Maybe Strictness
+sizeLessIn :: (Function f, Sized f, Weighted f) => Model f -> Term f -> Term f -> Maybe Strictness
 sizeLessIn model t u =
   case minimumIn model m of
     Just l
@@ -59,13 +93,13 @@
     _ -> Nothing
   where
     (k, m) =
-      foldr (addSize id)
-        (foldr (addSize negate) (0, Map.empty) (subterms t))
-        (subterms u)
-    addSize op (App f _) (k, m) = (k + op (size f), m)
-    addSize op (Var x) (k, m) = (k, Map.insertWith (+) x (op 1) m)
+      add 1 u (add (-1) t (0, Map.empty))
 
-minimumIn :: Function f => Model f -> Map Var Int -> Maybe Int
+    add a (App f ts) (k, m) =
+      foldr (add (a * argWeight f)) (k + a * size f, m) (unpack ts)
+    add a (Var x) (k, m) = (k, Map.insertWith (+) x a m)
+
+minimumIn :: (Function f, Sized f) => Model f -> Map Var Integer -> Maybe Integer
 minimumIn model t =
   liftM2 (+)
     (fmap sum (mapM minGroup (varGroups model)))
@@ -90,7 +124,7 @@
       | k < 0 = Nothing
       | otherwise = Just k
 
-lexLessIn :: Function f => Model f -> Term f -> Term f -> Maybe Strictness
+lexLessIn :: (Function f, Sized f, Weighted f) => Model f -> Term f -> Term f -> Maybe Strictness
 lexLessIn _ t u | t == u = Just Nonstrict
 lexLessIn cond t u
   | Just a <- fromTerm t,
@@ -119,3 +153,38 @@
     loop _ _ = error "incorrect function arity"
 lexLessIn _ t _ | isMinimal t = Just Nonstrict
 lexLessIn _ _ _ = Nothing
+
+class Sized a where
+  -- | Compute the size.
+  size  :: a -> Integer
+
+class Weighted f where
+  argWeight :: f -> Integer
+
+instance (Weighted f, Labelled f) => Weighted (Fun f) where
+  argWeight = argWeight . fun_value
+
+weightedVars :: (Weighted f, Labelled f) => Term f -> [(Var, Integer)]
+weightedVars t = collate sum (loop 1 t)
+  where
+    loop k (Var x) = [(x, k)]
+    loop k (App f ts) =
+      concatMap (loop (k * argWeight f)) (unpack ts)
+
+instance (Labelled f, Sized f) => Sized (Fun f) where
+  size = size . fun_value
+
+instance (Labelled f, Sized f, Weighted f) => Sized (TermList f) where
+  size = aux 0
+    where
+      aux n Empty = n
+      aux n (Cons (App f t) u) =
+        aux (n + size f + argWeight f * size t) u
+      aux n (Cons (Var _) t) = aux (n+1) t
+
+instance (Labelled f, Sized f, Weighted f) => Sized (Term f) where
+  size = size . singleton
+
+instance (Labelled f, Sized f, Weighted f) => Sized (Equation f) where
+  size (x :=: y) = size x + size y
+
diff --git a/Twee/Label.hs b/Twee/Label.hs
--- a/Twee/Label.hs
+++ b/Twee/Label.hs
@@ -8,8 +8,8 @@
 import System.IO.Unsafe
 import qualified Data.Map.Strict as Map
 import Data.Map.Strict(Map)
-import qualified Data.IntMap.Strict as IntMap
-import Data.IntMap.Strict(IntMap)
+import qualified Data.DynamicArray as DynamicArray
+import Data.DynamicArray(Array)
 import Data.Typeable
 import GHC.Exts
 import Unsafe.Coerce
@@ -30,7 +30,7 @@
 -- The global cache of labels.
 {-# NOINLINE cachesRef #-}
 cachesRef :: IORef Caches
-cachesRef = unsafePerformIO (newIORef (Caches 0 Map.empty IntMap.empty))
+cachesRef = unsafePerformIO (newIORef (Caches 0 Map.empty DynamicArray.newArray))
 
 data Caches =
   Caches {
@@ -39,7 +39,7 @@
     -- A map from values to labels.
     caches_from   :: !(Map TypeRep (Cache Any)),
     -- The reverse map from labels to values.
-    caches_to     :: !(IntMap Any) }
+    caches_to     :: !(Array Any) }
 
 type Cache a = Map a Int32
 
@@ -105,7 +105,7 @@
       (caches {
          caches_nextId = n+1,
          caches_from = Map.insert ty (toAnyCache (Map.insert x n cache)) caches_from,
-         caches_to = IntMap.insert (fromIntegral n) (toAny x) caches_to },
+         caches_to = DynamicArray.updateWithDefault undefined (fromIntegral n) (toAny x) caches_to },
        Label n)
       where
         n = caches_nextId
@@ -121,5 +121,5 @@
 -- doesn't work.
 find (Label !n) = unsafeDupablePerformIO $ do
   Caches{..} <- readIORef cachesRef
-  x <- return $! fromAny (IntMap.findWithDefault undefined (fromIntegral n) caches_to)
+  x <- return $! fromAny (DynamicArray.getWithDefault undefined (fromIntegral n) caches_to)
   return x
diff --git a/Twee/PassiveQueue.hs b/Twee/PassiveQueue.hs
--- a/Twee/PassiveQueue.hs
+++ b/Twee/PassiveQueue.hs
@@ -4,7 +4,7 @@
   Params(..),
   Queue,
   Passive(..),
-  empty, insert, removeMin, mapMaybe) where
+  empty, insert, removeMin, mapMaybe, toList, queueSize) where
 
 import qualified Data.Heap as Heap
 import qualified Data.Vector.Unboxed as Vector
@@ -115,6 +115,18 @@
        packId proxy (if isLeft then passive_rule2 else passive_rule1),
        fromIntegral passive_pos)
 
+-- Convert a PassiveSet back into a list of Passives.
+{-# INLINEABLE unpackPassiveSet #-}
+unpackPassiveSet :: forall params.Params params => PassiveSet params -> (Int, [Passive params])
+unpackPassiveSet PassiveSet{..} =
+  (1 + Vector.length passiveset_left + Vector.length passiveset_right,
+   passiveset_best:
+   map (unpack proxy passiveset_rule True) (Vector.toList passiveset_left) ++
+   map (unpack proxy passiveset_rule False) (Vector.toList passiveset_right))
+  where
+    proxy :: Proxy params
+    proxy = Proxy
+
 -- Find and remove the best element from a PassiveSet.
 {-# INLINEABLE unconsPassiveSet #-}
 unconsPassiveSet :: forall params. Params params => PassiveSet params -> (Passive params, Maybe (PassiveSet params))
@@ -174,10 +186,16 @@
 mapMaybe :: Params params => (Passive params -> Maybe (Passive params)) -> Queue params -> Queue params
 mapMaybe f (Queue q) = Queue (Heap.mapMaybe g q)
   where
-    g PassiveSet{..} =
+    g passiveSet@PassiveSet{..} =
       makePassiveSet passiveset_rule $ Data.Maybe.mapMaybe f $
-        passiveset_best:
-        map (unpack proxy passiveset_rule True) (Vector.toList passiveset_left) ++
-        map (unpack proxy passiveset_rule False) (Vector.toList passiveset_right)
-    proxy :: Proxy params
-    proxy = Proxy
+        snd (unpackPassiveSet passiveSet)
+
+-- | Convert a queue into a list of 'Passive's.
+-- The 'Passive's are produced in batches, with each batch labelled
+-- with its size.
+{-# INLINEABLE toList #-}
+toList :: Params params => Queue params -> [(Int, [Passive params])]
+toList (Queue h) = map unpackPassiveSet (Heap.toList h)
+
+queueSize :: Params params => Queue params -> Int
+queueSize = sum . map fst . toList
diff --git a/Twee/Pretty.hs b/Twee/Pretty.hs
--- a/Twee/Pretty.hs
+++ b/Twee/Pretty.hs
@@ -69,17 +69,44 @@
 
 -- * Pretty-printing of terms.
 
-instance Pretty f => Pretty (Fun f) where
+instance (Pretty f, Labelled f) => Pretty (Fun f) where
   pPrintPrec l p = pPrintPrec l p . fun_value
 
-instance PrettyTerm f => PrettyTerm (Fun f) where
-  termStyle f = termStyle (fun_value f)
-
 instance PrettyTerm f => Pretty (Term f) where
   pPrintPrec l p (Var x) = pPrintPrec l p x
   pPrintPrec l p (App f xs) =
-    pPrintTerm (termStyle f) l p (pPrint f) (unpack xs)
+    pPrintTerm (termStyle (fun_value f)) l p (pPrint f) (unpack xs)
 
+data HighlightedTerm f = HighlightedTerm [ANSICode] (Maybe [Int]) (Term f)
+
+type ANSICode = String
+green, bold :: ANSICode
+green = "32"
+bold = "1"
+
+highlight :: [ANSICode] -> Doc -> Doc
+highlight cs d =
+  hsep (map escape cs) <#> d <#> hsep [escape "" | not (null cs)]
+  where
+    escape s = zeroWidthText ("\027[" ++ s ++ "m")
+
+maybeHighlight :: [ANSICode] -> Maybe [Int] -> Doc -> Doc
+maybeHighlight cs (Just []) d = highlight cs d
+maybeHighlight _ _ d = d
+
+instance PrettyTerm f => Pretty (HighlightedTerm f) where
+  pPrintPrec l p (HighlightedTerm cs h (Var x)) =
+    maybeHighlight cs h (pPrintPrec l p x)
+  pPrintPrec l p (HighlightedTerm cs h (App f xs)) =
+    maybeHighlight cs h $
+    pPrintTerm (termStyle (fun_value f)) l p (pPrint f)
+      (zipWith annotate [0..] (unpack xs))
+    where
+      annotate i t =
+        case h of
+          Just (n:ns) | i == n -> HighlightedTerm cs (Just ns) t
+          _ -> HighlightedTerm cs Nothing t
+
 instance PrettyTerm f => Pretty (TermList f) where
   pPrintPrec _ _ = pPrint . unpack
 
@@ -91,7 +118,7 @@
         | (x, t) <- substToList sub ]
 
 -- | A class for customising the printing of function symbols.
-class Pretty f => PrettyTerm f where
+class (Pretty f, Labelled f) => PrettyTerm f where
   -- | The style of the function symbol. Defaults to 'curried'.
   termStyle :: f -> TermStyle
   termStyle _ = curried
diff --git a/Twee/Proof.hs b/Twee/Proof.hs
--- a/Twee/Proof.hs
+++ b/Twee/Proof.hs
@@ -5,10 +5,11 @@
   Proof, Derivation(..), Axiom(..),
   certify, equation, derivation,
   -- ** Smart constructors for derivations
-  lemma, axiom, symm, trans, cong, congPath,
+  lemma, autoSubst, simpleLemma, axiom, symm, trans, cong, congPath,
 
   -- * Analysing proofs
-  simplify, usedLemmas, usedAxioms, usedLemmasAndSubsts, usedAxiomsAndSubsts,
+  simplify, steps, usedLemmas, usedAxioms, usedLemmasAndSubsts, usedAxiomsAndSubsts,
+  groundAxiomsAndSubsts, eliminateDefinitions, eliminateDefinitionsFromGoal,
 
   -- * Pretty-printing proofs
   Config(..), defaultConfig, Presentation(..),
@@ -18,12 +19,18 @@
 import Twee.Base hiding (invisible)
 import Twee.Equation
 import Twee.Utils
+import qualified Twee.Index as Index
 import Control.Monad
 import Data.Maybe
 import Data.List
 import Data.Ord
 import qualified Data.Set as Set
+import Data.Set(Set)
 import qualified Data.Map.Strict as Map
+import Data.Map(Map)
+import qualified Data.IntMap.Strict as IntMap
+import Control.Monad.Trans.State.Strict
+import Data.Graph
 
 ----------------------------------------------------------------------
 -- Equational proofs. Only valid proofs can be constructed.
@@ -78,9 +85,9 @@
 
 -- This is the trusted core of the module.
 {-# INLINEABLE certify #-}
+{-# SCC certify #-}
 certify :: PrettyTerm f => Derivation f -> Proof f
 certify p =
-  {-# SCC certify #-}
   case check p of
     Nothing -> error ("Invalid proof created!\n" ++ prettyShow p)
     Just eqn -> Proof eqn p
@@ -141,7 +148,7 @@
   pPrint = pPrintLemma defaultConfig (prettyShow . axiom_number) (prettyShow . equation)
 instance PrettyTerm f => Pretty (Derivation f) where
   pPrint (UseLemma lemma sub) =
-    text "subst" <#> pPrintTuple [text "lemma" <#> pPrint (equation lemma), pPrint sub]
+    text "subst" <#> pPrintTuple [text "lemma" <+> pPrint (equation lemma), pPrint sub]
   pPrint (UseAxiom axiom sub) =
     text "subst" <#> pPrintTuple [pPrint axiom, pPrint sub]
   pPrint (Refl t) =
@@ -158,39 +165,69 @@
     text "axiom" <#>
     pPrintTuple [pPrint axiom_number, text axiom_name, pPrint axiom_eqn]
 
--- | Simplify a derivation.
---
--- After simplification, a derivation has the following properties:
---
---   * 'Symm' is pushed down next to 'Lemma' and 'Axiom'
---   * 'Refl' only occurs inside 'Cong' or at the top level
---   * 'Trans' is right-associated and is pushed inside 'Cong' if possible
-simplify :: Minimal f => (Proof f -> Maybe (Derivation f)) -> Derivation f -> Derivation f
-simplify lem p = simp p
+foldLemmas :: PrettyTerm f => (Map (Proof f) a -> Derivation f -> a) -> [Derivation f] -> Map (Proof f) a
+foldLemmas op ds =
+  execState (mapM_ foldGoal ds) Map.empty
   where
-    simp p@(UseLemma q sub) =
-      case lem q of
-        Nothing -> p
-        Just r ->
-          let
-            -- Get rid of any variables that are not bound by sub
-            -- (e.g., ones which only occur internally in q)
-            dead = usort (vars r) \\ substDomain sub
-          in simp (subst sub (erase dead r))
-    simp (Symm p) = symm (simp p)
-    simp (Trans p q) = trans (simp p) (simp q)
-    simp (Cong f ps) = cong f (map simp ps)
-    simp p = p
+    foldGoal p = mapM_ foldLemma (usedLemmas p)
+    foldLemma p = do
+      m <- get
+      case Map.lookup p m of
+        Just x -> return x
+        Nothing -> do
+          mapM_ foldLemma (usedLemmas (derivation p))
+          m <- get
+          case Map.lookup p m of
+            Just x  -> return x
+            Nothing -> do
+              let x = op m (derivation p)
+              put (Map.insert p x m)
+              return x
 
+mapLemmas :: Function f => (Derivation f -> Derivation f) -> [Derivation f] -> [Derivation f]
+mapLemmas f ds = map (derivation . op lem) ds
+  where
+    op lem = certify . f . unfoldLemmas (\pf -> Just (simpleLemma (lem Map.! pf)))
+    lem = foldLemmas op ds
+
+allLemmas :: PrettyTerm f => [Derivation f] -> [Proof f]
+allLemmas ds =
+  reverse [p | (_, p, _) <- map vertex (topSort graph)]
+  where
+    used = foldLemmas (\_ p -> usedLemmas p) ds
+    (graph, vertex, _) =
+      graphFromEdges
+        [((), p, ps) | (p, ps) <- Map.toList used]
+
+unfoldLemmas :: Minimal f => (Proof f -> Maybe (Derivation f)) -> Derivation f -> Derivation f
+unfoldLemmas lem p@(UseLemma q sub) =
+  case lem q of
+    Nothing -> p
+    Just r ->
+      -- Get rid of any variables that are not bound by sub
+      -- (e.g., ones which only occur internally in q)
+      subst sub (eraseExcept (substDomain sub) r)
+unfoldLemmas lem (Symm p) = symm (unfoldLemmas lem p)
+unfoldLemmas lem (Trans p q) = trans (unfoldLemmas lem p) (unfoldLemmas lem q)
+unfoldLemmas lem (Cong f ps) = cong f (map (unfoldLemmas lem) ps)
+unfoldLemmas _ p = p
+
 lemma :: Proof f -> Subst f -> Derivation f
 lemma p sub = UseLemma p sub
 
+simpleLemma :: PrettyTerm f => Proof f -> Derivation f
+simpleLemma p =
+  UseLemma p (autoSubst (equation p))
+
 axiom :: Axiom f -> Derivation f
 axiom ax@Axiom{..} =
-  UseAxiom ax $
-    fromJust $
-    listToSubst [(x, build (var x)) | x <- vars axiom_eqn]
+  UseAxiom ax (autoSubst axiom_eqn)
 
+autoSubst :: Equation f -> Subst f
+autoSubst eqn =
+  fromJust $
+  listToSubst [(x, build (var x)) | x <- vars eqn]
+
 symm :: Derivation f -> Derivation f
 symm (Refl t) = Refl t
 symm (Symm p) = p
@@ -206,17 +243,8 @@
   -- p cannot be a Trans (if it was created with the smart
   -- constructors) but q could be.
   Trans p (trans q r)
--- Collect adjacent uses of congruence.
-trans (Cong f ps) (Cong g qs) | f == g =
-  transCong f ps qs
-trans (Cong f ps) (Trans (Cong g qs) r) | f == g =
-  trans (transCong f ps qs) r
 trans p q = Trans p q
 
-transCong :: Fun f -> [Derivation f] -> [Derivation f] -> Derivation f
-transCong f ps qs =
-  cong f (zipWith trans ps qs)
-
 cong :: Fun f -> [Derivation f] -> Derivation f
 cong f ps =
   case sequence (map unRefl ps) of
@@ -226,6 +254,62 @@
     unRefl (Refl t) = Just t
     unRefl _ = Nothing
 
+-- Transform a proof so that each step uses exactly one axiom
+-- or lemma. The proof will have the following form afterwards:
+--   * Trans only occurs at the outermost level and is right-associated
+--   * Each Cong has exactly one non-Refl argument (no parallel rewriting)
+--   * Symm only occurs innermost, i.e., next to UseLemma or UseAxiom
+--   * Refl only occurs as an argument to Cong, or outermost if the
+--     whole proof is a single reflexivity step
+flattenDerivation :: Function f => Derivation f -> Derivation f
+flattenDerivation p =
+  fromSteps (equation (certify p)) (steps p)
+
+-- | Simplify a derivation so that:
+--   * Symm occurs innermost
+--   * Trans is right-associated
+--   * Each Cong has at least one non-Refl argument
+--   * Refl is not used unnecessarily
+simplify :: PrettyTerm f => Derivation f -> Derivation f
+simplify (Symm p) = symm (simplify p)
+simplify (Trans p q) = trans (simplify p) (simplify q)
+simplify (Cong f ps) = cong f (map simplify ps)
+simplify p
+  | t == u = Refl t
+  | otherwise = p
+  where
+    t :=: u = equation (certify p)
+
+-- | Transform a derivation into a list of single steps.
+--   Each step has the following form:
+--     * Trans does not occur
+--     * Symm only occurs innermost, i.e., next to UseLemma or UseAxiom
+--     * Each Cong has exactly one non-Refl argument (no parallel rewriting)
+--     * Refl only occurs as an argument to Cong
+steps :: Function f => Derivation f -> [Derivation f]
+steps = steps1 . simplify
+  where
+    steps1 p@UseAxiom{} = [p]
+    steps1 p@UseLemma{} = [p]
+    steps1 (Refl _) = []
+    steps1 (Symm p) = map symm (reverse (steps1 p))
+    steps1 (Trans p q) = steps1 p ++ steps1 q
+    steps1 p@(Cong f qs) =
+      concat [ map (inside i) (steps1 q) | (i, q) <- zip [0..] qs ]
+      where
+        App _ ts :=: App _ us = equation (certify p)
+        inside i p =
+          Cong f $
+            map Refl (take i (unpack us)) ++
+            [p] ++
+            map Refl (drop (i+1) (unpack ts))
+
+-- | Convert a list of steps (plus the equation it is proving)
+-- back to a derivation.
+fromSteps :: Equation f -> [Derivation f] -> Derivation f
+fromSteps (t :=: _) [] = Refl t
+fromSteps _ ps = foldr1 Trans ps
+
 -- | Find all lemmas which are used in a derivation.
 usedLemmas :: Derivation f -> [Proof f]
 usedLemmas p = map fst (usedLemmasAndSubsts p)
@@ -256,6 +340,79 @@
     ax (Cong _ ps) = foldr (.) id (map ax ps)
     ax _ = id
 
+-- | Find all ground instances of axioms which are used in the
+-- expanded form of a derivation (no lemmas).
+groundAxiomsAndSubsts :: Function f => Derivation f -> Map (Axiom f) (Set (Subst f))
+groundAxiomsAndSubsts p = ax lem p
+  where
+    lem = foldLemmas ax [p]
+
+    ax _ (UseAxiom axiom sub) =
+      Map.singleton axiom (Set.singleton sub)
+    ax lem (UseLemma lemma sub) =
+      Map.map (Set.map substAndErase) (lem Map.! lemma)
+      where
+        substAndErase sub' =
+          eraseExcept (vars sub) (subst sub sub')
+    ax lem (Symm p) = ax lem p
+    ax lem (Trans p q) = Map.unionWith Set.union (ax lem p) (ax lem q)
+    ax lem (Cong _ ps) = Map.unionsWith Set.union (map (ax lem) ps)
+    ax _ _ = Map.empty
+
+eliminateDefinitionsFromGoal :: Function f => [Axiom f] -> ProvedGoal f -> ProvedGoal f
+eliminateDefinitionsFromGoal axioms pg =
+  pg {
+    pg_proof = certify (eliminateDefinitions axioms (derivation (pg_proof pg))) }
+
+eliminateDefinitions :: Function f => [Axiom f] -> Derivation f -> Derivation f
+eliminateDefinitions [] p = p
+eliminateDefinitions axioms p = head (mapLemmas elim [p])
+  where
+    elim (UseAxiom axiom sub)
+      | axiom `Set.member` axSet =
+        Refl (term (subst sub (eqn_rhs (axiom_eqn axiom))))
+      | otherwise = UseAxiom axiom (elimSubst sub)
+    elim (UseLemma lemma sub) =
+      UseLemma lemma (elimSubst sub)
+    elim (Refl t) = Refl (term t)
+    elim (Symm p) = Symm (elim p)
+    elim (Trans p q) = Trans (elim p) (elim q)
+    elim (Cong f ps) =
+      case find (build (app f (map var vs))) of
+        Nothing -> Cong f (map elim ps)
+        Just (rhs, Subst sub) ->
+          let proof (Cons (Var (V x)) Empty) = qs !! x in
+          replace (proof <$> sub) rhs
+      where
+        vs = map V [0..length ps-1]
+        qs = map (simpleLemma . certify . elim) ps -- avoid duplicating proofs of ts
+
+    elimSubst (Subst sub) = Subst (singleton <$> term <$> unsingleton <$> sub)
+      where
+        unsingleton (Cons t Empty) = t
+
+    term = build . term'
+    term' (Var x) = var x
+    term' t@(App f ts) =
+      case find t of
+        Nothing -> app f (map term' (unpack ts))
+        Just (rhs, sub) ->
+          term' (subst sub rhs)
+
+    find t =
+      listToMaybe $ do
+        Axiom{axiom_eqn = l :=: r} <- Index.approxMatches t idx
+        sub <- maybeToList (match l t)
+        return (r, sub)
+
+    replace sub (Var (V x)) =
+      IntMap.findWithDefault undefined x sub
+    replace sub (App f ts) =
+      cong f (map (replace sub) (unpack ts))
+
+    axSet = Set.fromList axioms
+    idx = Index.fromListWith (eqn_lhs . axiom_eqn) axioms
+
 -- | Applies a derivation at a particular path in a term.
 congPath :: [Int] -> Term f -> Derivation f -> Derivation f
 congPath [] _ p = p
@@ -273,22 +430,31 @@
 ----------------------------------------------------------------------
 
 -- | Options for proof presentation.
-data Config =
+data Config f =
   Config {
     -- | Never inline lemmas.
     cfg_all_lemmas :: !Bool,
     -- | Inline all lemmas.
     cfg_no_lemmas :: !Bool,
+    -- | Make the proof ground.
+    cfg_ground_proof :: !Bool,
     -- | Print out explicit substitutions.
-    cfg_show_instances :: !Bool }
+    cfg_show_instances :: !Bool,
+    -- | Print out proofs in colour.
+    cfg_use_colour :: !Bool,
+    -- | Print out which instances of some axioms were used.
+    cfg_show_uses_of_axioms :: Axiom f -> Bool }
 
 -- | The default configuration.
-defaultConfig :: Config
+defaultConfig :: Config f
 defaultConfig =
   Config {
     cfg_all_lemmas = False,
     cfg_no_lemmas = False,
-    cfg_show_instances = False }
+    cfg_ground_proof = False,
+    cfg_show_instances = False,
+    cfg_use_colour = False,
+    cfg_show_uses_of_axioms = const False }
 
 -- | A proof, with all axioms and lemmas explicitly listed.
 data Presentation f =
@@ -346,220 +512,290 @@
   pPrint = pPrintPresentation defaultConfig
 
 -- | Simplify and present a proof.
-present :: Function f => Config -> [ProvedGoal f] -> Presentation f
-present config goals =
-  -- First find all the used lemmas, then hand off to presentWithGoals
-  presentWithGoals config goals
-    (snd (used Set.empty (concatMap (usedLemmas . derivation . pg_proof) goals)))
+present :: Function f => Config f -> [ProvedGoal f] -> Presentation f
+present config@Config{..} goals =
+  Presentation axioms lemmas goals'
   where
-    used lems [] = (lems, [])
-    used lems (x:xs)
-      | x `Set.member` lems = used lems xs
-      | otherwise =
-        let (lems1, ys) = used (Set.insert x lems) (usedLemmas (derivation x))
-            (lems2, zs) = used lems1 xs
-        in (lems2, ys ++ [x] ++ zs)
+    ps =
+      mapLemmas flattenDerivation $
+      simplifyProof config $ map (derivation . pg_proof) goals
 
-presentWithGoals ::
+    goals' =
+      [ decodeGoal (goal{pg_proof = certify p})
+      | (goal, p) <- zip goals ps ]
+
+    axioms = usort $
+      concatMap (usedAxioms . derivation . pg_proof) goals' ++
+      concatMap (usedAxioms . derivation) lemmas
+
+    lemmas = allLemmas (map (derivation . pg_proof) goals')
+
+groundProof :: Function f => [Derivation f] -> [Derivation f]
+groundProof ds
+  | all (isGround . equation) (allLemmas ds) = ds
+  | otherwise = groundProof (mapLemmas f ds)
+  where
+    f (UseLemma lemma sub) =
+      simpleLemma $ certify $
+      eraseExcept (vars sub) $
+      subst sub $
+      derivation lemma
+    f p@UseAxiom{} = p
+    f p@Refl{} = p
+    f (Symm p) = Symm (f p)
+    f (Trans p q) = Trans (f p) (f q)
+    f (Cong fun ps) = Cong fun (map f ps)
+
+simplifyProof :: Function f => Config f -> [Derivation f] -> [Derivation f]
+simplifyProof config@Config{..} goals =
+  canonicaliseLemmas (fixpointOn key simp' (fixpointOn key simp goals))
+  where
+    simpCore =
+      (inlineUsedOnceLemmas `onlyIf` not cfg_all_lemmas) .
+      inlineTrivialLemmas config .
+      tightenProof
+
+    simp = simpCore . generaliseProof
+    -- generaliseProof undoes the effect of groundProof!
+    -- But we still want to run generaliseProof first, to simplify the proof
+    simp' = (simpCore . groundProof) `onlyIf` cfg_ground_proof
+
+    key ds =
+      (ds, [(equation p, derivation p) | p <- allLemmas ds])
+
+    pass `onlyIf` True = pass
+    _    `onlyIf` False = id
+
+simplificationPass ::
   Function f =>
-  Config -> [ProvedGoal f] -> [Proof f] -> Presentation f
-presentWithGoals config@Config{..} goals lemmas
-  -- We inline a lemma if one of the following holds:
+  -- A transformation on lemmas
+  (Map (Proof f) (Derivation f) -> Proof f -> Derivation f) ->
+  -- A transformation on goals
+  (Map (Proof f) (Derivation f) -> Derivation f -> Derivation f) ->
+  [Derivation f] -> [Derivation f]
+simplificationPass lemma goal p = map (op goal lem) p
+  where
+    lem = foldLemmas (op (\lem -> lemma lem . certify)) p
+    op f lem p =
+      f lem (unfoldLemmas (\lemma -> Just (lem Map.! lemma)) p)
+
+inlineTrivialLemmas :: Function f => Config f -> [Derivation f] -> [Derivation f]
+inlineTrivialLemmas Config{..} =
+  -- A lemma is trivial if one of the following holds:
   --   * It only has one step
   --   * It is subsumed by an earlier lemma
-  --   * It is only used once
   --   * It has to do with $equals (for printing of the goal proof)
   --   * The option cfg_no_lemmas is true
-  -- First we compute all inlinings, then apply simplify to remove them,
-  -- then repeat if any lemma was inlined
-  | Map.null inlinings =
-    let
-      axioms = usort $
-        concatMap (usedAxioms . derivation . pg_proof) goals ++
-        concatMap (usedAxioms . derivation) lemmas
-    in
-      Presentation axioms
-        (map flattenProof lemmas)
-        [ decodeGoal (goal { pg_proof = flattenProof pg_proof })
-        | goal@ProvedGoal{..} <- goals ]
+  simplificationPass inlineTrivial (const id)
+  where
+    inlineTrivial lem p
+      | shouldInline p = derivation p
+      | (q:_) <- subsuming lem (equation p) = q
+      | otherwise = simpleLemma p
 
-  | otherwise =
-    let
-      inline lemma = Map.lookup lemma inlinings
+    shouldInline p =
+      cfg_no_lemmas ||
+      length (filter (not . invisible) (map (equation . certify) (steps (derivation p)))) <= 1 ||
+      (not cfg_all_lemmas &&
+       (isJust (decodeEquality (eqn_lhs (equation p))) ||
+        isJust (decodeEquality (eqn_rhs (equation p)))))
 
-      goals' =
-        [ decodeGoal (goal { pg_proof = certify $ simplify inline (derivation pg_proof) })
-        | goal@ProvedGoal{..} <- goals ]
-      lemmas' =
-        [ certify $ simplify inline (derivation lemma)
-        | lemma <- lemmas, not (lemma `Map.member` inlinings) ]
-    in
-      presentWithGoals config goals' lemmas'
+    subsuming lem (t :=: u) =
+      subsuming1 lem (t :=: u) ++
+      map symm (subsuming1 lem (u :=: t))
+    subsuming1 lem eq =
+      [ subst sub d
+      | (q, d) <- Map.toList lem,
+        sub <- maybeToList (matchEquation (equation q) eq) ]
 
+inlineUsedOnceLemmas :: Function f => [Derivation f] -> [Derivation f]
+inlineUsedOnceLemmas ds =
+  -- Inline any lemma that's only used once in the proof
+  simplificationPass (const inlineOnce) (const id) ds
   where
-    inlinings =
-      Map.fromList
-        [ (lemma, p)
-        | lemma <- lemmas, Just p <- [tryInline lemma]]
+    uses = Map.unionsWith (+) $
+      map countUses ds ++ Map.elems (foldLemmas (const countUses) ds)
 
-    tryInline p
-      | shouldInline p = Just (derivation p)
-    tryInline p
-      -- Check for subsumption by an earlier lemma
-      | Just (m, q) <- Map.lookup (canonicalise (t :=: u)) equations, m < n =
-        Just (subsume p (derivation q))
-      | Just (m, q) <- Map.lookup (canonicalise (u :=: t)) equations, m < n =
-        Just (subsume p (Symm (derivation q)))
+    countUses p =
+      Map.fromListWith (+) (zip (usedLemmas p) (repeat (1 :: Int)))
+
+    inlineOnce p
+      | usedOnce p = derivation p
+      | otherwise = simpleLemma p
       where
-        t :=: u = equation p
-        Just (n, _) = Map.lookup (canonicalise (equation p)) equations
-    tryInline _ = Nothing
+        usedOnce p =
+          case Map.lookup p uses of
+            Just 1 -> True
+            _ -> False
 
-    shouldInline p =
-      cfg_no_lemmas ||
-      oneStep (derivation p) ||
-      (not cfg_all_lemmas &&
-       (isJust (decodeEquality (eqn_lhs (equation p))) ||
-        isJust (decodeEquality (eqn_rhs (equation p))) ||
-        Map.lookup p uses == Just 1))
-  
-    subsume p q =
-      -- Rename q so its variables match p's
-      subst sub q
+tightenProof :: Function f => [Derivation f] -> [Derivation f]
+tightenProof = mapLemmas tightenLemma
+  where
+    tightenLemma p =
+      fromSteps eq (map fst (fixpointOn length (tightenSteps eq) (zip ps eqs)))
       where
-        t  :=: u  = equation p
-        t' :=: u' = equation (certify q)
-        Just sub  = matchList (buildList [t', u']) (buildList [t, u])
+        eq = equation (certify p)
+        ps = steps p
+        eqs = map (equation . certify) ps
 
-    -- Record which lemma proves each equation
-    equations =
-      Map.fromList
-        [ (canonicalise (equation p), (i, p))
-        | (i, p) <- zip [0..] lemmas]
+    tightenSteps eq steps = head (cands ++ [steps])
+      where
+        -- Look for a segment of ps which can be removed, in the
+        -- sense that the terms at both ends of the segment are
+        -- unifiable without altering eq.
+        cands =
+          [ subst sub (before ++ after)
+          | (before, mid1) <- splits steps,
+            -- 'reverse' means we start with big segments.
+            (mid@(_:_), after) <- reverse (splits mid1),
+            let t :=: _ = snd (head mid)
+                _ :=: u = snd (last mid),
+            sub <- maybeToList (unify t u),
+            subst sub eq == eq ] ++
+          [ subst sub before
+          | (before, after@(_:_)) <- splits steps,
+            let t :=: _ = snd (head after)
+                _ :=: u = snd (last after),
+            sub <- maybeToList (match t u),
+            subst sub (eqn_lhs eq) == eqn_lhs eq ] ++
+          [ subst sub after
+          | (before@(_:_), after) <- reverse (splits steps),
+            let t :=: _ = snd (head before)
+                _ :=: u = snd (last before),
+            sub <- maybeToList (match u t),
+            subst sub (eqn_rhs eq) == eqn_rhs eq ]
 
-    -- Count how many times each lemma is used
-    uses =
-      Map.fromListWith (+)
-        [ (p, 1)
-        | p <-
-            concatMap usedLemmas
-              (map (derivation . pg_proof) goals ++
-               map derivation lemmas) ]
+generaliseProof :: Function f => [Derivation f] -> [Derivation f]
+generaliseProof =
+  simplificationPass (const generaliseLemma) (const generaliseGoal)
+  where
+    generaliseLemma p = lemma (certify q) sub
+      where
+        (q, sub) = generalise p
+    generaliseGoal p = subst sub q
+      where
+        (q, sub) = generalise (certify p)
 
-    -- Check if a proof only has one step.
-    -- Trans only occurs at the top level by this point.
-    oneStep Trans{} = False
-    oneStep _ = True
+    generalise p = (q, sub)
+      where
+        eq = equation p
+        n = freshVar eq
+        qs = evalState (mapM generaliseStep (steps (derivation p))) n
+        Just sub1 = unifyMany (stepsConstraints qs)
+        q = canonicalise (fromSteps eq (subst sub1 qs))
+        Just sub = matchEquation (equation (certify q)) eq
 
+    generaliseStep (UseAxiom axiom _) =
+      freshen (vars (axiom_eqn axiom)) (UseAxiom axiom)
+    generaliseStep (UseLemma lemma _) =
+      freshen (vars (equation lemma)) (UseLemma lemma)
+    generaliseStep (Refl _) = do
+      n <- get
+      put (n+1)
+      return (Refl (build (var (V n))))
+    generaliseStep (Symm p) =
+      Symm <$> generaliseStep p
+    generaliseStep (Trans p q) =
+      liftM2 Trans (generaliseStep p) (generaliseStep q)
+    generaliseStep (Cong f ps) =
+      Cong f <$> mapM generaliseStep ps
+
+    freshen xs f = do
+      n <- get
+      put (n + length xs)
+      let Just sub = listToSubst [(x, build (var (V i))) | (x, i) <- zip (usort xs) [n..]]
+      return (f sub)
+
+    stepsConstraints ps = zipWith combine eqs (tail eqs)
+      where
+        eqs = map (equation . certify) ps
+        combine (_ :=: t) (u :=: _) = (t, u)
+
+canonicaliseLemmas :: Function f => [Derivation f] -> [Derivation f]
+canonicaliseLemmas =
+  simplificationPass (const canonicaliseLemma) (const canonicalise)
+  where
+    -- Present the equation left-to-right, and with variables
+    -- named canonically
+    canonicaliseLemma p
+      | u `lessEqSkolem` t = canon (derivation p)
+      | otherwise = symm (canon (symm (derivation p)))
+      where
+        t :=: u = equation p
+        -- This ensures that we also renumber variables in the derivation that
+        -- do not occur in the equation, but that variables in the equation
+        -- get priority.
+        symbolic p = (equation p, derivation p)
+        before = symbolic p
+        after = canonicalise (symbolic p)
+        Just sub1 = matchManyList (terms before) (terms after)
+        Just sub2 = matchManyList (terms after) (terms before)
+        canon p = subst sub2 (simpleLemma (certify (subst sub1 p)))
+
 invisible :: Function f => Equation f -> Bool
 invisible (t :=: u) = show (pPrint t) == show (pPrint u)
 
 -- Pretty-print the proof of a single lemma.
-pPrintLemma :: Function f => Config -> (Axiom f -> String) -> (Proof f -> String) -> Proof f -> Doc
-pPrintLemma Config{..} axiomNum lemmaNum p =
-  ppTerm (eqn_lhs (equation q)) $$ pp (derivation q)
+pPrintLemma :: Function f => Config f -> (Axiom f -> String) -> (Proof f -> String) -> Proof f -> Doc
+pPrintLemma Config{..} axiomNum lemmaNum p
+  | null qs = text "Reflexivity."
+  | equation (certify (fromSteps (equation p) qs)) == equation p =
+    vcat (zipWith pp hl qs) $$ ppTerm (eqn_rhs (equation p))
+  | otherwise = error "lemma changed by pretty-printing!"
   where
-    q = flattenProof p
+    qs = steps (derivation p)
+    hl = map highlightStep qs
 
-    pp (Trans p q) = pp p $$ pp q
-    pp p | invisible (equation (certify p)) = pPrintEmpty
-    pp p =
-      (text "= { by" <+>
-       ppStep
-         (nub (map (show . ppLemma) (usedLemmasAndSubsts p)) ++
-          nub (map (show . ppAxiom) (usedAxiomsAndSubsts p))) <+>
-       text "}" $$
-       ppTerm (eqn_rhs (equation (certify p))))
+    pp _ p | invisible (equation (certify p)) = pPrintEmpty
+    pp h p =
+      ppTerm (HighlightedTerm [green | cfg_use_colour] (Just h) (eqn_lhs (equation (certify p)))) $$
+      text "=" <+> highlight [bold | cfg_use_colour] (text "{ by" <+> ppStep p <+> text "}")
 
+    highlightStep UseAxiom{} = []
+    highlightStep UseLemma{} = []
+    highlightStep (Symm p) = highlightStep p
+    highlightStep (Cong _ ps) = i:highlightStep p
+      where
+        [(i, p)] = filter (not . isRefl . snd) (zip [0..] ps)
+
     ppTerm t = text "  " <#> pPrint t
 
-    ppStep [] = text "reflexivity" -- ??
-    ppStep [x] = text x
-    ppStep xs =
-      hcat (punctuate (text ", ") (map text (init xs))) <+>
-      text "and" <+>
-      text (last xs)
+    ppStep = pp True
+      where
+        pp dir (UseAxiom axiom@Axiom{..} sub) =
+          text "axiom" <+> text (axiomNum axiom) <+> parens (text axiom_name) <+> ppDir dir <#> showSubst sub
+        pp dir (UseLemma lemma sub) =
+          text "lemma" <+> text (lemmaNum lemma) <+> ppDir dir <#> showSubst sub
+        pp dir (Symm p) =
+          pp (not dir) p
+        pp dir (Cong _ ps) = pp dir p
+          where
+            [p] = filter (not . isRefl) ps
 
-    ppLemma (p, sub) =
-      text "lemma" <+> text (lemmaNum p) <#> showSubst sub
-    ppAxiom (axiom@Axiom{..}, sub) =
-      text "axiom" <+> text (axiomNum axiom) <+> parens (text axiom_name) <#> showSubst sub
+    ppDir True = pPrintEmpty
+    ppDir False = text "R->L"
 
     showSubst sub
       | cfg_show_instances && not (null (substToList sub)) =
-        text " with " <#>
-        fsep (punctuate comma
-          [ pPrint x <+> text "->" <+> pPrint t
-          | (x, t) <- substToList sub ])
+        text " with " <#> pPrintSubst sub
       | otherwise = pPrintEmpty
 
--- Transform a proof so that each step uses exactly one axiom
--- or lemma. The proof will have the following form afterwards:
---   * Trans only occurs at the outermost level and is right-associated
---   * Each Cong has exactly one non-Refl argument (no parallel rewriting)
---   * Symm only occurs innermost, i.e., next to UseLemma or UseAxiom
---   * Refl only occurs as an argument to Cong, or outermost if the
---     whole proof is a single reflexivity step
-flattenProof :: Function f => Proof f -> Proof f
-flattenProof =
-  certify . flat . simplify (const Nothing) . derivation
-  where
-    flat (Trans p q) = trans (flat p) (flat q)
-    flat p@(Cong f ps) =
-      foldr trans (reflAfter p)
-        [ Cong f $
-            map reflAfter (take i ps) ++
-            [p] ++
-            map reflBefore (drop (i+1) ps)
-        | (i, q) <- zip [0..] qs,
-          p <- steps q ]
-      where
-        qs = map flat ps
-    flat p = p
-
-    reflBefore p = Refl (eqn_lhs (equation (certify p)))
-    reflAfter p  = Refl (eqn_rhs (equation (certify p)))
-
-    steps Refl{} = []
-    steps (Trans p q) = steps p ++ steps q
-    steps p = [p]
-
-    trans (Trans p q) r = trans p (trans q r)
-    trans Refl{} p = p
-    trans p Refl{} = p
-    trans p q =
-      case strip q of
-        Nothing -> Trans p q
-        Just q' -> trans p q'
-
-    strip p
-      | t == u = Just (Refl t)
-      | otherwise = strip' t p
-      where
-        t :=: u = equation (certify p)
-    strip' t (Trans _ q)
-      | eqn_lhs (equation (certify q)) == t = Just q
-      | otherwise = strip' t q
-    strip' _ _ = Nothing
+    isRefl Refl{} = True
+    isRefl _ = False
 
--- Transform a derivation into a list of single steps.
--- Each step has the following form:
---   * Trans does not occur
---   * Symm only occurs innermost, i.e., next to UseLemma or UseAxiom
---   * Each Cong has exactly one non-Refl argument (no parallel rewriting)
---   * Refl only occurs as an argument to Cong
-derivSteps :: Function f => Derivation f -> [Derivation f]
-derivSteps = steps . derivation . flattenProof . certify
-  where
-    steps Refl{} = []
-    steps (Trans p q) = steps p ++ steps q
-    steps p = [p]
+-- Pretty-print a substitution.
+pPrintSubst :: Function f => Subst f -> Doc
+pPrintSubst sub =
+  fsep (punctuate comma
+    [ pPrint x <+> text "->" <+> pPrint t
+    | (x, t) <- substToList sub ])
 
 -- | Print a presented proof.
-pPrintPresentation :: forall f. Function f => Config -> Presentation f -> Doc
+pPrintPresentation :: forall f. Function f => Config f -> Presentation f -> Doc
 pPrintPresentation config (Presentation axioms lemmas goals) =
   vcat $ intersperse (text "") $
-    vcat [ describeEquation "Axiom" (axiomNum axiom) (Just name) eqn
+    vcat [ describeEquation "Axiom" (axiomNum axiom) (Just name) eqn $$
+           ppAxiomUses axiom
          | axiom@(Axiom _ name eqn) <- axioms,
            not (invisible eqn) ]:
     [ pp "Lemma" (lemmaNum p) Nothing (equation p) emptySubst p
@@ -593,6 +829,24 @@
             else pPrintEmpty,
             text ""]
 
+    ppAxiomUses axiom
+      | cfg_show_uses_of_axioms config axiom && not (null uses) =
+        text "Used with:" $$
+        nest 2 (vcat
+          [ pPrint i <#> text "." <+> pPrintSubst sub
+          | (i, sub) <- zip [1 :: Int ..] uses ])
+      | otherwise = pPrintEmpty
+      where
+        uses = Set.toList (axiomUses axiom)
+
+    axiomUses axiom = Map.findWithDefault Set.empty axiom usesMap
+    usesMap =
+      Map.unionsWith Set.union
+        [ Map.map (Set.delete emptySubst . Set.map ground)
+            (groundAxiomsAndSubsts p)
+        | goal <- goals,
+          let p = derivation (pg_proof goal) ]
+
 -- | Format an equation nicely.
 --
 -- Used both here and in the main file.
@@ -662,9 +916,9 @@
 maybeDecodeGoal ProvedGoal{..}
   --  N.B. presentWithGoals takes care of expanding any lemma which mentions
   --  $equals, and flattening the proof.
-  | isFalseTerm u = extract (derivSteps deriv)
+  | isFalseTerm u = extract (steps deriv)
     -- Orient the equation so that $false is the RHS.
-  | isFalseTerm t = extract (derivSteps (symm deriv))
+  | isFalseTerm t = extract (steps (symm deriv))
   | otherwise = Nothing
   where
     isFalseTerm, isTrueTerm :: Term f -> Bool
diff --git a/Twee/Rule.hs b/Twee/Rule.hs
--- a/Twee/Rule.hs
+++ b/Twee/Rule.hs
@@ -12,8 +12,8 @@
 import Data.Maybe
 import Data.List
 import Twee.Utils
-import qualified Data.Set as Set
-import Data.Set(Set)
+import qualified Data.Map.Strict as Map
+import Data.Map(Map)
 import qualified Twee.Term as Term
 import Data.Ord
 import Twee.Equation
@@ -29,18 +29,33 @@
 data Rule f =
   Rule {
     -- | Information about whether and how the rule is oriented.
-    orientation :: !(Orientation f),
+    orientation :: Orientation f,
     -- Invariant:
     -- For oriented rules: vars rhs `isSubsetOf` vars lhs
     -- For unoriented rules: vars lhs == vars rhs
+    
+    -- | A proof that the rule holds.
+    rule_proof :: !(Proof f),
 
     -- | The left-hand side of the rule.
     lhs :: {-# UNPACK #-} !(Term f),
     -- | The right-hand side of the rule.
     rhs :: {-# UNPACK #-} !(Term f) }
-  deriving (Eq, Ord, Show)
+  deriving Show
+instance Eq (Rule f) where
+  x == y = compare x y == EQ
+instance Ord (Rule f) where
+  compare = comparing (\rule -> (lhs rule, rhs rule))
 type RuleOf a = Rule (ConstantOf a)
 
+ruleDerivation :: Rule f -> Derivation f
+ruleDerivation r =
+  case (matchEquation (Proof.equation (rule_proof r)) (lhs r :=: rhs r),
+        matchEquation (Proof.equation (rule_proof r)) (rhs r :=: lhs r)) of
+    (Just sub, _) -> Proof.lemma (rule_proof r) sub
+    (_, Just sub) -> Proof.symm (Proof.lemma (rule_proof r) sub)
+    _ -> error "rule out of sync with proof"
+
 -- | A rule's orientation.
 --
 -- 'Oriented' and 'WeaklyOriented' rules are used only left-to-right.
@@ -73,15 +88,10 @@
 oriented WeaklyOriented{} = True
 oriented _ = False
 
--- | Is a rule weakly oriented?
-weaklyOriented :: Orientation f -> Bool
-weaklyOriented WeaklyOriented{} = True
-weaklyOriented _ = False
-
 instance Symbolic (Rule f) where
   type ConstantOf (Rule f) = f
-  termsDL (Rule or t u) = termsDL or `mplus` termsDL t `mplus` termsDL u
-  subst_ sub (Rule or t u) = Rule (subst_ sub or) (subst_ sub t) (subst_ sub u)
+  termsDL (Rule _ _ t _) = termsDL t
+  subst_ sub (Rule or pf t u) = Rule (subst_ sub or) pf (subst_ sub t) (subst_ sub u)
 
 instance f ~ g => Has (Rule f) (Term g) where
   the = lhs
@@ -100,7 +110,7 @@
   subst_ _   Unoriented = Unoriented
 
 instance PrettyTerm f => Pretty (Rule f) where
-  pPrint (Rule or l r) =
+  pPrint (Rule or _ l r) =
     pPrint l <+> text (showOrientation or) <+> pPrint r
     where
       showOrientation Oriented = "->"
@@ -110,16 +120,15 @@
 
 -- | Turn a rule into an equation.
 unorient :: Rule f -> Equation f
-unorient (Rule _ l r) = l :=: r
+unorient (Rule _ _ l r) = l :=: r
 
 -- | Turn an equation t :=: u into a rule t -> u by computing the
 -- orientation info (e.g. oriented, permutative or unoriented).
 --
--- Crashes if t -> u is not a valid rule, for example if there is
--- a variable in @u@ which is not in @t@. To prevent this happening,
--- combine with 'Twee.CP.split'.
-orient :: Function f => Equation f -> Rule f
-orient (t :=: u) = Rule o t u
+-- Crashes if either @t < u@, or there is a variable in @u@ which is
+-- not in @t@. To avoid this problem, combine with 'Twee.CP.split'.
+orient :: Function f => Equation f -> Proof f -> Rule f
+orient (t :=: u) pf = Rule o pf t u
   where
     o | lessEq u t =
         case unify t u of
@@ -165,7 +174,7 @@
 
 -- | Flip an unoriented rule so that it goes right-to-left.
 backwards :: Rule f -> Rule f
-backwards (Rule or t u) = Rule (back or) u t
+backwards (Rule or pf t u) = Rule (back or) pf u t
   where
     back (Permutative xs) = Permutative (map swap xs)
     back Unoriented = Unoriented
@@ -177,12 +186,9 @@
 
 -- | Compute the normal form of a term wrt only oriented rules.
 {-# INLINEABLE simplify #-}
+{-# SCC simplify #-}
 simplify :: (Function f, Has a (Rule f)) => Index f a -> Term f -> Term f
-simplify !idx !t = {-# SCC simplify #-} simplify1 idx t
-
-{-# INLINEABLE simplify1 #-}
-simplify1 :: (Function f, Has a (Rule f)) => Index f a -> Term f -> Term f
-simplify1 idx t
+simplify !idx !t
   | t == u = t
   | otherwise = simplify idx u
   where
@@ -196,19 +202,9 @@
     simp (Cons (App f ts) us) =
       app f (simp ts) `mappend` simp us
 
--- | Check if a term can be simplified.
-{-# INLINEABLE canSimplify #-}
-canSimplify :: (Function f, Has a (Rule f)) => Index f a -> Term f -> Bool
-canSimplify idx t = canSimplifyList idx (singleton t)
-
-{-# INLINEABLE canSimplifyList #-}
-canSimplifyList :: (Function f, Has a (Rule f)) => Index f a -> TermList f -> Bool
-canSimplifyList idx t =
-  {-# SCC canSimplifyList #-}
-  any (isJust . simpleRewrite idx) (filter isApp (subtermsList t))
-
 -- | Find a simplification step that applies to a term.
 {-# INLINEABLE simpleRewrite #-}
+{-# SCC simpleRewrite #-}
 simpleRewrite :: (Function f, Has a (Rule f)) => Index f a -> Term f -> Maybe (Rule f, Subst f)
 simpleRewrite idx t =
   -- Use instead of maybeToList to make fusion work
@@ -226,193 +222,101 @@
 -- | A strategy gives a set of possible reductions for a term.
 type Strategy f = Term f -> [Reduction f]
 
--- | A multi-step rewrite proof @t ->* u@
-data Reduction f =
-    -- | Apply a single rewrite rule to the root of a term
-    Step {-# UNPACK #-} !(Proof f) !(Rule f) !(Subst f)
-    -- | Reflexivity
-  | Refl {-# UNPACK #-} !(Term f)
-    -- | Transivitity
-  | Trans !(Reduction f) !(Reduction f)
-    -- | Congruence
-  | Cong {-# UNPACK #-} !(Fun f) ![Reduction f]
-  deriving Show
-
-instance Symbolic (Reduction f) where
-  type ConstantOf (Reduction f) = f
-  termsDL (Step _ _ sub) = termsDL sub
-  termsDL (Refl t) = termsDL t
-  termsDL (Trans p q) = termsDL p `mplus` termsDL q
-  termsDL (Cong _ ps) = termsDL ps
-
-  subst_ sub (Step lemma rule s) = Step lemma rule (subst_ sub s)
-  subst_ sub (Refl t) = Refl (subst_ sub t)
-  subst_ sub (Trans p q) = Trans (subst_ sub p) (subst_ sub q)
-  subst_ sub (Cong f ps) = Cong f (subst_ sub ps)
-
-instance Function f => Pretty (Reduction f) where
-  pPrint = pPrint . reductionProof
+-- | A reduction proof is just a sequence of rewrite steps, stored
+-- as a list in reverse order. In each rewrite step, all subterms that
+-- are exactly equal to the LHS of the rule are replaced by the RHS,
+-- i.e. the rewrite step is performed as a parallel rewrite without
+-- matching.
+type Reduction f = [Rule f]
 
--- | A smart constructor for Trans which simplifies Refl.
+-- | Transitivity for reduction sequences.
 trans :: Reduction f -> Reduction f -> Reduction f
-trans Refl{} p = p
-trans p Refl{} = p
--- Make right-associative to improve performance of 'result'
-trans p (Trans q r) = Trans (Trans p q) r
-trans p q = Trans p q
-
--- | A smart constructor for Cong which simplifies Refl.
-cong :: Fun f -> [Reduction f] -> Reduction f
-cong f ps
-  | all isRefl ps = Refl (result (reduce (Cong f ps)))
-  | otherwise = Cong f ps
-  where
-    isRefl Refl{} = True
-    isRefl _ = False
+trans p q = q ++ p
 
--- | The list of all rewrite rules used in a rewrite proof.
-steps :: Reduction f -> [Reduction f]
-steps r = aux r []
+-- | Compute the final term resulting from a reduction, given the
+-- starting term.
+result :: Term f -> Reduction f -> Term f
+result t [] = t
+result t (r:rs) = ruleResult u r
   where
-    aux step@Step{} = (step:)
-    aux (Refl _) = id
-    aux (Trans p q) = aux p . aux q
-    aux (Cong _ ps) = foldr (.) id (map aux ps)
+    u = result t rs
 
 -- | Turn a reduction into a proof.
-reductionProof :: Reduction f -> Derivation f
-reductionProof (Step lemma _ sub) =
-  Proof.lemma lemma sub
-reductionProof (Refl t) = Proof.Refl t
-reductionProof (Trans p q) =
-  Proof.trans (reductionProof p) (reductionProof q)
-reductionProof (Cong f ps) = Proof.cong f (map reductionProof ps)
-
--- | Construct a basic rewrite step.
-{-# INLINE step #-}
-step :: (Has a (Rule f), Has a (Proof f)) => a -> Subst f -> Reduction f
-step x sub = Step (the x) (the x) sub
-
-----------------------------------------------------------------------
--- | A rewrite proof with the final term attached.
--- Has an @Ord@ instance which compares the final term.
-----------------------------------------------------------------------
-
-data Resulting f =
-  Resulting {
-    result :: {-# UNPACK #-} !(Term f),
-    reduction :: !(Reduction f) }
-  deriving Show
-
-instance Eq (Resulting f) where x == y = compare x y == EQ
-instance Ord (Resulting f) where compare = comparing result
-
-instance Symbolic (Resulting f) where
-  type ConstantOf (Resulting f) = f
-  termsDL (Resulting t red) =
-    termsDL t `mplus` termsDL red
-  subst_ sub (Resulting t red) =
-    Resulting (subst_ sub t) (subst_ sub red)
-
-instance Function f => Pretty (Resulting f) where
-  pPrint = pPrint . reduction
-
--- | Construct a 'Resulting' from a 'Reduction'.
-reduce :: Reduction f -> Resulting f
-reduce p =
-  Resulting (res p) p
+reductionProof :: PrettyTerm f => Term f -> Reduction f -> Derivation f
+reductionProof t ps = red t (Proof.Refl t) (reverse ps)
   where
-    res (Trans _ q) = res q
-    res (Refl t) = t
-    res p = {-# SCC res_emitRes #-} build (emitResult p)
+    red _ p [] = p
+    red t p (q:qs) =
+      red (ruleResult t q) (p `Proof.trans` ruleProof t q) qs
 
-    emitResult (Step _ r sub) = Term.subst sub (rhs r)
-    emitResult (Refl t) = builder t
-    emitResult (Trans _ q) = emitResult q
-    emitResult (Cong f ps) = app f (map emitResult ps)
+-- Helpers for result and reductionProof.
+ruleResult :: Term f -> Rule f -> Term f
+ruleResult t r = build (replace (lhs r) (rhs r) (singleton t))
 
+ruleProof :: PrettyTerm f => Term f -> Rule f -> Derivation f
+ruleProof t r@(Rule _ _ lhs _)
+  | t == lhs = ruleDerivation r
+  | len t < len lhs = Proof.Refl t
+ruleProof (App f ts) rule =
+  Proof.cong f [ruleProof u rule | u <- unpack ts]
+ruleProof t _ = Proof.Refl t
+
 --------------------------------------------------------------------------------
--- * Strategy combinators.
+-- * Normalisation.
 --------------------------------------------------------------------------------
 
 -- | Normalise a term wrt a particular strategy.
 {-# INLINE normaliseWith #-}
-normaliseWith :: Function f => (Term f -> Bool) -> Strategy f -> Term f -> Resulting f
-normaliseWith ok strat t = {-# SCC normaliseWith #-} res
+{-# SCC normaliseWith #-}
+normaliseWith :: Function f => (Term f -> Bool) -> Strategy f -> Term f -> Reduction f
+normaliseWith ok strat t = res
   where
-    res = aux 0 (Refl t) t
+    res = aux 0 [] t
     aux 1000 p _ =
       error $
         "Possibly nonterminating rewrite:\n" ++ prettyShow p
     aux n p t =
-      case parallel strat t of
-        (q:_) | u <- result (reduce q), ok u ->
+      case anywhere strat t of
+        (q:_) | u <- result t q, ok u ->
           aux (n+1) (p `trans` q) u
-        _ -> Resulting t p
+        _ -> p
 
 -- | Compute all normal forms of a set of terms wrt a particular strategy.
-normalForms :: Function f => Strategy f -> Set (Resulting f) -> Set (Resulting f)
+normalForms :: Function f => Strategy f -> Map (Term f) (Reduction f) -> Map (Term f) (Term f, Reduction f)
 normalForms strat ps = snd (successorsAndNormalForms strat ps)
 
 -- | Compute all successors of a set of terms (a successor of a term @t@
 -- is a term @u@ such that @t ->* u@).
-successors :: Function f => Strategy f -> Set (Resulting f) -> Set (Resulting f)
+successors :: Function f => Strategy f -> Map (Term f) (Reduction f) -> Map (Term f) (Term f, Reduction f)
 successors strat ps =
-  Set.union qs rs
+  Map.union qs rs
   where
     (qs, rs) = successorsAndNormalForms strat ps
 
 {-# INLINEABLE successorsAndNormalForms #-}
-successorsAndNormalForms :: Function f => Strategy f -> Set (Resulting f) ->
-  (Set (Resulting f), Set (Resulting f))
+{-# SCC successorsAndNormalForms #-}
+successorsAndNormalForms :: Function f => Strategy f -> Map (Term f) (Reduction f) ->
+  (Map (Term f) (Term f, Reduction f), Map (Term f) (Term f, Reduction f))
 successorsAndNormalForms strat ps =
-  {-# SCC successorsAndNormalForms #-} go Set.empty Set.empty ps
+  go Map.empty Map.empty (Map.mapWithKey (\t red -> (t, red)) ps)
   where
     go dead norm ps =
-      case Set.minView ps of
+      case Map.minViewWithKey ps of
         Nothing -> (dead, norm)
-        Just (p, ps)
-          | p `Set.member` dead -> go dead norm ps
-          | p `Set.member` norm -> go dead norm ps
-          | null qs -> go dead (Set.insert p norm) ps
+        Just ((t, p), ps)
+          | t `Map.member` dead -> go dead norm ps
+          | t `Map.member` norm -> go dead norm ps
+          | null qs -> go dead (Map.insert t p norm) ps
           | otherwise ->
-            go (Set.insert p dead) norm (Set.fromList qs `Set.union` ps)
+            go (Map.insert t p dead) norm (Map.fromList qs `Map.union` ps)
           where
             qs =
-              [ reduce (reduction p `trans` q)
-              | q <- anywhere strat (result p) ]
+              [ (result t q, (fst p, (snd p `trans` q)))
+              | q <- anywhere strat t ]
 
 -- | Apply a strategy anywhere in a term.
 anywhere :: Strategy f -> Strategy f
-anywhere strat t = strat t ++ nested (anywhere strat) t
-
--- | Apply a strategy to some child of the root function.
-nested :: Strategy f -> Strategy f
-nested _ Var{} = []
-nested strat (App f ts) =
-  cong f <$> inner [] ts
-  where
-    inner _ Empty = []
-    inner before (Cons t u) =
-      [ reverse before ++ [p] ++ map Refl (unpack u)
-      | p <- strat t ] ++
-      inner (Refl t:before) u
-
--- | Apply a strategy in parallel in as many places as possible.
--- Takes only the first rewrite of each strategy.
-{-# INLINE parallel #-}
-parallel :: PrettyTerm f => Strategy f -> Strategy f
-parallel strat t =
-  case par t of
-    Refl{} -> []
-    p -> [p]
-  where
-    par t | p:_ <- strat t = p
-    par (App f ts) = cong f (inner [] ts)
-    par t = Refl t
-
-    inner before Empty = reverse before
-    inner before (Cons t u) = inner (par t:before) u
+anywhere strat t = concatMap strat (subterms t)
 
 --------------------------------------------------------------------------------
 -- * Basic strategies. These only apply at the root of the term.
@@ -420,24 +324,24 @@
 
 -- | A strategy which rewrites using an index.
 {-# INLINE rewrite #-}
-rewrite :: (Function f, Has a (Rule f), Has a (Proof f)) => (Rule f -> Subst f -> Bool) -> Index f a -> Strategy f
+rewrite :: (Function f, Has a (Rule f)) => (Rule f -> Subst f -> Bool) -> Index f a -> Strategy f
 rewrite p rules t = do
   rule <- Index.approxMatches t rules
   tryRule p rule t
 
 -- | A strategy which applies one rule only.
 {-# INLINEABLE tryRule #-}
-tryRule :: (Function f, Has a (Rule f), Has a (Proof f)) => (Rule f -> Subst f -> Bool) -> a -> Strategy f
+tryRule :: (Function f, Has a (Rule f)) => (Rule f -> Subst f -> Bool) -> a -> Strategy f
 tryRule p rule t = do
   sub <- maybeToList (match (lhs (the rule)) t)
   guard (p (the rule) sub)
-  return (step rule sub)
+  return [subst sub (the rule)]
 
 -- | Check if a rule can be applied, given an ordering <= on terms.
 {-# INLINEABLE reducesWith #-}
 reducesWith :: Function f => (Term f -> Term f -> Bool) -> Rule f -> Subst f -> Bool
-reducesWith _ (Rule Oriented _ _) _ = True
-reducesWith _ (Rule (WeaklyOriented min ts) _ _) sub =
+reducesWith _ (Rule Oriented _ _ _) _ = True
+reducesWith _ (Rule (WeaklyOriented min ts) _ _ _) sub =
   -- Be a bit careful here not to build new terms
   -- (reducesWith is used in simplify).
   -- This is the same as:
@@ -449,7 +353,7 @@
 
     isMinimal (App f Empty) = f == min
     isMinimal _ = False
-reducesWith p (Rule (Permutative ts) _ _) sub =
+reducesWith p (Rule (Permutative ts) _ _ _) sub =
   aux ts
   where
     aux [] = False
@@ -459,8 +363,8 @@
       where
         t' = subst sub t
         u' = subst sub u
-reducesWith p (Rule Unoriented t u) sub =
-  p u' t' && u' /= t'
+reducesWith p (Rule Unoriented _ t u) sub =
+  t' /= u' && p u' t'
   where
     t' = subst sub t
     u' = subst sub u
@@ -486,9 +390,4 @@
 {-# INLINEABLE reducesSkolem #-}
 reducesSkolem :: Function f => Rule f -> Subst f -> Bool
 reducesSkolem rule sub =
-  reducesWith (\t u -> lessEq (subst skolemise t) (subst skolemise u)) rule sub
-  where
-    skolemise (V x) = con (skolem (V (x + k)))
-    -- Make sure the Skolem constants we choose don't overlap with any
-    -- already in the rule
-    V k = maximum (V 0:map succ (catMaybes (map getSkolem (funs rule))))
+  reducesWith (\t u -> lessEqSkolem t u) rule sub
diff --git a/Twee/Rule/Index.hs b/Twee/Rule/Index.hs
--- a/Twee/Rule/Index.hs
+++ b/Twee/Rule/Index.hs
@@ -13,21 +13,19 @@
 data RuleIndex f a =
   RuleIndex {
     index_oriented :: !(Index f a),
-    index_weak     :: !(Index f a),
     index_all      :: !(Index f a) }
   deriving Show
 
 empty :: RuleIndex f a
-empty = RuleIndex Index.empty Index.empty Index.empty
+empty = RuleIndex Index.empty Index.empty
 
 insert :: forall f a. Has a (Rule f) => Term f -> a -> RuleIndex f a -> RuleIndex f a
 insert t x RuleIndex{..} =
   RuleIndex {
     index_oriented = insertWhen (oriented or) index_oriented,
-    index_weak = insertWhen (weaklyOriented or) index_weak,
     index_all = insertWhen True index_all }
   where
-    Rule or _ _ = the x :: Rule f
+    Rule or _ _ _ = the x :: Rule f
 
     insertWhen False idx = idx
     insertWhen True idx = Index.insert t x idx
@@ -36,10 +34,9 @@
 delete t x RuleIndex{..} =
   RuleIndex {
     index_oriented = deleteWhen (oriented or) index_oriented,
-    index_weak = deleteWhen (weaklyOriented or) index_weak,
     index_all = deleteWhen True index_all }
   where
-    Rule or _ _ = the x :: Rule f
+    Rule or _ _ _ = the x :: Rule f
 
     deleteWhen False idx = idx
     deleteWhen True idx = Index.delete t x idx
diff --git a/Twee/Term.hs b/Twee/Term.hs
--- a/Twee/Term.hs
+++ b/Twee/Term.hs
@@ -22,11 +22,11 @@
   -- * Terms
   Term, pattern Var, pattern App, isApp, isVar, singleton, len,
   -- * Termlists
-  TermList, pattern Empty, pattern Cons, pattern ConsSym,
-  pattern UnsafeCons, pattern UnsafeConsSym,
+  TermList, pattern Empty, pattern Cons, pattern ConsSym, hd, tl, rest,
+  pattern UnsafeCons, pattern UnsafeConsSym, uhd, utl, urest,
   empty, unpack, lenList,
   -- * Function symbols and variables
-  Fun, fun, fun_id, fun_value, pattern F, Var(..), 
+  Fun, fun, fun_id, fun_value, pattern F, Var(..), Labelled(..),
   -- * Building terms
   Build(..),
   Builder,
@@ -46,29 +46,35 @@
   -- ** Other operations on substitutions
   foldSubst, allSubst, substDomain,
   substSize,
-  substCompose, substCompatible, substUnion, idempotent, idempotentOn,
+  substCompatible, substUnion, idempotent, idempotentOn,
   canonicalise,
   -- * Matching
-  match, matchIn, matchList, matchListIn, isInstanceOf, isVariantOf,
+  match, matchIn, matchList, matchListIn,
+  matchMany, matchManyIn, matchManyList, matchManyListIn,
+  isInstanceOf, isVariantOf,
   -- * Unification
-  unify, unifyList,
-  unifyTri, unifyListTri, unifyListTriFrom,
-  TriangleSubst(..),
+  unify, unifyList, unifyMany,
+  unifyTri, unifyTriFrom, unifyListTri, unifyListTriFrom,
+  TriangleSubst(..), emptyTriangleSubst,
   close,
   -- * Positions in terms
   positionToPath, pathToPosition,
   replacePosition,
   replacePositionSub,
+  replace,
   -- * Miscellaneous functions
   bound, boundList, boundLists, mapFun, mapFunList, (<<)) where
 
 import Prelude hiding (lookup)
 import Twee.Term.Core hiding (F)
+import qualified Twee.Term.Core as Core
 import Data.List hiding (lookup, find)
 import Data.Maybe
 import Data.Semigroup(Semigroup(..))
 import Data.IntMap.Strict(IntMap)
 import qualified Data.IntMap.Strict as IntMap
+import Control.Arrow((&&&))
+import Twee.Utils
 
 --------------------------------------------------------------------------------
 -- * A type class for builders
@@ -109,8 +115,9 @@
 
 -- | Build a termlist.
 {-# INLINE buildList #-}
+{-# SCC buildList #-}
 buildList :: Build a => a -> TermList (BuildFun a)
-buildList x = {-# SCC buildList #-} buildTermList (builder x)
+buildList x = buildTermList (builder x)
 
 -- | Build a constant (a function with no arguments).
 {-# INLINE con #-}
@@ -202,7 +209,7 @@
 newtype Subst f =
   Subst {
     unSubst :: IntMap (TermList f) }
-  deriving Eq
+  deriving (Eq, Ord)
 
 -- | Return the highest-number variable in a substitution plus 1.
 {-# INLINE substSize #-}
@@ -237,11 +244,6 @@
 unsafeExtendList :: Var -> TermList f -> Subst f -> Subst f
 unsafeExtendList x !t (Subst sub) = Subst (IntMap.insert (var_id x) t sub)
 
--- | Compose two substitutions.
-substCompose :: Substitution s => Subst (SubstFun s) -> s -> Subst (SubstFun s)
-substCompose (Subst !sub1) !sub2 =
-  Subst (IntMap.map (buildList . substList sub2) sub1)
-
 -- | Check if two substitutions are compatible (they do not send the same
 -- variable to different terms).
 substCompatible :: Subst f -> Subst f -> Bool
@@ -272,14 +274,17 @@
 idempotentOn !sub = aux
   where
     aux Empty = True
-    aux (ConsSym App{} t) = aux t
+    aux ConsSym{hd = App{}, rest = t} = aux t
     aux (Cons (Var x) t) = isNothing (lookupList x sub) && aux t
 
 -- | Iterate a triangle substitution to make it idempotent.
 close :: TriangleSubst f -> Subst f
 close (Triangle sub)
   | idempotent sub = sub
-  | otherwise      = close (Triangle (substCompose sub sub))
+  | otherwise      = close (Triangle (compose sub sub))
+  where
+    compose (Subst !sub1) !sub2 =
+      Subst (IntMap.map (buildList . substList sub2) sub1)
 
 -- | Return a substitution which renames the variables of a list of terms to put
 -- them in a canonical order.
@@ -298,16 +303,20 @@
     loop sub _ Empty [] = sub
     loop sub Empty _ _ = sub
     loop sub vs Empty (t:ts) = loop sub vs t ts
-    loop sub vs (ConsSym App{} t) ts = loop sub vs t ts
+    loop sub vs ConsSym{hd = App{}, rest = t} ts = loop sub vs t ts
     loop sub vs0@(Cons v vs) (Cons (Var x) t) ts =
       case extend x v sub of
         Just sub -> loop sub vs  t ts
         Nothing  -> loop sub vs0 t ts
 
 -- | The empty substitution.
-{-# NOINLINE emptySubst #-}
+emptySubst :: Subst f
 emptySubst = Subst IntMap.empty
 
+-- | The empty triangle substitution.
+emptyTriangleSubst :: TriangleSubst f
+emptyTriangleSubst = Triangle emptySubst
+
 -- | Construct a substitution from a list.
 -- Returns @Nothing@ if a variable is bound to several different terms.
 listToSubst :: [(Var, Term f)] -> Maybe (Subst f)
@@ -337,20 +346,47 @@
 
 -- | A variant of 'match' which works on termlists
 -- and extends an existing substitution.
+{-# SCC matchListIn #-}
 matchListIn :: Subst f -> TermList f -> TermList f -> Maybe (Subst f)
 matchListIn !sub !pat !t
   | lenList t < lenList pat = Nothing
   | otherwise =
-    let loop !_ !_ !_ | False = undefined
-        loop sub Empty Empty = Just sub
-        loop sub (ConsSym (App f _) pat) (ConsSym (App g _) t)
-          | f == g = loop sub pat t
-        loop sub (Cons (Var x) pat) (Cons t u) = do
-          sub <- extend x t sub
-          loop sub pat u
+    let 
+        loop !sub ConsSym{hd = pat, tl = pats, rest = pats1} !ts = do
+          ConsSym{hd = t, tl = ts, rest = ts1} <- Just ts
+          case (pat, t) of
+            (App f _, App g _) | f == g ->
+              loop sub pats1 ts1
+            (Var x, _) -> do
+              sub <- extend x t sub
+              loop sub pats ts
+            _ -> Nothing
+        loop sub _ Empty = Just sub
         loop _ _ _ = Nothing
-    in {-# SCC match #-} loop sub pat t
+    in loop sub pat t
 
+-- | A variant of 'match' which works on lists of terms.
+matchMany :: [Term f] -> [Term f] -> Maybe (Subst f)
+matchMany pat t = matchManyIn emptySubst pat t
+
+-- | A variant of 'match' which works on lists of terms,
+-- and extends an existing substitution.
+matchManyIn :: Subst f -> [Term f] -> [Term f] -> Maybe (Subst f)
+matchManyIn sub ts us = matchManyListIn sub (map singleton ts) (map singleton us)
+
+-- | A variant of 'match' which works on lists of termlists.
+matchManyList :: [TermList f] -> [TermList f] -> Maybe (Subst f)
+matchManyList pat t = matchManyListIn emptySubst pat t
+
+-- | A variant of 'match' which works on lists of termlists,
+-- and extends an existing substitution.
+matchManyListIn :: Subst f -> [TermList f] -> [TermList f] -> Maybe (Subst f)
+matchManyListIn !sub [] [] = return sub
+matchManyListIn sub (t:ts) (u:us) = do
+  sub <- matchListIn sub t u
+  matchManyListIn sub ts us
+matchManyListIn _ _ _ = Nothing
+
 --------------------------------------------------------------------------------
 -- Unification.
 --------------------------------------------------------------------------------
@@ -397,32 +433,49 @@
   -- Not strict so that isJust (unify t u) doesn't force the substitution
   return (close sub)
 
+-- | Unify a collection of pairs of terms.
+unifyMany :: [(Term f, Term f)] -> Maybe (Subst f)
+unifyMany ts = unifyList us vs
+  where
+    us = buildList (map fst ts)
+    vs = buildList (map snd ts)
+
 -- | Unify two terms, returning a triangle substitution.
 -- This is slightly faster than 'unify'.
 unifyTri :: Term f -> Term f -> Maybe (TriangleSubst f)
 unifyTri t u = unifyListTri (singleton t) (singleton u)
 
+-- | Unify two terms, starting from an existing substitution.
+unifyTriFrom :: Term f -> Term f -> TriangleSubst f -> Maybe (TriangleSubst f)
+unifyTriFrom t u sub = unifyListTriFrom (singleton t) (singleton u) sub
+
 -- | Unify two termlists, returning a triangle substitution.
 -- This is slightly faster than 'unify'.
 unifyListTri :: TermList f -> TermList f -> Maybe (TriangleSubst f)
 unifyListTri t u = unifyListTriFrom t u (Triangle emptySubst)
 
+{-# SCC unifyListTriFrom #-}
 unifyListTriFrom :: TermList f -> TermList f -> TriangleSubst f -> Maybe (TriangleSubst f)
 unifyListTriFrom !t !u (Triangle !sub) =
-  fmap Triangle ({-# SCC unify #-} loop sub t u)
+  fmap Triangle (loop sub t u)
   where
     loop !_ !_ !_ | False = undefined
-    loop sub Empty Empty = Just sub
-    loop sub (ConsSym (App f _) t) (ConsSym (App g _) u)
-      | f == g = loop sub t u
-    loop sub (Cons (Var x) t) (Cons u v) = do
-      sub <- var sub x u
-      loop sub t v
-    loop sub (Cons t u) (Cons (Var x) v) = do
-      sub <- var sub x t
-      loop sub u v
+    loop sub (ConsSym{hd = t, tl = ts, rest = ts1}) u = do
+      ConsSym{hd = u, tl = us, rest =  us1} <- Just u
+      case (t, u) of
+        (App f _, App g _) | f == g ->
+          loop sub ts1 us1
+        (Var x, _) -> do
+          sub <- var sub x u
+          loop sub ts us
+        (_, Var x) -> do
+          sub <- var sub x t
+          loop sub ts us
+        _ -> Nothing
+    loop sub _ Empty = Just sub
     loop _ _ _ = Nothing
 
+    {-# INLINE var #-}
     var sub x t =
       case lookupList x sub of
         Just u -> loop sub u (singleton t)
@@ -439,15 +492,17 @@
       occurs sub x (singleton t)
       extend x t sub
 
-    occurs !_ !_ Empty = Just ()
-    occurs sub x (ConsSym App{} t) = occurs sub x t
-    occurs sub x (ConsSym (Var y) t)
-      | x == y = Nothing
-      | otherwise = do
-          occurs sub x t
-          case lookupList y sub of
-            Nothing -> Just ()
-            Just u  -> occurs sub x u
+    occurs !sub !x (ConsSym{hd = t, rest = ts}) =
+      case t of
+        App{} -> occurs sub x ts
+        Var y
+          | x == y -> Nothing
+          | otherwise -> do
+            occurs sub x ts
+            case lookupList y sub of
+              Nothing -> Just ()
+              Just u  -> occurs sub x u
+    occurs _ _ _ = Just ()
 
 --------------------------------------------------------------------------------
 -- Miscellaneous stuff.
@@ -462,7 +517,7 @@
 children :: Term f -> TermList f
 children t =
   case singleton t of
-    UnsafeConsSym _ ts -> ts
+    UnsafeConsSym{urest = ts} -> ts
 
 -- | Convert a termlist into an ordinary list of terms.
 unpack :: TermList f -> [Term f]
@@ -471,15 +526,15 @@
     op Empty = Nothing
     op (Cons t ts) = Just (t, ts)
 
-instance Show (Term f) where
+instance (Labelled f, Show f) => Show (Term f) where
   show (Var x) = show x
   show (App f Empty) = show f
   show (App f ts) = show f ++ "(" ++ intercalate "," (map show (unpack ts)) ++ ")"
 
-instance Show (TermList f) where
+instance (Labelled f, Show f) => Show (TermList f) where
   show = show . unpack
 
-instance Show (Subst f) where
+instance (Labelled f, Show f) => Show (Subst f) where
   show subst =
     show
       [ (i, t)
@@ -509,27 +564,19 @@
 bound t = boundList (singleton t)
 
 -- | Return the lowest- and highest-numbered variables in a termlist.
-{-# INLINE boundList #-}
 boundList :: TermList f -> (Var, Var)
-boundList t = boundListFrom (V maxBound) (V minBound) t
-
-boundListFrom :: Var -> Var -> TermList f -> (Var, Var)
-boundListFrom !m !n Empty = (m, n)
-boundListFrom m n (ConsSym App{} t) = boundListFrom m n t
-boundListFrom m n (ConsSym (Var x) t) =
-  boundListFrom (m `min` x) (n `max` x) t
+boundList t = boundListFrom (V maxBound, V minBound) t
 
 -- | Return the lowest- and highest-numbered variables in a list of termlists.
 boundLists :: [TermList f] -> (Var, Var)
-boundLists t = boundListsFrom (V maxBound) (V minBound) t
+boundLists ts = foldl' boundListFrom (V maxBound, V minBound) ts
 
-boundListsFrom :: Var -> Var -> [TermList f] -> (Var, Var)
-boundListsFrom !m !n [] = (m, n)
-boundListsFrom m n (t:ts) =
-  let
-    (m', n') = boundListFrom m n t
-  in
-    boundListsFrom m' n' ts
+{-# INLINE boundListFrom #-}
+boundListFrom :: (Var, Var) -> TermList f -> (Var, Var)
+boundListFrom (V !ex, V !ey) ts = (V x, V y)
+  where
+    !(!x, !y) = foldl' op (ex, ey) [x | Var (V x) <- subtermsList ts]
+    op (!mn, !mx) x = (mn `intMin` x, mx `intMax` x)
 
 -- | Check if a variable occurs in a term.
 {-# INLINE occurs #-}
@@ -542,7 +589,7 @@
 subtermsList t = unfoldr op t
   where
     op Empty = Nothing
-    op (ConsSym t u) = Just (t, u)
+    op ConsSym{hd = t, rest = u} = Just (t, u)
 
 -- | Find all subterms of a term.
 {-# INLINE subterms #-}
@@ -588,6 +635,16 @@
     aux (Cons (Var x) ts) = var x `mappend` aux ts
     aux (Cons (App ff ts) us) = app (f ff) (aux ts) `mappend` aux us
 
+{-# INLINE replace #-}
+replace :: (Build a, BuildFun a ~ f) => Term f -> a -> TermList f -> Builder f
+replace !_ !_ Empty = mempty
+replace t u (Cons v vs)
+  | t == v = builder u `mappend` replace t u vs
+  | otherwise =
+    case v of
+      Var x -> var x `mappend` replace t u vs
+      App f ts -> app f (replace t u ts) `mappend` replace t u vs
+
 -- | Replace the term at a given position in a term with a different term.
 {-# INLINE replacePosition #-}
 replacePosition :: (Build a, BuildFun a ~ f) => Int -> a -> TermList f -> Builder f
@@ -645,11 +702,28 @@
     list k (Cons t u) n ns =
       list (k+len t) u (n-1) ns
 
+class Labelled f where
+  -- | Labels should be small positive integers!
+  label :: f -> Int
+  find :: Int -> f
+
+instance (Labelled f, Show f) => Show (Fun f) where show = show . fun_value
+
 -- | A pattern which extracts the 'fun_value' from a 'Fun'.
-pattern F :: f -> Fun f
-pattern F x <- (fun_value -> x)
+pattern F :: Labelled f => Int -> f -> Fun f
+pattern F x y <- (fun_id &&& fun_value -> (x, y))
 {-# COMPLETE F #-}
 
 -- | Compare the 'fun_value's of two 'Fun's.
-(<<) :: Ord f => Fun f -> Fun f -> Bool
+(<<) :: (Labelled f, Ord f) => Fun f -> Fun f -> Bool
 f << g = fun_value f < fun_value g
+
+-- | Construct a 'Fun' from a function symbol.
+{-# INLINEABLE fun #-}
+fun :: Labelled f => f -> Fun f
+fun f = Core.F (fromIntegral (label f))
+
+-- | The underlying function symbol of a 'Fun'.
+{-# INLINEABLE fun_value #-}
+fun_value :: Labelled f => Fun f -> f
+fun_value x = find (fun_id x)
diff --git a/Twee/Term/Core.hs b/Twee/Term/Core.hs
--- a/Twee/Term/Core.hs
+++ b/Twee/Term/Core.hs
@@ -23,8 +23,6 @@
 import GHC.Prim
 import GHC.ST hiding (liftST)
 import Data.Ord
-import Twee.Label
-import Data.Typeable
 import Data.Semigroup(Semigroup(..))
 
 --------------------------------------------------------------------------------
@@ -42,7 +40,7 @@
 
 instance Show Symbol where
   show Symbol{..}
-    | isFun = show (F index) ++ "=" ++ show size
+    | isFun = "f" ++ show index ++ "=" ++ show size
     | otherwise = show (V index)
 
 -- Convert symbols to/from Int64 for storage in flatterms.
@@ -133,7 +131,7 @@
 -- | Like 'Cons', but does not check that the termlist is non-empty. Use only if
 -- you are sure the termlist is non-empty.
 pattern UnsafeCons :: Term f -> TermList f -> TermList f
-pattern UnsafeCons t ts <- (unsafePatHead -> Just (t, _, ts))
+pattern UnsafeCons t ts <- (unsafePatHead -> (t, _, ts))
 
 -- | Matches a non-empty termlist, unpacking it into head and
 -- /everything except the root symbol of the head/.
@@ -144,21 +142,21 @@
 --
 -- > u  = f(x,y)
 -- > us = [x, y, g(z)]
-pattern ConsSym :: Term f -> TermList f -> TermList f
-pattern ConsSym t ts <- (patHead -> Just (t, ts, _))
+pattern ConsSym :: Term f -> TermList f -> TermList f -> TermList f
+pattern ConsSym{hd, tl, rest} <- (patHead -> Just (hd, rest, tl))
 
 -- | Like 'ConsSym', but does not check that the termlist is non-empty. Use only
 -- if you are sure the termlist is non-empty.
-pattern UnsafeConsSym :: Term f -> TermList f -> TermList f
-pattern UnsafeConsSym t ts <- (unsafePatHead -> Just (t, ts, _))
+pattern UnsafeConsSym :: Term f -> TermList f -> TermList f -> TermList f
+pattern UnsafeConsSym{uhd, utl, urest} <- (unsafePatHead -> (uhd, urest, utl))
 
 -- A helper for UnsafeCons/UnsafeConsSym.
 {-# INLINE unsafePatHead #-}
-unsafePatHead :: TermList f -> Maybe (Term f, TermList f, TermList f)
+unsafePatHead :: TermList f -> (Term f, TermList f, TermList f)
 unsafePatHead TermList{..} =
-  Just (Term x (TermList low (low+size) array),
-        TermList (low+1) high array,
-        TermList (low+size) high array)
+  (Term x (TermList low (low+size) array),
+   TermList (low+1) high array,
+   TermList (low+size) high array)
   where
     !x = indexByteArray array low
     Symbol{..} = toSymbol x
@@ -168,7 +166,7 @@
 patHead :: TermList f -> Maybe (Term f, TermList f, TermList f)
 patHead t@TermList{..}
   | low == high = Nothing
-  | otherwise = unsafePatHead t
+  | otherwise = Just (unsafePatHead t)
 
 -- Pattern synonyms for single terms.
 -- * Var :: Var -> Term f
@@ -178,20 +176,9 @@
 -- by the user; @'Fun' f@ is an @f@ together with an automatically-generated unique number.
 newtype Fun f =
   F {
-    -- | The unique number of a 'Fun'.
+    -- | The unique number of a 'Fun'. Must fit in 32 bits.
     fun_id :: Int }
-instance Eq (Fun f) where
-  f == g = fun_id f == fun_id g
-instance Ord (Fun f) where
-  compare = comparing fun_id
-
--- | Construct a 'Fun' from a function symbol.
-fun :: (Ord f, Typeable f) => f -> Fun f
-fun f = F (fromIntegral (labelNum (label f)))
-
--- | The underlying function symbol of a 'Fun'.
-fun_value :: Fun f -> f
-fun_value f = find (unsafeMkLabel (fromIntegral (fun_id f)))
+  deriving (Eq, Ord)
 
 -- | A variable.
 newtype Var =
@@ -200,8 +187,8 @@
     -- Don't use huge variable numbers:
     -- they will be truncated to 32 bits when stored in a term.
     var_id :: Int } deriving (Eq, Ord, Enum)
-instance Show (Fun f) where show f = "f" ++ show (fun_id f)
-instance Show Var     where show x = "x" ++ show (var_id x)
+instance Show Var where
+  show x = "x" ++ show (var_id x)
 
 -- | Matches a variable.
 pattern Var :: Var -> Term f
@@ -216,60 +203,29 @@
 -- A helper function for Var and App.
 {-# INLINE patTerm #-}
 patTerm :: Term f -> Either Var (Fun f, TermList f)
-patTerm t@Term{..}
+patTerm Term{..}
   | isFun     = Right (F index, ts)
   | otherwise = Left (V index)
   where
     Symbol{..} = toSymbol root
-    !(UnsafeConsSym _ ts) = singleton t
+    !UnsafeConsSym{urest = ts} = termlist
 
 -- | Convert a term to a termlist.
 {-# INLINE singleton #-}
 singleton :: Term f -> TermList f
 singleton Term{..} = termlist
 
--- We can implement equality almost without access to the
--- internal representation of the termlists, but we cheat by
--- comparing Int64s instead of Symbols.
 instance Eq (TermList f) where
-  -- Manual worker-wrapper to prevent too much from being inlined.
-  t == u = eqTermList t u
-
-{-# INLINE eqTermList #-}
-eqTermList :: TermList f -> TermList f -> Bool
-eqTermList
-  (TermList (I# low1) (I# high1) (ByteArray array1))
-  (TermList (I# low2) (I# high2) (ByteArray array2)) =
-    weqTermList low1 high1 array1 low2 high2 array2
-
--- Manually worker-wrapper transform the thing, ugh...
-{-# NOINLINE weqTermList #-}
-weqTermList ::
-  Int# -> Int# -> ByteArray# ->
-  Int# -> Int# -> ByteArray# ->
-  Bool
-weqTermList low1 high1 array1 low2 high2 array2 =
-  lenList t == lenList u && eqSameLength t u
-  where
-    t = TermList (I# low1) (I# high1) (ByteArray array1)
-    u = TermList (I# low2) (I# high2) (ByteArray array2)
-    eqSameLength Empty !_ = True
-    eqSameLength (ConsSym s1 t) (UnsafeConsSym s2 u) =
-      root s1 == root s2 && eqSameLength t u
+  t == u = compare t u == EQ
 
 instance Ord (TermList f) where
   {-# INLINE compare #-}
   compare t u =
-    case compare (lenList t) (lenList u) of
-      EQ -> compareContents t u
-      x  -> x
-
-compareContents :: TermList f -> TermList f -> Ordering
-compareContents Empty !_ = EQ
-compareContents (ConsSym s1 t) (UnsafeConsSym s2 u) =
-  case compare (root s1) (root s2) of
-    EQ -> compareContents t u
-    x  -> x
+    compare (lenList t) (lenList u) `mappend`
+    compareByteArrays (array t) (low t * k)
+      (array u) (low u * k) ((high t - low t) * k)
+    where
+      k = sizeOf (fromSymbol undefined)
 
 --------------------------------------------------------------------------------
 -- Building terms.
@@ -310,10 +266,11 @@
           case m s mbytearray# n# 0# of
             (# s, n# #) -> (# s, I# n# #)
       if n' <= n then do
+        resizeMutableByteArray (MutableByteArray mbytearray#) (n' * sizeOf (fromSymbol undefined))
         !bytearray <- unsafeFreezeByteArray (MutableByteArray mbytearray#)
         return (TermList 0 n' bytearray)
        else loop (n'*2)
-  loop 32
+  loop 128
 
 -- Get at the term array.
 {-# INLINE getByteArray #-}
@@ -408,22 +365,10 @@
 {-# INLINE isSubtermOfList #-}
 isSubtermOfList :: Term f -> TermList f -> Bool
 isSubtermOfList t u =
-  isSubArrayOf (singleton t) u
-
--- N.B. this one should not be exported from Twee.Term
--- because subarray is not the same as subterm if t is not
--- a singleton
-isSubArrayOf :: TermList f -> TermList f -> Bool
-isSubArrayOf t u =
-  lenList t <= lenList u && (here t u || next t u)
+  or [ singleton t == u{low = low u + i, high = low u + i + n}
+     | i <- [0..lenList u - n]]
   where
-    here Empty _ = True
-    here (ConsSym s1 t) (UnsafeConsSym s2 u) =
-      root s1 == root s2 && here t u
-
-    -- This is safe because lenList t <= lenList u
-    -- so if u = Empty, then t = Empty and here t u = True.
-    next t (UnsafeConsSym _ u) = isSubArrayOf t u
+    n = lenList (singleton t)
 
 -- | Check if a variable occurs in a termlist.
 {-# INLINE occursList #-}
@@ -432,4 +377,4 @@
 
 symbolOccursList :: Int64 -> TermList f -> Bool
 symbolOccursList !_ Empty = False
-symbolOccursList n (ConsSym t ts) = root t == n || symbolOccursList n ts
+symbolOccursList n ConsSym{hd = t, rest = ts} = root t == n || symbolOccursList n ts
diff --git a/Twee/Utils.hs b/Twee/Utils.hs
--- a/Twee/Utils.hs
+++ b/Twee/Utils.hs
@@ -11,6 +11,7 @@
 import GHC.Prim
 import GHC.Types
 import Data.Bits
+import System.Random
 --import Test.QuickCheck hiding ((.&.))
 
 repeatM :: Monad m => m a -> m [a]
@@ -70,12 +71,15 @@
   | otherwise = (x:xs) `isSubsequenceOf` ys
 #endif
 
-{-# INLINE fixpoint #-}
 fixpoint :: Eq a => (a -> a) -> a -> a
-fixpoint f x = fxp x
+fixpoint = fixpointOn id
+
+{-# INLINE fixpoint #-}
+fixpointOn :: Eq b => (a -> b) -> (a -> a) -> a -> a
+fixpointOn key f x = fxp x
   where
     fxp x
-      | x == y = x
+      | key x == key y = x
       | otherwise = fxp y
       where
         y = f x
@@ -121,3 +125,30 @@
   where
     splits = splitInterval k (lo, hi)
 -}
+
+reservoir :: Int -> [(Integer, Int)]
+reservoir k =
+  zip (map fromIntegral prefix) prefix ++
+  zip (map (+fromIntegral k) (scanl1 (+) is)) ks
+  where
+    xs, ys :: [Double]
+    xs = randomRs (0, 1) (mkStdGen 314159265)
+    ys = randomRs (0, 1) (mkStdGen 358979323)
+    ks = randomRs (0, k-1) (mkStdGen 846264338)
+
+    ws = scanl1 (*) [ x ** (1 / fromIntegral k) | x <- xs ]
+    is = zipWith gen ws ys
+    gen w y = floor (log y / log (1-w)) + 1
+    prefix = [0..k-1]
+
+-- A combined inits/tails.
+splits :: [a] -> [([a], [a])]
+splits [] = [([], [])]
+splits (x:xs) =
+  [([], x:xs)] ++
+  [(x:ys, zs) | (ys, zs) <- splits xs]
+
+-- Fold over the natural numbers.
+foldn :: (a -> a) -> a -> Int -> a
+foldn _ e 0 = e
+foldn op e n | n > 0 = op (foldn op e (n-1))
diff --git a/twee-lib.cabal b/twee-lib.cabal
--- a/twee-lib.cabal
+++ b/twee-lib.cabal
@@ -1,5 +1,5 @@
 name:                twee-lib
-version:             2.2
+version:             2.3
 synopsis:            An equational theorem prover
 homepage:            http://github.com/nick8325/twee
 license:             BSD3
@@ -69,8 +69,10 @@
     dlist,
     pretty >= 1.1.2.0,
     ghc-prim,
-    primitive >= 0.6.2.0,
-    vector
+    primitive >= 0.7.1.0,
+    vector,
+    uglymemo,
+    random
   hs-source-dirs:      .
   ghc-options:         -W -fno-warn-incomplete-patterns
   default-language:    Haskell2010