Quelle  Examples.thy

theory Examples
imports "../IMP2" "../lib/IMP2_Aux_Lemmas"
begin

section ‹Examples›  
lemmas nat_distribs = nat_add_distrib nat_diff_distrib Suc_diff_le nat_mult_distrib nat_div_distrib

subsection ‹Common Loop Patterns›

subsubsection ‹Count Up›
text ‹
 Counter ‹c› counts from ‹0› to ‹n›, such that loop is executed ‹n› times.
 The result is computed in an accumulator ‹a›.
 ›
text ‹The invariant states that we have computed the function for the counter value ‹c››
text ‹The variant is the difference between ‹n› and ‹c›, i.e., the number of
 loop iterations that we still have to do›  


program_spec exp_count_up
  assumes "0≤n"
  ensures "a = 2^nat n₀"
  defines ‹
 a = 1;
 c = 0;
 while (c<n)
 @variant ‹n-c›
 @invariant ‹0≤c ∧ c≤n ∧ a=2^nat c›
 {
 G_par = a; scope { a = G_par; a=2*a; G_ret = a }; a = G_ret;
 c=c+1
 }
 ›
  apply vcg
  by (auto simp: algebra_simps nat_distribs)

program_spec sum_prog
  assumes "n ≥ 0" ensures "s = ∑{0..n₀}"
  defines ‹
 s = 0;
 i = 0;
 while (i < n)
 @variant ‹n₀ - i›
 @invariant ‹n₀ = n ∧ 0 ≤ i ∧ i ≤ n ∧ s = ∑{0..i}›
 {
 i = i + 1;
 s = s + i
 }
 ›
  apply vcg_cs
  by (simp_all add: intvs_incdec)
  
program_spec sq_prog
  assumes "n ≥ 0" ensures "a = n₀ * n₀"
  defines ‹
 a = 0;
 z = 1;
 i = 0;
 while (i < n)
 @variant ‹n₀ - i›
 @invariant ‹n₀ = n ∧ 0 ≤ i ∧ i ≤ n ∧ a = i * i ∧ z = 2 * i + 1›
 {
 a = a + z;
 z = z + 2;
 i = i + 1
 }
 ›
  by vcg_cs (simp add: algebra_simps)
  
fun factorial :: "int ==> int" where
  "factorial i = (if i ≤ 0 then 1 else i * factorial (i - 1))"
  
program_spec factorial_prog
  assumes "n ≥ 0" ensures "a = factorial n₀"
  defines ‹
 a = 1;
 i = 1;
 while (i ≤ n)
 @variant ‹n₀ + 1 - i›
 @invariant ‹n₀ = n ∧ 1 ≤ i ∧ i ≤ n + 1 ∧ a = factorial (i - 1)›
 {
 a = a * i;
 i = i + 1
 }
 ›
  by vcg (simp add: antisym_conv)+

  
fun fib :: "int ==> int" where
  "fib i = (if i ≤ 0 then 0 else if i = 1 then 1 else fib (i - 2) + fib (i - 1))"

lemma fib_simps[simp]:
  "i ≤ 0 ==> fib i = 0"
  "i = 1 ==> fib i = 1"
  "i > 1 ==> fib i = fib (i - 2) + fib (i - 1)"
  by simp+

lemmas [simp del] = fib.simps

text ‹With precondition›
program_spec fib_prog
  assumes "n ≥ 0" ensures "a = fib n"
  defines ‹
 a = 0; b = 1;
 i = 0;
 while (i < n)
 @variant ‹n₀ - i›
 @invariant ‹n = n₀ ∧ 0 ≤ i ∧ i ≤ n ∧ a = fib i ∧ b = fib (i + 1)›
 {
 c = b;
 b = a + b;
 a = c;
 i = i + 1
 }
 ›
  by vcg_cs (simp add: algebra_simps)

text ‹Without precondition, returning ‹0› for negative numbers›
program_spec fib_prog'
  assumes True ensures "a = fib n₀"
  defines ‹
 a = 0; b = 1;
 i = 0;
 while (i < n)
 @variant ‹n₀ - i›
 @invariant ‹n = n₀ ∧ (0 ≤ i ∧ i ≤ n ∨ n₀ < 0 ∧ i = 0) ∧ a = fib i ∧ b = fib (i+ 1)›
 {
 c = b;
 b = a + b;
 a = c;
 i = i + 1
 }
 ›
  by vcg_cs (auto simp: algebra_simps)



subsubsection ‹Count down›  
text ‹Essentially the same as count up, but we (ab)use the input variable as a counter.›
  
text ‹The invariant is the same as for count-up.
 Only that we have to compute the actual number
 of loop iterations by ‹n₀ - n›. We locally introduce the name ‹c› for that.
 ›  
  

program_spec exp_count_down
  assumes "0≤n"
  ensures "a = 2^nat n₀"
  defines ‹
 a = 1;
 while (n>0)
 @variant ‹n›
 @invariant ‹let c = n₀-n in 0≤n ∧ n≤n₀ ∧ a=2^nat c›
 {
 a=2*a;
 n=n-1
 }
 ›
  apply vcg_cs
  by (auto simp: algebra_simps nat_distribs)

  
subsubsection ‹Approximate from Below›
text ‹Used to invert a monotonic function.
 We count up, until we overshoot the desired result,
 then we subtract one.
 ›  
text ‹The invariant states that the ‹r-1› is not too big.
 When the loop terminates, ‹r-1› is not too big, but ‹r› is already too big,
 so ‹r-1› is the desired value (rounding down).
 ›
text ‹The variant measures the gap that we have to the correct result.
 Note that the loop will do a final iteration, when the result has been reached
 exactly. We account for that by adding one, such that the measure also decreases in this case.
 ›

program_spec sqr_approx_below
  assumes "0≤n"  
  ensures "0≤r ∧ r² ≤ n₀ ∧ n₀ < (r+1)²"
  defines ‹
 r=1;
 while (r*r ≤ n)
 @variant ‹n + 1 - r*r›
 @invariant ‹0≤r ∧ (r-1)² ≤ n₀›
 { r = r + 1 };
 r = r - 1
 ›
  apply vcg
  by (auto simp: algebra_simps power2_eq_square)


subsubsection ‹Bisection›
text ‹A more efficient way of inverting monotonic functions is by bisection,
 that is, one keeps track of a possible interval for the solution, and halfs the
 interval in each step. The program will need $O(\log n)$ iterations, and is
 thus very efficient in practice.

 Although the final algorithm looks quite simple, getting it right can be
 quite tricky.
 ›
text ‹The invariant is surprisingly easy, just stating that the solution
 is in the interval ‹l..<h\›. ›

lemma "∧h l n₀ :: int.
       [🍋''invar-final''; 0 ≤ n₀; ¬ 1 + l < h; 0 ≤ l; l < h; l * l ≤ n₀;
        n₀ < h * h]
       ==> n₀ < 1 + (l * l + l * 2)"  
  by (smt mult.commute semiring_normalization_rules(3))

program_spec sqr_bisect
  assumes "0≤n" ensures "r²≤n₀ ∧ n₀<(r+1)²"
  defines ‹
 l=0; h=n+1;
 while (l+1 < h)
 @variant ‹h-l›
 @invariant ‹0≤l ∧ l<h ∧ l²≤n ∧ n < h²›
 {
 m = (l + h) / 2;
 if (m*m ≤ n) l=m else h=m
 };
 r=l
 ›
  apply vcg
  text ‹We use quick-and-dirty apply style proof to discharge the VCs›
        apply (auto simp: power2_eq_square algebra_simps add_pos_pos)
   apply (smt not_sum_squares_lt_zero)
  by (smt mult.commute semiring_normalization_rules(3))

  
  
subsection ‹Debugging›

subsubsection ‹Testing Programs›

text ‹Stepwise›
schematic_goal "Map.empty: (sqr_approx_below,<''n'':=λ_. 4>) ==> ?s"
  unfolding sqr_approx_below_def
  apply big_step'
   apply big_step'
  apply big_step'
   apply big_step'
    apply big_step'
   apply big_step'
    apply big_step'
   apply big_step'
  apply big_step'
  done

text ‹Or all steps at once›
schematic_goal "Map.empty: (sqr_bisect,<''n'':=λ_. 4900000001>) ==> ?s"
  unfolding sqr_bisect_def
  by big_step
  
  
subsection ‹More Numeric Algorithms›  

subsubsection ‹Euclid's Algorithm (with subtraction)›

(* Crucial Lemmas *)
thm gcd.commute gcd_diff1

