diff --git a/doc/manual/manual.md b/doc/manual/manual.md
index 15e184d248..709fe3d9ec 100644
--- a/doc/manual/manual.md
+++ b/doc/manual/manual.md
@@ -4055,39 +4055,45 @@ N.B: The following commands must first be enabled using `enable_experimental`.
 
 SAW contains a bisimulation prover to prove that two terms simulate each other.
 This prover allows users to prove that two terms executing in lockstep satisfy
-some relation over the state of each circuit and their outputs.  This type of
+some relations over the state of each circuit and their outputs.  This type of
 proof is useful in demonstrating the eventual equivalence of two circuits, or
 of a circuit and a functional specification.  SAW enables these proofs with the
 experimental `prove_bisim` command:
 
 ~~~~
-prove_bisim : ProofScript () -> Term -> Term -> Term -> TopLevel ProofResult
+prove_bisim : ProofScript () -> [BisimTheorem] -> Term -> Term -> Term -> Term -> TopLevel BisimTheorem
 ~~~~
 
-When invoking `prove_bisim strat relation lhs rhs`, the arguments represent the
-following:
+When invoking `prove_bisim strat theorems srel orel lhs rhs`, the arguments
+represent the following:
 
 1. `strat`: A proof strategy to use during verification.
-2. `relation`: A relation between `lhs` and `rhs`.  This relation must have the
-   type `(lhsState, output) -> (rhsState, output) -> Bit`. The relation's first
-   argument is a pair consisting of `lhs`'s state and output following
+2. `theorems`: A list of already proven bisimulation theorems.
+3. `srel`: A state relation between `lhs` and `rhs`.  This relation must have
+   the type `lhsState -> rhsState -> Bit`. The relation's first argument is
+   `lhs`'s state prior to execution. The relation's second argument is `rhs`'s
+   state prior to execution. `srel` then returns a `Bit` indicating whether
+   the two arguments satisfy the bisimulation's state relation.
+4. `orel`: An output relation between `lhs` and `rhs`.  This relation must have
+   the type `(lhsState, output) -> (rhsState, output) -> Bit`. The relation's
+   first argument is a pair consisting of `lhs`'s state and output following
    execution. The relation's second argument is a pair consisting of `rhs`'s
-   state and output following execution. `relation` then then returns a `Bit`
-   indicating whether the two arguments satisfy the bisimulation relation. That
-   is, whether the terms simulate each other.
-3. `lhs`: A term that simulates `rhs`. `lhs` must have the type
+   state and output following execution. `orel` then returns a `Bit` indicating
+   whether the two arguments satisfy the bisimulation's output relation.
+5. `lhs`: A term that simulates `rhs`. `lhs` must have the type
    `(lhsState, input) -> (lhsState, output)`.  The first argument to `lhs` is a
    tuple containing the internal state of `lhs`, as well as the input to `lhs`.
    `lhs` returns a tuple containing its internal state after execution, as well
    as its output.
-4. `rhs`: A term that simulates `lhs`. `rhs` must have the type
+6. `rhs`: A term that simulates `lhs`. `rhs` must have the type
    `(rhsState, input) -> (rhsState, output)`.  The first argument to `rhs` is a
    tuple containing the internal state of `rhs`, as well as the input to `rhs`.
    `rhs` returns a tuple containing its internal state after execution, as well
    as its output.
 
-`prove_bisim` returns a `ProofResult` indicating whether `lhs` and `rhs` are
-bisimilar, given the relation `relation`.
+On success, `prove_bisim` returns a `BisimTheorem` that can be used in future
+bisimulation proofs to enable compositional bisimulation proofs.  On failure,
+`prove_bisim` will abort.
 
 ## Bisimulation Example
 
@@ -4155,22 +4161,38 @@ is the result of the function's computation.
 It computes the logical AND directly on the function's inputs using Cryptol's
 `(&&)` operator.
 
+Next, we define a state relation over `andImp` and `andSpec`:
+
+~~~~
+andStateRel : andState -> () -> Bit
+andStateRel _ () = True
+~~~~
+
+`andStateRel` takes two arguments:
+
+1. An `andState` for `andImp`.
+2. An empty state (`()`) for `andSpec`.
+
+`andStateRel` returns a `Bit` indicating whether the relation is satisified.  In
+this case, `andStateRel` always returns `True` because `andSpec` is stateless
+and therefore the state relation permits `andImp` to accept any state.
+
 Lastly, we define a relation over `andImp` and `andSpec`:
 
 ~~~~
-andRel : (andState, (Bit, Bit)) -> ((), (Bit, Bit)) -> Bit
-andRel (s, (impReady, impO)) ((), (_, specO)) =
+andOutputRel : (andState, (Bit, Bit)) -> ((), (Bit, Bit)) -> Bit
+andOutputRel (s, (impReady, impO)) ((), (_, specO)) =
   if impReady then impO == specO else True
 ~~~~
 
-`andRel` takes two arguments:
+`andOutputRel` takes two arguments:
 
 1. A return value from `andImp`. Specifically, a pair consisting of an
    `andState` and a pair containing a ready bit and result of the logical AND.
 2. A return value from `andSpec`. Specifically, a pair consisting of an empty
    state `()` and a pair containing a ready bit and result of the logical AND.
 
-`andRel` returns a `Bit` indicating whether the relation is satisfied.  It
+`andOutputRel` returns a `Bit` indicating whether the relation is satisfied.  It
 considers the relation satisfied in two ways:
 
 1. If `andImp`'s ready bit is set, the relation is satisfied if the output
@@ -4188,12 +4210,14 @@ terms are bisimilar--by using `prove_bisim`:
 import "And.cry";
 enable_experimental;
 
-res <- prove_bisim z3 {{ andRel }} {{ andImp }} {{ andSpec }};
-print res;
+and_bisim <- prove_bisim z3 [] {{ andStateRel }} {{ andOutputRel }} {{ andImp }} {{ andSpec }};
 ~~~~
 
-Upon running this script, SAW prints `Valid` indicating that `andImp` and
-`andSpec` are bisimilar, given the relation `andRel`.
+Upon running this script, SAW prints:
+
+~~~~
+Successfully proved bisimulation between andImp and andSpec
+~~~~
 
 ### Building a NAND gate
 
@@ -4213,7 +4237,7 @@ result:
 nandImp : (nandState, (Bit, Bit)) -> (nandState, (Bit, Bit))
 nandImp x = (s, (andReady, ~andRes))
   where
-    (s ,(andReady, andRes)) = andImp x
+    (s, (andReady, andRes)) = andImp x
 ~~~~
 
 Note that `nandImp` is careful to preserve the ready status of `andImp`.
@@ -4230,41 +4254,110 @@ nandSpec (_, (x, y)) = ((), (True, ~ (andSpec ((), (x, y))).1.1))
 As with `andSpec`, `nandSpec` is pure and computes its result in a single
 cycle.
 
-Lastly, we define a relation indicating that `nandImp` and `nandSpec` produce
-equivalent results once `nandImp`'s ready bit is `1`:
+Next, we define a state relation over `nandImp` and `nandSpec`:
 
 ~~~~
-nandRel : (nandState, (Bit, Bit)) -> ((), (Bit, Bit)) -> Bit
-nandRel (s, (impReady, impO)) ((), (_, specO)) =
+nandStateRel : andState -> () -> Bit
+nandStateRel _ () = True
+~~~~
+
+As with `andStateRel`, this state relation is always `True` because `nandSpec`
+is stateless.
+
+Lastly, we define an output relation indicating that `nandImp` and `nandSpec`
+produce equivalent results once `nandImp`'s ready bit is `1`:
+
+~~~~
+nandOutputRel : (nandState, (Bit, Bit)) -> ((), (Bit, Bit)) -> Bit
+nandOutputRel (s, (impReady, impO)) ((), (_, specO)) =
   if impReady then impO == specO else True
 ~~~~
 
 To prove that `nandImp` and `nandSpec` are bisimilar, we again use
-`prove_bisim`:
+`prove_bisim`. This time however, we can reuse the bisimulation proof for the
+AND gate by including it in the `theorems` paramter for `prove_bisim`:
 
 ~~~~
-res2 <- prove_bisim z3 {{ nandRel }} {{ nandImp }} {{ nandSpec }};
-print res2;
+prove_bisim z3 [and_bisim] {{ nandStateRel }} {{ nandOutputRel }} {{ nandImp }} {{ nandSpec }};
 ~~~~
 
-Upon running this portion of the script, SAW prints `Valid` confirming that the
-two terms are bisimilar, given the relation `nandRel`.
+Upon running this script, SAW prints:
 
-## Understanding the proof goal
+~~~~
+Successfully proved bisimulation between nandImp and nandSpec
+~~~~
+
+## Understanding the proof goals
 
 While not necessary for simple proofs, more advanced proofs may require
