(* ========================================================================== *)
(* FLYSPECK - BOOK FORMALIZATION                                              *)
(*                                                                            *)
(* Chapter: nonlinear inequalities                                            *)
(* Author:  Thomas Hales     *)
(* Date: 2010-10-17                                                         *)
(* ========================================================================== *)

(*
Definitions and Lemmas used in the generation of interval arithmetic
code from the HOL Light specifications.

*)

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

module Nonlinear_lemma = struct

let ineq = Sphere.ineq;;

let NONLIN = new_definition `NONLIN = T`;;

let sqrt_x1 = define `sqrt_x1 x1 x2 x3 x4 x5 x6 = sqrt x1`;;

let sqrt_x2 = define `sqrt_x2 x1 x2 x3 x4 x5 x6 = sqrt x2`;;

let sqrt_x3 = define `sqrt_x3 x1 x2 x3 x4 x5 x6 = sqrt x3`;;

let sqrt_x4 = define `sqrt_x4 x1 x2 x3 x4 x5 x6 = sqrt x4`;;

let sqrt_x5 = define `sqrt_x5 x1 x2 x3 x4 x5 x6 = sqrt x5`;;

let sqrt_x6 = define `sqrt_x6 x1 x2 x3 x4 x5 x6 = sqrt x6`;;

let halfbump_x = new_definition `halfbump_x x = bump (sqrt x / &2)`;;

let halfbump_x1 = new_definition `halfbump_x1 x1 x2 x3 x4 x5 x6 = halfbump_x x1`;;

let halfbump_x4 = new_definition `halfbump_x4 x1 x2 x3 x4 x5 x6 = halfbump_x x4`;;

let unit6 = define `unit6 x1 x2 x3 x4 x5 x6 = &1`;;

let proj_x1 = define `proj_x1 (x1:A) (x2:B) (x3:C) (x4:D) (x5:E) (x6:F) = x1`;;

let proj_x2 = define `proj_x2 (x1:A) (x2:B) (x3:C) (x4:D) (x5:E) (x6:F) = x2`;;

let proj_x3 = define `proj_x3 (x1:A) (x2:B) (x3:C) (x4:D) (x5:E) (x6:F) = x3`;;

let proj_x4 = define `proj_x4 (x1:A) (x2:B) (x3:C) (x4:D) (x5:E) (x6:F) = x4`;;

let proj_x5 = define `proj_x5 (x1:A) (x2:B) (x3:C) (x4:D) (x5:E) (x6:F) = x5`;;

let proj_x6 = define `proj_x6 (x1:A) (x2:B) (x3:C) (x4:D) (x5:E) (x6:F) = x6`;;

let promote = define `promote f (x1:A) (x2:B) (x3:C) (x4:D) (x5:E) (x6:F) = f x1`;;

let unit0 = define `unit0 = &1`;; 

let pow1 = new_definition `pow1 y = y pow 1`;;

let pow2 = new_definition `pow2 y = y pow 2`;;

let pow3 = new_definition `pow3 y = y pow 3`;;

let pow4 = new_definition `pow4 y = y pow 4`;;

let promote_pow2 = new_definition `promote_pow2 x1 (x2:A) (x3:B) (x4:C) (x5:D) (x6:E) = x1 pow 2`;;

let promote_pow3 = new_definition `promote_pow3 x1 (x2:A) (x3:B) (x4:C) (x5:D) (x6:E) = x1 pow 3`;;

let compose6 = new_definition `compose6 f p1 p2 p3 p4 p5 p6 
  (x1:real) (x2:real) (x3:real) (x4:real) (x5:real) (x6:real) =
  (f:real->real->real->real->real->real->real)
  (p1 x1 x2 x3 x4 x5 x6)
    (p2 x1 x2 x3 x4 x5 x6)
    (p3 x1 x2 x3 x4 x5 x6)
    (p4 x1 x2 x3 x4 x5 x6)
    (p5 x1 x2 x3 x4 x5 x6)
    (p6 x1 x2 x3 x4 x5 x6)`;;

let scale6 = new_definition `scale6 f
   (r:real)   (x1:real) (x2:real) (x3:real) (x4:real) (x5:real) (x6:real) = 
   (f x1 x2 x3 x4 x5 x6) * r`;;

let quadratic_root_plus_curry = 
  new_definition `quadratic_root_plus_curry a b c = quadratic_root_plus (a,b,c)`;;

let sq_pow2 = prove_by_refinement(
  `!a x. a pow 2 <= x ==> (sqrt x * sqrt x = x)`,
  (* {{{ proof *)
  [
  REWRITE_TAC[GSYM REAL_POW_2;SQRT_POW2];
  MESON_TAC[REAL_LE_TRANS;Collect_geom.REAL_LE_SQUARE_POW];
  ]);;
  (* }}} *)

let sqrt2_sqrt2 = prove_by_refinement(
  `sqrt2 * sqrt2 = &2`,
  (* {{{ proof *)
  [
  REWRITE_TAC[Sphere.sqrt2];
    MATCH_MP_TAC sq_pow2;
  EXISTS_TAC`&0`;
  REAL_ARITH_TAC;
  ]);;
  (* }}} *)

(* SOME DEFINITIONS *)

let vol3f_sqrt2_lmplus = new_definition 
  `vol3f_sqrt2_lmplus y1 y2 (y3:real) (y4:real) (y5:real) y6 = 
    (&2 * mm1 / pi) *
             (&2 * dih_y y1 y2 sqrt2 sqrt2 sqrt2 y6 +
              &2 * dih2_y y1 y2 sqrt2 sqrt2 sqrt2 y6 +
              &2 * dih6_y y1 y2 sqrt2 sqrt2 sqrt2 y6 +
              dih3_y y1 y2 sqrt2 sqrt2 sqrt2 y6 +
              dih4_y y1 y2 sqrt2 sqrt2 sqrt2 y6 +
              dih5_y y1 y2 sqrt2 sqrt2 sqrt2 y6 - &3 * pi) -
             (&8 * mm2 / pi) *
             (
              lfun (y2 / &2) * dih2_y y1 y2 sqrt2 sqrt2 sqrt2 y6 +
              lfun (y6 / &2) * dih6_y y1 y2 sqrt2 sqrt2 sqrt2 y6)`;;

let vol3f_x_sqrt2_lmplus = new_definition
  `vol3f_x_sqrt2_lmplus x1 x2 x3 x4 x5 x6 =
   vol3f_sqrt2_lmplus (sqrt x1) (sqrt x2) (sqrt x3) (sqrt x4) (sqrt x5)      (sqrt x6)`;;

let vol3f_x_lfun = new_definition 
 `vol3f_x_lfun x1 x2 (x3:real) (x4:real) (x5:real) x6 = vol3f (sqrt x1) (sqrt x2) (sqrt x6)  sqrt2 lfun `;;

let vol3_x_sqrt = new_definition
 `vol3_x_sqrt x1 x2 (x3:real) (x4:real) (x5:real) x6  = vol_y sqrt2 sqrt2 sqrt2 (sqrt x1) (sqrt x2) (sqrt x6) `;;

let monomial = new_definition `monomial n1 n2 n3 n4 n5 n6 y1 y2 y3 y4 y5 y6 = 
   (y1 pow n1) * (y2 pow n2) * (y3 pow n3) * (y4 pow n4) * (y5 pow n5) * (y6 pow n6)`;;

let arclength_x_234 = new_definition `arclength_x_234  (x1:real) (x2:real) (x3:real) (x4:real) (x5:real) (x6:real) = arclength (sqrt x2) (sqrt x3) (sqrt x4)`;;

let arclength_x_126 = new_definition `arclength_x_126  (x1:real) (x2:real) (x3:real) (x4:real) (x5:real) (x6:real) = arclength (sqrt x1) (sqrt x2) (sqrt x6)`;;


(* ========================================================================== *)
(*    BASIC LEMMAS                                                            *)
(* ========================================================================== *)

let strip_let_tm t = snd(dest_eq(concl(REDEPTH_CONV let_CONV t)));;

let strip_let t = REWRITE_RULE[REDEPTH_CONV let_CONV (concl t )] t;;

let tame_table_d_values = prove_by_refinement(
  `tame_table_d 2 1 = #0.103 /\ tame_table_d 1 2 = #0.2759 /\ tame_table_d 0 3 = #0.4488 /\ tame_table_d 4 1 = #0.6548  /\
    tame_table_d 6 0 = #0.7578 /\
    tame_table_d 3 1 = #0.3789 /\
    tame_table_d 2 2 = #0.5518 /\
   tame_table_d 1 3 = #0.7247 /\
   tame_table_d 0 4 = #0.8976 /\
   tame_table_d 5 0 = #0.4819 /\
   tame_table_d 4 1 = #0.6548 /\
   tame_table_d 3 2 = #0.8277 /\
   tame_table_d 2 3 = #1.0006
 `,
  (* {{{ proof *)
  [
    REWRITE_TAC[Sphere.tame_table_d;ARITH_RULE `2 + 2 * 1 > 3 /\ 1 + 2 * 2 > 3 /\ 0 + 2 * 3 > 3 /\ 4 + 2 *1 > 3 /\ 6 + 2 * 0 > 3 /\ (2 + 2 * 3 > 3) /\ (3 + 2 * 2 > 3) /\ (5 + 2 * 0 > 3) /\ (0 + 2 * 4 > 3) /\ (1 + 2 * 3 > 3) /\ (2 + 2 * 2 > 3) /\ (3 + 2 * 1 > 3)` ];
    REAL_ARITH_TAC;
  ]);;
  (* }}} *)

let unit0f = prove_by_refinement(
  `f x1 x2 x3 x4 x5 x6 * unit0 = f x1 x2 x3 x4 x5 x6`,
  (* {{{ proof *)
  [
  REWRITE_TAC[unit0] THEN REAL_ARITH_TAC
  ]);;
  (* }}} *)

let sqrt8_sqrt2 = prove_by_refinement(
  `sqrt8 = &2 * sqrt2`,
  (* {{{ proof *)
  [
  SIMP_TAC[Sphere.sqrt8;Sphere.sqrt2;SQRT_MUL;
	   REAL_ARITH `&8 = &2 pow 2 * &2 /\ &0 <= &2 /\ &0 <= &2 pow 2 /\ 
      abs(&2) = &2`;POW_2_SQRT_ABS];
  ]);;
  (* }}} *)

let sqrt2_sqrt8 = prove_by_refinement(
  `sqrt2 = #0.5 * sqrt8`,
  (* {{{ proof *)
  [
  REWRITE_TAC[sqrt8_sqrt2];
  REAL_ARITH_TAC;
  ]);;
  (* }}} *)

let SQRT_MUL_POW_2= prove_by_refinement(`!(a:real) b. (a>= &0) /\ (b>= &0) ==> sqrt((a*a)*b)= a* sqrt(b)`,
[
SIMP_TAC[SQRT_MUL;REAL_LE_SQUARE;real_ge];
MESON_TAC[sq_pow2;REAL_ARITH `&0 pow 2 = &0`];
]);;

(* sol0 = const1 * pi, repeated from TameGeneral.hl *)
let sol0_EQ_sol_y = prove(`sol0 = sol_y (&2) (&2) (&2) (&2) (&2) (&2)`,
    REWRITE_TAC[Sphere.sol0; Sphere.sol_y; Sphere.dih_y; Sphere.dih_x; Sphere.delta_x4; Sphere.delta_x] THEN
      CONV_TAC (DEPTH_CONV let_CONV) THEN
      CONV_TAC (REAL_RAT_REDUCE_CONV) THEN
      MP_TAC (SPEC `&1 / &3` Trigonometry2.acs_atn2) THEN
      CONV_TAC (REAL_RAT_REDUCE_CONV) THEN
      DISCH_THEN (fun th -> REWRITE_TAC[th]) THEN
      REWRITE_TAC [REAL_ARITH `&3 * (a - b) - c = (a + d) + (a + d) + (a + d) - c <=> --b = d`] THEN
      MP_TAC (SPECL [`sqrt (&8 / &9)`; `&1 / &3`] Trigonometry1.ATN2_RNEG) THEN
      CONV_TAC (REAL_RAT_REDUCE_CONV) THEN
      DISCH_THEN (fun th -> REWRITE_TAC[SYM th]) THEN
      SUBGOAL_THEN `sqrt (&2048) = &48 * sqrt (&8 / &9) /\ -- &16 = &48 * (-- &1 / &3)` ASSUME_TAC THENL
      [
	MP_TAC (SPECL [`&48`; `&8 / &9`] SQRT_MUL_POW_2) THEN
	  CONV_TAC (REAL_RAT_REDUCE_CONV);
	ALL_TAC
      ] THEN
      ASM_REWRITE_TAC[] THEN
      MATCH_MP_TAC (GSYM Trigonometry1.ATN2_LMUL_EQ) THEN
      REAL_ARITH_TAC);;

