(* ========================================================================= *)
(*                FLYSPECK - BOOK FORMALIZATION                              *)
(*                                                                           *)
(*      Authour   : VU KHAC KY                                               *)
(*      Book lemma: LEPJBDJ                                                  *)
(*      Chaper    : Packing (Marchal cells)                                  *)
(*                                                                           *)
(* ========================================================================= *)

(* ========================================================================= *)
(*                     FILES NEED TO BE LOADED                               *)
     
(*  flyspeck_needs "packing/marchal_cells_2.hl";;                            *)

(* ========================================================================= *)
module Lepjbdj = struct

open Rogers;;
(* open Prove_by_refinement;; *)
open Vukhacky_tactics;;
open Pack_defs;;
open Pack_concl;; 
open Pack1;;
open Sphere;; 
open Marchal_cells;;
open Emnwuus;;
open Marchal_cells_2_new;;

let LEPJBDJ_concl = `!V ul k.  
  saturated V /\ packing V /\ barV V 3 ul /\ (1 <= k) /\ (k <= 4) /\ 
  ~(mcell k V ul = {})==>
  (V INTER mcell k V ul = set_of_list (truncate_simplex (k-1) ul))`;;

let LEPJBDJ_0_concl =
 `!V ul.
         saturated V /\
         packing V /\
         barV V 3 ul
        ==> V INTER mcell 0 V ul = {}`;;

(* ========================================================================= *)

