Theory Common

section "Common"

theory Common
  imports
  "HOL-Library.Going_To_Filter"
  "Akra_Bazzi.Akra_Bazzi_Method"
  "Akra_Bazzi.Akra_Bazzi_Approximation"
  "HOL-Library.Code_Target_Numeral"
  "Root_Balanced_Tree.Time_Monad"
begin

type_synonym point = "int * int"

subsection "Auxiliary Functions and Lemmas"

subsubsection "Time Monad"

lemma time_distrib_bind:
  "time (bind_tm tm f) = time tm + time (f (val tm))"
  unfolding bind_tm_def by (simp split: tm.split)

lemmas time_simps = time_distrib_bind tick_def

lemma bind_tm_cong[fundef_cong]:
  assumes "⋀v. v = val n ⟹ f v = g v" "m = n"
  shows "bind_tm m f = bind_tm n g"
  using assms unfolding bind_tm_def by (auto split: tm.split)

subsubsection "Landau Auxiliary"

text ‹
  The following lemma expresses a procedure for deriving complexity properties of
  the form @{prop"t ∈ O[m going_to at_top within A](f o m)"} where
    ▪ ‹t› is a (timing) function on same data domain (e.g. lists),
    ▪ ‹m› is a measure function on that data domain (e.g. length),
    ▪ ‹t'› is a function on @{typ nat},
    ▪ ‹A› is the set of valid inputs for the data domain.
  One needs to show that
    ▪ ‹t› is bounded by @{term "t' o m"} for valid inputs
    ▪ @{prop"t' ∈ O(f)"}
  to conclude the overall property @{prop"t ∈ O[m going_to at_top within A](f o m)"}.
›

lemma bigo_measure_trans:
  fixes t :: "'a ⇒ real" and t' :: "nat ⇒ real" and m :: "'a ⇒ nat" and f ::"nat ⇒ real"
  assumes "⋀x. x ∈ A ⟹ t x ≤ (t' o m) x"
      and "t' ∈ O(f)"
      and "⋀x. x ∈ A ⟹ 0 ≤ t x"
  shows "t ∈ O[m going_to at_top within A](f o m)"
proof -
  have 0: "⋀x. x ∈ A ⟹ 0 ≤ (t' o m) x" by (meson assms(1,3) order_trans)
  have 1: "t ∈ O[m going_to at_top within A](t' o m)"
    apply(rule bigoI[where c=1]) using assms 0
    by (simp add: eventually_inf_principal going_to_within_def)
  have 2: "t' o m ∈ O[m going_to at_top](f o m)"
    unfolding o_def going_to_def
    by(rule landau_o.big.filtercomap[OF assms(2)])
  have 3: "t' o m ∈ O[m going_to at_top within A](f o m)"
    using landau_o.big.filter_mono[OF _2] going_to_mono[OF _subset_UNIV] by blast
  show ?thesis by(rule landau_o.big_trans[OF 1 3])
qed

lemma const_1_bigo_n_ln_n:
  "(λ(n::nat). 1) ∈ O(λn. n * ln n)"
proof -
  have "∃N. ∀(n::nat) ≥ N. (λx. 1 ≤ x * ln x) n"
  proof -
    have "∀(n::nat) ≥ 3. (λx. 1 ≤ x * ln x) n"
    proof standard
      fix n
      show "3 ≤ n ⟶ 1 ≤ real n * ln (real n)"
      proof standard
        assume "3 ≤ n"
        hence A: "1 ≤ real n"
          by simp
        have B: "1 ≤ ln (real n)" 
          using ln_ln_nonneg' ‹3 ≤ n› by simp
        show "1 ≤ real n * ln (real n)"
          using mult_mono [OF A B] by simp
      qed
    qed
    thus ?thesis
      by blast
  qed
  thus ?thesis
    by auto
qed

subsubsection "Miscellaneous Lemmas"

lemma set_take_drop_i_le_j:
  "i ≤ j ⟹ set xs = set (take j xs) ∪ set (drop i xs)"