let sol0_over_pi_EQ_const1 = prove(`sol0 / pi = const1`,
   REWRITE_TAC[sol0_EQ_sol_y; Sphere.const1]);;

let sol0_const1 = prove_by_refinement(
  `sol0 = pi * const1`,
  (* {{{ proof *)
  [
  REWRITE_TAC[GSYM sol0_over_pi_EQ_const1];
  MP_TAC PI_POS;
  CONV_TAC REAL_FIELD;
  ]);;
  (* }}} *)

let ineq_lemma_b = prove_by_refinement(
  `!a y b. (&0 <= a  /\  &0 <= b /\  a <= y /\ y <= b)  ==> 
    a pow 2 <= y pow 2 /\ y pow 2 <= b pow 2 /\ (sqrt (y pow 2) = y)`,
  (* {{{ proof *)
  [
  REPEAT GEN_TAC;
    STRIP_TAC;
    SUBGOAL_THEN `&0 <= y` MP_TAC;
  ASM_MESON_TAC [REAL_LE_TRANS];
  ASM_MESON_TAC[Collect_geom.POW2_COND;POW_2_SQRT_ABS;
		REAL_ARITH `&0 <= x ==> (x = abs x)`];
  ]);;
  (* }}} *)

let ineq_square2  = prove_by_refinement(
  `(&0 <= a1 /\ &0 <= a2 /\ &0 <= a3 /\ &0 <= a4 /\ &0 <= a5 /\ &0 <= a6 /\
      &0 <= b1 /\ &0 <= b2 /\ &0 <= b3 /\ &0 <= b4 /\ &0 <= b5 /\ &0 <= b6 )
   /\ 
    (!x1 x2 x3 x4 x5 x6. 
       ineq [(a1 pow 2,x1,b1 pow 2);(a2 pow 2,x2,b2 pow 2);(a3 pow 2,x3,b3 pow 2);
        (a4 pow 2,x4,b4 pow 2);(a5 pow 2,x5,b5 pow 2);(a6 pow 2,x6,b6 pow 2)]
          (P (sqrt x1) (sqrt x2) (sqrt x3) (sqrt x4) (sqrt x5) (sqrt x6))) ==>
    (!y1 y2 y3 y4 y5 y6.
      ineq 
        [(a1,y1,b1);(a2,y2,b2);(a3,y3,b3);(a4,y4,b4);(a5,y5,b5);(a6,y6,b6)]
          (P y1 y2 y3 y4 y5 y6))`,
  (* {{{ proof *)
  [
  REWRITE_TAC[ineq];
  REPEAT STRIP_TAC;
  FIRST_X_ASSUM (fun t-> MP_TAC (SPECL [`y1 pow 2`;`y2 pow 2`;`y3 pow 2`;`y4 pow 2`;`y5 pow 2`;`y6 pow 2`] t));
  SUBGOAL_THEN `(sqrt (y1 pow 2) = y1) /\ (sqrt (y2 pow 2) = y2) /\  (sqrt (y3 pow 2) = y3) /\   (sqrt (y4 pow 2) = y4) /\   (sqrt (y5 pow 2) = y5) /\   (sqrt (y6 pow 2) = y6)` (fun t -> REWRITE_TAC[t]);
  ASM_MESON_TAC[ineq_lemma_b];
  REWRITE_TAC[TAUT `(a ==> b ==> c) <=> (a /\ b ==> c)`];
  DISCH_THEN MATCH_MP_TAC;
  ASM_MESON_TAC[ineq_lemma_b];
  ]);;
  (* }}} *)

let ineq_square2_9  = prove_by_refinement(
  `(&0 <= a1 /\ &0 <= a2 /\ &0 <= a3 /\ &0 <= a4 /\ &0 <= a5 /\ &0 <= a6 /\
     &0 <= a7 /\ &0 <= a8 /\ &0 <= a9 /\
      &0 <= b1 /\ &0 <= b2 /\ &0 <= b3 /\ &0 <= b4 /\ &0 <= b5 /\ &0 <= b6 /\
    &0 <= b7 /\ &0 <= b8 /\ &0 <= b9
)
   /\ 
    (!x1 x2 x3 x4 x5 x6 x7 x8 x9. 
       ineq [(a1 pow 2,x1,b1 pow 2);(a2 pow 2,x2,b2 pow 2);(a3 pow 2,x3,b3 pow 2);
        (a4 pow 2,x4,b4 pow 2);(a5 pow 2,x5,b5 pow 2);(a6 pow 2,x6,b6 pow 2);
        (a7 pow 2,x7,b7 pow 2);(a8 pow 2,x8,b8 pow 2);(a9 pow 2,x9,b9 pow 2)]
          (P (sqrt x1) (sqrt x2) (sqrt x3) (sqrt x4) (sqrt x5) (sqrt x6) (sqrt x7) (sqrt x8) (sqrt x9))) ==>
    (!y1 y2 y3 y4 y5 y6 y7 y8 y9.
      ineq 
        [(a1,y1,b1);(a2,y2,b2);(a3,y3,b3);(a4,y4,b4);(a5,y5,b5);(a6,y6,b6);(a7,y7,b7);(a8,y8,b8);(a9,y9,b9)]
          (P y1 y2 y3 y4 y5 y6 y7 y8 y9))`,
  (* {{{ proof *)
  [
  REWRITE_TAC[ineq];
  REPEAT STRIP_TAC;
  FIRST_X_ASSUM (fun t-> MP_TAC (SPECL [`y1 pow 2`;`y2 pow 2`;`y3 pow 2`;`y4 pow 2`;`y5 pow 2`;`y6 pow 2`;`y7 pow 2`;`y8 pow 2`;`y9 pow 2`] t));
  SUBGOAL_THEN `(sqrt (y1 pow 2) = y1) /\ (sqrt (y2 pow 2) = y2) /\  (sqrt (y3 pow 2) = y3) /\   (sqrt (y4 pow 2) = y4) /\   (sqrt (y5 pow 2) = y5) /\   (sqrt (y6 pow 2) = y6) /\  (sqrt (y7 pow 2) = y7) /\  (sqrt (y8 pow 2) = y8) /\  (sqrt (y9 pow 2) = y9)` (fun t -> REWRITE_TAC[t]);
  ASM_MESON_TAC[ineq_lemma_b];
  REWRITE_TAC[TAUT `(a ==> b ==> c) <=> (a /\ b ==> c)`];
  DISCH_THEN MATCH_MP_TAC;
  ASM_MESON_TAC[ineq_lemma_b];
  ]);;
  (* }}} *)

let sqrt8_nn = prove_by_refinement(
  `&0 <= sqrt8`,
  (* {{{ proof *)
  [
  REWRITE_TAC[Sphere.sqrt8];
  MATCH_MP_TAC SQRT_POS_LE;
  REAL_ARITH_TAC;
  ]);;
  (* }}} *)

let sqrt2_nn = prove_by_refinement(
  `&0 <= sqrt2`,
  (* {{{ proof *)
  [
  REWRITE_TAC[Sphere.sqrt2];
  MATCH_MP_TAC SQRT_POS_LE;
  REAL_ARITH_TAC;
  ]);;
  (* }}} *)

let sqrt3_nn = prove_by_refinement(
  `&0 <= sqrt(&3)`,
  (* {{{ proof *)
  [
  MATCH_MP_TAC SQRT_POS_LE;
  REAL_ARITH_TAC;
  ]);;
  (* }}} *)

let basic_constants_nn = [
 REAL_ARITH `&0 <= #2.18 /\ &0 <= &2 /\ &0 <= #2.52 /\ #2.0 = &2 /\ #2 = &2 /\ &0 <= #2.25 `;
 sqrt8_nn;sqrt2_nn;
 ];;

let abc_quadratic = prove (`abc_of_quadratic (\x. a * (x pow 2) + b * x + c) = (a,b,c)`,
 REWRITE_TAC[Sphere.abc_of_quadratic] THEN
 (REPEAT LET_TAC) THEN
 REWRITE_TAC[PAIR_EQ] THEN
 REPEAT(FIRST_X_ASSUM MP_TAC)THEN
 ARITH_TAC);;

let delta_quadratic = prove( `-- delta_x x1 x2 x3 x4 x5 x6 = 
  (x1) * (x4 pow 2) + (x1*x1 + (x3 - x5)*(x2 - x6) 
   - x1*(x2 + x3 + x5 + x6)) * x4 
   + ( x1*x3*x5 + x1*x2*x6 - x3*(x1 + x2 - x3 + x5 - x6)*x6 
    - x2*x5*(x1 - x2 + x3 - x5 + x6) ) `,
REWRITE_TAC[Sphere.delta_x] THEN
ARITH_TAC);;

let edge_flat_rewrite = 
 REWRITE_RULE[abc_quadratic;delta_quadratic] Sphere.edge_flat;;

let enclosed_rewrite = 
  REWRITE_RULE[abc_quadratic] 
  (strip_let(REWRITE_RULE[Mur.muRa;Cayleyr.cayleyR_quadratic] Enclosed.enclosed));;

let quad_root_plus_curry = 
  REWRITE_RULE[Sphere.quadratic_root_plus] quadratic_root_plus_curry;;

let y_of_x_e = prove(`!y1 y2 y3 y4 y5 y6. y_of_x f y1 y2 y3 y4 y5 y6 =
     f (y1*y1) (y2*y2) (y3*y3) (y4*y4) (y5*y5) (y6*y6)`,
     REWRITE_TAC[Sphere.y_of_x]);;

let vol_y_e = prove(`!y1 y2 y3 y4 y5 y6. vol_y y1 y2 y3 y4 y5 y6 = 
    y_of_x vol_x y1 y2 y3 y4 y5 y6`,
    REWRITE_TAC[Sphere.vol_y]);;

let rad2_y_e = prove(`!y1 y2 y3 y4 y5 y6. rad2_y y1 y2 y3 y4 y5 y6 = 
    y_of_x rad2_x y1 y2 y3 y4 y5 y6`,
    REWRITE_TAC[Sphere.rad2_y]);;

let rad2_x_y = prove_by_refinement(
  `!a1 a2 a3 a4 a5 a6 x1 x2 x3 x4 x5 x6.
    (a1 pow 2 <= x1) /\ (a2 pow 2 <= x2) /\ (a3 pow 2 <= x3) /\ (a4 pow 2 <= x4) /\
    (a5 pow 2 <= x5) /\ (a6 pow 2 <= x6) ==> 
    (rad2_y (sqrt x1) (sqrt x2) (sqrt x3) (sqrt x4) (sqrt x5) (sqrt x6) = 
         rad2_x x1 x2 x3 x4 x5 x6)`,
  (* {{{ proof *)
  [
  REWRITE_TAC[Sphere.rad2_y;y_of_x_e;LET_DEF;LET_END_DEF];
  ASM_MESON_TAC[sq_pow2];
  ]);;
  (* }}} *)

let delta_x4_delta4_y = prove_by_refinement(
  `!a1 a2 a3 a4 a5 a6 x1 x2 x3 x4 x5 x6.
    (a1 pow 2 <= x1) /\ (a2 pow 2 <= x2) /\ (a3 pow 2 <= x3) /\ (a4 pow 2 <= x4) /\
    (a5 pow 2 <= x5) /\ (a6 pow 2 <= x6) ==> 
  (delta4_y (sqrt x1) (sqrt x2) (sqrt x3) (sqrt x4) (sqrt x5) (sqrt x6) = delta_x4 x1 x2 x3 x4 x5 x6)`,
  (* {{{ proof *)
  [
  REWRITE_TAC[Sphere.delta4_y;y_of_x_e;LET_DEF;LET_END_DEF];
  ASM_MESON_TAC[sq_pow2];
  ]);;
  (* }}} *)

