diff --git a/Complexity.lean b/Complexity.lean index 6dea70cad3..0234d7c3ee 100644 --- a/Complexity.lean +++ b/Complexity.lean @@ -36,6 +36,16 @@ public import Complexity.Algorithms.Streaming.CountMin.Basic public import Complexity.Algorithms.Streaming.CountMin.LinearHash public import Complexity.Algorithms.Streaming.CountMin.Probability public import Complexity.Algorithms.Streaming.CountMin.Program +public import Complexity.Algorithms.WordRAM.Dijkstra.BinaryHeap +public import Complexity.Algorithms.WordRAM.Dijkstra.ExactNotation +public import Complexity.Algorithms.WordRAM.Dijkstra.HeapProgram +public import Complexity.Algorithms.WordRAM.Dijkstra.HeapRelax +public import Complexity.Algorithms.WordRAM.Dijkstra.HeapSift +public import Complexity.Algorithms.WordRAM.Dijkstra.HeapStatement +public import Complexity.Algorithms.WordRAM.Dijkstra.Invariant +public import Complexity.Algorithms.WordRAM.Dijkstra.ListEncode +public import Complexity.Algorithms.WordRAM.Dijkstra.LogLinear +public import Complexity.Algorithms.WordRAM.Dijkstra.Statement public import Complexity.Algorithms.WordRAM.DynamicArray public import Complexity.Algorithms.WordRAM.ElementDistinctness public import Complexity.Algorithms.WordRAM.ElementDistinctnessRAM diff --git a/Complexity/Algorithms/WordRAM/Dijkstra/BinaryHeap.lean b/Complexity/Algorithms/WordRAM/Dijkstra/BinaryHeap.lean new file mode 100644 index 0000000000..51023949d2 --- /dev/null +++ b/Complexity/Algorithms/WordRAM/Dijkstra/BinaryHeap.lean @@ -0,0 +1,384 @@ +/- +Copyright (c) 2026 Renfei Zhou. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: Renfei Zhou +-/ + +module +-- upstream: candidate + +public import Complexity.Algorithms.WordRAM.Dijkstra.HeapProgram + +/-! +# Dijkstra's algorithm with a binary heap: the heap + +This file states the invariant of the binary min-heap of +`Cslib.Algorithms.WordRAM.Dijkstra.heapDijkstra` and proves how its operations change it. The heap +lives in a memory `mem` with three arrays past a base `Db`, for `N` vertices: the keys +`mem (Db + v)`, the heap entries `mem (Db + N + i)` for the positions `1 ≤ i ≤ sz`, and the +positions `mem (Db + 2 N + 1 + v)`, with `0` for a vertex not in the heap. The key of a position +is the key of the vertex there (`Cslib.Algorithms.WordRAM.Dijkstra.heapKey`). + +The invariant (`Cslib.Algorithms.WordRAM.Dijkstra.HeapInv`) has a shape part +(`Cslib.Algorithms.WordRAM.Dijkstra.HeapShape`: the positions array is the inverse of the heap, +and vertices in the heap have nonzero keys) and an order part: the key of every position is at +least the key of its parent. A sift-up from position `k` repairs a heap ordered except between `k` +and its parent (`Cslib.Algorithms.WordRAM.Dijkstra.UpOrdered`), and a sift-down from `k` repairs a +heap ordered except between `k` and its children (`Cslib.Algorithms.WordRAM.Dijkstra.DownOrdered`). +Insertion and decrease-key leave the first kind of defect, and removal of the root the second. + +## Main definitions + +* `Cslib.Algorithms.WordRAM.Dijkstra.HeapShape`, `Cslib.Algorithms.WordRAM.Dijkstra.HeapInv`: the + heap invariant. +* `Cslib.Algorithms.WordRAM.Dijkstra.UpOrdered`, `Cslib.Algorithms.WordRAM.Dijkstra.DownOrdered`: + heap order with one defect. + +## Main results + +* `Cslib.Algorithms.WordRAM.Dijkstra.HeapInv.push`, + `Cslib.Algorithms.WordRAM.Dijkstra.HeapInv.decrease`, + `Cslib.Algorithms.WordRAM.Dijkstra.HeapInv.extract`: insertion, decrease-key and removal of the + root leave a heap with one defect. +* `Cslib.Algorithms.WordRAM.Dijkstra.HeapInv.root_le`: the root has the least key. + +## References + +* [T. H. Cormen, C. E. Leiserson, R. L. Rivest and C. Stein, *Introduction to + Algorithms*][CLRS2009] +-/ + +@[expose] public section + +namespace Cslib.Algorithms.WordRAM.Dijkstra + +open Cslib.WordRAM HeapReg + +/-- The key of heap position `i`: the key of the vertex at position `i`. -/ +def heapKey (N Db : ℕ) (mem : ℕ → ℕ) (i : ℕ) : ℕ := mem (Db + mem (Db + N + i)) + +/-- **The shape of the heap.** The positions `1, …, sz` hold vertices below `N`, the positions +array is the inverse of the heap and is `0` off the heap, and vertices in the heap have nonzero +keys. -/ +structure HeapShape (N Db sz : ℕ) (mem : ℕ → ℕ) : Prop where + /-- The heap holds at most `N` entries. -/ + size_le : sz ≤ N + /-- The heap holds vertices. -/ + vertex_lt : ∀ i, 1 ≤ i → i ≤ sz → mem (Db + N + i) < N + /-- The position of the vertex at position `i` is `i`. -/ + pos_vertex : ∀ i, 1 ≤ i → i ≤ sz → mem (Db + 2 * N + 1 + mem (Db + N + i)) = i + /-- Positions are at most the size. -/ + pos_le : ∀ v < N, mem (Db + 2 * N + 1 + v) ≤ sz + /-- A vertex with a nonzero position is at that position. -/ + vertex_pos : ∀ v < N, mem (Db + 2 * N + 1 + v) ≠ 0 → + mem (Db + N + mem (Db + 2 * N + 1 + v)) = v + /-- Vertices in the heap have nonzero keys. -/ + key_ne : ∀ v < N, mem (Db + 2 * N + 1 + v) ≠ 0 → mem (Db + v) ≠ 0 + +/-- **Heap order**: the key of every position `2 ≤ i ≤ sz` is at least that of its parent. -/ +def Ordered (f : ℕ → ℕ) (sz : ℕ) : Prop := ∀ i, 2 ≤ i → i ≤ sz → f (i / 2) ≤ f i + +/-- Heap order except between position `k` and its parent, where the parent of `k` is at most +the children of `k`: the state of a sift-up at `k`. -/ +def UpOrdered (f : ℕ → ℕ) (sz k : ℕ) : Prop := + ∀ i, 2 ≤ i → i ≤ sz → (i ≠ k → f (i / 2) ≤ f i) ∧ (i / 2 = k → 2 ≤ k → f (k / 2) ≤ f i) + +/-- Heap order except between position `k` and its children, where the parent of `k` is at most +the children of `k`: the state of a sift-down at `k`. -/ +def DownOrdered (f : ℕ → ℕ) (sz k : ℕ) : Prop := + ∀ i, 2 ≤ i → i ≤ sz → (i / 2 ≠ k → f (i / 2) ≤ f i) ∧ (i / 2 = k → 2 ≤ k → f (k / 2) ≤ f i) + +/-- **The heap invariant**: the shape and the heap order. -/ +structure HeapInv (N Db sz : ℕ) (mem : ℕ → ℕ) : Prop extends HeapShape N Db sz mem where + /-- Heap order. -/ + ordered : Ordered (heapKey N Db mem) sz + +variable {N Db S : ℕ} {mem : ℕ → ℕ} + +/-- A heap that misses a vertex has room for it. -/ +theorem HeapShape.size_lt (h : HeapShape N Db S mem) {v : ℕ} (hv : v < N) + (hp : mem (Db + 2 * N + 1 + v) = 0) : S < N := by + have hcard := Finset.card_le_card_of_injOn (fun i => mem (Db + N + i)) + (s := Finset.Icc 1 S) (t := (Finset.range N).erase v) (fun i hi => ?_) + (fun i hi j hj he => ?_) + · rw [Nat.card_Icc, Finset.card_erase_of_mem (Finset.mem_range.2 hv), Finset.card_range] + at hcard + omega + · rw [Finset.mem_coe, Finset.mem_Icc] at hi + rw [Finset.mem_coe, Finset.mem_erase, Finset.mem_range] + refine ⟨fun he => ?_, h.vertex_lt i hi.1 hi.2⟩ + have := h.pos_vertex i hi.1 hi.2 + simp only at he + rw [he, hp] at this + omega + · rw [Finset.mem_coe, Finset.mem_Icc] at hi hj + have h1 := h.pos_vertex i hi.1 hi.2 + have h2 := h.pos_vertex j hj.1 hj.2 + simp only at he + rw [he] at h1 + omega + +/-- **The root has the least key** among the vertices in the heap. -/ +theorem HeapInv.root_le (h : HeapInv N Db S mem) {v : ℕ} (hv : v < N) + (hp : mem (Db + 2 * N + 1 + v) ≠ 0) : mem (Db + mem (Db + N + 1)) ≤ mem (Db + v) := by + have key : ∀ i, 1 ≤ i → i ≤ S → heapKey N Db mem 1 ≤ heapKey N Db mem i := by + intro i + induction i using Nat.strong_induction_on with + | _ i ih => + intro h1 h2 + rcases Nat.lt_or_ge i 2 with hi | hi + · obtain rfl : i = 1 := by omega + exact le_rfl + · exact (ih (i / 2) (by omega) (by omega) (by omega)).trans (h.ordered i hi h2) + have := key _ (Nat.one_le_iff_ne_zero.2 hp) (h.pos_le v hv) + simpa only [heapKey, h.vertex_pos v hv hp] using this + +/-- **Insertion.** Giving a vertex `v` out of the heap the nonzero key `c` and appending it at +position `S + 1` leaves a heap ordered except at `S + 1`, which holds `v`. -/ +theorem HeapInv.push (h : HeapInv N Db S mem) {v c : ℕ} (hv : v < N) + (hp : mem (Db + 2 * N + 1 + v) = 0) (hc : c ≠ 0) : + HeapShape N Db (S + 1) (Function.update (Function.update (Function.update mem (Db + v) c) + (Db + N + (S + 1)) v) (Db + 2 * N + 1 + v) (S + 1)) ∧ + UpOrdered (heapKey N Db (Function.update (Function.update (Function.update mem (Db + v) c) + (Db + N + (S + 1)) v) (Db + 2 * N + 1 + v) (S + 1))) (S + 1) (S + 1) := by + have hSN := h.size_lt hv hp + have hne : ∀ i, 1 ≤ i → i ≤ S → mem (Db + N + i) ≠ v := fun i h1 h2 he => by + have := h.pos_vertex i h1 h2 + rw [he, hp] at this + omega + set m' := Function.update (Function.update (Function.update mem (Db + v) c) + (Db + N + (S + 1)) v) (Db + 2 * N + 1 + v) (S + 1) with hm' + have hkey : ∀ i, 1 ≤ i → i ≤ S → heapKey N Db m' i = heapKey N Db mem i := by + intro i h1 h2 + have := hne i h1 h2 + have := h.vertex_lt i h1 h2 + simp only [hm', heapKey, Function.update_apply] + split_ifs <;> omega + refine ⟨⟨by omega, fun i h1 h2 => ?_, fun i h1 h2 => ?_, fun w hw => ?_, fun w hw hw0 => ?_, + fun w hw hw0 => ?_⟩, fun i h1 h2 => ⟨fun hi => ?_, fun hi _ => by omega⟩⟩ + · rcases Nat.lt_or_ge i (S + 1) with hi | hi + · have := h.vertex_lt i h1 (by omega) + simp only [hm', Function.update_apply] + split_ifs <;> omega + · simp only [hm', Function.update_apply] + split_ifs <;> omega + · rcases Nat.lt_or_ge i (S + 1) with hi | hi + · have := h.vertex_lt i h1 (by omega) + have := hne i h1 (by omega) + have := h.pos_vertex i h1 (by omega) + simp only [hm', Function.update_apply] + split_ifs <;> omega + · simp only [hm', Function.update_apply] + split_ifs <;> omega + · have := h.pos_le w hw + simp only [hm', Function.update_apply] + split_ifs <;> omega + · by_cases hwv : w = v + · subst hwv + simp only [hm', Function.update_apply] + split_ifs <;> omega + · have hP : m' (Db + 2 * N + 1 + w) = mem (Db + 2 * N + 1 + w) := by + simp only [hm', Function.update_apply] + split_ifs <;> omega + rw [hP] at hw0 ⊢ + have := h.pos_le w hw + have := h.vertex_pos w hw hw0 + simp only [hm', Function.update_apply] + split_ifs <;> omega + · by_cases hwv : w = v + · subst hwv + simp only [hm', Function.update_apply] + split_ifs <;> omega + · have hP : m' (Db + 2 * N + 1 + w) = mem (Db + 2 * N + 1 + w) := by + simp only [hm', Function.update_apply] + split_ifs <;> omega + rw [hP] at hw0 + have := h.key_ne w hw hw0 + simp only [hm', Function.update_apply] + split_ifs <;> omega + · rw [hkey i (by omega) (by omega), hkey (i / 2) (by omega) (by omega)] + exact h.ordered i h1 (by omega) + +/-- **Decrease-key.** Lowering the key of a vertex `v` in the heap to the nonzero `c` leaves a +heap ordered except at the position of `v`. -/ +theorem HeapInv.decrease (h : HeapInv N Db S mem) {v c : ℕ} (hv : v < N) + (hp : mem (Db + 2 * N + 1 + v) ≠ 0) (hc : c ≠ 0) (hle : c ≤ mem (Db + v)) : + HeapShape N Db S (Function.update mem (Db + v) c) ∧ + UpOrdered (heapKey N Db (Function.update mem (Db + v) c)) S (mem (Db + 2 * N + 1 + v)) := by + have hk1 : 1 ≤ mem (Db + 2 * N + 1 + v) := Nat.one_le_iff_ne_zero.2 hp + have hkS := h.pos_le v hv + have hHk := h.vertex_pos v hv hp + have hval : ∀ j, heapKey N Db (Function.update mem (Db + v) c) j = + if mem (Db + N + j) = v then c else heapKey N Db mem j := by + intro j + have hin : Function.update mem (Db + v) c (Db + N + j) = mem (Db + N + j) := + Function.update_of_ne (by omega) _ _ + simp only [heapKey, hin, Function.update_apply] + split_ifs <;> first | rfl | omega + have hle' : ∀ j, heapKey N Db (Function.update mem (Db + v) c) j ≤ heapKey N Db mem j := by + intro j + rw [hval] + split_ifs with h1 + · simp only [heapKey, h1]; exact hle + · exact le_rfl + have heq : ∀ j, 1 ≤ j → j ≤ S → j ≠ mem (Db + 2 * N + 1 + v) → + heapKey N Db (Function.update mem (Db + v) c) j = heapKey N Db mem j := by + intro j h1 h2 hj + rw [hval] + split_ifs with he + · have := h.pos_vertex j h1 h2 + rw [he] at this + omega + · rfl + refine ⟨⟨h.size_le, fun i h1 h2 => ?_, fun i h1 h2 => ?_, fun w hw => ?_, fun w hw hw0 => ?_, + fun w hw hw0 => ?_⟩, fun i h1 h2 => ⟨fun hi => ?_, fun hi h2k => ?_⟩⟩ + · rw [Function.update_of_ne (by omega)] + exact h.vertex_lt i h1 h2 + · have := h.vertex_lt i h1 h2 + rw [Function.update_of_ne (by omega), Function.update_of_ne (by omega)] + exact h.pos_vertex i h1 h2 + · rw [Function.update_of_ne (by omega)] + exact h.pos_le w hw + · rw [Function.update_of_ne (by omega)] at hw0 ⊢ + have := h.pos_le w hw + rw [Function.update_of_ne (by omega)] + exact h.vertex_pos w hw hw0 + · rw [Function.update_of_ne (by omega)] at hw0 + rw [Function.update_apply] + split_ifs + · exact hc + · exact h.key_ne w hw hw0 + · rw [heq i (by omega) h2 hi] + exact (hle' _).trans (h.ordered i h1 h2) + · rw [heq i (by omega) h2 (by omega)] + exact (hle' _).trans ((h.ordered _ h2k hkS).trans (hi ▸ h.ordered i h1 h2)) + +/-- **Removal of the root**, from a heap of at least two entries: clearing the position of the +root `U` and moving the last entry `L` to the root leaves a heap of `S - 1` entries ordered except +between the root and its children, which holds exactly the vertices of the old heap but `U`. -/ +theorem HeapInv.extract (h : HeapInv N Db S mem) (hS : 2 ≤ S) : + HeapShape N Db (S - 1) (Function.update (Function.update + (Function.update mem (Db + 2 * N + 1 + mem (Db + N + 1)) 0) (Db + N + 1) + (mem (Db + N + S))) (Db + 2 * N + 1 + mem (Db + N + S)) 1) ∧ + DownOrdered (heapKey N Db (Function.update (Function.update + (Function.update mem (Db + 2 * N + 1 + mem (Db + N + 1)) 0) (Db + N + 1) + (mem (Db + N + S))) (Db + 2 * N + 1 + mem (Db + N + S)) 1)) (S - 1) 1 ∧ + ∀ v < N, (Function.update (Function.update + (Function.update mem (Db + 2 * N + 1 + mem (Db + N + 1)) 0) (Db + N + 1) + (mem (Db + N + S))) (Db + 2 * N + 1 + mem (Db + N + S)) 1 (Db + 2 * N + 1 + v) = 0 ↔ + mem (Db + 2 * N + 1 + v) = 0 ∨ v = mem (Db + N + 1)) := by + have hSN := h.size_le + have hU := h.vertex_lt 1 le_rfl (by omega) + have hL := h.vertex_lt S (by omega) le_rfl + have hPU := h.pos_vertex 1 le_rfl (by omega) + have hPL := h.pos_vertex S (by omega) le_rfl + have hUL : mem (Db + N + 1) ≠ mem (Db + N + S) := fun e => by + rw [e, hPL] at hPU + omega + -- The middle positions `2, …, S - 1` hold neither the root nor the last vertex. + have hmid : ∀ j, 2 ≤ j → j ≤ S - 1 → mem (Db + N + j) ≠ mem (Db + N + 1) ∧ + mem (Db + N + j) ≠ mem (Db + N + S) ∧ mem (Db + N + j) < N ∧ + mem (Db + 2 * N + 1 + mem (Db + N + j)) = j := by + intro j h1 h2 + have hj := h.pos_vertex j (by omega) (by omega) + refine ⟨fun e => ?_, fun e => ?_, h.vertex_lt j (by omega) (by omega), hj⟩ + · rw [e, hPU] at hj + omega + · rw [e, hPL] at hj + omega + -- Any other vertex has its position among the middle positions or is out of the heap. + have hother : ∀ w < N, w ≠ mem (Db + N + 1) → w ≠ mem (Db + N + S) → + mem (Db + 2 * N + 1 + w) ≤ S - 1 ∧ mem (Db + 2 * N + 1 + w) ≠ 1 := by + intro w hw h1 h2 + have := h.pos_le w hw + by_cases h0 : mem (Db + 2 * N + 1 + w) = 0 + · omega + have hv := h.vertex_pos w hw h0 + refine ⟨?_, fun e => h1 (by rw [e] at hv; exact hv.symm)⟩ + by_contra hc + exact h2 (by rw [show mem (Db + 2 * N + 1 + w) = S by omega] at hv; exact hv.symm) + set m' := Function.update (Function.update + (Function.update mem (Db + 2 * N + 1 + mem (Db + N + 1)) 0) (Db + N + 1) + (mem (Db + N + S))) (Db + 2 * N + 1 + mem (Db + N + S)) 1 with hm' + -- The position of a vertex in the new heap. + have hpos : ∀ w < N, m' (Db + 2 * N + 1 + w) = if w = mem (Db + N + S) then 1 + else if w = mem (Db + N + 1) then 0 else mem (Db + 2 * N + 1 + w) := by + intro w hw + simp only [hm', Function.update_apply] + split_ifs <;> omega + have hkey : ∀ j, 2 ≤ j → j ≤ S - 1 → heapKey N Db m' j = heapKey N Db mem j := by + intro j h1 h2 + have := hmid j h1 h2 + simp only [hm', heapKey, Function.update_apply] + split_ifs <;> omega + refine ⟨⟨by omega, fun i h1 h2 => ?_, fun i h1 h2 => ?_, fun w hw => ?_, fun w hw hw0 => ?_, + fun w hw hw0 => ?_⟩, fun i h1 h2 => ⟨fun hi => ?_, fun _ h => by omega⟩, fun w hw => ?_⟩ + · have := h.vertex_lt i h1 (by omega) + simp only [hm', Function.update_apply] + split_ifs <;> omega + · rcases Nat.lt_or_ge i 2 with hi | hi + · obtain rfl : i = 1 := by omega + simp only [hm', Function.update_apply] + split_ifs <;> omega + · have := hmid i hi h2 + simp only [hm', Function.update_apply] + split_ifs <;> omega + · rw [hpos w hw] + split_ifs with h1 h2 + · omega + · omega + · exact (hother w hw h2 h1).1 + · rw [hpos w hw] at hw0 ⊢ + split_ifs at hw0 ⊢ with h1 h2 + · simp only [hm', Function.update_apply] + split_ifs <;> omega + · exact absurd rfl hw0 + · have := hother w hw h2 h1 + have := h.vertex_pos w hw hw0 + simp only [hm', Function.update_apply] + split_ifs <;> omega + · rw [hpos w hw] at hw0 + have hk : m' (Db + w) = mem (Db + w) := by + simp only [hm', Function.update_apply] + split_ifs <;> omega + rw [hk] + split_ifs at hw0 with h1 h2 + · subst h1 + exact h.key_ne _ hw (by omega) + · exact absurd rfl hw0 + · exact h.key_ne w hw hw0 + · rw [hkey i h1 h2, hkey (i / 2) (by omega) (by omega)] + exact h.ordered i h1 (by omega) + · rw [hpos w hw] + split_ifs with h1 h2 + · subst h1 + refine iff_of_false id ?_ + rintro (h0 | h0) + · omega + · exact hUL h0.symm + · omega + · omega + +/-- **Removal of the only entry**: clearing the position of the root of a heap of one entry +leaves the empty heap. -/ +theorem HeapInv.extract_one (h : HeapInv N Db 1 mem) : + HeapInv N Db 0 (Function.update mem (Db + 2 * N + 1 + mem (Db + N + 1)) 0) := by + have hU := h.vertex_lt 1 le_rfl le_rfl + have hz : ∀ w < N, Function.update mem (Db + 2 * N + 1 + mem (Db + N + 1)) 0 + (Db + 2 * N + 1 + w) = 0 := by + intro w hw + rw [Function.update_apply] + split_ifs with he + · rfl + · have := h.pos_le w hw + by_contra h0 + have hv := h.vertex_pos w hw h0 + rw [show mem (Db + 2 * N + 1 + w) = 1 by omega] at hv + omega + exact ⟨⟨Nat.zero_le _, fun i h1 h2 => by omega, fun i h1 h2 => by omega, + fun w hw => (hz w hw).le, fun w hw hw0 => absurd (hz w hw) hw0, + fun w hw hw0 => absurd (hz w hw) hw0⟩, fun i h1 h2 => by omega⟩ + +end Cslib.Algorithms.WordRAM.Dijkstra + +end diff --git a/Complexity/Algorithms/WordRAM/Dijkstra/ExactNotation.lean b/Complexity/Algorithms/WordRAM/Dijkstra/ExactNotation.lean new file mode 100644 index 0000000000..5e844c6af2 --- /dev/null +++ b/Complexity/Algorithms/WordRAM/Dijkstra/ExactNotation.lean @@ -0,0 +1,150 @@ +/- +Copyright (c) 2026 Renfei Zhou. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: Renfei Zhou +-/ + +module +-- upstream: candidate + +public import Complexity.Computability.WordRAM.Notation + +/-! +# Notation for structured word-RAM programs over the exact basis + +`ram_exact[ s₁; s₂; … ]` is `ram[ s₁; s₂; … ]` over the exact basis +`Cslib.WordRAM.Basis.exact`: it has the same statements, but a binary operation `d ← i - j`, +`d ← i / j`, … is the instruction `Cslib.WordRAM.Op.exactBin`, which faults where its exact value +is undefined (subtraction below zero, division and remainder by zero), and `d ← width%` is +`Cslib.WordRAM.Op.exactWidth`. The tactic `ram_exact_eval` is `ram_eval` for such programs. + +The notation, with its expansions and the tactic, is scoped to +`Cslib.Algorithms.WordRAM.Dijkstra`. Its statements are those of the scoped `wramStmt` syntax of +`Cslib.WordRAM`, so they parse only where both namespaces are open. + +## Main definitions + +* `ram_exact[ s₁; s₂; … ]`: structured word-RAM programs over the exact basis. +* `ram_exact_eval`: symbolic execution of such programs. + +## Tactics + +* `ram_exact_eval`, `ram_exact_eval n`: as `ram_eval`. A subtraction runs when the simp lemmas + show that its right operand is at most its left one, and a division or remainder when they + show that its divisor is nonzero. + +## References + +The statements are those of the `ram[` notation (`WordRAM.Notation`), whose forms follow the +word-RAM notation of [Algolean][Algolean2026]. +-/ + +public section + +namespace Cslib.Algorithms.WordRAM.Dijkstra + +open scoped Cslib.WordRAM + +/-- A structured word-RAM program over the exact basis, written as a block of statements. -/ +scoped syntax "ram_exact[" sepBy1(wramStmt, "; ") "]" : term + +/-- Internal: a nested block of statements, as a structured program over the exact basis. -/ +scoped syntax "ram_exact_block%" Cslib.WordRAM.wramBlock : term + +/-- Internal: a single statement, as a structured program over the exact basis. -/ +scoped syntax "ram_exact_stmt%" wramStmt : term + +scoped macro_rules + | `(ram_exact[ $s:wramStmt ]) => `(ram_exact_stmt% $s) + | `(ram_exact[ $s:wramStmt; $ss;* ]) => + `(Cslib.WordRAM.Code.seq (ram_exact_stmt% $s) (ram_exact[ $ss;* ])) + | `(ram_exact_block% [ $ss;* ]) => `(ram_exact[ $ss;* ]) + +scoped macro_rules + | `(ram_exact_stmt% $d:ident ← #$v) => `(Cslib.WordRAM.Code.op (.const $d $v)) + | `(ram_exact_stmt% $d:ident ← $s:ident) => `(Cslib.WordRAM.Code.op (.copy $d $s)) + | `(ram_exact_stmt% $d:ident ← $i:ident + $j:ident) => + `(Cslib.WordRAM.Code.op (.exactBin .add $d $i $j)) + | `(ram_exact_stmt% $d:ident ← $i:ident - $j:ident) => + `(Cslib.WordRAM.Code.op (.exactBin .sub $d $i $j)) + | `(ram_exact_stmt% $d:ident ← $i:ident * $j:ident) => + `(Cslib.WordRAM.Code.op (.exactBin .mul $d $i $j)) + | `(ram_exact_stmt% $d:ident ← $i:ident / $j:ident) => + `(Cslib.WordRAM.Code.op (.exactBin .div $d $i $j)) + | `(ram_exact_stmt% $d:ident ← $i:ident % $j:ident) => + `(Cslib.WordRAM.Code.op (.exactBin .mod $d $i $j)) + | `(ram_exact_stmt% $d:ident ← $i:ident &&& $j:ident) => + `(Cslib.WordRAM.Code.op (.exactBin .and $d $i $j)) + | `(ram_exact_stmt% $d:ident ← $i:ident ||| $j:ident) => + `(Cslib.WordRAM.Code.op (.exactBin .or $d $i $j)) + | `(ram_exact_stmt% $d:ident ← $i:ident ^^^ $j:ident) => + `(Cslib.WordRAM.Code.op (.exactBin .xor $d $i $j)) + | `(ram_exact_stmt% $d:ident ← $i:ident <<< $j:ident) => + `(Cslib.WordRAM.Code.op (.exactBin .shiftLeft $d $i $j)) + | `(ram_exact_stmt% $d:ident ← $i:ident >>> $j:ident) => + `(Cslib.WordRAM.Code.op (.exactBin .shiftRight $d $i $j)) + | `(ram_exact_stmt% $d:ident ← *$a:ident) => `(Cslib.WordRAM.Code.op (.load $d $a)) + | `(ram_exact_stmt% *$a:ident ← $s:ident) => `(Cslib.WordRAM.Code.op (.store $a $s)) + | `(ram_exact_stmt% $d:ident ← width%) => `(Cslib.WordRAM.Code.op (.exactWidth $d)) + | `(ram_exact_stmt% if $i:ident = $j:ident then $p:wramBlock else $n:wramBlock) => + `(Cslib.WordRAM.Code.branch (.eq $i $j) (ram_exact_block% $p) (ram_exact_block% $n)) + | `(ram_exact_stmt% if $i:ident < $j:ident then $p:wramBlock else $n:wramBlock) => + `(Cslib.WordRAM.Code.branch (.lt $i $j) (ram_exact_block% $p) (ram_exact_block% $n)) + | `(ram_exact_stmt% while $i:ident = $j:ident do $b:wramBlock) => + `(Cslib.WordRAM.Code.loop (.eq $i $j) (ram_exact_block% $b)) + | `(ram_exact_stmt% while $i:ident < $j:ident do $b:wramBlock) => + `(Cslib.WordRAM.Code.loop (.lt $i $j) (ram_exact_block% $b)) + | `(ram_exact_stmt% ⟪$c⟫) => `(($c : Cslib.WordRAM.Code .exact _)) + +/-- Symbolically execute structured code over the exact basis on a goal `ExecWithin w c x y T` or +`Exec w c x y T`, allowing each loop the given number of iterations (default `0`), and deciding +with the given simp lemmas which written numbers fit a word, which subtractions stay at or above +zero, and which divisors are nonzero. -/ +scoped syntax (name := ramExactEval) + "ram_exact_eval" (ppSpace num)? (" [" Lean.Parser.Tactic.simpArg,* "]")? : tactic + +scoped macro_rules + | `(tactic| ram_exact_eval) => `(tactic| ram_exact_eval 0 []) + | `(tactic| ram_exact_eval [$args,*]) => `(tactic| ram_exact_eval 0 [$args,*]) + | `(tactic| ram_exact_eval $n:num) => `(tactic| ram_exact_eval $n []) + | `(tactic| ram_exact_eval $n:num [$args,*]) => `(tactic| + first + | (apply Cslib.WordRAM.ExecWithin.of_eval (fuel := $n) + · set_option linter.unusedSimpArgs false in + simp +decide only [Cslib.WordRAM.Code.eval, Cslib.WordRAM.Code.evalLoop, + Cslib.WordRAM.Op.eval_const, Cslib.WordRAM.Op.eval_copy, Cslib.WordRAM.Op.eval_alu, + Cslib.WordRAM.Op.eval_load, Cslib.WordRAM.Op.eval_store, + Cslib.WordRAM.Op.eval_exactBin, Cslib.WordRAM.Op.eval_exactWidth, + Cslib.WordRAM.Basis.exact_interp_bin, Cslib.WordRAM.Basis.exact_interp_width, + Cslib.WordRAM.Data.put, Cslib.WordRAM.fits_zero, Cslib.WordRAM.BinOp.evalExact, + Cslib.WordRAM.Cond.holds_eq, Cslib.WordRAM.Cond.holds_lt, + Cslib.WordRAM.Data.set_regs, Cslib.WordRAM.Data.set_mem, + Cslib.WordRAM.Data.write_regs, Cslib.WordRAM.Data.write_mem, Function.update_self, + Function.update_apply, bind, Option.map_some, Option.bind_some, Option.map_none, + Option.bind_none, ↓reduceIte, Nat.reduceAdd, Nat.reduceMul, Nat.reduceSub, + Nat.reduceDiv, Nat.reduceMod, Nat.reducePow, Nat.reduceAnd, Nat.reduceOr, + Nat.reduceXor, Nat.reduceShiftLeft, Nat.reduceShiftRight, Nat.reduceLT, + Nat.reduceLeDiff, Nat.reduceEqDiff, $args,*] + try rfl + · try omega) + | (apply Cslib.WordRAM.Exec.of_eval (fuel := $n) + set_option linter.unusedSimpArgs false in + simp +decide only [Cslib.WordRAM.Code.eval, Cslib.WordRAM.Code.evalLoop, + Cslib.WordRAM.Op.eval_const, Cslib.WordRAM.Op.eval_copy, Cslib.WordRAM.Op.eval_alu, + Cslib.WordRAM.Op.eval_load, Cslib.WordRAM.Op.eval_store, + Cslib.WordRAM.Op.eval_exactBin, Cslib.WordRAM.Op.eval_exactWidth, + Cslib.WordRAM.Basis.exact_interp_bin, Cslib.WordRAM.Basis.exact_interp_width, + Cslib.WordRAM.Data.put, Cslib.WordRAM.fits_zero, Cslib.WordRAM.BinOp.evalExact, + Cslib.WordRAM.Cond.holds_eq, Cslib.WordRAM.Cond.holds_lt, + Cslib.WordRAM.Data.set_regs, Cslib.WordRAM.Data.set_mem, + Cslib.WordRAM.Data.write_regs, Cslib.WordRAM.Data.write_mem, Function.update_self, + Function.update_apply, bind, Option.map_some, Option.bind_some, Option.map_none, + Option.bind_none, ↓reduceIte, Nat.reduceAdd, Nat.reduceMul, Nat.reduceSub, + Nat.reduceDiv, Nat.reduceMod, Nat.reducePow, Nat.reduceAnd, Nat.reduceOr, + Nat.reduceXor, Nat.reduceShiftLeft, Nat.reduceShiftRight, Nat.reduceLT, + Nat.reduceLeDiff, Nat.reduceEqDiff, $args,*] + try rfl)) + +end Cslib.Algorithms.WordRAM.Dijkstra + +end diff --git a/Complexity/Algorithms/WordRAM/Dijkstra/HeapProgram.lean b/Complexity/Algorithms/WordRAM/Dijkstra/HeapProgram.lean new file mode 100644 index 0000000000..3f98b4f5a7 --- /dev/null +++ b/Complexity/Algorithms/WordRAM/Dijkstra/HeapProgram.lean @@ -0,0 +1,237 @@ +/- +Copyright (c) 2026 Renfei Zhou. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: Renfei Zhou +-/ + +module +-- upstream: candidate + +public import Complexity.Algorithms.WordRAM.Dijkstra.HeapStatement +public import Complexity.Algorithms.WordRAM.Dijkstra.ExactNotation + +/-! +# Dijkstra's algorithm with a binary heap: the program + +This file writes Dijkstra's algorithm [Dijkstra][Dijkstra1959] [Cormen et al., Section +24.3][CLRS2009] on adjacency lists as a structured word-RAM program over the exact basis +`Cslib.WordRAM.Basis.exact` and 27 registers, with a binary min-heap [Cormen et al., Chapter +6][CLRS2009] that supports decrease-key through an array of positions. + +The program reads the input of `Cslib.Algorithms.WordRAM.Dijkstra.ListInstance.encode`: `n` at +address `1`, the source at address `2`, the offset of vertex `u` at address `u + 3`, and edge `j` +as its target at address `n + 4 + 2 j` and its weight at the next address. With `m` the number of +edges and `Db = n + 2 m + 4` the first address past the input, it keeps three arrays: + +* the distance words `D v` at `Db + v`, which are the output cells: `0` for an undiscovered + vertex and `d + 1` for tentative distance `d`; +* the heap `H i` at `Db + n + i` for the positions `1 ≤ i ≤ sz`, ordered by the keys `D (H i)`; +* the positions `P v` at `Db + 2 n + 1 + v`: the position of `v` in the heap, and `0` when `v` is + not in it. + +Memory past the input starts at `0`, so every vertex starts undiscovered and out of the heap, and +the arrays need no initialization. A vertex is visited once it has left the heap: its distance +word is nonzero and its position is `0`. + +Each round extracts the root `u` of the heap and scans the edges out of `u`. An edge to an +undiscovered vertex inserts it with its new key; an edge that lowers the key of a vertex in the +heap moves it up. The rounds stop when the heap is empty. + +On the exact basis a subtraction below zero or a division by zero faults. The program divides +only by the constant `2`, and subtracts only `1` from the heap size, which is positive then. + +## Main definitions + +* `Cslib.Algorithms.WordRAM.Dijkstra.siftUp`, `Cslib.Algorithms.WordRAM.Dijkstra.siftDown`: the + two heap repairs. +* `Cslib.Algorithms.WordRAM.Dijkstra.extractMin`: removal of the root of the heap. +* `Cslib.Algorithms.WordRAM.Dijkstra.pushNew`: insertion of a newly discovered vertex into the + heap. +* `Cslib.Algorithms.WordRAM.Dijkstra.relaxEdges`: the scan of the edges out of a vertex. +* `Cslib.Algorithms.WordRAM.Dijkstra.heapDijkstra`: Dijkstra's algorithm with a binary heap. + +## References + +* [E. W. Dijkstra, *A note on two problems in connexion with graphs*][Dijkstra1959] +* [T. H. Cormen, C. E. Leiserson, R. L. Rivest and C. Stein, *Introduction to + Algorithms*][CLRS2009] +-/ + +@[expose] public section + +namespace Cslib.Algorithms.WordRAM.Dijkstra + +open Cslib.WordRAM + +/-! ### Register names -/ + +namespace HeapReg + +/-- The constant `0`. -/ +abbrev zero : Fin 27 := 0 +/-- The constant `1`. -/ +abbrev one : Fin 27 := 1 +/-- The constant `2`. -/ +abbrev two : Fin 27 := 2 +/-- The constant `3`. -/ +abbrev three : Fin 27 := 3 +/-- The number `n` of vertices. -/ +abbrev n : Fin 27 := 4 +/-- The address `n + 4` of the first edge. -/ +abbrev eb : Fin 27 := 5 +/-- The base `Db` of the distance words. -/ +abbrev db : Fin 27 := 6 +/-- The base `Db + n` of the heap. -/ +abbrev hb : Fin 27 := 7 +/-- The base `Db + 2 n + 1` of the positions. -/ +abbrev pb : Fin 27 := 8 +/-- The size of the heap. -/ +abbrev sz : Fin 27 := 9 +/-- The vertex of the current round. -/ +abbrev u : Fin 27 := 10 +/-- The distance word of `u`. -/ +abbrev best : Fin 27 := 11 +/-- The index of the current edge. -/ +abbrev j : Fin 27 := 12 +/-- The index past the last edge out of `u`. -/ +abbrev jend : Fin 27 := 13 +/-- The target of the current edge. -/ +abbrev tv : Fin 27 := 14 +/-- The weight of the current edge. -/ +abbrev tc : Fin 27 := 15 +/-- The candidate distance word through the current edge. -/ +abbrev cand : Fin 27 := 16 +/-- The distance word of the target. -/ +abbrev dv : Fin 27 := 17 +/-- The heap position being repaired. -/ +abbrev k : Fin 27 := 18 +/-- The parent or the chosen child of position `k`. -/ +abbrev p : Fin 27 := 19 +/-- The vertex at position `k`. -/ +abbrev vk : Fin 27 := 20 +/-- The vertex at position `p`. -/ +abbrev vp : Fin 27 := 21 +/-- The key of `vk`. -/ +abbrev dk : Fin 27 := 22 +/-- The key of `vp`. -/ +abbrev dp : Fin 27 := 23 +/-- A scratch register, mostly for addresses. -/ +abbrev t : Fin 27 := 24 +/-- The heap size plus one, during a sift-down. -/ +abbrev szp : Fin 27 := 25 +/-- The left child `2 k` of position `k`, during a sift-down. -/ +abbrev l : Fin 27 := 26 + +end HeapReg + +open HeapReg + +/-- Swap the heap entries at positions `k` and `p`, which hold the vertices `vk` and `vp`, and +update their positions. -/ +def swapCode : Code .exact 27 := ram_exact[ + t ← hb + k; *t ← vp; + t ← hb + p; *t ← vk; + t ← pb + vp; *t ← k; + t ← pb + vk; *t ← p ] + +/-- One step of the sift-up: swap position `k` with its parent if the parent has the larger key, +and otherwise stop by setting `k` to `1`. -/ +def siftUpBody : Code .exact 27 := ram_exact[ + p ← k / two; + t ← hb + k; vk ← *t; t ← db + vk; dk ← *t; + t ← hb + p; vp ← *t; t ← db + vp; dp ← *t; + if dk < dp then [ ⟪swapCode⟫; k ← p ] else [ k ← one ] ] + +/-- **Sift-up** from position `k`. -/ +def siftUp : Code .exact 27 := ram_exact[ while one < k do [ ⟪siftUpBody⟫ ] ] + +/-- The first half of a sift-down step: put the child of least key of position `k`, among the +left child `l` and the right child `l + 1` if it is in the heap, in `p`, its vertex in `vp` and +its key in `dp`. -/ +def siftDownChoose : Code .exact 27 := ram_exact[ + p ← l; + t ← hb + l; vp ← *t; t ← db + vp; dp ← *t; + t ← l + one; + if t < szp then [ + t ← hb + t; vk ← *t; t ← db + vk; dk ← *t; + if dk < dp then [ p ← l + one; vp ← vk; dp ← dk ] else [ p ← p ] ] else [ p ← p ] ] + +/-- The second half of a sift-down step: swap position `k` with the position `p` chosen by +`siftDownChoose` and move to it if its key `dp` is smaller, and otherwise stop by setting `l` to +`szp`. -/ +def siftDownSink : Code .exact 27 := ram_exact[ + t ← hb + k; vk ← *t; t ← db + vk; dk ← *t; + if dp < dk then [ ⟪swapCode⟫; k ← p; l ← k + k ] else [ l ← szp ] ] + +/-- One step of the sift-down: swap position `k` with its child of least key if that key is +smaller, and otherwise stop by setting `l` to `szp`. -/ +def siftDownBody : Code .exact 27 := ram_exact[ ⟪siftDownChoose⟫; ⟪siftDownSink⟫ ] + +/-- **Sift-down** from position `k`. -/ +def siftDown : Code .exact 27 := ram_exact[ + l ← k + k; szp ← sz + one; + while l < szp do [ ⟪siftDownBody⟫ ] ] + +/-- **Extract-min**: remove the root `u` of the heap, with distance word `best`, move the last +entry to the root and sift it down. -/ +def extractMin : Code .exact 27 := ram_exact[ + t ← hb + one; u ← *t; + t ← db + u; best ← *t; + t ← pb + u; *t ← zero; + t ← hb + sz; vk ← *t; + sz ← sz - one; + if zero < sz then [ + t ← hb + one; *t ← vk; + t ← pb + vk; *t ← one; + k ← one; ⟪siftDown⟫ ] else [ k ← k ] ] + +/-- **Insert** the vertex `tv`, whose distance word is already set, as a new last entry, and sift +it up. -/ +def pushNew : Code .exact 27 := ram_exact[ + sz ← sz + one; + t ← hb + sz; *t ← tv; + t ← pb + tv; *t ← sz; + k ← sz; ⟪siftUp⟫ ] + +/-- Relax edge `j` out of `u`: read its target `tv` and weight `tc`, and if `best + tc` improves +the distance word of `tv`, store it and insert `tv` into the heap or move it up (decrease-key). -/ +def edgeBody : Code .exact 27 := ram_exact[ + t ← j + j; t ← t + eb; tv ← *t; t ← t + one; tc ← *t; j ← j + one; + cand ← best + tc; + t ← db + tv; dv ← *t; + if dv = zero then [ *t ← cand; ⟪pushNew⟫ ] else [ + if cand < dv then [ *t ← cand; t ← pb + tv; k ← *t; ⟪siftUp⟫ ] else [ t ← t ] ] ] + +/-- The scan of the edges `j, …, jend - 1`. -/ +def relaxEdges : Code .exact 27 := ram_exact[ while j < jend do [ ⟪edgeBody⟫ ] ] + +/-- One round: extract the root `u` of the heap and relax the edges out of it, which are those +between its offset, at address `u + 3`, and the next offset. -/ +def heapRound : Code .exact 27 := ram_exact[ + ⟪extractMin⟫; + t ← u + three; j ← *t; + t ← t + one; jend ← *t; + ⟪relaxEdges⟫ ] + +/-- The setup: load the constants and `n`, compute the array bases from `n` and the number `m` +of edges, which is the offset at address `n + 3`, give the source the distance word `1` and make +it the only heap entry. -/ +def heapSetup : Code .exact 27 := ram_exact[ + one ← #1; two ← #2; three ← #3; + t ← #1; n ← *t; + t ← #2; u ← *t; + t ← n + three; t ← *t; + db ← t + t; db ← db + n; db ← db + three; db ← db + one; + eb ← n + three; eb ← eb + one; + hb ← db + n; pb ← hb + n; pb ← pb + one; + t ← db + u; *t ← one; + t ← hb + one; *t ← u; + t ← pb + u; *t ← one; + sz ← one ] + +/-- **Dijkstra's algorithm with a binary heap**, on adjacency lists. -/ +def heapDijkstra : Code .exact 27 := ram_exact[ ⟪heapSetup⟫; while zero < sz do [ ⟪heapRound⟫ ] ] + +end Cslib.Algorithms.WordRAM.Dijkstra + +end diff --git a/Complexity/Algorithms/WordRAM/Dijkstra/HeapRelax.lean b/Complexity/Algorithms/WordRAM/Dijkstra/HeapRelax.lean new file mode 100644 index 0000000000..9ce5246682 --- /dev/null +++ b/Complexity/Algorithms/WordRAM/Dijkstra/HeapRelax.lean @@ -0,0 +1,362 @@ +/- +Copyright (c) 2026 Renfei Zhou. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: Renfei Zhou +-/ + +module +-- upstream: candidate + +public import Complexity.Algorithms.WordRAM.Dijkstra.HeapSift +public import Complexity.Algorithms.WordRAM.Dijkstra.ListEncode + +/-! +# Dijkstra's algorithm with a binary heap: the scan of the edges + +This file proves what one run of `Cslib.Algorithms.WordRAM.Dijkstra.relaxEdges` does on the +edges out of a vertex `U` of an adjacency-list instance: within `O(deg U · (⌊log₂ n⌋ + 1) + 1)` +instructions it applies `Cslib.Algorithms.WordRAM.Dijkstra.edgesWord` to every distance word, and +keeps the heap invariant, with the heap holding the discovered vertices that are not visited. + +## Main results + +* `Cslib.Algorithms.WordRAM.Dijkstra.relaxEdges_spec`: the scan of the edges out of a vertex. + +## References + +* [E. W. Dijkstra, *A note on two problems in connexion with graphs*][Dijkstra1959] +* [T. H. Cormen, C. E. Leiserson, R. L. Rivest and C. Stein, *Introduction to + Algorithms*][CLRS2009] +-/ + +@[expose] public section + +namespace Cslib.Algorithms.WordRAM.Dijkstra + +open Cslib.WordRAM HeapReg + +/-- The straight-line start of `edgeBody`: read the edge and the distance word of its target. -/ +private def edgePrefix : Code .exact 27 := ram_exact[ + t ← j + j; t ← t + eb; tv ← *t; t ← t + one; tc ← *t; j ← j + one; + cand ← best + tc; + t ← db + tv; dv ← *t ] + +/-- The straight-line start of an insertion: store the key and append the vertex to the heap. -/ +private def insPrefix : Code .exact 27 := ram_exact[ + *t ← cand; sz ← sz + one; t ← hb + sz; *t ← tv; t ← pb + tv; *t ← sz; k ← sz ] + +/-- The straight-line start of a decrease-key: store the key and read the position. -/ +private def decPrefix : Code .exact 27 := ram_exact[ *t ← cand; t ← pb + tv; k ← *t ] + +/-- `edgeBody` split into its straight-line parts and the sift-ups. -/ +private theorem edgeBody_eq : edgeBody = appendCode edgePrefix (.branch (.eq dv zero) + (appendCode insPrefix siftUp) + (.branch (.lt cand dv) (appendCode decPrefix siftUp) ram_exact[ t ← t ])) := rfl + +/-- `edgeBody` writes only its scratch registers and `sz`. -/ +private theorem edgeBody_preserves : + ∀ i : Fin 27, i ∉ [j, tv, tc, cand, dv, t, k, p, vk, vp, dk, dp, sz] → + edgeBody.Preserves i := by + decide + +/-- Relaxing through the first `i + 1` edges is relaxing through the first `i` and then +through edge `i`. -/ +private theorem edgesWord_take_succ {N b d i : ℕ} {es : List (Fin N × ℕ)} (v : Fin N) + (hi : i < es.length) : + edgesWord b (es.take (i + 1)) v d = if es[i].1 = v then + edgeWord b es[i].2 (edgesWord b (es.take i) v d) else edgesWord b (es.take i) v d := by + rw [List.take_add_one, List.getElem?_eq_getElem hi] + simp only [edgesWord, List.foldl_append, Option.toList_some, List.foldl_cons, List.foldl_nil] + +/-- A nonzero distance word that no edge of the list improves stays. -/ +private theorem edgesWord_of_le {N b d : ℕ} {es : List (Fin N × ℕ)} {v : Fin N} (hd : d ≠ 0) + (hle : ∀ c, (v, c) ∈ es → d ≤ b + c) : edgesWord b es v d = d := by + induction es with + | nil => rfl + | cons e es ih => + have h' : (if e.1 = v then edgeWord b e.2 d else d) = d := by + split_ifs with he + · have := hle e.2 (by rw [← he]; simp) + simp [edgeWord, hd, Nat.not_lt.2 this] + · rfl + simp only [edgesWord, List.foldl_cons] at ih ⊢ + rw [h'] + exact ih fun c hc => hle c (List.mem_cons_of_mem _ hc) + +/-- **One edge.** The relaxation of an edge to `V` of weight `C`, with `B ≠ 0` the distance word +of its tail, from a heap that holds every discovered vertex whose word the edge could improve: +it advances `j`, relaxes the word of `V` and keeps the heap invariant and the visited vertices. +-/ +private theorem edgeBody_spec {G : ListInstance} {w B J S C : ℕ} {V : Fin G.n} {y : Data 27} + (hzero : y.regs zero = 0) (hone : y.regs one = 1) (htwo : y.regs two = 2) + (heb : y.regs eb = G.n + 4) (hdb : y.regs db = G.dbase) (hhb : y.regs hb = G.dbase + G.n) + (hpb : y.regs pb = G.dbase + 2 * G.n + 1) (hsz : y.regs sz = S) (hbest : y.regs best = B) + (hB : B ≠ 0) (hj : y.regs j = J) (hJ : J < G.m) + (hfit : Fits w (G.dbase + 3 * G.n + 1)) (hfitC : Fits w (B + C)) + (hT : y.mem (J + J + (G.n + 4)) = V) (hW : y.mem (J + J + (G.n + 4) + 1) = C) + (hheap : HeapInv G.n G.dbase S y.mem) + (hvis : y.mem (G.dbase + V) ≠ 0 → y.mem (G.dbase + 2 * G.n + 1 + V) = 0 → + y.mem (G.dbase + V) ≤ B + C) : + ∃ z, ExecWithin w edgeBody y z (30 * Nat.log 2 G.n + 30) ∧ + z.regs j = J + 1 ∧ z.regs sz ≤ S + 1 ∧ HeapInv G.n G.dbase (z.regs sz) z.mem ∧ + (∀ v : Fin G.n, z.mem (G.dbase + v) = + if V = v then edgeWord B C (y.mem (G.dbase + v)) else y.mem (G.dbase + v)) ∧ + (∀ v : Fin G.n, (z.mem (G.dbase + v) ≠ 0 ∧ z.mem (G.dbase + 2 * G.n + 1 + v) = 0) ↔ + (y.mem (G.dbase + v) ≠ 0 ∧ y.mem (G.dbase + 2 * G.n + 1 + v) = 0)) ∧ + (∀ a < G.dbase, z.mem a = y.mem a) := by + have hDb : G.dbase = G.n + 2 * G.m + 4 := rfl + have hV := V.isLt + have f1 : Fits w (J + J) := hfit.of_le (by omega) + have f2 : Fits w (J + J + (G.n + 4)) := hfit.of_le (by omega) + have f3 : Fits w (J + J + (G.n + 4) + 1) := hfit.of_le (by omega) + have f4 : Fits w (J + 1) := hfit.of_le (by omega) + have f5 : Fits w (G.dbase + V) := hfit.of_le (by omega) + have f8 : Fits w (G.dbase + 2 * G.n + 1 + V) := hfit.of_le (by omega) + let y₁ := ((((((((y.set t (J + J)).set t (J + J + (G.n + 4))).set tv V).set t + (J + J + (G.n + 4) + 1)).set tc C).set j (J + 1)).set cand (B + C)).set t + (G.dbase + V)).set dv (y.mem (G.dbase + V)) + have e₁ : ExecWithin w edgePrefix y y₁ 9 := by + ram_exact_eval [edgePrefix, y₁, hone, heb, hdb, hbest, hj, hT, hW, hfitC, f1, f2, f3, f4, f5] + have hlog {K : ℕ} (hK : K ≤ G.n) : Nat.log 2 K ≤ Nat.log 2 G.n := Nat.log_mono_right hK + by_cases hD0 : y.mem (G.dbase + V) = 0 + · -- `V` is undiscovered: insert it + have hP : y.mem (G.dbase + 2 * G.n + 1 + V) = 0 := by + by_contra hP + exact hheap.key_ne V hV hP hD0 + have hSN : S < G.n := hheap.size_lt hV hP + have f6 : Fits w (S + 1) := hfit.of_le (by omega) + have f7 : Fits w (G.dbase + G.n + (S + 1)) := hfit.of_le (by omega) + let y₂ := ((((((y₁.write (G.dbase + V) (B + C)).set sz (S + 1)).set t + (G.dbase + G.n + (S + 1))).write (G.dbase + G.n + (S + 1)) V).set t + (G.dbase + 2 * G.n + 1 + V)).write (G.dbase + 2 * G.n + 1 + V) (S + 1)).set k (S + 1) + have e₂ : ExecWithin w insPrefix y₁ y₂ 7 := by + ram_exact_eval [insPrefix, y₁, y₂, hone, hhb, hpb, hsz, f6, f7, f8] + have hm₂ : y₂.mem = Function.update (Function.update (Function.update y.mem + (G.dbase + V) (B + C)) (G.dbase + G.n + (S + 1)) V) + (G.dbase + 2 * G.n + 1 + V) (S + 1) := rfl + obtain ⟨hsh, hord⟩ := hheap.push (c := B + C) hV hP (by omega) + rw [← hm₂] at hsh hord + obtain ⟨z, e₃, hzr, hzh, hzp, hzm⟩ := siftUp_spec (S := S + 1) (K := S + 1) (x := y₂) + (by simp +decide [y₂, y₁, hone]) (by simp +decide [y₂, y₁, htwo]) + (by simp +decide [y₂, y₁, hdb]) (by simp +decide [y₂, y₁, hhb]) + (by simp +decide [y₂, y₁, hpb]) (by simp +decide [y₂]) (by omega) le_rfl hfit hsh hord + have hzj : z.regs j = J + 1 := by rw [hzr j (by decide)]; simp +decide [y₂, y₁] + have hzsz : z.regs sz = S + 1 := by rw [hzr sz (by decide)]; simp +decide [y₂] + have hzw (v : Fin G.n) (hv : V ≠ v) : z.mem (G.dbase + v) = y.mem (G.dbase + v) := by + have hvv : (V : ℕ) ≠ v := fun h => hv (Fin.ext h) + have := v.isLt + rw [hzm _ (Or.inl (by omega)), hm₂, Function.update_of_ne (by omega), + Function.update_of_ne (by omega), Function.update_of_ne (by omega)] + refine ⟨z, ?_, hzj, by omega, hzsz ▸ hzh, ?_, ?_, ?_⟩ + · rw [edgeBody_eq] + refine (execWithin_append e₁ + (ExecWithin.branch_pos ?_ (execWithin_append e₂ e₃))).mono ?_ + · simp +decide [Cond.holds_eq, y₁, hzero, hD0] + · have := hlog (K := S + 1) (by omega) + omega + · intro v + by_cases hv : V = v + · subst hv + simp only [↓reduceIte] + rw [hzm _ (Or.inl (by omega)), hm₂, Function.update_of_ne (by omega), + Function.update_of_ne (by omega), Function.update_self, hD0] + simp [edgeWord] + · simp only [hv, ↓reduceIte] + rw [hzw v hv] + · intro v + by_cases hv : V = v + · subst hv + have h₁ : z.mem (G.dbase + 2 * G.n + 1 + V) ≠ 0 := by + rw [Ne, hzp V hV, hm₂, Function.update_self] + omega + exact ⟨fun h => absurd h.2 h₁, fun h => absurd hD0 h.1⟩ + · have hvv : (V : ℕ) ≠ v := fun h => hv (Fin.ext h) + have := v.isLt + rw [hzw v hv, hzp v v.isLt, hm₂, Function.update_of_ne (by omega), + Function.update_of_ne (by omega), Function.update_of_ne (by omega)] + · intro a ha + rw [hzm _ (Or.inl (by omega)), hm₂, Function.update_of_ne (by omega), + Function.update_of_ne (by omega), Function.update_of_ne (by omega)] + by_cases hlt : B + C < y.mem (G.dbase + V) + · -- `V` is in the heap and its key decreases + have hP : y.mem (G.dbase + 2 * G.n + 1 + V) ≠ 0 := fun hP => by + have := hvis hD0 hP + omega + have hPS : y.mem (G.dbase + 2 * G.n + 1 + V) ≤ S := hheap.pos_le V hV + have hne : G.dbase + 2 * G.n + 1 + V ≠ G.dbase + V := by omega + let y₂ := ((y₁.write (G.dbase + V) (B + C)).set t (G.dbase + 2 * G.n + 1 + V)).set k + (y.mem (G.dbase + 2 * G.n + 1 + V)) + have e₂ : ExecWithin w decPrefix y₁ y₂ 3 := by + ram_exact_eval [decPrefix, y₁, y₂, hpb, f8, hne] + have hm₂ : y₂.mem = Function.update y.mem (G.dbase + V) (B + C) := rfl + obtain ⟨hsh, hord⟩ := hheap.decrease (c := B + C) hV hP (by omega) hlt.le + rw [← hm₂] at hsh hord + obtain ⟨z, e₃, hzr, hzh, hzp, hzm⟩ := siftUp_spec (S := S) + (K := y.mem (G.dbase + 2 * G.n + 1 + V)) (x := y₂) + (by simp +decide [y₂, y₁, hone]) (by simp +decide [y₂, y₁, htwo]) + (by simp +decide [y₂, y₁, hdb]) (by simp +decide [y₂, y₁, hhb]) + (by simp +decide [y₂, y₁, hpb]) (by simp +decide [y₂]) (by omega) hPS hfit hsh hord + have hzj : z.regs j = J + 1 := by rw [hzr j (by decide)]; simp +decide [y₂, y₁] + have hzsz : z.regs sz = S := by rw [hzr sz (by decide)]; simp +decide [y₂, y₁, hsz] + have hzw (v : Fin G.n) (hv : V ≠ v) : z.mem (G.dbase + v) = y.mem (G.dbase + v) := by + have hvv : (V : ℕ) ≠ v := fun h => hv (Fin.ext h) + have := v.isLt + rw [hzm _ (Or.inl (by omega)), hm₂, Function.update_of_ne (by omega)] + refine ⟨z, ?_, hzj, hzsz.le.trans (Nat.le_succ S), hzsz ▸ hzh, ?_, ?_, ?_⟩ + · rw [edgeBody_eq] + refine (execWithin_append e₁ (ExecWithin.branch_neg ?_ (ExecWithin.branch_pos ?_ + (execWithin_append e₂ e₃)))).mono ?_ + · simp +decide [Cond.holds_eq, y₁, hzero, hD0] + · simp +decide [Cond.holds_lt, y₁, hlt] + · have := hlog (K := y.mem (G.dbase + 2 * G.n + 1 + V)) (hPS.trans hheap.size_le) + omega + · intro v + by_cases hv : V = v + · subst hv + simp only [↓reduceIte] + rw [hzm _ (Or.inl (by omega)), hm₂, Function.update_self] + simp [edgeWord, hD0, hlt] + · simp only [hv, ↓reduceIte] + rw [hzw v hv] + · intro v + by_cases hv : V = v + · subst hv + have h₁ : z.mem (G.dbase + 2 * G.n + 1 + V) ≠ 0 := by + rw [Ne, hzp V hV, hm₂, Function.update_of_ne hne] + exact hP + exact ⟨fun h => absurd h.2 h₁, fun h => absurd h.2 hP⟩ + · have hvv : (V : ℕ) ≠ v := fun h => hv (Fin.ext h) + have := v.isLt + rw [hzw v hv, hzp v v.isLt, hm₂, Function.update_of_ne (by omega)] + · intro a ha + rw [hzm _ (Or.inl (by omega)), hm₂, Function.update_of_ne (by omega)] + · -- no change + have e₂ : ExecWithin (B := .exact) w ram_exact[ t ← t ] y₁ (y₁.set t (G.dbase + V)) 1 := by + ram_exact_eval [y₁] + refine ⟨y₁.set t (G.dbase + V), ?_, by simp +decide [y₁], by simp +decide [y₁, hsz], + by simpa +decide [y₁, hsz] using hheap, ?_, fun v => Iff.rfl, fun a _ => rfl⟩ + · rw [edgeBody_eq] + refine (execWithin_append e₁ + (ExecWithin.branch_neg ?_ (ExecWithin.branch_neg ?_ e₂))).mono ?_ + · simp +decide [Cond.holds_eq, y₁, hzero, hD0] + · simp +decide [Cond.holds_lt, y₁, hlt] + · omega + · intro v + by_cases hv : V = v + · subst hv + simp [y₁, edgeWord, hD0, hlt] + · simp [y₁, hv] + +/-- The invariant of the scan of the edges out of `U`, started from `x`: the edges before `j` +have been relaxed. -/ +private structure ScanInv (G : ListInstance) (U : Fin G.n) (B S₀ : ℕ) (x y : Data 27) : + Prop where + /-- The registers other than the scratch registers and `sz` keep their contents. -/ + frame : ∀ i ∉ [j, tv, tc, cand, dv, t, k, p, vk, vp, dk, dp, sz], y.regs i = x.regs i + /-- The index is at least the offset of `U`. -/ + j_ge : G.offset U ≤ y.regs j + /-- The index is at most the offset of `U + 1`. -/ + j_le : y.regs j ≤ G.offset ((U : ℕ) + 1) + /-- Each relaxed edge inserts at most one vertex. -/ + sz_le : y.regs sz ≤ S₀ + (y.regs j - G.offset U) + /-- The heap invariant. -/ + heap : HeapInv G.n G.dbase (y.regs sz) y.mem + /-- The distance words are relaxed through the edges before `j`. -/ + word : ∀ v : Fin G.n, y.mem (G.dbase + v) = + edgesWord B ((G.adj U).take (y.regs j - G.offset U)) v (x.mem (G.dbase + v)) + /-- The visited vertices stay. -/ + visited : ∀ v : Fin G.n, + (y.mem (G.dbase + v) ≠ 0 ∧ y.mem (G.dbase + 2 * G.n + 1 + v) = 0) ↔ + (x.mem (G.dbase + v) ≠ 0 ∧ x.mem (G.dbase + 2 * G.n + 1 + v) = 0) + /-- The input stays. -/ + input : ∀ a < G.dbase, y.mem a = x.mem a + +/-- **The scan of the edges out of `U`.** Started with the index `j` and the end `jend` at the +offsets of `U` and `U + 1`, the distance word `B ≠ 0` of `U` in `best`, the input intact below +`Db = G.dbase` and a heap that holds the discovered vertices other than the visited ones, the +scan runs within `deg U · (30 ⌊log₂ n⌋ + 32) + 2` instructions. It replaces every distance word +by its relaxation through the edges out of `U`, inserts at most `deg U` vertices into the heap and +keeps the heap invariant, the visited vertices, the input and the registers other than its +scratch registers and `sz`. + +A vertex is visited when it is discovered (its distance word is nonzero) but not in the heap (its +position is `0`). The scan assumes that no edge out of `U` improves a visited vertex. -/ +theorem relaxEdges_spec {G : ListInstance} {w S₀ B : ℕ} (U : Fin G.n) {x : Data 27} + (hzero : x.regs zero = 0) (hone : x.regs one = 1) (htwo : x.regs two = 2) + (heb : x.regs eb = G.n + 4) (hdb : x.regs db = G.dbase) (hhb : x.regs hb = G.dbase + G.n) + (hpb : x.regs pb = G.dbase + 2 * G.n + 1) (hsz : x.regs sz = S₀) (hbest : x.regs best = B) + (hB : B ≠ 0) (hj : x.regs j = G.offset U) (hjend : x.regs jend = G.offset ((U : ℕ) + 1)) + (hfit : Fits w (G.dbase + 3 * G.n + 1)) (hfitC : ∀ e ∈ G.adj U, Fits w (B + e.2)) + (hinput : ∀ a < G.dbase, x.mem a = (Data.init G.encode : Data 27).mem a) + (hheap : HeapInv G.n G.dbase S₀ x.mem) + (hvis : ∀ v : Fin G.n, x.mem (G.dbase + v) ≠ 0 → x.mem (G.dbase + 2 * G.n + 1 + v) = 0 → + ∀ c, (v, c) ∈ G.adj U → x.mem (G.dbase + v) ≤ B + c) : + ∃ y, ExecWithin w relaxEdges x y ((G.adj U).length * (30 * Nat.log 2 G.n + 32) + 2) ∧ + (∀ i ∉ [j, tv, tc, cand, dv, t, k, p, vk, vp, dk, dp, sz], y.regs i = x.regs i) ∧ + y.regs sz ≤ S₀ + (G.adj U).length ∧ + HeapInv G.n G.dbase (y.regs sz) y.mem ∧ + (∀ v : Fin G.n, y.mem (G.dbase + v) = edgesWord B (G.adj U) v (x.mem (G.dbase + v))) ∧ + (∀ v : Fin G.n, (y.mem (G.dbase + v) ≠ 0 ∧ y.mem (G.dbase + 2 * G.n + 1 + v) = 0) ↔ + (x.mem (G.dbase + v) ≠ 0 ∧ x.mem (G.dbase + 2 * G.n + 1 + v) = 0)) ∧ + (∀ a < G.dbase, y.mem a = x.mem a) := by + have hDb : G.dbase = G.n + 2 * G.m + 4 := rfl + have hdeg : G.offset ((U : ℕ) + 1) = G.offset U + (G.adj U).length := G.offset_succ U.isLt + have hm : G.offset ((U : ℕ) + 1) ≤ G.m := G.offset_n ▸ G.offset_mono (by omega) + have hI₀ : ScanInv G U B S₀ x x := by + refine ⟨fun _ _ => rfl, hj.ge, hj ▸ G.offset_mono (by omega), by omega, hsz ▸ hheap, + fun v => ?_, fun v => Iff.rfl, fun _ _ => rfl⟩ + rw [hj, Nat.sub_self, List.take_zero] + rfl + have hbody : ∀ y, ScanInv G U B S₀ x y → (Cond.lt j jend).Holds y → + ∃ z, ExecWithin w edgeBody y z (30 * Nat.log 2 G.n + 30) ∧ ScanInv G U B S₀ x z ∧ + G.offset ((U : ℕ) + 1) - z.regs j < G.offset ((U : ℕ) + 1) - y.regs j := by + intro y hy hc + have hjend' : y.regs jend = G.offset ((U : ℕ) + 1) := (hy.frame jend (by decide)).trans hjend + have hlt : y.regs j < G.offset ((U : ℕ) + 1) := by + simpa only [Cond.holds_lt, hjend'] using hc + obtain ⟨i, hji⟩ : ∃ i, y.regs j = G.offset U + i := + ⟨y.regs j - G.offset U, by have := hy.j_ge; omega⟩ + have hi : i < (G.adj U).length := by omega + have hword (v : Fin G.n) : y.mem (G.dbase + v) = + edgesWord B ((G.adj U).take i) v (x.mem (G.dbase + v)) := by + rw [hy.word v, hji, Nat.add_sub_cancel_left] + have haddr : y.regs j + y.regs j + (G.n + 4) = 2 * (G.offset U + i) + G.n + 4 := by omega + have hT : y.mem (y.regs j + y.regs j + (G.n + 4)) = ((G.adj U)[i].1 : ℕ) := by + rw [hy.input _ (by omega), hinput _ (by omega), haddr, G.init_target U hi] + have hW : y.mem (y.regs j + y.regs j + (G.n + 4) + 1) = (G.adj U)[i].2 := by + rw [hy.input _ (by omega), hinput _ (by omega), haddr, G.init_weight U hi] + have hmem := List.getElem_mem hi + have hvisy : y.mem (G.dbase + (G.adj U)[i].1) ≠ 0 → + y.mem (G.dbase + 2 * G.n + 1 + (G.adj U)[i].1) = 0 → + y.mem (G.dbase + (G.adj U)[i].1) ≤ B + (G.adj U)[i].2 := by + intro h₁ h₂ + obtain ⟨hx₁, hx₂⟩ := (hy.visited _).1 ⟨h₁, h₂⟩ + rw [hword, edgesWord_of_le hx₁ fun c hc => hvis _ hx₁ hx₂ c (List.mem_of_mem_take hc)] + exact hvis _ hx₁ hx₂ _ hmem + obtain ⟨z, he, hzj, hzsz, hzh, hzw, hzv, hzi⟩ := edgeBody_spec + ((hy.frame zero (by decide)).trans hzero) ((hy.frame one (by decide)).trans hone) + ((hy.frame two (by decide)).trans htwo) ((hy.frame eb (by decide)).trans heb) + ((hy.frame db (by decide)).trans hdb) ((hy.frame hb (by decide)).trans hhb) + ((hy.frame pb (by decide)).trans hpb) rfl ((hy.frame best (by decide)).trans hbest) hB + rfl (by omega) hfit (hfitC _ hmem) hT hW hy.heap hvisy + refine ⟨z, he, ⟨?_, by omega, by omega, by have := hy.sz_le; omega, hzh, fun v => ?_, + fun v => (hzv v).trans (hy.visited v), fun a ha => (hzi a ha).trans (hy.input a ha)⟩, + by omega⟩ + · intro r hr + rw [he.regs_eq (edgeBody_preserves r hr), hy.frame r hr] + · rw [hzw v, hword v, hzj, hji, Nat.add_assoc, Nat.add_sub_cancel_left, + edgesWord_take_succ v hi] + obtain ⟨y, he, hy, hexit⟩ := ExecWithin.loop_invariant (cond := .lt j jend) (body := edgeBody) + (ScanInv G U B S₀ x) (fun y => G.offset ((U : ℕ) + 1) - y.regs j) hbody hI₀ + have hjend' : y.regs jend = G.offset ((U : ℕ) + 1) := (hy.frame jend (by decide)).trans hjend + have hlen : y.regs j - G.offset U = (G.adj U).length := by + have := hy.j_le + simp only [Cond.holds_lt, hjend'] at hexit + omega + refine ⟨y, he.mono ?_, hy.frame, hlen ▸ hy.sz_le, hy.heap, fun v => ?_, hy.visited, hy.input⟩ + · rw [hj, show G.offset ((U : ℕ) + 1) - G.offset U = (G.adj U).length by omega] + · rw [hy.word v, hlen, List.take_length] + +end Cslib.Algorithms.WordRAM.Dijkstra + +end diff --git a/Complexity/Algorithms/WordRAM/Dijkstra/HeapSift.lean b/Complexity/Algorithms/WordRAM/Dijkstra/HeapSift.lean new file mode 100644 index 0000000000..28ad89341e --- /dev/null +++ b/Complexity/Algorithms/WordRAM/Dijkstra/HeapSift.lean @@ -0,0 +1,863 @@ +/- +Copyright (c) 2026 Renfei Zhou. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: Renfei Zhou +-/ + +module +-- upstream: candidate + +public import Complexity.Algorithms.WordRAM.Dijkstra.BinaryHeap + +/-! +# Dijkstra's algorithm with a binary heap: the heap operations + +This file proves what the heap operations of `Cslib.Algorithms.WordRAM.Dijkstra.heapDijkstra` +do on a heap with the invariant `Cslib.Algorithms.WordRAM.Dijkstra.HeapInv`, and how long they +take. A sift-up or a sift-down halves or doubles the repaired position at every step, so each +takes `O(⌊log₂ sz⌋ + 1)` instructions. + +## Main definitions + +* `Cslib.Algorithms.WordRAM.Dijkstra.appendCode`: one code block after another, nested as + `ram_exact[…]` nests sequences, so that a program splits into blocks proved one at a time. + +## Main results + +* `Cslib.Algorithms.WordRAM.Dijkstra.siftUp_spec`, + `Cslib.Algorithms.WordRAM.Dijkstra.siftDown_spec`: the two repairs. +* `Cslib.Algorithms.WordRAM.Dijkstra.extractMin_spec`: removal of the root. + +## References + +* [T. H. Cormen, C. E. Leiserson, R. L. Rivest and C. Stein, *Introduction to + Algorithms*][CLRS2009] +-/ + +@[expose] public section + +namespace Cslib.Algorithms.WordRAM.Dijkstra + +open Cslib.WordRAM HeapReg + +variable {N Db S : ℕ} + +/-! ### Sequencing -/ + +/-- `appendCode c d` runs `c` and then `d`, as a sequence nested to the right, as the notation +`ram_exact[…]` nests it. -/ +def appendCode {B : Basis} {r : ℕ} : Code B r → Code B r → Code B r + | .seq c₁ c₂, d => .seq c₁ (appendCode c₂ d) + | c, d => .seq c d + +/-- Running `c` and then `d`. -/ +theorem execWithin_append {B : Basis} {r w T₁ T₂ : ℕ} {c d : Code B r} {x y z : Data r} + (h₁ : ExecWithin w c x y T₁) (h₂ : ExecWithin w d y z T₂) : + ExecWithin w (appendCode c d) x z (T₁ + T₂) := by + induction c generalizing x T₁ with + | seq c₁ c₂ _ ih => + obtain ⟨n, hn, he⟩ := h₁ + cases he with + | seq he₁ he₂ => exact (he₁.within.seq (ih he₂.within)).mono (by omega) + | op o => exact h₁.seq h₂ + | branch c pos neg _ _ => exact h₁.seq h₂ + | loop c body _ => exact h₁.seq h₂ + | rand dst src => exact h₁.seq h₂ + +/-! ### Heap order with one defect -/ + +section Order + +variable {f g : ℕ → ℕ} {a c : ℕ} + +/-- A heap ordered except above `a`, where the parent of `a` is at most `a`, is ordered. -/ +private theorem UpOrdered.ordered (h : UpOrdered f S a) + (hle : 2 ≤ a → a ≤ S → f (a / 2) ≤ f a) : Ordered f S := by + intro i hi hiS + by_cases hia : i = a + · subst hia + exact hle hi hiS + · exact (h i hi hiS).1 hia + +/-- An ordered heap is ordered except above the root. -/ +private theorem Ordered.upOrdered_one (h : Ordered f S) : UpOrdered f S 1 := + fun i hi hiS => ⟨fun _ => h i hi hiS, fun _ h2 => absurd h2 (by omega)⟩ + +/-- **A sift-up step.** Swapping position `a ≥ 2` with its parent of larger key moves the defect +of a sift-up from `a` to its parent. Here `g` is `f` with the values at `a` and `a / 2` +exchanged. -/ +private theorem UpOrdered.swap (h : UpOrdered f S a) (ha : 2 ≤ a) (haS : a ≤ S) + (hlt : f a < f (a / 2)) (hga : g a = f (a / 2)) (hgp : g (a / 2) = f a) + (hgo : ∀ i, 1 ≤ i → i ≤ S → i ≠ a → i ≠ a / 2 → g i = f i) : + UpOrdered g S (a / 2) := by + intro i hi hiS + obtain ⟨h₁, h₂⟩ := h i hi hiS + constructor + · intro hne + rcases eq_or_ne i a with rfl | hia + · rw [hga, hgp] + exact hlt.le + rw [hgo i (by omega) hiS hia hne] + rcases eq_or_ne (i / 2) a with hia2 | hia2 + · rw [hia2, hga] + exact h₂ hia2 ha + rcases eq_or_ne (i / 2) (a / 2) with hip | hip + · have := h₁ hia + rw [hip] at this + rw [hip, hgp] + omega + · rw [hgo (i / 2) (by omega) (by omega) hia2 hip] + exact h₁ hia + · intro hip hp2 + have hpa := (h (a / 2) hp2 (by omega)).1 (by omega) + rw [hgo (a / 2 / 2) (by omega) (by omega) (by omega) (by omega)] + rcases eq_or_ne i a with rfl | hia + · rw [hga] + exact hpa + · have := h₁ hia + rw [hip] at this + rw [hgo i (by omega) hiS hia (by omega)] + omega + +/-- A heap ordered except below `a`, where `a` is at most its children, is ordered. -/ +private theorem DownOrdered.ordered (h : DownOrdered f S a) + (hle : ∀ i, 2 ≤ i → i ≤ S → i / 2 = a → f a ≤ f i) : Ordered f S := by + intro i hi hiS + by_cases hia : i / 2 = a + · rw [hia] + exact hle i hi hiS hia + · exact (h i hi hiS).1 hia + +/-- **A sift-down step.** Swapping position `a` with its child `c` of least key, when that key is +smaller than the key of `a`, moves the defect of a sift-down from `a` to `c`. Here `g` is `f` with +the values at `a` and `c` exchanged. -/ +private theorem DownOrdered.swap (h : DownOrdered f S a) (hc : c / 2 = a) (hcS : c ≤ S) + (hmin : ∀ i, 2 ≤ i → i ≤ S → i / 2 = a → f c ≤ f i) (hlt : f c < f a) + (hga : g a = f c) (hgc : g c = f a) + (hgo : ∀ i, 1 ≤ i → i ≤ S → i ≠ a → i ≠ c → g i = f i) : + DownOrdered g S c := by + intro i hi hiS + obtain ⟨h₁, -⟩ := h i hi hiS + constructor + · intro hic + rcases eq_or_ne i c with rfl | hic' + · rw [hc, hga, hgc] + exact hlt.le + rcases eq_or_ne (i / 2) a with hia | hia + · rw [hia, hga, hgo i (by omega) hiS (by omega) hic'] + exact hmin i hi hiS hia + rcases eq_or_ne i a with rfl | hia' + · rw [hga, hgo (i / 2) (by omega) (by omega) (by omega) hic] + exact (h c (by omega) hcS).2 hc hi + · rw [hgo i (by omega) hiS hia' hic', hgo (i / 2) (by omega) (by omega) hia hic] + exact h₁ hia + · intro hic _ + have := h₁ (by omega) + rw [hic] at this + rw [hc, hga, hgo i (by omega) hiS (by omega) (by omega)] + exact this + +end Order + +/-! ### Swapping two heap entries -/ + +/-- The memory after `swapCode` exchanges the vertices `va` at position `a` and `vb` at position +`b` and updates their positions. -/ +private def swapMem (N Db : ℕ) (mem : ℕ → ℕ) (a b va vb : ℕ) : ℕ → ℕ := + Function.update (Function.update (Function.update (Function.update mem (Db + N + a) vb) + (Db + N + b) va) (Db + 2 * N + 1 + vb) a) (Db + 2 * N + 1 + va) b + +/-- **Swapping two heap entries** keeps the shape of the heap and the set of vertices in it, +exchanges the keys of the two positions, and changes no memory outside the heap and positions +arrays. -/ +private theorem HeapShape.swap {mem : ℕ → ℕ} {a b : ℕ} (h : HeapShape N Db S mem) + (ha : 1 ≤ a) (haS : a ≤ S) (hb : 1 ≤ b) (hbS : b ≤ S) (hab : a ≠ b) : + let m' := swapMem N Db mem a b (mem (Db + N + a)) (mem (Db + N + b)) + HeapShape N Db S m' ∧ heapKey N Db m' a = heapKey N Db mem b ∧ + heapKey N Db m' b = heapKey N Db mem a ∧ + (∀ i, 1 ≤ i → i ≤ S → i ≠ a → i ≠ b → heapKey N Db m' i = heapKey N Db mem i) ∧ + (∀ v < N, m' (Db + 2 * N + 1 + v) = 0 ↔ mem (Db + 2 * N + 1 + v) = 0) ∧ + (∀ c, c ≤ Db + N ∨ Db + 3 * N < c → m' c = mem c) := by + intro m' + have hSN := h.size_le + have hva : mem (Db + N + a) < N := h.vertex_lt a ha haS + have hvb : mem (Db + N + b) < N := h.vertex_lt b hb hbS + have hpa := h.pos_vertex a ha haS + have hpb := h.pos_vertex b hb hbS + have hne : mem (Db + N + a) ≠ mem (Db + N + b) := fun he => hab (by rw [← hpa, ← hpb, he]) + have hH : ∀ i, 1 ≤ i → i ≤ S → m' (Db + N + i) = if i = b then mem (Db + N + a) + else if i = a then mem (Db + N + b) else mem (Db + N + i) := by + intro i hi hiS + simp only [m', swapMem, Function.update_apply] + split_ifs <;> omega + have hP : ∀ v < N, m' (Db + 2 * N + 1 + v) = if v = mem (Db + N + a) then b + else if v = mem (Db + N + b) then a else mem (Db + 2 * N + 1 + v) := by + intro v hv + simp only [m', swapMem, Function.update_apply] + split_ifs <;> omega + have hD : ∀ c, c ≤ Db + N ∨ Db + 3 * N < c → m' c = mem c := by + intro c hc + simp only [m', swapMem, Function.update_apply] + split_ifs <;> omega + have hK : ∀ v < N, m' (Db + v) = mem (Db + v) := fun v hv => hD _ (Or.inl (by omega)) + refine ⟨⟨hSN, ?_, ?_, ?_, ?_, ?_⟩, ?_, ?_, ?_, ?_, hD⟩ + · intro i hi hiS + rw [hH i hi hiS] + split_ifs + exacts [hva, hvb, h.vertex_lt i hi hiS] + · intro i hi hiS + rw [hH i hi hiS] + split_ifs with h1 h2 + · rw [hP _ hva, ite_eq_left rfl, h1] + · rw [hP _ hvb, ite_eq_right hne.symm, ite_eq_left rfl, h2] + · have hpi := h.pos_vertex i hi hiS + have h1' : mem (Db + N + i) ≠ mem (Db + N + a) := fun he => h2 (by rw [← hpi, he, hpa]) + have h2' : mem (Db + N + i) ≠ mem (Db + N + b) := fun he => h1 (by rw [← hpi, he, hpb]) + rw [hP _ (h.vertex_lt i hi hiS), ite_eq_right h1', ite_eq_right h2', hpi] + · intro v hv + rw [hP v hv] + split_ifs + exacts [hbS, haS, h.pos_le v hv] + · intro v hv hv0 + rw [hP v hv] at hv0 ⊢ + split_ifs at hv0 ⊢ with h1 h2 + · rw [hH b hb hbS, ite_eq_left rfl, h1] + · rw [hH a ha haS, ite_eq_right hab, ite_eq_left rfl, h2] + · have hvp := h.vertex_pos v hv hv0 + have hle := h.pos_le v hv + have h1' : mem (Db + 2 * N + 1 + v) ≠ b := fun he => h2 (by rw [← hvp, he]) + have h2' : mem (Db + 2 * N + 1 + v) ≠ a := fun he => h1 (by rw [← hvp, he]) + rw [hH _ (by omega) hle, ite_eq_right h1', ite_eq_right h2', hvp] + · intro v hv hv0 + rw [hK v hv] + rw [hP v hv] at hv0 + split_ifs at hv0 with h1 h2 + · subst h1 + exact h.key_ne _ hva (by omega) + · subst h2 + exact h.key_ne _ hvb (by omega) + · exact h.key_ne v hv hv0 + · simp only [heapKey] + rw [hH a ha haS, ite_eq_right hab, ite_eq_left rfl, hK _ hvb] + · simp only [heapKey] + rw [hH b hb hbS, ite_eq_left rfl, hK _ hva] + · intro i hi hiS hia hib + simp only [heapKey] + rw [hH i hi hiS, ite_eq_right hib, ite_eq_right hia, hK _ (h.vertex_lt i hi hiS)] + · intro v hv + rw [hP v hv] + split_ifs with h1 h2 + · subst h1 + rw [hpa] + omega + · subst h2 + rw [hpb] + omega + · rfl + + +/-! ### Sift-up -/ + +/-- **One sift-up step** from position `K ≥ 2`: if the key at `K` is smaller than the key of its +parent, it swaps the two entries and moves to the parent; otherwise it moves to the root. -/ +private theorem siftUpBody_spec {w K : ℕ} {y : Data 27} (hone : y.regs one = 1) + (htwo : y.regs two = 2) (hdb : y.regs db = Db) (hhb : y.regs hb = Db + N) + (hpb : y.regs pb = Db + 2 * N + 1) (hk : y.regs k = K) (hK : 2 ≤ K) (hKS : K ≤ S) + (hfit : Fits w (Db + 3 * N + 1)) (hshape : HeapShape N Db S y.mem) : + ∃ z, ExecWithin w siftUpBody y z 28 ∧ + (∀ i ∉ [k, p, vk, vp, dk, dp, t], z.regs i = y.regs i) ∧ + (heapKey N Db y.mem K < heapKey N Db y.mem (K / 2) ∧ z.regs k = K / 2 ∧ + z.mem = swapMem N Db y.mem K (K / 2) (y.mem (Db + N + K)) (y.mem (Db + N + K / 2)) ∨ + ¬heapKey N Db y.mem K < heapKey N Db y.mem (K / 2) ∧ z.regs k = 1 ∧ z.mem = y.mem) := by + have hSN := hshape.size_le + obtain ⟨VK, hvk⟩ : ∃ v, y.mem (Db + N + K) = v := ⟨_, rfl⟩ + obtain ⟨VP, hvp⟩ : ∃ v, y.mem (Db + N + K / 2) = v := ⟨_, rfl⟩ + have hVK : VK < N := hvk ▸ hshape.vertex_lt K (by omega) hKS + have hVP : VP < N := hvp ▸ hshape.vertex_lt (K / 2) (by omega) (by omega) + have f₁ : Fits w (K / 2) := hfit.of_le (by omega) + have f₂ : Fits w (Db + N + K) := hfit.of_le (by omega) + have f₃ : Fits w (Db + N + K / 2) := hfit.of_le (by omega) + have f₄ : Fits w (Db + VK) := hfit.of_le (by omega) + have f₅ : Fits w (Db + VP) := hfit.of_le (by omega) + have f₆ : Fits w (Db + 2 * N + 1 + VK) := hfit.of_le (by omega) + have f₇ : Fits w (Db + 2 * N + 1 + VP) := hfit.of_le (by omega) + have hframe : ∀ z T, ExecWithin w siftUpBody y z T → + ∀ i ∉ [k, p, vk, vp, dk, dp, t], z.regs i = y.regs i := by + intro z T he i hi + apply he.regs_eq + simp only [List.mem_cons, List.not_mem_nil, or_false, not_or] at hi + simp [siftUpBody, swapCode, Code.Preserves, Op.dest, Op.exactBin, eq_comm, hi.1, hi.2.1, + hi.2.2.1, hi.2.2.2.1, hi.2.2.2.2.1, hi.2.2.2.2.2.1, hi.2.2.2.2.2.2] + have hkey : ∀ a, heapKey N Db y.mem a = y.mem (Db + y.mem (Db + N + a)) := fun _ => rfl + rw [hkey, hkey, hvk, hvp] + by_cases hlt : y.mem (Db + VK) < y.mem (Db + VP) + · obtain ⟨z, he, hz⟩ : ∃ z, ExecWithin w siftUpBody y z 28 ∧ z.regs k = K / 2 ∧ + z.mem = swapMem N Db y.mem K (K / 2) VK VP := by + refine ⟨?z₁, ?e₁, ?_, ?_⟩ + case e₁ => + ram_exact_eval [siftUpBody, swapCode, hone, htwo, hdb, hhb, hpb, hk, hvk, hvp, f₁, f₂, f₃, + f₄, f₅, f₆, f₇, hlt] + · simp + · simp [swapMem] + exact ⟨z, he, hframe z _ he, Or.inl ⟨hlt, hz⟩⟩ + · obtain ⟨z, he, hz⟩ : ∃ z, ExecWithin w siftUpBody y z 28 ∧ z.regs k = 1 ∧ + z.mem = y.mem := by + refine ⟨?z₂, ?e₂, ?_, ?_⟩ + case e₂ => + ram_exact_eval [siftUpBody, swapCode, hone, htwo, hdb, hhb, hpb, hk, hvk, hvp, f₁, f₂, f₃, + f₄, f₅, f₆, f₇, hlt] + · simp + · simp + exact ⟨z, he, hframe z _ he, Or.inr ⟨hlt, hz⟩⟩ + +/-- The invariant of the sift-up loop started from the data `x`. -/ +private structure UpInv (N Db S : ℕ) (x y : Data 27) : Prop where + /-- Only the registers of the sift-up change. -/ + frame : ∀ i ∉ [k, p, vk, vp, dk, dp, t], y.regs i = x.regs i + /-- The position is in the heap. -/ + one_le : 1 ≤ y.regs k + /-- The position is in the heap. -/ + le_size : y.regs k ≤ S + /-- The heap keeps its shape. -/ + shape : HeapShape N Db S y.mem + /-- The heap is ordered except above the position. -/ + ordered : UpOrdered (heapKey N Db y.mem) S (y.regs k) + /-- The heap keeps its vertices. -/ + pos : ∀ v < N, y.mem (Db + 2 * N + 1 + v) = 0 ↔ x.mem (Db + 2 * N + 1 + v) = 0 + /-- Memory outside the heap and positions arrays is unchanged. -/ + outside : ∀ a, a ≤ Db + N ∨ Db + 3 * N < a → y.mem a = x.mem a + +/-- **Sift-up.** From a heap ordered except between position `K` and its parent, the sift-up +runs within `30 ⌊log₂ K⌋ + 2` instructions, writes only the registers `k`, `p`, `vk`, `vp`, `dk`, +`dp`, `t` and the heap and positions arrays, and leaves a heap with the same vertices. -/ +theorem siftUp_spec {w K : ℕ} {x : Data 27} (hone : x.regs one = 1) (htwo : x.regs two = 2) + (hdb : x.regs db = Db) (hhb : x.regs hb = Db + N) (hpb : x.regs pb = Db + 2 * N + 1) + (hk : x.regs k = K) (hK : 1 ≤ K) (hKS : K ≤ S) (hfit : Fits w (Db + 3 * N + 1)) + (hshape : HeapShape N Db S x.mem) (hord : UpOrdered (heapKey N Db x.mem) S K) : + ∃ y, ExecWithin w siftUp x y (30 * Nat.log 2 K + 2) ∧ + (∀ i ∉ [k, p, vk, vp, dk, dp, t], y.regs i = x.regs i) ∧ + HeapInv N Db S y.mem ∧ + (∀ v < N, y.mem (Db + 2 * N + 1 + v) = 0 ↔ x.mem (Db + 2 * N + 1 + v) = 0) ∧ + (∀ a, a ≤ Db + N ∨ Db + 3 * N < a → y.mem a = x.mem a) := by + have hinit : UpInv N Db S x x := + ⟨fun _ _ => rfl, hk ▸ hK, hk ▸ hKS, hshape, hk ▸ hord, fun _ _ => Iff.rfl, fun _ _ => rfl⟩ + have hbody : ∀ y, UpInv N Db S x y → (Cond.lt one k).Holds y → + ∃ z, ExecWithin w siftUpBody y z 28 ∧ UpInv N Db S x z ∧ + Nat.log 2 (z.regs k) < Nat.log 2 (y.regs k) := by + intro y hy hc + have ho : y.regs one = 1 := (hy.frame one (by decide)).trans hone + have hK2 : 2 ≤ y.regs k := by + simp only [Cond.holds_lt, ho] at hc + omega + have hlog : 0 < Nat.log 2 (y.regs k) := Nat.log_pos (by norm_num) hK2 + have hle := hy.le_size + obtain ⟨z, he, hfr, hz⟩ := siftUpBody_spec ho ((hy.frame two (by decide)).trans htwo) + ((hy.frame db (by decide)).trans hdb) ((hy.frame hb (by decide)).trans hhb) + ((hy.frame pb (by decide)).trans hpb) rfl hK2 hy.le_size hfit hy.shape + have hfr' : ∀ i ∉ [k, p, vk, vp, dk, dp, t], z.regs i = x.regs i := + fun i hi => (hfr i hi).trans (hy.frame i hi) + rcases hz with ⟨hlt, hzk, hzm⟩ | ⟨hlt, hzk, hzm⟩ + · obtain ⟨hs, hga, hgb, hgo, hpos, hout⟩ := HeapShape.swap hy.shape (a := y.regs k) + (b := y.regs k / 2) (by omega) hy.le_size (by omega) (by omega) (by omega) + rw [← hzm] at hs hga hgb hgo hpos hout + refine ⟨z, he, ⟨hfr', by omega, by omega, hs, ?_, ?_, ?_⟩, ?_⟩ + · rw [hzk] + exact hy.ordered.swap hK2 hy.le_size hlt hga hgb hgo + · intro v hv + rw [hpos v hv] + exact hy.pos v hv + · intro a ha + rw [hout a ha] + exact hy.outside a ha + · rw [hzk, Nat.log_div_base] + omega + · refine ⟨z, he, ⟨hfr', by omega, by omega, hzm ▸ hy.shape, ?_, ?_, ?_⟩, ?_⟩ + · rw [hzk, hzm] + exact (hy.ordered.ordered (fun _ _ => not_lt.1 hlt)).upOrdered_one + · intro v hv + rw [hzm] + exact hy.pos v hv + · intro a ha + rw [hzm] + exact hy.outside a ha + · rw [hzk, Nat.log_one_right] + exact hlog + obtain ⟨y, he, hy, hexit⟩ := ExecWithin.loop_invariant (cond := .lt one k) + (body := siftUpBody) (UpInv N Db S x) (fun y => Nat.log 2 (y.regs k)) hbody hinit + have ho : y.regs one = 1 := (hy.frame one (by decide)).trans hone + have hk1 : y.regs k = 1 := by + have := hy.one_le + simp only [Cond.holds_lt, ho] at hexit + omega + refine ⟨y, he.mono (by rw [hk]; omega), hy.frame, ⟨hy.shape, ?_⟩, hy.pos, hy.outside⟩ + have hord' := hy.ordered + rw [hk1] at hord' + exact hord'.ordered (fun h => absurd h (by omega)) + +/-! ### Sift-down -/ + +/-- The child of position `a` of least key, for a position `a` with a child in a heap of `S` +entries with keys `f`: the right child `2 a + 1` if it is in the heap and has the smaller key, +and otherwise the left child `2 a`. -/ +private def minChild (f : ℕ → ℕ) (S a : ℕ) : ℕ := + if 2 * a + 1 ≤ S ∧ f (2 * a + 1) < f (2 * a) then 2 * a + 1 else 2 * a + +/-- The child of least key is a child in the heap, and its key is at most the keys of the +children of `a`. -/ +private theorem minChild_spec {f : ℕ → ℕ} {a : ℕ} (ha : 2 * a ≤ S) : + (minChild f S a = 2 * a ∨ minChild f S a = 2 * a + 1) ∧ minChild f S a ≤ S ∧ + ∀ i, i ≤ S → i / 2 = a → f (minChild f S a) ≤ f i := by + unfold minChild + split_ifs with h + · refine ⟨Or.inr rfl, h.1, fun i hiS hia => ?_⟩ + rcases (by omega : i = 2 * a ∨ i = 2 * a + 1) with rfl | rfl + · exact h.2.le + · exact le_rfl + · refine ⟨Or.inl rfl, ha, fun i hiS hia => ?_⟩ + rcases (by omega : i = 2 * a ∨ i = 2 * a + 1) with rfl | rfl + · exact le_rfl + · exact not_lt.1 fun hlt => h ⟨hiS, hlt⟩ + +/-- `siftDownChoose` leaves the child of least key of position `K` in `p`, its vertex in `vp` and +its key in `dp`. -/ +private theorem siftDownChoose_spec {w K : ℕ} {y : Data 27} (hone : y.regs one = 1) + (hdb : y.regs db = Db) (hhb : y.regs hb = Db + N) (hl : y.regs l = 2 * K) + (hszp : y.regs szp = S + 1) (hK : 1 ≤ K) (hKS : 2 * K ≤ S) (hfit : Fits w (Db + 3 * N + 1)) + (hshape : HeapShape N Db S y.mem) : + ∃ z, ExecWithin w siftDownChoose y z 15 ∧ + (∀ i ∉ [p, t, vp, dp, vk, dk], z.regs i = y.regs i) ∧ z.mem = y.mem ∧ + z.regs p = minChild (heapKey N Db y.mem) S K ∧ + z.regs vp = y.mem (Db + N + minChild (heapKey N Db y.mem) S K) ∧ + z.regs dp = heapKey N Db y.mem (minChild (heapKey N Db y.mem) S K) := by + have hSN := hshape.size_le + obtain ⟨VL, hvl⟩ : ∃ v, y.mem (Db + N + 2 * K) = v := ⟨_, rfl⟩ + have hVL : VL < N := hvl ▸ hshape.vertex_lt (2 * K) (by omega) hKS + have f₁ : Fits w (Db + N + 2 * K) := hfit.of_le (by omega) + have f₂ : Fits w (Db + VL) := hfit.of_le (by omega) + have f₃ : Fits w (2 * K + 1) := hfit.of_le (by omega) + have hframe : ∀ z T, ExecWithin w siftDownChoose y z T → + ∀ i ∉ [p, t, vp, dp, vk, dk], z.regs i = y.regs i := by + intro z T he i hi + apply he.regs_eq + simp only [List.mem_cons, List.not_mem_nil, or_false, not_or] at hi + simp [siftDownChoose, Code.Preserves, Op.dest, Op.exactBin, eq_comm, hi.1, hi.2.1, + hi.2.2.1, hi.2.2.2.1, hi.2.2.2.2.1, hi.2.2.2.2.2] + have hkey : ∀ a, heapKey N Db y.mem a = y.mem (Db + y.mem (Db + N + a)) := fun _ => rfl + by_cases h₁ : 2 * K + 1 < S + 1 + · obtain ⟨VR, hvr⟩ : ∃ v, y.mem (Db + N + (2 * K + 1)) = v := ⟨_, rfl⟩ + have hVR : VR < N := hvr ▸ hshape.vertex_lt (2 * K + 1) (by omega) (by omega) + have f₄ : Fits w (Db + N + (2 * K + 1)) := hfit.of_le (by omega) + have f₅ : Fits w (Db + VR) := hfit.of_le (by omega) + by_cases h₂ : y.mem (Db + VR) < y.mem (Db + VL) + · have hC : minChild (heapKey N Db y.mem) S K = 2 * K + 1 := by + unfold minChild + exact ite_eq_left ⟨by omega, by rw [hkey, hkey, hvl, hvr]; exact h₂⟩ + rw [hC, hkey, hvr] + obtain ⟨z, he, hz⟩ : ∃ z, ExecWithin w siftDownChoose y z 15 ∧ z.mem = y.mem ∧ + z.regs p = 2 * K + 1 ∧ z.regs vp = VR ∧ z.regs dp = y.mem (Db + VR) := by + refine ⟨?z₁, ?e₁, ?_, ?_, ?_, ?_⟩ + case e₁ => + ram_exact_eval [siftDownChoose, hone, hdb, hhb, hl, hszp, hvl, hvr, f₁, f₂, f₃, f₄, f₅, + h₁, h₂] + all_goals simp +decide + exact ⟨z, he, hframe z _ he, hz⟩ + · have hC : minChild (heapKey N Db y.mem) S K = 2 * K := by + unfold minChild + exact ite_eq_right fun h => h₂ (by rw [hkey, hkey, hvl, hvr] at h; exact h.2) + rw [hC, hkey, hvl] + obtain ⟨z, he, hz⟩ : ∃ z, ExecWithin w siftDownChoose y z 15 ∧ z.mem = y.mem ∧ + z.regs p = 2 * K ∧ z.regs vp = VL ∧ z.regs dp = y.mem (Db + VL) := by + refine ⟨?z₂, ?e₂, ?_, ?_, ?_, ?_⟩ + case e₂ => + ram_exact_eval [siftDownChoose, hone, hdb, hhb, hl, hszp, hvl, hvr, f₁, f₂, f₃, f₄, f₅, + h₁, h₂] + all_goals simp +decide + exact ⟨z, he, hframe z _ he, hz⟩ + · have hC : minChild (heapKey N Db y.mem) S K = 2 * K := by + unfold minChild + exact ite_eq_right fun h => h₁ (by omega) + rw [hC, hkey, hvl] + obtain ⟨z, he, hz⟩ : ∃ z, ExecWithin w siftDownChoose y z 15 ∧ z.mem = y.mem ∧ + z.regs p = 2 * K ∧ z.regs vp = VL ∧ z.regs dp = y.mem (Db + VL) := by + refine ⟨?z₃, ?e₃, ?_, ?_, ?_, ?_⟩ + case e₃ => + ram_exact_eval [siftDownChoose, hone, hdb, hhb, hl, hszp, hvl, f₁, f₂, f₃, h₁] + all_goals simp +decide + exact ⟨z, he, hframe z _ he, hz⟩ + +/-- `siftDownSink` swaps position `K` with the position `C` in `p` if the key in `dp` is smaller +than the key at `K`, and moves to `C`; otherwise it stops the sift-down. -/ +private theorem siftDownSink_spec {w K C : ℕ} {y : Data 27} (hdb : y.regs db = Db) + (hhb : y.regs hb = Db + N) (hpb : y.regs pb = Db + 2 * N + 1) (hk : y.regs k = K) + (hp : y.regs p = C) (hvp : y.regs vp = y.mem (Db + N + C)) (hszp : y.regs szp = S + 1) + (hK : 1 ≤ K) (hKS : K ≤ S) (hC : 1 ≤ C) (hCS : C ≤ S) (hfit : Fits w (Db + 3 * N + 1)) + (hshape : HeapShape N Db S y.mem) : + ∃ z, ExecWithin w siftDownSink y z 15 ∧ (∀ i ∉ [t, vk, dk, k, l], z.regs i = y.regs i) ∧ + (y.regs dp < heapKey N Db y.mem K ∧ z.regs k = C ∧ z.regs l = 2 * C ∧ + z.mem = swapMem N Db y.mem K C (y.mem (Db + N + K)) (y.mem (Db + N + C)) ∨ + ¬y.regs dp < heapKey N Db y.mem K ∧ z.regs k = K ∧ z.regs l = S + 1 ∧ + z.mem = y.mem) := by + have hSN := hshape.size_le + obtain ⟨VK, hvk⟩ : ∃ v, y.mem (Db + N + K) = v := ⟨_, rfl⟩ + obtain ⟨VC, hvc⟩ : ∃ v, y.mem (Db + N + C) = v := ⟨_, rfl⟩ + rw [hvc] at hvp + have hVK : VK < N := hvk ▸ hshape.vertex_lt K hK hKS + have hVC : VC < N := hvc ▸ hshape.vertex_lt C hC hCS + have f₁ : Fits w (Db + N + K) := hfit.of_le (by omega) + have f₂ : Fits w (Db + VK) := hfit.of_le (by omega) + have f₃ : Fits w (Db + N + C) := hfit.of_le (by omega) + have f₄ : Fits w (Db + 2 * N + 1 + VK) := hfit.of_le (by omega) + have f₅ : Fits w (Db + 2 * N + 1 + VC) := hfit.of_le (by omega) + have f₆ : Fits w (C + C) := hfit.of_le (by omega) + have hframe : ∀ z T, ExecWithin w siftDownSink y z T → + ∀ i ∉ [t, vk, dk, k, l], z.regs i = y.regs i := by + intro z T he i hi + apply he.regs_eq + simp only [List.mem_cons, List.not_mem_nil, or_false, not_or] at hi + simp [siftDownSink, swapCode, Code.Preserves, Op.dest, Op.exactBin, eq_comm, hi.1, hi.2.1, + hi.2.2.1, hi.2.2.2.1, hi.2.2.2.2] + have hkey : heapKey N Db y.mem K = y.mem (Db + VK) := by simp only [heapKey, hvk] + rw [hkey, hvk, hvc] + by_cases hlt : y.regs dp < y.mem (Db + VK) + · obtain ⟨z, he, hz⟩ : ∃ z, ExecWithin w siftDownSink y z 15 ∧ z.regs k = C ∧ + z.regs l = 2 * C ∧ z.mem = swapMem N Db y.mem K C VK VC := by + refine ⟨?z₁, ?e₁, ?_, ?_, ?_⟩ + case e₁ => + ram_exact_eval [siftDownSink, swapCode, hdb, hhb, hpb, hk, hp, hvp, hszp, hvk, f₁, f₂, f₃, + f₄, f₅, f₆, hlt] + · simp +decide + · simp only [Data.set_regs, Data.write_regs, Function.update_idem, Function.update_self] + omega + · simp [swapMem] + exact ⟨z, he, hframe z _ he, Or.inl ⟨hlt, hz⟩⟩ + · obtain ⟨z, he, hz⟩ : ∃ z, ExecWithin w siftDownSink y z 15 ∧ z.regs k = K ∧ + z.regs l = S + 1 ∧ z.mem = y.mem := by + refine ⟨?z₂, ?e₂, ?_, ?_, ?_⟩ + case e₂ => + ram_exact_eval [siftDownSink, swapCode, hdb, hhb, hpb, hk, hp, hvp, hszp, hvk, f₁, f₂, hlt] + · simp +decide [hk] + · simp +decide + · simp + exact ⟨z, he, hframe z _ he, Or.inr ⟨hlt, hz⟩⟩ + +/-- One round of the sift-down at position `K` with a child in the heap: with `C` the child of +least key, it swaps `K` with `C` and moves to `C` if `C` has the smaller key, and otherwise it +stops the sift-down. -/ +private theorem siftDownBody_spec {w K C : ℕ} {y : Data 27} (hone : y.regs one = 1) + (hdb : y.regs db = Db) (hhb : y.regs hb = Db + N) (hpb : y.regs pb = Db + 2 * N + 1) + (hk : y.regs k = K) (hl : y.regs l = 2 * K) (hszp : y.regs szp = S + 1) (hK : 1 ≤ K) + (hKS : 2 * K ≤ S) (hfit : Fits w (Db + 3 * N + 1)) (hshape : HeapShape N Db S y.mem) + (hC : C = minChild (heapKey N Db y.mem) S K) : + ∃ z, ExecWithin w siftDownBody y z 30 ∧ + (∀ i ∉ [k, p, vk, vp, dk, dp, t, l], z.regs i = y.regs i) ∧ + (heapKey N Db y.mem C < heapKey N Db y.mem K ∧ z.regs k = C ∧ z.regs l = 2 * C ∧ + z.mem = swapMem N Db y.mem K C (y.mem (Db + N + K)) (y.mem (Db + N + C)) ∨ + ¬heapKey N Db y.mem C < heapKey N Db y.mem K ∧ z.regs k = K ∧ z.regs l = S + 1 ∧ + z.mem = y.mem) := by + obtain ⟨z₁, he₁, hfr₁, hm₁, hp₁, hvp₁, hdp₁⟩ := + siftDownChoose_spec hone hdb hhb hl hszp hK hKS hfit hshape + rw [← hC] at hp₁ hvp₁ hdp₁ + obtain ⟨hCc, hCS, -⟩ := minChild_spec (f := heapKey N Db y.mem) hKS + rw [← hC] at hCc hCS + have hr : ∀ i ∉ [p, t, vp, dp, vk, dk], z₁.regs i = y.regs i := hfr₁ + obtain ⟨z, he₂, hfr₂, hz⟩ := siftDownSink_spec (y := z₁) (K := K) (C := C) + ((hr db (by decide)).trans hdb) ((hr hb (by decide)).trans hhb) + ((hr pb (by decide)).trans hpb) ((hr k (by decide)).trans hk) hp₁ (hm₁ ▸ hvp₁) + ((hr szp (by decide)).trans hszp) hK (by omega) (by omega) hCS hfit (hm₁ ▸ hshape) + rw [hdp₁, hm₁] at hz + refine ⟨z, ?_, fun i hi => ?_, hz⟩ + · exact (he₁.seq he₂).mono (by norm_num) + · refine (hfr₂ i fun h => hi ?_).trans (hr i fun h => hi ?_) <;> + · simp only [List.mem_cons, List.not_mem_nil, or_false] at h ⊢ + tauto + +/-- The invariant of the sift-down loop started from the data `x`. -/ +private structure DownInv (N Db S : ℕ) (x y : Data 27) : Prop where + /-- Only the registers of the sift-down change. -/ + frame : ∀ i ∉ [k, p, vk, vp, dk, dp, t, l, szp], y.regs i = x.regs i + /-- The register `szp` holds `S + 1`. -/ + size_reg : y.regs szp = S + 1 + /-- The position is in the heap. -/ + one_le : 1 ≤ y.regs k + /-- The position is in the heap. -/ + le_size : y.regs k ≤ S + /-- The heap keeps its shape. -/ + shape : HeapShape N Db S y.mem + /-- Either the sift-down goes on at `k`, with `l = 2 k`, and the heap is ordered except below + `k`, or it has stopped, with `l = S + 1`, and the heap is ordered. -/ + order : y.regs l = 2 * y.regs k ∧ DownOrdered (heapKey N Db y.mem) S (y.regs k) ∨ + y.regs l = S + 1 ∧ Ordered (heapKey N Db y.mem) S + /-- The heap keeps its vertices. -/ + pos : ∀ v < N, y.mem (Db + 2 * N + 1 + v) = 0 ↔ x.mem (Db + 2 * N + 1 + v) = 0 + /-- Memory outside the heap and positions arrays is unchanged. -/ + outside : ∀ a, a ≤ Db + N ∨ Db + 3 * N < a → y.mem a = x.mem a + +/-- The set-up of the sift-down loop. -/ +private def downSetup : Code .exact 27 := ram_exact[ l ← k + k; szp ← sz + one ] + +/-- `siftDown` is `downSetup` followed by the loop. -/ +private theorem siftDown_eq : + siftDown = appendCode downSetup (.loop (.lt l szp) siftDownBody) := rfl + +/-- The measure of the sift-down loop: the number of rounds left is at most this. -/ +private def downMeasure (S : ℕ) (y : Data 27) : ℕ := + if y.regs l < S + 1 then Nat.log 2 (S + 1) + 1 - Nat.log 2 (y.regs l) else 0 + +/-- **Sift-down.** From a heap of `S` entries ordered except between position `K` and its +children, the sift-down runs within `45 ⌊log₂ (S + 1)⌋ + 4` instructions, writes only the +registers `k`, `p`, `vk`, `vp`, `dk`, `dp`, `t`, `l`, `szp` and the heap and positions arrays, and +leaves a heap with the same vertices. -/ +theorem siftDown_spec {w K : ℕ} {x : Data 27} (hone : x.regs one = 1) + (hdb : x.regs db = Db) (hhb : x.regs hb = Db + N) (hpb : x.regs pb = Db + 2 * N + 1) + (hsz : x.regs sz = S) (hk : x.regs k = K) (hK : 1 ≤ K) (hKS : K ≤ S) + (hfit : Fits w (Db + 3 * N + 1)) + (hshape : HeapShape N Db S x.mem) (hord : DownOrdered (heapKey N Db x.mem) S K) : + ∃ y, ExecWithin w siftDown x y (45 * Nat.log 2 (S + 1) + 4) ∧ + (∀ i ∉ [k, p, vk, vp, dk, dp, t, l, szp], y.regs i = x.regs i) ∧ + HeapInv N Db S y.mem ∧ + (∀ v < N, y.mem (Db + 2 * N + 1 + v) = 0 ↔ x.mem (Db + 2 * N + 1 + v) = 0) ∧ + (∀ a, a ≤ Db + N ∨ Db + 3 * N < a → y.mem a = x.mem a) := by + have hSN := hshape.size_le + have f₁ : Fits w (K + K) := hfit.of_le (by omega) + have f₂ : Fits w (S + 1) := hfit.of_le (by omega) + obtain ⟨x₁, he₁, hx₁l, hx₁s, hx₁m, hx₁r⟩ : ∃ x₁, ExecWithin w downSetup x x₁ 2 ∧ + x₁.regs l = 2 * K ∧ x₁.regs szp = S + 1 ∧ x₁.mem = x.mem ∧ + ∀ i, i ≠ l → i ≠ szp → x₁.regs i = x.regs i := by + refine ⟨?z₀, ?e₀, ?_, ?_, ?_, ?_⟩ + case e₀ => ram_exact_eval [downSetup, hone, hk, hsz, f₁, f₂] + · simp +decide only [Data.set_regs, ne_eq, Function.update_of_ne, Function.update_self] + omega + · simp + · rfl + · intro i h₁ h₂ + simp [h₁, h₂] + have hx₁k : x₁.regs k = K := (hx₁r k (by decide) (by decide)).trans hk + have hinit : DownInv N Db S x x₁ := by + refine ⟨fun i hi => hx₁r i ?_ ?_, hx₁s, hx₁k ▸ hK, hx₁k ▸ hKS, hx₁m ▸ hshape, + Or.inl ⟨by rw [hx₁l, hx₁k], by rw [hx₁m, hx₁k]; exact hord⟩, fun _ _ => by rw [hx₁m], + fun _ _ => by rw [hx₁m]⟩ + · rintro rfl + simp at hi + · rintro rfl + simp at hi + have hbody : ∀ y, DownInv N Db S x y → (Cond.lt l szp).Holds y → + ∃ z, ExecWithin w siftDownBody y z 30 ∧ DownInv N Db S x z ∧ + downMeasure S z < downMeasure S y := by + intro y hy hc + have hr : ∀ i ∉ [k, p, vk, vp, dk, dp, t, l, szp], y.regs i = x.regs i := hy.frame + have hsz' := hy.size_reg + have hk1 := hy.one_le + have hkS := hy.le_size + simp only [Cond.holds_lt, hsz'] at hc + obtain ⟨hl, hdo⟩ : y.regs l = 2 * y.regs k ∧ DownOrdered (heapKey N Db y.mem) S (y.regs k) := + hy.order.resolve_right fun h => by omega + obtain ⟨C, hC⟩ : ∃ C, C = minChild (heapKey N Db y.mem) S (y.regs k) := ⟨_, rfl⟩ + obtain ⟨hCc, hCS, hmin⟩ := minChild_spec (S := S) (f := heapKey N Db y.mem) + (a := y.regs k) (by omega) + rw [← hC] at hCc hCS hmin + obtain ⟨z, he, hfr, hz⟩ := siftDownBody_spec (hr one (by decide) ▸ hone) + (hr db (by decide) ▸ hdb) (hr hb (by decide) ▸ hhb) (hr pb (by decide) ▸ hpb) rfl hl hsz' + hk1 (by omega) hfit hy.shape hC + have hfr' : ∀ i ∉ [k, p, vk, vp, dk, dp, t, l, szp], z.regs i = x.regs i := by + intro i hi + refine (hfr i fun h => hi ?_).trans (hr i hi) + simp only [List.mem_cons, List.not_mem_nil, or_false] at h ⊢ + tauto + have hzs : z.regs szp = S + 1 := (hfr szp (by decide)).trans hsz' + have hlog : Nat.log 2 (y.regs l) ≤ Nat.log 2 (S + 1) := Nat.log_mono_right (by omega) + have hμ : downMeasure S y = Nat.log 2 (S + 1) + 1 - Nat.log 2 (y.regs l) := by + simp only [downMeasure] + exact ite_eq_left hc + rcases hz with ⟨hlt, hzk, hzl, hzm⟩ | ⟨hlt, hzk, hzl, hzm⟩ + · obtain ⟨hs, hga, hgb, hgo, hpos, hout⟩ := HeapShape.swap hy.shape (a := y.regs k) + (b := C) hk1 hkS (by omega) hCS (by omega) + rw [← hzm] at hs hga hgb hgo hpos hout + refine ⟨z, he, ⟨hfr', hzs, by omega, by omega, hs, Or.inl ⟨by rw [hzl, hzk], ?_⟩, ?_, ?_⟩, + ?_⟩ + · rw [hzk] + exact hdo.swap (by omega) hCS (fun i _ hiS hia => hmin i hiS hia) hlt hga hgb hgo + · intro v hv + rw [hpos v hv] + exact hy.pos v hv + · intro a ha + rw [hout a ha] + exact hy.outside a ha + · have h₁ : Nat.log 2 (y.regs l) ≤ Nat.log 2 C := Nat.log_mono_right (by omega) + have h₂ : Nat.log 2 (z.regs l) = Nat.log 2 C + 1 := by + rw [hzl, mul_comm] + exact Nat.log_mul_base (by norm_num) (by omega) + rw [hμ] + simp only [downMeasure] + split_ifs with h₃ + · have : Nat.log 2 (z.regs l) ≤ Nat.log 2 (S + 1) := Nat.log_mono_right (by omega) + omega + · omega + · refine ⟨z, he, ⟨hfr', hzs, by omega, by omega, hzm ▸ hy.shape, Or.inr ⟨hzl, ?_⟩, ?_, ?_⟩, + ?_⟩ + · rw [hzm] + exact hdo.ordered fun i _ hiS hia => (not_lt.1 hlt).trans (hmin i hiS hia) + · intro v hv + rw [hzm] + exact hy.pos v hv + · intro a ha + rw [hzm] + exact hy.outside a ha + · rw [hμ] + simp only [downMeasure, hzl, lt_irrefl, ite_false] + omega + obtain ⟨y, he, hy, hexit⟩ := ExecWithin.loop_invariant (cond := .lt l szp) + (body := siftDownBody) (DownInv N Db S x) (downMeasure S) hbody hinit + have hμ : downMeasure S x₁ ≤ Nat.log 2 (S + 1) := by + simp only [downMeasure, hx₁l] + split_ifs with h + · have : Nat.log 2 (2 * K) = Nat.log 2 K + 1 := by + rw [mul_comm] + exact Nat.log_mul_base (by norm_num) (by omega) + omega + · omega + refine ⟨y, ?_, hy.frame, ⟨hy.shape, ?_⟩, hy.pos, hy.outside⟩ + · rw [siftDown_eq] + refine (execWithin_append he₁ he).mono ?_ + have := Nat.mul_le_mul_right 32 hμ + omega + · have := hy.size_reg + simp only [Cond.holds_lt, this] at hexit + rcases hy.order with ⟨hl, hdo⟩ | ⟨-, hord'⟩ + · exact hdo.ordered fun i _ hiS hia => absurd hiS (by omega) + · exact hord' + +/-! ### Extract-min -/ + +/-- The part of `extractMin` before its test: it reads the root and removes it from the +positions array, reads the last entry and shrinks the heap. -/ +private def extractHead : Code .exact 27 := ram_exact[ + t ← hb + one; u ← *t; t ← db + u; best ← *t; + t ← pb + u; *t ← zero; + t ← hb + sz; vk ← *t; sz ← sz - one ] + +/-- The part of `extractMin` that moves the last entry to the root, before the sift-down. -/ +private def extractMove : Code .exact 27 := ram_exact[ + t ← hb + one; *t ← vk; t ← pb + vk; *t ← one; k ← one ] + +/-- The part of `extractMin` run when the heap becomes empty. -/ +private def extractStop : Code .exact 27 := ram_exact[ k ← k ] + +/-- `extractMin` is `extractHead` followed by the test on the new size. -/ +private theorem extractMin_eq : extractMin = appendCode extractHead + (.branch (.lt zero sz) (appendCode extractMove siftDown) extractStop) := rfl + +/-- **Extract-min.** From a heap of `S ≥ 1` entries, the removal of the root runs within +`45 ⌊log₂ S⌋ + 20` instructions. It leaves the root vertex in `u` and its key in `best`, and a +heap of `S - 1` entries that holds the vertices of the old heap but the root. -/ +theorem extractMin_spec {w : ℕ} {x : Data 27} (hzero : x.regs zero = 0) (hone : x.regs one = 1) + (hdb : x.regs db = Db) (hhb : x.regs hb = Db + N) (hpb : x.regs pb = Db + 2 * N + 1) + (hsz : x.regs sz = S) (hS : 1 ≤ S) (hfit : Fits w (Db + 3 * N + 1)) + (hheap : HeapInv N Db S x.mem) : + ∃ y, ExecWithin w extractMin x y (45 * Nat.log 2 S + 20) ∧ + (∀ i ∉ [t, u, best, vk, sz, k, p, vp, dk, dp, l, szp], y.regs i = x.regs i) ∧ + y.regs u = x.mem (Db + N + 1) ∧ y.regs best = x.mem (Db + x.mem (Db + N + 1)) ∧ + y.regs sz = S - 1 ∧ HeapInv N Db (S - 1) y.mem ∧ + (∀ v < N, y.mem (Db + 2 * N + 1 + v) = 0 ↔ + x.mem (Db + 2 * N + 1 + v) = 0 ∨ v = x.mem (Db + N + 1)) ∧ + (∀ a, a ≤ Db + N ∨ Db + 3 * N < a → y.mem a = x.mem a) := by + have hSN := hheap.size_le + obtain ⟨U, hU⟩ : ∃ v, x.mem (Db + N + 1) = v := ⟨_, rfl⟩ + obtain ⟨L, hL⟩ : ∃ v, x.mem (Db + N + S) = v := ⟨_, rfl⟩ + have hUN : U < N := hU ▸ hheap.vertex_lt 1 le_rfl hS + have hLN : L < N := hL ▸ hheap.vertex_lt S hS le_rfl + rw [hU] + have f₁ : Fits w (Db + N + 1) := hfit.of_le (by omega) + have f₂ : Fits w (Db + U) := hfit.of_le (by omega) + have f₃ : Fits w (Db + 2 * N + 1 + U) := hfit.of_le (by omega) + have f₄ : Fits w (Db + N + S) := hfit.of_le (by omega) + have f₅ : Fits w (S - 1) := hfit.of_le (by omega) + have hne : ¬Db + N + S = Db + 2 * N + 1 + U := by omega + obtain ⟨x₁, he₁, hx₁u, hx₁b, hx₁v, hx₁s, hx₁m, hx₁r⟩ : ∃ x₁, ExecWithin w extractHead x x₁ 9 ∧ + x₁.regs u = U ∧ x₁.regs best = x.mem (Db + U) ∧ x₁.regs vk = L ∧ x₁.regs sz = S - 1 ∧ + x₁.mem = Function.update x.mem (Db + 2 * N + 1 + U) 0 ∧ + ∀ i ∉ [t, u, best, vk, sz], x₁.regs i = x.regs i := by + refine ⟨?z₀, ?e₀, ?_, ?_, ?_, ?_, ?_, ?_⟩ + case e₀ => + ram_exact_eval [extractHead, hzero, hone, hdb, hhb, hpb, hsz, hU, hL, hne, hS, f₁, f₂, f₃, + f₄, f₅] + · simp +decide + · simp +decide + · simp +decide + · simp +decide + · simp + · intro i hi + simp only [List.mem_cons, List.not_mem_nil, or_false, not_or] at hi + simp [hi.1, hi.2.1, hi.2.2.1, hi.2.2.2.1, hi.2.2.2.2] + have hr₁ : ∀ i ∉ [t, u, best, vk, sz], x₁.regs i = x.regs i := hx₁r + have hU₀ : Function.update x.mem (Db + 2 * N + 1 + U) 0 (Db + N + 1) = U := by + rw [Function.update_of_ne (by omega), hU] + have hout₁ : ∀ a, a ≤ Db + N ∨ Db + 3 * N < a → x₁.mem a = x.mem a := by + intro a ha + rw [hx₁m, Function.update_of_ne (by omega)] + have hframe : ∀ (L₁ L₂ : List (Fin 27)), (∀ i ∈ L₁, i ∈ L₂) → ∀ i ∉ L₂, i ∉ L₁ := + fun _ _ h i hi hi' => hi (h i hi') + by_cases hS2 : 2 ≤ S + · obtain ⟨hs, hdo, hpos⟩ := hheap.extract hS2 + rw [hU, hL] at hs hdo hpos + have f₆ : Fits w (Db + 2 * N + 1 + L) := hfit.of_le (by omega) + obtain ⟨x₂, he₂, hx₂k, hx₂m, hx₂r⟩ : ∃ x₂, ExecWithin w extractMove x₁ x₂ 5 ∧ + x₂.regs k = 1 ∧ x₂.mem = Function.update (Function.update + (Function.update x.mem (Db + 2 * N + 1 + U) 0) (Db + N + 1) L) (Db + 2 * N + 1 + L) 1 ∧ + ∀ i ∉ [t, k], x₂.regs i = x₁.regs i := by + refine ⟨?z₁, ?e₁, ?_, ?_, ?_⟩ + case e₁ => + ram_exact_eval [extractMove, hr₁ one (by decide) ▸ hone, hr₁ hb (by decide) ▸ hhb, + hr₁ pb (by decide) ▸ hpb, hx₁v, hx₁m, f₁, f₆] + · simp +decide + · simp [hx₁m] + · intro i hi + simp only [List.mem_cons, List.not_mem_nil, or_false, not_or] at hi + simp [hi.1, hi.2] + have hr₂ : ∀ i ∉ [t, k], x₂.regs i = x₁.regs i := hx₂r + have hr₂₁ : ∀ i ∉ [t, u, best, vk, sz, k], x₂.regs i = x.regs i := fun i hi => + (hr₂ i (hframe _ _ (by simp) i hi)).trans (hr₁ i (hframe _ _ (by simp) i hi)) + obtain ⟨y, he₃, hfr₃, hheap', hpos₃, hout₃⟩ := siftDown_spec (S := S - 1) (K := 1) + (hr₂₁ one (by decide) ▸ hone) (hr₂₁ db (by decide) ▸ hdb) (hr₂₁ hb (by decide) ▸ hhb) + (hr₂₁ pb (by decide) ▸ hpb) ((hr₂ sz (by decide)).trans hx₁s) hx₂k le_rfl (by omega) + hfit (hx₂m ▸ hs) (hx₂m ▸ hdo) + rw [Nat.sub_add_cancel hS] at he₃ + have hr₃ : ∀ i ∉ [k, p, vk, vp, dk, dp, t, l, szp], y.regs i = x₂.regs i := hfr₃ + have hc : (Cond.lt zero sz).Holds x₁ := by + simp only [Cond.holds_lt, hr₁ zero (by decide) ▸ hzero, hx₁s] + omega + refine ⟨y, ?_, fun i hi => ?_, ?_, ?_, ?_, hheap', fun v hv => ?_, fun a ha => ?_⟩ + · rw [extractMin_eq] + refine (execWithin_append he₁ (ExecWithin.branch_pos hc + (execWithin_append he₂ he₃))).mono ?_ + omega + · exact (hr₃ i (hframe _ _ (by simp) i hi)).trans (hr₂₁ i (hframe _ _ (by simp) i hi)) + · rw [hr₃ u (by decide), hr₂ u (by decide), hx₁u] + · rw [hr₃ best (by decide), hr₂ best (by decide), hx₁b] + · rw [hr₃ sz (by decide), hr₂ sz (by decide), hx₁s] + · rw [hpos₃ v hv, hx₂m] + exact hpos v hv + · rw [hout₃ a ha, hx₂m, Function.update_of_ne (by omega), Function.update_of_ne (by omega), + Function.update_of_ne (by omega)] + · obtain rfl : S = 1 := by omega + have hc : ¬(Cond.lt zero sz).Holds x₁ := by + simp only [Cond.holds_lt, hr₁ zero (by decide) ▸ hzero, hx₁s] + omega + obtain ⟨x₂, he₂, hx₂⟩ : ∃ x₂, ExecWithin w extractStop x₁ x₂ 1 ∧ x₂ = x₁ := by + refine ⟨?z₂, ?e₂, ?_⟩ + case e₂ => ram_exact_eval [extractStop] + simp [Data.set] + subst hx₂ + have hheap' := hheap.extract_one + rw [hU] at hheap' + refine ⟨x₂, ?_, fun i hi => hr₁ i (hframe _ _ (by simp) i hi), hx₁u, hx₁b, hx₁s, + hx₁m ▸ hheap', fun v hv => ?_, hout₁⟩ + · rw [extractMin_eq] + exact (execWithin_append he₁ (ExecWithin.branch_neg hc he₂)).mono (by omega) + · rw [hx₁m, Function.update_apply] + split_ifs with h + · exact iff_of_true rfl (Or.inr (by omega)) + · exact ⟨Or.inl, fun h' => h'.resolve_right (by omega)⟩ + +end Cslib.Algorithms.WordRAM.Dijkstra + +end diff --git a/Complexity/Algorithms/WordRAM/Dijkstra/HeapStatement.lean b/Complexity/Algorithms/WordRAM/Dijkstra/HeapStatement.lean new file mode 100644 index 0000000000..d0cf5c57c7 --- /dev/null +++ b/Complexity/Algorithms/WordRAM/Dijkstra/HeapStatement.lean @@ -0,0 +1,169 @@ +/- +Copyright (c) 2026 Renfei Zhou. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: Renfei Zhou +-/ + +module +-- upstream: candidate + +public import Complexity.Algorithms.WordRAM.Dijkstra.Statement +public import Complexity.Computability.WordRAM.Time + +/-! +# Dijkstra's algorithm on the word RAM with adjacency lists: statements + +This file states, without proof, what a single-source shortest-path program on the word RAM +must satisfy to run in `O((m + 1) log n)` time on a graph with `n` vertices and `m` edges, the +bound of Dijkstra's algorithm with a binary heap [Dijkstra][Dijkstra1959] +[Cormen et al., Section 24.3][CLRS2009]. The bound has no `n` term: memory past the input starts +at `0`, so the per-vertex arrays need no initialization, and the offsets give the adjacency list +of any vertex directly, so the program touches only the vertices it reaches. + +An instance (`Cslib.Algorithms.WordRAM.Dijkstra.ListInstance`) is a directed graph on `Fin n` +given by its adjacency lists of `(target, weight)` pairs, with natural-number weights, and a +source vertex. The lists may repeat a pair of endpoints (parallel edges) and may contain +self-loops. The instance determines a weight-matrix instance of +`Cslib.Algorithms.WordRAM.Dijkstra.Instance` that keeps, for each pair of endpoints, the +lightest of its edges (`Cslib.Algorithms.WordRAM.Dijkstra.ListInstance.toInstance`); the distance +is the distance of that instance. + +The program reads the input words `n`, the source, the `n + 1` offsets of the adjacency lists in +the edge array, and the edges in order, each as two words, its target and its weight. When it +halts, the `n` memory cells just past the input hold the distances from the source, with `0` for +an unreachable vertex and `d + 1` for distance `d`. + +## Main definitions + +* `Cslib.Algorithms.WordRAM.Dijkstra.ListInstance`: a directed graph by adjacency lists, with a + source. +* `Cslib.Algorithms.WordRAM.Dijkstra.ListInstance.toInstance`: its weight-matrix instance. +* `Cslib.Algorithms.WordRAM.Dijkstra.ListInstance.encode`: its input words. +* `Cslib.Algorithms.WordRAM.Dijkstra.LogLinearSpec`: a flat program computes every distance from + the source within `K (m + 1) (⌊log₂ n⌋ + 1)` instructions. +* `Cslib.Algorithms.WordRAM.Dijkstra.LogLinearStatement`: some program meets `LogLinearSpec`. + +## Main results + +* `Cslib.Algorithms.WordRAM.Dijkstra.ListInstance.lightest_eq_some_iff`: the matrix entry of a + pair of endpoints is the least weight of the edges between them. +* `Cslib.Algorithms.WordRAM.Dijkstra.ListInstance.lightest_eq_none_iff`: it is absent exactly + when there is no such edge. + +## References + +* [E. W. Dijkstra, *A note on two problems in connexion with graphs*][Dijkstra1959] +* [T. H. Cormen, C. E. Leiserson, R. L. Rivest and C. Stein, *Introduction to + Algorithms*][CLRS2009] +-/ + +@[expose] public section + +namespace Cslib.Algorithms.WordRAM.Dijkstra + +open Cslib.WordRAM + +/-- An instance of the **single-source shortest-path problem** given by adjacency lists: a +directed graph on `Fin n` and a source vertex. The list `adj u` holds a pair `(v, c)` for each +edge from `u` to `v` of weight `c`; it may hold several pairs with the same `v`. -/ +structure ListInstance where + /-- The number of vertices. -/ + n : ℕ + /-- The adjacency lists. -/ + adj : Fin n → List (Fin n × ℕ) + /-- The source. -/ + source : Fin n + +namespace ListInstance + +variable (G : ListInstance) + +/-- The least weight of an edge from `u` to `v`, and `none` when there is no such edge. -/ +def lightest (u v : Fin G.n) : Option ℕ := + ((G.adj u).filterMap fun e => if e.1 = v then some e.2 else none).min? + +/-- The weight-matrix instance with the same vertices and source, whose entry for `u` and `v` is +the least weight of an edge from `u` to `v`. Keeping a lightest edge between each pair preserves +distances. -/ +def toInstance : Instance where + n := G.n + weight := G.lightest + source := G.source + +/-- The **shortest-path distance** from `u` to `v`: the distance of the weight-matrix instance, +`Cslib.Algorithms.WordRAM.Dijkstra.Instance.dist`. -/ +noncomputable def dist (u v : Fin G.n) : ℕ∞ := + G.toInstance.dist u v + +/-- The number of edges. -/ +def m : ℕ := ∑ u, (G.adj u).length + +/-- The largest edge weight, and `0` when there is no edge. -/ +def maxWeight : ℕ := + Finset.univ.sup fun u => ((G.adj u).map Prod.snd).foldr max 0 + +/-- The offset of the adjacency list of vertex `u` in the edge array: the number of edges out of +the vertices below `u`. The offset of `n` is the number `m` of edges. -/ +def offset (u : ℕ) : ℕ := + (((List.finRange G.n).take u).map fun x => (G.adj x).length).sum + +/-- The input words: `n`, the source, the offsets of the vertices `0, …, n`, and then the edges of +the lists `adj 0, …, adj (n - 1)` in order, each as its target and its weight. Edge `j` in this +order has its target at word `n + 3 + 2 j` and its weight at word `n + 4 + 2 j`, and the edges +out of `u` are those with `offset u ≤ j < offset (u + 1)`. -/ +def encode : List ℕ := + G.n :: (G.source : ℕ) :: + ((List.range (G.n + 1)).map G.offset ++ + (List.finRange G.n).flatMap fun u => (G.adj u).flatMap fun e => [(e.1 : ℕ), e.2]) + +/-- The output address of vertex `v`: the `n` cells just past the input, which has `n + 2 m + 3` +words and so ends at address `n + 2 m + 3`. -/ +def outAddr (v : Fin G.n) : ℕ := G.n + 2 * G.m + 4 + v + +/-- The matrix entry of `u` and `v` is `c` exactly when `c` is the least weight of an edge from +`u` to `v`. -/ +theorem lightest_eq_some_iff {u v : Fin G.n} {c : ℕ} : + G.lightest u v = some c ↔ (v, c) ∈ G.adj u ∧ ∀ c', (v, c') ∈ G.adj u → c ≤ c' := by + simp only [lightest, List.min?_eq_some_iff, List.mem_filterMap, Option.ite_none_right_eq_some, + Option.some.injEq, Prod.exists] + constructor + · rintro ⟨⟨a, b, hab, rfl, rfl⟩, hle⟩ + exact ⟨hab, fun c' hc' => hle c' ⟨a, c', hc', rfl, rfl⟩⟩ + · rintro ⟨hc, hle⟩ + refine ⟨⟨v, c, hc, rfl, rfl⟩, ?_⟩ + rintro b ⟨a, b', hab, rfl, rfl⟩ + exact hle b' hab + +/-- The matrix entry of `u` and `v` is absent exactly when there is no edge from `u` to `v`. -/ +theorem lightest_eq_none_iff {u v : Fin G.n} : + G.lightest u v = none ↔ ∀ c, (v, c) ∉ G.adj u := by + simp only [lightest, List.min?_eq_none_iff, List.filterMap_eq_nil_iff] + constructor + · intro h c hc + simpa using h _ hc + · rintro h ⟨a, b⟩ hab + simp only [ite_eq_right_iff, reduceCtorEq, imp_false] + rintro rfl + exact h b hab + +end ListInstance + +/-- **Shortest paths in `O((m + 1) log n)` time.** At every word size `w` that holds +`8 (n + m + 2) (W + 2)`, where `m` is the number of edges and `W` the largest weight, the flat +program `p`, started on the input words of an adjacency-list instance, halts within +`K (m + 1) (⌊log₂ n⌋ + 1)` instructions, the final `halt` included, and its memory then holds at +`G.outAddr v` the encoded distance from the source to `v`, for every vertex `v`. The register +value returned by `halt` is not part of the result. -/ +def LogLinearSpec {r : ℕ} (p : Program .exact r) (K : ℕ) : Prop := + ∀ (G : ListInstance) (w : ℕ), Fits w (8 * (G.n + G.m + 2) * (G.maxWeight + 2)) → + ∃ t steps, Runs w p (.init G.encode) t steps ∧ (∃ v, step w p t = .halt v) ∧ + steps + 1 ≤ K * (G.m + 1) * (Nat.log 2 G.n + 1) ∧ + ∀ v : Fin G.n, t.data.mem (G.outAddr v) = encodeDist (G.dist G.source v) + +/-- **Log-linear main theorem (statement).** Some word-RAM program over the exact basis computes +all shortest-path distances from the source of a weighted directed graph with `n` vertices and +`m` edges, given by adjacency lists, within `O((m + 1) log n)` instructions. -/ +def LogLinearStatement : Prop := + ∃ r, ∃ p : Program .exact r, ∃ K, LogLinearSpec p K + +end Cslib.Algorithms.WordRAM.Dijkstra diff --git a/Complexity/Algorithms/WordRAM/Dijkstra/Invariant.lean b/Complexity/Algorithms/WordRAM/Dijkstra/Invariant.lean new file mode 100644 index 0000000000..effd609f5a --- /dev/null +++ b/Complexity/Algorithms/WordRAM/Dijkstra/Invariant.lean @@ -0,0 +1,257 @@ +/- +Copyright (c) 2026 Renfei Zhou. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: Renfei Zhou +-/ + +module +-- upstream: candidate + +public import Complexity.Algorithms.WordRAM.Dijkstra.Statement + +/-! +# Dijkstra's algorithm: the invariant + +This file proves the graph theory behind Dijkstra's algorithm [Dijkstra][Dijkstra1959] +[Cormen et al., Section 24.3][CLRS2009], independently of the machine. The algorithm keeps a set +`S` of visited vertices and a tentative distance `D v` for every vertex, and +`Cslib.Algorithms.WordRAM.Dijkstra.Instance.Inv` states what it maintains: `D` bounds the +distance from above, is exact on `S`, and is already relaxed along every edge out of `S`. + +The key step (`Cslib.Algorithms.WordRAM.Dijkstra.Instance.Inv.step`) is the classical cut +argument: on every path from the source to an unvisited vertex, either the source is unvisited, +with tentative distance `0`, or the path has an edge from `S` to an unvisited vertex. Either way +its weight is at least the least tentative distance outside `S`. Hence an unvisited +vertex of least finite tentative distance has the exact distance, and relaxing its edges keeps +the invariant with that vertex visited. + +## Main definitions + +* `Cslib.Algorithms.WordRAM.Dijkstra.decodeDist`: the distance a word encodes. +* `Cslib.Algorithms.WordRAM.Dijkstra.Instance.relax`: relaxation of the edges out of a vertex. +* `Cslib.Algorithms.WordRAM.Dijkstra.Instance.Inv`: the invariant of Dijkstra's algorithm. + +## Main results + +* `Cslib.Algorithms.WordRAM.Dijkstra.decodeDist_encodeDist`: decoding inverts the encoding of + distances. +* `Cslib.Algorithms.WordRAM.Dijkstra.decodeDist_le`: decoding preserves order between nonzero + words. +* `Cslib.Algorithms.WordRAM.Dijkstra.Instance.le_maxWeight`: every edge weight is at most the + largest one. +* `Cslib.Algorithms.WordRAM.Dijkstra.Instance.dist_le_edge`: the distance to the head of an edge + is at most the distance to its tail plus its weight. +* `Cslib.Algorithms.WordRAM.Dijkstra.Instance.Inv.init`: the invariant holds initially. +* `Cslib.Algorithms.WordRAM.Dijkstra.Instance.Inv.step`: visiting an unvisited vertex of least + finite tentative distance and relaxing its edges keeps the invariant. +* `Cslib.Algorithms.WordRAM.Dijkstra.Instance.Inv.eq_dist`: once every unvisited vertex has + infinite tentative distance, the tentative distances are the distances. + +## References + +* [E. W. Dijkstra, *A note on two problems in connexion with graphs*][Dijkstra1959] +* [T. H. Cormen, C. E. Leiserson, R. L. Rivest and C. Stein, *Introduction to + Algorithms*][CLRS2009] +-/ + +@[expose] public section + +namespace Cslib.Algorithms.WordRAM.Dijkstra + +/-- The distance a word encodes: `⊤` for `0` and `d` for `d + 1`. -/ +def decodeDist (x : ℕ) : ℕ∞ := if x = 0 then ⊤ else ((x - 1 : ℕ) : ℕ∞) + +@[simp] +theorem encodeDist_decodeDist (x : ℕ) : encodeDist (decodeDist x) = x := by + cases x with + | zero => simp [decodeDist, encodeDist] + | succ x => simp [decodeDist, encodeDist] + +/-- Decoding inverts the encoding of distances: a distance word decodes to its distance. -/ +@[simp] +theorem decodeDist_encodeDist (d : ℕ∞) : decodeDist (encodeDist d) = d := by + induction d using ENat.recTopCoe with + | top => simp [decodeDist, encodeDist] + | coe d => simp [decodeDist, encodeDist] + +/-- A nonzero distance word decodes to a finite distance. -/ +theorem decodeDist_ne_top {b : ℕ} (hb : b ≠ 0) : decodeDist b ≠ ⊤ := by + simp [decodeDist, hb] + +/-- Decoding preserves order between nonzero distance words. -/ +theorem decodeDist_le {a b : ℕ} (ha : a ≠ 0) (hb : b ≠ 0) (h : a ≤ b) : + decodeDist a ≤ decodeDist b := by + simp only [decodeDist, ha, hb, ite_false] + exact_mod_cast Nat.sub_le_sub_right h 1 + +namespace Instance + +variable (I : Instance) + +/-- Every edge weight is at most the largest one. -/ +theorem le_maxWeight {u v : Fin I.n} {c : ℕ} (h : I.weight u v = some c) : c ≤ I.maxWeight := by + have hle := Finset.le_sup (f := fun uv : Fin I.n × Fin I.n => + (I.weight uv.1 uv.2).getD 0) (Finset.mem_univ (u, v)) + simpa only [h, Option.getD_some, maxWeight] using hle + +/-- Relaxation of the edges out of `u`: the tentative distance of `v` becomes at most that of `u` +plus the weight of the edge from `u` to `v`, if there is one. -/ +noncomputable def relax (D : Fin I.n → ℕ∞) (u v : Fin I.n) : ℕ∞ := + match I.weight u v with + | none => D v + | some c => min (D v) (D u + c) + +/-- **The invariant of Dijkstra's algorithm**, for the visited set `S` and the tentative +distances `D`. -/ +structure Inv (S : Finset (Fin I.n)) (D : Fin I.n → ℕ∞) : Prop where + /-- The source has tentative distance `0`. -/ + source : D I.source = 0 + /-- Tentative distances bound the distances from above. -/ + dist_le : ∀ v, I.dist I.source v ≤ D v + /-- Tentative distances are exact on visited vertices. -/ + eq_of_mem : ∀ v ∈ S, D v = I.dist I.source v + /-- Every edge from a visited vertex to an unvisited one is relaxed. -/ + relaxed : ∀ a ∈ S, ∀ b ∉ S, ∀ c, I.weight a b = some c → D b ≤ D a + c + /-- Finite tentative distances are at most `|S|` times the largest weight. -/ + bound : ∀ v, D v ≠ ⊤ → D v ≤ ((S.card * I.maxWeight : ℕ) : ℕ∞) + +/-- A distance is at most the weight of any path between its endpoints. -/ +private theorem dist_le_addWeight {a b : I.Vertex} (p : Quiver.Path a b) : + I.dist a b ≤ p.addWeight I.arrowWeight := + iInf_le _ p + +/-- Extending a path by an edge bounds the distance to the new endpoint. -/ +theorem dist_le_edge {a b : Fin I.n} {c : ℕ} (he : I.weight a b = some c) : + I.dist I.source b ≤ I.dist I.source a + c := by + change I.dist I.source b ≤ + (⨅ p : Quiver.Path (V := I.Vertex) I.source a, p.addWeight I.arrowWeight) + c + rw [ENat.iInf_add] + apply le_iInf + intro p + let e : Quiver.Hom (V := I.Vertex) a b := ⟨c, he⟩ + exact (I.dist_le_addWeight (p.cons e)).trans_eq + (Quiver.Path.addWeight_cons (V := I.Vertex) I.arrowWeight p e) + +/-- Relaxation never increases a tentative distance. -/ +private theorem relax_le (D : Fin I.n → ℕ∞) (u v : Fin I.n) : + I.relax D u v ≤ D v := by + cases he : I.weight u v with + | none => simp [relax, he] + | some c => simpa only [relax, he] using min_le_left (D v) (D u + c) + +/-- Every path to an unvisited vertex bounds an unvisited tentative distance from above. -/ +private theorem Inv.cut {S : Finset (Fin I.n)} {D : Fin I.n → ℕ∞} (h : I.Inv S D) + {v : I.Vertex} (p : Quiver.Path (V := I.Vertex) I.source v) (hv : v ∉ S) : + ∃ x ∉ S, D x ≤ p.addWeight I.arrowWeight := by + induction p with + | nil => + refine ⟨I.source, hv, ?_⟩ + exact h.source.le.trans_eq (Quiver.Path.addWeight_nil (V := I.Vertex) I.arrowWeight _).symm + | @cons b v p e ih => + rw [Quiver.Path.addWeight_cons (V := I.Vertex) I.arrowWeight p e] + by_cases hb : b ∈ S + · refine ⟨v, hv, ?_⟩ + calc + D v ≤ D b + (e : {c : ℕ // I.weight b v = some c}).1 := + h.relaxed b hb v hv _ e.property + _ = I.dist I.source b + I.arrowWeight e := by rw [h.eq_of_mem b hb]; rfl + _ ≤ p.addWeight I.arrowWeight + I.arrowWeight e := + add_le_add (I.dist_le_addWeight p) le_rfl + · obtain ⟨x, hx, hxp⟩ := ih hb + refine ⟨x, hx, ?_⟩ + exact hxp.trans (le_add_right le_rfl) + +/-- The invariant holds with no vertex visited, the source at tentative distance `0` and every +other vertex at `⊤`. -/ +theorem Inv.init : I.Inv ∅ (Function.update (fun _ => ⊤) I.source 0) := by + classical + refine ⟨by simp, ?_, ?_, ?_, ?_⟩ + · intro v + by_cases hv : v = I.source + · subst v + simp [I.dist_self] + · simp [Function.update_of_ne hv] + · simp + · simp + · intro v hv + by_cases heq : v = I.source + · subst v + simp + · simp [Function.update_of_ne heq] at hv + +/-- **Dijkstra's step.** Visiting an unvisited vertex `u` of least finite tentative distance and +relaxing its edges keeps the invariant. -/ +theorem Inv.step {S : Finset (Fin I.n)} {D : Fin I.n → ℕ∞} (h : I.Inv S D) {u : Fin I.n} + (hu : u ∉ S) (hmin : ∀ v ∉ S, D u ≤ D v) (hfin : D u ≠ ⊤) : + I.Inv (insert u S) (I.relax D u) := by + classical + have hexact : D u = I.dist I.source u := by + apply le_antisymm _ (h.dist_le u) + apply le_iInf + intro p + obtain ⟨x, hx, hxp⟩ := h.cut I p hu + exact (hmin x hx).trans hxp + have hlower (v) : I.dist I.source v ≤ I.relax D u v := by + cases he : I.weight u v with + | none => simpa only [relax, he] using h.dist_le v + | some c => + simp only [relax, he] + exact le_min (h.dist_le v) (hexact ▸ I.dist_le_edge he) + have hold (a) (ha : a ∈ S) : I.relax D u a = D a := by + apply le_antisymm (I.relax_le D u a) + rw [h.eq_of_mem a ha] + exact hlower a + have hself : I.relax D u u = D u := by + apply le_antisymm (I.relax_le D u u) + rw [hexact] + exact hlower u + have hcard : S.card * I.maxWeight ≤ (insert u S).card * I.maxWeight := + Nat.mul_le_mul_right _ (Finset.card_le_card (Finset.subset_insert _ _)) + refine ⟨?_, hlower, ?_, ?_, ?_⟩ + · exact le_antisymm ((I.relax_le D u I.source).trans h.source.le) zero_le + · intro a ha + rcases Finset.mem_insert.mp ha with rfl | ha + · exact hself.trans hexact + · exact (hold a ha).trans (h.eq_of_mem a ha) + · intro a ha b hb c he + rcases Finset.mem_insert.mp ha with rfl | ha + · rw [hself] + simp only [relax, he] + exact min_le_right _ _ + · rw [hold a ha] + exact (I.relax_le D u b).trans + (h.relaxed a ha b (fun hbS => hb (Finset.mem_insert_of_mem hbS)) c he) + · intro v hv + by_cases hd : D v = ⊤ + · cases he : I.weight u v with + | none => simp [relax, he, hd] at hv + | some c => + have hc : (c : ℕ∞) ≤ (I.maxWeight : ℕ∞) := by + exact_mod_cast I.le_maxWeight he + calc + I.relax D u v ≤ D u + c := by simp [relax, he] + _ ≤ ((S.card * I.maxWeight : ℕ) : ℕ∞) + I.maxWeight := + add_le_add (h.bound u hfin) hc + _ = (((insert u S).card * I.maxWeight : ℕ) : ℕ∞) := by + rw [Finset.card_insert_of_notMem hu, Nat.add_mul, Nat.one_mul, Nat.cast_add] + · exact (I.relax_le D u v).trans ((h.bound v hd).trans (by exact_mod_cast hcard)) + +/-- **Termination.** Once every unvisited vertex has infinite tentative distance, the tentative +distances are the distances. -/ +theorem Inv.eq_dist {S : Finset (Fin I.n)} {D : Fin I.n → ℕ∞} (h : I.Inv S D) + (hdone : ∀ v ∉ S, D v = ⊤) (v : Fin I.n) : D v = I.dist I.source v := by + by_cases hv : v ∈ S + · exact h.eq_of_mem v hv + · rw [hdone v hv] + apply le_antisymm _ le_top + apply le_iInf + intro p + obtain ⟨x, hx, hxp⟩ := h.cut I p hv + rw [hdone x hx] at hxp + exact hxp + +end Instance + +end Cslib.Algorithms.WordRAM.Dijkstra + +end diff --git a/Complexity/Algorithms/WordRAM/Dijkstra/ListEncode.lean b/Complexity/Algorithms/WordRAM/Dijkstra/ListEncode.lean new file mode 100644 index 0000000000..6e5dd81409 --- /dev/null +++ b/Complexity/Algorithms/WordRAM/Dijkstra/ListEncode.lean @@ -0,0 +1,257 @@ +/- +Copyright (c) 2026 Renfei Zhou. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: Renfei Zhou +-/ + +module +-- upstream: candidate + +public import Complexity.Algorithms.WordRAM.Dijkstra.HeapStatement +public import Complexity.Algorithms.WordRAM.Dijkstra.Invariant +import Mathlib.Algebra.BigOperators.Fin + +/-! +# Dijkstra's algorithm with a binary heap: the input words + +This file reads the input words `Cslib.Algorithms.WordRAM.Dijkstra.ListInstance.encode` of an +adjacency-list instance back from the initial memory of a word RAM: `n`, the source, the offsets +of the adjacency lists and the targets and weights of the edges each sit at the address that +`Cslib.Algorithms.WordRAM.Dijkstra.heapDijkstra` reads them from, and memory past the input is +`0`. + +It also defines the effect of scanning the edges out of a vertex on one distance word +(`Cslib.Algorithms.WordRAM.Dijkstra.edgesWord`), and shows that it is the relaxation of the +weight-matrix instance (`Cslib.Algorithms.WordRAM.Dijkstra.decode_edgesWord`): the scan of all +parallel edges to a vertex has the effect of the lightest one. + +## Main definitions + +* `Cslib.Algorithms.WordRAM.Dijkstra.ListInstance.dbase`: the first address past the input. +* `Cslib.Algorithms.WordRAM.Dijkstra.edgeWord`, `Cslib.Algorithms.WordRAM.Dijkstra.edgesWord`: + the relaxation of a distance word through one edge, and through a list of edges. + +## Main results + +* `Cslib.Algorithms.WordRAM.Dijkstra.ListInstance.encode_length`: the input has `n + 2 m + 3` + words. +* `Cslib.Algorithms.WordRAM.Dijkstra.ListInstance.init_offset`, + `Cslib.Algorithms.WordRAM.Dijkstra.ListInstance.init_target`, + `Cslib.Algorithms.WordRAM.Dijkstra.ListInstance.init_weight`: the input words in the initial + memory. +* `Cslib.Algorithms.WordRAM.Dijkstra.decode_edgesWord`: the scan of the edges out of a vertex + relaxes the tentative distances. + +## References + +* [E. W. Dijkstra, *A note on two problems in connexion with graphs*][Dijkstra1959] +* [T. H. Cormen, C. E. Leiserson, R. L. Rivest and C. Stein, *Introduction to + Algorithms*][CLRS2009] +-/ + +@[expose] public section + +namespace Cslib.Algorithms.WordRAM.Dijkstra + +open Cslib.WordRAM + +namespace ListInstance + +variable (G : ListInstance) {r : ℕ} + +/-- The first address past the input words: the base of the distance words. -/ +def dbase : ℕ := G.n + 2 * G.m + 4 + +/-- The edges of all lists, in order. -/ +private def edges : List (Fin G.n × ℕ) := (List.finRange G.n).flatMap G.adj + +/-- The offset of `u` is the number of edges out of the vertices below `u`. -/ +private theorem offset_eq_length (u : ℕ) : + G.offset u = (((List.finRange G.n).take u).flatMap G.adj).length := by + simp [offset, List.length_flatMap] + +/-- The edge list has `m` edges. -/ +private theorem length_edges : G.edges.length = G.m := by + simp [edges, m, List.length_flatMap, ← Fin.sum_univ_def] + +/-- Edge `offset u + i` of the edge list is edge `i` out of `u`. -/ +private theorem edges_getElem (u : Fin G.n) {i : ℕ} (hi : i < (G.adj u).length) : + G.edges[G.offset u + i]? = some (G.adj u)[i] := by + have hu : (u : ℕ) < (List.finRange G.n).length := by simp + have hs : List.finRange G.n = (List.finRange G.n).take u ++ + (u :: (List.finRange G.n).drop (u + 1)) := by + conv_lhs => rw [← List.take_append_drop u (List.finRange G.n)] + rw [List.drop_eq_getElem_cons hu] + simp + rw [edges, hs, List.flatMap_append, List.flatMap_cons, offset_eq_length, + List.getElem?_append_right (by omega), Nat.add_sub_cancel_left, + List.getElem?_append_left hi, List.getElem?_eq_getElem hi] + +/-- Word `2 j + k` of the edge words, for `k < 2`, is word `k` of edge `j`. -/ +private theorem getD_flatMap_pairs (l : List (Fin G.n × ℕ)) (j : ℕ) (hj : j < l.length) : + (l.flatMap fun e => [(e.1 : ℕ), e.2]).getD (2 * j) 0 = l[j].1 ∧ + (l.flatMap fun e => [(e.1 : ℕ), e.2]).getD (2 * j + 1) 0 = l[j].2 := by + induction l generalizing j with + | nil => simp at hj + | cons a l ih => + cases j with + | zero => simp + | succ j => + have h := ih j (by simpa using hj) + rw [show 2 * (j + 1) = 2 * j + 2 by ring, show 2 * j + 2 + 1 = 2 * j + 1 + 2 by ring] + simpa [List.getD_eq_getElem?_getD] using h + +/-- The edge words are the edge list, two words per edge. -/ +private theorem edgeWords_eq : ((List.finRange G.n).flatMap fun u => + (G.adj u).flatMap fun e => [(e.1 : ℕ), e.2]) = + G.edges.flatMap fun e => [(e.1 : ℕ), e.2] := by + rw [edges, List.flatMap_assoc] + +/-- Word `k + n + 3` of the input is word `k` of the edge words. -/ +private theorem encode_getD_edge (k : ℕ) : + G.encode.getD (k + G.n + 3) 0 = + (G.edges.flatMap fun e => [(e.1 : ℕ), e.2]).getD k 0 := by + rw [encode, show k + G.n + 3 = (k + G.n + 1) + 1 + 1 by ring, + List.getD_cons_succ, List.getD_cons_succ, edgeWords_eq, List.getD_eq_getElem?_getD, + List.getElem?_append_right (by simp), List.getD_eq_getElem?_getD] + simp + +/-- The input has `n + 2 m + 3` words. -/ +theorem encode_length : G.encode.length = G.n + 2 * G.m + 3 := by + have h : ∀ l : List (Fin G.n × ℕ), + (l.flatMap fun e => [(e.1 : ℕ), e.2]).length = 2 * l.length := by + intro l + induction l with + | nil => simp + | cons a l ih => + simp only [List.flatMap_cons, List.cons_append, List.nil_append, List.length_cons, ih] + ring + simp only [encode, edgeWords_eq, List.length_cons, List.length_append, List.length_map, + List.length_range, h, length_edges, Nat.add_right_cancel_iff] + ring + +/-- The offset of `u + 1` is the offset of `u` plus the length of the list of `u`. -/ +theorem offset_succ {u : ℕ} (hu : u < G.n) : + G.offset (u + 1) = G.offset u + (G.adj ⟨u, hu⟩).length := by + have hl : u < (List.finRange G.n).length := by simpa using hu + simp [offset, List.take_add_one, List.getElem?_eq_getElem hl] + +/-- The offset of `n` is the number of edges. -/ +theorem offset_n : G.offset G.n = G.m := by + rw [offset, List.take_of_length_le (by simp), ← Fin.sum_univ_def] + rfl + +/-- The offsets grow. -/ +theorem offset_mono {u v : ℕ} (h : u ≤ v) : G.offset u ≤ G.offset v := + ((List.take_sublist_take_left h).map _).sum_le_sum fun _ _ => Nat.zero_le _ + +/-- Address `1` initially holds `n`. -/ +theorem init_n : (Data.init G.encode : Data r).mem 1 = G.n := rfl + +/-- Address `2` initially holds the source. -/ +theorem init_source : (Data.init G.encode : Data r).mem 2 = G.source := rfl + +/-- Address `u + 3` initially holds the offset of `u`, for `u ≤ n`. -/ +theorem init_offset {u : ℕ} (hu : u ≤ G.n) : + (Data.init G.encode : Data r).mem (u + 3) = G.offset u := by + rw [Data.init_mem_succ, encode, List.getD_cons_succ, List.getD_cons_succ, + List.getD_eq_getElem?_getD, List.getElem?_append_left (by simp; omega)] + simp [Nat.lt_succ_of_le hu] + +/-- The target of edge `i` out of `u` is initially at address `2 (offset u + i) + n + 4`. -/ +theorem init_target (u : Fin G.n) {i : ℕ} (hi : i < (G.adj u).length) : + (Data.init G.encode : Data r).mem (2 * (G.offset u + i) + G.n + 4) = ((G.adj u)[i].1 : ℕ) := by + rw [show 2 * (G.offset u + i) + G.n + 4 = 2 * (G.offset u + i) + G.n + 3 + 1 by ring, + Data.init_mem_succ, encode_getD_edge] + obtain ⟨hl, he⟩ := List.getElem?_eq_some_iff.mp (G.edges_getElem u hi) + rw [(G.getD_flatMap_pairs _ _ hl).1, he] + +/-- The weight of edge `i` out of `u` is initially at address `2 (offset u + i) + n + 5`. -/ +theorem init_weight (u : Fin G.n) {i : ℕ} (hi : i < (G.adj u).length) : + (Data.init G.encode : Data r).mem (2 * (G.offset u + i) + G.n + 5) = (G.adj u)[i].2 := by + rw [show 2 * (G.offset u + i) + G.n + 5 = (2 * (G.offset u + i) + 1) + G.n + 3 + 1 by ring, + Data.init_mem_succ, encode_getD_edge] + obtain ⟨hl, he⟩ := List.getElem?_eq_some_iff.mp (G.edges_getElem u hi) + rw [(G.getD_flatMap_pairs _ _ hl).2, he] + +/-- Memory past the input is initially `0`. -/ +theorem init_past {a : ℕ} (ha : G.dbase ≤ a) : (Data.init G.encode : Data r).mem a = 0 := by + have hz : a ≠ 0 := by unfold dbase at ha; omega + have hl : G.encode.length ≤ a - 1 := by rw [encode_length]; unfold dbase at ha; omega + obtain ⟨a, rfl⟩ := Nat.exists_eq_succ_of_ne_zero hz + simp [List.getD_eq_getElem?_getD, List.getElem?_eq_none (by simpa using hl)] + +/-- Every edge weight is at most the largest weight. -/ +theorem le_maxWeight {u : Fin G.n} {e : Fin G.n × ℕ} (he : e ∈ G.adj u) : e.2 ≤ G.maxWeight := by + have h1 : e.2 ≤ ((G.adj u).map Prod.snd).foldr max 0 := + List.le_max_of_le (List.mem_map_of_mem he) le_rfl + exact h1.trans (Finset.le_sup (f := fun u => ((G.adj u).map Prod.snd).foldr max 0) + (Finset.mem_univ u)) + +/-- The largest weight of the weight-matrix instance is at most the largest edge weight. -/ +theorem toInstance_maxWeight_le : G.toInstance.maxWeight ≤ G.maxWeight := by + refine Finset.sup_le fun uv _ => ?_ + change (G.lightest uv.1 uv.2).getD 0 ≤ G.maxWeight + cases h : G.lightest uv.1 uv.2 with + | none => simp + | some c => exact G.le_maxWeight ((G.lightest_eq_some_iff.mp h).1) + +end ListInstance + +/-- The relaxation of the distance word `d` through an edge of weight `c` from a vertex with +distance word `b ≠ 0`. -/ +def edgeWord (b c d : ℕ) : ℕ := if d = 0 then b + c else if b + c < d then b + c else d + +/-- The relaxation of the distance word `d` of `v` through the edges to `v` in the list `es`, in +order, from a vertex with distance word `b`. -/ +def edgesWord {N : ℕ} (b : ℕ) (es : List (Fin N × ℕ)) (v : Fin N) (d : ℕ) : ℕ := + es.foldl (fun d e => if e.1 = v then edgeWord b e.2 d else d) d + +/-- The relaxation through one edge decodes to the relaxation of the tentative distance. -/ +private theorem decode_edgeWord {b : ℕ} (hb : b ≠ 0) (c d : ℕ) : + decodeDist (edgeWord b c d) = min (decodeDist d) (decodeDist b + c) := by + obtain ⟨b, rfl⟩ := Nat.exists_eq_succ_of_ne_zero hb + have he : b.succ + c - 1 = b + c := by omega + cases d with + | zero => simp [edgeWord, decodeDist, he] + | succ d => + simp only [edgeWord, Nat.succ_ne_zero, ite_false] + by_cases h : b.succ + c < d.succ + · have hc : (b : ℕ∞) + c ≤ d := by exact_mod_cast (by omega : b + c ≤ d) + simp [h, decodeDist, he, hc] + · have hc : (d : ℕ∞) ≤ (b : ℕ∞) + c := by exact_mod_cast (by omega : d ≤ b + c) + simp [h, decodeDist, hc] + +/-- The relaxation through the edges to `v` in `es` decodes to the relaxation through the +lightest of them. -/ +private theorem decode_edgesWord_list {N : ℕ} {b : ℕ} (hb : b ≠ 0) (es : List (Fin N × ℕ)) + (v : Fin N) (d : ℕ) : + decodeDist (edgesWord b es v d) = + match (es.filterMap fun e => if e.1 = v then some e.2 else none).min? with + | none => decodeDist d + | some c => min (decodeDist d) (decodeDist b + c) := by + induction es generalizing d with + | nil => rfl + | cons e es ih => + rw [edgesWord, List.foldl_cons, ← edgesWord, ih] + by_cases he : e.1 = v + · simp only [List.filterMap_cons, he, ↓reduceIte, List.min?_cons] + cases es.filterMap (fun e => if e.1 = v then some e.2 else none) |>.min? with + | none => simp [decode_edgeWord hb] + | some c => + simp only [decode_edgeWord hb, Option.elim_some] + rw [min_assoc, min_add_add_left, Nat.mono_cast.map_min] + · simp only [List.filterMap_cons, he, ↓reduceIte] + +/-- **The scan of the edges out of `u` relaxes the tentative distances.** If `b ≠ 0` is the +distance word of `u` and `d` that of `v`, the relaxation of `d` through the edges from `u` to `v` +decodes to the relaxation of the tentative distance of `v` in the weight-matrix instance. -/ +theorem decode_edgesWord (G : ListInstance) {D : Fin G.n → ℕ∞} {u v : Fin G.n} {b d : ℕ} + (hb : b ≠ 0) (hu : D u = decodeDist b) (hv : D v = decodeDist d) : + decodeDist (edgesWord b (G.adj u) v d) = G.toInstance.relax D u v := by + rw [decode_edgesWord_list hb, ← hu, ← hv] + rfl + +end Cslib.Algorithms.WordRAM.Dijkstra + +end diff --git a/Complexity/Algorithms/WordRAM/Dijkstra/LogLinear.lean b/Complexity/Algorithms/WordRAM/Dijkstra/LogLinear.lean new file mode 100644 index 0000000000..450a4af48d --- /dev/null +++ b/Complexity/Algorithms/WordRAM/Dijkstra/LogLinear.lean @@ -0,0 +1,445 @@ +/- +Copyright (c) 2026 Renfei Zhou. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: Renfei Zhou +-/ + +module +-- upstream: candidate + +public import Complexity.Algorithms.WordRAM.Dijkstra.HeapRelax +public import Complexity.Computability.WordRAM.Compile + +/-! +# Dijkstra's algorithm with a binary heap in `O((m + 1) (⌊log₂ n⌋ + 1))` time + +This file proves that `Cslib.Algorithms.WordRAM.Dijkstra.heapDijkstra` computes the distances +from the source of an adjacency-list instance within `O((m + 1) (⌊log₂ n⌋ + 1))` instructions, +and so proves `Cslib.Algorithms.WordRAM.Dijkstra.LogLinearStatement`. + +The rounds keep the invariant of Dijkstra's algorithm +(`Cslib.Algorithms.WordRAM.Dijkstra.Instance.Inv`) of the weight-matrix instance, for the visited +vertices (discovered and out of the heap) and the tentative distances in the output cells, along +with the heap invariant. A round costs `O((deg u + 1) (⌊log₂ n⌋ + 1))` for its vertex `u`, which +is paid by a potential that counts twice the edges out of the unvisited vertices plus the heap +size: the round visits `u`, so it removes `deg u` edges from the count and one entry from the +heap, and it inserts at most `deg u` entries. + +## Main definitions + +* `Cslib.Algorithms.WordRAM.Dijkstra.ListInstance.distOf`, + `Cslib.Algorithms.WordRAM.Dijkstra.ListInstance.visitedOf`: the tentative distances and the + visited vertices that a machine state holds. +* `Cslib.Algorithms.WordRAM.Dijkstra.ListInstance.remaining`: the number of edges out of the + unvisited vertices. +* `Cslib.Algorithms.WordRAM.Dijkstra.ListInstance.OuterInv`: the invariant of the rounds. +* `Cslib.Algorithms.WordRAM.Dijkstra.ListInstance.potential`: the potential that pays for them. + +## Main results + +* `Cslib.Algorithms.WordRAM.Dijkstra.heapSetup_spec`: the setup establishes the invariant of the + rounds within `25` instructions, with potential at most `100 (⌊log₂ n⌋ + 1) (2 m + 1)`. +* `Cslib.Algorithms.WordRAM.Dijkstra.heapRound_spec`: one round keeps the invariant of the + rounds, and the drop of the potential pays for its running time. +* `Cslib.Algorithms.WordRAM.Dijkstra.heapExit_spec`: when the rounds stop, with an empty heap, + the output cells hold the encoded distances. +* `Cslib.Algorithms.WordRAM.Dijkstra.heapDijkstra_spec`: the structured program leaves the + encoded distances in the output cells within `200 (m + 1) (⌊log₂ n⌋ + 1) + 27` instructions. +* `Cslib.Algorithms.WordRAM.Dijkstra.logLinearSpec`: its compiled program meets + `Cslib.Algorithms.WordRAM.Dijkstra.LogLinearSpec` with `K = 230`. +* `Cslib.Algorithms.WordRAM.Dijkstra.logLinearStatement`: the log-linear main theorem. + +## References + +* [E. W. Dijkstra, *A note on two problems in connexion with graphs*][Dijkstra1959] +* [T. H. Cormen, C. E. Leiserson, R. L. Rivest and C. Stein, *Introduction to + Algorithms*][CLRS2009] +-/ + +@[expose] public section + +namespace Cslib.Algorithms.WordRAM.Dijkstra + +open Cslib.WordRAM HeapReg + +namespace ListInstance + +variable (G : ListInstance) + +/-- The tentative distances held in the distance words of `y`. -/ +noncomputable def distOf (y : Data 27) (v : Fin G.n) : ℕ∞ := + decodeDist (y.mem (G.dbase + v)) + +/-- The visited vertices of `y`: those that are discovered, with a nonzero distance word, but +not in the heap, with position `0`. -/ +def visitedOf (y : Data 27) : Finset (Fin G.n) := + Finset.univ.filter fun v => y.mem (G.dbase + v) ≠ 0 ∧ y.mem (G.dbase + 2 * G.n + 1 + v) = 0 + +/-- The number of edges out of the unvisited vertices of `y`. -/ +def remaining (y : Data 27) : ℕ := + ∑ v ∈ Finset.univ.filter (· ∉ G.visitedOf y), (G.adj v).length + +/-- **The potential** of `y`, which pays for the remaining rounds: `100 (⌊log₂ n⌋ + 1)` times +twice the number of edges out of the unvisited vertices plus the heap size. -/ +def potential (y : Data 27) : ℕ := + 100 * (Nat.log 2 G.n + 1) * (2 * G.remaining y + y.regs sz) + +/-- **The invariant of the rounds.** The constant registers hold their values, the input is +intact, the heap invariant holds, and the memory satisfies the invariant of Dijkstra's algorithm +for the weight-matrix instance. -/ +structure OuterInv (y : Data 27) : Prop where + zero_eq : y.regs zero = 0 + one_eq : y.regs one = 1 + two_eq : y.regs two = 2 + three_eq : y.regs three = 3 + eb_eq : y.regs eb = G.n + 4 + db_eq : y.regs db = G.dbase + hb_eq : y.regs hb = G.dbase + G.n + pb_eq : y.regs pb = G.dbase + 2 * G.n + 1 + input : ∀ a < G.dbase, y.mem a = (Data.init G.encode : Data 27).mem a + heap : HeapInv G.n G.dbase (y.regs sz) y.mem + inv : G.toInstance.Inv (G.visitedOf y) (G.distOf y) + +end ListInstance + +variable {G : ListInstance} {w : ℕ} + +/-- The allowed word size holds every address of the three arrays. -/ +private theorem fits_arrays (hw : Fits w (8 * (G.n + G.m + 2) * (G.maxWeight + 2))) : + Fits w (G.dbase + 3 * G.n + 1) := + hw.of_le (by unfold ListInstance.dbase; nlinarith [Nat.zero_le G.maxWeight]) + +/-- The setup establishes the invariant of the rounds, with the source as the only heap entry. -/ +theorem heapSetup_spec (hw : Fits w (8 * (G.n + G.m + 2) * (G.maxWeight + 2))) : + ∃ y, ExecWithin w heapSetup (.init G.encode) y 25 ∧ G.OuterInv y ∧ + G.potential y ≤ 100 * (Nat.log 2 G.n + 1) * (2 * G.m + 1) := by + classical + have hf := fits_arrays hw + have hs := G.source.isLt + have hDb' : G.m + G.m + G.n + 3 + 1 = G.dbase := by unfold ListInstance.dbase; omega + have hpb' : G.dbase + G.n + G.n + 1 = G.dbase + 2 * G.n + 1 := by omega + have hoff : (Data.init G.encode : Data 27).mem (G.n + 3) = G.m := by + rw [G.init_offset le_rfl, G.offset_n] + have h1 : Fits w 1 := hf.of_le (by omega) + have h2 : Fits w 2 := hf.of_le (by omega) + have h3 : Fits w 3 := hf.of_le (by omega) + have hn : Fits w G.n := hf.of_le (by omega) + have hn3 : Fits w (G.n + 3) := hf.of_le (by unfold ListInstance.dbase; omega) + have hn4 : Fits w (G.n + 3 + 1) := hf.of_le (by unfold ListInstance.dbase; omega) + have hm2 : Fits w (G.m + G.m) := hf.of_le (by unfold ListInstance.dbase; omega) + have hm2n : Fits w (G.m + G.m + G.n) := hf.of_le (by unfold ListInstance.dbase; omega) + have hm2n3 : Fits w (G.m + G.m + G.n + 3) := hf.of_le (by unfold ListInstance.dbase; omega) + have hDb : Fits w G.dbase := hf.of_le (by omega) + have hhb : Fits w (G.dbase + G.n) := hf.of_le (by omega) + have hpb₁ : Fits w (G.dbase + G.n + G.n) := hf.of_le (by omega) + have hpb : Fits w (G.dbase + 2 * G.n + 1) := hf.of_le (by omega) + have hds : Fits w (G.dbase + G.source) := hf.of_le (by omega) + have hh1 : Fits w (G.dbase + G.n + 1) := hf.of_le (by omega) + have hps : Fits w (G.dbase + 2 * G.n + 1 + G.source) := hf.of_le (by omega) + let y := (((((((((((((((((((((((((Data.init G.encode : Data 27).set one 1).set two 2).set + three 3).set t 1).set n G.n).set t 2).set u G.source).set t (G.n + 3)).set t G.m).set db + (G.m + G.m)).set db (G.m + G.m + G.n)).set db (G.m + G.m + G.n + 3)).set db G.dbase).set eb + (G.n + 3)).set eb (G.n + 3 + 1)).set hb (G.dbase + G.n)).set pb (G.dbase + G.n + G.n)).set pb + (G.dbase + 2 * G.n + 1)).set t (G.dbase + G.source)).write (G.dbase + G.source) 1).set t + (G.dbase + G.n + 1)).write (G.dbase + G.n + 1) G.source).set t + (G.dbase + 2 * G.n + 1 + G.source)).write (G.dbase + 2 * G.n + 1 + G.source) 1).set sz 1 + have he : ExecWithin w heapSetup (.init G.encode) y 25 := by + ram_exact_eval [heapSetup, y, Data.init_regs, G.init_n, G.init_source, hoff, h1, h2, h3, hn, + hn3, hn4, hm2, hm2n, hm2n3, hDb', hDb, hhb, hpb₁, hpb', hpb, hds, hh1, hps] + have hm (a : ℕ) : y.mem a = Function.update (Function.update (Function.update + (Data.init G.encode : Data 27).mem (G.dbase + G.source) 1) (G.dbase + G.n + 1) G.source) + (G.dbase + 2 * G.n + 1 + G.source) 1 a := rfl + have hD (v : ℕ) (hv : v < G.n) : + y.mem (G.dbase + v) = if v = G.source then 1 else 0 := by + simp only [hm, Function.update_apply] + split_ifs <;> first | rfl | omega | exact G.init_past (by omega) + have hP (v : ℕ) (hv : v < G.n) : + y.mem (G.dbase + 2 * G.n + 1 + v) = if v = G.source then 1 else 0 := by + simp only [hm, Function.update_apply] + split_ifs <;> first | rfl | omega | exact G.init_past (by omega) + have hH : y.mem (G.dbase + G.n + 1) = G.source := by + simp only [hm, Function.update_apply] + split_ifs <;> first | rfl | omega + have hsz : y.regs sz = 1 := by simp [y] + have hvis : G.visitedOf y = ∅ := by + apply Finset.eq_empty_iff_forall_notMem.mpr + intro v + simp only [ListInstance.visitedOf, Finset.mem_filter, Finset.mem_univ, true_and, + hD v v.isLt, hP v v.isLt] + split_ifs <;> simp + have hdist : G.distOf y = Function.update (fun _ => ⊤) G.source 0 := by + funext v + simp only [ListInstance.distOf, hD v v.isLt, Fin.val_eq_val] + by_cases h : v = G.source + · subst v + simp [decodeDist] + · simp [h, decodeDist] + have hheap : HeapInv G.n G.dbase 1 y.mem := by + refine ⟨⟨by omega, ?_, ?_, ?_, ?_, ?_⟩, ?_⟩ + · intro i hi hi' + rw [show i = 1 by omega, hH] + exact hs + · intro i hi hi' + rw [show i = 1 by omega, hH, hP _ hs] + simp + · intro v hv + rw [hP v hv] + split_ifs <;> omega + · intro v hv hne + rw [hP v hv] at hne + split_ifs at hne with h + · subst h + simp [hP _ hv, hH] + · exact absurd rfl hne + · intro v hv hne + rw [hP v hv] at hne + split_ifs at hne with h + · subst h + simp [hD _ hv] + · exact absurd rfl hne + · intro i hi hi' + omega + refine ⟨y, he, ⟨by simp +decide [y], by simp +decide [y], by simp +decide [y], + by simp +decide [y], by simp +decide [y], by simp +decide [y], by simp +decide [y], + by simp +decide [y], ?_, hsz ▸ hheap, ?_⟩, ?_⟩ + · intro a ha + rw [hm, Function.update_of_ne (by omega), Function.update_of_ne (by omega), + Function.update_of_ne (by omega)] + · rw [hvis, hdist] + exact Instance.Inv.init G.toInstance + · have hrem : G.remaining y = G.m := by + simp [ListInstance.remaining, hvis, ListInstance.m] + simp [ListInstance.potential, hrem, hsz] + +/-- **One round** keeps the invariant of the rounds, and its running time and two more +instructions are paid by the drop of the potential. -/ +theorem heapRound_spec (hw : Fits w (8 * (G.n + G.m + 2) * (G.maxWeight + 2))) {y : Data 27} + (hy : G.OuterInv y) (hsz : (Cond.lt zero sz).Holds y) : + ∃ y' T, ExecWithin w heapRound y y' T ∧ G.OuterInv y' ∧ + T + 2 + G.potential y' ≤ G.potential y := by + classical + have hf := fits_arrays hw + have hS1 : 1 ≤ y.regs sz := by + simp only [Cond.holds_lt, hy.zero_eq] at hsz + omega + obtain ⟨x, ex, xr, xu, xbest, xsz, xheap, xpos, xmem⟩ := + extractMin_spec hy.zero_eq hy.one_eq hy.db_eq hy.hb_eq hy.pb_eq rfl hS1 hf hy.heap + have hU : y.mem (G.dbase + G.n + 1) < G.n := hy.heap.vertex_lt 1 le_rfl hS1 + have hPU : y.mem (G.dbase + 2 * G.n + 1 + y.mem (G.dbase + G.n + 1)) = 1 := + hy.heap.pos_vertex 1 le_rfl hS1 + have hB : y.mem (G.dbase + y.mem (G.dbase + G.n + 1)) ≠ 0 := + hy.heap.key_ne _ hU (by rw [hPU]; omega) + have hroot : ∀ v < G.n, y.mem (G.dbase + 2 * G.n + 1 + v) ≠ 0 → + y.mem (G.dbase + y.mem (G.dbase + G.n + 1)) ≤ y.mem (G.dbase + v) := + fun v hv hp => hy.heap.root_le hv hp + generalize hU0 : y.mem (G.dbase + G.n + 1) = U0 at xu xbest xpos hU hPU hB hroot + generalize hB0 : y.mem (G.dbase + U0) = B at xbest hB hroot + let U : Fin G.n := ⟨U0, hU⟩ + -- The four instructions that load the bounds of the edges out of `U`. + have xthree : x.regs three = 3 := (xr three (by decide)).trans hy.three_eq + have xone : x.regs one = 1 := (xr one (by decide)).trans hy.one_eq + have hlt : U0 + 3 + 1 < G.dbase := by unfold ListInstance.dbase; omega + have hf3 : Fits w (U0 + 3) := hf.of_le (by omega) + have hf4 : Fits w (U0 + 3 + 1) := hf.of_le (by omega) + have hx (a : ℕ) (ha : a < G.dbase) : x.mem a = (Data.init G.encode : Data 27).mem a := by + rw [xmem a (Or.inl (by omega)), hy.input a ha] + have hoff₀ : x.mem (U0 + 3) = G.offset U0 := by + rw [hx _ (by omega), G.init_offset (by omega)] + have hoff₁ : x.mem (U0 + 3 + 1) = G.offset (U0 + 1) := by + rw [hx _ (by omega), show U0 + 3 + 1 = U0 + 1 + 3 by omega, G.init_offset (by omega)] + let x₁ := x.set t (U0 + 3) + let x₂ := x₁.set j (G.offset U0) + let x₃ := x₂.set t (U0 + 3 + 1) + let x₄ := x₃.set jend (G.offset (U0 + 1)) + have e₁ : ExecWithin w ram_exact[t ← u + three] x x₁ 1 := by ram_exact_eval [x₁, xu, xthree, hf3] + have e₂ : ExecWithin (B := .exact) w ram_exact[j ← *t] x₁ x₂ 1 := by + ram_exact_eval [x₁, x₂, hoff₀] + have e₃ : ExecWithin w ram_exact[t ← t + one] x₂ x₃ 1 := by ram_exact_eval [x₁, x₂, x₃, xone, hf4] + have e₄ : ExecWithin (B := .exact) w ram_exact[jend ← *t] x₃ x₄ 1 := by + ram_exact_eval [x₁, x₂, x₃, x₄, hoff₁] + have r₄ (i : Fin 27) (hi : i ∉ [t, j, jend]) : x₄.regs i = x.regs i := by + simp only [List.mem_cons, List.not_mem_nil, or_false, not_or] at hi + simp [x₄, x₃, x₂, x₁, hi.1, hi.2.1, hi.2.2] + have reg (i : Fin 27) (h₁ : i ∉ [t, j, jend]) + (h₂ : i ∉ [t, u, best, vk, sz, k, p, vp, dk, dp, l, szp]) : x₄.regs i = y.regs i := + (r₄ i h₁).trans (xr i h₂) + have hxD (v : Fin G.n) : x₄.mem (G.dbase + v) = y.mem (G.dbase + v) := + xmem _ (Or.inl (by have := v.isLt; omega)) + -- Dijkstra's step for the vertex `U`. + have hdu : G.distOf y U = decodeDist B := by + change decodeDist (y.mem (G.dbase + U0)) = _ + rw [hB0] + have hfin : G.distOf y U ≠ ⊤ := hdu ▸ decodeDist_ne_top hB + have hu : U ∉ G.visitedOf y := by simp [ListInstance.visitedOf, U, hPU] + have hmin : ∀ v ∉ G.visitedOf y, G.distOf y U ≤ G.distOf y v := by + intro v hv + simp only [ListInstance.visitedOf, Finset.mem_filter, Finset.mem_univ, true_and, + not_and] at hv + rw [hdu] + by_cases hd : y.mem (G.dbase + v) = 0 + · simp [ListInstance.distOf, hd, decodeDist] + · exact decodeDist_le hB hd (hroot v v.isLt (hv hd)) + have hstep := Instance.Inv.step G.toInstance hy.inv hu hmin hfin + have hfitC : ∀ e ∈ G.adj U, Fits w (B + e.2) := by + intro e he + have hb := hy.inv.bound U hfin + rw [hdu] at hb + simp only [decodeDist, hB, ite_false] at hb + have hb' : B - 1 ≤ (G.visitedOf y).card * G.toInstance.maxWeight := by exact_mod_cast hb + have hc : (G.visitedOf y).card ≤ G.n := by simpa using Finset.card_le_univ (G.visitedOf y) + have hW := G.toInstance_maxWeight_le + have he2 := G.le_maxWeight he + have := Nat.mul_le_mul hc hW + have hB' : B ≤ G.n * G.maxWeight + 1 := by omega + have hle : G.n * G.maxWeight + G.maxWeight + 1 ≤ + 8 * (G.n + G.m + 2) * (G.maxWeight + 2) := by nlinarith + exact hw.of_le (by omega) + have hvis : ∀ v : Fin G.n, x₄.mem (G.dbase + v) ≠ 0 → + x₄.mem (G.dbase + 2 * G.n + 1 + v) = 0 → ∀ c, (v, c) ∈ G.adj U → + x₄.mem (G.dbase + v) ≤ B + c := by + intro v hd hp c hc + rw [hxD] at hd ⊢ + rcases (xpos v v.isLt).mp hp with hp | hv + · have hvV : v ∈ G.visitedOf y := by simp [ListInstance.visitedOf, hd, hp] + have h₁ := hstep.eq_of_mem v (Finset.mem_insert_of_mem hvV) + have h₂ := hy.inv.eq_of_mem v hvV + obtain ⟨c', hc'⟩ : ∃ c', G.lightest U v = some c' := by + cases h : G.lightest U v with + | none => exact absurd hc ((G.lightest_eq_none_iff.mp h) c) + | some c' => exact ⟨c', rfl⟩ + have hle := (G.lightest_eq_some_iff.mp hc').2 c hc + have hrel : G.toInstance.relax (G.distOf y) U v = + min (G.distOf y v) (G.distOf y U + c') := by + simp only [Instance.relax, show G.toInstance.weight U v = some c' from hc'] + have hmin' : G.distOf y v ≤ G.distOf y U + c' := by + rw [← min_eq_left_iff, ← hrel, h₁, h₂] + rw [hdu] at hmin' + simp only [ListInstance.distOf, decodeDist, hd, hB, ite_false] at hmin' + have : y.mem (G.dbase + v) - 1 ≤ B - 1 + c' := by exact_mod_cast hmin' + omega + · have : y.mem (G.dbase + v) = B := by rw [hv, hB0] + omega + obtain ⟨z, ez, zr, zsz, zheap, zD, zvis, zin⟩ := relaxEdges_spec (G := G) (w := w) U (x := x₄) + ((reg zero (by decide) (by decide)).trans hy.zero_eq) + ((reg one (by decide) (by decide)).trans hy.one_eq) + ((reg two (by decide) (by decide)).trans hy.two_eq) + ((reg eb (by decide) (by decide)).trans hy.eb_eq) + ((reg db (by decide) (by decide)).trans hy.db_eq) + ((reg hb (by decide) (by decide)).trans hy.hb_eq) + ((reg pb (by decide) (by decide)).trans hy.pb_eq) + ((r₄ sz (by decide)).trans xsz) ((r₄ best (by decide)).trans xbest) hB + (by simp +decide [x₄, x₃, x₂, U]) (by simp +decide [x₄, U]) hf hfitC hx xheap hvis + -- The invariant after the round. + have zreg (i : Fin 27) (h₁ : i ∉ [j, tv, tc, cand, dv, t, k, p, vk, vp, dk, dp, sz]) + (h₂ : i ∉ [t, j, jend]) (h₃ : i ∉ [t, u, best, vk, sz, k, p, vp, dk, dp, l, szp]) : + z.regs i = y.regs i := (zr i h₁).trans (reg i h₂ h₃) + have hvz : G.visitedOf z = insert U (G.visitedOf y) := by + ext v + simp only [ListInstance.visitedOf, Finset.mem_filter, Finset.mem_univ, true_and, + Finset.mem_insert] + rw [zvis v, hxD] + change (_ ∧ x.mem _ = 0) ↔ _ + rw [xpos v v.isLt] + by_cases hv : (v : ℕ) = U0 + · have hvU : v = U := Fin.ext hv + have hd : y.mem (G.dbase + v) ≠ 0 := by rw [hv, hB0]; exact hB + exact ⟨fun _ => Or.inl hvU, fun _ => ⟨hd, Or.inr hv⟩⟩ + · have hvU : v ≠ U := fun h => hv (by rw [h]) + simp [hv, hvU] + have hdz : G.distOf z = G.toInstance.relax (G.distOf y) U := by + funext v + change decodeDist (z.mem (G.dbase + v)) = _ + rw [zD v] + exact decode_edgesWord G hB hdu (by rw [hxD v]; rfl) + have hz : G.OuterInv z := by + refine ⟨?_, ?_, ?_, ?_, ?_, ?_, ?_, ?_, ?_, zheap, ?_⟩ + · exact (zreg zero (by decide) (by decide) (by decide)).trans hy.zero_eq + · exact (zreg one (by decide) (by decide) (by decide)).trans hy.one_eq + · exact (zreg two (by decide) (by decide) (by decide)).trans hy.two_eq + · exact (zreg three (by decide) (by decide) (by decide)).trans hy.three_eq + · exact (zreg eb (by decide) (by decide) (by decide)).trans hy.eb_eq + · exact (zreg db (by decide) (by decide) (by decide)).trans hy.db_eq + · exact (zreg hb (by decide) (by decide) (by decide)).trans hy.hb_eq + · exact (zreg pb (by decide) (by decide) (by decide)).trans hy.pb_eq + · intro a ha + rw [zin a ha] + exact hx a ha + · rw [hvz, hdz] + exact hstep + -- The potential pays for the round. + have hrem : G.remaining y = G.remaining z + (G.adj U).length := by + have hs : Finset.univ.filter (· ∉ G.visitedOf y) = + insert U (Finset.univ.filter (· ∉ insert U (G.visitedOf y))) := by + ext v + by_cases hv : v = U + · subst hv + simp [hu] + · simp [hv] + unfold ListInstance.remaining + rw [hvz, hs, Finset.sum_insert (by simp), add_comm] + have hlog : Nat.log 2 (y.regs sz) ≤ Nat.log 2 G.n := Nat.log_mono_right hy.heap.size_le + refine ⟨z, _, ex.seq (e₁.seq (e₂.seq (e₃.seq (e₄.seq ez)))), hz, ?_⟩ + unfold ListInstance.potential + rw [hrem] + generalize Nat.log 2 G.n = L at hlog ⊢ + generalize (G.adj U).length = d at zsz ⊢ + generalize G.remaining z = R + have h₁ : 100 * (L + 1) * (2 * R + z.regs sz) + 100 * (L + 1) * (d + 1) ≤ + 100 * (L + 1) * (2 * (R + d) + y.regs sz) := by + rw [← Nat.mul_add] + exact Nat.mul_le_mul_left _ (by omega) + have h₂ : d * (30 * L + 32) ≤ 100 * (L + 1) * d := by nlinarith + have h₃ : 100 * (L + 1) * (d + 1) = 100 * (L + 1) * d + 100 * L + 100 := by ring + omega + +/-- When the rounds stop, with an empty heap, the output cells hold the encoded distances. -/ +theorem heapExit_spec {y : Data 27} (hy : G.OuterInv y) (hsz : ¬(Cond.lt zero sz).Holds y) + (v : Fin G.n) : y.mem (G.outAddr v) = encodeDist (G.dist G.source v) := by + have hs : y.regs sz = 0 := by + simp only [Cond.holds_lt, hy.zero_eq] at hsz + omega + have hd : ∀ v ∉ G.visitedOf y, G.distOf y v = ⊤ := by + intro v hv + have hp : y.mem (G.dbase + 2 * G.n + 1 + v) = 0 := by + have := hy.heap.pos_le v v.isLt + omega + have h0 : y.mem (G.dbase + v) = 0 := by + simpa [ListInstance.visitedOf, hp] using hv + simp [ListInstance.distOf, h0, decodeDist] + rw [show G.dist G.source v = G.distOf y v from + (Instance.Inv.eq_dist G.toInstance hy.inv hd v).symm] + exact (encodeDist_decodeDist _).symm + +/-- **Correctness and running time of the structured program.** -/ +theorem heapDijkstra_spec (hw : Fits w (8 * (G.n + G.m + 2) * (G.maxWeight + 2))) : + ∃ y, ExecWithin w heapDijkstra (.init G.encode) y + (200 * (G.m + 1) * (Nat.log 2 G.n + 1) + 27) ∧ + ∀ v : Fin G.n, y.mem (G.outAddr v) = encodeDist (G.dist G.source v) := by + obtain ⟨x, hx, hix, hΦ⟩ := heapSetup_spec hw + obtain ⟨y, he, hi, hc⟩ := ExecWithin.loop_potential G.OuterInv G.potential + (fun _ h hc => heapRound_spec hw h hc) hix + refine ⟨y, (hx.seq he).mono ?_, fun v => heapExit_spec hi hc v⟩ + nlinarith + +/-- **Dijkstra's algorithm with a binary heap in `O((m + 1) (⌊log₂ n⌋ + 1))` time.** The +compiled program meets the specification `Cslib.Algorithms.WordRAM.Dijkstra.LogLinearSpec` with +`K = 230`. -/ +theorem logLinearSpec : LogLinearSpec (heapDijkstra.program zero) 230 := by + intro G w hw + obtain ⟨y, ⟨m, hm, he⟩, hd⟩ := heapDijkstra_spec hw + have hr := he.runs (pc := 0) (p := heapDijkstra.program zero) + (Program.Contains.self _).append_left + refine ⟨⟨heapDijkstra.size, y⟩, m, hr, ⟨y.regs zero, ?_⟩, ?_, hd⟩ + · exact step_halt (by simp [Code.program]) + · have : 1 ≤ (G.m + 1) * (Nat.log 2 G.n + 1) := Nat.one_le_iff_ne_zero.mpr (by positivity) + nlinarith + +/-- **Log-linear main theorem.** Some word-RAM program over the exact basis computes all +shortest-path distances from the source of a weighted directed graph with `n` vertices and `m` +edges, given by adjacency lists, within `O((m + 1) (⌊log₂ n⌋ + 1))` instructions. -/ +theorem logLinearStatement : LogLinearStatement := + ⟨27, _, 230, logLinearSpec⟩ + +end Cslib.Algorithms.WordRAM.Dijkstra + +end diff --git a/Complexity/Algorithms/WordRAM/Dijkstra/Statement.lean b/Complexity/Algorithms/WordRAM/Dijkstra/Statement.lean new file mode 100644 index 0000000000..038ebfbd74 --- /dev/null +++ b/Complexity/Algorithms/WordRAM/Dijkstra/Statement.lean @@ -0,0 +1,94 @@ +/- +Copyright (c) 2026 Renfei Zhou. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: Renfei Zhou +-/ + +module +-- upstream: candidate + +public import Complexity.Init +public import Mathlib.Combinatorics.Quiver.Path.Weight +public import Mathlib.Data.ENat.Lattice +public import Mathlib.Data.Fintype.Prod + +/-! +# Dijkstra's algorithm on the word RAM: statements + +This file defines the single-source shortest-path problem solved by Dijkstra's algorithm +[Dijkstra][Dijkstra1959] [Cormen et al., Section 24.3][CLRS2009]. + +An instance (`Cslib.Algorithms.WordRAM.Dijkstra.Instance`) is a directed graph on `Fin n` with +natural-number edge weights, given by its weight matrix, and a source vertex. The distance +`Cslib.Algorithms.WordRAM.Dijkstra.Instance.dist` is the infimum, in `ℕ∞`, of the weights of the +paths (`Quiver.Path`, `Quiver.Path.addWeight`) between two vertices; it is `⊤` when there is no +path. A program writes a distance as a word with `0` for an unreachable vertex and `d + 1` for +distance `d` (`Cslib.Algorithms.WordRAM.Dijkstra.encodeDist`). + +## Main definitions + +* `Cslib.Algorithms.WordRAM.Dijkstra.Instance`: a weighted directed graph with a source. +* `Cslib.Algorithms.WordRAM.Dijkstra.Instance.dist`: the shortest-path distance. +* `Cslib.Algorithms.WordRAM.Dijkstra.encodeDist`: a distance as a word. + +## Main results + +* `Cslib.Algorithms.WordRAM.Dijkstra.Instance.dist_self`: every vertex is at distance `0` from + itself. + +## References + +* [E. W. Dijkstra, *A note on two problems in connexion with graphs*][Dijkstra1959] +* [T. H. Cormen, C. E. Leiserson, R. L. Rivest and C. Stein, *Introduction to + Algorithms*][CLRS2009] +-/ + +@[expose] public section + +namespace Cslib.Algorithms.WordRAM.Dijkstra + +/-- An instance of the **single-source shortest-path problem**: a directed graph on `Fin n` with +natural-number edge weights, given by its weight matrix, and a source vertex. The entry +`weight u v` is `some c` for an edge from `u` to `v` of weight `c`, and `none` for no edge. -/ +structure Instance where + /-- The number of vertices. -/ + n : ℕ + /-- The weight matrix. -/ + weight : Fin n → Fin n → Option ℕ + /-- The source. -/ + source : Fin n + +namespace Instance + +variable (I : Instance) + +/-- The vertices of an instance, as a quiver: an arrow from `u` to `v` is a weight `c` with +`I.weight u v = some c`, so there is at most one. -/ +def Vertex : Type := Fin I.n + +instance : Quiver I.Vertex where + Hom u v := {c : ℕ // I.weight u v = some c} + +/-- The weight of an arrow. -/ +def arrowWeight {u v : I.Vertex} (e : u ⟶ v) : ℕ∞ := + (show {c : ℕ // I.weight u v = some c} from e).1 + +/-- The **shortest-path distance** from `u` to `v`: the infimum of the weights of the paths from +`u` to `v`, and `⊤` when there is none. -/ +noncomputable def dist (u v : Fin I.n) : ℕ∞ := + ⨅ p : Quiver.Path (V := I.Vertex) u v, p.addWeight I.arrowWeight + +/-- Every vertex is at distance `0` from itself, by the empty path. -/ +theorem dist_self (u : Fin I.n) : I.dist u u = 0 := + le_antisymm (iInf_le_of_le .nil (Quiver.Path.addWeight_nil _ _).le) (by simp) + +/-- The largest edge weight, and `0` when there is no edge. -/ +def maxWeight : ℕ := + Finset.univ.sup fun uv : Fin I.n × Fin I.n => (I.weight uv.1 uv.2).getD 0 + +end Instance + +/-- A distance as a word: `0` for an unreachable vertex and `d + 1` for distance `d`. -/ +def encodeDist (d : ℕ∞) : ℕ := if d = ⊤ then 0 else d.toNat + 1 + +end Cslib.Algorithms.WordRAM.Dijkstra diff --git a/Complexity/Computability/WordRAM/Notation.lean b/Complexity/Computability/WordRAM/Notation.lean index ff6a00d290..722c3e240b 100644 --- a/Complexity/Computability/WordRAM/Notation.lean +++ b/Complexity/Computability/WordRAM/Notation.lean @@ -30,10 +30,12 @@ for instance local abbreviations of elements of `Fin r`. The statements are: The notation unfolds to the constructors of `Cslib.WordRAM.Code`, so it adds nothing to the semantics. The binary operations and `width%` are instructions of the standard basis (`Cslib.WordRAM.Op.bin`, `Cslib.WordRAM.Op.width`); the other statements work over every basis. -It is scoped: `open Cslib.WordRAM` enables it. Its statement forms follow the word-RAM notation of -[Algolean][Algolean2026]. Its only new tokens are `ram[`, `width%` and the -internal `ram_stmt%` and `ram_block%`, so it reserves no identifiers: `ram` and `width` remain -ordinary names. +It is scoped: `open Cslib.WordRAM` enables it, with its statement forms, their expansions and the +tactic `ram_eval`. Only the category `wramStmt` and the syntax abbreviation `wramBlock` are +global, since Lean cannot scope them; without the open, `wramStmt` has no rules, so they parse +nothing. The statement forms follow the word-RAM notation of [Algolean][Algolean2026]. Its only +new tokens are `ram[`, `width%` and the internal `ram_stmt%` and `ram_block%`, so it reserves no +identifiers: `ram` and `width` remain ordinary names. ## Main definitions @@ -41,14 +43,14 @@ ordinary names. ## Tactics -* `ram_eval`, `ram_eval n`: symbolic execution. On a goal `Cslib.WordRAM.ExecWithin w c x y T` - or `Cslib.WordRAM.Exec w c x y T`, it runs the interpreter `Cslib.WordRAM.Code.eval` on `c` and - `x` with `simp`, allowing each loop `n` iterations (default `0`, for loop-free code), and - closes the goal when the result is `y` within (exactly, for `Exec`) `T` instructions; with `y` - a metavariable, it determines `y`. Which constants and results fit a word is decided for - numerals and otherwise taken from the given simp lemmas, stated with `Cslib.WordRAM.Fits`. - Programs given by name are unfolded by passing their definitions among the simp lemmas, as in - `ram_eval [myProgram, h]`. +* `ram_eval`, `ram_eval n` (scoped, like the notation): symbolic execution. On a goal + `Cslib.WordRAM.ExecWithin w c x y T` or `Cslib.WordRAM.Exec w c x y T`, it runs the interpreter + `Cslib.WordRAM.Code.eval` on `c` and `x` with `simp`, allowing each loop `n` iterations (default + `0`, for loop-free code), and closes the goal when the result is `y` within (exactly, for `Exec`) + `T` instructions; with `y` a metavariable, it determines `y`. Which constants and results fit a + word is decided for numerals and otherwise taken from the given simp lemmas, stated with + `Cslib.WordRAM.Fits`. Programs given by name are unfolded by passing their definitions among the + simp lemmas, as in `ram_eval [myProgram, h]`. ## References @@ -75,54 +77,54 @@ scoped syntax "ram_block%" wramBlock : term scoped syntax "ram_stmt%" wramStmt : term /-- Set a register to a constant. -/ -syntax ident " ← " "#" term:max : wramStmt +scoped syntax ident " ← " "#" term:max : wramStmt /-- Copy a register. -/ -syntax ident " ← " ident : wramStmt +scoped syntax ident " ← " ident : wramStmt /-- Addition. -/ -syntax ident " ← " ident " + " ident : wramStmt +scoped syntax ident " ← " ident " + " ident : wramStmt /-- Truncated subtraction. -/ -syntax ident " ← " ident " - " ident : wramStmt +scoped syntax ident " ← " ident " - " ident : wramStmt /-- Multiplication. -/ -syntax ident " ← " ident " * " ident : wramStmt +scoped syntax ident " ← " ident " * " ident : wramStmt /-- Division. -/ -syntax ident " ← " ident " / " ident : wramStmt +scoped syntax ident " ← " ident " / " ident : wramStmt /-- Remainder. -/ -syntax ident " ← " ident " % " ident : wramStmt +scoped syntax ident " ← " ident " % " ident : wramStmt /-- Bitwise and. -/ -syntax ident " ← " ident " &&& " ident : wramStmt +scoped syntax ident " ← " ident " &&& " ident : wramStmt /-- Bitwise or. -/ -syntax ident " ← " ident " ||| " ident : wramStmt +scoped syntax ident " ← " ident " ||| " ident : wramStmt /-- Bitwise exclusive or. -/ -syntax ident " ← " ident " ^^^ " ident : wramStmt +scoped syntax ident " ← " ident " ^^^ " ident : wramStmt /-- Shift left. -/ -syntax ident " ← " ident " <<< " ident : wramStmt +scoped syntax ident " ← " ident " <<< " ident : wramStmt /-- Shift right. -/ -syntax ident " ← " ident " >>> " ident : wramStmt +scoped syntax ident " ← " ident " >>> " ident : wramStmt /-- Load from the memory cell whose address is in a register. -/ -syntax ident " ← " "*" ident : wramStmt +scoped syntax ident " ← " "*" ident : wramStmt /-- Store to the memory cell whose address is in a register. -/ -syntax "*" ident " ← " ident : wramStmt +scoped syntax "*" ident " ← " ident : wramStmt /-- Read the word size. -/ -syntax ident " ← " "width%" : wramStmt +scoped syntax ident " ← " "width%" : wramStmt /-- Sample uniformly from `[0, regs s]` into a register. -/ -syntax ident " ← " &"rand" ppSpace ident : wramStmt +scoped syntax ident " ← " &"rand" ppSpace ident : wramStmt /-- A conditional on the equality of two registers. -/ -syntax "if " ident " = " ident " then " wramBlock " else " wramBlock : wramStmt +scoped syntax "if " ident " = " ident " then " wramBlock " else " wramBlock : wramStmt /-- A conditional on the comparison of two registers. -/ -syntax "if " ident " < " ident " then " wramBlock " else " wramBlock : wramStmt +scoped syntax "if " ident " < " ident " then " wramBlock " else " wramBlock : wramStmt /-- A loop while two registers are equal. -/ -syntax "while " ident " = " ident " do " wramBlock : wramStmt +scoped syntax "while " ident " = " ident " do " wramBlock : wramStmt /-- A loop while one register is smaller than another. -/ -syntax "while " ident " < " ident " do " wramBlock : wramStmt +scoped syntax "while " ident " < " ident " do " wramBlock : wramStmt /-- A structured program given as a term. -/ -syntax "⟪" term "⟫" : wramStmt +scoped syntax "⟪" term "⟫" : wramStmt -macro_rules +scoped macro_rules | `(ram[ $s:wramStmt ]) => `(ram_stmt% $s) | `(ram[ $s:wramStmt; $ss;* ]) => `(Cslib.WordRAM.Code.seq (ram_stmt% $s) (ram[ $ss;* ])) | `(ram_block% [ $ss;* ]) => `(ram[ $ss;* ]) -macro_rules +scoped macro_rules | `(ram_stmt% $d:ident ← #$v) => `(Cslib.WordRAM.Code.op (.const $d $v)) | `(ram_stmt% $d:ident ← $s:ident) => `(Cslib.WordRAM.Code.op (.copy $d $s)) | `(ram_stmt% $d:ident ← $i:ident + $j:ident) => @@ -162,10 +164,10 @@ macro_rules /-- Symbolically execute structured code on a goal `ExecWithin w c x y T` or `Exec w c x y T`, allowing each loop the given number of iterations (default `0`), and deciding with the given simp lemmas which written numbers fit a word. -/ -syntax (name := ramEval) "ram_eval" (ppSpace num)? (" [" Lean.Parser.Tactic.simpArg,* "]")? : +scoped syntax (name := ramEval) "ram_eval" (ppSpace num)? (" [" Lean.Parser.Tactic.simpArg,* "]")? : tactic -macro_rules +scoped macro_rules | `(tactic| ram_eval) => `(tactic| ram_eval 0 []) | `(tactic| ram_eval [$args,*]) => `(tactic| ram_eval 0 [$args,*]) | `(tactic| ram_eval $n:num) => `(tactic| ram_eval $n []) diff --git a/Complexity/Computability/WordRAM/Structured.lean b/Complexity/Computability/WordRAM/Structured.lean index a5dbf60ea0..a83d3abfb5 100644 --- a/Complexity/Computability/WordRAM/Structured.lean +++ b/Complexity/Computability/WordRAM/Structured.lean @@ -62,6 +62,8 @@ injectively chosen registers of a larger register file and leaves the others unt * `Cslib.WordRAM.Exec.unique`: executions are deterministic. * `Cslib.WordRAM.Exec.of_eval`: the interpreter only finds executions. * `Cslib.WordRAM.ExecWithin.loop_invariant`: the loop rule with an invariant and a measure. +* `Cslib.WordRAM.ExecWithin.loop_potential`: the loop rule with an invariant and a potential + that pays for iterations of different times. * `Cslib.WordRAM.Exec.regs_eq`: registers that a program never writes keep their contents. * `Cslib.WordRAM.Exec.rename`: renamed programs run on the chosen registers. * `Cslib.WordRAM.Exec.mono_width`, `Cslib.WordRAM.Exec.bounded`: a program that does not read the @@ -387,6 +389,26 @@ theorem loop_invariant (I : Data r → Prop) (μ : Data r → ℕ) {T : ℕ} omega · exact ⟨x, (loop_exit hc).mono (by omega), hx, hc⟩ +/-- **Loop rule, invariant and potential.** If every state `x` satisfying the invariant `I` and +the condition is led by the body to a state `y` satisfying `I`, within a number `T` of instructions +with `T + 2 + Φ y ≤ Φ x`, then from a state satisfying `I` the loop terminates within `Φ x + 2` +instructions, in a state satisfying `I` where the condition fails. The two extra instructions of +each iteration are the condition test and the jump back; the time of the body may differ between +iterations, as long as the potential `Φ` pays for it. -/ +theorem loop_potential (I : Data r → Prop) (Φ : Data r → ℕ) + (hbody : ∀ x, I x → cond.Holds x → + ∃ y T, ExecWithin w body x y T ∧ I y ∧ T + 2 + Φ y ≤ Φ x) + (hx : I x) : + ∃ y, ExecWithin w (.loop cond body) x y (Φ x + 2) ∧ I y ∧ ¬cond.Holds y := by + induction hm : Φ x using Nat.strong_induction_on generalizing x with + | _ m ih => + subst hm + by_cases hc : cond.Holds x + · obtain ⟨y, T, hy, hIy, hΦ⟩ := hbody x hx hc + obtain ⟨z, hz, hIz, hcz⟩ := ih (Φ y) (by omega) hIy rfl + exact ⟨z, (loop_iter hc hy hz).mono (by omega), hIz, hcz⟩ + · exact ⟨x, (loop_exit hc).mono (by omega), hx, hc⟩ + /-! ### Data instructions -/ /-- Setting register `d` to a fitting constant `v` executes within one instruction. -/ diff --git a/ComplexityTests/WordRAMDijkstra.lean b/ComplexityTests/WordRAMDijkstra.lean new file mode 100644 index 0000000000..5263307a61 --- /dev/null +++ b/ComplexityTests/WordRAMDijkstra.lean @@ -0,0 +1,73 @@ +/- +Copyright (c) 2026 Renfei Zhou. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: Renfei Zhou +-/ + +module +-- upstream: candidate + +public import Complexity.Algorithms.WordRAM.Dijkstra.LogLinear + +/-! Nonvacuity and interpreter checks for Dijkstra's algorithm with a binary heap. -/ + +@[expose] public section + +open Cslib Cslib.WordRAM Cslib.Algorithms.WordRAM.Dijkstra + +namespace WordRAMDijkstraTests + +/-- Three vertices, a parallel edge and a shorter path through a third vertex: edges `0 → 1` of +weights `4` and `5`, `0 → 2` of weight `1`, and `2 → 1` of weight `2`, from source `0`. The +distances are `0`, `3` and `1`. -/ +def graph : ListInstance where + n := 3 + adj := ![[(1, 4), (2, 1), (1, 5)], [], [(1, 2)]] + source := 0 + +example : graph.m = 4 ∧ graph.maxWeight = 5 := by decide + +/-- The input words: `n`, the source, the offsets `0, 3, 3, 4`, then the edges. -/ +example : graph.encode = [3, 0, 0, 3, 3, 4, 1, 4, 2, 1, 1, 5, 1, 2] := by decide + +/-- The word-size hypothesis of `LogLinearSpec` holds for the instance at word size `10`. -/ +theorem fits : Fits 10 (8 * (graph.n + graph.m + 2) * (graph.maxWeight + 2)) := by decide + +/-- The main theorem applies to the instance: the compiled program halts within +`230 (m + 1) (⌊log₂ n⌋ + 1)` instructions and leaves `1` at the source's output cell, the +encoded distance `0`. -/ +example : ∃ t steps, Runs 10 (heapDijkstra.program HeapReg.zero) (.init graph.encode) t steps ∧ + (∃ v, step 10 (heapDijkstra.program HeapReg.zero) t = .halt v) ∧ steps + 1 ≤ 230 * 5 * 2 ∧ + t.data.mem (graph.outAddr graph.source) = 1 := by + obtain ⟨t, steps, hr, hh, hle, hd⟩ := logLinearSpec graph 10 fits + have hm : graph.m = 4 := by decide + have hlog : Nat.log 2 graph.n = 1 := Nat.log_eq_of_pow_le_of_lt_pow (by decide) (by decide) + rw [hm, hlog] at hle + refine ⟨t, steps, hr, hh, hle, ?_⟩ + rw [hd] + exact (congrArg encodeDist (graph.toInstance.dist_self graph.toInstance.source)).trans + (by simp [encodeDist]) + +/-- The interpreter runs the structured program on the instance and finds the encoded distances +`1`, `4` and `2` in the output cells, within the bound `200 (m + 1) (⌊log₂ n⌋ + 1) + 27` of +`heapDijkstra_spec`. -/ +example : ((heapDijkstra.eval 10 8 (.init graph.encode)).map fun yn => + ((List.finRange 3).map fun v => yn.1.mem (graph.outAddr v), decide (yn.2 ≤ 2027))) = + some ([1, 4, 2], true) := by + decide +kernel + +/-- Two vertices and no edge, from source `0`: vertex `1` is unreachable. -/ +def edgeless : ListInstance where + n := 2 + adj := ![[], []] + source := 0 + +/-- With no edge the bound `200 (m + 1) (⌊log₂ n⌋ + 1) + 27` of `heapDijkstra_spec` is `427`. The +interpreter finds `1` (distance `0`) in the source's output cell and `0` (unreachable) in the other +one, within that bound. -/ +example : ((heapDijkstra.eval 10 8 (.init edgeless.encode)).map fun yn => + ((List.finRange 2).map fun v => yn.1.mem (edgeless.outAddr v), decide (yn.2 ≤ 427))) = + some ([1, 0], true) := by + decide +kernel + +end WordRAMDijkstraTests diff --git a/abstracts/ComplexityAbstracts.lean b/abstracts/ComplexityAbstracts.lean index caf5fa44df..a6b4d5d743 100644 --- a/abstracts/ComplexityAbstracts.lean +++ b/abstracts/ComplexityAbstracts.lean @@ -12,6 +12,8 @@ import ComplexityAbstracts.CookLevinAbstract import ComplexityAbstracts.CookLevinBridge import ComplexityAbstracts.CryptographyAbstract import ComplexityAbstracts.CryptographyBridge +import ComplexityAbstracts.DijkstraAbstract +import ComplexityAbstracts.DijkstraBridge import ComplexityAbstracts.Encoding import ComplexityAbstracts.FaginAbstract import ComplexityAbstracts.FaginBridge diff --git a/abstracts/ComplexityAbstracts/DijkstraAbstract.lean b/abstracts/ComplexityAbstracts/DijkstraAbstract.lean new file mode 100644 index 0000000000..dca8da617e --- /dev/null +++ b/abstracts/ComplexityAbstracts/DijkstraAbstract.lean @@ -0,0 +1,286 @@ +/- +Copyright (c) 2026 Renfei Zhou. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: Renfei Zhou +-/ + +-- upstream: local-only (formal abstract library) + +import ComplexityAbstracts.Check +import ComplexityAbstracts.DijkstraBridge + +/-! # Dijkstra's algorithm: shortest paths in `O((m + 1) log n)` time + +Given a directed graph with `n` vertices, `m` edges and natural-number edge weights, and a source +vertex `s`, the single-source shortest-path problem asks for the distance from `s` to every +vertex. Dijkstra's algorithm with a binary heap solves it in `O((m + 1) log n)` time. The theorem +below states this bound for a program on the word RAM, with the graph given by adjacency lists. + +The graph, its distances, the input and output layouts, the word-size condition and the machine are +all defined here from scratch. The machine, defined in the appendix, is the word RAM over the exact +basis of operations: each operation computes its exact value, and an instruction faults instead of +producing `x - y` for `y > x`, `x / 0`, `x % 0`, or a value that is not a word. Apart from `Runs`, +its definitions are those of `WordRAMAbstract`. To check that the theorem says what is claimed, a +reader only needs to inspect the definitions and the theorem statement in this file. Lean checks the +proof through `DijkstraBridge`, which proves that the distance here is the library's and that the +machine here runs as the library's word RAM `Cslib.WordRAM` over its exact basis, and applies the +library theorem `Cslib.Algorithms.WordRAM.Dijkstra.logLinearStatement`. `#check_abstract` checks +mechanically that every definition and theorem statement here depends only on this file, Lean, +Mathlib and upstream CSLib. + +## References + +* [E. W. Dijkstra, *A note on two problems in connexion with graphs*][Dijkstra1959] +* [T. H. Cormen, C. E. Leiserson, R. L. Rivest and C. Stein, *Introduction to + Algorithms*][CLRS2009], Section 24.3. +-/ + +namespace ComplexityAbstracts.DijkstraAbstract + +variable {n : ℕ} + +/-! ## Graphs and distances + +The vertices are `0, …, n - 1` (the type `Fin n`). A walk is a sequence of edges, each leaving +the vertex the previous one enters; vertices and edges may repeat. Its weight is the sum of the +weights of its edges. The distance from `u` to `v` is the least weight of a walk from `u` to `v`, +in `ℕ∞ = ℕ ∪ {⊤}`, and `⊤` when there is no such walk. Weights are natural numbers, so a least +weight exists whenever a walk does. +-/ + +/-- A directed graph on the vertices `Fin n` with natural-number edge weights, given by adjacency +lists: `G u` holds a pair `(v, c)` for each edge from `u` to `v` of weight `c`. The same pair of +endpoints may appear more than once (parallel edges), and `v` may be `u` (self-loops). -/ +abbrev Graph (n : ℕ) := Fin n → List (Fin n × ℕ) + +/-- `IsWalk G u l v`: the list `l` of edges `(target, weight)` is a walk from `u` to `v`. Its +first edge is in the list `G u` of edges out of `u`, each later edge is in the list of edges out +of the target of the edge before it, and the target of its last edge is `v`. The empty walk goes +from `u` to `u`. -/ +def IsWalk (G : Graph n) : Fin n → List (Fin n × ℕ) → Fin n → Prop + | u, [], v => u = v + | u, e :: l, v => e ∈ G u ∧ IsWalk G e.1 l v + +/-- The **distance** from `u` to `v`: the least weight of a walk from `u` to `v`, the weight of +a walk being the sum of the weights of its edges, and `⊤` when there is no walk. -/ +noncomputable def dist (G : Graph n) (u v : Fin n) : ℕ∞ := + ⨅ (l : List (Fin n × ℕ)) (_ : IsWalk G u l v), ((l.map Prod.snd).sum : ℕ∞) + +/-- The number `m` of edges. -/ +def numEdges (G : Graph n) : ℕ := ∑ u, (G u).length + +/-- The largest edge weight, and `0` when there is no edge. -/ +def maxWeight (G : Graph n) : ℕ := + Finset.univ.sup fun u => ((G u).map Prod.snd).foldr max 0 + +/-! ## Input, output and cost + +The program reads the graph as a list of natural numbers: `n`, the source `s`, the `n + 1` offsets +of the adjacency lists in an array of edges, and then the `m` edges, each as two numbers, its target +and its weight. The input and its length header are placed in memory before the run at no cost; +memory beyond them and all registers start at `0`. When it halts, the memory cell `outAddr G v` of +each vertex `v` holds `0` when `v` is unreachable from `s` and `d + 1` when the distance from `s` to +`v` is `d`. The machine is a word RAM with `w`-bit words, and the theorem +asks `8 (n + m + 2) (W + 2) < 2 ^ w`, where `W` is the largest weight, so a word size +`w = O(log (n + m + W + 2))` suffices. Each instruction costs one unit of time, whatever its +operands. The operations are `+`, `-`, `*`, `/`, `%`, bitwise operations and shifts, each with its +exact value; an instruction faults on `x - y` with `y > x`, on `x / 0` and `x % 0`, and on a +constant or operation value of `2 ^ w` or more. The bound counts the instructions executed, the +final `halt` included. The appendix gives the exact layout (`offset`, `encode`, `outAddr`, +`encodeDist`) and the machine (`Program`, `init`, `step`, `Runs`). +-/ + +/-! ## Representation details + +The input is a list of natural numbers: `n`, the source, the `n + 1` offsets of the adjacency +lists in the edge array, and the edges, in the order of the lists `G 0, …, G (n - 1)`, each as +two numbers, its target and its weight. This is the usual "compressed" layout: the edges out of +`u` are the edges `j` with `offset G u ≤ j < offset G (u + 1)`. + +The machine stores the input of length `k = n + 2 m + 3` at addresses `1, …, k` (see `init` +below). The output is the `n` memory cells just past it: the cell of vertex `v` holds `0` when `v` +is unreachable from the source and `d + 1` when its distance is `d`. +-/ + +/-- The offset of the adjacency list of vertex `u` in the edge array: the number of edges out of +the vertices below `u`. The offset of `n` is the number of edges. -/ +def offset (G : Graph n) (u : ℕ) : ℕ := + (((List.finRange n).take u).map fun x => (G x).length).sum + +/-- The input: `n`, the source `s`, the offsets of the vertices `0, …, n`, and then the edges of +`G 0, …, G (n - 1)` in order, each as its target and its weight. -/ +def encode (G : Graph n) (s : Fin n) : List ℕ := + n :: (s : ℕ) :: + ((List.range (n + 1)).map (offset G) ++ + (List.finRange n).flatMap fun u => (G u).flatMap fun e => [(e.1 : ℕ), e.2]) + +/-- The output address of vertex `v`: the input occupies addresses `0, …, n + 2 m + 3` (its length +at `0`), and the outputs follow it. -/ +def outAddr (G : Graph n) (v : Fin n) : ℕ := n + 2 * numEdges G + 4 + v + +/-- A distance as an output word: `0` for `⊤` (unreachable) and `d + 1` for distance `d`. -/ +def encodeDist (d : ℕ∞) : ℕ := if d = ⊤ then 0 else d.toNat + 1 + +/-! **The machine.** The word RAM over the exact basis, as in `WordRAMAbstract`. At word size +`w`, a *word* is a number below `2 ^ w`. A machine has `r` registers `R[0]`, …, `R[r - 1]` and a +memory cell `M[a]` at every address `a : ℕ`; all hold natural numbers. The machine stores arbitrary +natural numbers, but an instruction faults when it would create a value that is undefined or not a +word. There are ten binary operations, and each computes its exact value. Subtraction `x - y` is +undefined for `y > x`, and division `x / y` and remainder `x % y` are undefined for `y = 0`; +otherwise division rounds down. The operations `and`, `or` and `xor` act bitwise, and `shl` and +`shr` shift `x` left and right by `y` bits, that is, give `x * 2 ^ y` and `x / 2 ^ y` rounded +down. In `AluOp.eval`, a case value `n` stands for `some n`. +-/ + +variable {r : ℕ} + +/-- The binary operations. Division rounds down, `and`, `or`, `xor` are bitwise, and `shl`, +`shr` shift left and right by `y` bits. -/ +inductive AluOp + | add | sub | mul | div | mod | and | or | xor | shl | shr + +/-- The exact value of a binary operation, before any check that it is a word, or `none` where it +is undefined: `x - y` for `y > x`, and `x / 0` and `x % 0`. -/ +def AluOp.eval : AluOp → ℕ → ℕ → Option ℕ + | .add, x, y => x + y + | .sub, x, y => if y ≤ x then x - y else none + | .mul, x, y => x * y + | .div, x, y => if y = 0 then none else x / y + | .mod, x, y => if y = 0 then none else x % y + | .and, x, y => x &&& y + | .or, x, y => x ||| y + | .xor, x, y => x ^^^ y + | .shl, x, y => some (x * 2 ^ y) + | .shr, x, y => some (x / 2 ^ y) + +/-! The number `r` of registers is a parameter of the types below (the implicit variable `r`). +There are ten instructions; the comment on each gives its effect. Register arguments `d`, `i`, +`j`, `a` range over `Fin r`, that is, `0, …, r - 1` (in `load` and `store`, register `a` holds +the address); a constant `x` and a jump target `t` are arbitrary natural numbers. The +instruction `width` reads the word size `w`, so a program can depend on it. +-/ + +/-- The instructions over `r` registers. Jump targets are instruction numbers. -/ +inductive Instr (r : ℕ) + | const (d : Fin r) (x : ℕ) -- `R[d] := x` + | copy (d i : Fin r) -- `R[d] := R[i]` + | alu (f : AluOp) (d i j : Fin r) -- `R[d] := f(R[i], R[j])` + | width (d : Fin r) -- `R[d] := w`, the word size + | load (d a : Fin r) -- `R[d] := M[R[a]]` + | store (a i : Fin r) -- `M[R[a]] := R[i]` + | jump (t : ℕ) -- go to instruction `t` + | branchEq (i j : Fin r) (t : ℕ) -- go to `t` if `R[i] = R[j]`, else to the next one + | branchLt (i j : Fin r) (t : ℕ) -- go to `t` if `R[i] < R[j]`, else to the next one + | halt (i : Fin r) -- stop with output `R[i]` + +/-- A program is a list of instructions, numbered from `0`. -/ +abbrev Program (r : ℕ) := List (Instr r) + +/-! The memory is a function on all of `ℕ`: there is no allocation and no bound on addresses. One +step executes instruction number `pc`. `const` and `alu` write their value to `R[d]` if it is a +word, and fault otherwise; there is no wrap-around modulo `2 ^ w`, and an `alu` whose value is +undefined (`none`) also faults. `width` writes `w`, which is always a word, since `w < 2 ^ w`. +`copy`, `load` and `store` move a value without a check. `jump` and the two branches set `pc`; +every other instruction except `halt` moves on to instruction `pc + 1`. `halt i` stops with +output `R[i]`. If `pc` is not the number of an instruction, that is, `pc` is at least the length +of the program, the step faults. + +Copies, loads and stores need no check: when the input and its length are words, an induction +along the run shows that every register and memory cell always holds a word, and so every address +used is below `2 ^ w`. This induction is not stated here; it only explains why the check is not +needed. Under the word-size condition of the theorem, the input and its length are words. +-/ + +/-- A state: the number `pc` of the next instruction, the registers and the memory. -/ +structure State (r : ℕ) where + /-- The number of the next instruction. -/ + pc : ℕ + /-- `reg i` is `R[i]`. -/ + reg : Fin r → ℕ + /-- `mem a` is `M[a]`. -/ + mem : ℕ → ℕ + +/-- One step continues in a new state, halts with an output, or faults. -/ +inductive Outcome (r : ℕ) + | next (s : State r) + | halt (v : ℕ) + | fault + +/-- One step of program `p` at word size `w`. A constant or operation whose value is not a word +faults (no wrap-around), and so do an operation whose value is undefined and running off the end +of the program. Copies, loads and stores are not checked. -/ +def step (w : ℕ) (p : Program r) (s : State r) : Outcome r := + let set (d : Fin r) (x : ℕ) : Outcome r := -- `R[d] := x`, then the next instruction + .next { s with pc := s.pc + 1, reg := Function.update s.reg d x } + let setWord (d : Fin r) (x : ℕ) : Outcome r := if x < 2 ^ w then set d x else .fault + match p[s.pc]? with + | none => .fault + | some (.const d x) => setWord d x + | some (.copy d i) => set d (s.reg i) + | some (.alu f d i j) => + match f.eval (s.reg i) (s.reg j) with + | some x => setWord d x + | none => .fault + | some (.width d) => set d w -- `w < 2 ^ w`, so this never faults + | some (.load d a) => set d (s.mem (s.reg a)) + | some (.store a i) => + .next { s with pc := s.pc + 1, mem := Function.update s.mem (s.reg a) (s.reg i) } + | some (.jump t) => .next { s with pc := t } + | some (.branchEq i j t) => .next { s with pc := if s.reg i = s.reg j then t else s.pc + 1 } + | some (.branchLt i j t) => .next { s with pc := if s.reg i < s.reg j then t else s.pc + 1 } + | some (.halt i) => .halt (s.reg i) + +/-! A program reads the length `k` of its input `x₁, …, xₖ` at `M[0]` and the input at +`M[1], …, M[k]`. A run of `k` steps follows `step` one instruction at a time; each of its steps +continues to a next state. The theorem asks for a run to a state whose next step is a `halt`, so +the run executes one more instruction, the `halt`, and then stops. Every instruction takes one +unit of time, whatever its operands and the word size. +-/ + +/-- The start state on input `x₁, …, xₖ`: instruction `0`, registers `0`, and memory +`k, x₁, …, xₖ` at addresses `0, 1, …, k`, then `0` everywhere. -/ +def init (input : List ℕ) : State r := + ⟨0, fun _ => 0, fun a => (input.length :: input).getD a 0⟩ + +/-- `Runs w p s t k`: from state `s`, `k` steps of program `p` at word size `w` lead to state +`t`, each continuing to a next state, none halting or faulting. (`Relation.RelatesInSteps R s t k` +holds when `k` steps of the relation `R` lead from `s` to `t`: `s = t` for `k = 0`, and for +`k + 1`, `k` steps to some `u` and then `R u t`.) -/ +def Runs (w : ℕ) (p : Program r) : State r → State r → ℕ → Prop := + Relation.RelatesInSteps fun s t => step w p s = .next t + +/-! ## Main theorem -/ + +/-- **Dijkstra's algorithm.** Some program `p` of the word RAM defined in the representation +details (over the exact basis), with some number `r` of registers, and some constant `K` satisfy +the following. +For every graph `G` with `n` vertices and `m` edges and largest weight `W`, every source `s`, and +every word size `w` with `8 (n + m + 2) (W + 2) < 2 ^ w`: started on the input `encode G s`, the +program runs for `steps` steps to a state `t` whose next instruction is `halt`, executes +`steps + 1 ≤ K (m + 1) (Nat.log 2 n + 1)` instructions in all, and then holds at `outAddr G v` +the encoded distance from `s` to `v`, for every vertex `v`. Here `Nat.log 2 n` is the integer part +of `log₂ n` (and `0` for `n = 0`). The output of `halt` is not part of the result. -/ +theorem dijkstra_logLinear : + ∃ r, ∃ p : Program r, ∃ K, ∀ (n : ℕ) (G : Graph n) (s : Fin n) (w : ℕ), + 8 * (n + numEdges G + 2) * (maxWeight G + 2) < 2 ^ w → + ∃ t steps, Runs w p (init (encode G s)) t steps ∧ + (∃ v, step w p t = .halt v) ∧ + steps + 1 ≤ K * (numEdges G + 1) * (Nat.log 2 n + 1) ∧ + ∀ v : Fin n, t.mem (outAddr G v) = encodeDist (dist G s v) := + -- the bridge's halting runs and agreeing outcomes are unused here: `OutcomeAgrees` is `True` + DijkstraBridge.logLinear + ⟨AluOp, .add, .sub, .mul, .div, .mod, .and, .or, .xor, .shl, .shr, @AluOp.rec, AluOp.eval, + Instr, .const, .copy, .alu, .width, .load, .store, .jump, .branchEq, .branchLt, .halt, + @Instr.rec, State, .mk, Outcome, .next, .halt, .fault, step, + DijkstraBridge.haltsVia Instr State Outcome step fun f g a o => o.rec f g a, fun _ _ => True⟩ + (by constructor <;> simp_all [step, DijkstraBridge.haltsVia, AluOp.eval]) + State.mem (fun _ _ _ => rfl) (W := @IsWalk) ⟨fun _ _ _ => Iff.rfl, fun _ _ _ _ _ => Iff.rfl⟩ + +/-- +info: abstract audit passed for `ComplexityAbstracts.DijkstraAbstract`: 218 root +declarations, 229 dependencies checked +-/ +#guard_msgs (whitespace := lax) in +#check_abstract ComplexityAbstracts.DijkstraAbstract + +end ComplexityAbstracts.DijkstraAbstract diff --git a/abstracts/ComplexityAbstracts/DijkstraBridge.lean b/abstracts/ComplexityAbstracts/DijkstraBridge.lean new file mode 100644 index 0000000000..aba4d540af --- /dev/null +++ b/abstracts/ComplexityAbstracts/DijkstraBridge.lean @@ -0,0 +1,241 @@ +/- +Copyright (c) 2026 Renfei Zhou. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: Renfei Zhou +-/ + +-- upstream: local-only (formal abstract library) + +import ComplexityAbstracts.WordRAMBridge +import Complexity.Algorithms.WordRAM.Dijkstra.LogLinear + +/-! +# Bridge for the Dijkstra formal abstract + +`DijkstraAbstract.lean` defines the distance of a graph given by adjacency lists as the least +weight of a walk, a list of edges each of which leaves the target of the one before. The library +defines it (`Cslib.Algorithms.WordRAM.Dijkstra.ListInstance.dist`) as the least weight of a +`Quiver.Path` in the weight-matrix instance that keeps the lightest of the parallel edges. This +file proves that the two agree, for every walk predicate `W` that satisfies `WalkLaws`. + +The abstract also defines its own word RAM over the exact basis, the machine of +`WordRAMAbstract`. `WordRAMBridge` proves, for every presentation of that machine that satisfies +its `Laws`, that one step of a library program over `Cslib.WordRAM.Basis.exact` at an instruction +other than `rand` matches one step of its translation `WordRAMBridge.Presentation.programOf` +(`WordRAMBridge.Presentation.step_matches_at`). A step that does not fault is not at a `rand` +(`isRandFree_of_step`), so this file extends that to runs of a given number of steps (`runs_map`), +and transfers +`Cslib.Algorithms.WordRAM.Dijkstra.logLinearStatement` to the walk distance and the abstract's +machine. + +## Main results + +* `walkDist_eq`: the walk distance is the library's distance. +* `runs_map`: a run of a library program is a run of its translation. +* `logLinear`: the library's main theorem, over unfolded adjacency lists, the walk distance and + any presentation of the exact word RAM. + +## Syntax in the imports + +The library's Dijkstra program is written in the word-RAM program syntax of +`Complexity/Computability/WordRAM/Notation.lean` and its exact-basis variant +`Complexity/Algorithms/WordRAM/Dijkstra/ExactNotation.lean`, so the import of the library theorem +brings both files in (`LogLinear` → `HeapRelax` → `HeapSift` → `BinaryHeap` → `HeapProgram` → +`ExactNotation` → `Notation`). All their syntax rules, `macro_rules` and tactics are scoped, to +`Cslib.WordRAM` and to `Cslib.Algorithms.WordRAM.Dijkstra`. Only the syntax category `wramStmt` and +the syntax abbreviation `wramBlock` are global, since Lean cannot scope them, and outside those +namespaces `wramStmt` has no rules. `DijkstraAbstract.lean` opens neither namespace, so none of +this syntax is active there. +-/ + +namespace ComplexityAbstracts.DijkstraBridge + +open Cslib.WordRAM Cslib.Algorithms.WordRAM.Dijkstra ComplexityTests.Abstracts + +/-- The defining equations of a walk predicate: `W G u l v` says that the list `l` of edges +`(target, weight)` is a walk from `u` to `v` in the graph with adjacency lists `G`. -/ +structure WalkLaws + (W : ∀ {n : ℕ}, (Fin n → List (Fin n × ℕ)) → Fin n → List (Fin n × ℕ) → Fin n → Prop) : + Prop where + nil : ∀ {n : ℕ} (G : Fin n → List (Fin n × ℕ)) (u v : Fin n), W G u [] v ↔ u = v + cons : ∀ {n : ℕ} (G : Fin n → List (Fin n × ℕ)) (u : Fin n) (e : Fin n × ℕ) + (l : List (Fin n × ℕ)) (v : Fin n), W G u (e :: l) v ↔ e ∈ G u ∧ W G e.1 l v + +variable {W : ∀ {n : ℕ}, (Fin n → List (Fin n × ℕ)) → Fin n → List (Fin n × ℕ) → Fin n → Prop} + +/-- The least weight of a walk from `u` to `v`, and `⊤` when there is none. -/ +noncomputable def walkDist {n : ℕ} (W : (Fin n → List (Fin n × ℕ)) → Fin n → + List (Fin n × ℕ) → Fin n → Prop) (G : Fin n → List (Fin n × ℕ)) (u v : Fin n) : ℕ∞ := + ⨅ (l : List (Fin n × ℕ)) (_ : W G u l v), ((l.map Prod.snd).sum : ℕ∞) + +namespace ListInstance + +variable (hW : WalkLaws W) (G : ListInstance) +include hW + +/-- A walk followed by an edge out of its end is a walk. -/ +theorem walk_append {u b : Fin G.n} {l : List (Fin G.n × ℕ)} (hl : W G.adj u l b) + {e : Fin G.n × ℕ} (he : e ∈ G.adj b) : W G.adj u (l ++ [e]) e.1 := by + induction l generalizing u with + | nil => + rw [hW.nil] at hl + subst hl + exact (hW.cons _ _ _ _ _).2 ⟨he, (hW.nil _ _ _).2 rfl⟩ + | cons e' l ih => + rw [hW.cons] at hl + exact (hW.cons _ _ _ _ _).2 ⟨hl.1, ih hl.2⟩ + +/-- Every path of the weight-matrix instance is a walk of the same weight. -/ +theorem exists_walk {u v : G.toInstance.Vertex} (p : Quiver.Path u v) : + ∃ l, W G.adj u l v ∧ + ((l.map Prod.snd).sum : ℕ∞) = p.addWeight G.toInstance.arrowWeight := by + induction p with + | nil => + refine ⟨[], (hW.nil _ _ _).2 rfl, ?_⟩ + rw [Quiver.Path.addWeight_nil] + simp + | @cons b c p e ih => + obtain ⟨l, hl, hw⟩ := ih + have he : ((show Fin G.n from c), e.1) ∈ G.adj b := (G.lightest_eq_some_iff.1 e.2).1 + refine ⟨l ++ [((show Fin G.n from c), e.1)], walk_append hW G hl he, ?_⟩ + rw [Quiver.Path.addWeight_cons, ← hw] + simp only [List.map_append, List.map_cons, List.map_nil, List.sum_append, List.sum_cons, + List.sum_nil, add_zero, Nat.cast_add] + rfl + +/-- Every walk gives a path of the weight-matrix instance of at most the same weight: each edge +is replaced by the lightest edge between its endpoints. -/ +theorem exists_path {v : G.toInstance.Vertex} (l : List (Fin G.n × ℕ)) {u : G.toInstance.Vertex} + (hl : W G.adj u l v) : + ∃ p : Quiver.Path u v, + p.addWeight G.toInstance.arrowWeight ≤ ((l.map Prod.snd).sum : ℕ∞) := by + induction l generalizing u with + | nil => + obtain rfl : u = v := (hW.nil _ _ _).1 hl + refine ⟨.nil, ?_⟩ + rw [Quiver.Path.addWeight_nil] + exact zero_le + | cons e l ih => + replace hl := (hW.cons _ _ _ _ _).1 hl + obtain ⟨q, hq⟩ := ih hl.2 + obtain ⟨c, hc⟩ : ∃ c, G.lightest u e.1 = some c := by + cases h : G.lightest u e.1 with + | none => exact absurd hl.1 (G.lightest_eq_none_iff.1 h e.2) + | some c => exact ⟨c, rfl⟩ + have hce : c ≤ e.2 := (G.lightest_eq_some_iff.1 hc).2 e.2 hl.1 + let x : G.toInstance.Vertex := e.1 + let a : u ⟶ x := ⟨c, hc⟩ + let q' : Quiver.Path x v := q + refine ⟨(Quiver.Path.nil.cons a).comp q', ?_⟩ + rw [Quiver.Path.addWeight_comp, Quiver.Path.addWeight_cons, Quiver.Path.addWeight_nil, + zero_add] + simp only [List.map_cons, List.sum_cons, Nat.cast_add] + exact add_le_add (show ((c : ℕ) : ℕ∞) ≤ e.2 by exact_mod_cast hce) hq + +end ListInstance + +/-- The walk distance is the library's distance. -/ +theorem walkDist_eq (hW : WalkLaws W) (G : ListInstance) (u v : Fin G.n) : + walkDist W G.adj u v = G.dist u v := by + apply le_antisymm + · refine le_iInf fun p => ?_ + obtain ⟨l, hl, hw⟩ := ListInstance.exists_walk hW G p + exact (iInf₂_le l hl).trans hw.le + · refine le_iInf₂ fun l hl => ?_ + obtain ⟨p, hp⟩ := ListInstance.exists_path hW G l hl + exact (iInf_le _ p).trans hp + +/-! ## The machine -/ + +/-- Halting runs defined from one step by cases on its outcome, as the abstract `Halts` of +`WordRAMAbstract` is: `cases f g a o` is `f t`, `g x` or `a` when `o` continues in state `t`, +halts with output `x` or faults. The Dijkstra abstract states no halting runs; this fills the +field `Halts` of its `WordRAMBridge.Presentation`. -/ +def haltsVia (Instr State Outcome : ℕ → Type) + (step : ∀ {r : ℕ}, ℕ → List (Instr r) → State r → Outcome r) + (cases : ∀ {r : ℕ}, (State r → Prop) → (ℕ → Prop) → Prop → Outcome r → Prop) {r : ℕ} + (w : ℕ) (p : List (Instr r)) : State r → ℕ → ℕ → Prop + | _, _, 0 => False + | s, v, n + 1 => + cases (fun t => haltsVia Instr State Outcome step cases w p t v n) (fun x => x = v ∧ n = 0) + False (step w p s) + +namespace Machine + +variable {S : WordRAMBridge.Presentation} (h : WordRAMBridge.Laws S) {r : ℕ} +include h + +/-- The abstract state with the next instruction, registers and memory of a library state. -/ +def stateOf (S : WordRAMBridge.Presentation) (t : Cslib.WordRAM.State r) : S.State r := + S.state t.pc t.data.regs t.data.mem + +omit h in +theorem stateOf_init (input : List ℕ) : + stateOf S (.init input : Cslib.WordRAM.State r) = S.init input := by + rw [← WordRAMBridge.Presentation.init_eq] + rfl + +omit h in +/-- A step that does not fault is not at a `rand` instruction, which faults deterministically. -/ +theorem isRandFree_of_step {w : ℕ} {q : Program .exact r} {t : Cslib.WordRAM.State r} + (ht : Cslib.WordRAM.step w q t ≠ .fault) : ∀ j, q[t.pc]? = some j → j.IsRandFree := by + rintro (_ | _ | _ | _ | _) hj <;> first | trivial | exact ht (step_rand hj) + +/-- A step of a library program that continues is a step of its translation. -/ +theorem step_next {w : ℕ} {q : Program .exact r} {t t' : Cslib.WordRAM.State r} + (ht : Cslib.WordRAM.step w q t = .next t') : + S.step w (S.programOf q) (stateOf S t) = S.next (stateOf S t') := by + have hnr := isRandFree_of_step (w := w) (q := q) (t := t) (by simp [ht]) + obtain ⟨pc, reg, mem⟩ := t + rcases S.step_matches_at h w q pc reg mem hnr with ⟨_, _, _, h₁, h₂⟩ | ⟨_, _, h₂⟩ | ⟨_, h₂⟩ <;> + rw [ht] at h₂ <;> cases h₂ + exact h₁ + +/-- A step of a library program that halts is a step of its translation. -/ +theorem step_halt {w : ℕ} {q : Program .exact r} {t : Cslib.WordRAM.State r} {v : ℕ} + (ht : Cslib.WordRAM.step w q t = .halt v) : + S.step w (S.programOf q) (stateOf S t) = S.halted v := by + have hnr := isRandFree_of_step (w := w) (q := q) (t := t) (by simp [ht]) + obtain ⟨pc, reg, mem⟩ := t + rcases S.step_matches_at h w q pc reg mem hnr with ⟨_, _, _, _, h₂⟩ | ⟨_, h₁, h₂⟩ | ⟨_, h₂⟩ <;> + rw [ht] at h₂ <;> cases h₂ + exact h₁ + +/-- **Runs agree.** A run of `k` steps of a library program is a run of `k` steps of its +translation, between the corresponding states. -/ +theorem runs_map {w : ℕ} {q : Program .exact r} {s t : Cslib.WordRAM.State r} {k : ℕ} + (hr : Runs w q s t k) : + Relation.RelatesInSteps (fun a b => S.step w (S.programOf q) a = S.next b) (stateOf S s) + (stateOf S t) k := by + induction hr with + | refl => exact .refl _ + | tail _ _ _ _ hn ih => exact .tail _ _ _ _ ih (step_next h hn) + +end Machine + +/-- **Transfer of the main theorem.** `Cslib.Algorithms.WordRAM.Dijkstra.logLinearStatement`, +for the instance `⟨n, G, s⟩`, with the distance given by walks, and on any presentation `S` of the +exact word RAM; `mem t` is the memory of the abstract state `t`. -/ +theorem logLinear (S : WordRAMBridge.Presentation) (h : WordRAMBridge.Laws S) + (mem : ∀ {r : ℕ}, S.State r → ℕ → ℕ) + (hmem : ∀ {r : ℕ} (pc : ℕ) (reg : Fin r → ℕ) (m : ℕ → ℕ), mem (S.state pc reg m) = m) + (hW : WalkLaws W) : + ∃ r, ∃ p : List (S.Instr r), ∃ K, ∀ (n : ℕ) (G : Fin n → List (Fin n × ℕ)) (s : Fin n) + (w : ℕ), 8 * (n + (ListInstance.mk n G s).m + 2) * + ((ListInstance.mk n G s).maxWeight + 2) < 2 ^ w → + ∃ t steps, Relation.RelatesInSteps (fun a b => S.step w p a = S.next b) + (S.init (ListInstance.mk n G s).encode) t steps ∧ + (∃ v, S.step w p t = S.halted v) ∧ + steps + 1 ≤ K * ((ListInstance.mk n G s).m + 1) * (Nat.log 2 n + 1) ∧ + ∀ v : Fin n, mem t (n + 2 * (ListInstance.mk n G s).m + 4 + v) = + encodeDist (walkDist W G s v) := by + obtain ⟨r, q, K, hq⟩ := logLinearStatement + refine ⟨r, S.programOf q, K, fun n G s w hw => ?_⟩ + obtain ⟨t, steps, hr, ⟨v, hv⟩, hle, hd⟩ := hq ⟨n, G, s⟩ w hw + refine ⟨Machine.stateOf S t, steps, ?_, ⟨v, Machine.step_halt h hv⟩, hle, fun u => ?_⟩ + · rw [← Machine.stateOf_init] + exact Machine.runs_map h hr + · rw [Machine.stateOf, hmem, walkDist_eq hW ⟨n, G, s⟩] + exact hd u + +end ComplexityAbstracts.DijkstraBridge diff --git a/references.bib b/references.bib index 945d5311ae..12977eaa41 100644 --- a/references.bib +++ b/references.bib @@ -1177,6 +1177,16 @@ @article{DietzfelbingerHagerupKatajainenPenttonen1997 doi = {10.1006/jagm.1997.0873} } +@article{Dijkstra1959, + author = {Dijkstra, Edsger W.}, + title = {A Note on Two Problems in Connexion with Graphs}, + journal = {Numerische Mathematik}, + volume = {1}, + pages = {269--271}, + year = {1959}, + doi = {10.1007/BF01386390} +} + @article{DIL2007, author = {Demaine, Erik D. and Iacono, John and Langerman, Stefan}, title = {Retroactive Data Structures},