Theory IMP2_Binary_Heap

theory IMP2_Binary_Heap
  imports IMP2.IMP2 IMP2.IMP2_Aux_Lemmas
begin
section ‹Introduction›
text ‹In this submission imperative versions of the following array-based binary minimum heap
      functions are implemented and verified: insert, get-min, delete-min, make-heap.
      The latter three are then used to prove the correctness of an in-place heapsort, which sorts
      an array in descending order. To do that in Isabelle/HOL, the proof framework IMP2
      \<^cite>‹"IMP2-AFP"› is used. Here arrays are modeled by ‹int ⇒ int› functions. The imperative
      implementations are iterative versions of the partly recursive algorithms described in
      \<^cite>‹"MS"› and \<^cite>‹"CLRS"›.

      This submission starts with the basic definitions and lemmas, which are needed for array-based 
      binary heaps. These definitions and lemmas are parameterised with an arbitrary (transitive) 
      comparison function (where such a function is needed), so they are not only applicable to 
      minimum heaps. After some more general, useful lemmas on arrays, the imperative minimum heap 
      functions and the heapsort are implemented and verified.›

section ‹Heap Related Definitions and Theorems›
subsection ‹Array Bounds›
text ‹A small helper function is used to define valid array indices. Note that the lower index 
      bound ‹l› is arbitrary and not fixed to 0 or 1. The upper index bound ‹r› is not a valid
      index itself, so that the empty array can be denoted by having ‹l = r›.›
abbreviation bounded :: "int ⇒ int ⇒ int ⇒ bool" where
  "bounded l r x ≡ l ≤ x ∧ x < r"

subsection ‹Parent and Children›

subsubsection ‹Definitions›
text ‹For the notion of an array-based binary heap, the parent and child relations on the array 
      indices need to be defined.›
definition parent :: "int ⇒ int ⇒ int" where
  "parent l c = l + (c - l - 1) div 2"

definition l_child :: "int ⇒ int ⇒ int" where
  "l_child l p = 2 * p - l + 1"

definition r_child :: "int ⇒ int ⇒ int" where
  "r_child l p = 2 * p - l + 2"

subsubsection ‹Lemmas›
lemma parent_upper_bound: "parent l c < c ⟷ l ≤ c"
  unfolding parent_def by auto

lemma parent_upper_bound_alt: "l ≤ parent l c ⟹ parent l c < c"
  unfolding parent_def by simp

lemma parent_lower_bound: "l ≤ parent l c ⟷ l < c"
  unfolding parent_def by linarith

lemma grand_parent_upper_bound: "parent l (parent l c) < c ⟷ l ≤ c"
  unfolding parent_def by linarith

corollary parent_bounds: "l < x ⟹ x < r ⟹ bounded l r (parent l x)"
  using parent_lower_bound parent_upper_bound_alt by fastforce


lemma l_child_lower_bound: " p < l_child l p ⟷ l ≤ p"
  unfolding l_child_def by simp

corollary l_child_lower_bound_alt: "l ≤ x ⟹ x ≤ p ⟹ x < l_child l p"
  using l_child_lower_bound[of p l] by linarith

lemma parent_l_child[simp]: "parent l (l_child l n) = n"
  unfolding parent_def l_child_def by simp


lemma r_child_lower_bound: "l ≤ p ⟹ p < r_child l p"
  unfolding r_child_def by simp

corollary r_child_lower_bound_alt: "l ≤ x ⟹ x ≤ p ⟹ x < r_child l p"
  using r_child_lower_bound[of l p] by linarith

lemma parent_r_child[simp]: "parent l (r_child l n) = n"
  unfolding parent_def r_child_def by simp


lemma smaller_l_child: "l_child l x < r_child l x"
  unfolding l_child_def r_child_def by simp

lemma parent_two_children: 
  "(c = l_child l p ∨ c = r_child l p) ⟷ parent l c = p"
  unfolding parent_def l_child_def r_child_def by linarith

subsection ‹Heap Invariants›
subsubsection ‹Definitions›
text ‹The following heap invariants and the following lemmas are parameterised with an arbitrary 
      (transitive) comparison function. For the concrete function implementations at the end of 
      this submission ‹≤› on ints is used.›

text ‹For the ‹make_heap› function, which transforms an unordered array into a valid heap,
      the notion of a partial heap is needed. Here the heap invariant only holds for array indices
      between a certain valid array index ‹m› and ‹r›. The standard heap invariant is then
      simply the special case where ‹m = l›.›
definition is_partial_heap 
  :: "('a::order ⇒ 'a::order ⇒ bool) ⇒ (int ⇒ 'a::order) ⇒ int ⇒ int ⇒ int ⇒ bool" where
  "is_partial_heap cmp heap l m r = (∀ x. bounded m r x ⟶ 
               bounded m r (parent l x) ⟶ cmp (heap (parent l x)) (heap x))"

abbreviation is_heap 
  :: "('a::order ⇒ 'a::order ⇒ bool) ⇒ (int ⇒ 'a::order) ⇒ int ⇒ int ⇒ bool" where
  "is_heap cmp heap l r ≡ is_partial_heap cmp heap l l r"


text ‹During all of the modifying heap functions the heap invariant is temporarily violated at
      a single index ‹i› and it is then gradually restored by either ‹sift_down› or 
      ‹sift_up›. The following definitions formalize these weakened invariants.

      The second part of the conjunction in the following definitions states, that the comparison 
      between the parent of ‹i› and each of the children of ‹i› evaluates to ‹True› without 
      explicitly using the child relations.›