let dih_x_y = prove_by_refinement(
  `!a1 a2 a3 a4 a5 a6 x1 x2 x3 x4 x5 x6.
    (a1 pow 2 <= x1) /\ (a2 pow 2 <= x2) /\ (a3 pow 2 <= x3) /\ (a4 pow 2 <= x4) /\
    (a5 pow 2 <= x5) /\ (a6 pow 2 <= x6) ==> 
    (dih_y (sqrt x1) (sqrt x2) (sqrt x3) (sqrt x4) (sqrt x5) (sqrt x6) = 
         dih_x x1 x2 x3 x4 x5 x6)`,
  (* {{{ proof *)
  [
  REWRITE_TAC[Sphere.dih_y;LET_DEF;LET_END_DEF];
  ASM_MESON_TAC[sq_pow2];
  ]);;
  (* }}} *)

let dih2_x_y = prove_by_refinement(
  `!a1 a2 a3 a4 a5 a6 x1 x2 x3 x4 x5 x6.
    (a1 pow 2 <= x1) /\ (a2 pow 2 <= x2) /\ (a3 pow 2 <= x3) /\ (a4 pow 2 <= x4) /\
    (a5 pow 2 <= x5) /\ (a6 pow 2 <= x6) ==> 
    (dih2_y (sqrt x1) (sqrt x2) (sqrt x3) (sqrt x4) (sqrt x5) (sqrt x6) = 
         dih2_x x1 x2 x3 x4 x5 x6)`,
  (* {{{ proof *)
  [
  REWRITE_TAC[Sphere.dih2_y;Sphere.dih2_x;Sphere.dih_y;LET_DEF;LET_END_DEF];
  ASM_MESON_TAC[sq_pow2];
  ]);;
  (* }}} *)

let dih3_x_y = prove_by_refinement(
  `!a1 a2 a3 a4 a5 a6 x1 x2 x3 x4 x5 x6.
    (a1 pow 2 <= x1) /\ (a2 pow 2 <= x2) /\ (a3 pow 2 <= x3) /\ (a4 pow 2 <= x4) /\
    (a5 pow 2 <= x5) /\ (a6 pow 2 <= x6) ==> 
    (dih3_y (sqrt x1) (sqrt x2) (sqrt x3) (sqrt x4) (sqrt x5) (sqrt x6) = 
         dih3_x x1 x2 x3 x4 x5 x6)`,
  (* {{{ proof *)
  [
  REWRITE_TAC[Sphere.dih3_y;Sphere.dih3_x;Sphere.dih_y;LET_DEF;LET_END_DEF];
  ASM_MESON_TAC[sq_pow2];
  ]);;
  (* }}} *)

let delta_x_y = prove_by_refinement(
  `!a1 a2 a3 a4 a5 a6 x1 x2 x3 x4 x5 x6.
    (a1 pow 2 <= x1) /\ (a2 pow 2 <= x2) /\ (a3 pow 2 <= x3) /\ (a4 pow 2 <= x4) /\
    (a5 pow 2 <= x5) /\ (a6 pow 2 <= x6) ==> 
    (delta_y (sqrt x1) (sqrt x2) (sqrt x3) (sqrt x4) (sqrt x5) (sqrt x6) = 
         delta_x x1 x2 x3 x4 x5 x6)`,
  (* {{{ proof *)
  [
  REWRITE_TAC[Sphere.delta_y;LET_DEF;LET_END_DEF];
  ASM_MESON_TAC[sq_pow2];
  ]);;
  (* }}} *)

let vol_x_y = prove_by_refinement(
  `!a1 a2 a3 a4 a5 a6 x1 x2 x3 x4 x5 x6.
    (a1 pow 2 <= x1) /\ (a2 pow 2 <= x2) /\ (a3 pow 2 <= x3) /\ (a4 pow 2 <= x4) /\
    (a5 pow 2 <= x5) /\ (a6 pow 2 <= x6) ==> 
    (vol_y (sqrt x1) (sqrt x2) (sqrt x3) (sqrt x4) (sqrt x5) (sqrt x6) = 
         vol_x x1 x2 x3 x4 x5 x6)`,
  (* {{{ proof *)
  [
  REWRITE_TAC[Sphere.vol_y;Sphere.y_of_x;LET_DEF;LET_END_DEF];
  ASM_MESON_TAC[sq_pow2];
  ]);;
  (* }}} *)


let sqrt8_2 = prove_by_refinement(
  `sqrt8 * sqrt8 = #8.0`,
  (* {{{ proof *)
  [
  REWRITE_TAC[Sphere.sqrt8];
  MESON_TAC[REAL_POW_2;SQRT_WORKS;REAL_ARITH `&8 = #8.0 /\ &0 <= &8`];
  ]);;
  (* }}} *)

let dih_x_sym = prove_by_refinement(
  `!x1 x2 x3 x4 x5 x6. dih_x x1 x2 x3 x4 x5 x6 = dih_x x1 x3 x2 x4 x6 x5`,
  (* {{{ proof *)
  [
  REWRITE_TAC[Sphere.dih_x;LET_DEF;LET_END_DEF];
  REPEAT GEN_TAC;
  REPEAT AP_TERM_TAC;
  REWRITE_TAC[Sphere.delta_x;Sphere.delta_x4;PAIR_EQ];
  CONJ_TAC THEN AP_TERM_TAC THEN REAL_ARITH_TAC;
  ]);;
  (* }}} *)

let dih_x_sym2 = prove_by_refinement(
  `!x1 x2 x3 x4 x5 x6. dih_x x1 x2 x3 x4 x5 x6 = dih_x x1 x5 x6 x4 x2 x3`,
  (* {{{ proof *)
  [
  REWRITE_TAC[Sphere.dih_x;LET_DEF;LET_END_DEF];
  REPEAT GEN_TAC;
  REPEAT AP_TERM_TAC;
  REWRITE_TAC[Sphere.delta_x;Sphere.delta_x4;PAIR_EQ];
  CONJ_TAC THEN AP_TERM_TAC THEN REAL_ARITH_TAC;
  ]);;
  (* }}} *)


let dih_y_sym = prove_by_refinement(
  `!y1 y2 y3 y4 y5 y6. dih_y y1 y2 y3 y4 y5 y6 = dih_y y1 y3 y2 y4 y6 y5`,
  (* {{{ proof *)
  [
  REWRITE_TAC[Sphere.dih_y;LET_DEF;LET_END_DEF];
  MESON_TAC[dih_x_sym];
  ]);;
  (* }}} *)

let dih_y_sym2 = prove_by_refinement(
  `!y1 y2 y3 y4 y5 y6. dih_y y1 y2 y3 y4 y5 y6 = dih_y y1 y5 y6 y4 y2 y3`,
  (* {{{ proof *)
  [
  REWRITE_TAC[Sphere.dih_y;LET_DEF;LET_END_DEF];
  MESON_TAC[dih_x_sym2];
  ]);;
  (* }}} *)

let sol_y_123 = prove_by_refinement(
  `!y1 y2 y3 y4 y5 y6. sol_y y1 y2 y3 y4 y5 y6 = 
    dih_y y1 y2 y3 y4 y5 y6 + dih2_y y1 y2 y3 y4 y5 y6 + dih3_y y1  y2 y3 y4 y5 y6 - pi`,
  (* {{{ proof *)
  [
  REWRITE_TAC[Sphere.sol_y;Sphere.dih2_y;Sphere.dih3_y];
  REPEAT GEN_TAC;
  MATCH_MP_TAC (REAL_ARITH `(b = b') ==> (a + b + c -pi = a + b' + c - pi)`);
  REWRITE_TAC[Sphere.dih_y;LET_DEF;LET_END_DEF];
  MESON_TAC[dih_x_sym];
  ]);;
  (* }}} *)

let rhazim2_alt = prove_by_refinement(
  `!y1 y2 y3 y4 y5 y6. rhazim2 y1 y2 y3 y4 y5 y6 = rho y2 * dih2_y y1 y2 y3 y4 y5 y6`,
  (* {{{ proof *)
  [
  REWRITE_TAC[Sphere.rhazim2;Sphere.node2_y;Sphere.rhazim;Sphere.dih2_y;Sphere.dih_y;LET_DEF;LET_END_DEF];
  MESON_TAC[dih_x_sym];
  ]);;
  (* }}} *)

let rhazim3_alt = prove_by_refinement(
  `!y1 y2 y3 y4 y5 y6. rhazim3 y1 y2 y3 y4 y5 y6 = rho y3 * dih3_y y1 y2 y3 y4 y5 y6`,
  (* {{{ proof *)
  [
  REWRITE_TAC[Sphere.rhazim3;Sphere.node3_y;Sphere.rhazim;Sphere.dih3_y;Sphere.dih_y;LET_DEF;LET_END_DEF];
  ]);;
  (* }}} *)

let taum_123 = prove_by_refinement(
  `!y1 y2 y3 y4 y5 y6. taum y1 y2 y3 y4 y5 y6 = rhazim y1 y2 y3 y4 y5 y6 + rhazim2 y1 y2 y3 y4 y5 y6 + rhazim3 y1  y2 y3 y4 y5 y6 - (&1 + const1)* pi`,
  (* {{{ proof *)
  [
  REPEAT GEN_TAC;
  REWRITE_TAC[Sphere.taum;sol_y_123;Sphere.lnazim];
    SUBGOAL_THEN `dih_y y2 y3 y1 y5 y6 y4 = dih2_y y1 y2 y3 y4 y5 y6 /\ dih_y y3 y1 y2 y6 y4 y5 = dih3_y y1 y2 y3 y4 y5 y6` (fun t-> REWRITE_TAC[t]);
  REWRITE_TAC[Sphere.dih_y;Sphere.dih2_y;Sphere.dih3_y;LET_DEF;LET_END_DEF] THEN MESON_TAC[dih_x_sym];
  REWRITE_TAC[Sphere.rhazim;rhazim2_alt;rhazim3_alt;Sphere.rho];
  REAL_ARITH_TAC;
  ]);;
  (* }}} *)

let tauq_x_y = prove_by_refinement(
  `!x1 x2 x3 x4 x5 x6 x7 x8 x9.
   (tauq (sqrt x1) (sqrt x2) (sqrt x3) (sqrt x4) (sqrt x5) (sqrt x6) (sqrt x7) (sqrt x8) (sqrt x9) = 
    taum_x x1 x2 x3 x4 x5 x6 + taum_x x7 x2 x3 x4 x8 x9)`,
  (* {{{ proof *)
  [
  REWRITE_TAC[Sphere.tauq;Sphere.taum_x;taum_123;LET_DEF;LET_END_DEF;Sphere.rhazim_x;Sphere.rhazim2_x;Sphere.rhazim3_x];
  ]);;					   	  
  (* }}} *)					   


let rho_alt = prove_by_refinement(
  `!y. rho y = &1 + const1 *(y - &2) / (#0.52)`,
  (* {{{ proof *)
  [
  GEN_TAC;
  REWRITE_TAC[Sphere.rho;Sphere.ly;Sphere.interp;REAL_ARITH `#2.52 - &2 = #0.52`];
  REAL_ARITH_TAC;
  ]);;
  (* }}} *)

(* renamed from rho_x to avoid clash with rho_x in sphere.hl *)

let rho_sqrtx = prove_by_refinement(
  `!x. rho (sqrt x) = &1 + const1 * (sqrt x - &2) / (#0.52)`,
  (* {{{ proof *)
  [
  REWRITE_TAC[rho_alt];
  ]);;
  (* }}} *)

let lfun_ly = prove_by_refinement(
  `!h. lfun h = ly (&2 * h)`,
  (* {{{ proof *)
  [
  REWRITE_TAC[Sphere.lfun;Sphere.ly;Sphere.interp;Sphere.h0];
  REAL_ARITH_TAC;
  ]);;
  (* }}} *)

let lfun1 = prove_by_refinement(
  `lfun (&1) = (&1)`,
  (* {{{ proof *)
  [
  REWRITE_TAC[Sphere.lfun;Sphere.h0];
  REAL_ARITH_TAC;
  ]);;
  (* }}} *)

let rhazim2_alt = prove_by_refinement(
  `!y1 y2 y3 y4 y5 y6. rhazim2 y1 y2 y3 y4 y5 y6 = rho y2 * dih2_y y1 y2 y3 y4 y5 y6`,
  (* {{{ proof *)
  [
  REWRITE_TAC[Sphere.rhazim2;Sphere.node2_y;Sphere.rhazim;Sphere.dih2_y;];
  MESON_TAC[dih_y_sym];
  ]);;
  (* }}} *)