program_spec euclid1
  assumes "a>0 ∧ b>0"
  ensures "a = gcd a₀ b₀"
  defines ‹
 while (a≠b)
 @invariant ‹gcd a b = gcd a₀ b₀ ∧ (a>0 ∧ b>0)›
 @variant ‹a+b›
 {
 if (a<b) b = b-a
 else a = a-b
 }
 ›
  apply vcg_cs
   apply (metis gcd.commute gcd_diff1)
  apply (metis gcd_diff1)
  done


subsubsection ‹Euclid's Algorithm (with mod)›

(* Crucial Lemmas *)
thm gcd_red_int[symmetric]

program_spec euclid2
  assumes "a>0 ∧ b>0"
  ensures "a = gcd a₀ b₀"
  defines ‹
 while (b≠0)
 @invariant ‹gcd a b = gcd a₀ b₀ ∧ b≥0 ∧ a>0›
 @variant ‹b›
 {
 t = a;
 a = b;
 b = t mod b
 }
 ›
  apply vcg_cs
  by (simp add: gcd_red_int[symmetric])
  
subsubsection ‹Extended Euclid's Algorithm›
(* Provided by Simon Wimmer. *)

locale extended_euclid_aux_lemmas begin

lemma aux2:
  fixes a b :: int
  assumes "b = t * b₀ + s * a₀" "q = a div b" "gcd a b = gcd a₀ b₀"
  shows "gcd b (a - (a₀ * (s * q) + b₀ * (t * q))) = gcd a₀ b₀"
proof -
  have "a - (a₀ * (s * q) + b₀ * (t * q)) = a - b * q"
    unfolding ‹b = _› by (simp add: algebra_simps)
  also have "a - b * q = a mod b"
    unfolding ‹q = _› by (simp add: algebra_simps)
  finally show ?thesis
    unfolding ‹gcd a b = _›[symmetric] by (simp add: gcd_red_int[symmetric])
qed

lemma aux3:
  fixes a b :: int
  assumes "b = t * b₀ + s * a₀" "q = a div b" "b > 0"
  shows "t * (b₀ * q) + s * (a₀ * q) ≤ a"
proof -
  have "t * (b₀ * q) + s * (a₀ * q) = q * b"
    unfolding ‹b = _› by (simp add: algebra_simps)
  then show ?thesis
    using ‹b > 0›
    by (simp add: algebra_simps ‹q = a div b›)
      (simp flip: minus_mod_eq_mult_div)
qed

end


text ‹The following is a direct translation of the pseudocode for the Extended Euclidean algorithm
  described by the English version of Wikipedia
 🪙‹https://en.wikipedia.org/wiki/Extended_Euclidean_algorithm›):
 ›
program_spec euclid_extended
  assumes "a>0 ∧ b>0"
  ensures "old_r = gcd a₀ b₀ ∧ gcd a₀ b₀ = a₀ * old_s + b₀ * old_t"
defines ‹
 s = 0; old_s = 1;
 t = 1; old_t = 0;
 r = b; old_r = a;
 while (r≠0)
 @invariant ‹
 gcd old_r r = gcd a₀ b₀ ∧ r≥0 ∧ old_r>0
 ∧ a₀ * old_s + b₀ * old_t = old_r ∧ a₀ * s + b₀ * t = r
 ›
 @variant ‹r›
 {
 quotient = old_r / r;
 temp = old_r;
 old_r = r;
 r = temp - quotient * r;
 temp = old_s;
 old_s = s;
 s = temp - quotient * s;
 temp = old_t;
 old_t = t;
 t = temp - quotient * t
 }
 ›
proof -
  interpret extended_euclid_aux_lemmas .
  show ?thesis
    apply vcg_cs
     apply (simp add: algebra_simps)
     apply (simp add: aux2 aux3 minus_div_mult_eq_mod)+
    done  
    
qed


text ‹Non-Wikipedia version›
context extended_euclid_aux_lemmas begin
  lemma aux:
    fixes a b q x y:: int
    assumes "a = old_y * b₀ + old_x * a₀" "b = y * b₀ + x * a₀" "q = a div b"
    shows
    "a mod b + (a₀ * (x * q) + b₀ * (y * q)) = a"
  proof -
    have *: "a₀ * (x * q) + b₀ * (y * q) = q * b"
      unfolding ‹b = _› by (simp add: algebra_simps)
    show ?thesis
      unfolding * unfolding ‹q = _› by simp
  qed
end

program_spec euclid_extended'
  assumes "a>0 ∧ b>0"
  ensures "a = gcd a₀ b₀ ∧ gcd a₀ b₀ = a₀ * x + b₀ * y"
defines ‹
 x = 0;
 y = 1;
 old_x = 1;
 old_y = 0;
 while (b≠0)
 @invariant ‹
 gcd a b = gcd a₀ b₀ ∧ b≥0 ∧ a>0 ∧ a = a₀ * old_x + b₀ * old_y ∧ b = a₀ * x + b₀ * y
 ›
 @variant ‹b›
 {
 q = a / b;
 t = a;
 a = b;
 b = t mod b;
 t = x;
 x = old_x - q * x;
 old_x = t;
 t = y;
 y = old_y - q * y;
 old_y = t
 };
 x = old_x;
 y = old_y
 ›
proof -
  interpret extended_euclid_aux_lemmas .
  show ?thesis
    apply vcg_cs
    apply (simp add: gcd_red_int[symmetric])
    apply (simp add: algebra_simps)
    apply (rule aux; simp add: algebra_simps)
    done
qed
  
subsubsection ‹Exponentiation by Squaring›
(* Contributed by Simon Wimmer *)

lemma ex_by_sq_aux:
  fixes x :: int and n :: nat
  assumes "n mod 2 = 1"
  shows "x * (x * x) ^ (n div 2) = x ^ n"
proof -
  have "n > 0"
    using assms by presburger
  have "2 * (n div 2) = n - 1"
    using assms by presburger
  then have "(x * x) ^ (n div 2) = x ^ (n - 1)"
    by (simp add: semiring_normalization_rules)
  with ‹0 < n› show ?thesis
    by simp (metis Suc_pred power.simps(2))
qed

text ‹A classic algorithm for computing ‹xⁿ› works by repeated squaring,
  the following recurrence:
 ▪ ‹xⁿ = x * x^n-1)/22› if ‹n› is odd
 ▪ ‹xⁿ = x^/22› if ‹n› is even

 ›
program_spec ex_by_sq
  assumes "n≥0"
  ensures "r = x₀ ^ nat n₀"
defines ‹
 r = 1;
 while (n≠0)
 @invariant ‹
 n ≥ 0 ∧ r * x ^ nat n = x₀ ^ nat n₀
 ›
 @variant ‹n›
 {
 if (n mod 2 == 1) {
 r = r * x
 };
 x = x * x;
 n = n / 2
 }
 ›
  apply vcg_cs
   apply (auto simp flip: odd_iff_mod_2_eq_one
    simp add: nat_add_distrib nat_mult_distrib power2_eq_square power_mult elim!: oddE evenE)
  done


subsubsection ‹Power-Tower of 2s›

fun tower2 where
  "tower2 0 = 1"
| "tower2 (Suc n) = 2 ^ tower2 n"  

definition "tower2' n = int (tower2 (nat n))"

program_spec tower2_imp
  assumes ‹m>0›
  ensures ‹a = tower2' m₀›
  defines ‹
 a=1;
 while (m>0)
 @variant ‹m›
 @invariant ‹0≤m ∧ m≤m₀ ∧ a = tower2' (m₀-m)›
 {
 n=a;
 
 a = 1;
 while (n>0)
 @variant ‹n›
 @invariant ‹True› ― ‹This will get ugly, there is no ‹n₀› that we could use!›
 {
 a=2*a;
 n=n-1
 };
 
 m=m-1
 }
 ›
  oops

text ‹We prove the inner loop separately instead!
 (It happens to be exactly our const‹exp_count_down› program.)
 ›  
  
