(* ========================================================================== *)
(* BCC LATTICE                                                                *)
(*                                                                            *)
(* Nonlinear Inequalities for BCC lattice optimality                          *)
(* Author: Thomas C. Hales                                                    *)
(* Date: 2012-03-20                                                           *)
(* ========================================================================== *)

(*
The inequalities in this file are not part of the Flyspeck project.

Nonlinear inequalities for the optimality of the BCC with respect
to surface areas when volumes of Voronoi cells are normalized to be 1.

Nonlinear inequalities for the optimality of the FCC with respect
to surface areas of Voronoi cells when the lattice packs a unit ball.

Reference: "Voronoi polyhedra of unit ball packings with small surface area"
Bezdek, KIss, Liu,  Periodica Mathematica Hungarica, Vol 39 (1-3), 1999, pp. 107--118.

The bottom of the file contains the Mathematica code used to prove
the local optimality of BCC and FCC, respectively.

To verify in C++, load the Flyspeck project, then this file,
then run the scripts in scripts.hl.

*)

module Bcc_lattice = struct

open Ineq;;

let bcc_value = new_definition `bcc_value = 
  (&3 / &2 + &3 * sqrt(&3))/ root 3 (&2)`;;
let  selling_volume2 = new_definition `selling_volume2 p01 p02 p03 p12 p13 p23 = 
  p01*p02*p03 + p01*p03*p12 + p02*p03*p12 +
  p01*p02*p13 + p02*p03*p13 + 
  p01*p12*p13 + p02*p12*p13 + p03*p12*p13 + 
  p01*p02*p23 + p01*p03*p23 + 
  p01*p12*p23 + p02*p12*p23 + p03*p12*p23 + 
  p01*p13*p23 + p02*p13*p23 + p03*p13*p23`;;
let  selling_surface_num = new_definition 
  ` selling_surface_num p01 p02 p03 p12 p13 p23 = 
  let  p0_123 = p01 + p02 + p03 in
  let  p1_023 = p01 + p12 + p13 in
  let  p2_013 = p02 + p12 + p23 in
  let  p3_012 = p03 + p13 + p23 in
  let  p01_23 = p02 + p03 + p12 + p13 in
  let  p02_13 = p01 + p03 + p12 + p23 in
  let  p03_12 = p01 + p02 + p13 + p23 in
  let  F01_23 = p01 * p23 * sqrt(p01_23) in
  let  F02_13 = p02 * p13 * sqrt(p02_13) in
  let  F03_12 = p03 * p12 * sqrt(p03_12) in
  let  F0_123 = (p12*p13+p12*p23+p13*p23)*sqrt(p0_123) in
  let  F1_023 = (p02*p03+p02*p23+p03*p23)*sqrt(p1_023) in
  let  F2_013 = (p01*p03+p01*p13+p03*p13)*sqrt(p2_013) in
  let  F3_012 = (p01*p02+p01*p12+p02*p12)*sqrt(p3_012) in
   &2*(F01_23+F02_13+F03_12+F0_123+F1_023+F2_013+F3_012)`;;
let selling_surface_nn = new_definition
  `selling_surface_nn p01 p02 p03 p12 p13 p23 = 
  selling_surface_num p01 p02 p03 p12 p13 p23 - bcc_value`;;
let sqrt01 = new_definition
  `sqrt01 t = t* sqrt(#0.1) / (#0.1)`;;
 
(* lower approx to sqrt on [0,0.1]. *)

(* replace sqrt with analytic sqrt01, near sqrt(0), to maintain analyticity *)

let  selling_surface_num2_013_approx = new_definition 
  ` selling_surface_num2_013_approx p01 p02 p03 p12 p13 p23 = 
  let  p0_123 = p01 + p02 + p03 in
  let  p1_023 = p01 + p12 + p13 in
  let  p2_013 = p02 + p12 + p23 in
  let  p3_012 = p03 + p13 + p23 in
  let  p01_23 = p02 + p03 + p12 + p13 in
  let  p02_13 = p01 + p03 + p12 + p23 in
  let  p03_12 = p01 + p02 + p13 + p23 in
  let  F01_23 = p01 * p23 * sqrt(p01_23) in
  let  F02_13 = p02 * p13 * sqrt(p02_13) in
  let  F03_12 = p03 * p12 * sqrt(p03_12) in
  let  F0_123 = (p12*p13+p12*p23+p13*p23)*sqrt(p0_123) in
  let  F1_023 = (p02*p03+p02*p23+p03*p23)*sqrt(p1_023) in
  let  F2_013 = (p01*p03+p01*p13+p03*p13)*sqrt01(p2_013) in
  let  F3_012 = (p01*p02+p01*p12+p02*p12)*sqrt(p3_012) in
   &2*(F01_23+F02_13+F03_12+F0_123+F1_023+F2_013+F3_012)`;;