let rhazim3_alt = prove_by_refinement(
  `!y1 y2 y3 y4 y5 y6. rhazim3 y1 y2 y3 y4 y5 y6 = rho y3 * dih3_y y1 y2 y3 y4 y5 y6`,
  (* {{{ proof *)
  [
  REWRITE_TAC[Sphere.rhazim3;Sphere.node3_y;Sphere.rhazim;Sphere.dih3_y;];
  ]);;
  (* }}} *)

let beta_bump_force_x = prove_by_refinement(
  `!x1 x2 x3 x4 x5 x6. beta_bump_force_y (sqrt x1) (sqrt x2) (sqrt x3) (sqrt x4) (sqrt x5) (sqrt x6) =
    halfbump_x1 x1 x2 x3 x4 x5 x6 - halfbump_x4 x1 x2 x3 x4 x5 x6`,
  (* {{{ proof *)
  [
  REWRITE_TAC[Sphere.beta_bump_force_y;halfbump_x1;halfbump_x4;halfbump_x];
  ]);;
  (* }}} *)

let halfbump_x_expand = prove_by_refinement(
  `!x. &0 <= x ==> (halfbump_x x = 
   -- (&4398119 / &2376200) + (&17500 / &11881) * sqrt x - (&31250 / &106929) * x)`,
  (* {{{ proof *)
  [
  REWRITE_TAC[halfbump_x;Sphere.bump;Sphere.h0;Sphere.hplus];
  REPEAT STRIP_TAC;
  REWRITE_TAC[REAL_ARITH`(a/ &2 - b) pow 2 = (a pow 2) / &4 - a * b + b pow 2`];
  ASM_SIMP_TAC[SQRT_POW_2];
  REAL_ARITH_TAC;
  ]);;
  (* }}} *)

let vol4f_palt = prove_by_refinement(
  `!f y1 y2 y3 y4 y5 y6. vol4f y1 y2 y3 y4 y5 y6 f  = 
    (-- &8 * mm1)
    + 
    (&4 * mm1 / pi) * 
    (dih_y y1 y2 y3 y4 y5 y6 + dih2_y y1 y2 y3 y4 y5 y6 +
    dih3_y y1 y2 y3 y4 y5 y6 +
    dih4_y y1 y2 y3 y4 y5 y6 + dih5_y y1 y2 y3 y4 y5 y6 + 
     dih6_y y1 y2 y3 y4 y5 y6)
    +
    (-- &8 * mm2 / pi) * 
    (f (y1/ &2) * dih_y y1 y2 y3 y4 y5 y6 + 
    f(y2/ &2) * dih2_y y1 y2 y3 y4 y5 y6 +  
    f (y3/ &2) * dih3_y y1 y2 y3 y4 y5 y6 +
    f (y4/ &2) * dih4_y y1 y2 y3 y4 y5 y6 + 
    f(y5/ &2) * dih5_y y1 y2 y3 y4 y5 y6 + 
    f(y6/ &2) * dih6_y y1 y2 y3 y4 y5 y6)
    `,
  (* {{{ proof *)
  [
  REPEAT GEN_TAC;
  REWRITE_TAC[Sphere.vol4f;Sphere.sol_y;lfun_ly;REAL_ARITH `&2 * y / &2 = y`;Sphere.rhazim;rhazim2_alt;rhazim3_alt;Sphere.rhazim4;Sphere.rhazim5;Sphere.rhazim6;];
  SUBGOAL_THEN `dih_y y2 y3 y1 y5 y6 y4 = dih2_y y1 y2 y3 y4 y5 y6 /\ 
   dih_y y3 y1 y2 y6 y4 y5 = dih3_y y1 y2 y3 y4 y5 y6 /\ 
   dih_y y3 y1 y2 y6 y4 y5 = dih3_y y1 y2 y3 y4 y5 y6 /\ 
   dih_y y4 y3 y5 y1 y6 y2 = dih4_y y1 y2 y3 y4 y5 y6 /\ 
   dih_y y5 y1 y6 y2 y4 y3 = dih5_y y1 y2 y3 y4 y5 y6 /\ 
   dih_y y6 y1 y5 y3 y4 y2 = dih6_y y1 y2 y3 y4 y5 y6 /\ 
   dih_y y6 y4 y2 y3 y1 y5 = dih6_y y1 y2 y3 y4 y5 y6 /\ 
   dih_y y2 y6 y4 y5 y3 y1 = dih2_y y1 y2 y3 y4 y5 y6 /\ 
   dih_y y1 y5 y6  y4 y2 y3 = dih_y y1 y2 y3 y4 y5 y6 /\ 
   dih_y y5 y6 y1 y2 y3 y4 = dih5_y y1 y2 y3 y4 y5 y6 /\ 
   dih_y y4 y5 y3 y1 y2 y6 = dih4_y y1 y2 y3 y4 y5 y6 /\ 
   dih_y y5 y3 y4 y2 y6 y1 = dih5_y y1 y2 y3 y4 y5 y6 /\ 
   dih_y y3 y4 y5 y6 y1 y2 = dih3_y y1 y2 y3 y4 y5 y6 /\ 
   dih_y y4 y2 y6 y1 y5 y3 = dih4_y y1 y2 y3 y4 y5 y6` (fun t-> REWRITE_TAC[t]);
  REWRITE_TAC[Sphere.dih2_y;Sphere.dih3_y;Sphere.dih4_y;Sphere.dih5_y;Sphere.dih6_y];
  MESON_TAC[dih_y_sym;dih_y_sym2];
  ONCE_REWRITE_TAC[REAL_ARITH `x  = y <=> x - y = &0`];
  ABBREV_TAC `a = mm1/pi `;
  SUBGOAL_THEN `mm1 = a * pi ` (fun t->REWRITE_TAC[t]);
  POP_ASSUM MP_TAC;
  SUBGOAL_THEN  `~(pi = &0)` MP_TAC;
  SIMP_TAC[PI_POS;REAL_ARITH `&0 < x ==> ~(x= &0)`];
  CONV_TAC REAL_FIELD;
  REAL_ARITH_TAC;
  ]);;
  (* }}} *)


let edge_flat2_x_rewrite = prove_by_refinement(
  `!x1 x2 x3 x4 x5 x6. (&0 <= x1 /\ &0 <= x2 /\ &0 <= x3 /\ &0 <= x5 /\
    &0 <= x6) ==> (edge_flat2_x x1 x2 x3 x4 x5 x6 = 
    (sqrt (quadratic_root_plus
           (x1,
            (x1) * x1 +
            (x3 - x5) * (x2 - x6) -
            (x1) * (x2 + x3 + x5 + x6),
            (x1) * (x3) * x5 +
            (x1) * (x2) * x6 -
            (x3) *
            (x1 + x2 - x3 + x5 - x6) * x6 -
            (x2) *
            (x5) *
            (x1 - x2 + x3 - x5 + x6)))) pow 2)`,
  (* {{{ proof *)
  [
  REWRITE_TAC[Sphere.edge_flat2_x];
  REPEAT STRIP_TAC;
  AP_THM_TAC;
  AP_TERM_TAC;
  REWRITE_TAC[edge_flat_rewrite];
  ASM_SIMP_TAC[REAL_ARITH `sqrt x * sqrt x = (sqrt x) pow 2`;SQRT_POW2;SQRT_WORKS];
  ]);;
  (* }}} *)

let edge_quadratic = prove_by_refinement(
  `!x1 x2 x3 x5 x6. quadratic_root_plus (x1,
            x1 * x1 +
            (x3 - x5) * (x2 - x6) -
            x1 * (x2 + x3 + x5 + x6),
            x1 * x3 * x5 +
            (x1) * (x2) * x6 -  (x3) * (x1 + x2 - x3 + x5 - x6) * x6 -
	      (x2) * (x5) * (x1 - x2 + x3 - x5 + x6)) =
  (-- (x1 * x1) + x1*x2 + x1*x3 - x2*x3 + x1*x5 + x2*x5 + x1*x6 + x3*x6 - 
  x5*x6 + sqrt(ups_x x1 x3 x5 * ups_x x1 x2 x6))/(&2*x1)`,
  (* {{{ proof *)
  [
  REWRITE_TAC[Sphere.quadratic_root_plus;Sphere.ups_x];
  REPEAT STRIP_TAC;
  AP_THM_TAC;
  AP_TERM_TAC;
  REWRITE_TAC[REAL_ARITH `--(x1 * x1 + (x3 - x5) * (x2 - x6) - x1 * (x2 + x3 + x5 + x6)) + a = -- (x1 * x1) + x1 * x2 + x1 * x3 - x2 * x3 + x1 * x5 + x2 * x5 + x1 * x6 + x3 * x6 - x5 * x6 + a`];
  REPEAT (AP_THM_TAC ORELSE AP_TERM_TAC);
  REAL_ARITH_TAC;
  ]);;
  (* }}} *)

let lmfun0 = prove_by_refinement(
  `!y. &2 * h0 <= y ==> (lmfun (y/ &2)  = &0)`,
  (* {{{ proof *)
  [
  REWRITE_TAC[Sphere.lmfun;REAL_ARITH `&2 * h0 <= y <=> (~(y/ &2 <= h0) \/ (y/ &2 = h0))`];
  GEN_TAC;
  DISCH_THEN DISJ_CASES_TAC THEN ASM_REWRITE_TAC[REAL_ARITH `h0 <= h0`;REAL_FIELD `(h0 - h0)/(h0 - &1) = &0`];
  ]);;
  (* }}} *)

let lmfun_lfun = prove_by_refinement(
  `!y. y <= &2 * h0 ==> (lmfun (y/ &2) = lfun(y/ &2))`,
  (* {{{ proof *)
  [
  REWRITE_TAC[Sphere.lmfun;Sphere.lfun;REAL_ARITH `y <= &2 * h0 <=> y/ &2 <= h0`];
  REPEAT STRIP_TAC THEN ASM_REWRITE_TAC[];
  ]);;
  (* }}} *)

let lmfun_lfun2 = prove_by_refinement(
  `!y. y <=  h0 ==> (lmfun (y) = lfun(y))`,
  (* {{{ proof *)
  [
  REWRITE_TAC[Sphere.lmfun;Sphere.lfun];
  REPEAT STRIP_TAC THEN ASM_REWRITE_TAC[];
  ]);;
  (* }}} *)

(* compute hminus *)

let quartic_has_real_derivative = prove_by_refinement(
  `!x c0 c1 c2 c3 c4 . ((\x. c0 * &1 + c1 * x pow 1 + c2 * x pow 2 + c3 * x pow 3 + c4 * x pow 4) has_real_derivative (c0 * &0 + c1 * (&1 * x pow (1-1) * &1) + c2 * (&2 * x pow (2-1)) * &1 + c3 * (&3 * x pow (3-1)) * &1 + c4 * (&4 * x pow (4-1)) * &1)) (atreal x)`,
  (* {{{ proof *)
  [
  REPEAT STRIP_TAC THEN 
REPEAT (MATCH_MP_TAC HAS_REAL_DERIVATIVE_ADD THEN CONJ_TAC) THEN MATCH_MP_TAC (HAS_REAL_DERIVATIVE_LMUL_ATREAL) THEN REWRITE_TAC[HAS_REAL_DERIVATIVE_CONST]  THEN REWRITE_TAC[REAL_ARITH `(a * b) * c = a * b * (c:real)`] THEN MATCH_MP_TAC (HAS_REAL_DERIVATIVE_POW_ATREAL) THEN REWRITE_TAC[HAS_REAL_DERIVATIVE_ID];
  ]);;
  (* }}} *)

let POLY_CONTINUITY_TAC =  
 REPEAT (MATCH_MP_TAC REAL_CONTINUOUS_ON_ADD THEN 
 CONJ_TAC) THEN 
 REPEAT (MATCH_MP_TAC REAL_CONTINUOUS_ON_LMUL) THEN 
 (MATCH_MP_TAC REAL_CONTINUOUS_ON_POW) THEN 
 REWRITE_TAC[REAL_CONTINUOUS_ON_ID];;

let quartic_continuous_on = prove_by_refinement(
  `!s c0 c1 c2 c3 c4 .  (\x. c0 * x pow 0 + c1 * x pow 1 + c2 * x pow 2 + c3 * x pow 3 + c4 * x pow 4) real_continuous_on s`,
  (* {{{ proof *)
  [
  REPEAT STRIP_TAC;
  POLY_CONTINUITY_TAC;
  ]);;
  (* }}} *)