program_spec tower2_imp
  assumes ‹m>0›
  ensures ‹a = tower2' m₀›
  defines ‹
 a=1;
 while (m>0)
 @variant ‹m›
 @invariant ‹0≤m ∧ m≤m₀ ∧ a = tower2' (m₀-m)›
 {
 n=a;
 inline exp_count_down;
 m=m-1
 }
 ›
  apply vcg_cs
  by (auto simp: algebra_simps tower2'_def nat_distribs)



subsection ‹Array Algorithms›  
  
(*
  Presentation in lecture:

  introduce array syntax and semantics
    CONVENTIONS: 
      * Every variable is an array. 
      * Arrays are modeled as int\<Rightarrow>int. Negative indexes allowed.
      * Syntactic sugar to use it at index 0 by default
      
    V vname \<mapsto> Vidx vname aexp  Syntax: x[i]
    abbreviation "V x = Vidx x 0"
  
    Assign vname aexp \<mapsto> AssignIdx vname aexp aexp  Syntax: x[i] = a
    abbreviation: Assign x a = AssignIdx x 0 a
  
    NEW COMMAND: ArrayCpy vname vname Syntax: x[] = y
      copy two arrays as a whole. Assign x (V y) only copies index 0
  
  semantics: Obvious
    aval (Vidx x i) s = s x (aval i s)
    
    (x[i] = a, s) \<Rightarrow> s(x := (s x)(aval i s := aval a s))
    
    (x[] = y, s) \<Rightarrow> s(x := s y)
    
  weakest preconditions: Also obvious
    \<dots>
    
Reasoning with arrays:    
    
  Usually: Model as int\<Rightarrow>int is appropriate.
  
  Common idioms: 
    * {l..<h} the set of indexes i, l\<le>i<h
    * a`{l..<h} the set of values in range {l..<h}
  
    * ran_sorted a l h = \<dots>
    
    Warning: Sledgehammer does not perform well on these!

Data Refinement    
  Sometimes, abstraction to list may be useful.
  
  Idea: Model an algorithm on lists first. Prove it correct.
    Then implement algorithm. Show that it implements the model.
    By transitivity, get: Algorithm is correct!

*)


  
  
subsubsection ‹Summation›

program_spec array_sum
  assumes "l≤h"
  ensures "r = (∑i=l₀..<h₀. a₀ i)"
  defines ‹
 r = 0;
 while (l<h)
 @invariant ‹l₀≤l ∧ l≤h ∧ r = (∑i=l₀..<l. a₀ i)›
 @variant ‹h-l›
 {
 r = r + a[l];
 l=l+1
 }
 ›
  apply vcg_cs
  by (auto simp: intvs_incdec)
  
(* Exercise: Remove l\<le>h precondition! *)  

subsubsection ‹Finding Least Index of Element›  
  
program_spec find_least_idx
  assumes ‹l≤h›
  ensures ‹if l=h₀ then x₀ ∉ a₀`{l<sub>0</sub>..<h\<sub>0</sub>} else l∈{l<sub>0</sub>..<h\<sub>0</sub>} ∧ a<sub>0</sub> l = x<sub>0</sub> ∧ x<sub>0</sub>∉a<sub>0</sub>`{l₀..<l} ›
  defines ‹
 while (l<h ∧ a[l] ≠ x)
 @variant ‹h-l›
 @invariant ‹l₀≤l ∧ l≤h ∧ x∉a`{l₀..<l}›
 l=l+1
 ›
  apply vcg_cs
  by (smt atLeastLessThan_iff imageI)
  


  

subsubsection ‹Check for Sortedness›  
  
term ran_sorted

program_spec check_sorted
  assumes ‹l≤h›
  ensures ‹r≠0 ⟷ ran_sorted a₀ l₀ h₀›
  defines ‹
 if (l==h) r=1
 else {
 l=l+1;
 while (l<h ∧ a[l-1] ≤ a[l])
 @variant ‹h-l›
 @invariant ‹l₀<l ∧ l≤h ∧ ran_sorted a l₀ l›
 l=l+1;
 
 if (l==h) r=1 else r=0
 }
 ›
  apply vcg_cs
   apply (auto simp: ran_sorted_def) 
  by (smt atLeastLessThan_iff) 

  
subsubsection ‹Find Equilibrium Index›  
  
definition "is_equil a l h i ≡ l≤i ∧ i<h ∧ (∑j=l..<i. a j) = (∑j=i..<h. a j)"

program_spec equilibrium
  assumes ‹l≤h›
  ensures ‹is_equil a l h i ∨ i=h ∧ (∄i. is_equil a l h i)›
  defines ‹
 usum=0; i=l;
 while (i<h)
 @variant ‹h-i›
 @invariant ‹l≤i ∧ i≤h ∧ usum = (∑j=l..<i. a j)›
 {
 usum = usum + a[i]; i=i+1
 };
 i=l; lsum=0;
 while (usum ≠ lsum ∧ i<h)
 @variant ‹h-i›
 @invariant ‹l≤i ∧ i≤h
 ∧ lsum=(∑j=l..<i. a j)
 ∧ usum=(∑j=i..<h. a j)
 ∧ (∀j<i. ¬is_equil a l h j)
 ›
 {
 lsum = lsum + a[i];
 usum = usum - a[i];
 i=i+1
 }
 ›
  apply vcg_cs
     apply (auto simp: intvs_incdec is_equil_def)
    apply (metis atLeastLessThan_iff eq_iff finite_atLeastLessThan_int sum_diff1)
   apply force
  by force

  
subsubsection ‹Rotate Right›
program_spec rotate_right
  assumes "0<n"
  ensures "∀i∈{0..<n}. a i = a₀ ((i-1) mod n)"
  defines ‹
 i = 0;
 prev = a[n - 1];
 while (i < n)
 @invariant ‹
 0 ≤ i ∧ i ≤ n
 ∧ (∀j∈{0..<i}. a j = a₀ ((j-1) mod n))
 ∧ (∀j∈{i..<n}. a j = a₀ j)
 ∧ prev = a₀ ((i-1) mod n)
 ›
 @variant ‹n - i›
 {
 temp = a[i];
 a[i] = prev;
 prev = temp;
 i = i + 1
 }
 ›
  apply vcg_cs
  by (simp add: zmod_minus1)

  
subsubsection ‹Binary Search, Leftmost Element›  
  
text ‹We first specify the pre- and postcondition›
definition "bin_search_pre a l h ⟷ l≤h ∧ ran_sorted a l h"

definition "bin_search_post a l h x i ⟷
  l≤i ∧ i≤h ∧ (∀i∈{l..<i}. a i < x) ∧ (∀i∈{i..<h}. x ≤ a i)"    
  
text ‹Then we prove that the program is correct›  
program_spec binsearch 
  assumes ‹bin_search_pre a l h›
  ensures ‹bin_search_post a₀ l₀ h₀ x₀ l›
  defines ‹
 while (l < h)
 @variant ‹h-l›
 @invariant ‹l₀≤l ∧ l≤h ∧ h≤h₀ ∧ (∀i∈{l₀..<l}. a i < x) ∧ (∀i∈{h..<h\₀}. x ≤ a i)›
 {
 m = (l + h) / 2;
 if (a[m] < x) l = m + 1
 else h = m
 }
 ›
  apply vcg_cs
       apply (auto simp: algebra_simps bin_search_pre_def bin_search_post_def)
  txt ‹Driving sledgehammer to its limits ...›
   apply (smt div_add_self1 even_succ_div_two odd_two_times_div_two_succ ran_sorted_alt)
  by (smt div_add_self1 even_succ_div_two odd_two_times_div_two_succ ran_sorted_alt)

text ‹Next, we show that our postcondition (which was easy to prove)
 implies the expected properties of the algorithm.
 ›
lemma
  assumes "bin_search_pre a l h" "bin_search_post a l h x i"
  shows bin_search_decide_membership: "x∈a`{l..<h} ⟷ (i<h ∧ x = a i)"
    and bin_search_leftmost: "x∉a`{l..<i}"
  using assms apply -
   apply (auto simp: bin_search_post_def bin_search_pre_def)
  apply (smt atLeastLessThan_iff ran_sorted_alt)
  done  
  
subsubsection ‹Naive String Search›  
(* By Simon Wimmer *)

program_spec match_string
  assumes "l1 ≤ h1"
  ensures "(∀j ∈ {0..<i}. a (l + j) = b (l1 + j)) ∧ (i < h1 - l1 ⟶ a (l + i) ≠ b (l1 + i))
    ∧ 0 ≤ i ∧ i ≤ h1 - l1"
  defines ‹
 i = 0;
 while (l1 + i < h1 ∧ a[l + i] == b[l1 + i])
 @invariant ‹(∀j ∈ {0..<i}. a (l + j) = b (l1 + j)) ∧ 0 ≤ i ∧ i ≤ h1 - l1›
 @variant ‹(h1 - (l1 + i))›
 {
 i = i + 1
 }
 ›
  by vcg_cs auto

lemma lran_eq_iff': "lran a l1 (l1 + (h - l)) = lran a' l h
  ⟷ (∀i. 0 ≤ i ∧ i < h - l ⟶ a (l1 + i) = a' (l + i))" if "l ≤ h"
  using that
proof (induction "nat (h - l)" arbitrary: h)
  case 0
  then show ?case
    by auto
next
  case (Suc x)
  then have *: "x = nat (h - 1 - l)" "l ≤ h - 1"
    by auto
  note IH = Suc.hyps(1)[OF *]
  from * have 1:
    "lran a l1 (l1 + (h - l)) = lran a l1 (l1 + (h - 1 - l)) @ [a (l1 + (h - 1 - l))]"
    "lran a' l h = lran a' l (h - 1) @ [a' (h - 1)]"
    by (auto simp: lran_bwd_simp algebra_simps lran_butlast[symmetric])
  from IH * Suc.hyps(2) Suc.prems show ?case
    unfolding 1
    apply auto
    subgoal for i
      by (cases "i = h - 1 - l") auto
    done
qed

program_spec match_string'
  assumes "l1 ≤ h1"
ensures "i = h1 - l1 ⟷ lran a l (l + (h1 - l1)) = lran b l1 h1"
for i h1 l1 l a[] b[]
defines ‹inline match_string›
  by vcg_cs (auto simp: lran_eq_iff')

program_spec substring
  assumes "l ≤ h ∧ l1 ≤ h1"
  ensures "match = 1 ⟷ (∃ j ∈ {l₀..<h₀}. lran a j (j + (h1 - l1)) = lran b l1 h1)"
  for a[] b[]
  defines ‹
 match = 0;
 while (l < h ∧ match == 0)
 @invariant‹l₀ ≤ l ∧ l ≤ h ∧ match ∈ {0,1} ∧
 (if match = 1
 then lran a l (l + (h1 - l1)) = lran b l1 h1 ∧ l < h
 else (∀j ∈ {l₀..<l}. lran a j (j + (h1 - l1)) ≠ lran b l1 h1))›
 @variant‹(h - l) * (1 - match)›
 {
 inline match_string';
 if (i == h1 - l1) {match = 1}
 else {l = l + 1}
 }
 ›
  by vcg_cs auto

program_spec substring'
  assumes "l ≤ h ∧ l1 ≤ h1"
  ensures "match = 1 ⟷ (∃ j ∈ {l₀..h₀-(h1 - l1)}. lran a j (j + (h1 - l1)) = lran b l1 h1)"
  for a[] b[]
  defines ‹
 match = 0;
 if (l + (h1 - l1) ≤ h) {
 h = h - (h1 - l1) + 1;
 inline substring
 }
 ›
  by vcg_cs auto

(* Or combined: *)
program_spec substring''
  assumes "l ≤ h ∧ l1 ≤ h1"
  ensures "match = 1 ⟷ (∃ j ∈ {l₀..<h₀-(h1 - l1)}. lran a j (j + (h1 - l1)) = lran b l1 h1)"
  for a[] b[]
  defines ‹
 match = 0;
 if (l + (h1 - l1) ≤ h) {
 while (l + (h1 - l1) < h ∧ match == 0)
 @invariant‹l₀ ≤ l ∧ l ≤ h - (h1 - l1) ∧ match ∈ {0,1} ∧
 (if match = 1
 then lran a l (l + (h1 - l1)) = lran b l1 h1 ∧ l < h - (h1 - l1)
 else (∀j ∈ {l₀..<l}. lran a j (j + (h1 - l1)) ≠ lran b l1 h1))›
 @variant‹(h - l) * (1 - match)›
 {
 inline match_string';
 if (i == h1 - l1) {match = 1}
 else {l = l + 1}
 }
 }
 ›
  by vcg_cs auto

lemma lran_split:
  "lran a l h = lran a l p @ lran a p h" if "l ≤ p" "p ≤ h"
  using that by (induction p; simp; simp add: lran.simps)

lemma lran_eq_append_iff:
  "lran a l h = as @ bs ⟷ (∃ i. l ≤ i ∧ i ≤ h ∧ as = lran a l i ∧ bs = lran a i h)" if "l ≤ h"
  apply safe
  subgoal
    using that
  proof (induction as arbitrary: l)
    case Nil
    then show ?case
      by auto
  next
    case (Cons x as)
    from this(2-) have "lran a (l + 1) h = as @ bs" "a l = x" "l + 1 ≤ h"
        apply -
      subgoal
        by simp
      subgoal
        by (smt append_Cons list.inject lran.simps lran_append1)
      subgoal
        using add1_zle_eq by fastforce
      done
    from Cons.IH[OF this(1,3)] obtain i
      where IH: "l + 1 ≤ i" "i ≤ h" "as = lran a (l + 1) i" "bs = lran a i h"
      by auto
    with ‹a l = x› show ?case
      apply (intro exI[where x = i])
      apply auto
      by (smt IH(3) lran_prepend1)
  qed
  apply (subst lran_split; simp)
  done

lemma lran_split':
  "(∃j∈{l..h - (h1 - l1)}. lran a j (j + (h1 - l1)) = lran b l1 h1)
 = (∃as bs. lran a l h = as @ lran b l1 h1 @ bs)" if "l ≤ h" "l1 ≤ h1"
proof safe
  fix j
  assume j: "j ∈ {l..h - (h1 - l1)}" and match: "lran a j (j + (h1 - l1)) = lran b l1 h1"
  with ‹l1 ≤ h1› have "lran a l h = lran a l j @ lran a j (j + (h1 - l1)) @ lran a (j + (h1 - l1)) h"
    apply (subst lran_split[where p = j], simp, simp)
    apply (subst (2) lran_split[where p = "j + (h1 - l1)"]; simp)
    done
  then show "∃as bs. lran a l h = as @ lran b l1 h1 @ bs"
    by (auto simp: match)
next
  fix as bs
  assume "lran a l h = as @ lran b l1 h1 @ bs"
  with that lran_eq_append_iff[of l h a as "lran b l1 h1 @ bs"] obtain i where
    "as = lran a l i" "lran a i h = lran b l1 h1 @ bs" "l ≤ i" "i ≤ h"
    by auto
  with lran_eq_append_iff[of i h a "lran b l1 h1" bs] obtain j where j:
    "lran b l1 h1 = lran a i j" "bs = lran a j h" "i ≤ j" "j ≤ h"
    by auto
  moreover have "j = i + (h1 - l1)"
  proof -
    have "length (lran b l1 h1) = nat (h1 - l1)" "length (lran a i j) = nat (j - i)"
      by auto
    with j(1,3) that show ?thesis
      by auto
  qed
  ultimately show "∃j∈{l..h - (h1 - l1)}. lran a j (j + (h1 - l1)) = lran b l1 h1"
    using ‹l ≤ i› by auto
qed

program_spec substring_final
  assumes "l ≤ h ∧ 0 ≤ len"
  ensures "match = 1 ⟷ (∃as bs. lran a l₀ h₀ = as @ lran b 0 len @ bs)"
  for l h len match a[] b[]
  defines ‹l1 = 0; h1 = len; inline substring'›
  supply [simp] = lran_split'[symmetric]
  apply vcg_cs
  done


  
  
subsubsection ‹Insertion Sort›
    
text ‹
 We first prove the inner loop.
 The specification here specifies what the algorithm does as closely as possible,
 such that it becomes easier to prove. In this case, sortedness is not a precondition
 for the inner loop to move the key element backwards over all greater elements.
 ›
definition "insort_insert_post l j a₀ a i
  ⟷ (let key = a₀ j in
      i∈{l-1..<j} ― ‹‹i› is in range›
    ― ‹Content of new array›
    ∧ (∀k∈{l..i}. a k = a₀ k)
    ∧ a (i+1) = key
    ∧ (∀k∈{i+2..j}. a k = a₀ (k-1))
    ∧ a = a₀ on -{l..j}
    ― ‹Placement of ‹key››
    ∧ (i≥l ⟶ a i ≤ key) ― ‹Element at ‹i› smaller than ‹key›, if it exists ›
    ∧ (∀k∈{i+2..j}. a k > key) ― ‹Elements ‹≥i+2› greater than ‹key››
    )
    "
  for l j i :: int and a₀ a :: "int ==> int"
  
program_spec insort_insert 
  assumes "l<j"
  ensures "insort_insert_post l j a₀ a i"
  defines ‹
 key = a[j];
 i = j-1;
 while (i≥l ∧ a[i]>key)
 @variant ‹i-l+1›
 @invariant ‹l-1≤i ∧ i<j 
 ∧ (∀k∈{l..i}. a k = a₀ k)
 ∧ (∀k∈{i+2..j}. a k > key ∧ a k = a₀ (k-1))
 ∧ a = a₀ on -{l..j}
 ›
 {
 a[i+1] = a[i];
 i=i-1
 };
 a[i+1] = key
 ›
  unfolding insort_insert_post_def Let_def
  apply vcg
      apply auto
  by (smt atLeastAtMost_iff)
  
text ‹Next, we show that our specification that was easy to prove implies the
 specification that the outer loop expects:›

text ‹Invoking ‹insort_insert› will sort in the element›
lemma insort_insert_sorted:
  assumes "l<j"
  assumes "insort_insert_post l j a a' i'"
  assumes "ran_sorted a l j"
  shows "ran_sorted a' l (j + 1)"
  using assms unfolding insort_insert_post_def Let_def
  apply auto
  subgoal
    by (smt atLeastAtMost_iff ran_sorted_alt)
  subgoal 
    unfolding ran_sorted_alt Ball_def
    apply auto
    by smt
  done  
  
text ‹Invoking ‹insort_insert› will only mutate the elements›
lemma insort_insert_ran1:
  assumes "l<j"
  assumes "insort_insert_post l j a a' i"
  shows "mset_ran a' {l..j} = mset_ran a {l..j}"
proof -
  from assms have EQS:
    "a' = a on {l..i}"
    "a' (i+1) = a j"
    "a' = (a o (+)(-1)) on {i+2..j}"
    unfolding insort_insert_post_def eq_on_def Let_def by auto

  from assms have "l≤i+1" "i+1≤j" unfolding insort_insert_post_def Let_def by auto
    
  have ranshift: "mset_ran (a o (+)(-1)) {i+2..j} = mset_ran a {i+1..j-1}"  
    by (simp add: mset_ran_shift algebra_simps)
    
      
  have "mset_ran a' {l..j} = mset_ran a' {l..i} + {# a' (i+1) #} + mset_ran a' {i+2..j}"  
    using ‹l<j\› ‹l≤i+1› ‹i+1≤j›
    apply (simp add: mset_ran_combine)
    by (auto intro: arg_cong[where f="mset_ran a'"])
  also have "… = mset_ran a {l..i} + {# a j #} + mset_ran (a o (+)(-1)) {i+2..j}"
    using EQS(1,3)[THEN mset_ran_cong] EQS(2) by simp
  also have "… = mset_ran a {l..i} + mset_ran a {i+1..j-1} + {# a j #}"
    by (simp add: mset_ran_shift algebra_simps) 
  also have "… = mset_ran a {l..j}"
    using ‹l<j\› ‹l≤i+1› ‹i+1≤j›
    apply (simp add: mset_ran_combine)
    by (auto intro: arg_cong[where f="mset_ran a"])
  finally show ?thesis .
qed  

text ‹The property @{thm insort_insert_ran1} extends to the whole array to be sorted›
lemma insort_insert_ran2:
  assumes "l<j" "j<h"
  assumes "insort_insert_post l j a a' i"
  shows "mset_ran a' {l..<h} = mset_ran a {l..<h}" (is ?thesis1)
    and "a'=a on -{l..<h}" (is ?thesis2)
proof -
  from insort_insert_ran1 assms have "mset_ran a' {l..j} = mset_ran a {l..j}" by blast
  also from ‹insort_insert_post l j a a' i› have "a' = a on {j<..<h}"
    unfolding insort_insert_post_def Let_def by (auto simp: eq_on_def)
  hence "mset_ran a' {j<..<h} = mset_ran a {j<..<h}" by (rule mset_ran_cong)
  finally (mset_ran_combine_eqs) show ?thesis1
    by (simp add: assms ivl_disj_int_two(4) ivl_disj_un_two(4) le_less)
    
  from assms show ?thesis2   
    unfolding insort_insert_post_def Let_def eq_on_def 
    by auto
    
qed  

text ‹Finally, we specify and prove correct the outer loop›
program_spec insort
  assumes "l<h" 
  ensures "ran_sorted a l h ∧ mset_ran a {l..<h} = mset_ran a₀ {l..<h} ∧ a=a₀ on -{l..<h}"
  for a[]
  defines ‹
 j = l + 1;
 while (j<h)
 @variant ‹h-j›
 @invariant ‹
 l<j ∧ j≤h ― ‹‹j› in range›
 ∧ ran_sorted a l j ― ‹Array is sorted up to ‹j››
 ∧ mset_ran a {l..<h} = mset_ran a₀ {l..<h} ― ‹Elements in range only permuted›
 ∧ a=a₀ on -{l..<h}
 ›
 {
 inline insort_insert;
 j=j+1
 }
 ›
  apply vcg_cs
  apply (intro conjI)
  subgoal by (rule insort_insert_sorted)
  subgoal using insort_insert_ran2(1) by auto
  subgoal apply (frule (2) insort_insert_ran2(2)) by (auto simp: eq_on_def)
  done    
  

  
subsubsection ‹Quicksort›

procedure_spec partition_aux(a,l,h,p) returns (a,i)
  assumes "l≤h"
  ensures "mset_ran a₀ {l₀..<h₀} = mset_ran a {l₀..<h₀}
         ∧ (∀j∈{l₀..<i}. a j < p₀)
         ∧ (∀j∈{i..<h₀}. a j ≥ p₀)
         ∧ l₀≤i ∧ i≤h₀
         ∧ a₀ = a on -{l₀..<h₀}
         "
  defines ‹
 i=l; j=l;
 while (j<h)
 @invariant ‹
 l≤i ∧ i≤j ∧ j≤h
 ∧ mset_ran a₀ {l₀..<h\₀} = mset_ran a {l₀..<h\₀}
 ∧ (∀k∈{l..<i}. a k < p)
 ∧ (∀k∈{i..<j}. a k ≥ p)
 ∧ (∀k∈{j..<h}. a₀ k = a k)
 ∧ a₀ = a on -{l₀..<h\₀}
 ›
 @variant ‹(h-j)›
 {
 if (a[j]<p) {temp = a[i]; a[i] = a[j]; a[j] = temp; i=i+1};
 j=j+1
 }
 ›
  supply lran_eq_iff[simp] lran_tail[simp del]
  apply vcg_cs
  subgoal by (simp add: mset_ran_swap[unfolded swap_def])
  subgoal by auto
  done
  

procedure_spec partition(a,l,h,p) returns (a,i)
  assumes "l<h"
  ensures "mset_ran a₀ {l₀..<h₀} = mset_ran a {l₀..<h₀}
         ∧ (∀j∈{l₀..<i}. a j < a i)
         ∧ (∀j∈{i..<h₀}. a j ≥ a i)
         ∧ l₀≤i ∧ i<h₀ ∧ a₀ (h₀-1) = a i
         ∧ a₀ = a on -{l₀..<h₀}
         "
  defines ‹
 p = a[h-1];
 (a,i) = partition_aux(a,l,h-1,p);
 a[h-1] = a[i];
 a[i] = p
 ›
  apply vcg_cs
  apply (auto simp: eq_on_def mset_ran_swap[unfolded swap_def] intvs_incdec intro: mset_ran_combine_eq_diff)
  done
    
  
  
lemma quicksort_sorted_aux:  
  assumes BOUNDS: "l ≤ i" "i < h"
  
  assumes LESS: "∀j∈{l..<i}. a₁ j < a₁ i" 
  assumes GEQ: "∀j∈{i..<h}. a₁ i ≤ a₁ j"
  
  assumes R1: "mset_ran a₁ {l..<i} = mset_ran a₂ {l..<i}"
  assumes E1: "a₁ = a₂ on - {l..<i}"
  
  assumes SL: "ran_sorted a₂ l i"
  
  assumes R2: "mset_ran a₂ {i + 1..<h} = mset_ran a₃ {i + 1..<h}"
  assumes E2: "a₂ = a₃ on - {i + 1..<h}"

  assumes SH: "ran_sorted a₃ (i + 1) h"
  
  shows "ran_sorted a₃ l h"
proof -
  have [simp]: "{l..<i} ⊆ - {i + 1..<h}" by auto
  have [simp]: "a₁ i = a₃ i" using E1 E2 by (auto simp: eq_on_def)
  
  note X1 = mset_ran_xfer_pointwise[where P = ‹λx. x < p› for p, OF R1, simplified]
  note X2 = eq_on_xfer_pointwise[where P = ‹λx. x < p› for p, OF E2, of "{l..<i}", simplified]
  from LESS have LESS': "∀j∈{l..<i}. a₃ j < a₃ i"
    by (simp add: X1 X2)
    
  
  from GEQ have GEQ1: "∀j∈{i+1..<h}. a₁ i ≤ a₁ j" by auto
  have [simp]: "{i + 1..<h} ⊆ - {l..<i}" by auto
  note X3 = eq_on_xfer_pointwise[where P = ‹λx. x ≥ p› for p, OF E1, of "{i+1..<h}", simplified]
  note X4 = mset_ran_xfer_pointwise[where P = ‹λx. x ≥ p› for p, OF R2, simplified]  
  from GEQ1 have GEQ': "∀j∈{i+1..<h}. a₃ i ≤ a₃ j" by (simp add: X3 X4)

  
  from SL eq_on_xfer_ran_sorted[OF E2, of l i] have SL': "ran_sorted a₃ l i" by simp
  
  show ?thesis using combine_sorted_pivot[OF BOUNDS SL' SH LESS' GEQ'] .
qed  
  
lemma quicksort_mset_aux:
  assumes B: "l₀≤i" "i<h₀"
  assumes R1: "mset_ran a {l₀..<i} = mset_ran aa {l₀..<i}"
  assumes E1: "a = aa on - {l₀..<i}"
  assumes R2: "mset_ran aa {i + 1..<h₀} = mset_ran ab {i + 1..<h₀}"
  assumes E2: "aa = ab on - {i + 1..<h₀}"
  shows "mset_ran a {l₀..<h₀} = mset_ran ab {l₀..<h₀}"
  apply (rule trans)
   apply (rule mset_ran_eq_extend[OF R1 E1])
  using B apply auto [2]
  apply (rule mset_ran_eq_extend[OF R2 E2])
  using B apply auto [2]
  done 
  
  
recursive_spec quicksort(a,l,h) returns a 
  assumes "True"
  ensures "ran_sorted a l₀ h₀ ∧ mset_ran a₀ {l₀..<h₀} = mset_ran a {l₀..<h₀} ∧ a₀=a on -{l₀..<h₀}"
  variant "h-l"
  defines ‹
 if (l<h) {
 (a,i) = partition(a,l,h,a[l]);
 a = rec quicksort(a,l,i);
 a = rec quicksort(a,i+1,h)
 }
 ›
  apply (vcg_cs; (intro conjI)?)
  subgoal using quicksort_sorted_aux by metis
  subgoal using quicksort_mset_aux by metis
  subgoal by (smt ComplD ComplI atLeastLessThan_iff eq_on_def)
  subgoal by (auto simp: ran_sorted_def)
  done
  
  
subsection ‹Data Refinement›

subsubsection ‹Filtering›

program_spec array_filter_negative
  assumes "l≤h"
  ensures "lran a l₀ i = filter (λx. x≥0) (lran a₀ l₀ h₀)"
  defines ‹
 i=l; j=l;
 while (j<h)
 @invariant ‹
 l≤i ∧ i≤j ∧ j≤h
 ∧ lran a l i = filter (λx. x≥0) (lran a₀ l j)
 ∧ lran a j h = lran a₀ j h
 ›
 @variant ‹h-j›
 {
 if (a[j]≥0) {a[i] = a[j]; i=i+1};
 j=j+1
 }
 ›
  supply lran_eq_iff[simp] lran_tail[simp del]
  apply vcg_cs
  done
  
  (* Exercise: Use swap to modify this algorithm to partitioning! *)
  

  
subsubsection ‹Merge Two Sorted Lists›

text ‹
 We define the merge function abstractly first, as a functional program on lists.
 ›

fun merge where
  "merge [] ys = ys"
| "merge xs [] = xs"
| "merge (x#xs) (y#ys) = (if x<y then x#merge xs (y#ys) else y#merge (x#xs) ys)"   

lemma merge_add_simp[simp]: "merge xs [] = xs" by (cases xs) auto

text ‹It's straightforward to show that this produces a sorted list with the same elements.›
lemma merge_sorted:
  assumes "sorted xs" "sorted ys" 
  shows "sorted (merge xs ys) ∧ set (merge xs ys) = set xs ∪ set ys"
  using assms
  apply (induction xs ys rule: merge.induct)
    apply auto
  done
  
lemma merge_mset: "mset (merge xs ys) = mset xs + mset ys"  
  by (induction xs ys rule: merge.induct) auto
  
  
  
text ‹Next, we prove an equation that characterizes one step of the while loop,
 on the list level.›
lemma merge_eq: "xs≠[] ∨ ys≠[] ==> merge xs ys = (
  if ys=[] ∨ (xs≠[] ∧ hd xs < hd ys) then hd xs # merge (tl xs) ys
  else hd ys # merge xs (tl ys)
  )"
  by (cases xs; cases ys; simp)

text ‹
 We do a first proof that our merge implementation on the arrays and indexes
 behaves like the functional merge on the corresponding lists.
 
 The annotations use the @{const lran} function to map from the implementation level
 to the list level. Moreover, the invariant of the implementation, ‹l≤h›, is carried
 through explicitly.
 ›
program_spec merge_imp' 
  assumes "l1≤h1 ∧ l2≤h2"
  ensures "let ms = lran m 0 j; xs₀ = lran a1₀ l1₀ h1₀; ys₀ = lran a2₀ l2₀ h2₀ in
    j≥0 ∧ ms = merge xs₀ ys₀"
  defines ‹
 j=0;
 while (l1≠h1 ∨ l2≠h2)
 @variant ‹h1 + h2 - l1 - l2›
 @invariant ‹let
 xs=lran a1 l1 h1; ys = lran a2 l2 h2; ms = lran m 0 j;
 xs₀ = lran a1₀ l1₀ h1₀; ys₀ = lran a2₀ l2₀ h2₀
 in
 l1≤h1 ∧ l2≤h2 ∧ 0≤j ∧
 merge xs₀ ys₀ = ms@merge xs ys
 ›
 {
 if (l2==h2 ∨ (l1≠h1 ∧ a1[l1] < a2[l2])) {
 m[j] = a1[l1];
 l1=l1+1
 } else {
 m[j] = a2[l2];
 l2=l2+1
 };
 j=j+1
 }
 ›
  text ‹Given the @{thm [source] merge_eq} theorem, which captures the essence of a loop step,
 and the theorems @{thm lran_append1}, @{thm lran_tail}, and @{thm hd_lran}, which
 convert from the operations on arrays and indexes to operations on lists,
 the proof is straightforward
 ›
  apply vcg_cs
  subgoal apply (subst merge_eq) by auto
  subgoal by linarith
  subgoal apply (subst merge_eq) by auto
  done

text ‹
 In a next step, we refine our proof to
 combine it with the abstract properties we have proved about merge.
 The program does not change (we simply inline the original one here).
 ›
  
procedure_spec merge_imp (a1,l1,h1,a2,l2,h2) returns (m,j)
  assumes "l1≤h1 ∧ l2≤h2 ∧ sorted (lran a1 l1 h1) ∧ sorted (lran a2 l2 h2)"
  ensures "let ms = lran m 0 j in
      j≥0
    ∧ sorted ms
    ∧ mset ms = mset (lran a1₀ l1₀ h1₀) + mset (lran a2₀ l2₀ h2₀)"
  for l1 h1 l2 h2 a1[] a2[] m[] j
  defines ‹inline merge_imp'›
  apply vcg_cs
  apply (auto simp: Let_def merge_mset dest: merge_sorted)
  done
  
thm merge_imp_spec  
thm merge_imp_def
  
lemma [named_ss vcg_bb]:
  "UNIV ∪ a = UNIV"  
  "a ∪ UNIV = UNIV"  
  by auto
  
  
lemma merge_msets_aux: "[l≤m; m≤h] ==> mset (lran a l m) + mset (lran a m h) = mset (lran a l h)"  
  by (auto simp: mset_lran mset_ran_combine ivl_disj_un_two)
  
  
  
  
recursive_spec mergesort (a,l,h) returns (b,j)
  assumes "l≤h"
  ensures ‹0≤j ∧ sorted (lran b 0 j) ∧ mset (lran b 0 j) = mset (lran a₀ l₀ h₀)›
  variant ‹h-l›
  for a[] b[]
  defines ‹
 if (l==h) j=0
 else if (l+1==h) {
 b[0] = a[l];
 j=1
 } else {
 m = (h+l) / 2;
 (a1,h1) = rec mergesort (a,l,m);
 (a2,h2) = rec mergesort (a,m,h);
 (b,j) = merge_imp (a1,0,h1,a2,0,h2)
 }
 ›
  apply vcg
       apply auto []
      apply (auto simp: lran.simps) []
     apply auto []
    apply auto []
   apply auto []
  apply (auto simp: Let_def merge_msets_aux) []
  done
print_theorems    
  
  
subsubsection ‹Remove Duplicates from Array, using Bitvector Set›  
  
text ‹We use an array to represent a set of integers.

 If we only insert elements in range ‹{0..<n}›, this representation
 is called bit-vector (storing a single bit per index is enough).
 ›

definition set_of :: "(int ==> int) ==> int set" where "set_of a ≡ {i. a i ≠ 0}"

context notes [simp] = set_of_def begin
  lemma set_of_empty[simp]: "set_of (λ_. 0) = {}" by auto
  lemma set_of_insert[simp]: "x≠0 ==> set_of (a(i:=x)) = insert i (set_of a)" by auto
  lemma set_of_remove[simp]: "set_of (a(i:=0)) = set_of a - {i}" by auto
  lemma set_of_mem[simp]: "i∈set_of a ⟷ a i ≠ 0" by auto
end

program_spec dedup
  assumes ‹l≤h›
  ensures ‹set (lran a l i) = set (lran a₀ l h) ∧ distinct (lran a l i)›
  defines ‹
 i=l; j=l;
 clear b[];
 while (j<h)
 @variant ‹h-j›
 @invariant ‹l≤i ∧ i≤j ∧ j≤h
 ∧ set (lran a l i) = set (lran a₀ l j)
 ∧ distinct (lran a l i)
 ∧ lran a j h = lran a₀ j h
 ∧ set_of b = set (lran a l i)
 ›
 {
 if (b[a[j]] == 0) {
 a[i] = a[j]; i=i+1; b[a[j]] = 1
 };
 j=j+1
 }
 ›
  apply vcg_cs
   apply (auto simp: lran_eq_iff lran_upd_inside intro: arg_cong[where f=tl]) []
  apply (auto simp: lran_eq_iff) []
  done
  
  
procedure_spec bv_init () returns b 
  assumes True ensures ‹set_of b = {}›
  defines ‹clear b[]›
  by vcg_cs
  
procedure_spec bv_insert (x, b) returns b 
  assumes True ensures ‹set_of b = insert x₀ (set_of b₀)›
  defines ‹b[x] = 1›
  by vcg_cs
  
procedure_spec bv_remove (x, b) returns b
  assumes True ensures ‹set_of b = set_of b₀ - {x₀}›
  defines ‹b[x] = 0›
  by vcg_cs
  
procedure_spec bv_elem (x,b) returns r  
  assumes True ensures ‹r≠0 ⟷ x₀∈set_of b₀›
  defines ‹r = b[x]›
  by vcg_cs
  
  
procedure_spec dedup' (a,l,h) returns (a,l,i)
  assumes ‹l≤h› ensures ‹set (lran a l i) = set (lran a₀ l₀ h₀) ∧ distinct (lran a l i) ›
  for b[]
  defines ‹
 b = bv_init();
 
 i=l; j=l;
 
 while (j<h)
 @variant ‹h-j›
 @invariant ‹l≤i ∧ i≤j ∧ j≤h
 ∧ set (lran a l i) = set (lran a₀ l j)
 ∧ distinct (lran a l i)
 ∧ lran a j h = lran a₀ j h
 ∧ set_of b = set (lran a l i)
 ›
 {
 mem = bv_elem (a[j],b);
 if (mem == 0) {
 a[i] = a[j]; i=i+1; b = bv_insert(a[j],b)
 };
 j=j+1
 }
 ›
  apply vcg_cs
   apply (auto simp: lran_eq_iff lran_upd_inside intro: arg_cong[where f=tl])
  done
  
subsection ‹Recursion›  

subsubsection ‹Recursive Fibonacci›

recursive_spec fib_imp (i) returns r assumes True ensures ‹r = fib i₀› variant ‹i›
  defines ‹
 if (i≤0) r=0
 else if (i==1) r=1
 else {
 r1=rec fib_imp (i-2);
 r2=rec fib_imp (i-1);
 r = r1+r2
 }
 ›
  by vcg_cs    
    
subsubsection ‹Homeier's Cycling Termination›  

text ‹A contrived example from Homeier's thesis.
 Only the termination proof is done.›
(* TODO: Also show correctness: pedal (a,b) will multiply its arguments!
  correctness proof may be a good test how to handle specs with logical vars!
*)  
  
recursive_spec 
  pedal (n,m) returns () assumes ‹n≥0 ∧ m≥0› ensures True variant ‹n+m›
  defines ‹
 if (n≠0 ∧ m≠0) {
 G = G + m;
 if (n<m) rec coast (n-1,m-1) else rec pedal(n-1,m)
 }
 ›
and
  coast (n,m) returns () assumes ‹n≥0 ∧ m≥0› ensures True variant ‹n+m+1›
  defines ‹
 G = G + n;
 if (n<m) rec coast (n,m-1) else rec pedal (n,m)
 ›
  by vcg_cs  
  
  
subsubsection ‹Ackermann›  
  
fun ack :: "nat ==> nat ==> nat" where
  "ack 0 n = n+1"
| "ack m 0 = ack (m-1) 1"
| "ack m n = ack (m-1) (ack m (n-1))"
  
lemma ack_add_simps[simp]: 
  "m≠0 ==> ack m 0 = ack (m-1) 1"
  "[m≠0; n≠0] ==> ack m n = ack (m-1) (ack m (n-1))"
  subgoal by (cases m) auto
  subgoal by (cases "(m,n)" rule: ack.cases) (auto)
  done

recursive_spec relation "less_than <*lex*> less_than"
  ack_imp (m,n) returns r
    assumes "m≥0 ∧ n≥0" ensures "r=int (ack (nat m₀) (nat n₀))"
    variant "(nat m, nat n)" 
    defines ‹
 if (m==0) r = n+1
 else if (n==0) r = rec ack_imp (m-1,1)
 else {
 t = rec ack_imp (m,n-1);
 r = rec ack_imp (m-1,t)
 }
 ›
  supply nat_distribs[simp]
  by vcg_cs
  
  
subsubsection ‹McCarthy's 91 Function›
text ‹A standard benchmark for verification of recursive functions.
 We use Homeier's version with a global variable.›
  
recursive_spec p91(y) assumes True ensures "if 100<y₀ then G = y₀-10 else G = 91" variant "101-y" 
  for G
  defines ‹
 if (100<y) G=y-10
 else {
 rec p91 (y+11);
 rec p91 (G)
 }
 ›
  apply vcg_cs
   apply (auto split: if_splits)
  done
  
subsubsection ‹Odd/Even›  
  
recursive_spec 
  odd_imp (a) returns b 
    assumes True
    ensures "b≠0 ⟷ odd a₀"
    variant "∣a∣"
    defines ‹
 if (a==0) b=0
 else if (a<0) b = rec even_imp (a+1)
 else b = rec even_imp (a-1)
 ›
  and
  even_imp (a) returns b
    assumes True
    ensures "b≠0 ⟷ even a₀"
    variant "∣a∣"
    defines ‹
 if (a==0) b=1
 else if (a<0) b = rec odd_imp (a+1)
 else b = rec odd_imp (a-1)
 ›
  apply vcg  
           apply auto
  done  

thm even_imp_spec

  
subsubsection ‹Pandya and Joseph's Product Producers›
text ‹Again, taking the version from Homeier's thesis, but with a
 modification to also terminate for negative ‹y›.
 ›
  
recursive_spec relation ‹measure nat <*lex*> less_than›
  product () assumes True ensures ‹GZ = GZ₀ + GX₀*GY₀› variant "(∣GY∣,1::nat)" 
  for GX GY GZ
  defines
  ‹
 e = even_imp (GY);
 if (e≠0) rec evenproduct() else rec oddproduct()
 ›
and
  oddproduct() assumes ‹odd GY› ensures ‹GZ = GZ₀ + GX₀*GY₀› variant "(∣GY∣,0::nat)" 
  for GX GY GZ
  defines
  ‹
 if (GY<0) {
 GY = GY + 1;
 GZ = GZ - GX
 } else {
 GY = GY - 1;
 GZ = GZ + GX
 };
 rec evenproduct()
 ›
and
  evenproduct() assumes ‹even GY› ensures ‹GZ = GZ₀ + GX₀*GY₀› variant "(∣GY∣,0::nat)" 
  for GX GY GZ
  defines
  ‹
 if (GY≠0) {
 GX = 2*GX;
 GY = GY / 2;
 rec product()
 }
 ›
  apply vcg_cs
     apply (auto simp: algebra_simps)
   apply presburger+
  done
  
  
(*
  TODO: 
    * Infrastructure to prove different specification for same program
      (DONE) We can use inline, and just give program a new name.

    * ghost variables. Keep them existentially qualified, with naming tag.
        ABS_GHOST ''name'' (\<lambda>x. P x) \<equiv> \<exists>x. P x
        GHOST ''name'' value = True
        
        assumes ABS_GHOST name (\<lambda>x. P x) 
        obtains x where "GHOST name x" "P x"
        
        assumes "GHOST name x" and "P x"
        shows "ABS_GHOST name (\<lambda>x. P x)"
        
        Let_Ghost name (\<lambda>x s. C x s) = skip
        
        assumes "\<exists>x. C x s"
        assumes "\<And>x. \<lbrakk> GHOST name x; C x s \<rbrakk> \<Longrightarrow> Q s"
        shows "wp (Let_Ghost name (\<lambda>x s. C x s)) Q s"
        
      Let specification command also abstract over ghost variables in program.
        
*)  
  
(* IDEAS: 
  
  Extend arrays to multiple dimensions: int list \<Rightarrow> int
    Vidx vname "aexp list"
    AssignIdx vname "aexp list" aexp
    ArrayCpy vname vname  -- copy all dimensions (it's simpler)
    ArrayInit vname       -- Init all dimensions
    
    plain values on dimension []!

    then we can easily encode matrices and adjacency lists!
    
*)  

subsection ‹Graph Algorithms›

subsubsection ‹DFS›
(* By Simon Wimmer *)

text ‹A graph is stored as an array of integers.
 Each node is an index into this array, pointing to a size-prefixed list of successors.

 Example for node ‹i›, which has successors ‹s1... sn›:
 🍋‹
 Indexes: ... | i | i+1 | ... | i+n | ...
 Data: ... | n | s1 | ... | sn | ...
 ›
 ›

definition succs where
  "succs a i ≡ a ` {i+1..<a i}" for a :: "int ==> int"

definition Edges where
  "Edges a ≡ {(i, j). j ∈ succs a i}"

procedure_spec push' (x, stack, ptr) returns (stack, ptr)
  assumes "ptr ≥ 0" ensures ‹lran stack 0 ptr = lran stack₀ 0 ptr₀ @ [x₀] ∧ ptr = ptr₀ + 1›
  defines ‹stack[ptr] = x; ptr = ptr + 1›
  by vcg_cs

procedure_spec push (x, stack, ptr) returns (stack, ptr)
  assumes "ptr ≥ 0" ensures ‹stack ` {0..<ptr} = {x<sub>0</sub>} ∪ stack<sub>0</sub> ` {0..<ptr\₀} ∧ ptr = ptr₀ + 1›
  for stack[]
  defines ‹stack[ptr] = x; ptr = ptr + 1›
  by vcg_cs (auto simp: fun_upd_image)

program_spec get_succs
  assumes "j ≤ stop ∧ stop = a (j - 1) ∧ 0 ≤ i"
  ensures "
    stack ` {0..<i} = {x. (j<sub>0</sub> - 1, x) ∈ Edges a ∧ x ∉ set_of visited} ∪ stack<sub>0</sub> ` {0..<i₀}
    ∧ i ≥ i₀"
for i j stop stack[] a[] visited[]
defines
‹
 while (j < stop)
 @invariant ‹stack ` {0..<i} = {x. x ∈ a ` {j₀..<j} ∧ x ∉ set_of visited} ∪ stack₀ ` {0..<i\₀}
 ∧ j ≤ stop ∧ i₀ ≤ i ∧ j₀ ≤ j
 ›
 @variant ‹(stop - j)›
 {
 succ = a[j];
 is_elem = bv_elem(succ,visited);
 if (is_elem == 0) {
 (stack, i) = push (succ, stack, i)
 };
 j = j + 1
 }
 ›
  by vcg_cs (auto simp: intvs_incr_h Edges_def succs_def)

procedure_spec pop (stack, ptr) returns (x, ptr)
  assumes "ptr ≥ 1" ensures ‹stack₀ ` {0..<ptr\<sub>0</sub>} = stack<sub>0</sub> ` {0..<ptr} ∪ {x} ∧ ptr₀ = ptr + 1›
  for stack[]
  defines ‹ptr = ptr - 1; x = stack[ptr]›
  by vcg_cs (simp add: intvs_upper_decr)

procedure_spec stack_init () returns i
  assumes True ensures ‹i = 0›
  defines ‹i = 0›
  by vcg_cs

lemma Edges_empty:
  "Edges a `` {i} = {}" if "i + 1 ≥ a i"
  using that unfolding Edges_def succs_def by auto

text ‹This is one of the main insights of the algorithm: if a set of visited states is
  w.r.t.\ to the edge relation, then it is guaranteed to contain all the states that
  reachable from any state within the set.›
lemma reachability_invariant:
  assumes reachable: "(s, x) ∈ (Edges a)^*"
      and closed: "∀v∈visited. Edges a `` {v} ⊆ visited"
      and start: "s ∈ visited"
  shows "x ∈ visited"
  using reachable start closed by induction auto

program_spec (partial) dfs
  assumes "0 ≤ x ∧ 0 ≤ s"
  ensures "b = 1 ⟷ x ∈ (Edges a)^* `` {s}" defines ‹
 b = 0;
 clear stack[];
 i = stack_init();
 (stack, i) = push (s, stack, i);
 clear visited[];
 while (b == 0 ∧ i ≠ 0)
 @invariant ‹0 ≤ i ∧ (s ∈ stack ` {0..<i} ∨ s ∈ set_of visited) ∧ (b = 0 ∨ b = 1) ∧ (
 if b = 0 then
 stack ` {0..<i} ⊆ (Edges a)<sup>*</sup> `` {s}
 ∧ (∀v ∈ set_of visited. (Edges a) `` {v} ⊆ set_of visited ∪ stack ` {0..<i})
 ∧ (x ∉ set_of visited)
 else x ∈ (Edges a)<sup>*</sup> `` {s})
 ›
 {
 (next, i) = pop(stack, i); ―‹Take the top most element from the stack.›
 visited = bv_insert(next, visited); ―‹Mark it as visited,›
 if (next == x) {
 b = 1 ―‹If it is the target, we are done.›
 } else {
 ―‹Else, put its successors on the stack if they are not yet visited.›
 stop = a[next];
 j = next + 1;
 if (j ≤ stop) {
 inline get_succs
 }
 }
 }
 ›
  apply vcg_cs
  subgoal by (auto simp: set_of_def)
  subgoal using intvs_lower_incr by (auto simp: Edges_empty)
  subgoal by auto (fastforce simp: set_of_def dest!: reachability_invariant)
  done


text ‹Assuming that the input graph is finite, we can also prove that the algorithm terminates.
  will thus use an ‹Isabelle context› to fix a certain finite graph and a start state:
 ›

context
  fixes start :: int and edges
  assumes finite_graph[intro!]: "finite ((Edges edges)^* `` {start})"
begin

lemma sub_insert_same_iff: "s ⊂ insert x s ⟷ x∉s" by auto

program_spec dfs1
  assumes "0 ≤ x ∧ 0 ≤ s ∧ start = s ∧ edges = a"
  ensures "b = 1 ⟷ x ∈ (Edges a)^* `` {s}"
  for visited[]
  defines
‹
 b = 0;
 ―‹‹i› will point to the next free space in the stack (i.e. it is the size of the stack)›
 i = 1;
 ―‹Initially, we put ‹s› on the stack.›
 stack[0] = s;
 visited = bv_init();
 while (b == 0 ∧ i ≠ 0)
 @invariant ‹
 0 ≤ i ∧ (s ∈ stack ` {0..<i} ∨ s ∈ set_of visited) ∧ (b = 0 ∨ b = 1) ∧
 set_of visited ⊆ (Edges edges)<sup>*</sup> `` {start} ∧ (
 if b = 0 then
 stack ` {0..<i} ⊆ (Edges a)<sup>*</sup> `` {s}
 ∧ (∀v ∈ set_of visited. (Edges a) `` {v} ⊆ set_of visited ∪ stack ` {0..<i})
 ∧ (x ∉ set_of visited)
 else x ∈ (Edges a)^* `` {s})
 ›
 @relation ‹finite_psupset ((Edges edges)<sup>*</sup> `` {start}) <*lex*> less_than›
 @variant ‹ (set_of visited, nat i)›
 {
 ―‹Take the top most element from the stack.›
 (next, i) = pop(stack, i);
 if (next == x) {
 ―‹If it is the target, we are done.›
 visited = bv_insert(next, visited);
 b = 1
 } else {
 is_elem = bv_elem(next, visited);
 if (is_elem == 0) {
 visited = bv_insert(next, visited);
 ―‹Else, put its successors on the stack if they are not yet visited.›
 stop = a[next];
 j = next + 1;
 if (j ≤ stop) {
 inline get_succs
 }
 }
 }
 }
 ›
  apply vcg_cs
  subgoal by auto
  subgoal by (auto simp add: image_constant_conv)
  subgoal by (clarsimp simp: finite_psupset_def sub_insert_same_iff)
  subgoal by (auto simp: set_of_def)
  subgoal by (clarsimp simp: finite_psupset_def sub_insert_same_iff)
  subgoal by (auto simp: Edges_empty)
  subgoal by (clarsimp simp: finite_psupset_def sub_insert_same_iff)
  subgoal by (auto simp: set_of_def)
  subgoal by auto (fastforce simp: set_of_def dest!: reachability_invariant)
  done

end


end