let selling_surface_nn2_013 = new_definition
  `selling_surface_nn2_013 p01 p02 p03 p12 p13 p23 = 
  selling_surface_num2_013_approx p01 p02 p03 p12 p13 p23 - bcc_value`;;

let  selling_surface_num01_23_approx = new_definition 
  ` selling_surface_num01_23_approx p01 p02 p03 p12 p13 p23 = 
  let  p0_123 = p01 + p02 + p03 in
  let  p1_023 = p01 + p12 + p13 in
  let  p2_013 = p02 + p12 + p23 in
  let  p3_012 = p03 + p13 + p23 in
  let  p01_23 = p02 + p03 + p12 + p13 in
  let  p02_13 = p01 + p03 + p12 + p23 in
  let  p03_12 = p01 + p02 + p13 + p23 in
  let  F01_23 = p01 * p23 * sqrt01(p01_23) in
  let  F02_13 = p02 * p13 * sqrt(p02_13) in
  let  F03_12 = p03 * p12 * sqrt(p03_12) in
  let  F0_123 = (p12*p13+p12*p23+p13*p23)*sqrt(p0_123) in
  let  F1_023 = (p02*p03+p02*p23+p03*p23)*sqrt(p1_023) in
  let  F2_013 = (p01*p03+p01*p13+p03*p13)*sqrt(p2_013) in
  let  F3_012 = (p01*p02+p01*p12+p02*p12)*sqrt(p3_012) in
   &2*(F01_23+F02_13+F03_12+F0_123+F1_023+F2_013+F3_012)`;;
let selling_surface_nn01_23 = new_definition
  `selling_surface_nn01_23 p01 p02 p03 p12 p13 p23 = 
  selling_surface_num01_23_approx p01 p02 p03 p12 p13 p23 - bcc_value`;;
let selling_homog = new_definition
 `selling_homog y1 y2 y3 y4 y5 y6 = 
   (selling_surface_num y1 y2 y3 y4 y5 y6) -
      (bcc_value)*(root 6 (selling_volume2 y1 y2 y3 y4 y5 y6 pow 5))`;;
  

let fcc_ineq = new_definition
  `fcc_ineq y1 y2 y3 y4 y5 y6 = 
  selling_surface_num y1 y2 y3 y4 y5 y6 - 
    #12.0 * sqrt(&2) * sqrt(selling_volume2 y1 y2 y3 y4 y5 y6)`;;
let _ = 
  let v = map Parse_ineq.prep_autogen 
    [
      selling_volume2;
      bcc_value;
      selling_surface_num;
      selling_surface_nn;
      sqrt01;
      selling_surface_num2_013_approx;
      selling_surface_nn2_013;
      selling_surface_num01_23_approx;
      selling_surface_nn01_23;
      selling_homog;
      fcc_ineq;] in
    if (Lib.mem (hd v) !Parse_ineq.autogen) then () else 
      (Parse_ineq.autogen := !Parse_ineq.autogen @ v);;

  
let all_forall = Sphere.all_forall;;

Ineq.skip {
  idv="EIFIOKD";
  ineq = all_forall `ineq
   [
    (&0,y1,&81);
    (&0,y2,&81);
    (&0,y3,&81);
    (&0,y4,&81);
    (&0,y5,&81);
    (&0,y6,&81)
  ]
      ((selling_volume2 y1 y2 y3 y4 y5 y6 < &1) \/
   (selling_surface_num y1 y2 y3 y4 y5 y6 >= bcc_value)    )`;
  doc = "BCC cell main inequality. 
     This is the reference inequality, which is partitioned below into pieces.
     We must give special treatment at non-analytic points sqrt(0).
     BCC cell occurs when all variables are (16)^(-1/3) approx 0.39685.
   y1 = p01, y2 = p02, y3 = p03, y4 = p12, y5 = p13, y6 = p23";
  tags = [Cfsqp;Strongdodec];
};;