let marchal_minus_lfun = prove_by_refinement(
  `!h. marchal_quartic h - lfun h = 
   (inv(&65 * &1627 * (sqrt(&2) - &1))) * (h - &1)*
    (( -- &512505 + &770958*sqrt(&2)) * h pow 0 + 
    ( -- &364208 -  &1295359*sqrt(&2)) * h pow 1 +
   ( &1295359 + &585000*sqrt(&2))* h pow 2 + 
   (-- &585000) * h pow 3)`,
  (* {{{ proof *)
  [
  REWRITE_TAC[Sphere.marchal_quartic;Sphere.lfun;REAL_ARITH `x/y = x * inv y /\ #1.26 - &1 = &13/ &50 /\ #1.3254 - &1 = &1627/ &5000 /\ &65 = &5 * &13`;REAL_INV_INV;REAL_INV_MUL;Sphere.hplus;Sphere.h0];
  GEN_TAC;
  SUBGOAL_THEN `&0 < sqrt(&2) - &1 ` MP_TAC THENL [ALL_TAC;CONV_TAC REAL_FIELD];
  SUBGOAL_THEN `#1.414213 < sqrt(&2) ` MP_TAC THENL[ALL_TAC;REAL_ARITH_TAC];
  REWRITE_TAC[GSYM Sphere.sqrt2;Flyspeck_constants.bounds];
  ]);;
  (* }}} *)

let hminus_cont = prove_by_refinement(
  `!s. (\h. marchal_quartic h - lfun h) real_continuous_on s`,
  (* {{{ proof *)
  [
  GEN_TAC;
  REWRITE_TAC[marchal_minus_lfun];
  MATCH_MP_TAC REAL_CONTINUOUS_ON_LMUL;
 MATCH_MP_TAC REAL_CONTINUOUS_ON_MUL THEN REWRITE_TAC[REAL_ARITH `h - &1 = (-- &1) * h pow 0 + &1 * h pow 1`] THEN CONJ_TAC THEN POLY_CONTINUITY_TAC;
  ]);;
  (* }}} *)

let sqrt2_lb = prove_by_refinement(
  `#1.414213 < sqrt2 `,
  (* {{{ proof *)
  [
  REWRITE_TAC[Flyspeck_constants.bounds];
  ]);;
  (* }}} *)

let STRIP_NN_TAC = REPEAT (CONJ_TAC ORELSE MATCH_MP_TAC REAL_LE_MUL ORELSE CHANGED_TAC (REWRITE_TAC[REAL_LE_INV_EQ;REAL_ARITH `a >= &0 <=> &0 <= a /\ ((a*b)*c = a*b*c)`;REAL_INV_MUL]));;

let hminus_exists = prove_by_refinement(
  `?x. (#1.2 <= x /\ x < #1.3 /\ marchal_quartic x = lmfun x)`,
  (* {{{ proof *)
  [
    SUBGOAL_THEN `(?x. x IN real_interval [#1.2,#1.26] /\ (\x. marchal_quartic x - lfun x) x = &0)` MP_TAC;
    MATCH_MP_TAC REAL_IVT_INCREASING;
    BETA_TAC;
    REWRITE_TAC[hminus_cont];
    REWRITE_TAC[marchal_minus_lfun;GSYM Sphere.sqrt2];
   ASSUME_TAC sqrt2_lb;
   SUBGOAL_THEN `&0 < sqrt2 - &1` ASSUME_TAC;
   POP_ASSUM MP_TAC;
   REAL_ARITH_TAC;
   REWRITE_TAC[(* REAL_INV_MUL;REAL_INV_INV; *) REAL_ARITH `#1.2- &1 = &1 / &5 /\ #1.26 - &1 = &13 / &50`];
    CONJ_TAC THENL[REAL_ARITH_TAC;ALL_TAC];
    SUBGOAL_THEN `sqrt2 < #1.414214` ASSUME_TAC;
   REWRITE_TAC[Flyspeck_constants.bounds];
   SUBGOAL_THEN ` ((-- &512505 + &770958 * sqrt2) * #1.26 pow 0 +  (-- &364208 - &1295359 * sqrt2) * #1.26 pow 1 +  (&1295359 + &585000 * sqrt2) * #1.26 pow 2 +  -- &585000 * #1.26 pow 3) = -- &212787729/ &2500 + &3377583 *sqrt2 /(&50)` (fun t-> REWRITE_TAC[t]);
   REAL_ARITH_TAC;
   SUBGOAL_THEN `((-- &512505 + &770958 * sqrt2) * #1.2 pow 0 +  (-- &364208 - &1295359 * sqrt2) * #1.2 pow 1 +  (&1295359 + &585000 * sqrt2) * #1.2 pow 2 +  -- &585000 * #1.2 pow 3) = -- &2377941/ &25 + (& 294636* sqrt2)/ &5`  (fun t-> REWRITE_TAC[t]);
   REAL_ARITH_TAC;
   REWRITE_TAC[REAL_ARITH `a * b * c <= &0 <=> (&0 <= a * b * (-- c))`];
   STRIP_NN_TAC THEN (REPEAT (POP_ASSUM MP_TAC)) THEN TRY REAL_ARITH_TAC; 
   (* last subgoal *)
   REWRITE_TAC[real_interval;IN_ELIM_THM];
   REPEAT STRIP_TAC;
   EXISTS_TAC `x:real`;
   SUBGOAL_THEN `x <= &2 * h0` ASSUME_TAC;
   REWRITE_TAC[Sphere.h0];
   UNDISCH_TAC `x <= #1.26` THEN REAL_ARITH_TAC;
   ASM_SIMP_TAC [Sphere.h0;lmfun_lfun2];
   REPEAT (POP_ASSUM MP_TAC);
   REAL_ARITH_TAC;
  ]);;
  (* }}} *)

let hminus_prop = prove_by_refinement(
  `#1.2 <= hminus /\ hminus < #1.3 /\ marchal_quartic hminus = lmfun hminus`,
  (* {{{ proof *)
  [
  MP_TAC hminus_exists;
  MP_TAC Sphere.hminus;
  MESON_TAC[];
  ]);;
  (* }}} *)

let hminus_high = prove_by_refinement(
  `!h. (h0 <= h) ==> lmfun h = &0`,
  (* {{{ proof *)
  [
    REWRITE_TAC[Sphere.lmfun;REAL_ARITH `h0 <= h <=> (~(h <= h0) \/ (h = h0))`];
  GEN_TAC;
  DISCH_THEN DISJ_CASES_TAC THEN ASM_REWRITE_TAC[REAL_ARITH `h0 <= h0`;REAL_FIELD `(h0 - h0)/(h0 - &1) = &0`];
  ]);;
  (* }}} *)

let hminus_lt_h0 = prove_by_refinement(
  `!h. (&1 <= h) /\ (h < hplus) ==> (marchal_quartic h > &0)`,
  (* {{{ proof *)
  [
  REPEAT STRIP_TAC;
  REWRITE_TAC[Sphere.marchal_quartic];
  REWRITE_TAC[REAL_ARITH `(h-hplus)*(&9 * h pow 2 - &17 * h + &3)/u = (hplus - h)*(&17 * h - &9 * h pow 2 -  &3)/u /\ (u > &0 <=> &0 < u)`];
  MATCH_MP_TAC REAL_LT_MUL;
  CONJ_TAC;
  POP_ASSUM MP_TAC;
  SUBGOAL_THEN `&0 < sqrt(&2) - hplus` MP_TAC;
  REWRITE_TAC[Flyspeck_constants.bounds];
  REAL_ARITH_TAC;
MATCH_MP_TAC REAL_LT_MUL;
  CONJ_TAC;
  POP_ASSUM MP_TAC THEN REAL_ARITH_TAC;
  MATCH_MP_TAC REAL_LT_DIV;
  REWRITE_TAC[Flyspeck_constants.bounds];
  REPEAT (POP_ASSUM MP_TAC);
  REWRITE_TAC[Sphere.hplus];
  ABBREV_TAC `u = h - &1`;
  SUBGOAL_THEN `h = u + &1` (fun t->REWRITE_TAC[t]);
  POP_ASSUM MP_TAC THEN REAL_ARITH_TAC;
  REWRITE_TAC[REAL_FIELD `(&1 <= u + &1 <=> &0 <= u) /\ (u + &1 < #1.3254 <=> u < #0.3254) /\ (&17 * (u + &1) - &9 * (u+ &1) pow 2 - &3 = -- &9 * u pow 2 - u + &5)`];
  POP_ASSUM (fun t->ALL_TAC);
  REPEAT STRIP_TAC;
  SUBGOAL_THEN `u pow 2 < #0.3254 pow 2` MP_TAC;
  REWRITE_TAC[ GSYM REAL_LT_SQUARE_ABS];
  REPEAT (POP_ASSUM MP_TAC) THEN REAL_ARITH_TAC;
  REPEAT (POP_ASSUM MP_TAC);
  ABBREV_TAC `v = u pow 2`;
  REAL_ARITH_TAC;
  ]);;
  (* }}} *)

let hminus_lt_h0 = prove_by_refinement(
  `hminus < h0`,
  (* {{{ proof *)
  [
  REWRITE_TAC [REAL_ARITH `x < y <=> ~(x >= y)`];
  STRIP_TAC;
  MP_TAC hminus_prop;
  REPEAT STRIP_TAC;
  POP_ASSUM MP_TAC;
  MATCH_MP_TAC (REAL_ARITH `a > b ==>((a:real) = b ==> F)`);
  SUBGOAL_THEN `lmfun hminus = &0` (fun t -> REWRITE_TAC[t]);
  MATCH_MP_TAC hminus_high;
  UNDISCH_TAC `hminus >= h0` THEN REAL_ARITH_TAC;
  MATCH_MP_TAC   hminus_lt_h0;
  REPEAT (POP_ASSUM MP_TAC);
  REWRITE_TAC[Sphere.hplus;Sphere.h0];
  REAL_ARITH_TAC;
  ]);;
  (* }}} *)

let h0_lt_hplus = prove_by_refinement(
  `h0 < hplus`,
  (* {{{ proof *)
  [
  REWRITE_TAC[Sphere.h0;Sphere.hplus];
  REAL_ARITH_TAC;
  ]);;
  (* }}} *)

let hminus_gt = prove_by_refinement(
  `#1.2 <= hminus`,
  (* {{{ proof *)
  [
  REWRITE_TAC[hminus_prop];
  ]);;
  (* }}} *)

let lminus_ge_h0 = prove_by_refinement(
  `!h. (hplus <= h) /\ (h <= sqrt (&2)) ==> (marchal_quartic h <= &0)`,
  (* {{{ proof *)
  [
  REWRITE_TAC[Sphere.marchal_quartic];
  REWRITE_TAC[REAL_ARITH `a * b * c /d <= &0 <=> &0 <= a * b * (-- c)/d`];
  REPEAT STRIP_TAC;
  MATCH_MP_TAC REAL_LE_MUL;
  CONJ_TAC;
  POP_ASSUM MP_TAC;
  REAL_ARITH_TAC;
  MATCH_MP_TAC REAL_LE_MUL;
  CONJ_TAC;
  REPEAT (POP_ASSUM MP_TAC) THEN REAL_ARITH_TAC;
  MATCH_MP_TAC REAL_LE_DIV;
  CONJ_TAC;
  (* *)
    REPEAT (POP_ASSUM MP_TAC);
  REWRITE_TAC[Sphere.hplus];
  ABBREV_TAC `u = h - &1`;
  SUBGOAL_THEN `h = u + &1` (fun t->REWRITE_TAC[t]);
  POP_ASSUM MP_TAC THEN REAL_ARITH_TAC;
  REWRITE_TAC[REAL_FIELD `(&1 <= u + &1 <=> &0 <= u) /\ (#1.3254 <= u + &1 <=>  #0.3254 <= u) /\ (-- (&9* (u + &1) pow 2 - &17 * (u+ &1)  + &3) = -- &9 * u pow 2 - u + &5)`];
  REPEAT STRIP_TAC;
  SUBGOAL_THEN `u pow 2 < #0.42 pow 2` MP_TAC;
  REWRITE_TAC[ GSYM REAL_LT_SQUARE_ABS];
  POP_ASSUM MP_TAC;
  POP_ASSUM MP_TAC;
  EXPAND_TAC "u";
  SUBGOAL_THEN `sqrt(&2) - &1 < #0.42` (fun t-> (MP_TAC t) THEN REAL_ARITH_TAC);
  SUBGOAL_THEN `sqrt(&2) <  #1.414214` MP_TAC;
  REWRITE_TAC[GSYM Sphere.sqrt2;Flyspeck_constants.bounds];
  REAL_ARITH_TAC;
  SUBGOAL_THEN `u < #0.414214` MP_TAC;
  SUBGOAL_THEN `sqrt(&2) <  #1.414214` MP_TAC;
  REWRITE_TAC[GSYM Sphere.sqrt2;Flyspeck_constants.bounds];
  POP_ASSUM MP_TAC;
  REAL_ARITH_TAC;
  REAL_ARITH_TAC;
  MATCH_MP_TAC (REAL_ARITH `&0 < x ==> &0 <= x`);
  REWRITE_TAC[Flyspeck_constants.bounds];
  ]);;
  (* }}} *)