proof (induction xs arbitrary: i j)
  case (Cons x xs)
  show ?case
  proof (cases "i = 0")
    case True
    thus ?thesis
      using set_take_subset by force
  next
    case False
    hence "set xs = set (take (j - 1) xs) ∪ set (drop (i - 1) xs)"
      by (simp add: Cons diff_le_mono)
    moreover have "set (take j (x # xs)) = insert x (set (take (j - 1) xs))"
      using False Cons.prems by (auto simp: take_Cons')
    moreover have "set (drop i (x # xs)) = set (drop (i - 1) xs)"
      using False Cons.prems by (auto simp: drop_Cons')
    ultimately show ?thesis
      by auto
  qed
qed simp

lemma set_take_drop:
  "set xs = set (take n xs) ∪ set (drop n xs)"
  using set_take_drop_i_le_j by fast

lemma sorted_wrt_take_drop:
  "sorted_wrt f xs ⟹ ∀x ∈ set (take n xs). ∀y ∈ set (drop n xs). f x y"
  using sorted_wrt_append[of f "take n xs" "drop n xs"] by simp

lemma sorted_wrt_hd_less:
  assumes "sorted_wrt f xs" "⋀x. f x x"
  shows "∀x ∈ set xs. f (hd xs) x"
  using assms by (cases xs) auto

lemma sorted_wrt_hd_less_take:
  assumes "sorted_wrt f (x # xs)" "⋀x. f x x"
  shows "∀y ∈ set (take n (x # xs)). f x y"
  using assms sorted_wrt_hd_less [of f ‹x # xs›] in_set_takeD [of _ n ‹x # xs›]
  by auto

lemma sorted_wrt_take_less_hd_drop:
  assumes "sorted_wrt f xs" "n < length xs"
  shows "∀x ∈ set (take n xs). f x (hd (drop n xs))"
  using sorted_wrt_take_drop assms by fastforce

lemma sorted_wrt_hd_drop_less_drop:
  assumes "sorted_wrt f xs" "⋀x. f x x"
  shows "∀x ∈ set (drop n xs). f (hd (drop n xs)) x"
  using assms sorted_wrt_drop sorted_wrt_hd_less by blast

lemma length_filter_P_impl_Q:
   "(⋀x. P x ⟹ Q x) ⟹ length (filter P xs) ≤ length (filter Q xs)"
  by (induction xs) auto

lemma filter_Un:
  "set xs = A ∪ B ⟹ set (filter P xs) = { x ∈ A. P x } ∪ { x ∈ B. P x }"
  by (induction xs) (auto, metis UnI1 insert_iff, metis UnI2 insert_iff)

subsubsection ‹@{const length}›

fun length_tm :: "'a list ⇒ nat tm" where
  "length_tm [] =1 return 0"
| "length_tm (x # xs) =1
    do {
      l <- length_tm xs;
      return (1 + l)
    }"

lemma length_eq_val_length_tm:
  "val (length_tm xs) = length xs"
  by (induction xs) auto

lemma time_length_tm:
  "time (length_tm xs) = length xs + 1"
  by (induction xs) (auto simp: time_simps)

fun length_it' :: "nat ⇒ 'a list ⇒ nat" where
  "length_it' acc [] = acc"
| "length_it' acc (x#xs) = length_it' (acc+1) xs"

definition length_it :: "'a list ⇒ nat" where
  "length_it xs = length_it' 0 xs"

lemma length_conv_length_it':
  "length xs + acc = length_it' acc xs"
  by (induction acc xs rule: length_it'.induct) auto

lemma length_conv_length_it[code_unfold]:
  "length xs = length_it xs"
  unfolding length_it_def using length_conv_length_it' add_0_right by metis

subsubsection ‹@{const rev}›

fun rev_it' :: "'a list ⇒ 'a list ⇒ 'a list" where
  "rev_it' acc [] = acc"
| "rev_it' acc (x#xs) = rev_it' (x#acc) xs"

definition rev_it :: "'a list ⇒ 'a list" where
  "rev_it xs = rev_it' [] xs"

lemma rev_conv_rev_it':
  "rev xs @ acc = rev_it' acc xs"
  by (induction acc xs rule: rev_it'.induct) auto

lemma rev_conv_rev_it[code_unfold]:
  "rev xs = rev_it xs"
  unfolding rev_it_def using rev_conv_rev_it' append_Nil2 by metis

subsubsection ‹@{const take}›

fun take_tm :: "nat ⇒ 'a list ⇒ 'a list tm" where
  "take_tm n [] =1 return []"
| "take_tm n (x # xs) =1
    (case n of
       0 ⇒ return []
     | Suc m ⇒ do {
         ys <- take_tm m xs;
         return (x # ys)
       }
    )"

lemma take_eq_val_take_tm:
  "val (take_tm n xs) = take n xs"
  by (induction xs arbitrary: n) (auto split: nat.split)

lemma time_take_tm:
  "time (take_tm n xs) = min n (length xs) + 1"
  by (induction xs arbitrary: n) (auto simp: time_simps split: nat.split)

subsubsection ‹@{const filter}›

fun filter_tm :: "('a ⇒ bool) ⇒ 'a list ⇒ 'a list tm" where
  "filter_tm P [] =1 return []"
| "filter_tm P (x # xs) =1
    (if P x then
       do {
         ys <- filter_tm P xs;
         return (x # ys)
       }
     else
       filter_tm P xs
    )"

lemma filter_eq_val_filter_tm:
  "val (filter_tm P xs) = filter P xs"
  by (induction xs) auto

lemma time_filter_tm:
  "time (filter_tm P xs) = length xs + 1"
  by (induction xs) (auto simp: time_simps)

fun filter_it' :: "'a list ⇒ ('a ⇒ bool) ⇒ 'a list ⇒ 'a list" where
  "filter_it' acc P [] = rev acc"
| "filter_it' acc P (x#xs) = (
    if P x then
      filter_it' (x#acc) P xs
    else
      filter_it' acc P xs
  )"

definition filter_it :: "('a ⇒ bool) ⇒ 'a list ⇒ 'a list" where
  "filter_it P xs = filter_it' [] P xs"

lemma filter_conv_filter_it':
  "rev acc @ filter P xs = filter_it' acc P xs"
  by (induction acc P xs rule: filter_it'.induct) auto

lemma filter_conv_filter_it[code_unfold]:
  "filter P xs = filter_it P xs"
  unfolding filter_it_def using filter_conv_filter_it' append_Nil rev.simps(1) by metis

subsubsection ‹‹split_at››

fun split_at_tm :: "nat ⇒ 'a list ⇒ ('a list × 'a list) tm" where
  "split_at_tm n [] =1 return ([], [])"
| "split_at_tm n (x # xs) =1 (
    case n of
      0 ⇒ return ([], x # xs)
    | Suc m ⇒
      do {
        (xs', ys') <- split_at_tm m xs;
        return (x # xs', ys')
      }
  )"

fun split_at :: "nat ⇒ 'a list ⇒ 'a list × 'a list" where
  "split_at n [] = ([], [])"
| "split_at n (x # xs) = (
    case n of 
      0 ⇒ ([], x # xs)
    | Suc m ⇒ 
        let (xs', ys') = split_at m xs in
        (x # xs', ys')
  )"

lemma split_at_eq_val_split_at_tm:
  "val (split_at_tm n xs) = split_at n xs"
  by (induction xs arbitrary: n) (auto split: nat.split prod.split)

lemma split_at_take_drop_conv:
  "split_at n xs = (take n xs, drop n xs)"
  by (induction xs arbitrary: n) (auto simp: split: nat.split)

lemma time_split_at_tm:
  "time (split_at_tm n xs) = min n (length xs) + 1"
  by (induction xs arbitrary: n) (auto simp: time_simps split: nat.split prod.split)

fun split_at_it' :: "'a list ⇒ nat ⇒ 'a list ⇒ ('a list * 'a list)" where
  "split_at_it' acc n [] = (rev acc, [])"
| "split_at_it' acc n (x#xs) = (
    case n of
      0 ⇒ (rev acc, x#xs)
    | Suc m ⇒ split_at_it' (x#acc) m xs
  )"

definition split_at_it :: "nat ⇒ 'a list ⇒ ('a list * 'a list)" where
  "split_at_it n xs = split_at_it' [] n xs"

lemma split_at_conv_split_at_it':
  assumes "(ts, ds) = split_at n xs" "(ts', ds') = split_at_it' acc n xs"
  shows "rev acc @ ts = ts'"
    and "ds = ds'"
  using assms
  by (induction acc n xs arbitrary: ts rule: split_at_it'.induct)
     (auto simp: split: prod.splits nat.splits)

lemma split_at_conv_split_at_it_prod:
  assumes "(ts, ds) = split_at n xs" "(ts', ds') = split_at_it n xs"
  shows "(ts, ds) = (ts', ds')"
  using assms unfolding split_at_it_def 
  using split_at_conv_split_at_it' rev.simps(1) append_Nil by fast+

lemma split_at_conv_split_at_it[code_unfold]:
  "split_at n xs = split_at_it n xs"
  using split_at_conv_split_at_it_prod surj_pair by metis

declare split_at_tm.simps [simp del]
declare split_at.simps [simp del]


subsection "Mergesort"

subsubsection "Functional Correctness Proof"

definition sorted_fst :: "point list ⇒ bool" where
  "sorted_fst ps = sorted_wrt (λp⇩0 p⇩1. fst p⇩0 ≤ fst p⇩1) ps"

definition sorted_snd :: "point list ⇒ bool" where
  "sorted_snd ps = sorted_wrt (λp⇩0 p⇩1. snd p⇩0 ≤ snd p⇩1) ps"

fun merge_tm :: "('b ⇒ 'a::linorder) ⇒ 'b list ⇒ 'b list ⇒ 'b list tm" where
  "merge_tm f (x # xs) (y # ys) =1 (
    if f x ≤ f y then
      do {
        tl <- merge_tm f xs (y # ys);
        return (x # tl)
      }
    else
      do {
        tl <- merge_tm f (x # xs) ys;
        return (y # tl)
      }
  )"
| "merge_tm f [] ys =1 return ys"
| "merge_tm f xs [] =1 return xs"

fun merge :: "('b ⇒ 'a::linorder) ⇒ 'b list ⇒ 'b list ⇒ 'b list" where
  "merge f (x # xs) (y # ys) = (
    if f x ≤ f y then
      x # merge f xs (y # ys)
    else
      y # merge f (x # xs) ys
  )"
| "merge f [] ys = ys"
| "merge f xs [] = xs"

lemma merge_eq_val_merge_tm:
  "val (merge_tm f xs ys) = merge f xs ys"
  by (induction f xs ys rule: merge.induct) auto

lemma length_merge:
  "length (merge f xs ys) = length xs + length ys"
  by (induction f xs ys rule: merge.induct) auto

lemma set_merge:
  "set (merge f xs ys) = set xs ∪ set ys"
  by (induction f xs ys rule: merge.induct) auto

lemma distinct_merge:
  assumes "set xs ∩ set ys = {}" "distinct xs" "distinct ys"
  shows "distinct (merge f xs ys)"
  using assms by (induction f xs ys rule: merge.induct) (auto simp: set_merge)

lemma sorted_merge:
  assumes "P = (λx y. f x ≤ f y)"
  shows "sorted_wrt P (merge f xs ys) ⟷ sorted_wrt P xs ∧ sorted_wrt P ys"
  using assms by (induction f xs ys rule: merge.induct) (auto simp: set_merge)

declare split_at_take_drop_conv [simp]

function (sequential) mergesort_tm :: "('b ⇒ 'a::linorder) ⇒ 'b list ⇒ 'b list tm" where
  "mergesort_tm f [] =1 return []"
| "mergesort_tm f [x] =1 return [x]"
| "mergesort_tm f xs =1 (
    do {
      n <- length_tm xs;
      (xs⇩l, xs⇩r) <- split_at_tm (n div 2) xs;
      l <- mergesort_tm f xs⇩l;
      r <- mergesort_tm f xs⇩r;
      merge_tm f l r
    }
  )"
  by pat_completeness auto
termination mergesort_tm
  by (relation "Wellfounded.measure (λ(_, xs). length xs)")
     (auto simp add: length_eq_val_length_tm split_at_eq_val_split_at_tm)

fun mergesort :: "('b ⇒ 'a::linorder) ⇒ 'b list ⇒ 'b list" where
  "mergesort f [] = []"
| "mergesort f [x] = [x]"
| "mergesort f xs = ( 
    let n = length xs div 2 in
    let (l, r) = split_at n xs in
    merge f (mergesort f l) (mergesort f r)
  )"

declare split_at_take_drop_conv [simp del]

lemma mergesort_eq_val_mergesort_tm:
  "val (mergesort_tm f xs) = mergesort f xs"
  by (induction f xs rule: mergesort.induct)
     (auto simp add: length_eq_val_length_tm split_at_eq_val_split_at_tm merge_eq_val_merge_tm split: prod.split)

lemma sorted_wrt_mergesort:
  "sorted_wrt (λx y. f x ≤ f y) (mergesort f xs)"
  by (induction f xs rule: mergesort.induct) (auto simp: split_at_take_drop_conv sorted_merge)

lemma set_mergesort:
  "set (mergesort f xs) = set xs"
  by (induction f xs rule: mergesort.induct)
     (simp_all add: set_merge split_at_take_drop_conv, metis list.simps(15) set_take_drop)

lemma length_mergesort:
  "length (mergesort f xs) = length xs"
  by (induction f xs rule: mergesort.induct) (auto simp: length_merge split_at_take_drop_conv)

lemma distinct_mergesort:
  "distinct xs ⟹ distinct (mergesort f xs)"
proof (induction f xs rule: mergesort.induct)
  case (3 f x y xs)
  let ?xs' = "x # y # xs"
  obtain l r where lr_def: "(l, r) = split_at (length ?xs' div 2) ?xs'"
    by (metis surj_pair)
  have "distinct l" "distinct r"
    using "3.prems" split_at_take_drop_conv distinct_take distinct_drop lr_def by (metis prod.sel)+
  hence "distinct (mergesort f l)" "distinct (mergesort f r)"
    using "3.IH" lr_def by auto
  moreover have "set l ∩ set r = {}"
    using "3.prems" split_at_take_drop_conv lr_def by (metis append_take_drop_id distinct_append prod.sel)
  ultimately show ?case
    using lr_def by (auto simp: distinct_merge set_mergesort split: prod.splits)
qed auto

lemmas mergesort = sorted_wrt_mergesort set_mergesort length_mergesort distinct_mergesort

lemma sorted_fst_take_less_hd_drop:
  assumes "sorted_fst ps" "n < length ps"
  shows "∀p ∈ set (take n ps). fst p ≤ fst (hd (drop n ps))"
  using assms sorted_wrt_take_less_hd_drop[of "λp⇩0 p⇩1. fst p⇩0 ≤ fst p⇩1"] sorted_fst_def by fastforce

lemma sorted_fst_hd_drop_less_drop:
  assumes "sorted_fst ps"
  shows "∀p ∈ set (drop n ps). fst (hd (drop n ps)) ≤ fst p"
  using assms sorted_wrt_hd_drop_less_drop[of "λp⇩0 p⇩1. fst p⇩0 ≤ fst p⇩1"] sorted_fst_def by fastforce

subsubsection "Time Complexity Proof"

lemma time_merge_tm:
  "time (merge_tm f xs ys) ≤ length xs + length ys + 1"
  by (induction f xs ys rule: merge_tm.induct) (auto simp: time_simps)

function mergesort_recurrence :: "nat ⇒ real" where
  "mergesort_recurrence 0 = 1"
| "mergesort_recurrence 1 = 1"
| "2 ≤ n ⟹ mergesort_recurrence n = 4 + 3 * n + mergesort_recurrence (nat ⌊real n / 2⌋) + 
    mergesort_recurrence (nat ⌈real n / 2⌉)"
  by force simp_all
termination by akra_bazzi_termination simp_all

lemma mergesort_recurrence_nonneg[simp]:
  "0 ≤ mergesort_recurrence n"
  by (induction n rule: mergesort_recurrence.induct) (auto simp del: One_nat_def)

lemma time_mergesort_conv_mergesort_recurrence:
  "time (mergesort_tm f xs) ≤ mergesort_recurrence (length xs)"
proof (induction f xs rule: mergesort_tm.induct)
  case (1 f)
  thus ?case by (auto simp: time_simps)
next
  case (2 f x)
  thus ?case using mergesort_recurrence.simps(2) by (auto simp: time_simps)
next
  case (3 f x y xs')

  define xs where "xs = x # y # xs'"
  define n where "n = length xs"
  obtain l r where lr_def: "(l, r) = split_at (n div 2) xs"
    using prod.collapse by blast
  define l' where "l' = mergesort f l"
  define r' where "r' = mergesort f r"
  note defs = xs_def n_def lr_def l'_def r'_def

  have IHL: "time (mergesort_tm f l) ≤ mergesort_recurrence (length l)"
    using defs "3.IH"(1) by (auto simp: length_eq_val_length_tm split_at_eq_val_split_at_tm)
  have IHR: "time (mergesort_tm f r) ≤ mergesort_recurrence (length r)"
    using defs "3.IH"(2) by (auto simp: length_eq_val_length_tm split_at_eq_val_split_at_tm)

  have *: "length l = n div 2" "length r = n - n div 2"
    using defs by (auto simp: split_at_take_drop_conv)
  hence "(nat ⌊real n / 2⌋) = length l" "(nat ⌈real n / 2⌉) = length r"
    by linarith+
  hence IH: "time (mergesort_tm f l) ≤ mergesort_recurrence (nat ⌊real n / 2⌋)"
            "time (mergesort_tm f r) ≤ mergesort_recurrence (nat ⌈real n / 2⌉)"
    using IHL IHR by simp_all

  have "n = length l + length r"
    using * by linarith
  hence "time (merge_tm f l' r') ≤ n + 1"
    using time_merge_tm defs by (metis length_mergesort)
  
  have "time (mergesort_tm f xs) = 1 + time (length_tm xs) + time (split_at_tm (n div 2) xs) + 
          time (mergesort_tm f l) + time (mergesort_tm f r) + time (merge_tm f l' r')"
    using defs by (auto simp add: time_simps length_eq_val_length_tm  mergesort_eq_val_mergesort_tm 
                                  split_at_eq_val_split_at_tm 
                        split: prod.split)
  also have "... ≤ 4 + 3 * n + time (mergesort_tm f l) + time (mergesort_tm f r)"
    using time_length_tm[of xs] time_split_at_tm[of "n div 2" xs] n_def ‹time (merge_tm f l' r') ≤ n + 1› by simp
  also have "... ≤ 4 + 3 * n + mergesort_recurrence (nat ⌊real n / 2⌋) + mergesort_recurrence (nat ⌈real n / 2⌉)"
    using IH by simp
  also have "... = mergesort_recurrence n"
    using defs by simp
  finally show ?case
    using defs by simp
qed

theorem mergesort_recurrence:
  "mergesort_recurrence ∈ Θ(λn. n * ln n)"
  by (master_theorem) auto

theorem time_mergesort_tm_bigo:
  "(λxs. time (mergesort_tm f xs)) ∈ O[length going_to at_top]((λn. n * ln n) o length)"
proof -
  have 0: "⋀xs. time (mergesort_tm f xs) ≤ (mergesort_recurrence o length) xs"
    unfolding comp_def using time_mergesort_conv_mergesort_recurrence by blast
  show ?thesis
    using bigo_measure_trans[OF 0] by (simp add: bigthetaD1 mergesort_recurrence)
qed

subsubsection "Code Export"

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

fun merge_it' :: "('b ⇒ 'a::linorder) ⇒ 'b list ⇒ 'b list ⇒ 'b list ⇒ 'b list" where
  "merge_it' f acc [] [] = rev acc"
| "merge_it' f acc (x#xs) [] = merge_it' f (x#acc) xs []"
| "merge_it' f acc [] (y#ys) = merge_it' f (y#acc) ys []"
| "merge_it' f acc (x#xs) (y#ys) = (
    if f x ≤ f y then
      merge_it' f (x#acc) xs (y#ys)
    else
      merge_it' f (y#acc) (x#xs) ys
  )"

definition merge_it :: "('b ⇒ 'a::linorder) ⇒ 'b list ⇒ 'b list ⇒ 'b list" where
  "merge_it f xs ys = merge_it' f [] xs ys"

lemma merge_conv_merge_it':
  "rev acc @ merge f xs ys = merge_it' f acc xs ys"
  by (induction f acc xs ys rule: merge_it'.induct) auto

lemma merge_conv_merge_it[code_unfold]:
  "merge f xs ys = merge_it f xs ys"
  unfolding merge_it_def using merge_conv_merge_it' rev.simps(1) append_Nil by metis


subsection "Minimal Distance"

definition sparse :: "real ⇒ point set ⇒ bool" where
  "sparse δ ps ⟷ (∀p⇩0 ∈ ps. ∀p⇩1 ∈ ps. p⇩0 ≠ p⇩1 ⟶ δ ≤ dist p⇩0 p⇩1)"

lemma sparse_identity:
  assumes "sparse δ (set ps)" "∀p ∈ set ps. δ ≤ dist p⇩0 p"
  shows "sparse δ (set (p⇩0 # ps))"
  using assms by (simp add: dist_commute sparse_def)

lemma sparse_update:
  assumes "sparse δ (set ps)"
  assumes "dist p⇩0 p⇩1 ≤ δ" "∀p ∈ set ps. dist p⇩0 p⇩1 ≤ dist p⇩0 p"
  shows "sparse (dist p⇩0 p⇩1) (set (p⇩0 # ps))"
  using assms by (auto simp: dist_commute sparse_def, force+)

lemma sparse_mono:
  "sparse Δ P ⟹ δ ≤ Δ ⟹ sparse δ P"
  unfolding sparse_def by fastforce


subsection "Distance"

lemma dist_transform:
  fixes p :: point and δ :: real and l :: int
  shows "dist p (l, snd p) < δ ⟷ l - δ < fst p ∧ fst p < l + δ"
proof -
  have "dist p (l, snd p) = sqrt ((real_of_int (fst p) - l)⇧2)"
    by (auto simp add: dist_prod_def dist_real_def prod.case_eq_if)
  thus ?thesis
    by auto
qed

fun dist_code :: "point ⇒ point ⇒ int" where
  "dist_code p⇩0 p⇩1 = (fst p⇩0 - fst p⇩1)⇧2 + (snd p⇩0 - snd p⇩1)⇧2"

lemma dist_eq_sqrt_dist_code:
  fixes p⇩0 :: point
  shows "dist p⇩0 p⇩1 = sqrt (dist_code p⇩0 p⇩1)"
  by (auto simp: dist_prod_def dist_real_def split: prod.splits) 

lemma dist_eq_dist_code_lt:
  fixes p⇩0 :: point
  shows "dist p⇩0 p⇩1 < dist p⇩2 p⇩3 ⟷ dist_code p⇩0 p⇩1 < dist_code p⇩2 p⇩3"
  using dist_eq_sqrt_dist_code real_sqrt_less_iff by presburger

lemma dist_eq_dist_code_le:
  fixes p⇩0 :: point
  shows "dist p⇩0 p⇩1 ≤ dist p⇩2 p⇩3 ⟷ dist_code p⇩0 p⇩1 ≤ dist_code p⇩2 p⇩3"
  using dist_eq_sqrt_dist_code real_sqrt_le_iff by presburger

lemma dist_eq_dist_code_abs_lt:
  fixes p⇩0 :: point
  shows "¦c¦ < dist p⇩0 p⇩1 ⟷ c⇧2 < dist_code p⇩0 p⇩1"
  using dist_eq_sqrt_dist_code
  by (metis of_int_less_of_int_power_cancel_iff real_sqrt_abs real_sqrt_less_iff)

lemma dist_eq_dist_code_abs_le:
  fixes p⇩0 :: point
  shows "dist p⇩0 p⇩1 ≤ ¦c¦ ⟷ dist_code p⇩0 p⇩1 ≤ c⇧2"
  using dist_eq_sqrt_dist_code
  by (metis of_int_power_le_of_int_cancel_iff real_sqrt_abs real_sqrt_le_iff)

lemma dist_fst_abs:
  fixes p :: point and l :: int
  shows "dist p (l, snd p) = ¦fst p - l¦"
proof -
  have "dist p (l, snd p) = sqrt ((real_of_int (fst p) - l)⇧2)"
    by (simp add: dist_prod_def dist_real_def prod.case_eq_if)
  thus ?thesis
    by simp
qed

declare dist_code.simps [simp del]


subsection "Brute Force Closest Pair Algorithm"

subsubsection "Functional Correctness Proof"

fun find_closest_bf_tm :: "point ⇒ point list ⇒ point tm" where
  "find_closest_bf_tm _ [] =1 return undefined"
| "find_closest_bf_tm _ [p] =1 return p"
| "find_closest_bf_tm p (p⇩0 # ps) =1 (
    do {
      p⇩1 <- find_closest_bf_tm p ps;
      if dist p p⇩0 < dist p p⇩1 then
        return p⇩0
      else
        return p⇩1
    }
  )"

fun find_closest_bf :: "point ⇒ point list ⇒ point" where
  "find_closest_bf _ [] = undefined"
| "find_closest_bf _ [p] = p"
| "find_closest_bf p (p⇩0 # ps) = (
    let p⇩1 = find_closest_bf p ps in
    if dist p p⇩0 < dist p p⇩1 then
      p⇩0
    else
      p⇩1
  )"

lemma find_closest_bf_eq_val_find_closest_bf_tm:
  "val (find_closest_bf_tm p ps) = find_closest_bf p ps"
  by (induction p ps rule: find_closest_bf.induct) (auto simp: Let_def)

lemma find_closest_bf_set:
  "0 < length ps ⟹ find_closest_bf p ps ∈ set ps"
  by (induction p ps rule: find_closest_bf.induct)
     (auto simp: Let_def split: prod.splits if_splits)

lemma find_closest_bf_dist:
  "∀q ∈ set ps. dist p (find_closest_bf p ps) ≤ dist p q"
  by (induction p ps rule: find_closest_bf.induct)
     (auto split: prod.splits)

fun closest_pair_bf_tm :: "point list ⇒ (point × point) tm" where
  "closest_pair_bf_tm [] =1 return undefined"
| "closest_pair_bf_tm [_] =1 return undefined"
| "closest_pair_bf_tm [p⇩0, p⇩1] =1 return (p⇩0, p⇩1)"
| "closest_pair_bf_tm (p⇩0 # ps) =1 (
    do {
      (c⇩0::point, c⇩1::point) <- closest_pair_bf_tm ps;
      p⇩1 <- find_closest_bf_tm p⇩0 ps;
      if dist c⇩0 c⇩1 ≤ dist p⇩0 p⇩1 then
        return (c⇩0, c⇩1)
      else
        return (p⇩0, p⇩1)
    }
  )"

fun closest_pair_bf :: "point list ⇒ (point * point)" where
  "closest_pair_bf [] = undefined"
| "closest_pair_bf [_] = undefined"
| "closest_pair_bf [p⇩0, p⇩1] = (p⇩0, p⇩1)"
| "closest_pair_bf (p⇩0 # ps) = (
    let (c⇩0, c⇩1) = closest_pair_bf ps in
    let p⇩1 = find_closest_bf p⇩0 ps in
    if dist c⇩0 c⇩1 ≤ dist p⇩0 p⇩1 then
      (c⇩0, c⇩1)
    else
      (p⇩0, p⇩1) 
  )"

lemma closest_pair_bf_eq_val_closest_pair_bf_tm:
  "val (closest_pair_bf_tm ps) = closest_pair_bf ps"
  by (induction ps rule: closest_pair_bf.induct) 
     (auto simp: Let_def find_closest_bf_eq_val_find_closest_bf_tm split: prod.split)

lemma closest_pair_bf_c0:
  "1 < length ps ⟹ (c⇩0, c⇩1) = closest_pair_bf ps ⟹ c⇩0 ∈ set ps"
  by (induction ps arbitrary: c⇩0 c⇩1 rule: closest_pair_bf.induct)
     (auto simp: Let_def find_closest_bf_set split: if_splits prod.splits)

lemma closest_pair_bf_c1:
  "1 < length ps ⟹ (c⇩0, c⇩1) = closest_pair_bf ps ⟹ c⇩1 ∈ set ps" 
proof (induction ps arbitrary: c⇩0 c⇩1 rule: closest_pair_bf.induct)
  case (4 p⇩0 p⇩2 p⇩3 ps)
  let ?ps = "p⇩2 # p⇩3 # ps"
  obtain c⇩0 c⇩1 where c⇩0_def: "(c⇩0, c⇩1) = closest_pair_bf ?ps"
    using prod.collapse by blast
  define p⇩1 where p⇩1_def: "p⇩1 = find_closest_bf p⇩0 ?ps"
  note defs = c⇩0_def p⇩1_def
  have "c⇩1 ∈ set ?ps"
    using "4.IH" defs by simp
  moreover have "p⇩1 ∈ set ?ps"
    using find_closest_bf_set defs by blast
  ultimately show ?case
    using "4.prems"(2) defs by (auto simp: Let_def split: prod.splits if_splits)
qed auto

lemma closest_pair_bf_c0_ne_c1:
  "1 < length ps ⟹ distinct ps ⟹ (c⇩0, c⇩1) = closest_pair_bf ps ⟹ c⇩0 ≠ c⇩1"
proof (induction ps arbitrary: c⇩0 c⇩1 rule: closest_pair_bf.induct)
  case (4 p⇩0 p⇩2 p⇩3 ps)
  let ?ps = "p⇩2 # p⇩3 # ps"
  obtain c⇩0 c⇩1 where c⇩0_def: "(c⇩0, c⇩1) = closest_pair_bf ?ps"
    using prod.collapse by blast
  define p⇩1 where p⇩1_def: "p⇩1 = find_closest_bf p⇩0 ?ps"
  note defs = c⇩0_def p⇩1_def
  have "c⇩0 ≠ c⇩1"
    using "4.IH" "4.prems"(2) defs by simp
  moreover have "p⇩0 ≠ p⇩1"
    using find_closest_bf_set "4.prems"(2) defs
    by (metis distinct.simps(2) length_pos_if_in_set list.set_intros(1))
  ultimately show ?case
    using "4.prems"(3) defs by (auto simp: Let_def split: prod.splits if_splits)
qed auto

lemmas closest_pair_bf_c0_c1 = closest_pair_bf_c0 closest_pair_bf_c1 closest_pair_bf_c0_ne_c1

lemma closest_pair_bf_dist:
  assumes "1 < length ps" "(c⇩0, c⇩1) = closest_pair_bf ps"
  shows "sparse (dist c⇩0 c⇩1) (set ps)"
  using assms
proof (induction ps arbitrary: c⇩0 c⇩1 rule: closest_pair_bf.induct)
  case (4 p⇩0 p⇩2 p⇩3 ps)
  let ?ps = "p⇩2 # p⇩3 # ps"
  obtain c⇩0 c⇩1 where c⇩0_def: "(c⇩0, c⇩1) = closest_pair_bf ?ps"
    using prod.collapse by blast
  define p⇩1 where p⇩1_def: "p⇩1 = find_closest_bf p⇩0 ?ps"
  note defs = c⇩0_def p⇩1_def
  hence IH: "sparse (dist c⇩0 c⇩1) (set ?ps)"
    using 4 c⇩0_def by simp
  have *: "∀p ∈ set ?ps. (dist p⇩0 p⇩1) ≤ dist p⇩0 p"
    using find_closest_bf_dist defs by blast
  show ?case
  proof (cases "dist c⇩0 c⇩1 ≤ dist p⇩0 p⇩1")
    case True
    hence "∀p ∈ set ?ps. dist c⇩0 c⇩1 ≤ dist p⇩0 p"
      using * by auto
    hence "sparse (dist c⇩0 c⇩1) (set (p⇩0 # ?ps))"
      using sparse_identity IH by blast
    thus ?thesis
      using True "4.prems" defs by (auto split: prod.splits)
  next
    case False
    hence "sparse (dist p⇩0 p⇩1) (set (p⇩0 # ?ps))"
      using sparse_update[of "dist c⇩0 c⇩1" ?ps p⇩0 p⇩1] IH * defs by argo
    thus ?thesis
      using False "4.prems" defs by (auto split: prod.splits)
  qed
qed (auto simp: dist_commute sparse_def)

subsubsection "Time Complexity Proof"

lemma time_find_closest_bf_tm:
  "time (find_closest_bf_tm p ps) ≤ length ps + 1"
  by (induction p ps rule: find_closest_bf_tm.induct) (auto simp: time_simps)

lemma time_closest_pair_bf_tm:
  "time (closest_pair_bf_tm ps) ≤ length ps * length ps + 1"
proof (induction ps rule: closest_pair_bf_tm.induct)
  case (4 p⇩0 p⇩2 p⇩3 ps)
  let ?ps = "p⇩2 # p⇩3 # ps"
  have "time (closest_pair_bf_tm (p⇩0 # ?ps)) = 1 + time (find_closest_bf_tm p⇩0 ?ps) + time (closest_pair_bf_tm ?ps)"
    by (auto simp: time_simps split: prod.split)
  also have "... ≤ 2 + length ?ps + time (closest_pair_bf_tm ?ps)"
    using time_find_closest_bf_tm[of p⇩0 ?ps] by simp
  also have "... ≤ 2 + length ?ps + length ?ps * length ?ps + 1"
    using "4.IH" by simp
  also have "... ≤ length (p⇩0 # ?ps) * length (p⇩0 # ?ps) + 1"
    by auto
  finally show ?case
    by blast
qed (auto simp: time_simps)

subsubsection "Code Export"

fun find_closest_bf_code :: "point ⇒ point list ⇒ (int * point)" where
  "find_closest_bf_code p [] = undefined"
| "find_closest_bf_code p [p⇩0] = (dist_code p p⇩0, p⇩0)"
| "find_closest_bf_code p (p⇩0 # ps) = (
    let (δ⇩1, p⇩1) = find_closest_bf_code p ps in
    let δ⇩0 = dist_code p p⇩0 in
    if δ⇩0 < δ⇩1 then
      (δ⇩0, p⇩0)
    else
      (δ⇩1, p⇩1)
  )"

lemma find_closest_bf_code_dist_eq:
  "0 < length ps ⟹ (δ, c) = find_closest_bf_code p ps ⟹ δ = dist_code p c"
  by (induction p ps rule: find_closest_bf_code.induct)
     (auto simp: Let_def split: prod.splits if_splits)

lemma find_closest_bf_code_eq:
  "0 < length ps ⟹ c = find_closest_bf p ps ⟹ (δ', c') = find_closest_bf_code p ps ⟹ c = c'"
proof (induction p ps arbitrary: c δ' c' rule: find_closest_bf.induct)
  case (3 p p⇩0 p⇩2 ps)
  define δ⇩0 δ⇩0' where δ⇩0_def: "δ⇩0 = dist p p⇩0" "δ⇩0' = dist_code p p⇩0"
  obtain δ⇩1 p⇩1 δ⇩1' p⇩1' where δ⇩1_def: "δ⇩1 = dist p p⇩1" "p⇩1 = find_closest_bf p (p⇩2 # ps)"
    "(δ⇩1', p⇩1') = find_closest_bf_code p (p⇩2 # ps)"
    using prod.collapse by blast+
  note defs = δ⇩0_def δ⇩1_def
  have *: "p⇩1 = p⇩1'"
    using "3.IH" defs by simp
  hence "δ⇩0 < δ⇩1 ⟷ δ⇩0' < δ⇩1'"
    using find_closest_bf_code_dist_eq[of "p⇩2 # ps" δ⇩1' p⇩1' p]
          dist_eq_dist_code_lt defs
    by simp
  thus ?case
    using "3.prems"(2,3) * defs by (auto split: prod.splits)
qed auto

declare find_closest_bf_code.simps [simp del]

fun closest_pair_bf_code :: "point list ⇒ (int * point * point)" where
  "closest_pair_bf_code [] = undefined"
| "closest_pair_bf_code [p⇩0] = undefined"
| "closest_pair_bf_code [p⇩0, p⇩1] = (dist_code p⇩0 p⇩1, p⇩0, p⇩1)"
| "closest_pair_bf_code (p⇩0 # ps) = (
    let (δ⇩c, c⇩0, c⇩1) = closest_pair_bf_code ps in
    let (δ⇩p, p⇩1) = find_closest_bf_code p⇩0 ps in
    if δ⇩c ≤ δ⇩p then
      (δ⇩c, c⇩0, c⇩1)
    else
      (δ⇩p, p⇩0, p⇩1) 
  )"

lemma closest_pair_bf_code_dist_eq:
  "1 < length ps ⟹ (δ, c⇩0, c⇩1) = closest_pair_bf_code ps ⟹ δ = dist_code c⇩0 c⇩1"
proof (induction ps arbitrary: δ c⇩0 c⇩1 rule: closest_pair_bf_code.induct)
  case (4 p⇩0 p⇩2 p⇩3 ps)
  let ?ps = "p⇩2 # p⇩3 # ps"
  obtain δ⇩c c⇩0 c⇩1 where δ⇩c_def: "(δ⇩c, c⇩0, c⇩1) = closest_pair_bf_code ?ps"
    by (metis prod_cases3)
  obtain δ⇩p p⇩1 where δ⇩p_def: "(δ⇩p, p⇩1) = find_closest_bf_code p⇩0 ?ps"
    using prod.collapse by blast
  note defs = δ⇩c_def δ⇩p_def
  have "δ⇩c = dist_code c⇩0 c⇩1"
    using "4.IH" defs by simp
  moreover have "δ⇩p = dist_code p⇩0 p⇩1"
    using find_closest_bf_code_dist_eq defs by blast
  ultimately show ?case
    using "4.prems"(2) defs by (auto split: prod.splits if_splits)
qed auto

lemma closest_pair_bf_code_eq:
  assumes "1 < length ps" 
  assumes "(c⇩0, c⇩1) = closest_pair_bf ps" "(δ', c⇩0', c⇩1') = closest_pair_bf_code ps"
  shows "c⇩0 = c⇩0' ∧ c⇩1 = c⇩1'"
  using assms
proof (induction ps arbitrary: c⇩0 c⇩1 δ' c⇩0' c⇩1' rule: closest_pair_bf_code.induct)
  case (4 p⇩0 p⇩2 p⇩3 ps)
  let ?ps = "p⇩2 # p⇩3 # ps"
  obtain c⇩0 c⇩1 δ⇩c' c⇩0' c⇩1' where δ⇩c_def: "(c⇩0, c⇩1) = closest_pair_bf ?ps"
    "(δ⇩c', c⇩0', c⇩1') = closest_pair_bf_code ?ps"
    by (metis prod_cases3)
  obtain p⇩1 δ⇩p' p⇩1' where δ⇩p_def: "p⇩1 = find_closest_bf p⇩0 ?ps"
    "(δ⇩p', p⇩1') = find_closest_bf_code p⇩0 ?ps"
    using prod.collapse by blast
  note defs = δ⇩c_def δ⇩p_def
  have A: "c⇩0 = c⇩0' ∧ c⇩1 = c⇩1'"
    using "4.IH" defs by simp
  moreover have B: "p⇩1 = p⇩1'"
    using find_closest_bf_code_eq defs by blast
  moreover have "δ⇩c' = dist_code c⇩0' c⇩1'"
    using defs closest_pair_bf_code_dist_eq[of ?ps] by simp
  moreover have "δ⇩p' = dist_code p⇩0 p⇩1'"
    using defs find_closest_bf_code_dist_eq by blast
  ultimately have "dist c⇩0 c⇩1 ≤ dist p⇩0 p⇩1 ⟷ δ⇩c' ≤ δ⇩p'"
    by (simp add: dist_eq_dist_code_le)
  thus ?case
    using "4.prems"(2,3) defs A B by (auto simp: Let_def split: prod.splits)
qed auto


subsection "Geometry"

subsubsection "Band Filter"

lemma set_band_filter_aux:
  fixes δ :: real and ps :: "point list"
  assumes "p⇩0 ∈ ps⇩L" "p⇩1 ∈ ps⇩R" "p⇩0 ≠ p⇩1" "dist p⇩0 p⇩1 < δ" "set ps = ps⇩L ∪ ps⇩R"
  assumes "∀p ∈ ps⇩L. fst p ≤ l" "∀p ∈ ps⇩R. l ≤ fst p"
  assumes "ps' = filter (λp. l - δ < fst p ∧ fst p < l + δ) ps"
  shows "p⇩0 ∈ set ps' ∧ p⇩1 ∈ set ps'"
proof (rule ccontr)
  assume "¬ (p⇩0 ∈ set ps' ∧ p⇩1 ∈ set ps')"
  then consider (A) "p⇩0 ∉ set ps' ∧ p⇩1 ∉ set ps'"
              | (B) "p⇩0 ∈ set ps' ∧ p⇩1 ∉ set ps'"
              | (C) "p⇩0 ∉ set ps' ∧ p⇩1 ∈ set ps'"
    by blast
  thus False
  proof cases
    case A
    hence "fst p⇩0 ≤ l - δ ∨ l + δ ≤ fst p⇩0" "fst p⇩1 ≤ l - δ ∨ l + δ ≤ fst p⇩1"
      using assms(1,2,5,8) by auto
    hence "fst p⇩0 ≤ l - δ" "l + δ ≤ fst p⇩1"
      using assms(1,2,6,7) by force+
    hence "δ ≤ dist (fst p⇩0) (fst p⇩1)"
      using dist_real_def by simp
    hence "δ ≤ dist p⇩0 p⇩1"
      using dist_fst_le[of p⇩0 p⇩1] by (auto split: prod.splits)
    then show ?thesis
      using assms(4) by fastforce
  next
    case B
    hence "fst p⇩1 ≤ l - δ ∨ l + δ ≤ fst p⇩1"
      using assms(2,5,8) by auto
    hence "l + δ ≤ fst p⇩1"
      using assms(2,7) by auto
    moreover have "fst p⇩0 ≤ l"
      using assms(1,6) by simp
    ultimately have "δ ≤ dist (fst p⇩0) (fst p⇩1)"
      using dist_real_def by simp
    hence "δ ≤ dist p⇩0 p⇩1"
      using dist_fst_le[of p⇩0 p⇩1] less_le_trans by (auto split: prod.splits)
    thus ?thesis
      using assms(4) by simp
  next
    case C
    hence "fst p⇩0 ≤ l - δ ∨ l + δ ≤ fst p⇩0"
      using assms(1,2,5,8) by auto
    hence "fst p⇩0 ≤ l - δ"
      using assms(1,6) by auto
    moreover have "l ≤ fst p⇩1"
      using assms(2,7) by simp
    ultimately have "δ ≤ dist (fst p⇩0) (fst p⇩1)"
      using dist_real_def by simp
    hence "δ ≤ dist p⇩0 p⇩1"
      using dist_fst_le[of p⇩0 p⇩1] less_le_trans by (auto split: prod.splits)
    thus ?thesis
      using assms(4) by simp
  qed
qed
  
lemma set_band_filter:
  fixes δ :: real and ps :: "point list"
  assumes "p⇩0 ∈ set ps" "p⇩1 ∈ set ps" "p⇩0 ≠ p⇩1" "dist p⇩0 p⇩1 < δ" "set ps = ps⇩L ∪ ps⇩R"
  assumes "sparse δ ps⇩L" "sparse δ ps⇩R"
  assumes "∀p ∈ ps⇩L. fst p ≤ l" "∀p ∈ ps⇩R. l ≤ fst p"
  assumes "ps' = filter (λp. l - δ < fst p ∧ fst p < l + δ) ps"
  shows "p⇩0 ∈ set ps' ∧ p⇩1 ∈ set ps'"
proof -
  have "p⇩0 ∉ ps⇩L ∨ p⇩1 ∉ ps⇩L" "p⇩0 ∉ ps⇩R ∨ p⇩1 ∉ ps⇩R"
    using assms(3,4,6,7) sparse_def by force+
  then consider (A) "p⇩0 ∈ ps⇩L ∧ p⇩1 ∈ ps⇩R" | (B) "p⇩0 ∈ ps⇩R ∧ p⇩1 ∈ ps⇩L"
    using assms(1,2,5) by auto
  thus ?thesis
  proof cases
    case A
    thus ?thesis
      using set_band_filter_aux assms(3,4,5,8,9,10) by (auto split: prod.splits)
  next
    case B
    moreover have "dist p⇩1 p⇩0 < δ"
      using assms(4) dist_commute by metis
    ultimately show ?thesis
      using set_band_filter_aux assms(3)[symmetric] assms(5,8,9,10) by (auto split: prod.splits)
  qed
qed

subsubsection "2D-Boxes and Points"

lemma cbox_2D: 
  fixes x⇩0 :: real and y⇩0 :: real
  shows "cbox (x⇩0, y⇩0) (x⇩1, y⇩1) = { (x, y). x⇩0 ≤ x ∧ x ≤ x⇩1 ∧ y⇩0 ≤ y ∧ y ≤ y⇩1 }"
  by fastforce

lemma mem_cbox_2D:
  fixes x :: real and y :: real
  shows "x⇩0 ≤ x ∧ x ≤ x⇩1 ∧ y⇩0 ≤ y ∧ y ≤ y⇩1 ⟷ (x, y) ∈ cbox (x⇩0, y⇩0) (x⇩1, y⇩1)"
  by fastforce

lemma cbox_top_un:
  fixes x⇩0 :: real and y⇩0 :: real
  assumes "y⇩0 ≤ y⇩1" "y⇩1 ≤ y⇩2"
  shows "cbox (x⇩0, y⇩0) (x⇩1, y⇩1) ∪ cbox (x⇩0, y⇩1) (x⇩1, y⇩2) = cbox (x⇩0, y⇩0) (x⇩1, y⇩2)"
  using assms by auto

lemma cbox_right_un:
  fixes x⇩0 :: real and y⇩0 :: real
  assumes "x⇩0 ≤ x⇩1" "x⇩1 ≤ x⇩2"
  shows "cbox (x⇩0, y⇩0) (x⇩1, y⇩1) ∪ cbox (x⇩1, y⇩0) (x⇩2, y⇩1) = cbox (x⇩0, y⇩0) (x⇩2, y⇩1)"
  using assms by auto

lemma cbox_max_dist:
  assumes "p⇩0 = (x, y)" "p⇩1 = (x + δ, y + δ)"
  assumes "(x⇩0, y⇩0) ∈ cbox p⇩0 p⇩1" "(x⇩1, y⇩1) ∈ cbox p⇩0 p⇩1" "0 ≤ δ"
  shows "dist (x⇩0, y⇩0) (x⇩1, y⇩1) ≤ sqrt 2 * δ"
proof -
  have X: "dist x⇩0 x⇩1 ≤ δ"
    using assms dist_real_def by auto
  have Y: "dist y⇩0 y⇩1 ≤ δ"
    using assms dist_real_def by auto

  have "dist (x⇩0, y⇩0) (x⇩1, y⇩1) = sqrt ((dist x⇩0 x⇩1)⇧2 + (dist y⇩0 y⇩1)⇧2)"
    using dist_Pair_Pair by auto
  also have "... ≤ sqrt (δ⇧2 + (dist y⇩0 y⇩1)⇧2)"
    using X power_mono by fastforce
  also have "... ≤ sqrt (δ⇧2 + δ⇧2)"
    using Y power_mono by fastforce
  also have "... = sqrt 2 * sqrt (δ⇧2)"
    using real_sqrt_mult by simp
  also have "... = sqrt 2 * δ"
    using assms(5) by simp
  finally show ?thesis .
qed

subsubsection "Pigeonhole Argument"

lemma card_le_1_if_pairwise_eq:
  assumes "∀x ∈ S. ∀y ∈ S. x = y"
  shows "card S ≤ 1"
proof (rule ccontr)
  assume "¬ card S ≤ 1"
  hence "2 ≤ card S"
    by simp
  then obtain T where *: "T ⊆ S ∧ card T = 2"
    using ex_card by metis
  then obtain x y where "x ∈ T ∧ y ∈ T ∧ x ≠ y"
    by (meson card_2_iff')
  then show False
    using * assms by blast
qed

lemma card_Int_if_either_in:
  assumes "∀x ∈ S. ∀y ∈ S. x = y ∨ x ∉ T ∨ y ∉ T" 
  shows "card (S ∩ T) ≤ 1"
proof (rule ccontr)
  assume "¬ (card (S ∩ T) ≤ 1)"
  then obtain x y where *: "x ∈ S ∩ T ∧ y ∈ S ∩ T ∧ x ≠ y"
    by (meson card_le_1_if_pairwise_eq)
  hence "x ∈ T" "y ∈ T"
    by simp_all
  moreover have "x ∉ T ∨ y ∉ T"
    using assms * by auto
  ultimately show False
    by blast
qed

lemma card_Int_Un_le_Sum_card_Int:
  assumes "finite S"
  shows "card (A ∩ ⋃S) ≤ (∑B ∈ S. card (A ∩ B))"
  using assms
proof (induction "card S" arbitrary: S)
  case (Suc n)
  then obtain B T where *: "S = { B } ∪ T" "card T = n" "B ∉ T"
    by (metis card_Suc_eq Suc_eq_plus1 insert_is_Un)
  hence "card (A ∩ ⋃S) = card (A ∩ ⋃({ B } ∪ T))"
    by blast
  also have "... ≤ card (A ∩ B) + card (A ∩ ⋃T)"
    by (simp add: card_Un_le inf_sup_distrib1)
  also have "... ≤ card (A ∩ B) + (∑B ∈ T. card (A ∩ B))"
    using Suc * by simp
  also have "... ≤ (∑B ∈ S. card (A ∩ B))"
    using Suc.prems * by simp
  finally show ?case .
qed simp

lemma pigeonhole:
  assumes "finite T" "S ⊆ ⋃T" "card T < card S"
  shows "∃x ∈ S. ∃y ∈ S. ∃X ∈ T. x ≠ y ∧ x ∈ X ∧ y ∈ X"
proof (rule ccontr)
  assume "¬ (∃x ∈ S. ∃y ∈ S. ∃X ∈ T. x ≠ y ∧ x ∈ X ∧ y ∈ X)"
  hence *: "∀X ∈ T. card (S ∩ X) ≤ 1"
    using card_Int_if_either_in by metis
  have "card T < card (S ∩ ⋃T)"
    using Int_absorb2 assms by fastforce
  also have "... ≤ (∑X ∈ T. card (S ∩ X))"
    using assms(1) card_Int_Un_le_Sum_card_Int by blast
  also have "... ≤ card T"
    using * sum_mono by fastforce
  finally show False by simp
qed

subsubsection "Delta Sparse Points within a Square"

lemma max_points_square:
  assumes "∀p ∈ ps. p ∈ cbox (x, y) (x + δ, y + δ)" "sparse δ ps" "0 ≤ δ"
  shows "card ps ≤ 4"
proof (cases "δ = 0")
  case True
  hence "{ (x, y) } = cbox (x, y) (x + δ, y + δ)"
    using cbox_def by simp
  hence "∀p ∈ ps. p = (x, y)"
    using assms(1) by blast
  hence "∀p ∈ ps. ∀q ∈ ps. p = q"
    apply (auto split: prod.splits) by (metis of_int_eq_iff)+
  thus ?thesis
    using card_le_1_if_pairwise_eq by force
next
  case False
  hence δ: "0 < δ"
    using assms(3) by simp
  show ?thesis
  proof (rule ccontr)
    assume A: "¬ (card ps ≤ 4)"
    define PS where PS_def: "PS = (λ(x, y). (real_of_int x, real_of_int y)) ` ps"
    have "inj_on (λ(x, y). (real_of_int x, real_of_int y)) ps"
      using inj_on_def by fastforce
    hence *: "¬ (card PS ≤ 4)"
      using A PS_def by (simp add: card_image)

    let ?x' = "x + δ / 2"
    let ?y' = "y + δ / 2"

    let ?ll = "cbox ( x ,  y ) (?x'   , ?y'   )"
    let ?lu = "cbox ( x , ?y') (?x'   ,  y + δ)"
    let ?rl = "cbox (?x',  y ) ( x + δ, ?y'   )"
    let ?ru = "cbox (?x', ?y') ( x + δ,  y + δ)"

    let ?sq = "{ ?ll, ?lu, ?rl, ?ru }"

    have card_le_4: "card ?sq ≤ 4"
      by (simp add: card_insert_le_m1)

    have "cbox (x, y) (?x', y + δ) = ?ll ∪ ?lu"
      using cbox_top_un assms(3) by auto
    moreover have "cbox (?x', y) (x + δ, y + δ) = ?rl ∪ ?ru"
      using cbox_top_un assms(3) by auto
    moreover have "cbox (x, y) (?x', y + δ) ∪ cbox (?x', y) (x + δ, y + δ) = cbox (x, y) (x + δ, y + δ)"
      using cbox_right_un assms(3) by simp
    ultimately have "?ll ∪ ?lu ∪ ?rl ∪ ?ru = cbox (x, y) (x + δ, y + δ)"
      by blast

    hence "PS ⊆ ⋃(?sq)"
      using assms(1) PS_def by (auto split: prod.splits)
    moreover have "card ?sq < card PS"
      using * card_insert_le_m1 card_le_4 by linarith
    moreover have "finite ?sq"
      by simp
    ultimately have "∃p⇩0 ∈ PS. ∃p⇩1 ∈ PS. ∃S ∈ ?sq. (p⇩0 ≠ p⇩1 ∧ p⇩0 ∈ S ∧ p⇩1 ∈ S)"
      using pigeonhole[of ?sq PS] by blast
    then obtain S p⇩0 p⇩1 where #: "p⇩0 ∈ PS" "p⇩1 ∈ PS" "S ∈ ?sq" "p⇩0 ≠ p⇩1" "p⇩0 ∈ S" "p⇩1 ∈ S"
      by blast

    have D: "0 ≤ δ / 2"
      using assms(3) by simp
    have LL: "∀p⇩0 ∈ ?ll. ∀p⇩1 ∈ ?ll. dist p⇩0 p⇩1 ≤ sqrt 2 * (δ / 2)"
      using cbox_max_dist[of "(x, y)" x y "(?x', ?y')" "δ / 2"] D by auto
    have LU: "∀p⇩0 ∈ ?lu. ∀p⇩1 ∈ ?lu. dist p⇩0 p⇩1 ≤ sqrt 2 * (δ / 2)"
      using cbox_max_dist[of "(x, ?y')" x ?y' "(?x', y + δ)" "δ / 2"] D by auto
    have RL: "∀p⇩0 ∈ ?rl. ∀p⇩1 ∈ ?rl. dist p⇩0 p⇩1 ≤ sqrt 2 * (δ / 2)"
      using cbox_max_dist[of "(?x', y)" ?x' y "(x + δ, ?y')" "δ / 2"] D by auto
    have RU: "∀p⇩0 ∈ ?ru. ∀p⇩1 ∈ ?ru. dist p⇩0 p⇩1 ≤ sqrt 2 * (δ / 2)"
      using cbox_max_dist[of "(?x', ?y')" ?x' ?y' "(x + δ, y + δ)" "δ / 2"] D by auto

    have "∀p⇩0 ∈ S. ∀p⇩1 ∈ S. dist p⇩0 p⇩1 ≤ sqrt 2 * (δ / 2)"
      using # LL LU RL RU by blast
    hence "dist p⇩0 p⇩1 ≤ (sqrt 2 / 2) * δ"
      using # by simp
    moreover have "(sqrt 2 / 2) * δ < δ"
      using sqrt2_less_2 δ by simp
    ultimately have "dist p⇩0 p⇩1 < δ"
      by simp
    moreover have "δ ≤ dist p⇩0 p⇩1"
      using assms(2) # sparse_def PS_def by auto
    ultimately show False
      by simp
  qed
qed

end