Ineq.add {
  idv="EIFIOKD-a";
  ineq = all_forall `ineq
   [
    (&5,y1,&81);
    (&0,y2,&81);
    (&0,y3,&81);
    (&0,y4,&81);
    (&0,y5,&81);
    (&0,y6,&81)
  ]
      ( (selling_surface_nn y1 y2 y3 y4 y5 y6 > &0) \/
	  (  y2 + y3 + y4 + y5 < #0.1) \/
	 (  y2 + y4 + y6  < #0.1) \/
	 ( y3 + y5 + y6  < #0.1) \/
	 (selling_volume2 y1 y2 y3 y4 y5 y6 < &1) 
      )`;
  doc = "BCC cell main inequality. 
    Region: some variable at least 5, moved to the first slot.
    p01_23, p2_013, p3_012 >= 0.1.
     We are bounded away from sqrt(0) on this domain.
     BCC cell occurs when all variables are (16)^(-1/3) approx 0.39685.
   y1 = p01, y2 = p02, y3 = p03, y4 = p12, y5 = p13, y6 = p23";
  tags = [Cfsqp;Strongdodec;Widthcutoff 0.05];
};;





Ineq.add {
  idv="EIFIOKD-b";
  ineq = all_forall `ineq
   [
    (&5,y1,&81);
    (&0,y2,#0.1);
    (&0,y3,#0.1);
    (&0,y4,#0.1);
    (&0,y5,#0.1);
    (&0,y6,&81)
  ]
      (    (selling_surface_nn01_23 y1 y2 y3 y4 y5 y6 > &0) \/
	     (  y2 + y3 + y4 + y5  > #0.1) \/
	 (selling_volume2 y1 y2 y3 y4 y5 y6 < &1) 
    )`;
  doc = "BCC cell main inequality. 
    Region: some variable at least 5, moved to the first slot.
    p01_23 <= 0.1.
     We are bounded away from sqrt(0) on this domain, except for sqrt(p01_23).
     BCC cell occurs when all variables are (16)^(-1/3) approx 0.39685.
   y1 = p01, y2 = p02, y3 = p03, y4 = p12, y5 = p13, y6 = p23";
  tags = [Cfsqp;Strongdodec];
};;


Ineq.add {
  idv="EIFIOKD-c";
  ineq = all_forall `ineq
   [
    (&5,y1,&81);
    (&0,y2,#0.1);
    (&0,y3,&81);
    (&0,y4,#0.1);
    (&0,y5,&81);
    (&0,y6,#0.1)
  ]
      (  (selling_surface_nn2_013 y1 y2 y3 y4 y5 y6 > &0)  \/ (  y2 + y4 + y6  > #0.1) \/
	 (selling_volume2 y1 y2 y3 y4 y5 y6 < &1) 
     )`;
  doc = "BCC cell main inequality. 
    Region: some variable at least 5, moved to the first slot.
    p2_013 <= 0.1 or p3_012 < 0.1, by symmetry assume p2_013 < 0.1.
     We are bounded away from sqrt(0) on this domain, except for sqrt(p01_23).
     BCC cell occurs when all variables are (16)^(-1/3) approx 0.39685.
   y1 = p01, y2 = p02, y3 = p03, y4 = p12, y5 = p13, y6 = p23";
  tags = [Cfsqp;Strongdodec];
};;

Ineq.add {
  idv="EIFIOKD1";
  ineq = all_forall `ineq
   [
    (#0.4,y1,&5);
    (&0,y2,&5);
    (&0,y3,&5);
    (&0,y4,&5);
    (&0,y5,&5);
    (&0,y6,&5)
  ]
      ( (selling_surface_nn y1 y2 y3 y4 y5 y6 > &0) \/
   ((selling_volume2 y1 y2 y3 y4 y5 y6 < &1)     ))`;
  doc = "BCC cell main inequality. 
     Region: all 0-5, some variable at least 0.4
     We are bounded away from sqrt(0) on this domain.
     BCC cell occurs when all variables are (16)^(-1/3) approx 0.39685.
   y1 = p01, y2 = p02, y3 = p03, y4 = p12, y5 = p13, y6 = p23";
  tags = [Cfsqp;Strongdodec];
};;

Ineq.add {
  idv="EIFIOKD2";
  ineq = all_forall `ineq
   [
    (&0,y1,#0.39);
    (&0,y2,#0.4);
    (&0,y3,#0.4);
    (&0,y4,#0.4);
    (&0,y5,#0.4);
    (&0,y6,#0.4)
  ]
      ((selling_volume2 y1 y2 y3 y4 y5 y6 < &1) \/
   (selling_surface_nn y1 y2 y3 y4 y5 y6 > &0)    )`;
  doc = "BCC cell main inequality. 
     Region: all 0-0.4, some variable at most 0.39, 
     We are bounded away from sqrt(0) on this domain.
     BCC cell occurs when all variables are (16)^(-1/3) approx 0.39685.
   y1 = p01, y2 = p02, y3 = p03, y4 = p12, y5 = p13, y6 = p23";
  tags = [Cfsqp;Strongdodec];
};;