let gcy_high  = prove_by_refinement(
  `!y. (&2 * h0 <= y) ==> (gcy y = &4 * mm1/pi)`,
  (* {{{ proof *)
  [
  REWRITE_TAC[Sphere.gcy];
  SIMP_TAC[lmfun0];
  REAL_ARITH_TAC;
  ]);;
  (* }}} *)

let gcy_low = prove_by_refinement(
  `!y. (y <= &2 * h0) ==> (gcy y = gchi y)`,
  (* {{{ proof *)
  [
  REWRITE_TAC[Sphere.gcy;Sphere.gchi];
  SIMP_TAC[lmfun_lfun];
  REWRITE_TAC[Sphere.lfun;Sphere.h0];
  REPEAT STRIP_TAC;
  ABBREV_TAC `m = mm2/pi`;
  REAL_ARITH_TAC;
  ]);;
  (* }}} *)

let gcy_low_hminus = prove_by_refinement(
  `!y. (y <= &2 * hminus) ==> (gcy y = gchi y)`,
  (* {{{ proof *)
  [
  REPEAT STRIP_TAC;
  MATCH_MP_TAC gcy_low;
  POP_ASSUM MP_TAC;
  MP_TAC hminus_lt_h0;
  REAL_ARITH_TAC;
  ]);;
  (* }}} *)

let c2001 = prove_by_refinement(
  `!y. (y <= #2.001) ==> (y <= &2 * h0)`,
  [
    REWRITE_TAC[Sphere.h0];
  REAL_ARITH_TAC;
  ]
);;

let gcy_low_const = prove_by_refinement(
  `!y. (y <= #2.001) ==> (gcy y = gchi y)`,
  (* {{{ proof *)
  [
  REPEAT STRIP_TAC;
  MATCH_MP_TAC gcy_low;
  POP_ASSUM MP_TAC;
  REWRITE_TAC[Sphere.h0];
  REAL_ARITH_TAC;
  ]);;
  (* }}} *)

let gcy_high_hplus = prove_by_refinement(
  `!y. (&2 * hplus <= y) ==> (gcy y = &4 * mm1/pi)`,
  (* {{{ proof *)
  [
  REPEAT STRIP_TAC;
  MATCH_MP_TAC gcy_high;
  FIRST_X_ASSUM MP_TAC;
  MP_TAC h0_lt_hplus;
  REAL_ARITH_TAC;
  ]);;
  (* }}} *)

let hm0 = prove_by_refinement(
  `!y. ((y <= &2 * hminus) ==> (y <= &2 * h0))`,
  (* {{{ proof *)
  [
   MP_TAC hminus_lt_h0;
  CONV_TAC REAL_FIELD;
  ]);;
  (* }}} *)

let h0_lt_gt = prove_by_refinement(
  `((y <= #2.01) ==> (y <= &2 * h0)) /\
   ((#2.8 <= y) ==> (&2 * h0 <= y)) /\
  (( y <= &2) ==> (y <= &2 * h0)) /\
    ((sqrt8 <= y) ==> (&2 * h0 <= y)) /\   ((&2 * h0 <= y) ==> (&0 <= y)) /\
   ((&2 <= y) ==> (&0 <= y)) /\
   ((y <= &2 * hminus) ==> (y <= &2 * h0)) /\
   ((&2 * hminus <= y) ==> (&0 <= y))`,
  (* {{{ proof *)
  [
  REWRITE_TAC[Sphere.h0;sqrt8_sqrt2;hm0];
    MP_TAC sqrt2_lb;
  MP_TAC hminus_gt;
  REAL_ARITH_TAC;
  ]);;
  (* }}} *)

let sqrtxx = prove_by_refinement(
  `!x. &0 <= x ==> (sqrt(x * x) = x)`,
  (* {{{ proof *)
  [
    REWRITE_TAC[POW_2_SQRT_ABS;REAL_ARITH `x * x = x pow 2`];
  REAL_ARITH_TAC;
  ]);;
  (* }}} *)

let lmdih_ldih = prove_by_refinement(
  `!y1 y2 y3 y4 y5 y6. (&0 <= y1 /\ y1 <= &2 * h0) ==>(y_of_x lmdih_x_div_sqrtdelta_posbranch y1 y2 y3 y4 y5 y6 =  y_of_x ldih_x_div_sqrtdelta_posbranch y1 y2 y3 y4 y5 y6 )`,
  (* {{{ proof *)
  [
  REWRITE_TAC[Sphere.y_of_x;Nonlin_def.lmdih_x_div_sqrtdelta_posbranch;Nonlin_def.ldih_x_div_sqrtdelta_posbranch];
  MESON_TAC[sqrtxx;lmfun_lfun];
  ]);;
  (* }}} *)

let lmdih3_ldih3 = prove_by_refinement(
  `!y1 y2 y3 y4 y5 y6. (&0 <= y3 /\ y3 <= &2 * h0) ==>(y_of_x lmdih3_x_div_sqrtdelta_posbranch y1 y2 y3 y4 y5 y6 =  y_of_x ldih3_x_div_sqrtdelta_posbranch y1 y2 y3 y4 y5 y6 )`,
  (* {{{ proof *)
  [
  REWRITE_TAC[Sphere.rotate3;Nonlin_def.lmdih_x_div_sqrtdelta_posbranch;Nonlin_def.ldih_x_div_sqrtdelta_posbranch;Sphere.y_of_x;Nonlin_def.lmdih3_x_div_sqrtdelta_posbranch;Nonlin_def.ldih3_x_div_sqrtdelta_posbranch];
  MESON_TAC[sqrtxx;lmfun_lfun];
  ]);;
  (* }}} *)

let lmdih5_ldih5 = prove_by_refinement(
  `!y1 y2 y3 y4 y5 y6. (&0 <= y5 /\ y5 <= &2 * h0) ==>(y_of_x lmdih5_x_div_sqrtdelta_posbranch y1 y2 y3 y4 y5 y6 =  y_of_x ldih5_x_div_sqrtdelta_posbranch y1 y2 y3 y4 y5 y6 )`,
  (* {{{ proof *)
  [
  REWRITE_TAC[Sphere.rotate5;Nonlin_def.lmdih_x_div_sqrtdelta_posbranch;Nonlin_def.ldih_x_div_sqrtdelta_posbranch;Sphere.y_of_x;Nonlin_def.lmdih5_x_div_sqrtdelta_posbranch;Nonlin_def.ldih5_x_div_sqrtdelta_posbranch];
  MESON_TAC[sqrtxx;lmfun_lfun];
  ]);;
  (* }}} *)

let lmdih_ldih' = REWRITE_RULE[Sphere.y_of_x] lmdih_ldih;;
let lmdih3_ldih3' = REWRITE_RULE[Sphere.y_of_x] lmdih3_ldih3;;
let lmdih5_ldih5' = REWRITE_RULE[Sphere.y_of_x] lmdih5_ldih5;;

let lmdih0 = prove_by_refinement(
  `!y1 y2 y3 y4 y5 y6. (&2 * h0 <= y1 ) ==>(y_of_x lmdih_x_div_sqrtdelta_posbranch y1 y2 y3 y4 y5 y6 =  &0 )`,
  (* {{{ proof *)
  [
  REWRITE_TAC[Sphere.y_of_x;Nonlin_def.lmdih_x_div_sqrtdelta_posbranch;Nonlin_def.ldih_x_div_sqrtdelta_posbranch];
  MESON_TAC[sqrtxx;lmfun0;REAL_ARITH `&0 * x = &0 `;h0_lt_gt];
  ]);;
  (* }}} *)

let lmdih3_0 = prove_by_refinement(
  `!y1 y2 y3 y4 y5 y6. (&2 * h0 <= y3 ) ==>(y_of_x lmdih3_x_div_sqrtdelta_posbranch y1 y2 y3 y4 y5 y6 =  &0 )`,
  (* {{{ proof *)
  [
  REWRITE_TAC[Sphere.y_of_x;Nonlin_def.lmdih3_x_div_sqrtdelta_posbranch;Sphere.rotate3;Nonlin_def.lmdih_x_div_sqrtdelta_posbranch;Nonlin_def.ldih_x_div_sqrtdelta_posbranch];
  MESON_TAC[sqrtxx;lmfun0;REAL_ARITH `&0 * x = &0 `;h0_lt_gt];
  ]);;
  (* }}} *)

let lmdih5_0 = prove_by_refinement(
  `!y1 y2 y3 y4 y5 y6. (&2 * h0 <= y5 ) ==>(y_of_x lmdih5_x_div_sqrtdelta_posbranch y1 y2 y3 y4 y5 y6 =  &0 )`,
  (* {{{ proof *)
  [
  REWRITE_TAC[Sphere.y_of_x;Nonlin_def.lmdih5_x_div_sqrtdelta_posbranch;Sphere.rotate5;Nonlin_def.lmdih_x_div_sqrtdelta_posbranch;Nonlin_def.ldih_x_div_sqrtdelta_posbranch];
  MESON_TAC[sqrtxx;lmfun0;REAL_ARITH `&0 * x = &0 `;h0_lt_gt];
  ]);;
  (* }}} *)

let lmdih1_0' = REWRITE_RULE[Sphere.y_of_x] lmdih0;;
let lmdih3_0' = REWRITE_RULE[Sphere.y_of_x] lmdih3_0;;
let lmdih5_0' = REWRITE_RULE[Sphere.y_of_x] lmdih5_0;;

let vol4f_lmfun = prove_by_refinement(
  `! y1 y2 y3 y4 y5 y6. vol4f y1 y2 y3 y4 y5 y6 lmfun  = 
    (-- &8 * mm1)
    + 
    gcy y1 * dih_y y1 y2 y3 y4 y5 y6 + gcy y2 * dih2_y y1 y2 y3 y4 y5 y6 +
    gcy y3 * dih3_y y1 y2 y3 y4 y5 y6 +
    gcy y4 * dih4_y y1 y2 y3 y4 y5 y6 + 
      gcy y5 * dih5_y y1 y2 y3 y4 y5 y6 + 
    gcy y6 * dih6_y y1 y2 y3 y4 y5 y6`,
  (* {{{ proof *)
  [
    REWRITE_TAC[vol4f_palt;Sphere.gcy];
  REAL_ARITH_TAC;
  ]);;
  (* }}} *)

let gamma4fgcy_alt = prove_by_refinement(`gamma4fgcy y1 y2 y3 y4 y5 y6 lmfun =
    vol_y y1 y2 y3 y4 y5 y6 - 
    ((-- &8 * mm1)
    + 
    gcy y1 * dih_y y1 y2 y3 y4 y5 y6 + gcy y2 * dih2_y y1 y2 y3 y4 y5 y6 +
    gcy y3 * dih3_y y1 y2 y3 y4 y5 y6 +
    gcy y4 * dih4_y y1 y2 y3 y4 y5 y6 + 
      gcy y5 * dih5_y y1 y2 y3 y4 y5 y6 + 
    gcy y6 * dih6_y y1 y2 y3 y4 y5 y6)`,
[
REWRITE_TAC[Sphere.gamma4fgcy;Sphere.gamma4f;vol4f_lmfun];
]);;

