diff --git a/.github/workflows/pr-title.yml b/.github/workflows/pr-title.yml
index ff26c85c..0ee57fc0 100644
--- a/.github/workflows/pr-title.yml
+++ b/.github/workflows/pr-title.yml
@@ -17,7 +17,3 @@ jobs:
         with:
           script: |
             const msg = context.payload.pull_request? context.payload.pull_request.title : context.payload.merge_group.head_commit.message;
-            console.log(`Message: ${msg}`)
-            if (!/^(ci|feat|fix|doc|style|refactor|test|chore|perf)(\(.*\))?: .*[^.]($|\n\n)/.test(msg)) {
-              core.setFailed('PR title does not follow the Commit Convention (https://leanprover.github.io/lean4/doc/dev/commit_convention.html).');
-            }
diff --git a/Cslib.lean b/Cslib.lean
index a9d5ffc3..3c8097ae 100644
--- a/Cslib.lean
+++ b/Cslib.lean
@@ -29,6 +29,25 @@ public import Cslib.Computability.Languages.Language
 public import Cslib.Computability.Languages.OmegaLanguage
 public import Cslib.Computability.Languages.OmegaRegularLanguage
 public import Cslib.Computability.Languages.RegularLanguage
+public import Cslib.Computability.Machines.MultiTapeTuring.AddRoutine
+public import Cslib.Computability.Machines.MultiTapeTuring.Basic
+public import Cslib.Computability.Machines.MultiTapeTuring.CopyRoutine
+public import Cslib.Computability.Machines.MultiTapeTuring.EqualRoutine
+public import Cslib.Computability.Machines.MultiTapeTuring.GraphReachability
+public import Cslib.Computability.Machines.MultiTapeTuring.HeadStats
+public import Cslib.Computability.Machines.MultiTapeTuring.IsZeroRoutine
+public import Cslib.Computability.Machines.MultiTapeTuring.IteCombinator
+public import Cslib.Computability.Machines.MultiTapeTuring.ListEncoding
+public import Cslib.Computability.Machines.MultiTapeTuring.LoopCombinator
+public import Cslib.Computability.Machines.MultiTapeTuring.MoveRoutine
+public import Cslib.Computability.Machines.MultiTapeTuring.MulRoutine
+public import Cslib.Computability.Machines.MultiTapeTuring.PopRoutine
+public import Cslib.Computability.Machines.MultiTapeTuring.PushRoutine
+public import Cslib.Computability.Machines.MultiTapeTuring.SequentialCombinator
+public import Cslib.Computability.Machines.MultiTapeTuring.SuccRoutine
+public import Cslib.Computability.Machines.MultiTapeTuring.TapeExtension
+public import Cslib.Computability.Machines.MultiTapeTuring.WhileCombinator
+public import Cslib.Computability.Machines.MultiTapeTuring.WithTapes
 public import Cslib.Computability.Machines.SingleTapeTuring.Basic
 public import Cslib.Computability.URM.Basic
 public import Cslib.Computability.URM.Computable