definition is_partial_heap_except_down 
  :: "('a::order ⇒ 'a::order ⇒ bool) ⇒ (int ⇒ 'a::order) ⇒ int ⇒ int ⇒ int ⇒ int ⇒ bool" where
  "is_partial_heap_except_down cmp heap l m r i = (∀ x. bounded m r x ⟶
    ((parent l x ≠ i ⟶ bounded m r (parent l x) ⟶ cmp (heap (parent l x)) (heap x)) ∧ 
     (parent l x = i ⟶ bounded m r (parent l (parent l x))  
                     ⟶ cmp (heap (parent l (parent l x))) (heap x))))"

abbreviation is_heap_except_down 
  :: "('a::order ⇒ 'a::order ⇒ bool) ⇒ (int ⇒ 'a::order) ⇒ int ⇒ int ⇒ int ⇒ bool" where
  "is_heap_except_down cmp heap l r i ≡ is_partial_heap_except_down cmp heap l l r i"


text ‹As mentioned the notion of a partial heap is only needed for ‹make_heap›, 
      which only uses ‹sift_down› internally, so there doesn't need to be an additional 
      definition for the partial heap version of the ‹sift_up› invariant.›
definition is_heap_except_up 
  :: "('a::order ⇒ 'a::order ⇒ bool) ⇒ (int ⇒ 'a::order) ⇒ int ⇒ int ⇒ int ⇒ bool" where
  "is_heap_except_up cmp heap l r i = (∀ x. bounded l r x ⟶ 
    ((x ≠ i ⟶ bounded l r (parent l x) ⟶ cmp (heap (parent l x)) (heap x)) ∧ 
     (parent l x = i ⟶ bounded l r (parent l (parent l x)) 
                     ⟶ cmp (heap (parent l (parent l x))) (heap x))))"

subsubsection ‹Lemmas›

lemma empty_partial_heap[simp]: "is_partial_heap cmp heap l r r"
  unfolding is_partial_heap_def by linarith

lemma is_partial_heap_smaller_back: 
  "is_partial_heap cmp heap l m r ⟹ r' ≤ r ⟹ is_partial_heap cmp heap l m r'"
  unfolding is_partial_heap_def by simp

lemma is_partial_heap_smaller_front: 
  "is_partial_heap cmp heap l m r ⟹ m ≤ m' ⟹ is_partial_heap cmp heap l m' r"
  unfolding is_partial_heap_def by simp

text ‹The second half of each array is a is a partial binary heap, since it contains only leafs,
      which are all trivial binary heaps.›
lemma snd_half_is_partial_heap: 
  "(l + r) div 2 ≤ m ⟹ is_partial_heap cmp heap l m r" 
  unfolding is_partial_heap_def parent_def by linarith

lemma modify_outside_partial_heap: 
  assumes 
    "heap = heap' on {m..<r}" 
    "is_partial_heap cmp heap l m r"
  shows "is_partial_heap cmp heap' l m r"
  using assms eq_onD unfolding is_partial_heap_def by fastforce

text ‹The next few lemmas formalize how the heap invariant is weakened, when the heap is modified
      in a certain way.›

text ‹This lemma is used by ‹make_heap›.›
lemma partial_heap_added_first_el: 
  assumes 
    "l ≤ m" "m ≤ r"
    "is_partial_heap cmp heap l (m + 1) r"
  shows "is_partial_heap_except_down cmp heap l m r m"
 unfolding is_partial_heap_except_down_def
proof
  fix x
  let ?p_x = "parent l x"
  let ?gp_x = "parent l ?p_x"
  show "bounded m r x ⟶
        (?p_x ≠ m ⟶ bounded m r ?p_x ⟶ cmp (heap ?p_x) (heap x)) ∧
        (?p_x = m ⟶ bounded m r ?gp_x ⟶ cmp (heap ?gp_x) (heap x))"
  proof
    assume x_bound: "bounded m r x"
    have p_x_lower: "?p_x ≠ m ⟶ bounded m r ?p_x ⟶ ?p_x ≥ m + 1"
      by simp
    hence "?p_x ≠ m ⟶ bounded m r ?p_x ⟶ x ≥ m + 1"
      using parent_upper_bound[of l x] x_bound assms(1) by linarith
    hence p_invariant: "(?p_x ≠ m ⟶ bounded m r ?p_x ⟶ cmp (heap ?p_x) (heap x))"
      using assms(3) is_partial_heap_def p_x_lower x_bound by blast

    have gp_up_bound: "(?p_x = m ⟶ ?gp_x < m)"
      by (simp add: assms(1) parent_upper_bound)
    show "(?p_x ≠ m ⟶ bounded m r ?p_x ⟶ cmp (heap ?p_x) (heap x)) ∧
          (?p_x = m ⟶ bounded m r ?gp_x ⟶ cmp (heap ?gp_x) (heap x))"
      using gp_up_bound p_invariant by linarith
  qed
qed

text ‹This lemma is used by ‹del_min›.›
lemma heap_changed_first_el: 
  assumes "is_heap cmp heap l r" "l ≤ r"
  shows "is_heap_except_down cmp (heap(l := b)) l r l"
proof -
  have "is_partial_heap cmp heap l (l + 1) r"
    using assms(1) is_partial_heap_smaller_front by fastforce
  hence "is_partial_heap cmp (heap(l := b)) l (l + 1) r"
    using modify_outside_partial_heap[of heap] by simp
  thus ?thesis
    by (simp add: assms(2) partial_heap_added_first_el)
qed

text ‹This lemma is used by ‹insert›.›
lemma heap_appended_el: 
  assumes 
    "is_heap cmp heap l r" 
    "heap = heap' on {l..<r}"
  shows "is_heap_except_up cmp heap' l (r+1) r"
proof -
  have "is_heap cmp heap' l r" 
    using assms(1,2) modify_outside_partial_heap by blast 
  thus ?thesis unfolding is_partial_heap_def is_heap_except_up_def
    by (metis not_less_iff_gr_or_eq parent_upper_bound zless_add1_eq)