Ineq.add {
  idv="EIFIOKD3";
  ineq = all_forall `ineq
   [
    (&1,y1,&1);
    (#0.975,y2,#0.999991);
    (#0.975,y3,&1);
    (#0.975,y4,&1);
    (#0.975,y5,&1);
    (#0.975,y6,&1)
  ]
   (selling_homog  y1 y2 y3 y4 y5 y6 >= &0)`;
  doc = "BCC cell main inequality. 
     Local analysis near a critical point.  
   Remove constraint by forming homogeneous (deg 15 expression in p)
    0.975 = 39/40  40/39 < 1.026.
    Homogeneity used to make the largest variable 1.  Symmetry shifted into first position.
    This is the case that an adjacent variable is at most 0.999991.
     We are bounded away from sqrt(0) on this domain.
     BCC cell occurs when all variables are (16)^(-1/3) approx 0.39685.
   y1 = p01, y2 = p02, y3 = p03, y4 = p12, y5 = p13, y6 = p23";
  tags = [Cfsqp;Strongdodec];
};;

Ineq.add {
  idv="EIFIOKD4";
  ineq = all_forall `ineq
   [
    (&1,y1,&1);
    (#0.975,y2,&1);
    (#0.975,y3,&1);
    (#0.975,y4,&1);
    (#0.975,y5,&1);
    (#0.975,y6,#0.999991)
  ]
   (selling_homog  y1 y2 y3 y4 y5 y6 >= &0)`;
  doc = "BCC cell main inequality. 
     Local analysis near a critical point.  
   Remove constraint by forming homogeneous (deg 15 expression in p)
    0.975 = 39/40  40/39 < 1.026.
    Homogeneity used to make the largest variable 1.  Symmetry shifted into first position.
    This is the case that an adjacent variable is at most 0.999991.
     We are bounded away from sqrt(0) on this domain.
     BCC cell occurs when all variables are (16)^(-1/3) approx 0.39685.
   y1 = p01, y2 = p02, y3 = p03, y4 = p12, y5 = p13, y6 = p23";
  tags = [Cfsqp;Strongdodec];
};;

Ineq.skip {
  idv="EIFIOKD5";
  ineq = all_forall `ineq
   [
    (&1,y1,&1);
    (#0.99999,y2,&1);
    (#0.99999,y3,&1);
    (#0.99999,y4,&1);
    (#0.99999,y5,&1);
    (#0.99999,y6,&1)
  ]
   (selling_homog  y1 y2 y3 y4 y5 y6 >= &0)`;
  doc = "BCC cell main inequality. 
     Local analysis near a critical point.  
   Remove constraint by forming homogeneous (deg 15 expression in p)
    Homogeneity used to make the largest variable 1.  Symmetry shifted into first position.
    This is that all variables are at least 0.99999.
   This has been done with interval arithmetic in Mathematica, by checking
    local optimality conditions in a nbd of (1,1,1,1,1,1).  (See code below.)
   Thus, we can skip this case in the C++ interval arithmetic verifications.
     We are bounded away from sqrt(0) on this domain.
     BCC cell occurs when all variables are (16)^(-1/3) approx 0.39685.
   y1 = p01, y2 = p02, y3 = p03, y4 = p12, y5 = p13, y6 = p23";
  tags = [Cfsqp;Strongdodec];
};;



Ineq.skip {
  idv="FXZXPNS";
  ineq = all_forall `ineq
   [
    (#1.3,y1,&30);
    (&1,y2,&30);
    (&0,y3,&30);
    (&0,y4,&30);
    (&0,y5,&30);
    (&0,y6,&30)
  ]
      ( (fcc_ineq y1 y2 y3 y4 y5 y6 >= &0) \/
	  (  y1 + y2 + y3  < &4) \/
	  (  y1 + y4 + y5  < &4) \/
	  (  y2 + y4 + y6  < &4) \/
	  (  y3 + y5 + y6  < &4) \/
	  (  y2 + y3 + y4 + y5 < &4) \/
	  (  y1 + y3 + y4  + y6 < &4) \/
	  (  y1 + y2 + y5 + y6  < &4)
      )`;
  doc = "FCC main inequality
    Region: assume p01 largest wlog, then p0X123 >= 4 ==> p01 >= 4/3.
    Assume wlog p02 >= p03, p12, p13, ==> p01X23 >=4 ==> p02 >= 1.
     We are bounded away from sqrt(0) on this domain.
     FCC cell occurs when {y1,..,y6} = {2,2,0,0,2,2};
    Local analysis at FCC is done in mathematica below.  
   This allows us to assume that some constraint (1)..(6) goes out to 4.04.
   y1 = p01, y2 = p02, y3 = p03, y4 = p12, y5 = p13, y6 = p23";
  tags = [Cfsqp;Strongdodec];
};;