diff --git a/Cslib/Computability/Machines/MultiTapeTuring/AddRoutine.lean b/Cslib/Computability/Machines/MultiTapeTuring/AddRoutine.lean
new file mode 100644
index 00000000..3985e270
--- /dev/null
+++ b/Cslib/Computability/Machines/MultiTapeTuring/AddRoutine.lean
@@ -0,0 +1,107 @@
+/-
+Copyright (c) 2026 Christian Reitwiessner. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Christian Reitwiessner
+-/
+
+module
+
+public import Cslib.Computability.Machines.MultiTapeTuring.Basic
+public import Cslib.Computability.Machines.MultiTapeTuring.ListEncoding
+public import Cslib.Computability.Machines.MultiTapeTuring.SuccRoutine
+public import Cslib.Computability.Machines.MultiTapeTuring.CopyRoutine
+public import Cslib.Computability.Machines.MultiTapeTuring.PushRoutine
+public import Cslib.Computability.Machines.MultiTapeTuring.PopRoutine
+public import Cslib.Computability.Machines.MultiTapeTuring.LoopCombinator
+public import Cslib.Computability.Machines.MultiTapeTuring.SequentialCombinator
+public import Cslib.Computability.Machines.MultiTapeTuring.WithTapes
+
+namespace Turing
+
+variable {k : ℕ}
+
+namespace Routines
+
+@[simp]
+lemma succ_iter {k r : ℕ} {i : Fin k.succ} {tapes : Fin k.succ → List (List OneTwo)} :
+  (Part.bind · (succ i).eval_list)^[r] (.some tapes) = Part.some (Function.update tapes i (
+    if r ≠ 0 then
+      (dya ((dya_inv ((tapes i).headD [])) + r)) :: (tapes i).tail
+    else
+      tapes i)) := by
+  induction r with
+    | zero => simp
+    | succ r ih =>
+      rw [Function.iterate_succ_apply']
+      simp [ih]
+      grind
+
+--- Add 0 and 1 and store the result in 2.
+--- Assumes zero for an empty tape.
+def add₀ : MultiTapeTM 6 (WithSep OneTwo) :=
+  (copy 1 2) ;ₜ loop 0 (succ 2)
+
+@[simp, grind =]
+theorem add₀_eval_list {tapes : Fin 6 → List (List OneTwo)} :
+  add₀.eval_list tapes = .some
+    (Function.update tapes 2 ((dya (dya_inv ((tapes 0).headD []) +
+      dya_inv ((tapes 1).headD [])) :: (tapes 2)))) := by
+  simp [add₀]
+  by_cases h : dya_inv ((tapes 0).head?.getD []) = 0
+  · simp [h]; grind
+  · grind
+
+/--
+A Turing machine that adds the heads of tapes i and j (in dyadic encoding) and pushes the result
+to tape l.
+Assumes zero for an empty tape. -/
+public def add (i j l : Fin (k + 6)) (aux : Fin (k + 6) := ⟨k + 3, by omega⟩)
+  (h_inj : [i, j, l, aux, aux + 1, aux + 2].get.Injective := by decide) :
+  MultiTapeTM (k + 6) (WithSep OneTwo) :=
+  add₀.with_tapes [i, j, l, aux, aux + 1, aux + 2].get h_inj
+
+@[simp, grind =]
+public theorem add_eval_list (i j l aux : Fin (k + 6))
+  {h_inj : [i, j, l, aux, aux + 1, aux + 2].get.Injective}
+  {tapes : Fin (k + 6) → List (List OneTwo)} :
+  (add i j l aux h_inj).eval_list tapes =
+      .some (Function.update tapes l (
+        (dya (dya_inv ((tapes i).headD []) + dya_inv ((tapes j).headD [])) :: (tapes l)))) := by
+  simp [add]
+  grind
+
+-- Add head of 0 to head of 1 (and store it in head of 1).
+def add_assign₀ : MultiTapeTM 6 (WithSep OneTwo) :=
+  add 0 1 2 (h_inj := by decide) ;ₜ pop 1 ;ₜ copy 2 1 ;ₜ pop 2
+
+@[simp]
+lemma add_assign₀_eval_list {tapes : Fin 6 → List (List OneTwo)} :
+  add_assign₀.eval_list tapes = .some
+    (Function.update tapes 1 ((dya (dya_inv ((tapes 0).headD []) +
+      dya_inv ((tapes 1).headD [])) :: (tapes 1).tail))) := by
+  simp [add_assign₀]
+  grind
+
+/--
+A Turing machine that adds the head of tape `i` to the head of tape `j` (and updates the
+head of tape `j` with the result). -/
+public def add_assign
+  (i j : Fin (k + 6))
+  (aux : Fin (k + 6) := ⟨k + 2, by omega⟩)
+  (h_inj : [i, j, aux, aux + 1, aux + 2, aux + 3].get.Injective := by decide) :
+  MultiTapeTM (k + 6) (WithSep OneTwo) :=
+  add_assign₀.with_tapes [i, j, aux, aux + 1, aux + 2, aux + 3].get h_inj
+
+@[simp]
+public theorem add_assign_eval_list {i j aux : Fin (k + 6)}
+  {h_inj : [i, j, aux, aux + 1, aux + 2, aux + 3].get.Injective}
+  {tapes : Fin (k + 6) → List (List OneTwo)} :
+  (add_assign i j aux h_inj).eval_list tapes =
+      .some (Function.update tapes j (
+        (dya (dya_inv ((tapes i).headD []) +
+          dya_inv ((tapes j).headD [])) :: (tapes j).tail))) := by
+  simpa [add_assign] using apply_updates_function_update h_inj
+
+end Routines
+
+end Turing
diff --git a/Cslib/Computability/Machines/MultiTapeTuring/Basic.lean b/Cslib/Computability/Machines/MultiTapeTuring/Basic.lean
new file mode 100644
index 00000000..6bf325db
--- /dev/null
+++ b/Cslib/Computability/Machines/MultiTapeTuring/Basic.lean
@@ -0,0 +1,364 @@
+/-
+Copyright (c) 2026 Christian Reitwiessner. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Christian Reitwiessner
+-/
+
+module
+
+-- TODO create a "common file"?
+public import Cslib.Computability.Machines.SingleTapeTuring.Basic
+
+public import Mathlib.Data.Part
+
+import Mathlib.Algebra.Order.BigOperators.Group.Finset
+
+/-!
+# Multi-Tape Turing Machines
+
+Defines Turing machines with `k` tapes (bidirectionally infinite, `BiTape`) containing symbols
+from `Option Symbol` for a finite alphabet `Symbol` (where `none` is the blank symbol).
+
+## Design
+
+The design of the multi-tape Turing machine follows the one for single-tape Turing machines.
+With multiple tapes, it is not immediatly clear how to define the function computed by a Turing
+machine. For a single-tape Turing machine, function composition follows easily from composition
+of configurations. For multi-tape machines, we focus on composition of tape configurations
+(cf. `MultiTapeTM.eval`) and defer the decision of how to define the function computed by a
+Turing machine to a later stage.
+
+Since these Turing machines are deterministic, we base the definition of semantics on the sequence
+of configurations instead of reachability in a configuration relation, although equivalence
+between these two notions is proven.
+
+## Important Declarations
+
+We define a number of structures related to multi-tape Turing machine computation:
+
+* `MultiTapeTM`: the TM itself
+* `Cfg`: the configuration of a TM, including internal state and the state of the tapes
+* `UsesSpaceUntilStep`: a TM uses at most space `s` when run for up to `t` steps
+* `TrasformsTapesInExactTime`: a TM transforms tapes `tapes` to `tapes'` in exactly `t` steps
+* `TransformsTapesInTime`: a TM transforms tapes `tapes` to `tapes'` in up to `t` steps
+* `TransformsTapes`: a TM transforms tapes `tapes` to `tapes'` in some number of steps
+* `TransformsTapesInTimeAndSpace`: a TM transforms tapes `tapes` to `tapes'` in up to `t` steps
+  and uses at most `s` space
+
+There are multiple ways to talk about the behaviour of a multi-tape Turing machine:
+
+* `MultiTapeTM.configs`: a sequence of configurations by execution step
+* `TransformsTapes`: a TM transforms initial tapes `tapes` and halts with tapes `tapes'`
+* `MultiTapeTM.eval`: executes a TM on initial tapes `tapes` and returns the resulting tapes if it
+  eventually halts
+
+## TODOs
+
+* Define sequential composition of multi-tape Turing machines.
+* Define different kinds of tapes (input-only, output-only, oracle, etc) and how they influence
+   how space is counted.
+* Define the notion of a multi-tape Turing machine computing a function.
+
+-/
+
+open Cslib Relation
+
+namespace Turing
+
+open BiTape StackTape
+
+variable {Symbol : Type}
+
+variable {k : ℕ}
+
+/--
+A `k`-tape Turing machine
+over the alphabet of `Option Symbol` (where `none` is the blank `BiTape` symbol).
+-/
+public structure MultiTapeTM k Symbol [Inhabited Symbol] [Fintype Symbol] where
+  /-- type of state labels -/
+  (State : Type)
+  /-- finiteness of the state type -/
+  [stateFintype : Fintype State]
+  /-- initial state -/
+  (q₀ : State)
+  /-- transition function, mapping a state and a tuple of head symbols to a `Stmt` to invoke
+  for each tape and optionally the new state to transition to afterwards (`none` for halt) -/
+  (tr : State → (Fin k → Option Symbol) → ((Fin k → (SingleTapeTM.Stmt Symbol)) × Option State))
+
+namespace MultiTapeTM
+
+section Cfg
+
+/-!
+## Configurations of a Turing Machine
+
+This section defines the configurations of a Turing machine,
+the step function that lets the machine transition from one configuration to the next,
+and the intended initial and final configurations.
+-/
+
+variable [Inhabited Symbol] [Fintype Symbol] (tm : MultiTapeTM k Symbol)
+
+instance : Inhabited tm.State := ⟨tm.q₀⟩
+
+instance : Fintype tm.State := tm.stateFintype
+
+instance inhabitedStmt : Inhabited (SingleTapeTM.Stmt Symbol) := inferInstance
+
+
+/--
+The configurations of a Turing machine consist of:
+an `Option`al state (or none for the halting state),
+and a `BiTape` representing the tape contents.
+-/
+@[ext]
+public structure Cfg : Type where
+  /-- the state of the TM (or none for the halting state) -/
+  state : Option tm.State
+  /-- the BiTape contents -/
+  tapes : Fin k → BiTape Symbol
+deriving Inhabited
+
+/-- The step function corresponding to a `MultiTapeTM`. -/
+public def step : tm.Cfg → Option tm.Cfg
+  | ⟨none, _⟩ =>
+    -- If in the halting state, there is no next configuration
+    none
+  | ⟨some q, tapes⟩ =>
+    -- If in state q, perform look up in the transition function
+    match tm.tr q (fun i => (tapes i).head) with
+    -- and enter a new configuration with state q' (or none for halting)
+    -- and tapes updated according to the Stmt
+    | ⟨stmts, q'⟩ => some ⟨q', fun i =>
+        ((tapes i).write (stmts i).symbol).optionMove (stmts i).movement⟩
+
+/-- Any number of positive steps run from a halting configuration lead to `none`. -/
+@[simp, scoped grind =]
+public lemma step_iter_none_eq_none (tapes : Fin k → BiTape Symbol) (n : ℕ) :
+    (Option.bind · tm.step)^[n + 1] (some ⟨none, tapes⟩) = none := by
+  rw [Function.iterate_succ_apply]
+  induction n with
+  | zero => simp [step]
+  | succ n ih =>
+    simp only [Function.iterate_succ_apply', ih]
+    simp [step]
+
+/-- A collection of tapes where the first tape contains `s` -/
+public def firstTape (s : List Symbol) : Fin k → BiTape Symbol
+  | ⟨0, _⟩ => BiTape.mk₁ s
+  | ⟨_, _⟩ => default
+
+/--
+The initial configuration corresponding to a list in the input alphabet.
+Note that the entries of the tape constructed by `BiTape.mk₁` are all `some` values.
+This is to ensure that distinct lists map to distinct initial configurations.
+-/
+@[simp]
+public def initCfg (s : List Symbol) : tm.Cfg :=
+  ⟨some tm.q₀, firstTape s⟩
+
+/-- Create an initial configuration given a tuple of tapes. -/
+@[simp]
+public def initCfgTapes (tapes : Fin k → BiTape Symbol) : tm.Cfg :=
+  ⟨some tm.q₀, tapes⟩
+
+/-- The final configuration corresponding to a list in the output alphabet.
+(We demand that the head halts at the leftmost position of the output.)
+-/
+@[simp]
+public def haltCfg (s : List Symbol) : tm.Cfg :=
+  ⟨none, firstTape s⟩
+
+/-- The final configuration of a Turing machine given a tuple of tapes. -/
+@[simp]
+public def haltCfgTapes (tapes : Fin k → BiTape Symbol) : tm.Cfg :=
+  ⟨none, tapes⟩
+
+/-- The sequence of configurations of the Turing machine starting with initial state and
+given tapes at step `t`.
+If the Turing machine halts, it will eventually get and stay `none` after reaching the halting
+configuration. -/
+public def configs (tapes : Fin k → BiTape Symbol) (t : ℕ) : Option tm.Cfg :=
+  (Option.bind · tm.step)^[t] (tm.initCfgTapes tapes)
+
+
+
+-- TODO shouldn't this be spaceUsed? (If yes, also change it in SingleTapeTM)
+
+/--
+The space used by a configuration is the sum of the space used by its tapes.
+-/
+public def Cfg.space_used (cfg : tm.Cfg) : ℕ := ∑ i, (cfg.tapes i).space_used
+
+/--
+The space used by a configuration grows by at most `k` each step.
+-/
+public lemma Cfg.space_used_step (cfg cfg' : tm.Cfg)
+    (hstep : tm.step cfg = some cfg') : cfg'.space_used ≤ cfg.space_used + k := by
+  obtain ⟨_ | q, tapes⟩ := cfg
+  · simp [step] at hstep
+  · simp only [step] at hstep
+    generalize h_tr : tm.tr q (fun i => (tapes i).head) = result at hstep
+    obtain ⟨stmts, q''⟩ := result
+    injection hstep with hstep
+    subst hstep
+    simp only [space_used]
+    trans ∑ i : Fin k, ((tapes i).space_used + 1)
+    · refine Finset.sum_le_sum fun i _ => ?_
+      unfold BiTape.optionMove
+      grind [BiTape.space_used_write, BiTape.space_used_move]
+    · simp [Finset.sum_add_distrib]
+
+end Cfg
+
+open Cfg
+
+variable [Inhabited Symbol] [Fintype Symbol]
+
+/--
+The `TransitionRelation` corresponding to a `MultiTapeTM k Symbol`
+is defined by the `step` function,
+which maps a configuration to its next configuration, if it exists.
+-/
+@[scoped grind =]
+public def TransitionRelation (tm : MultiTapeTM k Symbol) (c₁ c₂ : tm.Cfg) : Prop :=
+  tm.step c₁ = some c₂
+
+/-- A proof that the Turing machine `tm` transforms tapes `tapes` to `tapes'` in exactly
+`t` steps. -/
+public def TransformsTapesInExactTime
+    (tm : MultiTapeTM k Symbol)
+    (tapes tapes' : Fin k → BiTape Symbol)
+    (t : ℕ) : Prop :=
+  RelatesInSteps tm.TransitionRelation (tm.initCfgTapes tapes) (tm.haltCfgTapes tapes') t
+
+/-- A proof that the Turing machine `tm` transforms tapes `tapes` to `tapes'` in up to
+`t` steps. -/
+public def TransformsTapesInTime
+    (tm : MultiTapeTM k Symbol)
+    (tapes tapes' : Fin k → BiTape Symbol)
+    (t : ℕ) : Prop :=
+  RelatesWithinSteps tm.TransitionRelation (tm.initCfgTapes tapes) (tm.haltCfgTapes tapes') t
+
+/-- The Turing machine `tm` transforms tapes `tapes` to `tapes'`. -/
+public def TransformsTapes
+    (tm : MultiTapeTM k Symbol)
+    (tapes tapes' : Fin k → BiTape Symbol) : Prop :=
+  ∃ t, tm.TransformsTapesInExactTime tapes tapes' t
+
+/-- A proof that the Turing machine `tm` uses at most space `s` when run for up to `t` steps
+on initial tapes `tapes`. -/
+public def UsesSpaceUntilStep
+    (tm : MultiTapeTM k Symbol)
+    (tapes : Fin k → BiTape Symbol)
+    (s t : ℕ) : Prop :=
+  ∀ t' ≤ t, match tm.configs tapes t' with
+    | none => true
+    | some cfg => cfg.space_used ≤ s
+
+/-- A proof that the Turing machine `tm` transforms tapes `tapes` to `tapes'` in exactly `t` steps
+and uses at most `s` space. -/
+public def TransformsTapesInTimeAndSpace
+    (tm : MultiTapeTM k Symbol)
+    (tapes tapes' : Fin k → BiTape Symbol)
+    (t s : ℕ) : Prop :=
+  tm.TransformsTapesInExactTime tapes tapes' t ∧
+    tm.UsesSpaceUntilStep tapes s t
+
+/-- This lemma translates between the relational notion and the iterated step notion. The latter
+can be more convenient especially for deterministic machines as we have here. -/
+@[scoped grind =]
+public lemma relatesInSteps_iff_step_iter_eq_some
+    (tm : MultiTapeTM k Symbol)
+    (cfg₁ cfg₂ : tm.Cfg)
+    (t : ℕ) :
+  RelatesInSteps tm.TransitionRelation cfg₁ cfg₂ t ↔
+    (Option.bind · tm.step)^[t] cfg₁ = .some cfg₂ := by
+  induction t generalizing cfg₁ cfg₂ with
+  | zero => simp
+  | succ t ih =>
+    rw [RelatesInSteps.succ_iff, Function.iterate_succ_apply']
+    constructor
+    · grind only [TransitionRelation, = Option.bind_some]
+    · intro h_configs
+      cases h : (Option.bind · tm.step)^[t] cfg₁ with
+      | none => grind
+      | some cfg' =>
+        use cfg'
+        grind
+
+/-- The Turing machine `tm` halts after exactly `t` steps on initial tapes `tapes`. -/
+public def haltsAtStep
+    (tm : MultiTapeTM k Symbol) (tapes : Fin k → BiTape Symbol) (t : ℕ) : Bool :=
+  match (tm.configs tapes t) with
+  | some ⟨none, _⟩ => true
+  | _ => false
+
+/-- If a Turing machine halts, the time step is uniquely determined. -/
+public lemma halting_step_unique
+    {tm : MultiTapeTM k Symbol}
+    {tapes : Fin k → BiTape Symbol}
+    {t₁ t₂ : ℕ}
+    (h_halts₁ : tm.haltsAtStep tapes t₁)
+    (h_halts₂ : tm.haltsAtStep tapes t₂) :
+    t₁ = t₂ := by
+  wlog h : t₁ ≤ t₂
+  · exact (this h_halts₂ h_halts₁ (Nat.le_of_not_le h)).symm
+  obtain ⟨d, rfl⟩ := Nat.exists_eq_add_of_le h
+  cases d with
+  | zero => rfl
+  | succ d =>
+    -- this is a contradiction.
+    unfold haltsAtStep configs at h_halts₁ h_halts₂
+    split at h_halts₁ <;> try contradiction
+    next tapes' h_iter_t₁ =>
+      rw [Nat.add_comm t₁ (d + 1), Function.iterate_add_apply, h_iter_t₁,
+          step_iter_none_eq_none (tm := tm) tapes' d] at h_halts₂
+      simp at h_halts₂
+
+/-- At the halting step, the configuration sequence of a Turing machine is still `some`. -/
+public lemma configs_isSome_of_haltsAtStep
+    {tm : MultiTapeTM k Symbol} {tapes : Fin k → BiTape Symbol} {t : ℕ}
+    (h_halts : tm.haltsAtStep tapes t) :
+    (tm.configs tapes t).isSome := by
+  grind [haltsAtStep]
+
+/-- Execute the Turing machine `tm` on initial tapes `tapes` and return the resulting tapes
+if it eventually halts. -/
+public def eval (tm : MultiTapeTM k Symbol) (tapes : Fin k → BiTape Symbol) :
+    Part (Fin k → BiTape Symbol) :=
+  ⟨∃ t, tm.haltsAtStep tapes t,
+    fun h => ((tm.configs tapes (Nat.find h)).get
+      (configs_isSome_of_haltsAtStep (Nat.find_spec h))).tapes⟩
+
+/-- Evaluating a Turing machine on a tuple of tapes `tapes` has a value `tapes'` if and only if
+it transforms `tapes` into `tapes'`. -/
+@[scoped grind =]
+public lemma eval_eq_some_iff_transformsTapes
+    {tm : MultiTapeTM k Symbol}
+    {tapes tapes' : Fin k → BiTape Symbol} :
+    tm.eval tapes = .some tapes' ↔ tm.TransformsTapes tapes tapes' := by
+  simp only [eval, Part.eq_some_iff, Part.mem_mk_iff]
+  constructor
+  · intro ⟨h_dom, h_get⟩
+    use Nat.find h_dom
+    rw [TransformsTapesInExactTime, relatesInSteps_iff_step_iter_eq_some]
+    rw [← configs, Option.eq_some_iff_get_eq]
+    use configs_isSome_of_haltsAtStep (Nat.find_spec h_dom)
+    ext1
+    · simp
+      grind [haltsAtStep, Nat.find_spec h_dom]
+    · exact h_get
+  · intro ⟨t, h_iter⟩
+    rw [TransformsTapesInExactTime, relatesInSteps_iff_step_iter_eq_some] at h_iter
+    rw [← configs] at h_iter
+    have h_halts_at_t : tm.haltsAtStep tapes t := by simp [haltsAtStep, h_iter]
+    let h_halts : ∃ t, tm.haltsAtStep tapes t := ⟨t, h_halts_at_t⟩
+    use h_halts
+    have h_eq : Nat.find h_halts = t := halting_step_unique (Nat.find_spec h_halts) h_halts_at_t
+    simp [h_eq, h_iter]
+
+end MultiTapeTM
+
+end Turing
diff --git a/Cslib/Computability/Machines/MultiTapeTuring/CopyRoutine.lean b/Cslib/Computability/Machines/MultiTapeTuring/CopyRoutine.lean
new file mode 100644
index 00000000..566a3aef
--- /dev/null
+++ b/Cslib/Computability/Machines/MultiTapeTuring/CopyRoutine.lean
@@ -0,0 +1,53 @@
+/-
+Copyright (c) 2026 Christian Reitwiessner. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Christian Reitwiessner
+-/
+
+module
+
+public import Cslib.Computability.Machines.MultiTapeTuring.Basic
+public import Cslib.Computability.Machines.MultiTapeTuring.ListEncoding
+public import Cslib.Computability.Machines.MultiTapeTuring.WithTapes
+
+namespace Turing
+
+namespace Routines
+
+variable [Inhabited α] [Fintype α]
+
+def copy₁ : MultiTapeTM 2 (WithSep α) where
+  State := PUnit
+  stateFintype := inferInstance
+  q₀ := PUnit.unit
+  tr _ syms := sorry
+
+@[simp]
+lemma copy₁_eval_list {tapes : Fin 2 → List (List α)} :
+  copy₁.eval_list tapes =
+    Part.some (Function.update tapes 1 (((tapes 0).headD []) :: tapes 1)) := by
+  sorry
+
+/--
+A Turing machine that copies the first word on tape `i` to tape `j`.
+If Tape `i` is empty, pushes the empty word to tape `j`.
+-/
+public def copy {k : ℕ} (i j : Fin k)
+  (h_inj : [i, j].get.Injective := by intro x y; grind) :
+    MultiTapeTM k (WithSep α) :=
+  copy₁.with_tapes [i, j].get h_inj
+
+@[simp, grind =]
+public lemma copy_eval_list
+  {k : ℕ}
+  {i j : Fin k}
+  (h_inj : [i, j].get.Injective)
+  {tapes : Fin k → List (List α)} :
+  (copy i j h_inj).eval_list tapes = Part.some
+    (Function.update tapes j (((tapes i).headD []) :: (tapes j))) := by
+  simp_all [copy]
+  grind
+
+end Routines
+
+end Turing
diff --git a/Cslib/Computability/Machines/MultiTapeTuring/EqualRoutine.lean b/Cslib/Computability/Machines/MultiTapeTuring/EqualRoutine.lean
new file mode 100644
index 00000000..b5c3d88b
--- /dev/null
+++ b/Cslib/Computability/Machines/MultiTapeTuring/EqualRoutine.lean
@@ -0,0 +1,65 @@
+/-
+Copyright (c) 2026 Christian Reitwiessner. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Christian Reitwiessner
+-/
+
+module
+
+public import Cslib.Computability.Machines.MultiTapeTuring.Basic
+public import Cslib.Computability.Machines.MultiTapeTuring.ListEncoding
+public import Cslib.Computability.Machines.MultiTapeTuring.WithTapes
+
+import Mathlib.Logic.Function.Basic
+
+namespace Turing
+
+namespace Routines
+
+/--
+A 3-tape Turing machine that pushes the new word "1"
+to the third tape if the first words on the first and second tape are the same
+and otherwise pushes the empty word to the third tape.
+If one of the first two tapes is empty, uses the empty word for comparison.
+-/
+def eq₀ : MultiTapeTM 3 (WithSep OneTwo) where
+  State := PUnit
+  stateFintype := inferInstance
+  q₀ := PUnit.unit
+  tr _ syms := sorry
+
+@[simp]
+theorem eq₀_eval_list {tapes : Fin 3 → List (List OneTwo)} :
+  eq₀.eval_list tapes =
+    Part.some (Function.update tapes 2 ((if (tapes 0).headD [] = (tapes 1).headD [] then
+      [.one]
+    else
+      []) :: (tapes 2))) := by
+  sorry
+
+/--
+A Turing machine that pushes the new word "1"
+to tape `t` if the first words on tape `q` and tape `s` are the same
+and otherwise pushes the empty word to tape `t`.
+If one of the tapes `q` or `s` are empty, uses the empty word for comparison.
+-/
+public def eq {k : ℕ} (q s t : Fin k)
+  (h_inj : [q, s, t].get.Injective := by intro x y; grind) :
+  MultiTapeTM k (WithSep OneTwo) :=
+  eq₀.with_tapes [q, s, t].get h_inj
+
+@[grind =]
+public theorem eq_eval_list {k : ℕ} {q s t : Fin k}
+  (h_inj : [q, s, t].get.Injective)
+  {tapes : Fin k → List (List OneTwo)} :
+  (eq q s t).eval_list tapes =
+    Part.some (Function.update tapes t ((if (tapes q).headD [] = (tapes s).headD [] then
+      [.one]
+    else
+      []) :: (tapes t))) := by
+  simp [eq]
+  grind
+
+end Routines
+
+end Turing
diff --git a/Cslib/Computability/Machines/MultiTapeTuring/GraphReachability.lean b/Cslib/Computability/Machines/MultiTapeTuring/GraphReachability.lean
new file mode 100644
index 00000000..508491ca
--- /dev/null
+++ b/Cslib/Computability/Machines/MultiTapeTuring/GraphReachability.lean
@@ -0,0 +1,74 @@
+/-
+Copyright (c) 2026 Christian Reitwiessner. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Christian Reitwiessner
+-/
+
+module
+
+import Cslib.Foundations.Data.BiTape
+import Cslib.Foundations.Data.RelatesInSteps
+
+public import Cslib.Computability.Machines.MultiTapeTuring.Basic
+public import Cslib.Computability.Machines.MultiTapeTuring.ListEncoding
+public import Cslib.Computability.Machines.MultiTapeTuring.HeadStats
+
+-- TODO create a "common file"
+import Cslib.Computability.Machines.SingleTapeTuring.Basic
+
+namespace Turing
+
+
+-- This is an attempt at proving Savitch's theorem. We start by stating a generic
+-- space-efficient graph reachability algorithm.
+
+/-
+
+General idea, assume distance is power of two:
+
+def reachable(a, b, t, r):
+  if t = 0:
+    return r(a, b)
+  else:
+    for c in vertices:
+      if reachable(a, c, t - 1, r) and reachable(c, b, t - 1, r):
+        return True
+    return False
+
+Until we have a generic mechanism for recursion, let's translate this into a program that
+uses "goto", and every variable is a stack:
+
+def reachable(a, b, t, r):
+  terminate = 0
+  result = 0
+  section = [0]
+  while !terminate:
+    match section.pop()
+    | 0 =>
+      if t = 0:
+        result = r(a, b)
+        terminate = 1
+        section.push(7)
+      else:
+        section.push(1)
+    | 1 =>
+      c.push(0)
+      section.push(2)
+    | 2 =>
+      if c.top() = num_vertices:
+        section.push(6)
+      else:
+        a.push(a.top())
+        b.push(c.top())
+        section.push(0)
+        t.push(t.top() - 1)
+        section.push(3)
+    | 3 =>
+      section.push(4)
+
+
+
+-/
+
+
+end Turing
diff --git a/Cslib/Computability/Machines/MultiTapeTuring/HeadStats.lean b/Cslib/Computability/Machines/MultiTapeTuring/HeadStats.lean
new file mode 100644
index 00000000..3dbba8fb
--- /dev/null
+++ b/Cslib/Computability/Machines/MultiTapeTuring/HeadStats.lean
@@ -0,0 +1,67 @@
+/-
+Copyright (c) 2026 Christian Reitwiessner. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Christian Reitwiessner
+-/
+
+module
+
+public import Cslib.Computability.Machines.MultiTapeTuring.Basic
+public import Cslib.Computability.Machines.MultiTapeTuring.TapeExtension
+
+
+namespace Turing
+
+variable [Inhabited Symbol] [Fintype Symbol]
+
+variable {k : ℕ}
+
+/-- Statistics on the tape head movements. -/
+public structure HeadStats where
+  /-- The minimal (left-most) position of the head during the computation,
+  relative to the starting position. -/
+  min : ℤ
+  /-- The maximal (right-most) position of the head during the computation,
+  relative to the starting position. -/
+  max : ℤ
+  /-- The final position of the head after the computation, relative to the
+  starting position. -/
+  final : ℤ
+  h_bounds : min ≤ final ∧ final ≤ max ∧ min ≤ 0 ∧ 0 ≤ max
+
+/-- The space required. -/
+public def HeadStats.space (hs : HeadStats) : ℕ :=
+  (1 + hs.max - hs.min).toNat -- TODO we know it is nonnegative, is there a way to make use of that?
+
+
+/-- Compute the head statistics for a turing machine starting with a certain tape configuration. -/
+public def headStats (tm : MultiTapeTM k Symbol) (tapes : Fin k → BiTape Symbol) :
+  Part (Fin k → HeadStats) := sorry
+
+/-- Execute a Turing machine and also compute head statistics. -/
+public def MultiTapeTM.evalWithStats (tm : MultiTapeTM k Symbol) (tapes : Fin k → BiTape Symbol) :
+  Part ((Fin k → BiTape Symbol) × (Fin k → HeadStats)) := sorry
+
+-- move this somewhere else
+def seq (tm₁ tm₂ : MultiTapeTM k Symbol) : MultiTapeTM k Symbol := sorry
+
+def seq_combine_stats (stats₁ stats₂ : Fin k → HeadStats) : Fin k → HeadStats :=
+  fun i => match (stats₁ i, stats₂ i) with
+  | (⟨min₁, max₁, final₁, h₁⟩, ⟨min₂, max₂, final₂, h₂⟩) =>
+    ⟨min min₁ (min₂ + final₁),
+    max max₁ (max₂ + final₁),
+    final₁ + final₂,
+    by omega⟩
+
+lemma seq_evalWithStats (tm₁ tm₂ : MultiTapeTM k Symbol) (tapes : Fin k → BiTape Symbol) (i : Fin k) :
+  (seq tm₁ tm₂).evalWithStats tapes = do
+      let (tapes', stats₁) ← tm₁.evalWithStats tapes
+      let (tapes'', stats₂) ← tm₂.evalWithStats tapes'
+      return (tapes'', seq_combine_stats stats₁ stats₂) := by sorry
+
+-- Next step: relate space requirements and head stats.
+
+theorem stats_and_space (tm : MultiTapeTM k Symbol) (tapes tapes' : Fin k → BiTape Symbol) (s : ℕ) :
+  (∃ t, tm.TransformsTapesInTimeAndSpace tapes tapes' t s) ↔
+    ∃ hs, (∑ i, (hs i).space) ≤ s ∧ tm.evalWithStats tapes = .some (tapes', hs) := by sorry
+end Turing
diff --git a/Cslib/Computability/Machines/MultiTapeTuring/IsZeroRoutine.lean b/Cslib/Computability/Machines/MultiTapeTuring/IsZeroRoutine.lean
new file mode 100644
index 00000000..265f9e42
--- /dev/null
+++ b/Cslib/Computability/Machines/MultiTapeTuring/IsZeroRoutine.lean
@@ -0,0 +1,35 @@
+/-
+Copyright (c) 2026 Christian Reitwiessner. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Christian Reitwiessner
+-/
+
+module
+
+public import Cslib.Computability.Machines.MultiTapeTuring.Basic
+public import Cslib.Computability.Machines.MultiTapeTuring.ListEncoding
+public import Cslib.Computability.Machines.MultiTapeTuring.PushRoutine
+public import Cslib.Computability.Machines.MultiTapeTuring.PopRoutine
+public import Cslib.Computability.Machines.MultiTapeTuring.SequentialCombinator
+public import Cslib.Computability.Machines.MultiTapeTuring.IteCombinator
+
+namespace Turing
+
+variable {k : ℕ}
+
+namespace Routines
+
+/--
+A Turing machine that computes the logical negation: It replaces an empty (or non-existing) head
+on tape `i` by the word "1" and everything else by the empty word. -/
+public def isZero (i : Fin k) := ite i (pop i ;ₜ push i []) (pop i ;ₜ push i [OneTwo.one])
+
+@[simp, grind =]
+public theorem isZero_eval_list {i : Fin k} {tapes : Fin k → List (List OneTwo)} :
+  (isZero i).eval_list tapes = .some (Function.update tapes i (
+    (if (tapes i).headD [] = [] then [OneTwo.one] else []) :: (tapes i).tail)) := by
+  simp [isZero]
+  grind
+
+end Routines
+end Turing
diff --git a/Cslib/Computability/Machines/MultiTapeTuring/IteCombinator.lean b/Cslib/Computability/Machines/MultiTapeTuring/IteCombinator.lean
new file mode 100644
index 00000000..9d9d41b9
--- /dev/null
+++ b/Cslib/Computability/Machines/MultiTapeTuring/IteCombinator.lean
@@ -0,0 +1,40 @@
+/-
+Copyright (c) 2026 Christian Reitwiessner. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Christian Reitwiessner
+-/
+
+module
+
+public import Cslib.Computability.Machines.MultiTapeTuring.Basic
+public import Cslib.Computability.Machines.MultiTapeTuring.ListEncoding
+
+namespace Turing
+
+namespace Routines
+
+variable [Inhabited α] [Fintype α]
+variable {k : ℕ}
+
+/--
+A Turing machine combinator that runs `tm₁` if the first word on tape `i` exists and is non-empty,
+otherwise it runs `tm₂`. -/
+public def ite (i : Fin k) (tm₁ tm₂ : MultiTapeTM k (WithSep α)) :
+    MultiTapeTM k (WithSep α) where
+  State := PUnit
+  stateFintype := inferInstance
+  q₀ := PUnit.unit
+  tr _ syms := sorry
+
+@[simp, grind =]
+public theorem ite_eval_list
+  {i : Fin k}
+  {tm₁ tm₂ : MultiTapeTM k (WithSep α)}
+  {tapes : Fin k → List (List α)} :
+  (ite i tm₁ tm₂).eval_list tapes =
+    if (tapes i).headD [] = [] then tm₂.eval_list tapes else tm₁.eval_list tapes := by
+  sorry
+
+end Routines
+
+end Turing
diff --git a/Cslib/Computability/Machines/MultiTapeTuring/ListEncoding.lean b/Cslib/Computability/Machines/MultiTapeTuring/ListEncoding.lean
new file mode 100644
index 00000000..f16a1984
--- /dev/null
+++ b/Cslib/Computability/Machines/MultiTapeTuring/ListEncoding.lean
@@ -0,0 +1,215 @@
+/-
+Copyright (c) 2026 Christian Reitwiessner. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Christian Reitwiessner
+-/
+
+module
+
+public import Cslib.Computability.Machines.MultiTapeTuring.Basic
+public import Cslib.Computability.Machines.MultiTapeTuring.WithTapes
+
+import Mathlib.Tactic.DeriveFintype
+
+namespace Turing
+
+variable {Symbol : Type}
+variable [Inhabited Symbol] [Fintype Symbol]
+
+/--
+An alphabet that contains exactly two symbols, 1 and 2.
+TODO use an embedding or something else that is more flexible
+-/
+public inductive OneTwo where
+  | one
+  | two
+deriving DecidableEq, Inhabited, Fintype
+
+
+/-- An alphabet for list encoding -/
+public inductive WithSep (Symbol : Type) where
+  | ofChar (c : Symbol)
+  | comma
+  -- TODO need to decide if we want to encode lists with parentheses or not.
+  -- Is annoying when pushing and popping from lists, but can be useful to avoid
+  -- running "off the tape"
+  -- | lParen
+  -- | rParen
+deriving Fintype, DecidableEq, Inhabited
+
+/-- A list of words is transformed by appending a comma after each word and concatenating.
+Note that the comma is not only a separator but also appears as the final character
+of the resulting string (if the list is non-empty). -/
+public def listToString (ls : List (List Symbol)) : List (WithSep Symbol) :=
+  (ls.map (fun w : List Symbol => (w.map .ofChar) ++ [.comma])).flatten
+
+public def stringToList (s : List (Option (WithSep Symbol))) : Option (List (List Symbol)) :=
+  s.foldr (fun c acc =>
+    match c with
+    | none => none
+    | some .comma => acc.map ([] :: ·)
+    | some (.ofChar c) => acc.bind (fun ws =>
+        match ws with
+        | [] => none
+        | w :: rest => some ((c :: w) :: rest)))
+  (some [])
+
+/-- Encodes a list of words into a tape. -/
+public def listToTape (ls : List (List Symbol)) : BiTape (WithSep Symbol) :=
+  BiTape.mk₁ (listToString ls)
+
+/-- Only tapes where the head is at the left end are accepted. -/
+public def tapeToList (tape : BiTape (WithSep Symbol)) : Option (List (List Symbol)) :=
+  match (tape.left.toList, tape.head) with
+  | ([], c) => stringToList (c :: tape.right.toList)
+  | _ => none
+
+public def tapesToLists (tapes : Fin k → BiTape (WithSep Symbol)) :
+    Option (Fin k → List (List Symbol)) :=
+  if h : ∀ i, (tapeToList (tapes i)).isSome then
+    some (fun i => (tapeToList (tapes i)).get (h i))
+  else
+    none
+
+/-- The Turing machine `tm` transforms the list-encoded tapes `tapes` into the list-encoded
+tapes `tapes'`. -/
+public def MultiTapeTM.TransformsLists
+    (tm : MultiTapeTM k (WithSep Symbol))
+    (tapes tapes' : Fin k → List (List Symbol)) : Prop :=
+  tm.TransformsTapes (listToTape ∘ tapes) (listToTape ∘ tapes')
+
+/-- The Turing machine `tm` halts starting with list-encoded tapes `tapes`
+and also always outputs properly list-encoded tapes. -/
+public def MultiTapeTM.HaltsOnLists
+    (tm : MultiTapeTM k (WithSep Symbol))
+    (tapes : Fin k → List (List Symbol)) : Prop :=
+  ∃ tapes', tm.TransformsLists tapes tapes'
+
+/-- Execute the Turing machine `tm` on the list-encoded tapes `tapes`. -/
+public def MultiTapeTM.eval_list
+    (tm : MultiTapeTM k (WithSep Symbol))
+    (tapes : Fin k → List (List Symbol)) :
+    Part (Fin k → List (List Symbol)) :=
+  (tm.eval (listToTape ∘ tapes)).bind fun tapes => tapesToLists tapes
+
+@[simp]
+public theorem MultiTapeTM.HaltsOnLists_of_eval_list
+    {tm : MultiTapeTM k (WithSep Symbol)}
+    {tapes : Fin k → List (List Symbol)}
+    (h_dom : (tm.eval_list tapes).Dom) :
+    tm.HaltsOnLists tapes := by
+  sorry
+
+@[simp]
+public lemma MultiTapeTM.eval_list_dom_of_halts_on_lists
+    {tm : MultiTapeTM k (WithSep Symbol)}
+    {tapes : Fin k → List (List Symbol)}
+    (h_halts : tm.HaltsOnLists tapes) :
+    (tm.eval_list tapes).Dom := by
+  sorry
+
+/-- Execute the Turing machine `tm` knowing that it always halts, thus yielding a total function
+on the tapes. -/
+@[simp]
+public def MultiTapeTM.eval_list_tot
+    (tm : MultiTapeTM k (WithSep Symbol))
+    (h_alwaysHalts : ∀ tapes, tm.HaltsOnLists tapes)
+    (tapes : Fin k → List (List Symbol)) :
+  Fin k → List (List Symbol) :=
+    (tm.eval_list tapes).get (tm.eval_list_dom_of_halts_on_lists (h_alwaysHalts tapes))
+
+@[simp, grind =]
+public theorem MultiTapeTM.extend_eval_list
+    {k₁ k₂ : ℕ} {h_le : k₁ ≤ k₂}
+    {tm : MultiTapeTM k₁ (WithSep Symbol)}
+    {tapes : Fin k₂ → List (List Symbol)} :
+  (tm.extend h_le).eval_list tapes =
+    (tm.eval_list (tapes ⟨·, by omega⟩)).map (fun tapes' =>
+      fun i : Fin k₂ => if h : i.val < k₁ then tapes' ⟨i, h⟩ else tapes i) := by
+  sorry
+
+@[simp, grind =]
+public theorem MultiTapeTM.permute_tapes_eval_list
+  (tm : MultiTapeTM k (WithSep Symbol)) (σ : Equiv.Perm (Fin k))
+  (tapes : Fin k → List (List Symbol)) :
+  (tm.permute_tapes σ).eval_list tapes =
+    (tm.eval_list (tapes ∘ σ)).map (fun tapes' => tapes' ∘ σ.symm) := by
+  sorry
+
+@[simp, grind =]
+public theorem MultiTapeTM.with_tapes_eval_list
+  {k₁ k₂ : ℕ}
+  {tm : MultiTapeTM k₁ (WithSep Symbol)} {f : Fin k₁ → Fin k₂} {h_inj : f.Injective}
+  {tapes : Fin k₂ → List (List Symbol)} :
+  (tm.with_tapes f h_inj).eval_list tapes =
+    (tm.eval_list (tapes ∘ f)).map
+      (fun tapes' => fun t => apply_updates tapes tapes' f t) := by
+   sorry
+
+-- def MultiTapeTM.TransformsListsWithStats
+--     (tm : MultiTapeTM k (WithSep Symbol))
+--     (tapes : Fin k → List (List Symbol))
+--     (ts : (Fin k → List (List Symbol)) × (Fin k → HeadStats)) : Prop :=
+--     tm.evalWithStats (listToTape ∘ tapes) = .some (listToTape ∘ ts.1, ts.2)
+
+-- /--
+-- Evaluate the Turing machine `tm` on the list-encoded tapes `tapes` and also return the head
+-- statistics of the computation.
+-- -/
+-- public noncomputable def MultiTapeTM.evalWithStats_list
+--     (tm : MultiTapeTM k (WithSep Symbol))
+--     (tapes : Fin k → List (List Symbol)) :
+--     Part ((Fin k → List (List Symbol)) × (Fin k → HeadStats)) :=
+--   ⟨∃ ts, tm.TransformsListsWithStats tapes ts, fun h => h.choose⟩
+
+-- TODO for machines running on lists, we can actually have more precise head stats:
+-- we know (and should enforce) that the head never moves to the right of the rightmost symbol
+-- and always starts and ends on the leftmost symbol (and if the tape is empty, it never moves
+-- to the right of the starting position).
+-- So the minimal information we need is (per tape) the amount of symbols that were used beyond
+-- the max of the ones in the initial and final configuration.
+-- TODO it also makes sense to allow upper bounds on that.
+
+/--
+The Turing machine `tm` computes a total function on lists and this uniquely
+determined function is `f`.
+-/
+public def MultiTapeTM.computes
+    (tm : MultiTapeTM k (WithSep Symbol))
+    (f : (Fin k → List (List Symbol)) → (Fin k → List (List Symbol))) : Prop :=
+  ∀ tapes, tm.eval_list tapes = .some (f tapes)
+
+public theorem MultiTapeTM.eval_of_computes
+    {Symbol : Type}
+    [Fintype Symbol]
+    {tm : MultiTapeTM k (WithSep Symbol)}
+    {f : (Fin k → List (List Symbol)) → (Fin k → List (List Symbol))}
+    (h_computes : tm.computes f)
+    {tapes : Fin k → List (List Symbol)} :
+    tm.eval_list tapes = .some (f tapes) := by
+  specialize h_computes tapes
+  simpa [MultiTapeTM.computes] using h_computes
+
+
+/-- Dyadic encoding of natural numbers. -/
+public def dya (n : ℕ) : List OneTwo :=
+  if n = 0 then []
+  else if Even n then
+    dya (n / 2 - 1) ++ [.two]
+  else
+    dya ((n - 1) / 2) ++ [.one]
+
+/-- Dyadic decoding of natural numbers. -/
+public def dya_inv : List OneTwo → ℕ := sorry
+
+@[simp, grind =]
+public lemma dya_inv_zero : dya_inv [] = 0 := by
+  sorry
+
+@[simp, grind =]
+public lemma dya_inv_dya (n : ℕ) : dya_inv (dya n) = n := by sorry
+
+@[simp, grind =]
+public lemma dya_dya_inv (w : List OneTwo) : dya (dya_inv w) = w := by sorry
+
+end Turing
diff --git a/Cslib/Computability/Machines/MultiTapeTuring/LoopCombinator.lean b/Cslib/Computability/Machines/MultiTapeTuring/LoopCombinator.lean
new file mode 100644
index 00000000..217f62a3
--- /dev/null
+++ b/Cslib/Computability/Machines/MultiTapeTuring/LoopCombinator.lean
@@ -0,0 +1,59 @@
+/-
+Copyright (c) 2026 Christian Reitwiessner. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Christian Reitwiessner
+-/
+
+module
+
+public import Cslib.Computability.Machines.MultiTapeTuring.Basic
+public import Cslib.Computability.Machines.MultiTapeTuring.ListEncoding
+
+namespace Turing
+
+namespace Routines
+
+variable {k : ℕ}
+
+/--
+A Turing machine that executes `tm` a number of times as given by the first word on tape `i`.
+If tape `i` is empty, do not execute the TM.
+Note that the iteration counter is not directly available to `tm`. -/
+public def loop (i : Fin k)
+  (tm : MultiTapeTM k (WithSep OneTwo)) : MultiTapeTM (k + 3) (WithSep OneTwo) :=
+  sorry
+  -- let target : Fin (k + (aux + 3)) := ⟨aux, by omega⟩
+  -- let counter : Fin (k + (aux + 3)) := ⟨aux + 1, by omega⟩
+  -- let cond : Fin (k + (aux + 3)) := ⟨aux + 2, by omega⟩
+  -- (copy (k := k + (aux + 3)) i target (h_neq := by grind) <;>
+  -- push counter [] <;>
+  -- neq target counter cond (by grind) (by grind) (by grind) <;>
+  -- doWhile cond (
+  --   pop cond <;>
+  --   tm.toMultiTapeTM <;>
+  --   succ counter <;>
+  --   neq target counter cond (by grind) (by grind) (by grind)) <;>
+  -- pop cond <;>
+  -- pop counter <;>
+  -- pop target).set_aux_tapes (aux + 3)
+
+
+@[simp]
+public theorem loop_eval_list {i : Fin k}
+  {tm : MultiTapeTM k (WithSep OneTwo)}
+  {tapes : Fin (k + 3) → List (List OneTwo)} :
+  (loop i tm).eval_list tapes =
+      (((Part.bind · tm.eval_list)^[dya_inv ((tapes ⟨i, by omega⟩).headD [])]
+        (Part.some (tapes_take tapes k (by omega))))).map
+          fun tapes' => tapes_extend_by tapes' tapes := by
+  sorry
+
+@[simp]
+public theorem loop_halts_of_halts {i : Fin k}
+  {tm : MultiTapeTM k (WithSep OneTwo)}
+  (h_halts : ∀ tapes, tm.HaltsOnLists tapes) :
+  ∀ tapes, (loop i tm).HaltsOnLists tapes := by
+  sorry
+
+end Routines
+end Turing
diff --git a/Cslib/Computability/Machines/MultiTapeTuring/MoveRoutine.lean b/Cslib/Computability/Machines/MultiTapeTuring/MoveRoutine.lean
new file mode 100644
index 00000000..443eb44a
--- /dev/null
+++ b/Cslib/Computability/Machines/MultiTapeTuring/MoveRoutine.lean
@@ -0,0 +1,102 @@
+/-
+Copyright (c) 2026 Christian Reitwiessner. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Christian Reitwiessner
+-/
+
+module
+
+public import Cslib.Computability.Machines.MultiTapeTuring.Basic
+
+namespace Turing
+
+namespace Routines
+
+variable [Inhabited α] [Fintype α]
+
+-- /-- A 1-tape Turing machine that moves its head in a given direction
+-- once and then halts. -/
+-- public def move (dir : Dir) : MultiTapeTM 1 α where
+--   Λ := PUnit
+--   q₀ := 0
+--   M _ syms := (fun i => ⟨syms i, some dir⟩, none)
+
+-- @[simp]
+-- public lemma move_eval (tape : BiTape α) (dir : Turing.Dir) :
+--   (move dir).eval (fun _ => tape) = .some (fun _ => tape.move dir) := by
+--   sorry
+--   rw [MultiTapeTM.eval_iff_exists_steps_iter_eq_some]
+--   use 1
+--   rfl
+
+-- /-- A 1-tape Turing machine that moves its head in a given direction until a condition
+-- on the read symbol is met. -/
+-- public def move_until (dir : Turing.Dir) (cond : (Option α) → Bool) : MultiTapeTM 1 α where
+--   Λ := PUnit
+--   q₀ := PUnit.unit
+--   M q syms := match cond (syms 0) with
+--     | false  => (fun _ => ⟨syms 0, some dir⟩, some q)
+--     | true => (fun _ => ⟨syms 0, none⟩, none)
+
+-- lemma move_until_step_cond_false
+--   {tape : BiTape α}
+--   {stop_condition : Option α → Bool}
+--   (h_neg_stop : ¬ stop_condition tape.head) :
+--   (move_until .right stop_condition).step
+--     ⟨some (move_until .right stop_condition).q₀, (fun _ => tape)⟩ =
+--     some ⟨some (move_until .right stop_condition).q₀, (fun _ => tape.move .right)⟩ := by
+--   simp [move_until, h_neg_stop, BiTape.optionMove, MultiTapeTM.step]
+
+-- lemma move_until_step_cond_true
+--   {tape : BiTape α}
+--   {stop_condition : Option α → Bool}
+--   (h_neg_stop : stop_condition tape.head) :
+--   (move_until .right stop_condition).step
+--     ⟨some (move_until .right stop_condition).q₀, (fun _ => tape)⟩ =
+--     some ⟨none, (fun _ => tape)⟩ := by
+--   simp [move_until, h_neg_stop, BiTape.optionMove, MultiTapeTM.step]
+
+-- public theorem move_until.right_semantics
+--   (tape : BiTape α)
+--   (stop_condition : Option α → Bool)
+--   (h_stop : ∃ n : ℕ, stop_condition (tape.nth n)) :
+--   (move_until .right stop_condition).eval (fun _ => tape) =
+--     .some (fun _ => tape.move_int (Nat.find h_stop)) := by
+--   rw [MultiTapeTM.eval_iff_exists_steps_iter_eq_some]
+--   let n := Nat.find h_stop
+--   use n.succ
+--   have h_not_stop_of_lt : ∀ k < n, ¬ stop_condition (tape.move_int k).head := by
+--     intro k hk
+--     simp [Nat.find_min h_stop hk]
+--   have h_iter : ∀ k < n, (Option.bind · (move_until .right stop_condition).step)^[k]
+--       (some ⟨some (move_until .right stop_condition).q₀, fun _ => tape⟩) =
+--       some ⟨some (move_until .right stop_condition).q₀, fun _ => tape.move_int k⟩ := by
+--     intro k hk
+--     induction k with
+--     | zero =>
+--       simp [BiTape.move_int]
+--     | succ k ih =>
+--       have hk' : k < n := Nat.lt_of_succ_lt hk
+--       rw [Function.iterate_succ_apply', ih hk']
+--       simp only [Option.bind_some, move_until_step_cond_false (h_not_stop_of_lt k hk')]
+--       simp [BiTape.move, ← BiTape.move_int_one_eq_move_right, BiTape.move_int_move_int]
+--   have h_n_eq : n = Nat.find h_stop := by grind
+--   by_cases h_n_zero : n = 0
+--   · have h_stop_cond : stop_condition (tape.head) := by simp_all [n]
+--     let h_step := move_until_step_cond_true h_stop_cond
+--     simp [h_step, ← h_n_eq, h_n_zero]
+--   · obtain ⟨n', h_n'_eq_n_succ⟩ := Nat.exists_eq_add_one_of_ne_zero h_n_zero
+--     rw [h_n'_eq_n_succ, Function.iterate_succ_apply', Function.iterate_succ_apply']
+--     have h_n'_lt_n : n' < n := by omega
+--     simp only [MultiTapeTM.initCfgTapes, MultiTapeTM.haltCfgTapes]
+--     rw [h_iter n' h_n'_lt_n]
+--     simp only [Option.bind_some, move_until_step_cond_false (h_not_stop_of_lt n' h_n'_lt_n)]
+--     simp only [BiTape.move, ← BiTape.move_int_one_eq_move_right, BiTape.move_int_move_int]
+--     rw [show (n' : ℤ) + 1 = n by omega]
+--     have h_n_stop : stop_condition ((tape.move_int n).head) := by
+--       simpa [n] using Nat.find_spec h_stop
+--     simpa using move_until_step_cond_true h_n_stop
+
+end Routines
+
+end Turing
diff --git a/Cslib/Computability/Machines/MultiTapeTuring/MulRoutine.lean b/Cslib/Computability/Machines/MultiTapeTuring/MulRoutine.lean
new file mode 100644
index 00000000..a5e4f72d
--- /dev/null
+++ b/Cslib/Computability/Machines/MultiTapeTuring/MulRoutine.lean
@@ -0,0 +1,78 @@
+/-
+Copyright (c) 2026 Christian Reitwiessner. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Christian Reitwiessner
+-/
+
+module
+
+public import Cslib.Computability.Machines.MultiTapeTuring.Basic
+public import Cslib.Computability.Machines.MultiTapeTuring.ListEncoding
+public import Cslib.Computability.Machines.MultiTapeTuring.AddRoutine
+public import Cslib.Computability.Machines.MultiTapeTuring.PushRoutine
+public import Cslib.Computability.Machines.MultiTapeTuring.LoopCombinator
+public import Cslib.Computability.Machines.MultiTapeTuring.SequentialCombinator
+public import Cslib.Computability.Machines.MultiTapeTuring.WithTapes
+
+namespace Turing
+
+variable {k : ℕ}
+
+namespace Routines
+
+-- Multiplies the heads of 0 and 1 and stores the result in 2.
+def mul₀ : MultiTapeTM 9 (WithSep OneTwo) :=
+  push 2 [] ;ₜ loop 0 (add_assign 1 2 3)
+
+@[simp]
+lemma add_assign_iter {i j aux : Fin (k + 6)} {r : ℕ}
+  (h_inj : [i, j, aux, aux + 1, aux + 2, aux + 3].get.Injective)
+  {tapes : Fin (k + 6) → List (List OneTwo)} :
+  (Part.bind · (add_assign i j aux h_inj).eval_list)^[r] (.some tapes) =
+    Part.some (Function.update tapes j (
+      if r = 0 then
+        tapes j
+      else
+        (dya ((dya_inv ((tapes j).headD [])) + r * (dya_inv ((tapes i).headD [])))) ::
+          (tapes j).tail)) := by
+  induction r with
+  | zero => simp
+  | succ r ih =>
+    rw [Function.iterate_succ_apply']
+    simp only [ih, Part.bind_some]
+    rw [add_assign_eval_list]
+    have h_neq : i ≠ j := Function.Injective.ne h_inj (a₁ := 0) (a₂ := 1) (by grind)
+    simp
+    grind
+
+@[simp]
+theorem mul₀_eval_list {tapes : Fin 9 → List (List OneTwo)} :
+  mul₀.eval_list tapes = .some
+    (Function.update tapes 2 (
+      (dya (dya_inv ((tapes 0).headD []) * dya_inv ((tapes 1).headD [])) :: (tapes 2)))) := by
+  by_cases h_zero: dya_inv ((tapes 0).head?.getD []) = 0
+  <;> simp [mul₀, h_zero] <;> grind
+
+/--
+A Turing machine that multiplies the heads of tapes i and j and pushes the result to tape l.
+If tapes are empty, their heads are assumed to be zero.
+-/
+public def mul
+  (i j l : Fin (k + 9))
+  (aux : Fin (k + 9) := ⟨k + 3, by omega⟩)
+  (h_inj : [i, j, l, aux, aux + 1, aux + 2, aux + 3, aux + 4, aux + 5].get.Injective := by decide) :
+  MultiTapeTM (k + 9) (WithSep OneTwo) :=
+  mul₀.with_tapes [i, j, l, aux, aux + 1, aux + 2, aux + 3, aux + 4, aux + 5].get h_inj
+
+@[simp, grind =]
+public theorem mul_eval_list (i j l aux : Fin (k + 9))
+  {h_inj : [i, j, l, aux, aux + 1, aux + 2, aux + 3, aux + 4, aux + 5].get.Injective}
+  {tapes : Fin (k + 9) → List (List OneTwo)} :
+  (mul i j l aux h_inj).eval_list tapes =
+      .some (Function.update tapes l (
+        (dya (dya_inv ((tapes i).headD []) * dya_inv ((tapes j).headD [])) :: (tapes l)))) := by
+  simpa [mul] using apply_updates_function_update h_inj
+
+end Routines
+
+end Turing
diff --git a/Cslib/Computability/Machines/MultiTapeTuring/PopRoutine.lean b/Cslib/Computability/Machines/MultiTapeTuring/PopRoutine.lean
new file mode 100644
index 00000000..a18b72b9
--- /dev/null
+++ b/Cslib/Computability/Machines/MultiTapeTuring/PopRoutine.lean
@@ -0,0 +1,45 @@
+/-
+Copyright (c) 2026 Christian Reitwiessner. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Christian Reitwiessner
+-/
+
+module
+
+public import Cslib.Computability.Machines.MultiTapeTuring.Basic
+public import Cslib.Computability.Machines.MultiTapeTuring.ListEncoding
+public import Cslib.Computability.Machines.MultiTapeTuring.WithTapes
+
+namespace Turing
+
+namespace Routines
+
+variable [Inhabited α] [Fintype α]
+
+def pop₁ : MultiTapeTM 1 (WithSep α) where
+  State := PUnit
+  stateFintype := inferInstance
+  q₀ := PUnit.unit
+  tr _ syms := sorry
+
+@[simp]
+lemma pop₁_eval_list {tapes : Fin 1 → List (List α)} :
+  pop₁.eval_list tapes = .some (Function.update tapes 0 (tapes 0).tail) := by
+  sorry
+
+/--
+A Turing machine that removes the first word on tape `i`. If the tape is empty, does nothing.
+-/
+public def pop {k : ℕ} (i : Fin k) : MultiTapeTM k (WithSep α) :=
+  pop₁.with_tapes [i].get (by intro x y; grind)
+
+@[simp, grind =]
+public theorem pop_eval_list {k : ℕ} {i : Fin k}
+  {tapes : Fin k → List (List α)} :
+  (pop i).eval_list tapes = .some (Function.update tapes i (tapes i).tail) := by
+  have h_inj : [i].get.Injective := by intro x y; grind
+  simp_all [pop]
+
+end Routines
+
+end Turing
diff --git a/Cslib/Computability/Machines/MultiTapeTuring/PushRoutine.lean b/Cslib/Computability/Machines/MultiTapeTuring/PushRoutine.lean
new file mode 100644
index 00000000..d6c0bac7
--- /dev/null
+++ b/Cslib/Computability/Machines/MultiTapeTuring/PushRoutine.lean
@@ -0,0 +1,46 @@
+/-
+Copyright (c) 2026 Christian Reitwiessner. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Christian Reitwiessner
+-/
+
+module
+
+public import Cslib.Computability.Machines.MultiTapeTuring.Basic
+public import Cslib.Computability.Machines.MultiTapeTuring.ListEncoding
+public import Cslib.Computability.Machines.MultiTapeTuring.WithTapes
+
+namespace Turing
+
+namespace Routines
+
+variable [Inhabited α] [Fintype α]
+
+def push₁ (w : List α) : MultiTapeTM 1 (WithSep α) where
+  State := PUnit
+  stateFintype := inferInstance
+  q₀ := PUnit.unit
+  tr _ syms := sorry
+
+@[simp]
+lemma push₁_eval_list {w : List α} {tapes : Fin 1 → List (List α)} :
+  (push₁ w).eval_list tapes = .some (Function.update tapes 0 (w :: (tapes 0))) := by
+  sorry
+
+/--
+A Turing machine that pushes the word `w` to tape `i`.
+-/
+public def push {k : ℕ} (i : Fin k) (w : List α) : MultiTapeTM k (WithSep α) :=
+  (push₁ w).with_tapes [i].get (by intro x y; grind)
+
+@[simp, grind =]
+public theorem push_eval_list {k : ℕ}
+  {i : Fin k} {w : List α} {tapes : Fin k → List (List α)} :
+  (push i w).eval_list tapes =
+    .some (Function.update tapes i (w :: (tapes i))) := by
+  have h_inj : [i].get.Injective := by intro x y; grind
+  simp_all [push]
+
+end Routines
+
+end Turing
diff --git a/Cslib/Computability/Machines/MultiTapeTuring/SequentialCombinator.lean b/Cslib/Computability/Machines/MultiTapeTuring/SequentialCombinator.lean
new file mode 100644
index 00000000..e39726df
--- /dev/null
+++ b/Cslib/Computability/Machines/MultiTapeTuring/SequentialCombinator.lean
@@ -0,0 +1,57 @@
+/-
+Copyright (c) 2026 Christian Reitwiessner. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Christian Reitwiessner
+-/
+
+module
+
+import Cslib.Foundations.Data.BiTape
+import Cslib.Foundations.Data.RelatesInSteps
+
+public import Cslib.Computability.Machines.MultiTapeTuring.Basic
+public import Cslib.Computability.Machines.MultiTapeTuring.ListEncoding
+
+namespace Turing
+
+namespace MultiTapeTM
+
+variable [Inhabited Symbol] [Fintype Symbol]
+variable {k : ℕ}
+
+/--
+Sequential combination of Turing machines. Runs `tm₁` and then `tm₂` on the resulting tapes
+(if the first one halts).
+-/
+public def seq (tm₁ tm₂ : MultiTapeTM k Symbol) : MultiTapeTM k Symbol := sorry
+
+public theorem seq_eval
+  (tm₁ tm₂ : MultiTapeTM k Symbol)
+  (tapes₀ : Fin k → BiTape Symbol) :
+  (seq tm₁ tm₂).eval tapes₀ =
+    tm₁.eval tapes₀ >>= fun tape₁ => tm₂.eval tape₁ := by
+  sorry
+
+@[simp, grind =]
+public theorem seq_eval_list
+  {tm₁ tm₂ : MultiTapeTM k (WithSep Symbol)}
+  {tapes₀ : Fin k → List (List Symbol)} :
+  (seq tm₁ tm₂).eval_list tapes₀ =
+    tm₁.eval_list tapes₀ >>= fun tape₁ => tm₂.eval_list tape₁ := by
+  sorry
+
+public theorem seq_associative
+  (tm₁ tm₂ tm₃ : MultiTapeTM k Symbol)
+  (tapes₀ : Fin k → List (List Symbol)) :
+  (seq (seq tm₁ tm₂) tm₃).eval = (seq tm₁ (seq tm₂ tm₃)).eval := by
+  sorry
+
+/--
+Sequential combination of Turing machines.
+-/
+infixl:90 " ;ₜ " => seq
+
+
+end MultiTapeTM
+
+end Turing
diff --git a/Cslib/Computability/Machines/MultiTapeTuring/SuccRoutine.lean b/Cslib/Computability/Machines/MultiTapeTuring/SuccRoutine.lean
new file mode 100644
index 00000000..0ed3fc41
--- /dev/null
+++ b/Cslib/Computability/Machines/MultiTapeTuring/SuccRoutine.lean
@@ -0,0 +1,60 @@
+/-
+Copyright (c) 2026 Christian Reitwiessner. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Christian Reitwiessner
+-/
+
+module
+
+import Cslib.Foundations.Data.BiTape
+import Cslib.Foundations.Data.RelatesInSteps
+
+public import Cslib.Computability.Machines.MultiTapeTuring.Basic
+public import Cslib.Computability.Machines.MultiTapeTuring.ListEncoding
+public import Cslib.Computability.Machines.MultiTapeTuring.WithTapes
+
+import Mathlib.Data.Nat.Bits
+
+namespace Turing
+
+namespace Routines
+
+def succ₀ : MultiTapeTM 1 (WithSep OneTwo) where
+  State := PUnit
+  stateFintype := inferInstance
+  q₀ := PUnit.unit
+  tr _ syms := sorry
+
+@[simp]
+lemma succ₀_eval_list {tapes : Fin 1 → List (List OneTwo)} :
+  succ₀.eval_list tapes = .some (Function.update tapes 0
+    ((dya (dya_inv ((tapes 0).headD [])).succ) :: (tapes 0).tail)) := by
+  sorry
+
+/--
+A Turing machine that increments the head of tape `i` by one (in dyadic encoding).
+Pushes zero if the tape is empty. -/
+public def succ {k : ℕ} (i : Fin k) : MultiTapeTM k (WithSep OneTwo) :=
+  succ₀.with_tapes [i].get (by intro x y; grind)
+
+@[simp]
+public theorem succ_eval_list {k : ℕ} {i : Fin k} {tapes : Fin k → List (List OneTwo)} :
+  (succ i).eval_list tapes = .some (Function.update tapes i
+    ((dya (dya_inv ((tapes i).headD [])).succ) :: (tapes i).tail)) := by
+  simpa [succ] using apply_updates_function_update (by intro x y; grind)
+
+-- lemma succ₀_evalWithStats_list {n : ℕ} {ls : List (List OneTwo)} :
+--   succ₀.evalWithStats_list [(dya n) :: ls].get =
+--     .some (
+--       [(dya n.succ) :: ls].get,
+--       -- this depends on if we have overflow on the highest dyadic character or not.
+--       if (dya n.succ).length = (dya n).length then
+--         [⟨0, (dya n).length, 0, by omega⟩].get
+--       else
+--         [⟨-1, (dya n).length, -1, by omega⟩].get) := by
+--   sorry
+
+
+end Routines
+
+end Turing
diff --git a/Cslib/Computability/Machines/MultiTapeTuring/TapeExtension.lean b/Cslib/Computability/Machines/MultiTapeTuring/TapeExtension.lean
new file mode 100644
index 00000000..7e414cdb
--- /dev/null
+++ b/Cslib/Computability/Machines/MultiTapeTuring/TapeExtension.lean
@@ -0,0 +1,72 @@
+/-
+Copyright (c) 2026 Christian Reitwiessner. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Christian Reitwiessner
+-/
+
+module
+
+public import Cslib.Computability.Machines.MultiTapeTuring.Basic
+
+namespace Turing
+
+variable [Inhabited Symbol] [Fintype Symbol]
+
+/-- Extend a Turing machine to work with more tapes.
+The added tapes are not acted upon. -/
+public def MultiTapeTM.extend {k₁ k₂ : ℕ} (h_le : k₁ ≤ k₂)
+    (tm : MultiTapeTM k₁ Symbol) : MultiTapeTM k₂ Symbol where
+  State := tm.State
+  stateFintype := tm.stateFintype
+  q₀ := tm.q₀
+  tr := fun q syms => match tm.tr q (fun i => syms ⟨i, by omega⟩) with
+    | (stmts, q') =>
+      (fun i => if h : i < k₁ then stmts ⟨i, h⟩ else default, q')
+
+/--
+Restrict a sequence of tapes to the first `k'` tapes.
+-/
+@[simp]
+public abbrev tapes_take
+  {γ : Type}
+  (tapes : Fin k → γ)
+  (k' : ℕ)
+  (h_le : k' ≤ k)
+  (i : Fin k') : γ :=
+  tapes ⟨i, by omega⟩
+
+@[simp]
+public lemma Function.update_tapes_take
+  {γ : Type}
+  (k : ℕ)
+  {k' : ℕ}
+  {h_le : k' ≤ k}
+  {tapes : Fin k → γ}
+  {p : Fin k'}
+  {v : γ} :
+  Function.update (tapes_take tapes k' h_le) p v =
+    tapes_take (Function.update tapes ⟨p, by omega⟩ v) k' h_le := by
+  sorry
+
+/--
+Extend a sequence of tapes by adding more tapes at the end.
+Ignores the first `k₁` tapes of `extend_by` and uses the rest.
+-/
+@[simp]
+public abbrev tapes_extend_by
+  {γ : Type}
+  {k₁ k₂ : ℕ}
+  (tapes : Fin k₁ → γ)
+  (extend_by : Fin k₂ → γ)
+  (i : Fin k₂) : γ :=
+  if h : i < k₁ then tapes ⟨i, h⟩ else extend_by i
+
+@[simp, grind =]
+public lemma MultiTapeTM.extend_eval {k₁ k₂ : ℕ} (h_le : k₁ ≤ k₂)
+  (tm : MultiTapeTM k₁ Symbol)
+  {tapes : Fin k₂ → BiTape Symbol} :
+  (tm.extend h_le).eval tapes =
+    (tm.eval (tapes ⟨·, by omega⟩)).map (fun tapes' => tapes_extend_by tapes' tapes) := by
+  sorry
+
+end Turing
diff --git a/Cslib/Computability/Machines/MultiTapeTuring/WhileCombinator.lean b/Cslib/Computability/Machines/MultiTapeTuring/WhileCombinator.lean
new file mode 100644
index 00000000..97d46133
--- /dev/null
+++ b/Cslib/Computability/Machines/MultiTapeTuring/WhileCombinator.lean
@@ -0,0 +1,67 @@
+/-
+Copyright (c) 2026 Christian Reitwiessner. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Christian Reitwiessner
+-/
+
+module
+
+import Cslib.Foundations.Data.BiTape
+import Cslib.Foundations.Data.RelatesInSteps
+
+public import Cslib.Computability.Machines.MultiTapeTuring.Basic
+public import Cslib.Computability.Machines.MultiTapeTuring.ListEncoding
+public import Cslib.Computability.Machines.MultiTapeTuring.HeadStats
+
+namespace Turing
+
+namespace Routines
+
+variable [Inhabited α] [Fintype α]
+variable {k : ℕ}
+
+/--
+Repeatedly run a sub routine as long as a condition on the symbol
+at the head of tape `i` is true.
+-/
+public def doWhileSymbol (cond : Option α → Bool) (i : Fin k) (tm : MultiTapeTM k α) :
+    MultiTapeTM k α where
+  State := PUnit
+  stateFintype := inferInstance
+  q₀ := PUnit.unit
+  tr _ syms := sorry
+
+@[simp]
+public theorem doWhileSymbol_eval
+  (cond : Option α → Bool)
+  (i : Fin k)
+  (tm : MultiTapeTM k α)
+  (tapes_seq : ℕ → Fin k → BiTape α)
+  (h_transform : ∀ j, tm.eval (tapes_seq j) = .some (tapes_seq j.succ))
+  (h_stops : ∃ m, cond (tapes_seq m i).head = false) :
+  (doWhileSymbol cond i tm).eval (tapes_seq 0) = .some (tapes_seq (Nat.find h_stops)) := by
+  sorry
+
+/-- Repeatedly run a sub routine as long as the first word on tape `i` is non-empty.
+-/
+public def doWhile (i : Fin k) (tm : MultiTapeTM k (WithSep α)) :
+    MultiTapeTM k (WithSep α) where
+  State := PUnit
+  stateFintype := inferInstance
+  q₀ := PUnit.unit
+  tr _ syms := sorry
+
+@[simp]
+public theorem doWhile_eval_list
+  {i : Fin k}
+  {tm : MultiTapeTM k (WithSep α)}
+  {tapes : Fin k → List (List α)}
+  (h_halts : ∀ tapes, tm.HaltsOnLists tapes) :
+  (doWhile i tm).eval_list tapes =
+    ⟨∃ n, ((tm.eval_list_tot h_halts)^[n] tapes i).head?.getD [] = [],
+      fun h_loopEnds => (tm.eval_list_tot h_halts)^[Nat.find h_loopEnds] tapes⟩ := by
+  sorry
+
+end Routines
+
+end Turing
diff --git a/Cslib/Computability/Machines/MultiTapeTuring/WithTapes.lean b/Cslib/Computability/Machines/MultiTapeTuring/WithTapes.lean
new file mode 100644
index 00000000..209179db
--- /dev/null
+++ b/Cslib/Computability/Machines/MultiTapeTuring/WithTapes.lean
@@ -0,0 +1,127 @@
+/-
+Copyright (c) 2026 Christian Reitwiessner. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Christian Reitwiessner
+-/
+
+module
+
+public import Cslib.Computability.Machines.MultiTapeTuring.Basic
+public import Cslib.Computability.Machines.MultiTapeTuring.TapeExtension
+
+public import Mathlib.Logic.Equiv.Fintype
+public import Mathlib.Data.Finset.Sort
+
+namespace Turing
+
+variable [Inhabited α] [Fintype α]
+
+variable {k : ℕ}
+
+/--
+Permute tapes according to a bijection.
+-/
+public def MultiTapeTM.permute_tapes
+  (tm : MultiTapeTM k α) (σ : Equiv.Perm (Fin k)) : MultiTapeTM k α where
+  State := tm.State
+  stateFintype := tm.stateFintype
+  q₀ := tm.q₀
+  tr := fun q syms => match tm.tr q (syms ∘ σ) with
+    | (stmts, q') => (stmts ∘ σ.symm, q')
+
+--- General theorem: permuting tapes commutes with evaluation
+@[simp, grind =]
+public theorem MultiTapeTM.permute_tapes_eval
+  (tm : MultiTapeTM k α) (σ : Equiv.Perm (Fin k)) (tapes : Fin k → BiTape α) :
+  (tm.permute_tapes σ).eval tapes =
+    (tm.eval (tapes ∘ σ)).map (fun tapes' => tapes' ∘ σ.symm) := by
+  sorry
+
+private def findFin {k : ℕ} (p : Fin k → Prop) [DecidablePred p] (h : ∃ i, p i) : Fin k :=
+  (Finset.univ.filter p).min' (by
+    simp only [Finset.Nonempty, Finset.mem_filter, Finset.mem_univ, true_and]
+    exact h)
+
+def inj_to_perm {k₁ k₂ : ℕ} (f : Fin k₁ → Fin k₂) (h_inj : f.Injective) :
+  Equiv.Perm (Fin k₂)
+  -- non-computable version
+  --  let f' : {i : Fin k₂ // i < k₁} → Fin k₂ := fun ⟨i, _⟩ => f ⟨i, by omega⟩
+  --  have h_f'_inj : f'.Injective := by intro a b h; grind
+  --  (Equiv.ofInjective f' h_f'_inj).extendSubtype
+  where
+  toFun := sorry
+  invFun := sorry
+  left_inv := by sorry
+  right_inv := by sorry
+
+/--
+Change the order of the tapes of a Turing machine.
+Example: For a 2-tape Turing machine tm,
+`tm.with_tapes [2, 4].get (by grind)` is a 5-tape Turing machine whose tape 2 is
+the original machine's tape 0 and whose tape 4 is the original
+machine's tape 1
+Note that `f` is a function to `Fin k₂`, which means that integer literals
+automatically wrap. You have to be careful to make sure that the target machine
+has the right amount of tapes.
+-/
+public def MultiTapeTM.with_tapes {k₁ k₂ : ℕ}
+-- TODO use embedding instead?
+  (tm : MultiTapeTM k₁ α) (f : Fin k₁ → Fin k₂) (h_inj : f.Injective) : MultiTapeTM k₂ α :=
+  (tm.extend
+    (by simpa using Fintype.card_le_of_injective f h_inj)).permute_tapes (inj_to_perm f h_inj)
+
+-- TODO do not use `h.choose` here but rather assume that `f` is injective.
+
+/--
+Updates `tapes` by choosing elements from `tapes'` according to (the partial inverse of) `f`.
+-/
+public noncomputable def apply_updates
+  {γ : Type}
+  {k₁ k₂ : ℕ}
+  (tapes : Fin k₂ → γ)
+  (tapes' : Fin k₁ → γ)
+  (f : Fin k₁ → Fin k₂)
+  (i : Fin k₂) : γ :=
+  if h : ∃ j, f j = i then tapes' h.choose else tapes i
+
+@[simp, grind =]
+public lemma apply_updates_function_update_apply
+  {γ : Type}
+  {k₁ k₂ : ℕ}
+  {tapes : Fin k₂ → γ}
+  {f : Fin k₁ → Fin k₂}
+  (h_inj : f.Injective)
+  {t : Fin k₁}
+  {new_val : γ}
+  {i : Fin k₂} :
+  apply_updates tapes (Function.update (tapes ∘ f) t new_val) f i =
+    Function.update tapes (f t) new_val i := by
+  sorry
+
+@[simp, grind =]
+public lemma apply_updates_function_update
+  {γ : Type}
+  {k₁ k₂ : ℕ}
+  {tapes : Fin k₂ → γ}
+  {f : Fin k₁ → Fin k₂}
+  (h_inj : f.Injective)
+  {t : Fin k₁}
+  {new_val : γ} :
+  apply_updates tapes (Function.update (tapes ∘ f) t new_val) f =
+    Function.update tapes (f t) new_val := by
+  funext i
+  apply apply_updates_function_update_apply h_inj
+
+@[simp, grind =]
+public theorem MultiTapeTM.with_tapes_eval
+  {k₁ k₂ : ℕ}
+  {tm : MultiTapeTM k₁ α} {f : Fin k₁ → Fin k₂} {h_inj : f.Injective}
+  {tapes : Fin k₂ → BiTape α} :
+  (tm.with_tapes f h_inj).eval tapes =
+    (tm.eval (tapes ∘ f)).map
+      (fun tapes' => fun t => apply_updates tapes tapes' f t) := by
+  simp [with_tapes]
+  sorry
+
+
+end Turing