let vol3f_palt = prove_by_refinement(
  `!y1 y2 y3 y4 y5 y6 r f. 
      (y3 = r) /\ (y4 = r) /\ (y5 = r) ==>
     (vol3f y1 y2 y6 r f =       (&2 * mm1 / pi) *
         (&2 * dih_y y1 y2 y3 y4 y5 y6 + &2 * dih2_y y1 y2 y3 y4 y5 y6 +
          &2 * dih6_y y1 y2 y3 y4 y5 y6 + dih3_y y1 y2 y3 y4 y5 y6 +
	  dih4_y y1 y2 y3 y4 y5 y6 + dih5_y y1 y2 y3 y4 y5 y6 - &3 * pi) -
         (&8 * mm2 / pi) *
         (f (y1 / &2) * dih_y y1 y2 y3 y4 y5 y6 +
          f (y2 / &2) * dih2_y y1 y2 y3 y4 y5 y6 +
          f (y6 / &2) * dih6_y y1 y2 y3 y4 y5 y6))`,
  (* {{{ proof *)
  [
  REPEAT GEN_TAC;
  REWRITE_TAC[Sphere.vol3f;Sphere.sol_y];
  DISCH_THEN (fun t->REWRITE_TAC[t]);
  REWRITE_TAC[Sphere.dih2_y;Sphere.dih3_y;Sphere.dih4_y;Sphere.dih5_y;Sphere.dih6_y];
  ABBREV_TAC `a = &8 * mm2/pi`;
  ABBREV_TAC `b = &2 * mm1/pi`;
  SUBGOAL_THEN `dih_y y2 r y1 r y6 r = dih_y y2 y1 r r r y6 /\ dih_y y2 y6 r r r y1 = dih_y y2 y1 r r r y6 /\ dih_y y6 r y2 r y1 r = dih_y y6 y1 r r r y2 /\  dih_y y1 r y6 r y2 r = dih_y y1 y2 r r r y6 /\ dih_y r y6 y1 y2 r r = dih_y r y1 y6 y2 r r` (fun t->REWRITE_TAC[t]);
 MESON_TAC[dih_y_sym2;dih_y_sym];
 MATCH_MP_TAC (REAL_ARITH `(a = a') ==> (a - c = a' - c)`);
 REAL_ARITH_TAC;
  ]);;
  (* }}} *)

let vol3f_135_palt = prove_by_refinement(
  `!y1 y2 y3 y4 y5 y6 r f. 
      (y2 = r) /\ (y4 = r) /\ (y6 = r) ==>
     (vol3f y1 y3 y5 r f =       (&2 * mm1 / pi) *
         (&2 * dih_y y1 y2 y3 y4 y5 y6 + &2 * dih3_y y1 y2 y3 y4 y5 y6 +
          &2 * dih5_y y1 y2 y3 y4 y5 y6 + dih2_y y1 y2 y3 y4 y5 y6 +
	  dih4_y y1 y2 y3 y4 y5 y6 + dih6_y y1 y2 y3 y4 y5 y6 - &3 * pi) -
         (&8 * mm2 / pi) *
         (f (y1 / &2) * dih_y y1 y2 y3 y4 y5 y6 +
          f (y3 / &2) * dih3_y y1 y2 y3 y4 y5 y6 +
          f (y5 / &2) * dih5_y y1 y2 y3 y4 y5 y6))`,
  (* {{{ proof *)
  [
  REPEAT GEN_TAC;
  REWRITE_TAC[Sphere.vol3f;Sphere.sol_y];
  DISCH_THEN (fun t->REWRITE_TAC[t]);
  REWRITE_TAC[Sphere.dih2_y;Sphere.dih3_y;Sphere.dih4_y;Sphere.dih5_y;Sphere.dih6_y];
  ABBREV_TAC `a = &8 * mm2/pi`;
  ABBREV_TAC `b = &2 * mm1/pi`;
  SUBGOAL_THEN `dih_y y3 r y1 r y5 r = dih_y y3 y5 r r r y1 /\ dih_y y5 r y3 r y1 r = dih_y y5 y1 r r r y3 /\ dih_y y1 r y5 r y3 r = dih_y y1 y3 r r r y5 /\ dih_y y1 r y3 r y5 r = dih_y y1 y3 r r r y5 /\ dih_y y3 y1 r r r y5 = dih_y y3 y5 r r r y1 /\ dih_y r r r y1 y5 y3 = dih_y r y3 y5 y1 r r /\ dih_y r y1 y5 y3 r r = dih_y r y5 y1 y3 r r` (fun t->REWRITE_TAC[t]);
 MESON_TAC[dih_y_sym2;dih_y_sym];
 MATCH_MP_TAC (REAL_ARITH `(a = a') ==> (a - c = a' - c)`);
 REAL_ARITH_TAC;
  ]);;
  (* }}} *)

let vol3r_alt = prove_by_refinement(
  `! y1 y2 y3 y4 y5 y6 r. (y3 = r ) /\ (y4 = r) /\ (y5 = r) ==> 
   (vol3r y1 y2 y6 r  = vol_y y1 y2 y3 y4 y5 y6)`,
  (* {{{ proof *)
  [
  REWRITE_TAC[Sphere.vol3r;Sphere.vol_y;Sphere.y_of_x;Sphere.vol_x;Sphere.delta_x];
  REPEAT GEN_TAC;
  DISCH_THEN (fun t->REWRITE_TAC[t]);
  AP_THM_TAC;
  AP_TERM_TAC;
  AP_TERM_TAC;
  REAL_ARITH_TAC;
  ]);;
  (* }}} *)

(* a few lemmas copied  from TameGeneral *)

let COS_PI3 = prove(`cos (pi / &3) = &1 / &2`,
   REWRITE_TAC[COS_SIN] THEN
     REWRITE_TAC[REAL_ARITH `pi / &2 - pi / &3 = pi / &6`] THEN
     REWRITE_TAC[SIN_PI6]);;

let ACS_2 = prove(`acs (&1 / &2) = pi / &3`,
   REWRITE_TAC[SYM COS_PI3] THEN
     MATCH_MP_TAC ACS_COS THEN
     REWRITE_TAC[REAL_ARITH `(&0 <= a / &3 <=> &0 <= a) /\ (a / &3 <= a <=> &0 <= a)`] THEN
     MATCH_MP_TAC REAL_LT_IMP_LE THEN
     REWRITE_TAC[PI_POS]);;

let sol0_POS = prove(`&0 < sol0`,
   REWRITE_TAC[Sphere.sol0] THEN
     REWRITE_TAC[REAL_ARITH `&0 < &3 * a - pi <=> pi / &3 < a`] THEN
     REWRITE_TAC[SYM ACS_2] THEN
     MATCH_MP_TAC ACS_MONO_LT THEN
     REAL_ARITH_TAC);;

let vol4f_alt = prove_by_refinement(
  `!y1 y2 y3 y4 y5 y6. vol4f y1 y2 y3 y4 y5 y6 lfun  = 
    (-- &8 * mm1)
    + 
    (&4 * mm1/pi - &8 * mm2 *(&1+const1)/(pi * const1) ) *
    (dih_y y1 y2 y3 y4 y5 y6 + dih2_y y1 y2 y3 y4 y5 y6 + dih3_y y1 y2 y3 y4 y5 y6 +
    dih4_y y1 y2 y3 y4 y5 y6 + dih5_y y1 y2 y3 y4 y5 y6 + dih6_y y1 y2 y3 y4 y5 y6)
    +
    (&8 * mm2 / (pi * const1)) * 
    (rhazim y1 y2 y3 y4 y5 y6 + rhazim2 y1 y2 y3 y4 y5 y6 + rhazim3 y1 y2 y3 y4 y5 y6 +
    rhazim4 y1 y2 y3 y4 y5 y6 + rhazim5 y1 y2 y3 y4 y5 y6 + rhazim6 y1 y2 y3 y4 y5 y6)
    `,
  (* {{{ proof *)
  [
  REPEAT GEN_TAC;
    REWRITE_TAC[vol4f_palt];
  REWRITE_TAC[Sphere.sol_y;lfun_ly;REAL_ARITH `&2 * y / &2 = y`;Sphere.rhazim;rhazim2_alt;rhazim3_alt;Sphere.rhazim4;Sphere.rhazim5;Sphere.rhazim6;];
  REWRITE_TAC[Sphere.rho;Sphere.node2_y;Sphere.node3_y;Sphere.rhazim];
  ONCE_REWRITE_TAC[REAL_ARITH `x  = y <=> x - y = &0`];
  SUBGOAL_THEN  `~(pi = &0)` ASSUME_TAC;
  SIMP_TAC[PI_POS;REAL_ARITH `&0 < x ==> ~(x= &0)`];
  SUBGOAL_THEN `~(const1 = &0)`   MP_TAC;
  REWRITE_TAC[GSYM sol0_over_pi_EQ_const1];
  MP_TAC sol0_POS;
  FIRST_X_ASSUM MP_TAC;
  CONV_TAC REAL_FIELD;
  FIRST_X_ASSUM MP_TAC;
  CONV_TAC REAL_FIELD;
  ]);;
  (* }}} *)

let vol2f_marchal_pow_y = prove_by_refinement(
  `!r y. vol2f y r marchal_quartic = 
       let fac = (-- (&8 * mm2/pi) * &2 * pi * inv ( #1.627 * (sqrt2 - &1)))  in 
    (&2 * mm1 / pi)  * &2 * pi 
     - (&2 * mm1 /pi) * &2 * pi * inv(r * &2) * y pow 1
   - fac * &3 * sqrt2 *hplus 
    + fac *(#1.5 * sqrt2 + #1.5 * hplus + #8.5 * sqrt2 * hplus) * y pow 1
   + fac * (-- #0.75  - #8.5 * sqrt2 * inv(&2) - #8.5 * hplus * inv(&2) - &9 * hplus * sqrt2 * inv (&4)) * y pow 2 
    + fac* ( #17.0 * inv (&8) + #9.0 * sqrt2 * inv(&8) + #9.0 * hplus * inv(&8)) * y pow 3 
    - fac * #9.0 * inv(&16) * y pow 4`,
  (* {{{ proof *)
  [
  REWRITE_TAC[Sphere.vol2f;Sphere.marchal_quartic;LET_DEF;LET_END_DEF;Sphere.hplus;GSYM Sphere.sqrt2;REAL_ARITH `(sqrt2 - &1) * &5 * (#1.3254 - &1) = #1.627 * (sqrt2- &1)`];
    REPEAT GEN_TAC;
  REWRITE_TAC[GSYM Sphere.hplus;real_div];
  REAL_ARITH_TAC;
  ]);;
  (* }}} *)

let vol2r_y = prove_by_refinement(
  `!y r. vol2r y r = &2 * pi * r * r * inv (&3) - #0.5 * pi * inv(&3) * y pow 2`,
  (* {{{ proof *)
  [
  REWRITE_TAC[Sphere.vol2r;real_div];
  REAL_ARITH_TAC;
  ]);;
  (* }}} *)

let ineq_expand6 = prove_by_refinement(
  `!a1 a2 a3 a4 a5 a6 b1 b2 b3 b4 b5 b6 x1 x2 x3 x4 x5 x6 P. (ineq [(a1,y1,b1);(a2,y2,b2);(a3,y3,b3);(a4,y4,b4);(a5,y5,b5);(a6,y6,b6)] P) <=> 
  (a1 <= y1 /\ y1 <= b1
     ==> a2 <= y2 /\ y2 <= b2
     ==> a3 <= y3 /\ y3 <= b3
     ==> a4 <= y4 /\ y4 <= b4
     ==> a5 <= y5 /\ y5 <= b5
     ==> a6 <= y6 /\ y6 <= b6
     ==> P)`,
  (* {{{ proof *)
  [
  REWRITE_TAC[ineq];
  ]);;
  (* }}} *)

let ineq_expand9 = prove_by_refinement(
  `!a1 a2 a3 a4 a5 a6 a7 a8 a9 b1 b2 b3 b4 b5 b6 b7 b8 b9 x1 x2 x3 x4 x5 x6 x7 x8 x9 P. (ineq [(a1,y1,b1);(a2,y2,b2);(a3,y3,b3);(a4,y4,b4);(a5,y5,b5);(a6,y6,b6);(a7,y7,b7);(a8,y8,b8);(a9,y9,b9)] P) <=> 
  (a1 <= y1 /\ y1 <= b1
     ==> a2 <= y2 /\ y2 <= b2
     ==> a3 <= y3 /\ y3 <= b3
     ==> a4 <= y4 /\ y4 <= b4
     ==> a5 <= y5 /\ y5 <= b5
     ==> a6 <= y6 /\ y6 <= b6
     ==> a7 <= y7 /\ y7 <= b7
     ==> a8 <= y8 /\ y8 <= b8
     ==> a9 <= y9 /\ y9 <= b9
     ==> P)`,
  (* {{{ proof *)
  [
  REWRITE_TAC[ineq];
  ]);;
  (* }}} *)

