Theory Closest_Pair_Points.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 (λp0 p1. fst p0  fst p1) ps"

definition sorted_snd :: "point list  bool" where
  "sorted_snd ps = sorted_wrt (λp0 p1. snd p0  snd p1) 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;
      (xsl, xsr) <- split_at_tm (n div 2) xs;
      l <- mergesort_tm f xsl;
      r <- mergesort_tm f xsr;
      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 "λp0 p1. fst p0  fst p1"] 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 "λp0 p1. fst p0  fst p1"] 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  (p0  ps. p1  ps. p0  p1  δ  dist p0 p1)"

lemma sparse_identity:
  assumes "sparse δ (set ps)" "p  set ps. δ  dist p0 p"
  shows "sparse δ (set (p0 # ps))"
  using assms by (simp add: dist_commute sparse_def)

lemma sparse_update:
  assumes "sparse δ (set ps)"
  assumes "dist p0 p1  δ" "p  set ps. dist p0 p1  dist p0 p"
  shows "sparse (dist p0 p1) (set (p0 # 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 p0 p1 = (fst p0 - fst p1)2 + (snd p0 - snd p1)2"

lemma dist_eq_sqrt_dist_code:
  fixes p0 :: point
  shows "dist p0 p1 = sqrt (dist_code p0 p1)"
  by (auto simp: dist_prod_def dist_real_def split: prod.splits) 

lemma dist_eq_dist_code_lt:
  fixes p0 :: point
  shows "dist p0 p1 < dist p2 p3  dist_code p0 p1 < dist_code p2 p3"
  using dist_eq_sqrt_dist_code real_sqrt_less_iff by presburger

lemma dist_eq_dist_code_le:
  fixes p0 :: point
  shows "dist p0 p1  dist p2 p3  dist_code p0 p1  dist_code p2 p3"
  using dist_eq_sqrt_dist_code real_sqrt_le_iff by presburger

lemma dist_eq_dist_code_abs_lt:
  fixes p0 :: point
  shows "¦c¦ < dist p0 p1  c2 < dist_code p0 p1"
  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 p0 :: point
  shows "dist p0 p1  ¦c¦  dist_code p0 p1  c2"
  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 (p0 # ps) =1 (
    do {
      p1 <- find_closest_bf_tm p ps;
      if dist p p0 < dist p p1 then
        return p0
      else
        return p1
    }
  )"

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

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 [p0, p1] =1 return (p0, p1)"
| "closest_pair_bf_tm (p0 # ps) =1 (
    do {
      (c0::point, c1::point) <- closest_pair_bf_tm ps;
      p1 <- find_closest_bf_tm p0 ps;
      if dist c0 c1  dist p0 p1 then
        return (c0, c1)
      else
        return (p0, p1)
    }
  )"

fun closest_pair_bf :: "point list  (point * point)" where
  "closest_pair_bf [] = undefined"
| "closest_pair_bf [_] = undefined"
| "closest_pair_bf [p0, p1] = (p0, p1)"
| "closest_pair_bf (p0 # ps) = (
    let (c0, c1) = closest_pair_bf ps in
    let p1 = find_closest_bf p0 ps in
    if dist c0 c1  dist p0 p1 then
      (c0, c1)
    else
      (p0, p1) 
  )"

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  (c0, c1) = closest_pair_bf ps  c0  set ps"
  by (induction ps arbitrary: c0 c1 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  (c0, c1) = closest_pair_bf ps  c1  set ps" 
proof (induction ps arbitrary: c0 c1 rule: closest_pair_bf.induct)
  case (4 p0 p2 p3 ps)
  let ?ps = "p2 # p3 # ps"
  obtain c0 c1 where c0_def: "(c0, c1) = closest_pair_bf ?ps"
    using prod.collapse by blast
  define p1 where p1_def: "p1 = find_closest_bf p0 ?ps"
  note defs = c0_def p1_def
  have "c1  set ?ps"
    using "4.IH" defs by simp
  moreover have "p1  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  (c0, c1) = closest_pair_bf ps  c0  c1"
proof (induction ps arbitrary: c0 c1 rule: closest_pair_bf.induct)
  case (4 p0 p2 p3 ps)
  let ?ps = "p2 # p3 # ps"
  obtain c0 c1 where c0_def: "(c0, c1) = closest_pair_bf ?ps"
    using prod.collapse by blast
  define p1 where p1_def: "p1 = find_closest_bf p0 ?ps"
  note defs = c0_def p1_def
  have "c0  c1"
    using "4.IH" "4.prems"(2) defs by simp
  moreover have "p0  p1"
    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" "(c0, c1) = closest_pair_bf ps"
  shows "sparse (dist c0 c1) (set ps)"
  using assms
proof (induction ps arbitrary: c0 c1 rule: closest_pair_bf.induct)
  case (4 p0 p2 p3 ps)
  let ?ps = "p2 # p3 # ps"
  obtain c0 c1 where c0_def: "(c0, c1) = closest_pair_bf ?ps"
    using prod.collapse by blast
  define p1 where p1_def: "p1 = find_closest_bf p0 ?ps"
  note defs = c0_def p1_def
  hence IH: "sparse (dist c0 c1) (set ?ps)"
    using 4 c0_def by simp
  have *: "p  set ?ps. (dist p0 p1)  dist p0 p"
    using find_closest_bf_dist defs by blast
  show ?case
  proof (cases "dist c0 c1  dist p0 p1")
    case True
    hence "p  set ?ps. dist c0 c1  dist p0 p"
      using * by auto
    hence "sparse (dist c0 c1) (set (p0 # ?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 p0 p1) (set (p0 # ?ps))"
      using sparse_update[of "dist c0 c1" ?ps p0 p1] 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 p0 p2 p3 ps)
  let ?ps = "p2 # p3 # ps"
  have "time (closest_pair_bf_tm (p0 # ?ps)) = 1 + time (find_closest_bf_tm p0 ?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 p0 ?ps] by simp
  also have "...  2 + length ?ps + length ?ps * length ?ps + 1"
    using "4.IH" by simp
  also have "...  length (p0 # ?ps) * length (p0 # ?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 [p0] = (dist_code p p0, p0)"