Ineq.add {
  idv="FXZXPNS-1a";
  ineq = all_forall `ineq
   [
    (#1.3,y1,&30);
    (&1,y2,&30);
    (&1,y3,&30);
    (&0,y4,&30);
    (&0,y5,&30);
    (&0,y6,&30)
  ]
      ( 
	  (  y1 + y2 + y3  < #4.04) \/
	  (  y1 + y4 + y5  < &4) \/
	  (  y2 + y4 + y6  < &4) \/
	  (  y3 + y5 + y6  < &4) \/
	  (  y2 + y3 + y4 + y5 < &4) \/
	  (  y1 + y3 + y4  + y6 < &4) \/
	  (  y1 + y2 + y5 + y6  < &4) \/ (fcc_ineq y1 y2 y3 y4 y5 y6 > &0) 
      )`;
  doc = "FCC main inequality, case 1a.   a: 1<=y3
    Region: assume p01 largest wlog, then p0X123 >= 4 ==> p01 >= 4/3.
    Assume wlog p02 >= p03, p12, p13, ==> p01X23 >=4 ==> p02 >= 1.
     We are bounded away from sqrt(0) on this domain.
     FCC cell occurs when {y1,..,y6} = {2,2,0,0,2,2};
    Local analysis at FCC is done in mathematica below.  
   This allows us to assume that some constraint (1)..(6) goes out to 4.04.
   y1 = p01, y2 = p02, y3 = p03, y4 = p12, y5 = p13, y6 = p23";
  tags = [Cfsqp;Strongdodec;Cfsqp_branch 7];
};;


Ineq.add {
  idv="FXZXPNS-1b";
  ineq = all_forall `ineq
   [
    (#1.3,y1,&30);
    (&1,y2,&30);
    (&0,y3,&1);
    (&0,y4,&30);
    (&1,y5,&30);
    (&0,y6,&30)
  ]
      ( 
	  (  y1 + y2 + y3  < #4.04) \/
	  (  y1 + y4 + y5  < &4) \/
	  (  y2 + y4 + y6  < &4) \/
	  (  y3 + y5 + y6  < &4) \/
	  (  y2 + y3 + y4 + y5 < &4) \/
	  (  y1 + y3 + y4  + y6 < &4) \/
	  (  y1 + y2 + y5 + y6  < &4) \/ (fcc_ineq y1 y2 y3 y4 y5 y6 > &0) 
      )`;
  doc = "FCC main inequality, case 1b.   b: 1<=y5
    Region: assume p01 largest wlog, then p0X123 >= 4 ==> p01 >= 4/3.
    Assume wlog p02 >= p03, p12, p13, ==> p01X23 >=4 ==> p02 >= 1.
     We are bounded away from sqrt(0) on this domain.
     FCC cell occurs when {y1,..,y6} = {2,2,0,0,2,2};
    Local analysis at FCC is done in mathematica below.  
   This allows us to assume that some constraint (1)..(6) goes out to 4.04.
   y1 = p01, y2 = p02, y3 = p03, y4 = p12, y5 = p13, y6 = p23";
  tags = [Cfsqp;Strongdodec];
};;


Ineq.add {
  idv="FXZXPNS-1c";
  ineq = all_forall `ineq
   [
    (#1.3,y1,&30);
    (&1,y2,&30);
    (&0,y3,&1);
    (&0,y4,&30);
    (&0,y5,&1);
    (&1,y6,&30)
  ]
      ( 
	  (  y1 + y2 + y3  < #4.04) \/
	  (  y1 + y4 + y5  < &4) \/
	  (  y2 + y4 + y6  < &4) \/
	  (  y3 + y5 + y6  < &4) \/
	  (  y2 + y3 + y4 + y5 < &4) \/
	  (  y1 + y3 + y4  + y6 < &4) \/
	  (  y1 + y2 + y5 + y6  < &4) \/ (fcc_ineq y1 y2 y3 y4 y5 y6 > &0) 
      )`;
  doc = "FCC main inequality, case 1c.   c: 1<=y6.   
   (Note: if y3,y5,y6<1, then y3+y5+y6<4 and we are done.)
    Region: assume p01 largest wlog, then p0X123 >= 4 ==> p01 >= 4/3.
    Assume wlog p02 >= p03, p12, p13, ==> p01X23 >=4 ==> p02 >= 1.
     We are bounded away from sqrt(0) on this domain.
     FCC cell occurs when {y1,..,y6} = {2,2,0,0,2,2};
    Local analysis at FCC is done in mathematica below.  
   This allows us to assume that some constraint (1)..(6) goes out to 4.04.
   y1 = p01, y2 = p02, y3 = p03, y4 = p12, y5 = p13, y6 = p23";
  tags = [Cfsqp;Cfsqp_branch 7;Strongdodec];
};;