let pathL_bound = prove_by_refinement(
  `!y a.  FST (pathL a) <= y /\ y <= SND (pathL a) ==>
    FST a <= y /\ y <= SND a`,
  (* {{{ proof *)
  [
  REPEAT GEN_TAC;
    SUBGOAL_THEN `pathL a = pathL (FST a, SND a)` (fun t->ONCE_REWRITE_TAC[t]);
  REWRITE_TAC[PAIR];
    REWRITE_TAC[Sphere.pathL];
    REAL_ARITH_TAC;
  ]);;
  (* }}} *)

let pathR_bound = prove_by_refinement(
  `!y a. FST (pathR a) <= y /\ y <= SND (pathR a) ==>
    FST a <= y /\ y <= SND a`,
  (* {{{ proof *)
  [
  REPEAT GEN_TAC;
    SUBGOAL_THEN `pathR a = pathR (FST a, SND a)` (fun t->ONCE_REWRITE_TAC[t]);
  REWRITE_TAC[PAIR];
    REWRITE_TAC[Sphere.pathR];
    REAL_ARITH_TAC;
  ]);;
  (* }}} *)

let delta_x_eq0 = prove_by_refinement(
  `delta_x (&8) (&4) (&4) (&8) (&4) (&4) = &0 /\
   delta_x (&4) (&8) (&4) (&4) (&8) (&4) = &0
  `,
  (* {{{ proof *)
  [
  REWRITE_TAC[Sphere.delta_x];
    REAL_ARITH_TAC;
  ]);;
  (* }}} *)

let delta_x4_eq64 = prove_by_refinement(
  `delta_x4 (&8) (&4) (&4) (&8) (&4) (&4) = -- &64 /\
   delta_x4 (&4) (&8) (&4) (&4) (&8) (&4) = &64
  `,
  (* {{{ proof *)
  [
  REWRITE_TAC[Sphere.delta_x4];
    REAL_ARITH_TAC;
  ]);;
  (* }}} *)

let atn2_0y = prove_by_refinement(
  `atn2 (&0,&64) = pi / &2 /\ atn2 (&0 ,-- &64) = -- pi/ &2`,
  (* {{{ proof *)
  [
  REWRITE_TAC[Sphere.atn2;REAL_ARITH `~(abs(&64) < &0) /\ &0 < &64 /\ ~(abs (-- &64) < &0) /\ ~(&0 < -- &64) /\ (-- &64 < &0) /\ (&0 / &64 = &0) /\ (&0 / -- &64 = &0)`;ATN_0];
  REAL_ARITH_TAC;
  ]);;
  (* }}} *)

let OR_RULE rule1 rule2 th = try rule1 th with _ -> rule2 th;;

let rec REPEAT_RULE rule =
  fun t -> if can rule t then REPEAT_RULE rule (rule t) else t;;

let vol3f_lmln = prove_by_refinement(
  `!y1 y2 y3 y4 y5 y6. (y1 <= &2 * h0) /\ (y2 <= &2 * h0) /\ (y6 <= &2 * h0) ==> (vol3f y1 y2 y6 sqrt2 lmfun = vol3f y1 y2 y6 sqrt2 lfun)`,
  (* {{{ proof *)
  [
  REWRITE_TAC[Sphere.vol3f];
  MESON_TAC[lmfun_lfun];
  ]);;
  (* }}} *)

let vol3_vol_x = prove_by_refinement(
  `!x1 x2 x3 x4 x5 x6.  
   &0 <= x1 /\ &0 <= x2 /\ &0 <= x6 ==>
    (vol3_x_sqrt x1 x2 x3 x4 x5 x6 = vol_x x1 x2 (&2) (&2) (&2) x6)`,
  (* {{{ proof *)
  [
   MP_TAC (REAL_ARITH `&0 <= &2`);
  ONCE_REWRITE_TAC[REAL_ARITH `&0 = &0 pow 2`];
  REPEAT STRIP_TAC;
  REWRITE_TAC[vol3_x_sqrt;Sphere.vol_x;Sphere.vol_y;Sphere.y_of_x;Sphere.sqrt2];
  REPEAT AP_THM_TAC;
  AP_TERM_TAC ;
  AP_TERM_TAC;
  ASM_SIMP_TAC[SPEC `&0` sq_pow2];
  REWRITE_TAC[Sphere.delta_x];
  REAL_ARITH_TAC;
  ]);;
  (* }}} *)

let spec0 = SPECL [`&0`;`&0`;`&0`;`&0`;`&0`;`&0`;];;

let vol3f_x_lfun_alt = prove_by_refinement(
  `!x1 x2 x3 x4 x5 x6.  
   &0 <= x1 /\ &0 <= x2 /\ &0 <= x6 ==>
   vol3f_x_lfun x1 x2 x3 x4 x5 x6 =   (&2 * mm1 / pi) *
             (&2 * dih_x x1 x2 (&2) (&2) (&2) x6 +
              &2 * dih2_x x1 x2 (&2) (&2) (&2) x6 +
              &2 * dih6_x x1 x2 (&2) (&2) (&2) x6 +
              dih3_x x1 x2 (&2) (&2) (&2) x6 +
              dih4_x x1 x2 (&2) (&2) (&2) x6 +
              dih5_x x1 x2 (&2) (&2) (&2) x6 - &3 * pi) -
             (&8 * mm2 / pi) *
             ( ldih_x x1 x2 (&2) (&2) (&2) x6 +
              ldih2_x x1 x2 (&2) (&2) (&2) x6 +
              ldih6_x x1 x2 (&2) (&2) (&2) x6)`,
  (* {{{ proof *)
  [
   MP_TAC (REAL_ARITH `&0 <= &2`);
  ONCE_REWRITE_TAC[REAL_ARITH `&0 = &0 pow 2`];
  REPEAT STRIP_TAC;
  REWRITE_TAC[vol3f_x_lfun;];
  MP_TAC (SPECL [`sqrt x1`;`sqrt x2`;`sqrt2`;`sqrt2`;`sqrt2`;`sqrt x6`;`sqrt2`;`lfun`] vol3f_palt);
  REWRITE_TAC[];
  DISCH_THEN (fun t->REWRITE_TAC[t]);
  REWRITE_TAC[Sphere.sqrt2;Sphere.ldih_x;Sphere.ldih2_x;Sphere.ldih6_x;Sphere.dih4_x;Sphere.dih5_x;Sphere.dih6_x];
  ASM_SIMP_TAC[spec0 dih_x_y;spec0 dih2_x_y;spec0 dih3_x_y];
  ]);;
  (* }}} *)

let spech0 = SPECL [`&2 * h0`;`&0`;`&0`;`&0`;`&0`;`&0`;];;

let vol3f_x_sqrt2_lmplus_alt = prove_by_refinement(
  `!x1 x2 x3 x4 x5 x6.  
   ((&2 * h0) pow 2 <= x1) /\ &0 <= x2 /\ &0 <= x6 ==>
   vol3f_x_sqrt2_lmplus x1 x2 x3 x4 x5 x6 =   (&2 * mm1 / pi) *
             (&2 * dih_x x1 x2 (&2) (&2) (&2) x6 +
              &2 * dih2_x x1 x2 (&2) (&2) (&2) x6 +
              &2 * dih6_x x1 x2 (&2) (&2) (&2) x6 +
              dih3_x x1 x2 (&2) (&2) (&2) x6 +
              dih4_x x1 x2 (&2) (&2) (&2) x6 +
              dih5_x x1 x2 (&2) (&2) (&2) x6 - &3 * pi) -
             (&8 * mm2 / pi) *
             (
              ldih2_x x1 x2 (&2) (&2) (&2) x6 +
              ldih6_x x1 x2 (&2) (&2) (&2) x6)`,
  (* {{{ proof *)
  [
   MP_TAC (REAL_ARITH `&0 <= &2`);
  ONCE_REWRITE_TAC[REAL_ARITH `&0 = &0 pow 2`];
  REPEAT STRIP_TAC;
  REWRITE_TAC[vol3f_x_sqrt2_lmplus;vol3f_sqrt2_lmplus];
  REWRITE_TAC[Sphere.sqrt2;Sphere.ldih_x;Sphere.ldih2_x;Sphere.ldih6_x;Sphere.dih4_x;Sphere.dih5_x;Sphere.dih6_x];
  ASM_SIMP_TAC[spech0 dih_x_y;spech0 dih2_x_y;spech0 dih3_x_y];
  ]);;
  (* }}} *)

let vol3f_lm0 = prove_by_refinement(
  `!y1 y2 y3 y4 y5 y6. 
  ( &2 * h0 <= y1) /\ (y2 <= &2 * h0) /\ (y6 <= &2 * h0) ==>
  (vol3f y1 y2 y6 sqrt2 lmfun = vol3f_sqrt2_lmplus y1 y2 y3 y4 y5 y6)`,
  (* {{{ proof *)
  [
  REPEAT STRIP_TAC;
  REWRITE_TAC[vol3f_sqrt2_lmplus];
  MP_TAC (SPECL [`y1:real`;`y2:real`;`sqrt2`;`sqrt2`;`sqrt2`;`y6:real`; `sqrt2`;`lmfun`] vol3f_palt);
  REWRITE_TAC[];
  DISCH_THEN (fun t-> REWRITE_TAC[t]);
  ASM_SIMP_TAC[lmfun_lfun;lmfun0];
  REWRITE_TAC[REAL_ARITH `&0 * x = &0 /\ &0 + y = y`];
  ]);;
  (* }}} *)

let vol3r_126_x = prove_by_refinement(
  `vol3r (sqrt x1) (sqrt x2) (sqrt x6) sqrt2 = vol3_x_sqrt x1 x2 (x3:real) (x4:real) (x5:real) x6`,
  (* {{{ proof *)
  [
  REWRITE_TAC[vol3_x_sqrt;Sphere.vol3r;Sphere.vol_y];
  ]);;
  (* }}} *)

let num1_poly = prove_by_refinement(
  `! x1 x2 x3 x4 x5 x6. num1 x1 x2 x3 x4 x5 x6 =
     &64*x1*x4 - &32*x2*x4 - &32*x3*x4 - &4*x1*(x4 pow 2) - &32*x2*x5 + &32*x3*x5 + 
   &4*x2*x4*x5 + &32*x2*x6 - &32*x3*x6 + &4*x3*x4*x6 `,
  (* {{{ proof *)
  [
  REWRITE_TAC[Sphere.num1];
  REAL_ARITH_TAC;
  ]);;
  (* }}} *)

let ineq6_of_ineq5 = prove_by_refinement(
  `!a1 a2 a3 a4 a5 y1 y2 y3 y4 y5 b1 b2 b3 b4 b5 P. 
((!x1 x2 x3 x4 x5 x6. ineq[(a1,x1,b1);(a2,x2,b2);(a3,x3,b3);(a4,x4,b4);(a5,x5,b5);(&1,x6,&1)] 
   (P x1 x2 x3 x4 x5)) ==> 
    ineq[(a1,y1,b1);(a2,y2,b2);(a3,y3,b3);(a4,y4,b4);(a5,y5,b5)]
   (P y1 y2 y3 y4 y5))`,
  (* {{{ proof *)
  [
  REWRITE_TAC[ineq];
  MESON_TAC[REAL_ARITH `&1 <= &1`];
  ]);;
  (* }}} *)

let ineq6_of_ineq1 = prove_by_refinement(
  `!a1 y1 b1 P. 
((!x1 x2 x3 x4 x5 x6. ineq[(a1,x1,b1);(&1,x2,&1);(&1,x3,&1);(&1,x4,&1);(&1,x5,&1);(&1,x6,&1)] 
   (P x1)) ==> ineq[(a1,y1,b1)]  (P y1))`,
  (* {{{ proof *)
  [
  REWRITE_TAC[ineq];
  MESON_TAC[REAL_ARITH `&1 <= &1`];
  ]);;
  (* }}} *)

let dih_x_alt = prove_by_refinement(
  `!x1 x2 x3 x4 x5 x6. dih_x x1 x2 x3 x4 x5 x6 = pi / &2 +
     atn2
     (sqrt (&4 * x1 * delta_x x1 x2 x3 x4 x5 x6),
      --delta_x4 x1 x2 x3 x4 x5 x6) `,
  (* {{{ proof *)
  [
  REWRITE_TAC[Sphere.dih_x];
  REPEAT LET_TAC;
  REWRITE_TAC[];
  ]);;
  (* }}} *)

end;;