| "find_closest_bf_code p (p0 # ps) = (
    let (δ1, p1) = find_closest_bf_code p ps in
    let δ0 = dist_code p p0 in
    if δ0 < δ1 then
      (δ0, p0)
    else
      (δ1, p1)
  )"

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 p0 p2 ps)
  define δ0 δ0' where δ0_def: "δ0 = dist p p0" "δ0' = dist_code p p0"
  obtain δ1 p1 δ1' p1' where δ1_def: "δ1 = dist p p1" "p1 = find_closest_bf p (p2 # ps)"
    "(δ1', p1') = find_closest_bf_code p (p2 # ps)"
    using prod.collapse by blast+
  note defs = δ0_def δ1_def
  have *: "p1 = p1'"
    using "3.IH" defs by simp
  hence "δ0 < δ1  δ0' < δ1'"
    using find_closest_bf_code_dist_eq[of "p2 # ps" δ1' p1' 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 [p0] = undefined"
| "closest_pair_bf_code [p0, p1] = (dist_code p0 p1, p0, p1)"
| "closest_pair_bf_code (p0 # ps) = (
    let (δc, c0, c1) = closest_pair_bf_code ps in
    let (δp, p1) = find_closest_bf_code p0 ps in
    if δc  δp then
      (δc, c0, c1)
    else
      (δp, p0, p1) 
  )"

lemma closest_pair_bf_code_dist_eq:
  "1 < length ps  (δ, c0, c1) = closest_pair_bf_code ps  δ = dist_code c0 c1"
proof (induction ps arbitrary: δ c0 c1 rule: closest_pair_bf_code.induct)
  case (4 p0 p2 p3 ps)
  let ?ps = "p2 # p3 # ps"
  obtain δc c0 c1 where δc_def: "(δc, c0, c1) = closest_pair_bf_code ?ps"
    by (metis prod_cases3)
  obtain δp p1 where δp_def: "(δp, p1) = find_closest_bf_code p0 ?ps"
    using prod.collapse by blast
  note defs = δc_def δp_def
  have "δc = dist_code c0 c1"
    using "4.IH" defs by simp
  moreover have "δp = dist_code p0 p1"
    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 "(c0, c1) = closest_pair_bf ps" "(δ', c0', c1') = closest_pair_bf_code ps"
  shows "c0 = c0'  c1 = c1'"
  using assms