qed

subsubsection ‹First Heap Element›
text ‹The next step is to show that the first element of the heap is always the ``smallest'' 
      according to the given comparison function. For the proof a rule for strong induction on lower
      bounded integers is needed. Its proof is based on the proof of strong induction on natural 
      numbers found in \<^cite>‹"Str_Ind"›.›
lemma strong_int_gr_induct_helper: 
  assumes "k < (i::int)" "(⋀i. k < i ⟹ (⋀j. k < j ⟹ j < i ⟹ P j) ⟹ P i)"
  shows "⋀j. k < j ⟹ j < i ⟹ P j"
  using assms
proof(induction i rule: int_gr_induct)
  case base 
  then show ?case by linarith 
next
  case (step i) 
  then show ?case 
  proof(cases "j = i")
    case True
    then show ?thesis using step.IH step.prems(1,3) by blast 
  next
    case False
    hence "j < i" using step.prems(2) by simp
    then show ?thesis using step.IH step.prems(1,3) by blast 
  qed
qed

theorem strong_int_gr_induct:
  assumes 
    "k < (i::int)" 
    "(⋀i. k < i ⟹ (⋀j. k < j ⟹ j < i ⟹ P j) ⟹ P i)" 
  shows "P i"
  using assms less_induct strong_int_gr_induct_helper by blast

text ‹Now the main theorem, that the first heap element is the ``smallest'' according to the
      given comparison function, can be proven.›
theorem heap_first_el:
  assumes 
    "is_heap cmp heap l r"
    "transp cmp"
    "l < x" "x < r"
  shows "cmp (heap l) (heap x)"
  using assms unfolding is_partial_heap_def
proof(induction x rule: strong_int_gr_induct[of l])
  case 1
  then show ?case using assms(3) by simp
next
  case (2 i)
  have cmp_pi_i: "cmp (heap (parent l i)) (heap i)"
    using "2.hyps" "2.prems"(1,4) parent_bounds by simp
  then show ?case 
  proof(cases)
    assume "parent l i > l"
    then have "cmp (heap l) (heap (parent l i))"
      using "2.IH" "2.prems"(1,2,4) parent_upper_bound_alt by simp
    then show ?thesis
      using "2.prems"(2) cmp_pi_i transpE by metis
  next
    assume "¬ parent l i > l"
    then have "parent l i = l"
      using "2.hyps" dual_order.order_iff_strict parent_lower_bound by metis
    then show ?thesis
      using cmp_pi_i by simp
  qed
qed



section ‹General Lemmas on Arrays›
text ‹Some additional lemmas on @{const "mset_ran"}, @{const "swap"} and @{const "eq_on"} are needed
      for the final proofs.›

subsection ‹Lemmas on @{const "mset_ran"}›
abbreviation arr_mset :: "(int ⇒ 'a) ⇒ int ⇒ int ⇒ 'a multiset" where
  "arr_mset arr l r ≡ mset_ran arr {l..<r}"

lemma in_mset_imp_in_array: 
  "x ∈# (arr_mset arr l r) ⟷ (∃i. bounded l r i ∧ arr i = x)"
  unfolding mset_ran_def by fastforce

lemma arr_mset_remove_last: 
  "l ≤ r ⟹ arr_mset arr l r = arr_mset arr l (r + 1) - {#arr r#}"
  by (simp add: intvs_upper_decr mset_ran_def)

lemma arr_mset_append: 
  "l ≤ r ⟹ arr_mset arr l (r + 1) = arr_mset arr l r + {#arr r#}"
  using arr_mset_remove_last[of l r arr] by simp

corollary arr_mset_append_alt: 
  "l ≤ r ⟹ arr_mset (arr(r := b)) l (r + 1) = arr_mset arr l r + {#b#}"
  by (simp add: arr_mset_append mset_ran_subst_outside)

lemma arr_mset_remove_first: 
  "i ≤ r ⟹ arr_mset arr (i - 1) r = arr_mset arr i r + {#arr (i - 1)#}"
  by(induction r rule: int_ge_induct) (auto simp add: arr_mset_append)

lemma arr_mset_split: 
  assumes "l ≤ m" "m ≤ r"
  shows "arr_mset arr l r = arr_mset arr l m + arr_mset arr m r"
  using assms