Ineq.add {
  idv="FXZXPNS-2";
  ineq = all_forall `ineq
   [
    (#1.3,y1,#4.04);
    (&1,y2,#4.04);
    (&0,y3,#4.04);
    (&0,y4,#4.04);
    (&0,y5,#4.04);
    (&0,y6,#4.04)
  ]
      ( 
	  (  y1 + y2 + y3  < &4) \/
	  (  y1 + y2 + y3  > #4.04) \/
	  (  y1 + y4 + y5  < #4.04) \/
	  (  y2 + y4 + y6  < &4) \/
	  (  y3 + y5 + y6  < &4) \/
	  (  y2 + y3 + y4 + y5 < &4) \/
	  (  y1 + y3 + y4  + y6 < &4) \/
	  (  y1 + y2 + y5 + y6  < &4) \/ (fcc_ineq y1 y2 y3 y4 y5 y6 > &0) 
      )`;
  doc = "FCC main inequality, case 2.  second constraint.
    Region: assume p01 largest wlog, then p0X123 >= 4 ==> p01 >= 4/3.
    Assume wlog p02 >= p03, p12, p13, ==> p01X23 >=4 ==> p02 >= 1.
     We are bounded away from sqrt(0) on this domain.
     FCC cell occurs when {y1,..,y6} = {2,2,0,0,2,2};
    Local analysis at FCC is done in mathematica below.  
   This allows us to assume that some constraint (1)..(6) goes out to 4.04.
   y1 = p01, y2 = p02, y3 = p03, y4 = p12, y5 = p13, y6 = p23";
  tags = [Cfsqp;Strongdodec];
};;


Ineq.add {
  idv="FXZXPNS-3";
  ineq = all_forall `ineq
   [
    (#1.3,y1,#4.04);
    (&1,y2,#4.04);
    (&0,y3,#4.04);
    (&0,y4,#4.04);
    (&0,y5,#4.04);
    (&0,y6,#4.04)
  ]
      ( 
	  (  y1 + y2 + y3  < &4) \/
	  (  y1 + y2 + y3  > #4.04) \/
	  (  y1 + y4 + y5  < &4) \/
	  (  y1 + y4 + y5  > #4.04) \/
	  (  y2 + y4 + y6  < #4.04) \/
	  (  y3 + y5 + y6  < &4) \/
	  (  y2 + y3 + y4 + y5 < &4) \/
	  (  y1 + y3 + y4  + y6 < &4) \/
	  (  y1 + y2 + y5 + y6  < &4) \/ (fcc_ineq y1 y2 y3 y4 y5 y6 > &0) 
      )`;
  doc = "FCC main inequality, case 3.  third constraint.
    Region: assume p01 largest wlog, then p0X123 >= 4 ==> p01 >= 4/3.
    Assume wlog p02 >= p03, p12, p13, ==> p01X23 >=4 ==> p02 >= 1.
     We are bounded away from sqrt(0) on this domain.
     FCC cell occurs when {y1,..,y6} = {2,2,0,0,2,2};
    Local analysis at FCC is done in mathematica below.  
   This allows us to assume that some constraint (1)..(6) goes out to 4.04.
   y1 = p01, y2 = p02, y3 = p03, y4 = p12, y5 = p13, y6 = p23";
  tags = [Cfsqp;Strongdodec];
};;