let LEPJBDJ = prove_by_refinement (
 LEPJBDJ_concl,
[(REPEAT STRIP_TAC);
 (NEW_GOAL `?u0 u1 u2 u3. ul = [u0;u1;u2;u3:real^3]`);
 (MATCH_MP_TAC BARV_3_EXPLICIT);
 (EXISTS_TAC `V:real^3->bool` THEN ASM_REWRITE_TAC[]);
 (UP_ASM_TAC THEN STRIP_TAC);
 (NEW_GOAL `~(u0 = u1:real^3)`);
 (STRIP_TAC);
 (NEW_GOAL `CARD {u0,u1,u2,u3:real^3} = 4`);
 (REWRITE_WITH `{u0,u1,u2,u3:real^3} = set_of_list ul`);
 (ASM_REWRITE_TAC[set_of_list]);
 (REWRITE_WITH `LENGTH (ul:(real^3)list) = 3 + 1 /\ CARD (set_of_list ul) = 3 + 1`);
 (MATCH_MP_TAC Rogers.BARV_IMP_LENGTH_EQ_CARD);
 (EXISTS_TAC `V:real^3->bool` THEN ASM_REWRITE_TAC[]);
 (ARITH_TAC);
 (UP_ASM_TAC THEN ASM_REWRITE_TAC[SET_RULE `{a,a,b,c} = {a,b,c}`]);
 (NEW_GOAL `CARD {u1,u2,u3:real^3} <= 3`);
 (REWRITE_TAC[Geomdetail.CARD3]);
 (ASM_ARITH_TAC);

(*==== Next, we will consider 4 case of k = 1,2,3 and k >= 4   ======= *)
(* -------------------------------------------------------------------- *)
 (ASM_CASES_TAC `k = 1`);
 (ASM_REWRITE_TAC[ARITH_RULE `1 - 1 = 0`; TRUNCATE_SIMPLEX_EXPLICIT_0; set_of_list; MCELL_EXPLICIT; mcell1; HD;TL; TRUNCATE_SIMPLEX_EXPLICIT_1]);
 (COND_CASES_TAC);
 (REWRITE_TAC[SET_EQ_LEMMA; IN_SING]);
 (REPEAT STRIP_TAC);
 (NEW_GOAL `x IN V INTER cball (u0:real^3,sqrt (&2))`);
 (UP_ASM_TAC THEN SET_TAC[]);
 (UP_ASM_TAC THEN REWRITE_TAC[IN_INTER; cball]);
 (STRIP_TAC THEN UP_ASM_TAC THEN REWRITE_TAC[IN;IN_ELIM_THM] THEN STRIP_TAC);
 (UNDISCH_TAC `packing (V:real^3->bool)` THEN REWRITE_TAC[packing]);
 (STRIP_TAC);
 (ASM_CASES_TAC `x = u0:real^3`);
 (ASM_MESON_TAC[]);
 (NEW_GOAL `F`);
 (NEW_GOAL `&2 <= dist (u0,x:real^3)`);
 (FIRST_ASSUM MATCH_MP_TAC);
 (ASM_REWRITE_TAC[]);
 (STRIP_TAC);
 (ONCE_REWRITE_TAC[GSYM IN] THEN MATCH_MP_TAC (SET_RULE `(?b. a IN b /\ b SUBSET c) ==> a IN c`));
 (EXISTS_TAC `set_of_list (ul:(real^3)list)`);
 (STRIP_TAC);
 (ASM_REWRITE_TAC[set_of_list] THEN SET_TAC[]);
 (MATCH_MP_TAC Packing3.BARV_SUBSET);
 (EXISTS_TAC `3` THEN ASM_REWRITE_TAC[]);
 (UNDISCH_TAC `x:real^3 IN V` THEN REWRITE_TAC[IN]);
 (NEW_GOAL `sqrt (&2) < &2`);
 (REWRITE_TAC[SQRT2_LT_2]);
 (ASM_REAL_ARITH_TAC);
 (ASM_MESON_TAC[]);

 (REWRITE_TAC[IN_INTER; IN_DIFF]);
 (REPEAT STRIP_TAC); (* break into 3 subgoals *)

 (ASM_REWRITE_TAC[]);
 (MATCH_MP_TAC (SET_RULE `(?b. a IN b /\ b SUBSET c) ==> a IN c`));
 (EXISTS_TAC `set_of_list (ul:(real^3)list)`);
 (STRIP_TAC);
 (ASM_REWRITE_TAC[set_of_list] THEN SET_TAC[]);
 (MATCH_MP_TAC Packing3.BARV_SUBSET);
 (EXISTS_TAC `3` THEN ASM_REWRITE_TAC[]);
 (REWRITE_TAC[GSYM (ASSUME `ul = [u0; u1; u2; u3:real^3]`)]);
 (REWRITE_WITH `rogers V ul = convex hull
   {HD ul, omega_list_n V ul 1, omega_list_n V ul 2, omega_list_n V ul 3}`);
 (MATCH_MP_TAC ROGERS_EXPLICIT);
 (ASM_REWRITE_TAC[]);
 (ASM_REWRITE_TAC[HD]);
 (MATCH_MP_TAC IN_SET_IMP_IN_CONVEX_HULL_SET);
 (SET_TAC[]);
 (ASM_REWRITE_TAC[cball; IN; IN_ELIM_THM; DIST_REFL]);
 (MATCH_MP_TAC SQRT_POS_LE);
 (REAL_ARITH_TAC);
 (UP_ASM_TAC THEN ASM_REWRITE_TAC[rcone_gt; rconesgn; IN; IN_ELIM_THM;  
   VECTOR_SUB_REFL; DOT_LZERO; ARITH_RULE `~(a > b) <=> a <= b`; DIST_REFL;  
   REAL_MUL_LZERO; REAL_LE_REFL]);
 (NEW_GOAL `F`);
 (UNDISCH_TAC `~(mcell k V ul = {})`);
 (ASM_REWRITE_TAC[MCELL_EXPLICIT; mcell1]);
 (ASM_MESON_TAC[]);

(* -------------------------------------------------------------------- *)
 (ASM_CASES_TAC `k = 2`);
 (ASM_REWRITE_TAC[ARITH_RULE `2 - 1 = 1`; TRUNCATE_SIMPLEX_EXPLICIT_1;
   MCELL_EXPLICIT; mcell2;HD;TL]);
 (COND_CASES_TAC);
 (LET_TAC);
 (NEW_GOAL `hl [u0;u1:real^3] = inv(&2) * dist (u1,u0)`);
 (REWRITE_WITH `hl [u0;u1:real^3] = 
   dist (circumcenter (set_of_list [u0;u1]),HD [u0;u1])`);
 (MATCH_MP_TAC Rogers.HL_EQ_DIST0);
 (EXISTS_TAC `V:real^3->bool` THEN EXISTS_TAC `1`);
 (ASM_REWRITE_TAC[]);
 (REWRITE_WITH `[u0;u1:real^3] = truncate_simplex 1 ul`);
 (ASM_REWRITE_TAC[TRUNCATE_SIMPLEX_EXPLICIT_1]);
 (MATCH_MP_TAC Packing3.TRUNCATE_SIMPLEX_BARV);
 (EXISTS_TAC `3` THEN ASM_REWRITE_TAC[ARITH_RULE `1 <= 3`]);
 (REWRITE_TAC[set_of_list; Rogers.CIRCUMCENTER_2; midpoint; HD; dist]);
 (NORM_ARITH_TAC);

 (NEW_GOAL `&0 < a /\ a < &1`);
 (EXPAND_TAC "a");
 (STRIP_TAC);
 (MATCH_MP_TAC REAL_LT_DIV);
 (SIMP_TAC[SQRT_POS_LT; ARITH_RULE `&0 < &2`]);
 (ASM_REWRITE_TAC[]);
 (MATCH_MP_TAC REAL_LT_MUL);
 (ASM_SIMP_TAC[REAL_LT_INV; ARITH_RULE `&0 < &2`]);
 (MATCH_MP_TAC DIST_POS_LT THEN ASM_MESON_TAC[]);
 (REWRITE_WITH `hl [u0; u1:real^3] / sqrt (&2) < &1 <=> 
                 hl [u0; u1] < &1 * sqrt (&2)`);
 (MATCH_MP_TAC REAL_LT_LDIV_EQ);
 (SIMP_TAC[SQRT_POS_LT; ARITH_RULE `&0 < &2`]);
 (REWRITE_TAC[REAL_MUL_LID]);
 (ASM_MESON_TAC[]);


 (REWRITE_TAC[set_of_list; SET_EQ_LEMMA]);
 (REPEAT STRIP_TAC);
 (ASM_CASES_TAC `x IN {u0,u1:real^3}`);
 (ASM_REWRITE_TAC[]);
 (NEW_GOAL `F`);
 (SWITCH_TAC THEN UP_ASM_TAC THEN REWRITE_TAC[IN_INTER] THEN STRIP_TAC);


 (ABBREV_TAC `v = proj_point (u1 - u0:real^3) (x - u0) `);
 (NEW_GOAL `(x - u0 - v:real^3) dot (u1 - u0) = &0`);
 (EXPAND_TAC "v");
 (REWRITE_TAC[GSYM projection_proj_point]);
 (REWRITE_TAC[Packing3.PROJECTION_ORTHOGONAL]);

 (NEW_GOAL `between (v+u0) (u0,u1:real^3)`);
 (ASM_CASES_TAC `between (v+u0) (u0,u1:real^3)`);
 (ASM_REWRITE_TAC[]);
 (NEW_GOAL `between (v+u0) (u0,u1:real^3) \/ between u0 (u1,v+u0) \/ between u1 (v+u0,u0)`);
 (REWRITE_TAC[GSYM COLLINEAR_BETWEEN_CASES]);
 (MATCH_MP_TAC Collect_geom.IN_AFFINE_HULL_IMP_COLLINEAR);
 (REWRITE_TAC[AFFINE_HULL_2; IN_ELIM_THM; IN]);
 (EXPAND_TAC "v");
 (REWRITE_TAC[PRO_EXP]);
 (ABBREV_TAC `t = ((x - u0) dot (u1 - u0)) / ((u1 - u0) dot (u1 - u0):real^3)`);
 (EXISTS_TAC `(&1 - t)` THEN EXISTS_TAC `t:real`);
 (REWRITE_TAC[REAL_ARITH `(&1 - t) + t = &1`]);
 (VECTOR_ARITH_TAC);
 (NEW_GOAL `F`);
 (UP_ASM_TAC THEN ASM_REWRITE_TAC[] THEN STRIP_TAC);

 (UNDISCH_TAC `x:real^3 IN rcone_ge u0 u1 a`);
 (REWRITE_TAC[rcone_ge;rconesgn;IN;IN_ELIM_THM; ARITH_RULE `~(a >= b) <=> a < b`]);
 (NEW_GOAL `(x - u0) dot (u1 - u0:real^3) <= &0`);
 (UP_ASM_TAC THEN REWRITE_TAC[BETWEEN_IN_CONVEX_HULL; CONVEX_HULL_2;  
   IN;IN_ELIM_THM] THEN STRIP_TAC);
 (REWRITE_WITH `(x - u0) dot (u1 - u0:real^3) = (x - u0 - v) dot (u1 - u0) + 
   v dot (u1 - u0)`);
 (VECTOR_ARITH_TAC);
 (REWRITE_TAC[ASSUME `(x - u0 - v) dot (u1 - u0:real^3) = &0`; REAL_ADD_LID]);
 (MP_TAC (ASSUME `u0 = u % u1 + v' % (v + u0:real^3)`) THEN 
   REWRITE_WITH `u0 = u % u1 + v' % (v + u0:real^3) 
   <=> (u + v') % u0 = u % u1 + v' % (v + u0)`);
 (ASM_REWRITE_TAC[]);
 (VECTOR_ARITH_TAC);
 (REWRITE_TAC[VECTOR_ARITH `(u + v') % u0 = u % u1 + v' % (v + u0) 
  <=> u % (u0 - u1) = v' % v`]);
 (STRIP_TAC);
 (UP_ASM_TAC);
 (REWRITE_WITH `u % (u0 - u1) = v' % v <=> u / v' % (u0 - u1:real^3) = v`);
 (ONCE_REWRITE_TAC[EQ_SYM_EQ]);
 (MATCH_MP_TAC Trigonometry2.VEC_DIV_MOV);
 (STRIP_TAC);
 (MP_TAC (ASSUME `u0 = u % u1 + v' % (v + u0:real^3)`));
 (REWRITE_WITH `u = &1 /\ v' = &0`);
 (ASM_REAL_ARITH_TAC);
 (ASM_REWRITE_TAC[VECTOR_MUL_LID; VECTOR_MUL_LZERO; VECTOR_ADD_RID]);
 (ONCE_REWRITE_TAC[EQ_SYM_EQ] THEN STRIP_TAC);
 (ASM_REWRITE_TAC[]);
 (REWRITE_TAC[DOT_LMUL]);
 (REWRITE_WITH `(u0 - u1) dot (u1 - u0) = --((u0 - u1:real^3) dot (u0 - u1))`);
 (VECTOR_ARITH_TAC);
 (REWRITE_TAC[DOT_SQUARE_NORM]);
 (REWRITE_TAC[REAL_ARITH `a * -- b <= &0 <=> &0 <= a * b`]);
 (MATCH_MP_TAC REAL_LE_MUL);
 (REWRITE_TAC[REAL_LE_POW_2]);
 (MATCH_MP_TAC REAL_LE_DIV);
 (ASM_REWRITE_TAC[]);
 (NEW_GOAL `&0 < dist (x,u0) * dist (u1,u0:real^3) * a`);
 (MATCH_MP_TAC REAL_LT_MUL);
 (STRIP_TAC);
 (MATCH_MP_TAC DIST_POS_LT);
 (ASM_SET_TAC[]);
 (MATCH_MP_TAC REAL_LT_MUL);
 (STRIP_TAC);
 (MATCH_MP_TAC DIST_POS_LT);
 (ASM_SET_TAC[]);
 (ASM_REWRITE_TAC[]);
 (ASM_REAL_ARITH_TAC);

(* ==== *)

 (UNDISCH_TAC `x:real^3 IN rcone_ge u1 u0 a`);
 (REWRITE_TAC[rcone_ge;rconesgn;IN;IN_ELIM_THM; ARITH_RULE `~(a >= b) <=> a < b`]);
 (NEW_GOAL `(x - u1) dot (u0 - u1:real^3) <= &0`);
 (UP_ASM_TAC THEN REWRITE_TAC[BETWEEN_IN_CONVEX_HULL; CONVEX_HULL_2;  
   IN;IN_ELIM_THM] THEN STRIP_TAC);
 (REWRITE_WITH `(x - u1) dot (u0 - u1:real^3) = 
  --((x - u0 - v) dot (u1 - u0)) +  (v + u0 - u1) dot (u0 - u1)`);
 (VECTOR_ARITH_TAC);
 (REWRITE_TAC[ASSUME `(x - u0 - v) dot (u1 - u0:real^3) = &0`; 
   REAL_NEG_0; REAL_ADD_LID]);


 (MP_TAC (ASSUME `u1 = u % (v + u0) + v' % u0:real^3`) THEN 
   REWRITE_WITH `u1 = u % (v + u0) + v' % u0 
   <=> u1 = u % v + (u + v') % u0:real^3`);
 (VECTOR_ARITH_TAC);
 (REWRITE_TAC[ASSUME `u + v' = &1`; VECTOR_ARITH `&1 % a = a`]);
 (REWRITE_TAC[VECTOR_ARITH `u1 = u % v + u0 <=> &1 % (u1 - u0:real^3) = u % v`]);
 (REWRITE_WITH `&1 % (u1 - u0) = u % v <=> &1 / u % (u1 - u0:real^3) = v`);

 (ONCE_REWRITE_TAC[EQ_SYM_EQ]);
 (MATCH_MP_TAC Trigonometry2.VEC_DIV_MOV);
 (STRIP_TAC);
 (MP_TAC (ASSUME `u1 = u % (v + u0) + v' % u0:real^3`));
 (REWRITE_WITH `u = &0 /\ v' = &1`);
 (ASM_REAL_ARITH_TAC);
 (REWRITE_TAC[VECTOR_MUL_LID; VECTOR_MUL_LZERO; VECTOR_ADD_LID]);
 (ASM_MESON_TAC[]);
 (ONCE_REWRITE_TAC[EQ_SYM_EQ] THEN STRIP_TAC);
 (REWRITE_WITH `v + u0 - u1 = (&1 / u - &1) % (u1 - u0:real^3)`);

 (REWRITE_TAC[ASSUME `v = &1 / u % (u1 - u0:real^3)`]);
 (VECTOR_ARITH_TAC);
 (REWRITE_TAC[DOT_LMUL]);

 (REWRITE_WITH `(u1 - u0) dot (u0 - u1) = --((u0 - u1:real^3) dot (u0 - u1))`);
 (VECTOR_ARITH_TAC);
 (REWRITE_TAC[DOT_SQUARE_NORM]);

 (REWRITE_TAC[REAL_ARITH `a * -- b <= &0 <=> &0 <= a * b`]);

 (MATCH_MP_TAC REAL_LE_MUL);
 (REWRITE_TAC[REAL_LE_POW_2]);

 (REWRITE_TAC[REAL_ARITH `&0 <= a - b <=> b <= a`]);

 (REWRITE_WITH `&1 <= &1 / u <=> &1 * u <= &1`);
 (MATCH_MP_TAC REAL_LE_RDIV_EQ);
 (ASM_CASES_TAC `u = &0`);
 (NEW_GOAL `F`);
 (MP_TAC (ASSUME `u1 = u % (v + u0) + v' % u0:real^3`));
 (REWRITE_WITH `u = &0 /\ v' = &1`);
 (ASM_REAL_ARITH_TAC);
 (REWRITE_TAC[VECTOR_MUL_LID; VECTOR_MUL_LZERO; VECTOR_ADD_LID]);
 (ASM_MESON_TAC[]);
 (ASM_MESON_TAC[]);
 (ASM_REAL_ARITH_TAC);
 (REWRITE_TAC[REAL_MUL_LID]);
 (ASM_REAL_ARITH_TAC);

 (NEW_GOAL `&0 < dist (x,u1) * dist (u0,u1:real^3) * a`);
 (MATCH_MP_TAC REAL_LT_MUL);
 (STRIP_TAC);
 (MATCH_MP_TAC DIST_POS_LT);
 (ASM_SET_TAC[]);
 (MATCH_MP_TAC REAL_LT_MUL);
 (STRIP_TAC);
 (MATCH_MP_TAC DIST_POS_LT);
 (ASM_SET_TAC[]);
 (ASM_REWRITE_TAC[]);
 (ASM_REAL_ARITH_TAC);
 (ASM_MESON_TAC[]);
(* Finish proving that (v + u0) is between u0 and u1 *)
(* ================================================= *)

 (NEW_GOAL `norm (x - u0) pow 2 = norm ((u0 + v) - u0) pow 2 + norm (x - (u0 + v):real^3) pow 2`);
 (MATCH_MP_TAC PYTHAGORAS);
 (REWRITE_TAC[orthogonal]);
 (REWRITE_TAC[VECTOR_ARITH `a:real^3 - (a + b) = --b`]);
 (REWRITE_TAC[VECTOR_ARITH `--v dot (x - (u0 + v)) = -- ((x - u0 - v) dot v)`]);
 (REWRITE_WITH `(x - u0 - v) dot v = (x - u0 - v) dot proj_point (u1 - u0) (x - u0:real^3)`);
 (ASM_REWRITE_TAC[]);
 (REWRITE_TAC[PRO_EXP]);
 (REWRITE_TAC[DOT_RMUL]);
 (ASM_REWRITE_TAC[] THEN REAL_ARITH_TAC);

 (NEW_GOAL `norm (x - u1) pow 2 = norm ((u0 + v) - u1) pow 2 + norm (x - (u0 + v):real^3) pow 2`);
 (MATCH_MP_TAC PYTHAGORAS);
 (REWRITE_TAC[orthogonal]);
 (REWRITE_TAC[VECTOR_ARITH `(u1 - (u0 + v)) dot (x - (u0 + v)) = ((x - u0 - v) dot (u1 - (u0 + v)))`]);
 (REWRITE_WITH `(x - u0 - v) dot (u1 - (u0 + v)) = 
  (x - u0 - v) dot (u1 - u0 - proj_point (u1 - u0) (x - u0:real^3))`);
 (ASM_REWRITE_TAC[]);
 (VECTOR_ARITH_TAC);
 (REWRITE_TAC[PRO_EXP]);
 (ABBREV_TAC `t = ((x - u0) dot (u1 - u0)) / ((u1 - u0) dot (u1 - u0:real^3))`);
 (REWRITE_TAC[VECTOR_ARITH `a - b - t% (a - b) = (&1 - t) % (a - b)`]);
 (REWRITE_TAC[DOT_RMUL]);
 (ASM_REWRITE_TAC[] THEN REAL_ARITH_TAC);

(* important lemma *)
 (NEW_GOAL `sqrt (&2) <= dist (u0 + v,u0:real^3)`);
 (UNDISCH_TAC `x:real^3 IN rcone_ge u0 u1 a`);
 (REWRITE_TAC[rcone_ge; rconesgn; IN; IN_ELIM_THM]);
 (REWRITE_WITH `(x - u0) dot (u1 - u0) = (x - u0 - v:real^3) dot (u1 - u0) + ((v+u0) - u0) dot (u1 - u0)`);
 (VECTOR_ARITH_TAC);
 (ASM_REWRITE_TAC[REAL_ADD_LID]);
 (EXPAND_TAC "a");
 (REWRITE_TAC[ASSUME `hl [u0;u1:real^3] = inv (&2) * dist (u1,u0)`]);

 (REWRITE_WITH `((v + u0) - u0) dot (u1 - u0) = dist (u0 + v, u0) * dist (u1,u0:real^3)`);
 (REWRITE_WITH `v + u0 = u0 + v:real^3`);
 (VECTOR_ARITH_TAC);
 (REWRITE_TAC[dist; NORM_CAUCHY_SCHWARZ_EQ]);
 (UNDISCH_TAC `between (v + u0) (u0,u1:real^3)` THEN  
   REWRITE_TAC[BETWEEN_IN_CONVEX_HULL; CONVEX_HULL_2; IN; 
   IN_ELIM_THM] THEN STRIP_TAC);
 (UP_ASM_TAC THEN REWRITE_WITH `v + u0 = u0 + v:real^3`);
 (VECTOR_ARITH_TAC);
 (STRIP_TAC);
 (ASM_REWRITE_TAC[]);
 (REWRITE_WITH `(u % u0 + v' % u1) - u0:real^3 = (u % u0 + v' % u1) - (u + v') % u0`);
 (ASM_REWRITE_TAC[VECTOR_MUL_LID]);
 (REWRITE_WITH `(u % u0 + v' % u1) - (u + v') % u0 = v' % (u1 - u0:real^3)`);
 (VECTOR_ARITH_TAC);
 (REWRITE_TAC[NORM_MUL; VECTOR_MUL_ASSOC]);
 (AP_THM_TAC THEN AP_TERM_TAC);
 (REWRITE_WITH `abs v' = v'`);
 (ASM_REWRITE_TAC[REAL_ABS_REFL]);
 (REAL_ARITH_TAC);
 (REWRITE_TAC[REAL_ARITH `a * b >= c * b * d <=> &0 <= (a - c* d) * b`]);
 (REWRITE_WITH `&0 <= (dist (u0 + v,u0) - dist (x,u0) * (inv (&2) * dist (u1,u0)) /   sqrt (&2)) *  dist (u1,u0) <=>  &0 <= dist (u0 + v,u0:real^3) -  dist (x,u0) *   (inv (&2) * dist (u1,u0:real^3)) / sqrt (&2)`);
 (MP_TAC REAL_LE_MUL_EQ);
 (STRIP_TAC);
 (FIRST_ASSUM MATCH_MP_TAC);
 (MATCH_MP_TAC DIST_POS_LT);
 (ASM_MESON_TAC[]);
 (STRIP_TAC);

 (NEW_GOAL `&0 <= dist (x,u0:real^3) * (inv (&2) * dist (u1,u0)) / sqrt (&2) - dist (u1,u0) / sqrt (&2)`);
 (REWRITE_TAC[REAL_ARITH `a * (b * c)/d - c/d = (a * b - &1) * (c/d)`]);
 (MATCH_MP_TAC REAL_LE_MUL);
 (ASM_SIMP_TAC[REAL_LE_DIV; DIST_POS_LE; SQRT_POS_LE; REAL_ARITH `&0 <= &2`]);
 (NEW_GOAL `&2 <= dist (x,u0:real^3)`);
 (UNDISCH_TAC `packing (V:real^3->bool)` );
 (REWRITE_TAC[packing]);
 (STRIP_TAC THEN FIRST_ASSUM MATCH_MP_TAC);
 (REPEAT STRIP_TAC);
 (ONCE_REWRITE_TAC[GSYM IN]);
 (ASM_REWRITE_TAC[]);
 (ONCE_REWRITE_TAC[GSYM IN]);
 (MATCH_MP_TAC (SET_RULE `(?b. a IN b /\ b SUBSET c) ==> a IN c`));
 (EXISTS_TAC `set_of_list (ul:(real^3)list)`);
 (STRIP_TAC);
 (ASM_REWRITE_TAC[set_of_list] THEN SET_TAC[]);
 (MATCH_MP_TAC Packing3.BARV_SUBSET);
 (EXISTS_TAC `3` THEN ASM_REWRITE_TAC[]);
 (ASM_SET_TAC[]);
 (NEW_GOAL `&0 <= dist (x,u0:real^3) * inv (&2) - &2 * inv (&2)`);
 (REWRITE_TAC[REAL_ARITH `a * x - b * x = (a - b) * x`]);
 (MATCH_MP_TAC REAL_LE_MUL);
 (STRIP_TAC);
 (ASM_REAL_ARITH_TAC);
 (ASM_SIMP_TAC[REAL_LE_INV; REAL_ARITH `&0 <= &2`]);
 (NEW_GOAL `&2 * inv (&2) = &1`);
 (ASM_SIMP_TAC [REAL_ARITH `~(&2 = &0)`; REAL_MUL_RINV]);
 (ASM_REAL_ARITH_TAC);
 (NEW_GOAL `dist (u1,u0) / sqrt (&2) <= dist (u0 + v,u0:real^3)`);
 (ASM_REAL_ARITH_TAC);
 (NEW_GOAL `sqrt (&2) <= dist (u1,u0:real^3) / sqrt (&2)`);
 (REWRITE_WITH `sqrt (&2) <= dist (u1,u0:real^3) / sqrt (&2) 
  <=> sqrt (&2) * sqrt (&2) <= dist (u1,u0:real^3)`);
 (MATCH_MP_TAC REAL_LE_RDIV_EQ);
 (ASM_SIMP_TAC[SQRT_POS_LT; REAL_ARITH `&0 < &2`]);
 (SIMP_TAC[GSYM REAL_POW_2; SQRT_POW_2; REAL_ARITH `&0 <= &2`]);
 (UNDISCH_TAC `packing (V:real^3->bool)` );
 (REWRITE_TAC[packing]);
 (STRIP_TAC THEN FIRST_ASSUM MATCH_MP_TAC);
 (NEW_GOAL `{u0,u1:real^3} SUBSET V`);
 (MATCH_MP_TAC (SET_RULE `(?b. a SUBSET b /\ b SUBSET c) ==> a SUBSET c`));
 (EXISTS_TAC `set_of_list (ul:(real^3)list)`);
 (STRIP_TAC);
 (ASM_REWRITE_TAC[set_of_list] THEN SET_TAC[]);
 (MATCH_MP_TAC Packing3.BARV_SUBSET);
 (EXISTS_TAC `3` THEN ASM_REWRITE_TAC[]);
 (ASM_SET_TAC[]);
 (ASM_REAL_ARITH_TAC);

(* continue the next important lemma, similarly as above *)

 (NEW_GOAL `sqrt (&2) <= dist (u0 + v,u1:real^3)`);
 (UNDISCH_TAC `x:real^3 IN rcone_ge u1 u0 a`);
 (REWRITE_TAC[rcone_ge; rconesgn; IN; IN_ELIM_THM]);
 (REWRITE_WITH `(x - u1) dot (u0 - u1) = --((x - u0 - v:real^3) dot (u1 - u0)) + ((v+u0) - u1) dot (u0 - u1)`);
 (VECTOR_ARITH_TAC);
 (ASM_REWRITE_TAC[REAL_ADD_LID;REAL_NEG_0]);
 (EXPAND_TAC "a");
 (REWRITE_TAC[ASSUME `hl [u0;u1:real^3] = inv (&2) * dist (u1,u0)`]);

 (REWRITE_WITH `((v + u0) - u1) dot (u0 - u1) = dist (u0 + v, u1) * dist (u1,u0:real^3)`);
 (REWRITE_WITH `v + u0 = u0 + v:real^3 /\ dist (u1,u0) = dist (u0,u1)`);
 (REWRITE_TAC[DIST_SYM]);
 (VECTOR_ARITH_TAC);
 (REWRITE_TAC[dist; NORM_CAUCHY_SCHWARZ_EQ]);
 (UNDISCH_TAC `between (v + u0) (u0,u1:real^3)` THEN  
   REWRITE_TAC[BETWEEN_IN_CONVEX_HULL; CONVEX_HULL_2; IN; 
   IN_ELIM_THM] THEN STRIP_TAC);
 (UP_ASM_TAC THEN REWRITE_WITH `v + u0 = u0 + v:real^3`);
 (VECTOR_ARITH_TAC);
 (STRIP_TAC);
 (ASM_REWRITE_TAC[]);

 (REWRITE_WITH `(u % u0 + v' % u1) - u1:real^3 = (u % u0 + v' % u1) - (u + v') % u1`);
 (ASM_REWRITE_TAC[VECTOR_MUL_LID]);
 (REWRITE_WITH `(u % u0 + v' % u1) - (u + v') % u1 = u % (u0 - u1:real^3)`);
 (VECTOR_ARITH_TAC);
 (REWRITE_TAC[NORM_MUL; VECTOR_MUL_ASSOC]);
 (AP_THM_TAC THEN AP_TERM_TAC);
 (REWRITE_WITH `abs u = u`);
 (ASM_REWRITE_TAC[REAL_ABS_REFL]);
 (REAL_ARITH_TAC);
 (REWRITE_WITH `dist (u0,u1:real^3) = dist (u1,u0)`);
 (REWRITE_TAC[DIST_SYM]);
 (REWRITE_TAC[REAL_ARITH `a * b >= c * b * d <=> &0 <= (a - c* d) * b`]);
 (REWRITE_WITH `&0 <= (dist (u0 + v,u1) - dist (x,u1) * (inv (&2) * dist (u1,u0)) /   sqrt (&2)) *  dist (u1,u0) <=>  &0 <= dist (u0 + v,u1:real^3) -  dist (x,u1) *   (inv (&2) * dist (u1,u0:real^3)) / sqrt (&2)`);
 (MP_TAC REAL_LE_MUL_EQ);
 (STRIP_TAC);
 (FIRST_ASSUM MATCH_MP_TAC);
 (MATCH_MP_TAC DIST_POS_LT);
 (ASM_MESON_TAC[]);
 (STRIP_TAC);

 (NEW_GOAL `&0 <= dist (x,u1:real^3) * (inv (&2) * dist (u1,u0)) / sqrt (&2) - dist (u1,u0) / sqrt (&2)`);
 (REWRITE_TAC[REAL_ARITH `a * (b * c)/d - c/d = (a * b - &1) * (c/d)`]);
 (MATCH_MP_TAC REAL_LE_MUL);
 (ASM_SIMP_TAC[REAL_LE_DIV; DIST_POS_LE; SQRT_POS_LE; REAL_ARITH `&0 <= &2`]);
 (NEW_GOAL `&2 <= dist (x,u1:real^3)`);
 (UNDISCH_TAC `packing (V:real^3->bool)` );
 (REWRITE_TAC[packing]);
 (STRIP_TAC THEN FIRST_ASSUM MATCH_MP_TAC);
 (REPEAT STRIP_TAC);
 (ONCE_REWRITE_TAC[GSYM IN]);
 (ASM_REWRITE_TAC[]);
 (ONCE_REWRITE_TAC[GSYM IN]);
 (MATCH_MP_TAC (SET_RULE `(?b. a IN b /\ b SUBSET c) ==> a IN c`));
 (EXISTS_TAC `set_of_list (ul:(real^3)list)`);
 (STRIP_TAC);
 (ASM_REWRITE_TAC[set_of_list] THEN SET_TAC[]);
 (MATCH_MP_TAC Packing3.BARV_SUBSET);
 (EXISTS_TAC `3` THEN ASM_REWRITE_TAC[]);
 (ASM_SET_TAC[]);
 (NEW_GOAL `&0 <= dist (x,u1:real^3) * inv (&2) - &2 * inv (&2)`);
 (REWRITE_TAC[REAL_ARITH `a * x - b * x = (a - b) * x`]);
 (MATCH_MP_TAC REAL_LE_MUL);
 (STRIP_TAC);
 (ASM_REAL_ARITH_TAC);
 (ASM_SIMP_TAC[REAL_LE_INV; REAL_ARITH `&0 <= &2`]);
 (NEW_GOAL `&2 * inv (&2) = &1`);
 (ASM_SIMP_TAC [REAL_ARITH `~(&2 = &0)`; REAL_MUL_RINV]);
 (ASM_REAL_ARITH_TAC);
 (NEW_GOAL `dist (u1,u0) / sqrt (&2) <= dist (u0 + v,u1:real^3)`);
 (ASM_REAL_ARITH_TAC);
 (NEW_GOAL `sqrt (&2) <= dist (u1,u0:real^3) / sqrt (&2)`);
 (REWRITE_WITH `sqrt (&2) <= dist (u1,u0:real^3) / sqrt (&2) 
  <=> sqrt (&2) * sqrt (&2) <= dist (u1,u0:real^3)`);
 (MATCH_MP_TAC REAL_LE_RDIV_EQ);
 (ASM_SIMP_TAC[SQRT_POS_LT; REAL_ARITH `&0 < &2`]);
 (SIMP_TAC[GSYM REAL_POW_2; SQRT_POW_2; REAL_ARITH `&0 <= &2`]);
 (UNDISCH_TAC `packing (V:real^3->bool)` );
 (REWRITE_TAC[packing]);
 (STRIP_TAC THEN FIRST_ASSUM MATCH_MP_TAC);
 (NEW_GOAL `{u0,u1:real^3} SUBSET V`);
 (MATCH_MP_TAC (SET_RULE `(?b. a SUBSET b /\ b SUBSET c) ==> a SUBSET c`));
 (EXISTS_TAC `set_of_list (ul:(real^3)list)`);
 (STRIP_TAC);
 (ASM_REWRITE_TAC[set_of_list] THEN SET_TAC[]);
 (MATCH_MP_TAC Packing3.BARV_SUBSET);
 (EXISTS_TAC `3` THEN ASM_REWRITE_TAC[]);
 (ASM_SET_TAC[]);
 (ASM_REAL_ARITH_TAC);
 (UNDISCH_TAC `between (v + u0) (u0,u1:real^3)`);
 (REWRITE_WITH `v + u0 = u0 + v:real^3`);
 (VECTOR_ARITH_TAC);
 (REWRITE_TAC[between] THEN STRIP_TAC);
 (NEW_GOAL `dist (u0,u1:real^3) = dist (u0 + v,u0) + dist (u0 + v,u1)`);
 (ASM_MESON_TAC[DIST_SYM]);
 (NEW_GOAL `&2 * sqrt (&2) <= dist (u0,u1:real^3)`);
 (NEW_GOAL `&2 * sqrt (&2) <= dist (u0 + v,u0) + dist (u0 + v,u1:real^3)`);
 (ASM_REAL_ARITH_TAC);
 (ASM_REAL_ARITH_TAC);
 (UP_ASM_TAC THEN REWRITE_WITH `dist (u0,u1) = &2 * hl [u0; u1:real^3]`);
 (REWRITE_TAC[ASSUME `hl [u0; u1] = inv (&2) * dist (u1,u0:real^3)`]);
 (REWRITE_TAC[REAL_ARITH `a * b * c = (a * b) * c`]);
 (REWRITE_WITH `&2 * inv (&2) = &1`);
 (ASM_SIMP_TAC[REAL_MUL_RINV; REAL_ARITH `~(&2 = &0)`]);
 (REWRITE_TAC[DIST_SYM;REAL_MUL_LID]);
 (REWRITE_TAC[REAL_ARITH `~(a * b <= a * c) <=> &0 < a * (b - c)`]);
 (MATCH_MP_TAC REAL_LT_MUL);
 (ASM_REAL_ARITH_TAC);
 (ASM_MESON_TAC[]);
(* Finish the most difficult part in case k = 2 *)

 (REWRITE_TAC[IN_INTER]);
 (STRIP_TAC);
 (NEW_GOAL `{u0,u1:real^3} SUBSET V`);
 (MATCH_MP_TAC (SET_RULE `(?b. a SUBSET b /\ b SUBSET c) ==> a SUBSET c`));
 (EXISTS_TAC `set_of_list (ul:(real^3)list)`);
 (STRIP_TAC);
 (ASM_REWRITE_TAC[set_of_list] THEN SET_TAC[]);
 (MATCH_MP_TAC Packing3.BARV_SUBSET);
 (EXISTS_TAC `3` THEN ASM_REWRITE_TAC[]);
 (ASM_SET_TAC[]);
 (REPEAT STRIP_TAC);
 (REWRITE_TAC[rcone_ge; rconesgn; IN; IN_ELIM_THM]);
 (ASM_CASES_TAC `x = u0:real^3`);
 (ASM_REWRITE_TAC[VECTOR_ARITH `(a - a) dot b = &0`; DIST_REFL; REAL_MUL_LZERO; REAL_ARITH `&0 >= &0`]);
 (ASM_CASES_TAC `x = u1:real^3`);
 (ASM_REWRITE_TAC[GSYM NORM_POW_2; GSYM dist; REAL_POW_2;
REAL_ARITH `a * a >= a * a * c <=> &0 <= a * a * (&1 - c)`]);
 (MATCH_MP_TAC REAL_LE_MUL);
 (ASM_SIMP_TAC[DIST_POS_LE]);
 (MATCH_MP_TAC REAL_LE_MUL);
 (ASM_SIMP_TAC[DIST_POS_LE]);
 (ASM_REAL_ARITH_TAC);
 (NEW_GOAL `F`);
 (ASM_SET_TAC[]);
 (ASM_MESON_TAC[]);

 (REWRITE_TAC[rcone_ge; rconesgn; IN; IN_ELIM_THM]);
 (ASM_CASES_TAC `x = u1:real^3`);
 (ASM_REWRITE_TAC[VECTOR_ARITH `(a - a) dot b = &0`; DIST_REFL; REAL_MUL_LZERO; REAL_ARITH `&0 >= &0`]);
 (ASM_CASES_TAC `x = u0:real^3`);
 (ASM_REWRITE_TAC[GSYM NORM_POW_2; GSYM dist; REAL_POW_2;
REAL_ARITH `a * a >= a * a * c <=> &0 <= a * a * (&1 - c)`]);
 (MATCH_MP_TAC REAL_LE_MUL);
 (ASM_SIMP_TAC[DIST_POS_LE]);
 (MATCH_MP_TAC REAL_LE_MUL);
 (ASM_SIMP_TAC[DIST_POS_LE]);
 (ASM_REAL_ARITH_TAC);
 (NEW_GOAL `F`);
 (ASM_SET_TAC[]);
 (ASM_MESON_TAC[]);

 (ABBREV_TAC `m = mxi V [u0;u1;u2;u3:real^3]`);
 (ABBREV_TAC `s = omega_list_n V [u0;u1;u2;u3:real^3] 3`);
 (REWRITE_TAC[aff_ge_def; IN_AFFSIGN; sgn_ge]);

 (ABBREV_TAC `f = (\t:real^3. if t = x then &1 else &0)`);
 (EXISTS_TAC `f:real^3->real`);
 (REPEAT STRIP_TAC);
 (EXPAND_TAC "f");
 (COND_CASES_TAC);
 (REAL_ARITH_TAC);
 (REAL_ARITH_TAC);
 (REWRITE_WITH `{u0, u1} UNION {m, s} = {m,s} UNION {u0,u1:real^3}`);
 (SET_TAC[]);
 (REWRITE_WITH `sum ({m, s} UNION {u0, u1}) f = sum ({u0, u1:real^3}) f`);
 (MATCH_MP_TAC SUM_UNION_LZERO);
 (REWRITE_TAC[Geomdetail.FINITE6]);
 (REPEAT STRIP_TAC THEN EXPAND_TAC "f");
 (COND_CASES_TAC);
 (NEW_GOAL `F`);
 (ASM_SET_TAC[]);
 (ASM_MESON_TAC[]);
 (MESON_TAC[]);
 (ASM_SIMP_TAC[Collect_geom.SUM_DIS2]);
 (EXPAND_TAC "f");
 (COND_CASES_TAC);
 (COND_CASES_TAC);
 (NEW_GOAL `F`);
 (ASM_MESON_TAC[]);
 (ASM_MESON_TAC[]);
 (REAL_ARITH_TAC);
 (COND_CASES_TAC);
 (REAL_ARITH_TAC);
 (NEW_GOAL `F`);
 (ASM_SET_TAC[]);
 (ASM_MESON_TAC[]);

 (REWRITE_WITH `{u0, u1} UNION {m, s} = {m,s} UNION {u0,u1:real^3}`);
 (SET_TAC[]);
 (REWRITE_WITH `vsum ({m:real^3, s} UNION {u0, u1}) (\x. f x % x) 
 = vsum ({u0,u1}) (\x. f x % x)`);
 (MATCH_MP_TAC VSUM_UNION_LZERO);
 (REWRITE_TAC[Geomdetail.FINITE6]);
 (REPEAT STRIP_TAC THEN EXPAND_TAC "f");
 (COND_CASES_TAC);
 (NEW_GOAL `F`);
 (ASM_SET_TAC[]);
 (ASM_MESON_TAC[]);
 (VECTOR_ARITH_TAC);
 (ASM_SIMP_TAC[Collect_geom.VSUM_DIS2]);
 (EXPAND_TAC "f");
 (COND_CASES_TAC);
 (COND_CASES_TAC);
 (NEW_GOAL `F`);
 (ASM_MESON_TAC[]);
 (ASM_MESON_TAC[]);
 (ASM_REWRITE_TAC[] THEN VECTOR_ARITH_TAC);
 (COND_CASES_TAC);
 (ASM_REWRITE_TAC[] THEN VECTOR_ARITH_TAC);
 (NEW_GOAL `F`);
 (ASM_SET_TAC[]);
 (ASM_MESON_TAC[]);

 (NEW_GOAL `F`);
 (UNDISCH_TAC `~(mcell k V ul = {})`);
 (ASM_REWRITE_TAC[MCELL_EXPLICIT;mcell2]);
 (LET_TAC);
 (COND_CASES_TAC);
 (NEW_GOAL `F`);
 (UP_ASM_TAC THEN ASM_REWRITE_TAC[TRUNCATE_SIMPLEX_EXPLICIT_1]);
 (ASM_MESON_TAC[]);
 (REWRITE_TAC[]);
 (ASM_MESON_TAC[]);

(* ========================================================================= *)
(*  -------------------- k = 4  --------------------------------------------*)

 (ASM_CASES_TAC `k = 4`);
 (ASM_REWRITE_TAC[TRUNCATE_SIMPLEX_EXPLICIT_3; ARITH_RULE `4 - 1 = 3`;
   set_of_list]);
 (ASM_SIMP_TAC[MCELL_EXPLICIT; ARITH_RULE `4 >= 4`]);
 (REWRITE_TAC[SET_EQ_LEMMA; IN_INTER]);
 (REWRITE_TAC[mcell4]);
 (COND_CASES_TAC);
 (REWRITE_TAC[set_of_list]);
 (REPEAT STRIP_TAC);

 (ABBREV_TAC `S = {u0, u1,u2, u3:real^3}`);
 (ABBREV_TAC `s3 = omega_list V (ul:(real^3)list)`);

 (ASM_CASES_TAC `(x:real^3) IN S`);
 (ASM_MESON_TAC[]);
 (NEW_GOAL `!x:real^3. x IN convex hull S /\ ~(x IN S)
                  ==> (?y. y IN S /\ norm (x - s3) < norm (y - s3))`);
 (MATCH_MP_TAC SIMPLEX_FURTHEST_LT_2);
 (EXPAND_TAC "S");
 (REWRITE_TAC[Geomdetail.FINITE6]);
 (NEW_GOAL `?y:real^3. y IN S /\ norm (x - s3) < norm (y - s3)`);
 (FIRST_ASSUM MATCH_MP_TAC);
 (ASM_REWRITE_TAC[IN]);
 (UP_ASM_TAC THEN EXPAND_TAC "S" THEN STRIP_TAC);
 (NEW_GOAL `s3 IN (voronoi_closed V (y:real^3))`);
 (EXPAND_TAC "s3");
 (MATCH_MP_TAC (SET_RULE `(?b. a IN b /\ b SUBSET c) ==> a IN c`));
 (EXISTS_TAC `voronoi_list V (ul:(real^3)list)`);
 (STRIP_TAC);
 (MATCH_MP_TAC Packing3.OMEGA_LIST_IN_VORONOI_LIST);
 (EXISTS_TAC `3` THEN ASM_REWRITE_TAC[]);
 (ASM_REWRITE_TAC[VORONOI_LIST; VORONOI_SET; set_of_list]);
 (NEW_GOAL `y:real^3 IN S`);
 (ASM_MESON_TAC[]);
 (ASM_SET_TAC[]);
 (UP_ASM_TAC THEN REWRITE_TAC[voronoi_closed; IN; IN_ELIM_THM]);
 (ONCE_REWRITE_TAC[DIST_SYM]);
 (REWRITE_TAC[dist]);
 (STRIP_TAC);
 (NEW_GOAL `F`);
 (NEW_GOAL `norm (y - s3) <= norm (x - s3:real^3)`);
 (FIRST_ASSUM MATCH_MP_TAC);
 (UNDISCH_TAC `x:real^3 IN V` THEN MESON_TAC[IN]);
 (ASM_REAL_ARITH_TAC);
 (ASM_MESON_TAC[]);
 (MATCH_MP_TAC (SET_RULE `(?b. a IN b /\ b SUBSET c ) ==> a IN c`));
 (EXISTS_TAC `set_of_list (ul:(real^3)list)`);
 (STRIP_TAC);
 (ASM_REWRITE_TAC[set_of_list]);
 (MATCH_MP_TAC Packing3.BARV_SUBSET);
 (EXISTS_TAC `3` THEN ASM_REWRITE_TAC[]);
 (UP_ASM_TAC THEN REWRITE_TAC[IN_SET_IMP_IN_CONVEX_HULL_SET]);
 (NEW_GOAL `F`);
 (UNDISCH_TAC `~(mcell k V ul = {})`);
 (ASM_SIMP_TAC[MCELL_EXPLICIT; ARITH_RULE `4 >= 4`]);
 (REWRITE_TAC[mcell4]);
 (COND_CASES_TAC);
 (NEW_GOAL `F`);
 (ASM_MESON_TAC[]);
 (ASM_MESON_TAC[]);
 (REFL_TAC);
 (ASM_MESON_TAC[]);

(* ========================================================================= *)
(*  -------------------- k = 3  --------------------------------------------*)


 (ASM_CASES_TAC `k = 3`);
 (NEW_GOAL `{u0,u1,u2,u3:real^3} SUBSET V`);
 (REWRITE_WITH `{u0,u1,u2,u3:real^3} = set_of_list ul`);
 (ASM_REWRITE_TAC[set_of_list]);
 (MATCH_MP_TAC Packing3.BARV_SUBSET);
 (EXISTS_TAC `3` THEN ASM_REWRITE_TAC[]);
 (ASM_REWRITE_TAC[ARITH_RULE `3 - 1 = 2`; TRUNCATE_SIMPLEX_EXPLICIT_2;
   MCELL_EXPLICIT; mcell3; set_of_list]);
 (ABBREV_TAC `m = mxi V [u0; u1; u2; u3:real^3]`);
 (REWRITE_TAC[SET_RULE `{a,b,c} UNION {d} = {a,b,c,d}`]);
 (COND_CASES_TAC);

 (REWRITE_TAC[set_of_list; SET_EQ_LEMMA; IN_INTER]);
 (REPEAT STRIP_TAC);
 (ASM_CASES_TAC `x IN {u0,u1,u2:real^3}`);
 (ASM_REWRITE_TAC[]);
 (NEW_GOAL `F`);
 (ABBREV_TAC `S = {u0, u1,u2, m:real^3}`);
 (ABBREV_TAC `s2 = omega_list_n V (ul:(real^3)list) 2`);
 (ABBREV_TAC `s3 = omega_list V (ul:(real^3)list)`);

 (NEW_GOAL `s2 IN voronoi_list V [u0;u1;u2:real^3]`);
 (EXPAND_TAC "s2");
 (REWRITE_TAC[OMEGA_LIST_TRUNCATE_2; ASSUME `ul = [u0; u1; u2; u3:real^3]`]);
 (MATCH_MP_TAC Packing3.OMEGA_LIST_IN_VORONOI_LIST);
 (EXISTS_TAC `2`);
 (REWRITE_WITH `[u0;u1;u2:real^3] = truncate_simplex 2 ul`);
 (ASM_REWRITE_TAC[TRUNCATE_SIMPLEX_EXPLICIT_2]);
 (MATCH_MP_TAC Packing3.TRUNCATE_SIMPLEX_BARV);
 (EXISTS_TAC `3` THEN ASM_REWRITE_TAC[ARITH_RULE `2 <= 3`]);

 (NEW_GOAL `s3:real^3 IN voronoi_list V ul`);
 (EXPAND_TAC "s3");
 (MATCH_MP_TAC Packing3.OMEGA_LIST_IN_VORONOI_LIST);
 (EXISTS_TAC `3` THEN ASM_REWRITE_TAC[]);

 (NEW_GOAL `s3:real^3 IN voronoi_list V [u0;u1;u2:real^3]`);
 (NEW_GOAL `voronoi_list V ul SUBSET voronoi_list V [u0;u1;u2:real^3]`);
 (ASM_REWRITE_TAC[VORONOI_LIST; VORONOI_SET; set_of_list]);
 (SET_TAC[]);
 (ASM_SET_TAC[]);

 (NEW_GOAL `?s. between s (omega_list_n V ul 2,omega_list_n V ul 3) /\
                  dist (u0,s) = sqrt (&2) /\
                  mxi V ul = s`);
 (MATCH_MP_TAC MXI_EXPLICIT);
 (EXISTS_TAC `u1:real^3` THEN EXISTS_TAC `u2:real^3` THEN EXISTS_TAC `u3:real^3`);
 (ASM_REWRITE_TAC[TRUNCATE_SIMPLEX_EXPLICIT_2]);
 (UP_ASM_TAC THEN ASM_REWRITE_TAC[] THEN STRIP_TAC);
 (NEW_GOAL `omega_list_n V [u0; u1; u2; u3] 3 = s3:real^3`);
 (EXPAND_TAC "s3" THEN ONCE_REWRITE_TAC[EQ_SYM_EQ]);
 (REWRITE_WITH `ul = truncate_simplex 3 [u0; u1; u2; u3:real^3]`);
 (ASM_REWRITE_TAC[TRUNCATE_SIMPLEX_EXPLICIT_3]);
 (MATCH_MP_TAC Packing3.OMEGA_LIST_LEMMA);
 (REWRITE_TAC[LENGTH] THEN ARITH_TAC);
 (MP_TAC (ASSUME `between s (s2,omega_list_n V [u0; u1; u2; u3:real^3] 3)`));
 (REWRITE_TAC[ASSUME `omega_list_n V [u0; u1; u2; u3] 3 = s3:real^3`]);
 (REWRITE_TAC[BETWEEN_IN_CONVEX_HULL]);
 (STRIP_TAC);
 (NEW_GOAL `convex hull {s2,s3:real^3} SUBSET voronoi_list V [u0; u1; u2]`);
 (REWRITE_WITH `convex hull {s2, s3} SUBSET voronoi_list V [u0; u1; u2:real^3]    <=> convex hull {s2, s3} SUBSET convex hull (voronoi_list V [u0; u1; u2])`);
 (ONCE_REWRITE_TAC[EQ_SYM_EQ]);
 (REWRITE_WITH `convex hull (voronoi_list V [u0; u1; u2]) = voronoi_list V [u0; u1; u2:real^3]`);
 (REWRITE_TAC[CONVEX_HULL_EQ; Packing3.CONVEX_VORONOI_LIST]);
 (MATCH_MP_TAC Marchal_cells.CONVEX_HULL_SUBSET);
 (ASM_SET_TAC[]);
 (NEW_GOAL `s IN voronoi_list V [u0;u1;u2:real^3]`);
 (ASM_SET_TAC[]);

 (NEW_GOAL `~(m:real^3 IN V)`);
 (UP_ASM_TAC THEN ASM_REWRITE_TAC[VORONOI_LIST;VORONOI_SET;set_of_list]);
 (REPEAT STRIP_TAC);
 (NEW_GOAL `s IN voronoi_closed V (u0:real^3)`);
 (ASM_SET_TAC[]);
 (UP_ASM_TAC THEN REWRITE_TAC[voronoi_closed; IN; IN_ELIM_THM]);
 (STRIP_TAC);
 (NEW_GOAL `dist (s, u0:real^3) <= dist (s, s)`);
 (FIRST_ASSUM MATCH_MP_TAC);
 (ASM_SET_TAC[]);
 (UP_ASM_TAC THEN REWRITE_TAC[DIST_REFL] THEN ONCE_REWRITE_TAC[DIST_SYM] 
   THEN ASM_REWRITE_TAC[REAL_ARITH `~(a <= b) <=> b < a`]);
 (ASM_SIMP_TAC[SQRT_POS_LT; REAL_ARITH `&0 < &2`]);

 (NEW_GOAL `!x:real^3. x IN convex hull S /\ ~(x IN S)
                  ==> (?y. y IN S /\ norm (x - m) < norm (y - m))`);
 (MATCH_MP_TAC SIMPLEX_FURTHEST_LT_2);
 (EXPAND_TAC "S");
 (REWRITE_TAC[Geomdetail.FINITE6]);
 (NEW_GOAL `?y:real^3. y IN S /\ norm (x - m) < norm (y - m)`);
 (FIRST_ASSUM MATCH_MP_TAC);
 (ASM_REWRITE_TAC[]);
 (EXPAND_TAC "S");
 (STRIP_TAC);
 (NEW_GOAL `x = m:real^3`);
 (ASM_SET_TAC[]);
 (ASM_MESON_TAC[]);
 (UP_ASM_TAC THEN EXPAND_TAC "S" THEN STRIP_TAC);

 (NEW_GOAL `norm (u0 - m:real^3) = sqrt (&2)`);
 (ASM_REWRITE_TAC[GSYM dist]);
 (UNDISCH_TAC `s IN voronoi_list V [u0;u1;u2:real^3]`);
 (ASM_REWRITE_TAC[VORONOI_LIST;VORONOI_SET;set_of_list]);
 (REPEAT STRIP_TAC);

 (NEW_GOAL `s IN voronoi_closed V (u0:real^3)`);
 (ASM_SET_TAC[]);
 (NEW_GOAL `s IN voronoi_closed V (u1:real^3)`);
 (ASM_SET_TAC[]);
 (NEW_GOAL `s IN voronoi_closed V (u2:real^3)`);
 (ASM_SET_TAC[]);
 (NEW_GOAL `norm (u0 - m:real^3) <= norm (u1 - m)`);
 (REWRITE_TAC[GSYM dist] THEN ONCE_REWRITE_TAC[DIST_SYM]);
 (ASM_REWRITE_TAC[]);
 (UNDISCH_TAC `s IN voronoi_closed V (u0:real^3)` THEN 
  REWRITE_TAC[voronoi_closed; IN; IN_ELIM_THM]);
 (STRIP_TAC);
 (FIRST_ASSUM MATCH_MP_TAC);
 (ASM_SET_TAC[]);

 (NEW_GOAL `norm (u1 - m:real^3) <= norm (u2 - m)`);
 (REWRITE_TAC[GSYM dist] THEN ONCE_REWRITE_TAC[DIST_SYM]);
 (ASM_REWRITE_TAC[]);
 (UNDISCH_TAC `s IN voronoi_closed V (u1:real^3)` THEN 
  REWRITE_TAC[voronoi_closed; IN; IN_ELIM_THM]);
 (STRIP_TAC);
 (FIRST_ASSUM MATCH_MP_TAC);
 (ASM_SET_TAC[]);

 (NEW_GOAL `norm (u2 - m:real^3) <= norm (u0 - m)`);
 (REWRITE_TAC[GSYM dist] THEN ONCE_REWRITE_TAC[DIST_SYM]);
 (ASM_REWRITE_TAC[]);
 (UNDISCH_TAC `s IN voronoi_closed V (u2:real^3)` THEN 
  REWRITE_TAC[voronoi_closed; IN; IN_ELIM_THM]);
 (STRIP_TAC);
 (FIRST_ASSUM MATCH_MP_TAC);
 (ASM_SET_TAC[]);

 (NEW_GOAL `norm (m - m:real^3) = &0`);
 (NORM_ARITH_TAC);
 (NEW_GOAL `norm (y - m:real^3) <= sqrt (&2)`);
 (UNDISCH_TAC `y IN {u0,u1,u2,m:real^3}` THEN 
   REWRITE_TAC[SET_RULE `(y IN {u0,u1,u2,m:real^3}) 
  <=> (y = u0 \/ y = u1 \/ y = u2 \/ y = m)`]);
 (NEW_GOAL `&0 <= sqrt (&2)`);
 (SIMP_TAC[REAL_ARITH `&0 <= &2`; SQRT_POS_LE]);
 (NEW_GOAL `norm (u1 - m) = sqrt (&2) /\ norm (u2 - m:real^3) = sqrt (&2)`);
 (ASM_REAL_ARITH_TAC);
 (STRIP_TAC);
 (REWRITE_TAC[ASSUME `y = u0:real^3`; ASSUME `norm (u0 - m:real^3) = sqrt (&2)`] THEN REAL_ARITH_TAC);
 (REWRITE_TAC[ASSUME `y = u1:real^3`; ASSUME `norm (u1 - m) = sqrt (&2) /\ norm (u2 - m:real^3) = sqrt (&2)`] THEN REAL_ARITH_TAC);
 (REWRITE_TAC[ASSUME `y = u2:real^3`; ASSUME `norm (u1 - m) = sqrt (&2) /\ norm (u2 - m:real^3) = sqrt (&2)`] THEN REAL_ARITH_TAC);
 (REWRITE_TAC[ASSUME `y = m:real^3`; ASSUME `norm (m - m:real^3) = &0`; 
  ASSUME `&0 <= sqrt (&2)`] THEN REAL_ARITH_TAC);
 (NEW_GOAL `norm (x - m:real^3) < norm (u0 -m)`);
 (ASM_REAL_ARITH_TAC);
 (UP_ASM_TAC THEN REWRITE_TAC[ASSUME `m = s:real^3`; GSYM dist]);
 (STRIP_TAC);
 (UNDISCH_TAC `s:real^3 IN voronoi_closed V u0`);
 (REWRITE_TAC[voronoi_closed; IN; IN_ELIM_THM]);
 (ONCE_REWRITE_TAC[DIST_SYM]);
 (STRIP_TAC);
 (NEW_GOAL `dist (u0,s:real^3) <= dist (x,s)`);
 (FIRST_ASSUM MATCH_MP_TAC);
 (ASM_SET_TAC[]);
 (ASM_REAL_ARITH_TAC);
 (ASM_MESON_TAC[]);
 (ASM_SET_TAC[]);
 (MATCH_MP_TAC IN_SET_IMP_IN_CONVEX_HULL_SET);
 (ASM_SET_TAC[]);
 (NEW_GOAL `F`);
 (UNDISCH_TAC `~( mcell k V (ul:(real^3)list) = {})`);
 (ASM_REWRITE_TAC[ARITH_RULE `3 - 1 = 2`; TRUNCATE_SIMPLEX_EXPLICIT_2;
   MCELL_EXPLICIT; mcell3; set_of_list]);
 (ASM_MESON_TAC[]);
 (NEW_GOAL `F`);
 (NEW_GOAL `k = 1 \/ k = 2 \/ k = 3 \/ k = 4`);
 (UNDISCH_TAC `1 <= k` THEN UNDISCH_TAC `k <= 4`);
 (ARITH_TAC);
 (ASM_ARITH_TAC);
 (ASM_MESON_TAC[])]);;

(* ========================================================================== *)
(*                 Lemma for the case mcell 0                                 *)
(* =========================================================================  *)

let LEPJBDJ_0 = prove_by_refinement 
(LEPJBDJ_0_concl,
[(REPEAT STRIP_TAC);
 (REWRITE_TAC[MCELL_EXPLICIT; mcell0]);
 (REWRITE_TAC[SET_RULE `a INTER b = {} <=> (!v. v IN a INTER b ==> F)`]);      
 (GEN_TAC THEN REWRITE_TAC[IN_INTER; IN_DIFF] THEN REPEAT STRIP_TAC );
 (NEW_GOAL `v:real^3 = HD ul`);
 (MATCH_MP_TAC ROGERS_INTER_V_LEMMA);
 (EXISTS_TAC `V:real^3->bool` THEN ASM_REWRITE_TAC[]);
 (ASM_SET_TAC[]);
 (UNDISCH_TAC `~(v:real^3 IN ball (HD ul, sqrt (&2)))`);
 (ASM_REWRITE_TAC[ball; IN; IN_ELIM_THM; DIST_REFL]);
 (ASM_SIMP_TAC[SQRT_POS_LT; REAL_ARITH `&0 < &2`])]);;

end;;