Theory Examples

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 "0n"
  ensures "a = 2^nat n0"
  defines ‹
    a = 1;
    c = 0;
    while (c<n) 
      @variant n-c 
      @invariant 0c  cn  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..n0}"
  defines ‹
    s = 0;
    i = 0;
    while (i < n)
      @variant n0 - i
      @invariant n0 = 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 = n0 * n0"
  defines ‹
    a = 0;
    z = 1;
    i = 0;
    while (i < n)
      @variant n0 - i
      @invariant n0 = 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 n0"
  defines ‹
    a = 1;
    i = 1;
    while (i  n)
      @variant n0 + 1 - i
      @invariant n0 = 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 n0 - i
      @invariant n = n0  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 n0"
  defines ‹
    a = 0; b = 1;
    i = 0;
    while (i < n) 
      @variant n0 - i
      @invariant n = n0  (0  i  i  n  n0 < 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 n0 - n›. We locally introduce the name c› for that.
›  
  

program_spec exp_count_down
  assumes "0n"
  ensures "a = 2^nat n0"
  defines ‹
    a = 1;
    while (n>0) 
      @variant n 
      @invariant let c = n0-n in 0n  nn0  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 "0n"  
  ensures "0r  r2  n0  n0 < (r+1)2"
  defines ‹
    r=1;
    while (r*r  n)
      @variant n + 1 - r*r 
      @invariant 0r  (r-1)2  n0    
      { 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 n0 :: int.
       ''invar-final''; 0  n0; ¬ 1 + l < h; 0  l; l < h; l * l  n0;
        n0 < h * h
        n0 < 1 + (l * l + l * 2)"  
  by (smt mult.commute semiring_normalization_rules(3))

program_spec sqr_bisect
  assumes "0n" ensures "r2n0  n0<(r+1)2"
  defines ‹
    l=0; h=n+1;
    while (l+1 < h) 
      @variant h-l
      @invariant 0l  l<h  l2n  n < h2    
    {
      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 a0 b0"
  defines while (ab) 
      @invariant gcd a b = gcd a0 b0  (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 a0 b0"
  defines while (b0) 
      @invariant gcd a b = gcd a0 b0  b0  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 * b0 + s * a0" "q = a div b" "gcd a b = gcd a0 b0"
  shows "gcd b (a - (a0 * (s * q) + b0 * (t * q))) = gcd a0 b0"
proof -
  have "a - (a0 * (s * q) + b0 * (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 * b0 + s * a0" "q = a div b" "b > 0"
  shows "t * (b0 * q) + s * (a0 * q)  a"
proof -
  have "t * (b0 * q) + s * (a0 * 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
as 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 a0 b0  gcd a0 b0 = a0 * old_s + b0 * old_t"
defines ‹
    s = 0;    old_s = 1;
    t = 1;    old_t = 0;
    r = b;    old_r = a;
    while (r0) 
      @invariant gcd old_r r = gcd a0 b0  r0  old_r>0
       a0 * old_s + b0 * old_t = old_r  a0 * s + b0 * 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 * b0 + old_x * a0" "b = y * b0 + x * a0" "q = a div b"
    shows
    "a mod b + (a0 * (x * q) + b0 * (y * q)) = a"
  proof -
    have *: "a0 * (x * q) + b0 * (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 a0 b0  gcd a0 b0 = a0 * x + b0 * y"
defines ‹
    x = 0;
    y = 1;
    old_x = 1;
    old_y = 0;
    while (b0) 
      @invariant gcd a b = gcd a0 b0  b0  a>0  a = a0 * old_x + b0 * old_y  b = a0 * x + b0 * 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 xn works by repeated squaring,
using the following recurrence:
   xn = x * x(n-1)/22 if n› is odd
   xn = xn/22 if n› is even

›
program_spec ex_by_sq
  assumes "n0"
  ensures "r = x0 ^ nat n0"
defines ‹
    r = 1;
    while (n0) 
      @invariant n  0  r * x ^ nat n = x0 ^ nat n0
      @variant n
    {
      if (n mod 2 == 1) {
        r = r * x
      };
      x = x * x;
      n = n / 2
    }
  apply vcg_cs
  subgoal
    apply (subst (asm) ex_by_sq_aux[symmetric])
     apply (smt Suc_1 nat_1 nat_2 nat_mod_distrib nat_one_as_int)
    by (simp add: div_int_pos_iff int_div_less_self nat_div_distrib)
   apply (simp add: algebra_simps semiring_normalization_rules nat_mult_distrib; fail)
  apply (simp add: algebra_simps semiring_normalization_rules nat_mult_distrib; fail)
  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' m0
  defines ‹
    a=1;
    while (m>0) 
      @variant m
      @invariant 0m  mm0  a = tower2' (m0-m)
    {
      n=a;
      
      a = 1;
      while (n>0) 
        @variant n 
        @invariant True ― ‹This will get ugly, there is no n0 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 constexp_count_down program.)
›  
  
program_spec tower2_imp
  assumes m>0
  ensures a = tower2' m0
  defines ‹
    a=1;
    while (m>0) 
      @variant m
      @invariant 0m  mm0  a = tower2' (m0-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⇒int. Negative indexes allowed.
      * Syntactic sugar to use it at index 0 by default
      
    V vname ↦ Vidx vname aexp  Syntax: x[i]
    abbreviation "V x = Vidx x 0"
  
    Assign vname aexp ↦ 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) ⇒ s(x := (s x)(aval i s := aval a s))
    
    (x[] = y, s) ⇒ s(x := s y)
    
  weakest preconditions: Also obvious
    …
    
Reasoning with arrays:    
    
  Usually: Model as int⇒int is appropriate.
  
  Common idioms: 
    * {l..<h} the set of indexes i, l≤i<h
    * a`{l..<h} the set of values in range {l..<h}
  
    * ran_sorted a l h = …
    
    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 "lh"
  ensures "r = (i=l0..<h0. a0 i)"
  defines ‹
    r = 0;
    while (l<h) 
      @invariant l0l  lh  r = (i=l0..<l. a0 i)
      @variant h-l
    {
      r = r + a[l];
      l=l+1
    }
  apply vcg_cs
  by (auto simp: intvs_incdec)
  
(* Exercise: Remove l≤h precondition! *)  

subsubsection ‹Finding Least Index of Element›  
  
program_spec find_least_idx
  assumes lh
  ensures if l=h0 then x0  a0`{l0..<h0} else l{l0..<h0}  a0 l = x0  x0a0`{l0..<l}
  defines while (l<h  a[l]  x) 
      @variant h-l
      @invariant l0l  lh  xa`{l0..<l}
      l=l+1
  ›
  apply vcg_cs
  by (smt atLeastLessThan_iff imageI)
  


  

subsubsection ‹Check for Sortedness›  
  
term ran_sorted

program_spec check_sorted
  assumes lh
  ensures r0  ran_sorted a0 l0 h0
  defines if (l==h) r=1
    else {
      l=l+1;
      while (l<h  a[l-1]  a[l]) 
        @variant h-l
        @invariant l0<l  lh  ran_sorted a l0 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  li  i<h  (j=l..<i. a j) = (j=i..<h. a j)"

program_spec equilibrium
  assumes lh
  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 li  ih  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 li  ih 
         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 = a0 ((i-1) mod n)"
  defines ‹
    i = 0;
    prev = a[n - 1];
    while (i < n)
      @invariant 0  i  i  n 
         (j{0..<i}. a j = a0 ((j-1) mod n))
         (j{i..<n}. a j = a0 j)
         prev = a0 ((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  lh  ran_sorted a l h"

definition "bin_search_post a l h x i  
  li  ih  (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 a0 l0 h0 x0 l
  defines while (l < h) 
      @variant h-l
      @invariant l0l  lh  hh0  (i{l0..<l}. a i < x)  (i{h..<h0}. 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: "xa`{l..<h}  (i<h  x = a i)"
    and bin_search_leftmost: "xa`{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  {l0..<h0}. lran a j (j + (h1 - l1)) = lran b l1 h1)"
  for a[] b[]
  defines ‹
  match = 0;
  while (l < h  match == 0)
    @invariantl0  l  l  h  match  {0,1} 
    (if match = 1
     then lran a l (l + (h1 - l1)) = lran b l1 h1  l < h
     else (j  {l0..<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  {l0..h0-(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  {l0..<h0-(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)
      @invariantl0  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  {l0..<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 l0 h0 = 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 a0 a i
   (let key = a0 j in
      i{l-1..<j}            ― ‹i› is in range›
    ― ‹Content of new array›
     (k{l..i}. a k = a0 k) 
     a (i+1) = key
     (k{i+2..j}. a k = a0 (k-1))
     a = a0 on -{l..j}
    ― ‹Placement of key›  
     (il  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 a0 a :: "int  int"
  
program_spec insort_insert 
  assumes "l<j"
  ensures "insort_insert_post l j a0 a i"
  defines ‹
    key = a[j];
    i = j-1;
    while (il  a[i]>key) 
      @variant i-l+1
      @invariant l-1i  i<j 
         (k{l..i}. a k = a0 k) 
         (k{i+2..j}. a k > key  a k = a0 (k-1))
         a = a0 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 "li+1" "i+1j" 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 li+1 i+1j
    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 li+1 i+1j
    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 a0 {l..<h}  a=a0 on -{l..<h}"
  for a[]
  defines ‹
    j = l + 1;
    while (j<h) 
      @variant h-j
      @invariant l<j  jh           ― ‹j› in range›
         ran_sorted a l j     ― ‹Array is sorted up to j›
         mset_ran a {l..<h} = mset_ran a0 {l..<h} ― ‹Elements in range only permuted›
         a=a0 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 "lh"
  ensures "mset_ran a0 {l0..<h0} = mset_ran a {l0..<h0}
          (j{l0..<i}. a j < p0)
          (j{i..<h0}. a j  p0)
          l0i  ih0
          a0 = a on -{l0..<h0}
         "
  defines ‹
  i=l; j=l;
  while (j<h) 
    @invariant li  ij  jh     
       mset_ran a0 {l0..<h0} = mset_ran a {l0..<h0}
       (k{l..<i}. a k < p)  
       (k{i..<j}. a k  p)  
       (k{j..<h}. a0 k = a k)  
       a0 = a on -{l0..<h0}
    @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 a0 {l0..<h0} = mset_ran a {l0..<h0}
          (j{l0..<i}. a j < a i) 
          (j{i..<h0}. a j  a i)
          l0i  i<h0  a0 (h0-1) = a i
          a0 = a on -{l0..<h0}
         "
  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}. a1 j < a1 i" 
  assumes GEQ: "j{i..<h}. a1 i  a1 j"
  
  assumes R1: "mset_ran a1 {l..<i} = mset_ran a2 {l..<i}"
  assumes E1: "a1 = a2 on - {l..<i}"
  
  assumes SL: "ran_sorted a2 l i"
  
  assumes R2: "mset_ran a2 {i + 1..<h} = mset_ran a3 {i + 1..<h}"
  assumes E2: "a2 = a3 on - {i + 1..<h}"

  assumes SH: "ran_sorted a3 (i + 1) h"
  
  shows "ran_sorted a3 l h"
proof -
  have [simp]: "{l..<i}  - {i + 1..<h}" by auto
  have [simp]: "a1 i = a3 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}. a3 j < a3 i"
    by (simp add: X1 X2)
    
  
  from GEQ have GEQ1: "j{i+1..<h}. a1 i  a1 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}. a3 i  a3 j" by (simp add: X3 X4)

  
  from SL eq_on_xfer_ran_sorted[OF E2, of l i] have SL': "ran_sorted a3 l i" by simp
  
  show ?thesis using combine_sorted_pivot[OF BOUNDS SL' SH LESS' GEQ'] .
qed  
  
lemma quicksort_mset_aux:
  assumes B: "l0i" "i<h0"
  assumes R1: "mset_ran a {l0..<i} = mset_ran aa {l0..<i}"
  assumes E1: "a = aa on - {l0..<i}"
  assumes R2: "mset_ran aa {i + 1..<h0} = mset_ran ab {i + 1..<h0}"
  assumes E2: "aa = ab on - {i + 1..<h0}"
  shows "mset_ran a {l0..<h0} = mset_ran ab {l0..<h0}"
  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 l0 h0  mset_ran a0 {l0..<h0} = mset_ran a {l0..<h0}  a0=a on -{l0..<h0}"
  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 "lh"
  ensures "lran a l0 i = filter (λx. x0) (lran a0 l0 h0)"
  defines ‹
    i=l; j=l;
    while (j<h) 
      @invariant li  ij  jh     
         lran a l i = filter (λx. x0) (lran a0 l j)
         lran a j h = lran a0 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 "l1h1  l2h2"
  ensures "let ms = lran m 0 j; xs0 = lran a10 l10 h10; ys0 = lran a20 l20 h20 in  
    j0  ms = merge xs0 ys0"
  defines ‹
    j=0;
    while (l1h1  l2h2) 
      @variant h1 + h2 - l1 - l2
      @invariant let 
          xs=lran a1 l1 h1; ys = lran a2 l2 h2; ms = lran m 0 j;
          xs0 = lran a10 l10 h10; ys0 = lran a20 l20 h20
        in 
          l1h1  l2h2  0j 
          merge xs0 ys0 = ms@merge xs ys
    {
      if (l2==h2  (l1h1  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 "l1h1  l2h2  sorted (lran a1 l1 h1)  sorted (lran a2 l2 h2)"
  ensures "let ms = lran m 0 j in 
      j0 
     sorted ms 
     mset ms = mset (lran a10 l10 h10) + mset (lran a20 l20 h20)"
  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: "lm; mh  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 "lh"
  ensures 0j  sorted (lran b 0 j)  mset (lran b 0 j) = mset (lran a0 l0 h0)
  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]: "x0  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]: "iset_of a  a i  0" by auto
end

program_spec dedup
  assumes lh
  ensures set (lran a l i) = set (lran a0 l h)  distinct (lran a l i)
  defines ‹
    i=l; j=l;
    clear b[];
    while (j<h) 
      @variant h-j
      @invariant li  ij  jh 
         set (lran a l i) = set (lran a0 l j)
         distinct (lran a l i)
         lran a j h = lran a0 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 x0 (set_of b0)
  defines ‹b[x] = 1›
  by vcg_cs
  
procedure_spec bv_remove (x, b) returns b
  assumes True ensures set_of b = set_of b0 - {x0}
  defines ‹b[x] = 0›
  by vcg_cs
  
procedure_spec bv_elem (x,b) returns r  
  assumes True ensures r0  x0set_of b0
  defines ‹r = b[x]
  by vcg_cs
  
  
procedure_spec dedup' (a,l,h) returns (a,l,i)
  assumes lh ensures set (lran a l i) = set (lran a0 l0 h0)  distinct (lran a l i)
  for b[]
  defines ‹
    b = bv_init();
  
    i=l; j=l;
    
    while (j<h) 
      @variant h-j
      @invariant li  ij  jh 
         set (lran a l i) = set (lran a0 l j)
         distinct (lran a l i)
         lran a j h = lran a0 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 i0 variant i
  defines if (i0) 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 n0  m0 ensures True variant n+m
  defines if (n0  m0) {
      G = G + m;
      if (n<m) rec coast (n-1,m-1) else rec pedal(n-1,m)
    }
and
  coast (n,m) returns () assumes n0  m0 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]: 
  "m0  ack m 0 = ack (m-1) 1"
  "m0; n0  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 "m0  n0" ensures "r=int (ack (nat m0) (nat n0))"
    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<y0 then G = y0-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 "b0  odd a0"
    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 "b0  even a0"
    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 = GZ0 + GX0*GY0 variant "(¦GY¦,1::nat)" 
  for GX GY GZ
  defines
  ‹
    e = even_imp (GY);
    if (e0) rec evenproduct() else rec oddproduct()
and
  oddproduct() assumes odd GY ensures GZ = GZ0 + GX0*GY0 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 = GZ0 + GX0*GY0 variant "(¦GY¦,0::nat)" 
  for GX GY GZ
  defines
  if (GY0) {
      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'' (λx. P x) ≡ ∃x. P x
        GHOST ''name'' value = True
        
        assumes ABS_GHOST name (λx. P x) 
        obtains x where "GHOST name x" "P x"
        
        assumes "GHOST name x" and "P x"
        shows "ABS_GHOST name (λx. P x)"
        
        Let_Ghost name (λx s. C x s) = skip
        
        assumes "∃x. C x s"
        assumes "⋀x. ⟦ GHOST name x; C x s ⟧ ⟹ Q s"
        shows "wp (Let_Ghost name (λ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 ⇒ 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 stack0 0 ptr0 @ [x0]  ptr = ptr0 + 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} = {x0}  stack0 ` {0..<ptr0}  ptr = ptr0 + 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. (j0 - 1, x)  Edges a  x  set_of visited}  stack0 ` {0..<i0}
     i  i0"
for i j stop stack[] a[] visited[]
defines
while (j < stop)
  @invariant stack ` {0..<i} = {x. x  a ` {j0..<j}  x  set_of visited}  stack0 ` {0..<i0}
     j  stop  i0  i  j0  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 stack0 ` {0..<ptr0} = stack0 ` {0..<ptr}  {x}  ptr0 = 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
closed w.r.t.\ to the edge relation, then it is guaranteed to contain all the states that
are reachable from any state within the set.›
lemma reachability_invariant:
  assumes reachable: "(s, x)  (Edges a)*"
      and closed: "vvisited. 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)* `` {s}
       (v  set_of visited. (Edges a) `` {v}  set_of visited  stack ` {0..<i})
       (x  set_of visited)
    else x  (Edges a)* `` {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.
We 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  xs" 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)* `` {start}  (
    if b = 0 then
      stack ` {0..<i}  (Edges a)* `` {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)* `` {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