-inspecting the proof goal.  `prove_bisim` generates and attempts to solve the
-proof goal:
+inspecting proof goals.  `prove_bisim` generates and attempts to solve the
+following proof goals:
 
 ~~~~
-forall s1 s2 in out1 out2.
-    relation (s1, out1) (s2, out2) -> relation (lhs (s1, in)) (rhs (s2, in))
+OUTPUT RELATION THEOREM:
+  forall s1 s2 in.
+    srel s1 s2 -> orel (lhs (s1, in)) (rhs (s2, in))
+
+STATE RELATION THEOREM:
+  forall s1 s2 out1 out2.
+    orel (s1, out1) (s2, out2) -> srel s1 s2
 ~~~~
 
-where the variables in the `forall` are:
+where the variables in the `forall`s are:
 
 * `s1`: Initial state for `lhs`
 * `s2`: Initial state for `rhs`
 * `in`: Input value to `lhs` and `rhs`
 * `out1`: Initial output value for `lhs`
 * `out2`: Initial output value for `rhs`
+
+The `STATE RELATION THEOREM` verifies that the output relation properly captures
+the guarantees of the state relation.  The `OUTPUT RELATION THEOREM` verifies
+that if `lhs` and `rhs` are executed with related states, then the result of
+that execution is also related.  These two theorems together guarantee that the
+terms simulate each other.
+
+When using composition, `prove_bisim` also generates and attempts to solve
+the proof goal below for any successfully applied `BisimTheorem` in the
+`theorems` list:
+
+~~~~
+COMPOSITION SIDE CONDITION:
+  forall g_lhs_s g_rhs_s.
+    g_srel g_lhs_s g_rhs_s -> f_srel f_lhs_s f_rhs_s
+    where
+      f_lhs_s = extract_inner_state g_lhs g_lhs_s f_lhs
+      f_rhs_s = extract_inner_state g_rhs g_rhs_s f_rhs
+~~~~
+
+where `g_lhs` is an outer term containing a call to an inner term `f_lhs`
+represented by a `BisimTheorem` and `g_rhs` is an outer term containing a call
+to an inner term `f_rhs` represented by the same `BisimTheorem`. The variables
+in `COMPOSITION SIDE CONDITION` are:
+
+* `extract_inner_state x x_s y`: A helper function that takes an outer term `x`, an
+  outer state `x_s`, and an inner term `y`, and returns the inner state of `x_s`
+  that `x` passes to `y`.
+* `g_lhs_s`: The state for `g_lhs`
+* `g_rhs_s`: The state for `g_rhs`
+* `g_srel`: The state relation for `g_lhs` and `g_rhs`
+* `f_srel`: The state relation for `f_lhs` and `f_rhs`
+* `f_lhs_s`: The state for `f_lhs`, as represented in `g_lhs_s` (extracted using
+  `extract_inner_state`).
+* `f_rhs_s`: The state for `f_rhs`, as represented in `g_rhs_s` (extracted using
+  `extract_inner_state`).
+
+The `COMPOSITION SIDE CONDITION` exists to verify that the terms in the
+bisimulation relation properly set up valid states for subterms they contain.
+
+## Limitiations
+
+For now, the `prove_bisim` command has a couple limitations:
+
+* `lhs` and `rhs` must be named functions.  This is because `prove_bisim` uses
+  these names to perform substitution when making use of compositionality.
+* Each subterm present in the list of bisimulation theorems already
+  proven may be invoked at most once in `lhs` or `rhs`.  That is, if some
+  function `g_lhs` calls `f_lhs`, and `prove_bisim` is invoked with a
+  `BisimTheorem` proving that `f_lhs` is bisimilar to `f_rhs`, then `g_lhs` may
+  call `f_lhs` at most once.
diff --git a/doc/manual/manual.pdf b/doc/manual/manual.pdf
index 631a6ee872..b5a529af27 100644
Binary files a/doc/manual/manual.pdf and b/doc/manual/manual.pdf differ
diff --git a/intTests/test_bisimulation/Comp.cry b/intTests/test_bisimulation/Comp.cry
index a66ca3dbd9..a866c63a4a 100644
--- a/intTests/test_bisimulation/Comp.cry
+++ b/intTests/test_bisimulation/Comp.cry
@@ -24,25 +24,35 @@ andImp (s, (x, y)) =
   then (s, (True, s.origX && s.origY))
   else ({ loaded = True, origX = x, origY = y }, (False, 0))
 
-// Relation between `andImp` and `andSpec`.  Checks that outputs are
+// Output relation between `andImp` and `andSpec`.  Checks that outputs are
 // equivalent, given that `andImp` is ready to be read.
-andRel : (andState, (Bit, Bit)) -> ((), (Bit, Bit)) -> Bit
-andRel (s, (impReady, impO)) ((), (_, specO)) =
+andOutputRel : (andState, (Bit, Bit)) -> ((), (Bit, Bit)) -> Bit
+andOutputRel (s, (impReady, impO)) ((), (_, specO)) =
   if impReady then impO == specO else True
 
-// Stateful NAND gate holds `andState`. It has no additional state.
-type nandState = andState
+// State relation between `andImp` and `andSpec`.  Trivial in this case because
+// `andSpec` is stateless.
+andStateRel : andState -> () -> Bit
+andStateRel _ () = True
+
+// Stateful NAND gate holds `andState` as well as an unused bit.
+type nandState = { as : andState, unused : Bit }
 
 // A stateful NAND gate.  First output bit is the ready bit, second bit is the
 // result.  This gate takes 2 cycles to compute.  It restarts when the input
 // changes.
 nandImp : (nandState, (Bit, Bit)) -> (nandState, (Bit, Bit))