Ineq.add {
  idv="FXZXPNS-4";
  ineq = all_forall `ineq
   [
    (#1.3,y1,#4.04);
    (&1,y2,#4.04);
    (&0,y3,#4.04);
    (&0,y4,#4.04);
    (&0,y5,#4.04);
    (&0,y6,#4.04)
  ]
      ( 
	  (  y1 + y2 + y3  < &4) \/
	  (  y1 + y2 + y3  > #4.04) \/
	  (  y1 + y4 + y5  < &4) \/
	  (  y1 + y4 + y5  > #4.04) \/
	  (  y2 + y4 + y6  < &4) \/
	  (  y2 + y4 + y6  > #4.04) \/
	  (  y3 + y5 + y6  < #4.04) \/
	  (  y2 + y3 + y4 + y5 < &4) \/
	  (  y1 + y3 + y4  + y6 < &4) \/
	  (  y1 + y2 + y5 + y6  < &4) \/ (fcc_ineq y1 y2 y3 y4 y5 y6 > &0) 
      )`;
  doc = "FCC main inequality, case 4.  fourth constraint.
    Region: assume p01 largest wlog, then p0X123 >= 4 ==> p01 >= 4/3.
    Assume wlog p02 >= p03, p12, p13, ==> p01X23 >=4 ==> p02 >= 1.
     We are bounded away from sqrt(0) on this domain.
     FCC cell occurs when {y1,..,y6} = {2,2,0,0,2,2};
    Local analysis at FCC is done in mathematica below.  
   This allows us to assume that some constraint (1)..(6) goes out to 4.04.
   y1 = p01, y2 = p02, y3 = p03, y4 = p12, y5 = p13, y6 = p23";
  tags = [Cfsqp;Strongdodec];
};;


Ineq.add {
  idv="FXZXPNS-5";
  ineq = all_forall `ineq
   [
    (#1.3,y1,#4.04);
    (&1,y2,#4.04);
    (&0,y3,#4.04);
    (&0,y4,#4.04);
    (&0,y5,#4.04);
    (&0,y6,#4.04)
  ]
      ( 
	  (  y1 + y2 + y3  < &4) \/
	  (  y1 + y2 + y3  > #4.04) \/
	  (  y1 + y4 + y5  < &4) \/
	  (  y1 + y4 + y5  > #4.04) \/
	  (  y2 + y4 + y6  < &4) \/
	  (  y2 + y4 + y6  > #4.04) \/
	  (  y3 + y5 + y6  < &4) \/
	  (  y3 + y5 + y6  > #4.04) \/
	  (  y2 + y3 + y4 + y5 < #4.04) \/
	  (  y1 + y3 + y4  + y6 < &4) \/
	  (  y1 + y2 + y5 + y6  < &4) \/ (fcc_ineq y1 y2 y3 y4 y5 y6 > &0) 
      )`;
  doc = "FCC main inequality, case 5. 
    Region: assume p01 largest wlog, then p0X123 >= 4 ==> p01 >= 4/3.
    Assume wlog p02 >= p03, p12, p13, ==> p01X23 >=4 ==> p02 >= 1.
     We are bounded away from sqrt(0) on this domain.
     FCC cell occurs when {y1,..,y6} = {2,2,0,0,2,2};
    Local analysis at FCC is done in mathematica below.  
   This allows us to assume that some constraint (1)..(6) goes out to 4.04.
   y1 = p01, y2 = p02, y3 = p03, y4 = p12, y5 = p13, y6 = p23";
  tags = [Cfsqp;Strongdodec];
};;


Ineq.add {
  idv="FXZXPNS-6";
  ineq = all_forall `ineq
   [
    (#1.3,y1,#4.04);
    (&1,y2,#4.04);
    (&0,y3,#4.04);
    (&0,y4,#4.04);
    (&0,y5,#4.04);
    (&0,y6,#4.04)
  ]
      ( 
	  (  y1 + y2 + y3  < &4) \/
	  (  y1 + y2 + y3  > #4.04) \/
	  (  y1 + y4 + y5  < &4) \/
	  (  y1 + y4 + y5  > #4.04) \/
	  (  y2 + y4 + y6  < &4) \/
	  (  y2 + y4 + y6  > #4.04) \/
	  (  y3 + y5 + y6  < &4) \/
	  (  y3 + y5 + y6  > #4.04) \/
	  (  y2 + y3 + y4 + y5 < &4) \/
	  (  y2 + y3 + y4 + y5 > #4.04) \/
	  (  y1 + y3 + y4  + y6 < #4.04) \/
	  (  y1 + y2 + y5 + y6  < &4) \/ (fcc_ineq y1 y2 y3 y4 y5 y6 > &0) 
      )`;
  doc = "FCC main inequality, case 6. Final case. 
    Region: assume p01 largest wlog, then p0X123 >= 4 ==> p01 >= 4/3.
    Assume wlog p02 >= p03, p12, p13, ==> p01X23 >=4 ==> p02 >= 1.
     We are bounded away from sqrt(0) on this domain.
     FCC cell occurs when {y1,..,y6} = {2,2,0,0,2,2};
    Local analysis at FCC is done in mathematica below.  
   This allows us to assume that some constraint (1)..(6) goes out to 4.04.
   y1 = p01, y2 = p02, y3 = p03, y4 = p12, y5 = p13, y6 = p23";
  tags = [Cfsqp;Strongdodec];
};;