proof(induction m rule: int_ge_induct[of l])
  case (step i)
  have add_last: "arr_mset arr l (i + 1) = arr_mset arr l i + {#arr i#}"
    using step arr_mset_append by blast
  have rem_first: "arr_mset arr (i+1) r = arr_mset arr i r - {#arr i#}"
    by (metis step.prems arr_mset_remove_first add_diff_cancel_right') 
  show ?case 
    using step add_last rem_first by fastforce
qed (simp)

text ‹That the first element in a heap is the ``smallest'', can now be expressed using multisets.›
corollary heap_first_el_alt:                                                   
  assumes 
    "transp cmp" 
    "is_heap cmp heap l r" 
    "x ∈# (arr_mset heap l r)" 
    "heap l ≠ x"
  shows "cmp (heap l) x"
  by(metis assms heap_first_el in_mset_imp_in_array le_less)


subsection ‹Lemmas on @{term "swap"} and @{term "eq_on"}›

lemma eq_on_subset: 
  "arr1 = arr2 on R ⟹ S ⊆ R ⟹ arr1 = arr2 on S"
  by (simp add: eq_on_def set_mp)

lemma swap_swaps: 
  "arr' = swap arr x y ⟹ arr' y = arr x ∧ arr' x = arr y"
  unfolding swap_def by simp

lemma swap_only_swaps: 
  "arr' = swap arr x y ⟹ z ≠ x ⟹ z ≠ y ⟹ arr' z = arr z"
  unfolding swap_def by simp

lemma swap_commute: "swap arr x y = swap arr y x"
  unfolding swap_def by fastforce

lemma swap_eq_on: 
  "arr1 = arr2 on S ⟹ x ∉ S ⟹ y ∉ S ⟹ arr1 = swap arr2 x y on S"
  unfolding swap_def by simp

corollary swap_parent_eq_on: 
  assumes 
    "arr1 = arr2 on - {l..<r}" 
    "l < c" "c < r" 
  shows "arr1 = swap arr2 (parent l c) c on - {l..<r} "
  using parent_bounds swap_eq_on assms by fastforce

corollary swap_child_eq_on: 
  assumes 
    "arr1 = arr2 on - {l..<r}" 
    "c = l_child l p ∨ c = r_child l p" 
    "l ≤ p" "c < r" 
  shows "arr1 = swap arr2 p c on - {l..<r} "
  by (metis assms parent_lower_bound parent_two_children swap_parent_eq_on)

lemma swap_child_mset: 
  assumes 
    "arr_mset arr1 l r = arr_mset arr2 l r" 
    "c = l_child l p ∨ c = r_child l p" 
    "l ≤ p" "c < r" 
  shows "arr_mset arr1 l r = arr_mset (swap arr2 p c) l r"
proof -
  have child_bounded: "l < c ∧ c < r"
    by (metis assms(2-4) parent_lower_bound parent_two_children)
  have parent_bounded: "bounded l r p"
    by (metis assms(2-4) dual_order.strict_trans parent_two_children parent_upper_bound_alt)
  thus ?thesis
    using assms(1) child_bounded mset_ran_swap[of p "{l..<r}" c arr2] atLeastLessThan_iff
    by simp
qed

text ‹The following lemma shows, which propositions have to hold on the pre-swap array, so that
      a comparison between two elements holds on the post-swap array. This is useful for the 
      proofs of the loop invariants of ‹sift_up› and ‹sift_down›. The lemma is kept 
      quite general (except for the argument order) and could probably be more closely
      related to the parent relation for more concise proofs.›
lemma cmp_swapI: 
  fixes arr::"'a::order ⇒ 'a::order"
  assumes
    "m < n ∧ x < y"
    "m < n ∧ x < y ⟹ x = m ⟹ y = n ⟹ P (arr n) (arr m)"
    "m < n ∧ x < y ⟹ x ≠ m ⟹ x ≠ n ⟹ y ≠ m ⟹ y ≠ n ⟹ P (arr m) (arr n)"
    "m < n ∧ x < y ⟹ x = m ⟹ y ≠ m ⟹ y ≠ n ⟹  P (arr y) (arr n)"
    "m < n ∧ x < y ⟹ x = n ⟹ y ≠ m ⟹ y ≠ n ⟹  P (arr m) (arr y)"
    "m < n ∧ x < y ⟹ x ≠ m ⟹ x ≠ n ⟹ y = n ⟹ P (arr m) (arr x)"
    "m < n ∧ x < y ⟹ x ≠ m ⟹ x ≠ n ⟹ y = m ⟹ P (arr x) (arr n)"
  shows "P (swap arr x y m) (swap arr x y n)"
  by (metis assms order.asym swap_only_swaps swap_swaps)

section ‹Imperative Heap Implementation›
text ‹The following imperative heap functions are based on \<^cite>‹"MS"› and \<^cite>‹"CLRS"›. All functions,
      that are recursive in these books, are iterative in the following implementations. The 
      function definitions are done with IMP2 \<^cite>‹"IMP2-AFP"›. From now on the heaps only contain 
      ints and only use ‹≤› as comparison function. The auxiliary lemmas used from now on are 
      heavily modeled after the proof goals, that are generated by the vcg tool (also part of IMP2).›

subsection ‹Simple Functions›
subsubsection ‹Parent, Children and Swap›
text ‹In this section the parent and children relations are expressed as IMP2 procedures. 
      Additionally a simple procedure, that swaps two array elements, is defined.›
procedure_spec prnt (l, x) returns p
  assumes True
    ensures "p = parent l⇩0 x⇩0"
  defines ‹p = ((x - l - 1) / 2 + l)›
  by vcg (simp add: parent_def)

procedure_spec left_child (l, x) returns lc
  assumes True 
    ensures "lc = l_child l⇩0 x⇩0"
  defines ‹lc = 2 * x - l + 1›
  by vcg (simp add: l_child_def)

procedure_spec right_child (l, x) returns rc
  assumes True
    ensures "rc = r_child l⇩0 x⇩0"
  defines ‹rc = 2 * x - l + 2›
  by vcg (simp add: r_child_def)

procedure_spec swp (heap, x, y) returns heap
  assumes True
    ensures "heap = swap heap⇩0 x⇩0 y⇩0 "
  defines ‹tmp = heap[x]; heap[x] = heap[y]; heap[y] = tmp›
  by vcg (simp add: swap_def)

subsubsection ‹‹get_min››
text ‹In this section ‹get_min› is defined, which simply returns the first element (the minimum) of 
      the heap. For this definition an additional theorem is proven, which enables the use of
      @{const "Min_mset"} in the postcondition.›

theorem heap_minimum: 
  assumes 
    "l < r" 
    "is_heap (≤) heap l r"
  shows "heap l = Min_mset (arr_mset heap l r)"
proof -
  have "(∀x ∈# (arr_mset heap l r). (heap l) ≤ x)"
    using assms(2) heap_first_el_alt transp_on_le by blast
  thus ?thesis
    by (simp add: assms(1) dual_order.antisym)
qed

procedure_spec get_min (heap, l, r) returns min
  assumes "l < r ∧ is_heap (≤) heap l r"  
    ensures "min = Min_mset (arr_mset heap⇩0 l⇩0 r⇩0)"
  for heap[] l r
  defines ‹min = heap[l]›
  by vcg (simp add: heap_minimum)

subsection ‹Modifying Functions›

subsubsection ‹‹sift_up› and ‹insert››
text ‹The next heap function is ‹insert›, which internally uses ‹sift_up›. In the beginning of 
      this section ‹sift_up_step› is proven, which states that each ‹sift_up› loop iteration 
      correctly transforms the weakened heap invariant. For its proof two additional
      auxiliary lemmas are used. After ‹sift_up_step› ‹sift_up› and then ‹insert› are verified.›

text ‹‹sift_up_step› can be proven directly by the smt-solver without auxiliary lemmas, but they
      were introduced to show the proof details. The analogous proofs for ‹sift_down› were 
      just solved with smt, since the proof structure should be very similar, even though the 
      ‹sift_down› proof goals are slightly more complex.›
lemma sift_up_step_aux1:
  fixes heap::"int ⇒ int"
  assumes  
    "is_heap_except_up (≤) heap l r x"
    "parent l x ≥ l" 
    "(heap x) ≤ (heap (parent l x))" 
    "bounded l r k" 
    "k ≠ (parent l x)" 
    "bounded l r (parent l k)"
  shows "(swap heap (parent l x) x (parent l k)) ≤ (swap heap (parent l x) x k)"
  apply(intro cmp_swapI[of "(parent l k)" k "(parent l x)" x "(≤)" heap])
  subgoal using assms(2,6) parent_upper_bound_alt by blast
  subgoal using assms(3) by blast
  subgoal using assms(1,4,6) unfolding is_heap_except_up_def by blast
  subgoal using assms(1,3,4,6) unfolding is_heap_except_up_def by fastforce
  subgoal using assms(5) by blast
  subgoal by blast 
  subgoal using assms(1,2,4) unfolding is_heap_except_up_def by simp
  done

lemma sift_up_step_aux2: 
  fixes heap::"int ⇒ int"
  assumes
    "is_heap_except_up (≤) heap l r x"
    "parent l x ≥ l"
    "heap x ≤ (heap (parent l x))"
    "bounded l r k"
    "parent l k = parent l x"
    "bounded l r (parent l (parent l k))"
  shows 
    "swap heap (parent l x) x (parent l (parent l k)) ≤ swap heap (parent l x) x k"
  using assms unfolding is_heap_except_up_def
proof-
  let ?gp_k = "parent l (parent l k)"
  let ?gp_x = "parent l (parent l x)"
  have gp_k_eq_gp_x: "swap heap (parent l x) x ?gp_k = heap ?gp_x"
    by (metis assms(2,5) grand_parent_upper_bound less_irrefl swap_only_swaps)
  show ?thesis
      using assms unfolding is_heap_except_up_def
  proof(cases)
    assume k_eq_x: "k = x"
    have "swap heap (parent l x) x k = heap (parent l x)"
      by (metis k_eq_x swap_swaps)
    then show ?thesis
      using assms(1,2,4,6) unfolding is_heap_except_up_def
      by (metis gp_k_eq_gp_x k_eq_x parent_bounds parent_lower_bound)
  next
    assume k_neq_x: "k ≠ x"
    have "swap heap (parent l x) x k = heap k"
      by (metis assms(5) gp_k_eq_gp_x k_neq_x swap_only_swaps)
    then show ?thesis using assms unfolding is_heap_except_up_def
      by (metis gp_k_eq_gp_x k_neq_x order_trans parent_bounds parent_lower_bound)
  qed
qed

lemma sift_up_step:
  fixes heap::"int ⇒ int"
  assumes  
    "is_heap_except_up (≤) heap l r x"
    "parent l x ≥ l" 
    "(heap x) ≤ (heap (parent l x))"
  shows "is_heap_except_up (≤) (swap heap (parent l x) x) l r (parent l x)"
  using assms sift_up_step_aux1 sift_up_step_aux2
  unfolding is_heap_except_up_def by blast

text ‹‹sift_up› restores the heap invariant, that is only violated at the current position, by 
      iteratively swapping the current element with its parent until the beginning of the array is 
      reached or the current element is bigger than its parent.›
procedure_spec sift_up (heap, l, r, x) returns heap
  assumes "is_heap_except_up (≤) heap l r x ∧ bounded l r x"
    ensures "is_heap (≤) heap l⇩0 r⇩0 ∧
             arr_mset heap⇩0 l⇩0 r⇩0 = arr_mset heap l⇩0 r⇩0 ∧
             heap⇩0 = heap on - {l⇩0..<r⇩0}"
  for heap[] l x r
  defines ‹
    p = prnt(l, x);
    while (x > l ∧ heap[x] ≤ heap[p]) 
      @variant ‹x - l›
      @invariant ‹is_heap_except_up (≤) heap l r x ∧ p = parent l x ∧ 
                  bounded l r x ∧ arr_mset heap⇩0 l⇩0 r⇩0 = arr_mset heap l r ∧
                  heap⇩0 = heap  on - {l..<r}› 
    {
        heap = swp(heap, p, x);
        x = p;
        p = prnt(l, x)
    }›
  apply vcg_cs
  apply(intro conjI)
  subgoal using parent_lower_bound sift_up_step by blast
  subgoal using parent_lower_bound by blast
  subgoal using parent_bounds by blast
  subgoal using parent_bounds by (simp add: mset_ran_swap)
  subgoal using swap_parent_eq_on by blast
  subgoal using parent_upper_bound by simp
  subgoal unfolding is_heap_except_up_def is_partial_heap_def
    by (metis le_less not_less parent_lower_bound)
  done

text ‹‹insert› inserts an element into a heap by appending it to the heap and restoring the heap 
      invariant with @{const "sift_up"}.›
procedure_spec insert (heap, l, r, el) returns (heap, l, r)
  assumes "is_heap (≤) heap l r ∧ l ≤ r"  
    ensures "is_heap (≤) heap l r ∧
             arr_mset heap l r = arr_mset heap⇩0 l⇩0 r⇩0 + {#el⇩0#} ∧
             l = l⇩0 ∧ r = r⇩0 + 1  ∧ heap⇩0 = heap on - {l..<r}"
  for heap l r el
  defines ‹
    heap[r] = el;
    x = r;
    r = r + 1;
    heap = sift_up(heap, l, r, x)
  ›
  apply vcg_cs 
  subgoal by (simp add: heap_appended_el)
  subgoal by (metis arr_mset_append_alt add_mset_add_single)
  done

subsubsection ‹‹sift_down›, ‹del_min› and ‹make_heap››
text ‹The next heap functions are ‹del_min› and ‹make_heap›, which both use ‹sift_down› to 
      restore/establish the heap invariant. ‹sift_down› is proven first (this time without 
      additional auxiliary lemmas) followed by ‹del_min› and ‹make_heap›.›

text ‹‹sift_down› restores the heap invariant, that is only violated at the current position, by 
      iteratively swapping the current element with its smallest child until the end of 
      the array is reached or the current element is smaller than its children.›
procedure_spec sift_down(heap, l, r, x) returns heap 
  assumes "is_partial_heap_except_down (≤) heap l x r x ∧ l ≤ x ∧ x ≤ r"
    ensures "is_partial_heap (≤) heap l⇩0 x⇩0 r⇩0 ∧ 
             arr_mset heap⇩0 l⇩0 r⇩0 = arr_mset heap l⇩0 r⇩0 ∧ 
             heap⇩0 = heap on - {l⇩0..<r⇩0}"
  defines ‹
   lc = left_child(l, x);
   rc = right_child(l, x);
    while (lc < r ∧ (heap[lc] < heap[x] ∨ (rc < r ∧ heap[rc] < heap[x]))) 
      @variant ‹r - x›
      @invariant ‹is_partial_heap_except_down (≤) heap l x⇩0 r x ∧
                  x⇩0 ≤ x ∧ x ≤ r ∧ lc = l_child l x ∧ rc = r_child l x ∧
                  arr_mset heap⇩0 l r = arr_mset heap l r ∧
                  heap⇩0 = heap on - {l..<r}›
  { 
    smallest = lc;
    if (rc < r ∧ heap[rc] < heap[lc]) {
      smallest = rc
    };
    heap = swp(heap, x, smallest);
    x = smallest;
    lc = left_child(l, x);
    rc = right_child(l, x)
  }›
  apply vcg_cs
  subgoal
    apply(intro conjI)
    subgoal unfolding is_partial_heap_except_down_def
      by (smt parent_two_children swap_swaps swap_only_swaps
          swap_commute parent_upper_bound_alt)
    subgoal using r_child_lower_bound_alt by fastforce
    subgoal using swap_child_mset order_trans by blast
    subgoal using swap_child_eq_on by fastforce
    done
  subgoal
    by (meson less_le_trans not_le order.asym r_child_lower_bound) 
  subgoal
    apply(intro conjI)
    subgoal unfolding is_partial_heap_except_down_def
      by (smt parent_two_children swap_swaps swap_only_swaps
          swap_commute parent_upper_bound_alt)
    subgoal using l_child_lower_bound_alt by fastforce
    subgoal using swap_child_mset order_trans by blast
    subgoal using swap_child_eq_on by fastforce
    done
  subgoal 
    by (meson less_le_trans not_le order.asym l_child_lower_bound)
  subgoal unfolding is_partial_heap_except_down_def is_partial_heap_def
    by (metis dual_order.strict_trans not_less parent_two_children smaller_l_child)
  done


text ‹‹del_min› needs an additional lemma which shows, that it actually removes (only) the minimum
      from the heap.›
lemma del_min_mset: 
  fixes heap::"int ⇒ int"
  assumes 
    "l < r"
    "is_heap (≤) heap l r"
    "mod_heap = heap(l := heap (r - 1))"
    "arr_mset mod_heap l (r - 1) = arr_mset new_heap l (r - 1)"
  shows 
    "arr_mset new_heap l (r - 1) = arr_mset heap l r - {#Min_mset (arr_mset heap l r)#}"
proof -
  let ?heap_mset = "arr_mset heap l r"
  have l_is_min: "heap l = Min_mset ?heap_mset"
    using assms(1,2) heap_minimum by blast
  have "(arr_mset mod_heap l r) = ?heap_mset + {#heap (r-1)#} - {#heap l#}"
    by (simp add: assms(1,3) mset_ran_subst_inside)
  hence "(arr_mset mod_heap l (r - 1)) = ?heap_mset - {#heap l#}"
    by (simp add: assms(1,3) arr_mset_remove_last) 
  thus ?thesis
    using assms(4) l_is_min by simp
qed

text ‹‹del_min› removes the minimum element from the heap by replacing the first element with the 
      last element, shrinking the array by one and subsequently restoring the heap invariant 
      with @{const "sift_down"}.›
procedure_spec del_min (heap, l, r) returns (heap, l, r)
  assumes "l < r ∧ is_heap (≤) heap l r"  
    ensures "is_heap (≤) heap l r ∧ 
             arr_mset heap l r = arr_mset heap⇩0 l⇩0 r⇩0 - {#Min_mset (arr_mset heap⇩0 l⇩0 r⇩0)#} ∧ 
             l = l⇩0 ∧ r = r⇩0 - 1 ∧ 
             heap⇩0 = heap on - {l⇩0..<r⇩0}"
  for heap l r
  defines ‹
    r = r - 1;
    heap[l] = heap[r];
    heap = sift_down(heap, l, r, l)
  ›
  apply vcg_cs
  subgoal by (simp add: heap_changed_first_el is_partial_heap_smaller_back)
  subgoal
    apply(rule conjI)
    subgoal using del_min_mset by blast
    subgoal by (simp add: eq_on_def intvs_incdec(3) intvs_lower_incr)
    done
  done


text ‹‹make_heap› transforms an arbitrary array into a heap by iterating through all array 
      positions from the middle of the array up to the beginning of the array and calling 
      @{const "sift_down"} for each one.›
procedure_spec make_heap (heap, l, r) returns heap
  assumes "l ≤ r"  
    ensures "is_heap (≤) heap l⇩0 r⇩0 ∧ 
             arr_mset heap l⇩0 r⇩0 = arr_mset heap⇩0 l⇩0 r⇩0 ∧ 
             heap⇩0 = heap on - {l⇩0..< r⇩0}"
  for heap[] l r
  defines ‹
    y = (r + l)/2 - 1;
    while (y ≥ l)
          @variant ‹y - l + 1›
          @invariant ‹is_partial_heap (≤) heap l (y + 1) r ∧
                      arr_mset heap l r = arr_mset heap⇩0 l⇩0 r⇩0 ∧
                      l - 1 ≤ y ∧ y < r ∧ heap⇩0 = heap on - {l..<r}›
    {
      heap = sift_down(heap, l, r, y);
      y = y - 1
    }› 
  apply(vcg_cs)
  subgoal 
    apply(rule conjI)
    subgoal by (simp add: snd_half_is_partial_heap add.commute)
    subgoal by linarith 
    done
  subgoal using partial_heap_added_first_el le_less by blast 
  subgoal using eq_on_trans by blast 
  subgoal using dual_order.antisym by fastforce
  done

subsection ‹Heapsort Implementation›

text ‹The final part of this submission is the implementation of the in-place heapsort. Firstly it 
      builds the ‹≤›-heap and then it iteratively removes the minimum of the heap, which is put at 
      the now vacant end of the shrinking heap. This is done until the heap is empty, which leaves 
      the array sorted in descending order.›
subsubsection ‹Auxiliary Lemmas›
text ‹Firstly the notion of a sorted array is needed. This is more or less the same as 
      @{const "ran_sorted"} generalized for arbitrary comparison functions.›
definition array_is_sorted :: "(int ⇒ int ⇒ bool) ⇒  (int ⇒ int) ⇒ int ⇒ int ⇒ bool" where
  "array_is_sorted cmp a l r ≡ ∀i. ∀j. bounded l r i ⟶ bounded l r j ⟶ i < j ⟶ cmp (a i) (a j)"

text ‹This lemma states, that the heapsort doesn't change the elements contained in the array during
      the loop iterations.›
lemma heap_sort_mset_step:
  fixes arr::"int ⇒ int"
  assumes 
    "l < m" "m ≤ r"
    "arr_mset arr' l (m - 1) = arr_mset arr l m - {#Min_mset (arr_mset arr l m)#}"
    "arr = arr' on - {l..<m}"
    "mod_arr = arr'(m - 1 := Min_mset (arr_mset arr l m))"
  shows "arr_mset arr l r = arr_mset mod_arr l r"
proof -
  let ?min = "{#Min_mset (arr_mset arr l m)#}"
  let ?new_arr_mset = "arr_mset mod_arr"
  have middle: "?new_arr_mset (m - 1) m = ?min"
    by (simp add: assms(5))
  have first_half: "?new_arr_mset l (m - 1) = arr_mset arr l m - ?min"
    by (simp add: assms(3,5) mset_ran_subst_outside)
  hence "?new_arr_mset l m = ?new_arr_mset l (m - 1) + ?new_arr_mset (m - 1) m"
    by (metis assms(1,3,5) diff_add_cancel middle arr_mset_append_alt zle_diff1_eq)
  hence first_half_middle: "?new_arr_mset l m = arr_mset arr l m"
    using middle first_half assms(1) by simp

  hence "mod_arr = arr on - {l..<m}"
    using assms(1,4,5) eq_on_sym eq_on_trans by auto
  then have second_half: "arr_mset arr m r = arr_mset mod_arr m r"
    by (simp add: eq_on_def mset_ran_cong)

  then show ?thesis
    by (metis arr_mset_split assms(1,2) first_half_middle le_less)
qed

text ‹This lemma states, that each loop iteration leaves the growing second half of the array 
      sorted in descending order.›
lemma heap_sort_second_half_sorted_step: 
  fixes arr::"int ⇒ int"
  assumes
    "l⇩0 < m" "m ≤ r⇩0"
    "arr = arr' on - {l⇩0..<m}"
    "∀i. ∀j. bounded m r⇩0 i ⟶ bounded m r⇩0 j ⟶  i < j ⟶ arr j ≤ arr i"
    "∀x∈#arr_mset arr l⇩0 m. ∀y∈#arr_mset arr m r⇩0. ¬ x < y"
    "bounded (m - 1) r⇩0 i" 
    "bounded (m - 1) r⇩0 j" 
    "i < j"
    "mod_arr = (arr'(m - 1 := Min_mset (arr_mset arr l⇩0 m)))"
  shows "mod_arr j ≤ mod_arr i"
proof -
  have second_half_eq: "mod_arr = arr on {m..< r⇩0}" 
    using assms(3, 9) unfolding eq_on_def by simp
  have j_stricter_bound: "bounded m r⇩0 j"
    using assms(6-8) by simp
  then have el_at_j: "mod_arr j ∈# arr_mset arr m r⇩0"
    using eq_onD second_half_eq by fastforce
  then show ?thesis 
  proof(cases)
    assume "i = (m-1)"
    then have "mod_arr i ∈# arr_mset arr l⇩0 m"
      by (simp add: assms(1, 9))
    then show ?thesis
      using assms(5) el_at_j not_less by blast 
  next
    assume "i ≠ (m-1)"  
    then have "bounded m r⇩0 i"
      using assms(6) by simp
    then show ?thesis 
      using assms(4, 8) eq_on_def j_stricter_bound second_half_eq by force
  qed
qed


text ‹The following lemma shows that all elements in the first part of the array (the binary heap) 
      are bigger than the elements in the second part (the sorted part) after every iteration. This 
      lemma and the invariant of the ‹heap_sort› loop use ‹¬ x < y› instead of ‹x ≥ y› since 
      ‹vcg_cs› doesn't terminate in the latter case.›
lemma heap_sort_fst_part_bigger_snd_part_step:
  fixes arr::"int ⇒ int" 
  assumes
    "l⇩0 < m"
    "m ≤ r⇩0"
    "arr_mset arr' l⇩0 (m - 1) = arr_mset arr l⇩0 m - {#Min_mset (arr_mset arr l⇩0 m)#}"
    "arr = arr' on - {l⇩0..<m}"
    "∀x∈#arr_mset arr l⇩0 m. ∀y∈#arr_mset arr m r⇩0. ¬ x < y"
    "mod_arr = arr'(m - 1 := Min_mset (arr_mset arr l⇩0 m))"
    "x∈#arr_mset mod_arr l⇩0 (m - 1)"
    "y∈#arr_mset mod_arr (m - 1) r⇩0"
  shows "¬ x < y"
proof -
  have "{m..<r⇩0} ⊆ - {l⇩0..<m}" 
    by auto
  hence "arr' = arr on {m..<r⇩0}" 
    using assms(4) eq_on_sym eq_on_subset by blast
  hence arr_eq_on: "mod_arr = arr on {m..<r⇩0}"
    by (simp add: assms(6))
  hence same_mset: "arr_mset mod_arr m r⇩0 = arr_mset arr m r⇩0"
    using mset_ran_cong by blast
  have "x ∈# arr_mset arr l⇩0 m" using same_mset assms
    by (metis assms(3,6,7) add_mset_remove_trivial_eq lran_upd_outside(2)
        mset_lran cancel_ab_semigroup_add_class.diff_right_commute
        diff_single_trivial multi_self_add_other_not_self order_refl)
  then have x_bigger_min: "x ≥ Min_mset (arr_mset arr l⇩0 m)"
    using Min_le by blast 
  have y_smaller_min: "y ≤ Min_mset (arr_mset arr l⇩0 m)" 
  proof(cases "y = mod_arr (m - 1)")
    case False
    hence "y∈#arr_mset mod_arr (m - 1) r⇩0 - {#mod_arr (m - 1)#}"
      by (metis assms(8) diff_single_trivial insert_DiffM insert_noteq_member)
    then have "y∈#arr_mset arr m r⇩0"
      by (simp add: assms(2) intvs_decr_l mset_ran_insert same_mset)
    then show ?thesis
      using assms(1) assms(5) by fastforce 
  qed (simp add: assms(6))
  then show ?thesis
    using x_bigger_min by linarith 
qed

subsubsection ‹Implementation›
text ‹Now finally the correctness of the ‹heap_sort› is shown. As mentioned, it starts by 
      transforming the array into a minimum heap using @{const "make_heap"}. Then in each 
      iteration it removes the first element from the heap with @{const "del_min"} after its value 
      was retrieved with @{const "get_min"}. This value is then put at the position freed by 
      @{const "del_min"}.›
program_spec heap_sort
  assumes "l ≤ r"  
    ensures "array_is_sorted (≥) arr l⇩0 r⇩0 ∧ 
             arr_mset arr⇩0 l⇩0 r⇩0 = arr_mset arr l⇩0 r⇩0 ∧ 
             arr⇩0 = arr on - {l⇩0 ..<r⇩0 } ∧ l = l⇩0 ∧ r = r⇩0"
  for l r arr[]
  defines ‹
    arr = make_heap(arr, l, r);
    m = r;
    while (m > l)
      @variant ‹m - l + 1› 
      @invariant ‹is_heap (≤) arr l m ∧
        array_is_sorted (≥) arr m r⇩0 ∧
        (∀x ∈# arr_mset arr l⇩0 m. ∀y ∈# arr_mset arr m r⇩0. ¬ x < y) ∧
        arr_mset arr⇩0 l⇩0 r⇩0 = arr_mset arr l⇩0 r⇩0 ∧ 
        l ≤ m ∧ m ≤ r⇩0 ∧ l = l⇩0 ∧ arr⇩0 = arr on - {l⇩0 ..<r⇩0}›
    {
      min = get_min(arr, l, m);
      (arr, l, m) = del_min(arr, l, m);
      arr[m] = min 
    }
  ›
  apply vcg_cs
  subgoal unfolding array_is_sorted_def by simp
  subgoal
    apply(intro conjI)
    subgoal unfolding is_partial_heap_def by simp
    subgoal unfolding array_is_sorted_def using heap_sort_second_half_sorted_step
      by blast
    subgoal using heap_sort_fst_part_bigger_snd_part_step by blast
    subgoal using heap_sort_mset_step by blast
    subgoal unfolding eq_on_def
      by (metis ComplD ComplI atLeastLessThan_iff less_le_trans)
    done
  done
end