-nandImp x = (s, (andReady, ~andRes))
+nandImp (s, in) = ({ as = as' , unused = False }, (andReady, ~andRes))
   where
-    (s ,(andReady, andRes)) = andImp x
+    (as', (andReady, andRes)) = andImp (s.as, in)
 
 // Relation between `nandImp` and `nandSpec`.  Checks that outputs are
 // equivalent, given that `nandImp` is ready to be read.
-nandRel : (nandState, (Bit, Bit)) -> ((), (Bit, Bit)) -> Bit
-nandRel (s, (impReady, impO)) ((), (_, specO)) =
+nandOutputRel : (nandState, (Bit, Bit)) -> ((), (Bit, Bit)) -> Bit
+nandOutputRel (s, (impReady, impO)) ((), (_, specO)) =
   if impReady then impO == specO else True
+
+// State relation between `nandImp` and `nandSpec`.  Trivial in this case
+// because `nandSpec` is stateless.
+nandStateRel : nandState -> () -> Bit
+nandStateRel _ () = True
diff --git a/intTests/test_bisimulation/Test.cry b/intTests/test_bisimulation/Test.cry
index 7bab304282..93874bbf80 100644
--- a/intTests/test_bisimulation/Test.cry
+++ b/intTests/test_bisimulation/Test.cry
@@ -20,6 +20,10 @@ f2_bad_input_type (ctr, inp) = (ctr+1, (drop inp)+(drop ctr))
 f2_bad_output_type : (S2, [8]) -> (S2, [9])
 f2_bad_output_type (ctr, inp) = (ctr+1, zext (inp+(drop ctr)))
 
-// Bisimulation relation for f1 and f2
-rel : (S1, [8]) -> (S2, [8]) -> Bit
-rel (s1, o1) (s2, o2) = s1 == drop s2 /\ o1 == o2
+// Output bisimulation relation for f1 and f2
+orel : (S1, [8]) -> (S2, [8]) -> Bit
+orel (s1, o1) (s2, o2) = s1 == drop s2 /\ o1 == o2
+
+// State bisimulation relation for f1 and f2
+srel : S1 -> S2 -> Bit
+srel s1 s2 = s1 == drop s2
diff --git a/intTests/test_bisimulation/comp.saw b/intTests/test_bisimulation/comp.saw
index fd2ee53cae..b9a96720e9 100644
--- a/intTests/test_bisimulation/comp.saw
+++ b/intTests/test_bisimulation/comp.saw
@@ -1,13 +1,29 @@
-/* Test the prove_bisim command with a compositional circuit. At the moment,
- * prove_bisim does not have any notion of compositionality, but this example
- * is simple enough to work as-is */
+/* Test the prove_bisim command with a compositional circuit */
 
 import "Comp.cry";
 
 enable_experimental;
 
-res <- prove_bisim z3 {{ andRel }} {{ andImp }} {{ andSpec }};
-print res;
+// Prove 'andImp' and 'andSpec' simulate each other
+and_res <- prove_bisim z3
+                       []
+                       {{ andStateRel }}
+                       {{ andOutputRel }}
+                       {{ andImp }}
+                       {{ andSpec }};
 
-res2 <- prove_bisim z3 {{ nandRel }} {{ nandImp }} {{ nandSpec }};
-print res2;
+// Prove 'nandImp' and 'nandSpec' simulate each other, using 'and_res' in the
+// process.
+nand_res <- prove_bisim (do {
+    unfolding ["nandImp", "nandSpec"];
+    w4_unint_z3 ["andImp", "andSpec" ];
+  })
+  [and_res] {{ \x -> nandStateRel x }} {{ \x -> nandOutputRel x }} {{ nandImp }} {{ nandSpec }};
+
+// Test using theorem that doesn't apply
+prove_bisim z3
+            [nand_res]
+            {{ andStateRel }}
+            {{ andOutputRel }}
+            {{ andImp }}
+            {{ andSpec }};
diff --git a/intTests/test_bisimulation/test.saw b/intTests/test_bisimulation/test.saw
index 86b1393d8f..d0b75feec0 100644
--- a/intTests/test_bisimulation/test.saw
+++ b/intTests/test_bisimulation/test.saw
@@ -4,17 +4,22 @@ import "Test.cry";
 
 enable_experimental;
 
-res <- prove_bisim z3 {{ rel }} {{ f1 }} {{ f2 }};
-print res;
+res <- prove_bisim z3 [] {{ srel }} {{ orel }} {{ f1 }} {{ f2 }};
+
+// Test non-named constant in bisimulation term
+fails (prove_bisim z3 [] {{ srel }} {{ orel }} {{ f1 }} {{ \x -> f2 x }});
 
 // Test incompatable input types
-fails (prove_bisim z3 {{ rel }} {{ f1 }} {{ f2_bad_input_type }});
+fails (prove_bisim z3 [] {{ srel }} {{ orel }} {{ f1 }} {{ f2_bad_input_type }});
 
 // Test incompatable output types
-fails (prove_bisim z3 {{ rel }} {{ f1 }} {{ f2_bad_output_type }});
+fails (prove_bisim z3 [] {{ srel }} {{ orel }} {{ f1 }} {{ f2_bad_output_type }});
+
+// Test bad output relation type
+fails (prove_bisim z3 [] {{ srel }} {{ True }} {{ f1 }} {{ f2 }});
 
-// Test bad relation type
-fails (prove_bisim z3 {{ True }} {{ f1 }} {{ f2 }});
+// Test bad state relation type
+fails (prove_bisim z3 [] {{ True }} {{ orel }} {{ f1 }} {{ f2 }});
 
 // Test wrong state type
-fails (prove_bisim z3 {{ rel }} {{ f2 }} {{ f2 }});
+fails (prove_bisim z3 [] {{ srel }} {{ orel }} {{ f2 }} {{ f2 }});
diff --git a/saw-core/src/Verifier/SAW/SharedTerm.hs b/saw-core/src/Verifier/SAW/SharedTerm.hs
index c573994c76..7a3fb99580 100644
--- a/saw-core/src/Verifier/SAW/SharedTerm.hs
+++ b/saw-core/src/Verifier/SAW/SharedTerm.hs
@@ -253,6 +253,7 @@ module Verifier.SAW.SharedTerm
 --  , scFalse
    , scOpenTerm
    , scCloseTerm
+   , scLambdaBody
     -- ** Variable substitution
   , instantiateVar
   , instantiateVarList
@@ -2589,7 +2590,6 @@ scGeneralizeExts sc exts x = loop (zip (inits exts) exts)
     -- base case, convert all the exts in the body of x into deBruijn variables
     loop [] = scExtsToLocals sc exts x
 
-
 scUnfoldConstants :: SharedContext -> [VarIndex] -> Term -> IO Term
 scUnfoldConstants sc names t0 = scUnfoldConstantSet sc True (Set.fromList names) t0
 
@@ -2708,3 +2708,14 @@ scCloseTerm close sc ec body = do
     lv <- scLocalVar sc 0
     body' <- scInstantiateExt sc (Map.insert (ecVarIndex ec) lv Map.empty) =<< incVars sc 0 1 body
     close sc (toShortName (ecName ec)) (ecType ec) body'
+
+-- | Compute the body of 0 or more nested lambda-abstractions by applying the
+-- lambdas to fresh 'ExtCns's. Note that we do this lambda-by-lambda, rather
+-- than generating all of the 'ExtCns's at the same time, because later
+-- variables could have types that depend on earlier variables, and so the
+-- substitution of earlier 'ExtCns's has to happen before we generate later
+-- ones.
+scLambdaBody :: SharedContext -> Term -> IO Term
+scLambdaBody sc (asLambda -> Just (nm, tp, body)) =
+  scOpenTerm sc nm tp 0 body >>= scLambdaBody sc . snd
+scLambdaBody _sc t = return t
diff --git a/saw-script.cabal b/saw-script.cabal
index 8d75960379..4d596c0018 100644
--- a/saw-script.cabal
+++ b/saw-script.cabal
@@ -114,6 +114,7 @@ library
     SAWScript.AutoMatch.JVM
     SAWScript.AutoMatch.Util
     SAWScript.Bisimulation
+    SAWScript.Bisimulation.BisimTheorem
     SAWScript.Builtins
     SAWScript.CongruenceClosure
     SAWScript.Exceptions
diff --git a/src/SAWScript/Bisimulation.hs b/src/SAWScript/Bisimulation.hs
index f2204ad6b8..726ea59409 100644
--- a/src/SAWScript/Bisimulation.hs
+++ b/src/SAWScript/Bisimulation.hs
@@ -13,144 +13,575 @@ functionality in the future.
 At its core, we want to prove that two circuits executing in lockstep satisfy
 some relation over the state of each circuit and their outputs. To achieve this,
 the 'proveBisimulation' command takes:
-* A relation @rel : (lhsState, output) -> (rhsState, output) -> Bit@
+* A list of already proved bisimulation theorems
+* A state relation @srel : lhsState -> rhsState -> Bit@
+* An output relation @orel : (lhsState, output) -> (rhsState, output) -> Bit@
 * A term @lhs : (lhsState, input) -> (lhsState, output)@
 * A term @rhs : (rhsState, input) -> (rhsState, output)@
-and considers @lhs@ and @rhs@ bisimilar when:
-  forall s1 s2 in out1 out2.
-    rel (s1, out1) (s2, out2) -> rel (lhs (s1, in)) (rhs (s2, in))
-
-One natural question is why the relation has the type
-@(lhsState, output) -> (rhsState, output) -> Bit@ instead of something simpler
-like @lhsState -> rhsState -> Bit@. We require the user to specify when outputs
-satisfy the relation because it is not always the case that circuit outputs
-agree exactly after each simulation step. Specifically, one circuit may complete
-some calculation in @N@ steps and another in @M@ steps where @N != M@. In this
-case, the user may not want to check the circuit outputs for equality until the
-greater of @N@ and @M@ steps have passed. The simpler relation type would not
-enable specification of this property.
-
-The main downside of this approach is that the resulting forall in the formula
-sent to the SMT solver quantifies over the initial output of the circuits prior
-to simulating a step. This can be confusing when reading the SAW source code,
-and could be resolved by requiring the user to provide two different relations
-(one over states, and one over states and outputs), but this puts more burden on
-the end user who would have to write two relations rather than just one. As
-such, we've chosen the approach that's easier on the end user, rather than the
-one that's easier on the SAW implementer.
+and considers @lhs@ and @rhs@ bisimilar when the following two theorems hold:
+* OUTPUT RELATION THEOREM:
+    forall s1 s2 in.
+      srel s1 s2 -> orel (lhs (s1, in)) (rhs (s2, in))
+* STATE RELATION THEOREM:
+    forall s1 s2 out1 out2.
+      orel (s1, out1) (s2, out2) -> srel s1 s2
+
+The OUTPUT RELATION THEOREM ensures that if the state relation holds prior to
+executing @lhs@ and @rhs@, then the output relation holds after executing @lhs@
+and @rhs@.  That is, the two terms step together.
+
+The STATE RELATION THEOREM ensures that if the output relation holds over some
+state and output value, then the state relation also holds over that same state.
+This ensures that the output relation captures the state relation.
+
+To enable compositionality, 'proveBisimulation' takes a list of proven
+bisimulation theorems as 'BisimTheorem's.  For each 'BisimTheorem', this
+implementation searches for invocations of the functions in the 'BisimTheorem'
+and replaces them with uninterpreted values that satisfy the 'BisimTheorem''s
+output relation.  Put another way, if @g_lhs@ contains some function @f_lhs@ and
+@g_rhs@ contains some function @f_rhs@ where @f_lhs@ and @f_rhs@ satisfy some
+output relation @f_orel@, then in the proof of some OUTPUT RELATION THEOREM for
+@g_lhs@ and @g_rhs@, the prover can replace @f_lhs@ with an uninterpreted value
+@v_lhs@ and @f_rhs@ with an uninterpreted value @v_rhs@.  Lastly, the prover
+assumes @v_lhs@ and @v_rhs@ satisfy @f_orel@.
+
+In order to make this assumption about @v_lhs@ and @v_rhs@, the state relation
+of the outer terms must ensure that the inner terms' state relation holds.  As
+such, the prover generates and proves a side condition for each applied
+'BisimTheorem'.  Let @g_lhs_s@ and @g_rhs_s@ be the states for @g_lhs@ and
+@g_rhs@ respectively.  Additionally, let there be a function
+@extract_inner_state x x_s y@ that takes an outer term @x@, an outer state
+@x_s@, and an inner term @y@, and returns the inner state of @x_s@ that @x@
+passes to @y@.  Lastly, let @g_srel@ be the state relation for the @g@ terms,
+and @f_srel@ be the state relation for the @f@ terms.  The prover then checks:
+  COMPOSITION SIDE CONDITION:
+    forall g_lhs_s g_rhs_s.
+      g_srel g_lhs_s g_rhs_s -> f_srel f_lhs_s f_rhs_s
+      where
+        f_lhs_s = extract_inner_state g_lhs g_lhs_s f_lhs
+        f_rhs_s = extract_inner_state g_rhs g_rhs_s f_rhs
+
+The reason for all of this complexity around composition is ultimately to reduce
+the burden on the SMT solver by uninterpreting functions so that the SMT solver
+may handle more complex proofs.
 -}
+
+{-# LANGUAGE LambdaCase #-}
 {-# LANGUAGE OverloadedStrings #-}
 {-# LANGUAGE TupleSections #-}
 
 module SAWScript.Bisimulation
-  ( proveBisimulation )
+  ( BisimTheorem, proveBisimulation )
   where
 
-import Control.Monad (unless)
+import Control.Monad.State.Strict as State
+import Data.Foldable (foldl')
+import qualified Data.IntSet as IntSet
+import qualified Data.Map.Strict as Map
+import qualified Data.Text as Text
 
 import qualified Cryptol.TypeCheck.Type as C
 import qualified Cryptol.Utils.PP as C
 
-import SAWScript.Builtins (provePrim)
-import SAWScript.Crucible.Common.MethodSpec (ppTypedTermType)
+import SAWScript.Bisimulation.BisimTheorem
+import SAWScript.Builtins (unfold_term)
+import SAWScript.Crucible.Common.MethodSpec (ppTypedTermType, ppTypedTerm)
+import SAWScript.Options (Verbosity(..))
+import SAWScript.Panic (panic)
 import SAWScript.Proof
+import SAWScript.Prover.Util (checkBooleanSchema)
 import SAWScript.Value
 
 import qualified Verifier.SAW.Cryptol as C
 import Verifier.SAW.SharedTerm
+import Verifier.SAW.Term.Functor
 import Verifier.SAW.TypedTerm
+import Verifier.SAW.Recognizer
+
+-- State used to facilitate the replacement of a 'Constant' application in a
+-- 'Term' with an 'ExtCns'.  Used in 'replaceConstantTerm' and
+-- 'replaceConstantTermF'
+data ReplaceState = ReplaceState {
+    rsMemo :: Map.Map TermIndex Term
+ -- ^ Memoization table to avoid re-visiting the same shared term
+  , rsExtCns :: Maybe Term
+ -- ^ ExtCns that replaces the 'Constant' application, if the constant could be
+ -- located.
+  , rsApp :: Maybe (TermF Term)
+ -- ^ Application that was replaced, if it could be located.
+  }
+
+-- Components needed for a bisimulation proof.
+data BisimComponents = BisimComponents {
+    bcTheorem :: BisimTheorem
+ -- ^ Bisimulation theorem capturing relations and outputs
+  , bcInputType :: C.Type
+ -- ^ Input type of the bisimilar functions
+  }
 
--- | Use bisimulation to prove that two terms simulate each other.
+-- An initial 'ReplaceState'
+emptyReplaceState :: ReplaceState
+emptyReplaceState = ReplaceState Map.empty Nothing Nothing
+
+-- | Translate a list of 'TypedTerm' representing theorems for input to the
+-- given validity-checking script and attempt to prove them.
+proveAll :: ProofScript () -> [TypedTerm] -> TopLevel ()
+proveAll script ts = do
+  sc <- getSharedContext
+  pos <- getPosition
+  forM_ (zip [0..] ts) $ \(n, t) -> do
+    io $ checkBooleanSchema (ttType t)
+    prop <- io $ predicateToProp sc Universal (ttTerm t)
+    let goal = ProofGoal
+              { goalNum  = n
+              , goalType = "prove"
+              , goalName = "prove_bisim"
+              , goalLoc  = show pos
+              , goalDesc = ""
+              , goalSequent = propToSequent prop
+              , goalTags = mempty
+              }
+    res <- runProofScript script prop goal Nothing "prove_bisim" True False
+    case res of
+      UnfinishedProof {} -> failProof res
+      ValidProof _ thm -> recordProof thm
+      InvalidProof {} ->  failProof res
+  where
+    failProof :: ProofResult -> TopLevel ()
+    failProof res = do
+      opts <- rwPPOpts <$> getTopLevelRW
+      fail $ "prove_bisim failed.\n" ++ showsProofResult opts res ""
+
+-- | Generate 'Term' for application of a relation
+scRelation :: TypedTerm -> Term -> Term -> TopLevel Term
+scRelation rel relLhs relRhs = do
+  sc <- getSharedContext
+  io $ scApplyAll sc (ttTerm rel) [relLhs, relRhs]
+
+-- | Build the COMPOSITION SIDE CONDITION for 'bc' and 'bt'.  See the
+-- documentation at the top of this file for information on the COMPOSITION SIDE
+-- CONDITION.
 --
--- Given the following:
--- * A relation @rel : (lhsState, output) -> (rhsState, output) -> Bit@
--- * A term @lhs : (lhsState, input) -> (lhsState, output)@
--- * A term @rhs : (rhsState, input) -> (rhsState, output)@
--- the prover considers @lhs@ and @rhs@ bisimilar when:
---   forall s1 s2 in out1 out2.
---     rel (s1, out1) (s2, out2) -> rel (lhs (s1, in)) (rhs (s2, in))
-proveBisimulation ::
-  ProofScript () ->
-  -- ^ Proof script to use over generated bisimulation term
-  TypedTerm ->
-  -- ^ Relation over states and outputs for terms to prove bisimilar. Must have
-  -- type @(lhsState, output) -> (rhsState, output) -> Bit@.
-  TypedTerm ->
-  -- ^ LHS of bisimulation. Must have type
-  -- @(lhsState, input) -> (lhsState, output)@
-  TypedTerm ->
-  -- ^ RHS of bisimulation. Must have type
-  -- @(rhsState, input) -> (rhsState, output)@
-  TopLevel ProofResult
-proveBisimulation script relation lhs rhs = do
+-- This function works by examining the specific calls to @f_lhs@ in @g_lhs@ and
+-- @f_rhs@ in @g_rhs@ to deduce how these subfunctions access state in their
+-- respective super functions.  It then extracts these accessors and
+-- instantiates them with the specific @g_lhs_s@ and @g_rhs_s@ values passed
+-- into @g_lhs@ and @g_rhs@ respectively.
+buildCompositionSideCondition
+  :: BisimComponents
+  -- ^ Components of the outer bisimulation under verification
+  -> BisimTheorem
+  -- ^ Bisimulation theorem concerning inner function
+  -> TopLevel TypedTerm
+buildCompositionSideCondition bc innerBt = do
   sc <- getSharedContext
+  let outerBt = bcTheorem bc
 
-  -- Typechecking
-  (lhsStateType, rhsStateType, outputType) <- typecheckRelation
-  lhsInputType <- typecheckSide lhs lhsStateType outputType
-  rhsInputType <- typecheckSide rhs rhsStateType outputType
-  unless (lhsInputType == rhsInputType) $
-    fail $ unlines [ "Error: Mismatched input types in bisimulation terms."
-                   , "  LHS input type: " ++ C.pretty lhsInputType
-                   , "  RHS input type: " ++ C.pretty rhsInputType ]
+  lhsOuterState <- io $ scLocalVar sc 0        -- g_lhs_s
+  rhsOuterState <- io $ scLocalVar sc 1        -- g_rhs_s
+
+  -- NOTE: Although not used in the final formula, we need to capture the input
+  -- to the outer functions because the extracted inner function applications
+  -- depend on it.  Therefore, it is necessary to match the expected form of the
+  -- inner ExtCns that this function instantiates.
+  input    <- io $ scLocalVar sc 2        -- in
+
+  -- Locate inner function calls on each side and replace their arguments with
+  -- 'ExtCns's
+  (lhsInnerEc, lhsInnerApp) <-
+    openConstantApp (bisimTheoremLhs innerBt) (bisimTheoremLhs outerBt)
+  (rhsInnerEc, rhsInnerApp) <-
+    openConstantApp (bisimTheoremRhs innerBt) (bisimTheoremRhs outerBt)
+
+  -- Extract state accessors from each inner function
+  lhsInnerState <- stateFromApp lhsInnerApp  -- f_lhs_s
+  rhsInnerState <- stateFromApp rhsInnerApp  -- f_rhs_s
+
+  -- Outer state relation
+  -- g_srel g_lhs_s g_rhs_s
+  outerRel <-
+    scRelation (bisimTheoremStateRelation outerBt) lhsOuterState rhsOuterState
+
+  -- Inner state relation
+  -- f_srel f_lhs_s f_rhs_s
+  innerRel <-
+    scRelation (bisimTheoremStateRelation innerBt) lhsInnerState rhsInnerState
+
+  -- Replace extcns in inner relation with outer inputs
+  lhsTuple <- io $ scTuple sc [lhsOuterState, input]  -- (f_lhs_s, in)
+  rhsTuple <- io $ scTuple sc [rhsOuterState, input]  -- (f_rhs_s, in)
+  innerRel' <- io $
+    scInstantiateExt sc
+                     (Map.fromList [ (ecVarIndex lhsInnerEc, lhsTuple)
+                                   , (ecVarIndex rhsInnerEc, rhsTuple)])
+                     innerRel
+
+  -- outer state relation implies inner state relation
+  -- g_srel g_lhs_s g_rhs_s -> f_srel f_lhs_s f_rhs_s
+  implication <- io $ scImplies sc outerRel innerRel'
+
+  -- Theorem to prove. Note that the 'input' is ultimately unused (see NOTE on
+  -- 'input' at the top of this function).
+  -- forall g_lhs_s g_rhs_s. g_srel g_lhs_s g_rhs_s -> f_srel f_lhs_s f_rhs_s
+  args <- io $ mapM
+    (\(name, t) -> (name,) <$> C.importType sc C.emptyEnv t)
+    [ ("input", bcInputType bc)
+    , ("rhsState", bisimTheoremRhsStateType outerBt)
+    , ("lhsState", bisimTheoremLhsStateType outerBt) ]
+  theorem <- io $ scLambdaList sc args implication
+  io $ mkTypedTerm sc theorem
+
+-- | Extract the state from the 'App' within a bisimulation side. Fails if 'app'
+-- is not an 'App'.
+stateFromApp :: TermF Term -> TopLevel Term
+stateFromApp app = do
+  sc <- getSharedContext
+  case app of
+    App _ arg -> io $ scFlatTermF sc $ PairLeft arg
+    _ -> do
+      term <- io $ scTermF sc app
+      fail $ "Error: " ++ showTerm term ++ " is not an App"
+
+-- | Given a term containing the application of a 'Constant', locate this
+-- application and replace its argument with an 'ExtCns'.  Returns the inserted
+-- 'ExtCns' and the updated 'App'.  Fails if 'constant' is not a 'Constant'.
+openConstantApp :: TypedTerm
+                -- ^ 'Constant' to search for.
+                -> TypedTerm
+                -- ^ 'TypeTerm' to locate 'Constant' in.  Must itself be a
+                -- 'Constant' (will be unfolded).
+                -> TopLevel (ExtCns Term, TermF Term)
+openConstantApp constant t = do
+  sc <- getSharedContext
+  -- Unfold constant
+  name <- Text.unpack <$> constantName (unwrapTermF (ttTerm t))
+  tUnfolded <- unfold_term [name] t
+
+  -- Break down lambda into its component parts.
+  (nm, tp, body) <- lambdaOrFail tUnfolded
+
+  -- Replace outer function's argument with an 'ExtCns'
+  -- NOTE: The bisimulation relation type ensures this is a single argument
+  -- lambda, so it's OK to apply scOpenTerm once and not recurse
+  (ec, tExtconsified) <- io $ scOpenTerm sc nm tp 0 body
+  extractedF <- extractApp constant tExtconsified
+  pure (ec, extractedF)
+
+  where
+    -- Break down lambda into its component parts.  Fails if 'tt' is not a
+    -- lambda.
+    lambdaOrFail :: TypedTerm -> TopLevel (LocalName, Term, Term)
+    lambdaOrFail tt =
+      case asLambda (ttTerm tt) of
+        Just lambda -> return lambda
+        Nothing ->
+          fail $ "Error: Expected a lambda term, got " ++ show (ppTypedTerm tt)
+
+-- Traverses 'term' and extracts the application of 'constant' within it.  Fails
+-- if 'constant' cannot be found within 'term'.
+extractApp :: TypedTerm
+           -- ^ 'Constant' to search for.
+           -> Term
+           -- ^ 'Term' to search in.
+           -> TopLevel (TermF Term)
+extractApp constant term =
+  case snd $ go (IntSet.empty, Nothing) term of
+    Just app -> return app
+    Nothing -> fail $ unlines [ "Error: Failed to locate constant in term."
+                              , "  Constant: " ++ show (ppTypedTerm constant)
+                              , "  Term: " ++ showTerm term ]
+  where
+    go :: (IntSet.IntSet, Maybe (TermF Term))
+       -> Term
+       -> (IntSet.IntSet, Maybe (TermF Term))
+    go (seen, acc) t =
+      case acc of
+        Just res -> (seen, Just res)
+        Nothing ->
+          case t of
+            Unshared tf -> termf (seen, acc) tf
+            STApp{ stAppIndex = i, stAppTermF = tf } ->
+              if IntSet.member i seen
+              then (seen, acc)
+              else termf (IntSet.insert i seen, acc) tf
+
+    termf :: (IntSet.IntSet, Maybe (TermF Term))
+          -> TermF Term
+          -> (IntSet.IntSet, Maybe (TermF Term))
+    termf (seen, acc) tf =
+      case tf of
+        App fn _ | unwrapTermF fn == unwrapTermF (ttTerm constant) ->
+          (seen, Just tf)
+        _ -> foldl' go (seen, acc) tf
+
+-- | Apply a prior proof result to simplify a bisimulation 'Term'.  Returns an
+-- updated 'Term', and a potential side condition that must be proven (or
+-- 'Nothing' if simplification had no effect).
+--
+-- This function works by searching for the functions from 'bt' in 'term'.  If
+-- it finds them, it replaces the calls in 'term' with 'ExtCns's and adds an
+-- assumption that the output relation holds over those 'ExtCns's.  This has the
+-- effect of saving the SMT solver from having to reason again about
+-- subrelations that have already been proven.  More formally, if @g_1@ and
+-- @g_2@ contain subfunctions @f_1@ and @f_2@ that satisfy some output relation
+-- @orel@, then this function replaces @f_1@ and @f_2@ with the uninterpreted
+-- values @v_1@ and @v_2@ respectively and adds the assumption @orel v_1 v_2@.
+--
+-- When simplification succeeds, this function also returns the side condition
+-- that the outer state relation from @g_1@ and @g_2@ implies that the inner
+-- state relation on @f_1@ and @f_2@ holds.  See the documentation on
+-- 'buildSuperStateRelation' for more information on the side condition.
+applyTheorem :: BisimTheorem
+             -- ^ Theorem to apply
+             -> BisimComponents
+             -- ^ Components of bisimulation under verification
+             -> Term
+             -- ^ Term to simplify
+             -> TopLevel (Term, Maybe TypedTerm)
+applyTheorem bt bc term0 = do
+  sc <- getSharedContext
+
+  -- Attempt to replace calls in 'term0' to the left side of 'bt' with an
+  -- 'ExtCns'.  In other words, replace @f_1@ with @v_1@.
+  let tpLhs = C.TCon (C.TC (C.TCTuple 2)) [ bisimTheoremLhsStateType bt
+                                          , bisimTheoremOutputType bt ]
+  (term1, lhsRs) <-
+    State.runStateT (replaceConstantTerm (bisimTheoremLhs bt) tpLhs term0)
+                    emptyReplaceState
+
+  -- Attempt to replace calls in 'term1' to the right side of 'bt' with an
+  -- 'ExtCns'.  In other words, replace @f_2@ with @v_2@.
+  let tpRhs = C.TCon (C.TC (C.TCTuple 2)) [ bisimTheoremRhsStateType bt
+                                          , bisimTheoremOutputType bt ]
+  (term2, rhsRs) <-
+    State.runStateT (replaceConstantTerm (bisimTheoremRhs bt) tpRhs term1)
+                    emptyReplaceState
+
+  case (rsExtCns lhsRs, rsExtCns rhsRs) of
+    (Just lhsEc, Just rhsEc) -> do
+      -- Simplification succeeded! Add assumption and generate side condition
+
+      -- Apply inner output relation to 'ExtCns's
+      -- orel v_1 v_2
+      app <- io $ scApplyAll sc
+                              (ttTerm (bisimTheoremOutputRelation bt))
+                              [ lhsEc , rhsEc ]
+
+      -- Use implication to assume relation holds in the main condition
+      -- orel v_1 v_2 -> term2
+      main <- io $ scImplies sc app term2
+
+      -- Build side condition
+      side <- buildCompositionSideCondition bc bt
+
+      pure (main, Just side)
+    _ ->
+      -- Application failed to match on one or both sides.  Return the Term
+      -- unchanged.
+      pure (term0, Nothing)
+
+-- | Apply multiple proof results to simplify a bisimulation 'Term'.  Returns an
+-- updated 'Term' and a list of side conditions that must be proven.  See the
+-- documentation for 'applyTheorem' for more info.
+applyAllTheorems :: [BisimTheorem]
+                 -- ^ Theorems to apply
+                 -> BisimComponents
+                 -- ^ Components of bisimulation under verification
+                 -> Term
+                 -- ^ Term to simplify
+                 -> TopLevel (Term, [TypedTerm])
+applyAllTheorems bthms bc term = foldM go (term, []) bthms
+  where
+    go :: (Term, [TypedTerm]) -> BisimTheorem -> TopLevel (Term, [TypedTerm])
+    go (curTerm, sides) bt = do
+      (nextTerm, mSideCondition) <- applyTheorem bt bc curTerm
+      case mSideCondition of
+        Nothing ->
+          -- Application failed
+          pure (nextTerm, sides)
+        Just sideCondition -> pure (nextTerm, sideCondition : sides)
+
+
+-- | Build the OUTPUT RELATION THEOREM for 'bc', as well as COMPOSITION SIDE
+-- CONDITION terms for any theorems used in the OUTPUT RELATION THEOREM.
+-- Returns a list of TypedTerms that must be proved. See the documentation at
+-- the top of this file for information on OUTPUT RELATION THEOREM and
+-- COMPOSITION SIDE CONDITION.
+buildOutputRelationTheorem :: [BisimTheorem]
+                           -- ^ Theorems available to use in proof
+                           -> BisimComponents
+                           -- ^ Components of bisimulation under verification
+                           -> TopLevel [TypedTerm]
+buildOutputRelationTheorem bthms bc = do
+  sc <- getSharedContext
+  let outerBt = bcTheorem bc
 
   -- Outer function inputs. See comments to the right of each line to see how
-  -- they line up with the @forall@ in the haddocs for this function.
-  input <- io $ scLocalVar sc 0           -- in
-  lhsState <- io $ scLocalVar sc 1        -- s1
-  rhsState <- io $ scLocalVar sc 2        -- s2
-  initLhsOutput <- io $ scLocalVar sc 3   -- out1
-  initRhsOutput <- io $ scLocalVar sc 4   -- out2
+  -- they line up with the documentation at the top of this file.
+  lhsState <- io $ scLocalVar sc 0        -- s1
+  rhsState <- io $ scLocalVar sc 1        -- s2
+  input <- io $ scLocalVar sc 2           -- in
 
   -- LHS/RHS inputs
   lhsTuple <- io $ scTuple sc [lhsState, input]  -- (s1, in)
   rhsTuple <- io $ scTuple sc [rhsState, input]  -- (s2, in)
 
   -- LHS/RHS outputs
-  lhsOutput <- io $ scApply sc (ttTerm lhs) lhsTuple  -- lhs (s1, in)
-  rhsOutput <- io $ scApply sc (ttTerm rhs) rhsTuple  -- rhs (s2, in)
-
-  -- Initial relation inputs
-  initRelArg1 <- io $ scTuple sc [lhsState, initLhsOutput]  -- (s1, out1)
-  initRelArg2 <- io $ scTuple sc [rhsState, initRhsOutput]  -- (s2, out2)
+  -- lhs (s1, in)
+  lhsOutput <- io $ scApply sc (ttTerm (bisimTheoremLhs outerBt)) lhsTuple
+  -- rhs (s2, in)
+  rhsOutput <- io $ scApply sc (ttTerm (bisimTheoremRhs outerBt)) rhsTuple
 
   -- Initial relation result
-  -- rel (s1, out1) (s2, out2)
-  initRelation <- scRelation sc initRelArg1 initRelArg2
+  -- srel s1 s2
+  initRelation <-
+    scRelation (bisimTheoremStateRelation outerBt) lhsState rhsState
 
   -- Relation over outputs
-  -- rel (lhs (s1, in)) (rhs (s2, in))
-  relationRes <- scRelation sc lhsOutput rhsOutput
+  -- orel (lhs (s1, in)) (rhs (s2, in))
+  relationRes <-
+    scRelation (bisimTheoremOutputRelation outerBt) lhsOutput rhsOutput
 
   -- initRelation implies relationRes
-  -- rel (s1, out1) (s2, out2) -> rel (lhs (s1, in)) (rhs (s2, in))
+  -- srel s1 s2 -> orel (lhs (s1, in)) (rhs (s2, in))
   implication <- io $ scImplies sc initRelation relationRes
 
+  -- Simplify Term with any theorems that apply
+  (implication', sideConditions) <- applyAllTheorems bthms bc implication
+
   -- Function to prove
   -- forall s1 s2 in out1 out2.
-  --   rel (s1, out1) (s2, out2) -> rel (lhs (s1, in)) (rhs (s2, in))
+  --   srel s1 s2 -> orel (lhs (s1, in)) (rhs (s2, in))
   args <- io $ mapM
     (\(name, t) -> (name,) <$> C.importType sc C.emptyEnv t)
-    [ ("initRhsOutput", outputType)
-    , ("initLhsOutput", outputType)
-    , ("rhsState", rhsStateType)
-    , ("lhsState", lhsStateType)
-    , ("input", lhsInputType) ]
+    [ ("input", bcInputType bc)
+    , ("rhsState", bisimTheoremRhsStateType outerBt)
+    , ("lhsState", bisimTheoremLhsStateType outerBt) ]
+  theorem <- io $ scLambdaList sc args implication'
+
+  tt <- io $ mkTypedTerm sc theorem
+
+  pure $ tt : sideConditions
+
+-- | Build the STATE RELATION THEOREM for 'bc'. See the documentation at the top
+-- of this file for information on the STATE RELATION THEOREM.
+buildStateRelationTheorem :: BisimComponents -> TopLevel TypedTerm
+buildStateRelationTheorem bc = do
+  sc <- getSharedContext
+  let outerBt = bcTheorem bc
+
+  -- Outer function inputs. See comments to the right of each line to see how
+  -- they line up with the documentation at the top of this file.
+  lhsState <- io $ scLocalVar sc 0        -- s1
+  rhsState <- io $ scLocalVar sc 1        -- s2
+  initLhsOutput <- io $ scLocalVar sc 2   -- out1
+  initRhsOutput <- io $ scLocalVar sc 3   -- out2
+
+  -- LHS/RHS initial outputs
+  lhsTuple <- io $ scTuple sc [lhsState, initLhsOutput]  -- (s1, out1)
+  rhsTuple <- io $ scTuple sc [rhsState, initRhsOutput]  -- (s2, out2)
+
+  -- LHS of implication
+  -- orel (s1, out1) (s2, out2)
+  lhsRelation <-
+    scRelation (bisimTheoremOutputRelation outerBt) lhsTuple rhsTuple
+
+  -- RHS of implication
+  -- srel s1 s2
+  rhsRelation <-
+    scRelation (bisimTheoremStateRelation outerBt) lhsState rhsState
+
+  -- lhsRelation implies rhsRelation
+  -- orel (s1, out1) (s2, out2) -> srel s1 s2
+  implication <- io $ scImplies sc lhsRelation rhsRelation
+
+  -- Function to prove
+  -- forall s1 s2 in out1 out2.
+  --   orel (s1, out1) (s2, out2) -> srel s1 s2
+  args <- io $ mapM
+    (\(name, t) -> (name,) <$> C.importType sc C.emptyEnv t)
+    [ ("initRhsOutput", bisimTheoremOutputType outerBt)
+    , ("initLhsOutput", bisimTheoremOutputType outerBt)
+    , ("rhsState", bisimTheoremRhsStateType outerBt)
+    , ("lhsState", bisimTheoremLhsStateType outerBt) ]
   theorem <- io $ scLambdaList sc args implication
 
-  io (mkTypedTerm sc theorem) >>= provePrim script
+  io $ mkTypedTerm sc theorem
+
+-- | Use bisimulation to prove that two terms simulate each other.  Proves the
+-- OUTPUT RELATION THEOREM and STATE RELATION THEOREM for the bisimilar terms.
+-- Additionally enables compositionality through the 'bthms' parameter.  This
+-- function also proves the COMPOSITION SIDE CONDITION for every successfully
+-- applied theorem in 'bthms'.  See the documentation at the top of this file
+-- for the definitions of bisimulation and these theorems.
+proveBisimulation ::
+  ProofScript () ->
+  -- ^ Proof script to use over generated bisimulation term
+  [BisimTheorem] ->
+  -- ^ Previously proved theorems that may be applied
+  TypedTerm ->
+  -- ^ Relation over states for terms to prove bisimilar. Must have type
+  -- @lhsState -> rhsState -> Bit@  Also called the "state relation".
+  TypedTerm ->
+  -- ^ Relation over states and outputs for terms to prove bisimilar. Must have
+  -- type @(lhsState, output) -> (rhsState, output) -> Bit@.  Also called the
+  -- "output relation".
+  TypedTerm ->
+  -- ^ LHS of bisimulation. Must have type
+  -- @(lhsState, input) -> (lhsState, output)@
+  TypedTerm ->
+  -- ^ RHS of bisimulation. Must have type
+  -- @(rhsState, input) -> (rhsState, output)@
+  TopLevel BisimTheorem
+proveBisimulation script bthms srel orel lhs rhs = do
+  -- Typechecking
+  (lhsStateType, rhsStateType, outputType) <- typecheckOutputRelation
+  typecheckStateRelation lhsStateType rhsStateType
+  (lhsName, lhsInputType) <- typecheckSide lhs lhsStateType outputType
+  (rhsName, rhsInputType) <- typecheckSide rhs rhsStateType outputType
+  unless (lhsInputType == rhsInputType) $
+    fail $ unlines [ "Error: Mismatched input types in bisimulation terms."
+                   , "  LHS input type: " ++ C.pretty lhsInputType
+                   , "  RHS input type: " ++ C.pretty rhsInputType ]
+
+  let bt = BisimTheorem
+           { bisimTheoremStateRelation = srel
+           , bisimTheoremOutputRelation = orel
+           , bisimTheoremLhs = lhs
+           , bisimTheoremRhs = rhs
+           , bisimTheoremOutputType = outputType
+           , bisimTheoremLhsStateType = lhsStateType
+           , bisimTheoremRhsStateType = rhsStateType
+           }
+
+  let bc = BisimComponents
+           { bcTheorem = bt
+           , bcInputType = lhsInputType
+           }
+  outputTheorems <- buildOutputRelationTheorem bthms bc
+  stateTheorem <- buildStateRelationTheorem bc
+  proveAll script $ stateTheorem : outputTheorems
+
+  -- Proof succeeded!
+  printOutLnTop Info $ concat [ "Successfully proved bisimulation between "
+                              , Text.unpack lhsName
+                              , " and "
+                              , Text.unpack rhsName ]
+
+  return $ BisimTheorem srel orel lhs rhs outputType lhsStateType rhsStateType
 
   where
-    -- Typecheck relation. The expected type for a relation is:
+    -- Typecheck output relation. The expected type is:
     -- @(lhsStateType, outputType) -> (rhsStateType, outputType) -> Bit@
     --
-    -- If the relation typechecks, 'typecheckRelation' evaluates to a tuple of:
+    -- If the relation typechecks, 'typecheckOutputRelation' evaluates to a
+    -- tuple of:
     -- @(lhsStateType, rhsStateType, outputType)@
     -- Otherwise, this invokes 'fail' with a description of the specific
     -- typechecking error.
-    typecheckRelation :: TopLevel (C.Type, C.Type, C.Type)
-    typecheckRelation =
-      case ttType relation of
+    typecheckOutputRelation :: TopLevel (C.Type, C.Type, C.Type)
+    typecheckOutputRelation =
+      case ttType orel of
         TypedTermSchema
           (C.Forall
             []
@@ -168,18 +599,51 @@ proveBisimulation script relation lhs rhs = do
             , "RHS output type: " ++ C.pretty o2 ]
 
           return (s1, s2, o1)
-        _ -> fail $ "Error: Unexpected relation type: "
-                 ++ show (ppTypedTermType (ttType relation))
+        _ -> fail $ "Error: Unexpected output relation type: "
+                 ++ show (ppTypedTermType (ttType orel))
+
+    -- Check that 'lhsStateType' and 'rhsStateType' match the extracted types
+    -- from 'typecheckOutputRelation'.  Invokes 'fail' if the types do not
+    -- match.
+    typecheckStateRelation :: C.Type -> C.Type -> TopLevel ()
+    typecheckStateRelation lhsStateType rhsStateType =
+      case ttType srel of
+        TypedTermSchema
+          (C.Forall
+            []
+            []
+            (C.TCon
+              (C.TC C.TCFun)
+              [ s1
+              , C.TCon
+                (C.TC C.TCFun)
+                [ s2, C.TCon (C.TC C.TCBit) []]])) -> do
+          unless (s1 == lhsStateType) $ fail $ unlines
+            [ "Error: LHS of state relation and output relations have incompatible state types:"
+            , "  State relation LHS state type: " ++ C.pretty s1
+            , "  Output relation LHS state type: " ++ C.pretty lhsStateType ]
+
+          unless (s2 == rhsStateType) $ fail $ unlines
+            [ "Error: RHS of state relation and output relations have incompatible state types:"
+            , "  State relation RHS state type: " ++ C.pretty s2
+            , "  Output relation RHS state type: " ++ C.pretty rhsStateType ]
+        _ -> fail $ "Error: Unexpected state relation type: "
+                 ++ show (ppTypedTermType (ttType srel))
 
     -- Typecheck bisimulation term. The expected type for a bisimulation term
     -- is:
     -- @(stateType, inputType) -> (stateType, outputType)@
     --
-    -- If the term typechecks, this function returns @inputType@.  Otherwise,
-    -- this funciton invokes 'fail' with a description of the specific
-    -- typechecking error.
-    typecheckSide :: TypedTerm -> C.Type -> C.Type -> TopLevel C.Type
-    typecheckSide side stateType outputType =
+    -- If the term typechecks, this function returns a pair containing the name
+    -- of the 'Constant' and @inputType@.  Otherwise, this function invokes
+    -- 'fail' with a description of the specific typechecking error.
+    typecheckSide
+      :: TypedTerm -> C.Type -> C.Type -> TopLevel (Text.Text, C.Type)
+    typecheckSide side stateType outputType = do
+      -- Check that 'side' is a 'Constant' (necessary for composition)
+      name <- constantName $ unwrapTermF $ ttTerm side
+
+      -- Check function arguments
       case ttType side of
         TypedTermSchema
           (C.Forall
@@ -201,10 +665,10 @@ proveBisimulation script relation lhs rhs = do
 
           unless (o == outputType) $ fail $ unlines
             [ "Error: Output type in bisimulation term does not match output type in relation."
-            , "  Expected: " ++ C.pretty outputType
+            ,"  Expected: " ++ C.pretty outputType
             , "  Actual: " ++ C.pretty o ]
 
-          return i
+          return (name, i)
         _ ->
           let stStr = C.pretty stateType in
           fail $ unlines
@@ -212,7 +676,79 @@ proveBisimulation script relation lhs rhs = do
           , "  Expected: (" ++ stStr ++ ", inputType) -> (" ++ stStr ++ ", outputType)"
           , "  Actual: " ++ show (ppTypedTermType (ttType side)) ]
 
-    -- Generate 'Term' for application of a relation
-    scRelation :: SharedContext -> Term -> Term -> TopLevel Term
-    scRelation sc relLhs relRhs = io $
-      scApplyAll sc (ttTerm relation) [relLhs, relRhs]
+-- | Replace the invocation of a specific 'Constant' with an 'ExtCns'.  The
+-- function returns the resulting 'Term' and updates a 'ReplaceState' to hold
+-- the generated 'ExtCns' and the specific 'App' that was replaced.
+replaceConstantTerm :: TypedTerm
+                    -- ^ 'Constant' to replace application of
+                    -> C.Type
+                    -- ^ 'constant's return type
+                    -> Term
+                    -- ^ 'Term' to perform replacement in
+                    -> State.StateT ReplaceState TopLevel Term
+replaceConstantTerm constant constantRetType term = do
+  sc <- lift getSharedContext
+  case term of
+    Unshared termF -> do
+      termF' <- replaceConstantTermF termF
+      liftIO $ scTermF sc termF'
+    STApp{ stAppIndex = i, stAppTermF = termF } -> do
+      table <- State.gets rsMemo
+      case Map.lookup i table of
+        Just x -> return x
+        Nothing -> do
+          termF' <- replaceConstantTermF termF
+          term' <- liftIO $ scTermF sc termF'
+          let table' = Map.insert i term' table
+          State.modify$ \st -> st { rsMemo = table' }
+          return term'
+  where
+    -- | Partner function to 'replaceConstantTerm' that operates over 'TermF's.
+    replaceConstantTermF
+      :: TermF Term -> State.StateT ReplaceState TopLevel (TermF Term)
+    replaceConstantTermF termF = do
+      case termF of
+        App x _ | unwrapTermF x == unwrapTermF (ttTerm constant) ->
+          State.gets rsExtCns >>= \case
+            Just v ->
+              State.gets rsApp >>= \case
+                Just a | a == termF ->
+                  -- Encountered another call to the function under replacement
+                  -- that matches the replaced function.  This can happen even
+                  -- when the underlying Cryptol does not explicitly make
+                  -- multiple calls because translation to SAWCore can insert
+                  -- additional function calls with the same arguments.  In this
+                  -- case, simply return the same 'ExtCns' already generated.
+                  return $ unwrapTermF v
+                Just _ -> do
+                  -- Encountered a call to the function under replacement with
+                  -- different arguments.  This isn't yet supported.
+                  name <- Text.unpack <$> lift (constantName (unwrapTermF x))
+                  fail $ concat ["Error: Encountered multiple calls to "
+                                , name
+                                , ". Use of composition with multiple "
+                                , "subfunction calls is not yet supported."
+                                ]
+                _ -> panic "replaceConstantTermF"
+                           ["rsApp should always exist when rsExtCns exists"]
+            Nothing -> do
+              sc <- lift getSharedContext
+
+              -- Generate an 'ExtCns' and return it, thereby replacing 'termF'
+              -- with it.
+              tp <- liftIO $ C.importType sc C.emptyEnv constantRetType
+              name <- lift $ constantName $ unwrapTermF x
+              v <- liftIO $ scFreshGlobal sc name tp
+              State.modify $ \st -> st { rsExtCns = Just v, rsApp = Just termF }
+              return $ unwrapTermF v
+        _ ->
+          -- Recurse
+          mapM (replaceConstantTerm constant constantRetType) termF
+
+-- Extract the name from a 'Constant'. Fails if provided another kind of 'TermF'
+constantName :: TermF Term -> TopLevel Text.Text
+constantName (Constant e _) = return $ toShortName $ ecName e
+constantName tf = do
+  sc <- getSharedContext
+  term <- io $ scTermF sc tf
+  fail $ "Error: Expected a constant, but got: " ++ showTerm term
diff --git a/src/SAWScript/Bisimulation/BisimTheorem.hs b/src/SAWScript/Bisimulation/BisimTheorem.hs
new file mode 100644
index 0000000000..fd95b5dedd
--- /dev/null
+++ b/src/SAWScript/Bisimulation/BisimTheorem.hs
@@ -0,0 +1,32 @@
+{- |
+Module      : SAWScript.Bisimulation.BisimTheorem
+License     : BSD3
+Maintainer  : bboston7
+Stability   : experimental
+-}
+
+module SAWScript.Bisimulation.BisimTheorem
+  ( BisimTheorem(..) )
+  where
+
+import qualified Cryptol.TypeCheck.Type as C
+import Verifier.SAW.TypedTerm ( TypedTerm )
+
+-- | A proved bisimulation theorem.  See the comment at the top of
+-- "SAWScript.Bisimulation" for an explanation of some of the terms used here.
+data BisimTheorem = BisimTheorem {
+    bisimTheoremStateRelation :: TypedTerm
+ -- ^ State relation
+  , bisimTheoremOutputRelation :: TypedTerm
+ -- ^ Output relation
+  , bisimTheoremLhs :: TypedTerm
+ -- ^ Left hand side of the bisimulation
+  , bisimTheoremRhs :: TypedTerm
+ -- ^ Right hand side of the bisimulation
+  , bisimTheoremOutputType :: C.Type
+ -- ^ Output type for the bisimilar terms
+  , bisimTheoremLhsStateType :: C.Type
+ -- ^ State type for the left term
+  , bisimTheoremRhsStateType :: C.Type
+ -- ^ State type for the right term
+  }
diff --git a/src/SAWScript/Interpreter.hs b/src/SAWScript/Interpreter.hs
index ea80c7a308..619605a321 100644
--- a/src/SAWScript/Interpreter.hs
+++ b/src/SAWScript/Interpreter.hs
@@ -1665,22 +1665,29 @@ primitives = Map.fromList
     , "successful, and aborts if unsuccessful."
     ]
 
-  , prim "prove_bisim"         "ProofScript () -> Term -> Term -> Term -> TopLevel ProofResult"
+  , prim "prove_bisim"         "ProofScript () -> [BisimTheorem] -> Term -> Term -> Term -> Term -> TopLevel BisimTheorem"
     (pureVal proveBisimulation)
     Experimental
-    [ "Use bisimulation to prove that two terms simulate each other.  The first"
-    , "argument is the proof strategy to use.  The second argument is a"
-    , "relation over the states and outputs for the third and fourth"
-    , "arguments. The relation must have the type"
-    , "'(lhsState, output) -> (rhsState, output) -> Bit'. The third and fourth"
-    , "arguments are the two terms to prove bisimilar. They must have the types"
-    , "'(lhsState, input) -> (lhsState, output)' and"
-    , "'(rhsState, input) -> (rhsState, output)' respectively."
+    [ "Use bisimulation to prove that two terms simulate each other.  The "
+    , "command takes the following arguments: "
+    , "1. The proof strategy to use"
+    , "2. A list of already proven bisimulation theorems"
+    , "3. A state relation `srel : lhsState -> rhsState -> Bit`"
+    , "4. An output relation `orel : (lhsState, output) -> (rhsState, output) -> Bit`"
+    , "5. A term `lhs : (lhsState, input) -> (lhsState, output)`"
+    , "6. A term `rhs : (rhsState, input) -> (rhsState, output)`"
+    , "and considers `lhs` and `rhs` bisimilar when the following two theorems hold:"
+    , "* OUTPUT RELATION THEOREM:"
+    , "   forall s1 s2 in."
+    , "     srel s1 s2 -> orel (lhs (s1, in)) (rhs (s2, in))"
+    , "* STATE RELATION THEOREM:"
+    , "   forall s1 s2 out1 out2."
+    , "     orel (s1, out1) (s2, out2) -> srel s1 s2"
     , ""
-    , "Let the second argument be called 'rel', the third 'lhs', and the"
-    , "fourth 'rhs'. The prover considers 'lhs' and 'rhs' bisimilar when:"
-    , "  forall s1 s2 in out1 out2."
-    , "    rel (s1, out1) (s2, out2) -> rel (lhs (s1, in)) (rhs (s2, in))"
+    , "LIMITATIONS: For now, the prove_bisim command has a couple limitations:"
+    , "* `lhs` and `rhs` (arguments 5 and 6) must be named functions."
+    , "* Each subterm present in the list of bisimulation theorems already"
+    , "  proven (argument 2) may be invoked at most once in `lhs` or `rhs`."
     ]
 
   , prim "sat"                 "ProofScript () -> Term -> TopLevel SatResult"
diff --git a/src/SAWScript/Value.hs b/src/SAWScript/Value.hs
index 25a8c46421..7875521e1e 100644
--- a/src/SAWScript/Value.hs
+++ b/src/SAWScript/Value.hs
@@ -62,6 +62,7 @@ import qualified Prettyprinter as PP
 import qualified Data.AIG as AIG
 
 import qualified SAWScript.AST as SS
+import SAWScript.Bisimulation.BisimTheorem (BisimTheorem)
 import qualified SAWScript.Exceptions as SS
 import qualified SAWScript.Position as SS
 import qualified SAWScript.Crucible.Common as Common
@@ -148,6 +149,7 @@ data Value
   | VSimpset SAWSimpset
   | VRefnset SAWRefnset
   | VTheorem Theorem
+  | VBisimTheorem BisimTheorem
   -----
   | VLLVMCrucibleSetup !(LLVMCrucibleSetupM Value)
   | VLLVMCrucibleMethodSpec (CMSLLVM.SomeLLVM CMS.ProvedSpec)
@@ -361,6 +363,7 @@ showsPrecValue opts nenv p v =
     VTheorem thm ->
       showString "Theorem " .
       showParen True (showString (prettyProp opts' nenv (thmProp thm)))
+    VBisimTheorem _ -> showString "<<Bisimulation theorem>>"
     VLLVMCrucibleSetup{} -> showString "<<Crucible Setup>>"
     VLLVMCrucibleSetupValue{} -> showString "<<Crucible SetupValue>>"
     VLLVMCrucibleMethodSpec{} -> showString "<<Crucible MethodSpec>>"
@@ -1263,6 +1266,13 @@ instance FromValue Theorem where
     fromValue (VTheorem t) = t
     fromValue _ = error "fromValue Theorem"
 
+instance IsValue BisimTheorem where
+    toValue = VBisimTheorem
+
+instance FromValue BisimTheorem where
+    fromValue (VBisimTheorem t) = t
+    fromValue _ = error "fromValue BisimTheorem"
+
 instance IsValue JavaType where
     toValue t = VJavaType t