(* 

"Mathematica code used to prove local optimality in an explicit epsilon neighborhood
  of BCC";

"definition of objective function mu:";

SellingVolume2 = p01*p02*p03 + p01*p03*p12 + p02*p03*p12 + 
      p01*p02*p13 + p02*p03*p13 + p01*p12*p13 + p02*p12*p13 + 
      p03*p12*p13 + p01*p02*p23 + p01*p03*p23 + p01*p12*p23 + 
      p02*p12*p23 + p03*p12*p23 + p01*p13*p23 + p02*p13*p23 + p03*p13*p23;

p0X123 = p01 + p02 + p03;
p1X023 = p01 + p12 + p13;
p2X013 = p02 + p12 + p23;
p3X012 = p03 + p13 + p23;
p01X23 = p02 + p03 + p12 + p13;
p02X13 = p01 + p03 + p12 + p23;
p03X12 = p01 + p02 + p13 + p23;
F01X23 = p01*p23*Sqrt[p01X23];
F02X13 = p02*p13*Sqrt[p02X13];
F03X12 = p03*p12*Sqrt[p03X12];
F0X123 = (p12*p13 + p12*p23 + p13*p23)*Sqrt[p0X123];
F1X023 = (p02*p03 + p02*p23 + p03*p23)*Sqrt[p1X023];
F2X013 = (p01*p03 + p01*p13 + p03*p13)*Sqrt[p2X013];
F3X012 = (p01*p02 + p01*p12 + p02*p12)*Sqrt[p3X012];

SellingSurfaceNum = 2*(
    F01X23 + F02X13 + F03X12 + F0X123 + F1X023 + F2X013 + F3X012);

muSelling = ((SellingSurfaceNum/2)^3/(SellingVolume2)^(5/2));

"Interval Arithmetic Evaluation";
eval[r_, ii_] := (r /. {p01 -> 1, p02 -> ii, p03 -> ii, p12 -> 
    ii, p13 -> ii, p23 -> ii});
ii = Interval[{0.99999, 1}];

"Gradient";
dd1[i_] := D[muSelling1, {pv[[i]]}];
dds = Table[dd1[i], {i, 1, 5}];
{"Gradient", tab1i = Table[eval[dds[[i]], 1], {i, 1, 5}] // Simplify}

"Principal Minors";
pv = {p02, p03, p12, p13, p23};
dd[i_, j_] :=  D[muSelling1, {pv[[i]]}, {pv[[j]]}];
ddij = Table[eval[dd[i, j], ii], {i, 1, 5}, {j, 1, 5}];
minor[k_] := Det[Table[ddij[[i]][[j]], {i, 1, k}, {j, 1, k}]];
{"Table of Principal Minors", Table[minor[k], {k, 1, 5}]  }

*)

(*
(* Local Optimality of fcc lattice with respect to surface area of Voronoi
cell among those lattices containing a unit ball packing. *)
(* Corollary : Gauss on fcc lattice optimality.  *)

surfcc = SellingSurfaceNum/(SellingVolume2)^(1/2);
  (* at fcc, we get value 12 Sqrt[2], or about 16.9705627484771405856202646905 *)

epssub = {p01 -> 2 + eps01, p02 -> 2 + eps02, p03 -> 0 + 
    eps03, p12 -> 0 + eps12, p13 -> 2 + eps13, p23 -> 2 + eps23};
systemFcc = {t1 == p0X123 - 4, t2 == p1X023 - 4, t3 == p2X013 - 
          4, t4 == p3X012 - 4,
        t5 == p01X23 - 4, t6 == p02X13 - 4} /. epssub;
tsub = Solve[systemFcc, {eps01, eps02, eps03, eps12, eps13, eps23}];
surt = (surfcc /. epssub) /. tsub;
tvar = {t1, t2, t3, t4, t5, t6}; (* positive coordinates on constrained domain *)
grad = Table[D[surt, tvar[[i]]], {i, 1, 6}];
evalt[f_, a_] := f /. {t1 -> a, t2 -> a, t3 -> a, t4 -> a, t5 -> a, t6 -> a};
{"gradient is pos", evalt[grad, Interval[{0, 0.04}]]}
*)


end;;