proof (induction ps arbitrary: c0 c1 δ' c0' c1' rule: closest_pair_bf_code.induct)
  case (4 p0 p2 p3 ps)
  let ?ps = "p2 # p3 # ps"
  obtain c0 c1 δc' c0' c1' where δc_def: "(c0, c1) = closest_pair_bf ?ps"
    "(δc', c0', c1') = closest_pair_bf_code ?ps"
    by (metis prod_cases3)
  obtain p1 δp' p1' where δp_def: "p1 = find_closest_bf p0 ?ps"
    "(δp', p1') = find_closest_bf_code p0 ?ps"
    using prod.collapse by blast
  note defs = δc_def δp_def
  have A: "c0 = c0'  c1 = c1'"
    using "4.IH" defs by simp
  moreover have B: "p1 = p1'"
    using find_closest_bf_code_eq defs by blast
  moreover have "δc' = dist_code c0' c1'"
    using defs closest_pair_bf_code_dist_eq[of ?ps] by simp
  moreover have "δp' = dist_code p0 p1'"
    using defs find_closest_bf_code_dist_eq by blast
  ultimately have "dist c0 c1  dist p0 p1  δ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 "p0  psL" "p1  psR" "p0  p1" "dist p0 p1 < δ" "set ps = psL  psR"
  assumes "p  psL. fst p  l" "p  psR. l  fst p"
  assumes "ps' = filter (λp. l - δ < fst p  fst p < l + δ) ps"
  shows "p0  set ps'  p1  set ps'"
proof (rule ccontr)
  assume "¬ (p0  set ps'  p1  set ps')"
  then consider (A) "p0  set ps'  p1  set ps'"
              | (B) "p0  set ps'  p1  set ps'"
              | (C) "p0  set ps'  p1  set ps'"
    by blast
  thus False
  proof cases
    case A
    hence "fst p0  l - δ  l + δ  fst p0" "fst p1  l - δ  l + δ  fst p1"
      using assms(1,2,5,8) by auto
    hence "fst p0  l - δ" "l + δ  fst p1"
      using assms(1,2,6,7) by force+
    hence "δ  dist (fst p0) (fst p1)"
      using dist_real_def by simp
    hence "δ  dist p0 p1"
      using dist_fst_le[of p0 p1] by (auto split: prod.splits)
    then show ?thesis
      using assms(4) by fastforce
  next
    case B
    hence "fst p1  l - δ  l + δ  fst p1"
      using assms(2,5,8) by auto
    hence "l + δ  fst p1"
      using assms(2,7) by auto
    moreover have "fst p0  l"
      using assms(1,6) by simp
    ultimately have "δ  dist (fst p0) (fst p1)"
      using dist_real_def by simp
    hence "δ  dist p0 p1"
      using dist_fst_le[of p0 p1] less_le_trans by (auto split: prod.splits)
    thus ?thesis
      using assms(4) by simp
  next
    case C
    hence "fst p0  l - δ  l + δ  fst p0"
      using assms(1,2,5,8) by auto
    hence "fst p0  l - δ"
      using assms(1,6) by auto
    moreover have "l  fst p1"
      using assms(2,7) by simp
    ultimately have "δ  dist (fst p0) (fst p1)"
      using dist_real_def by simp
    hence "δ  dist p0 p1"
      using dist_fst_le[of p0 p1] 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 "p0  set ps" "p1  set ps" "p0  p1" "dist p0 p1 < δ" "set ps = psL  psR"
  assumes "sparse δ psL" "sparse δ psR"
  assumes "p  psL. fst p  l" "p  psR. l  fst p"
  assumes "ps' = filter (λp. l - δ < fst p  fst p < l + δ) ps"
  shows "p0  set ps'  p1  set ps'"
proof -
  have "p0  psL  p1  psL" "p0  psR  p1  psR"
    using assms(3,4,6,7) sparse_def by force+
  then consider (A) "p0  psL  p1  psR" | (B) "p0  psR  p1  psL"
    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 p1 p0 < δ"
      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 x0 :: real and y0 :: real
  shows "cbox (x0, y0) (x1, y1) = { (x, y). x0  x  x  x1  y0  y  y  y1 }"
  by fastforce

lemma mem_cbox_2D:
  fixes x :: real and y :: real
  shows "x0  x  x  x1  y0  y  y  y1  (x, y)  cbox (x0, y0) (x1, y1)"
  by fastforce

lemma cbox_top_un:
  fixes x0 :: real and y0 :: real
  assumes "y0  y1" "y1  y2"
  shows "cbox (x0, y0) (x1, y1)  cbox (x0, y1) (x1, y2) = cbox (x0, y0) (x1, y2)"
  using assms by auto

lemma cbox_right_un:
  fixes x0 :: real and y0 :: real
  assumes "x0  x1" "x1  x2"
  shows "cbox (x0, y0) (x1, y1)  cbox (x1, y0) (x2, y1) = cbox (x0, y0) (x2, y1)"
  using assms by auto

lemma cbox_max_dist:
  assumes "p0 = (x, y)" "p1 = (x + δ, y + δ)"
  assumes "(x0, y0)  cbox p0 p1" "(x1, y1)  cbox p0 p1" "0  δ"
  shows "dist (x0, y0) (x1, y1)  sqrt 2 * δ"
proof -
  have X: "dist x0 x1  δ"
    using assms dist_real_def by auto
  have Y: "dist y0 y1  δ"
    using assms dist_real_def by auto

  have "dist (x0, y0) (x1, y1) = sqrt ((dist x0 x1)2 + (dist y0 y1)2)"
    using dist_Pair_Pair by auto
  also have "...  sqrt (δ2 + (dist y0 y1)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 "p0  PS. p1  PS. S  ?sq. (p0  p1  p0  S  p1  S)"
      using pigeonhole[of ?sq PS] by blast
    then obtain S p0 p1 where #: "p0  PS" "p1  PS" "S  ?sq" "p0  p1" "p0  S" "p1  S"
      by blast

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

    have "p0  S. p1  S. dist p0 p1  sqrt 2 * (δ / 2)"
      using # LL LU RL RU by blast
    hence "dist p0 p1  (sqrt 2 / 2) * δ"
      using # by simp
    moreover have "(sqrt 2 / 2) * δ < δ"
      using sqrt2_less_2 δ by simp
    ultimately have "dist p0 p1 < δ"
      by simp
    moreover have "δ  dist p0 p1"
      using assms(2) # sparse_def PS_def by auto
    ultimately show False
      by simp
  qed
qed

end