projects_discrete_geom/bcc_lattice.hl

   1 (* ========================================================================== *)
   2 (* BCC LATTICE                                                                *)
   3 (*                                                                            *)
   4 (* Nonlinear Inequalities for BCC lattice optimality                          *)
   5 (* Author: Thomas C. Hales                                                    *)
   6 (* Date: 2012-03-20                                                           *)
   7 (* ========================================================================== *)
   8
   9 (*
  10 The inequalities in this file are not part of the Flyspeck project.
  11
  12 Nonlinear inequalities for the optimality of the BCC with respect
  13 to surface areas when volumes of Voronoi cells are normalized to be 1.
  14
  15 Nonlinear inequalities for the optimality of the FCC with respect
  16 to surface areas of Voronoi cells when the lattice packs a unit ball.
  17
  18 Reference: "Voronoi polyhedra of unit ball packings with small surface area"
  19 Bezdek, KIss, Liu,  Periodica Mathematica Hungarica, Vol 39 (1-3), 1999, pp. 107--118.
  20
  21 The bottom of the file contains the Mathematica code used to prove
  22 the local optimality of BCC and FCC, respectively.
  23
  24 To verify in C++, load the Flyspeck project, then this file,
  25 then run the scripts in scripts.hl.
  26
  27 *)
  28
  29 module Bcc_lattice = struct
  30
  31 open Ineq;;
  32
  33 let bcc_value = new_definition `bcc_value =
  34   (&3 / &2 + &3 * sqrt(&3))/ root 3 (&2)`;;
  35
  36 let  selling_volume2 = new_definition `selling_volume2 p01 p02 p03 p12 p13 p23 =
  37   p01*p02*p03 + p01*p03*p12 + p02*p03*p12 +
  38   p01*p02*p13 + p02*p03*p13 +
  39   p01*p12*p13 + p02*p12*p13 + p03*p12*p13 +
  40   p01*p02*p23 + p01*p03*p23 +
  41   p01*p12*p23 + p02*p12*p23 + p03*p12*p23 +
  42   p01*p13*p23 + p02*p13*p23 + p03*p13*p23`;;
  43
  44 let  selling_surface_num = new_definition
  45   ` selling_surface_num p01 p02 p03 p12 p13 p23 =
  46   let  p0_123 = p01 + p02 + p03 in
  47   let  p1_023 = p01 + p12 + p13 in
  48   let  p2_013 = p02 + p12 + p23 in
  49   let  p3_012 = p03 + p13 + p23 in
  50   let  p01_23 = p02 + p03 + p12 + p13 in
  51   let  p02_13 = p01 + p03 + p12 + p23 in
  52   let  p03_12 = p01 + p02 + p13 + p23 in
  53   let  F01_23 = p01 * p23 * sqrt(p01_23) in
  54   let  F02_13 = p02 * p13 * sqrt(p02_13) in
  55   let  F03_12 = p03 * p12 * sqrt(p03_12) in
  56   let  F0_123 = (p12*p13+p12*p23+p13*p23)*sqrt(p0_123) in
  57   let  F1_023 = (p02*p03+p02*p23+p03*p23)*sqrt(p1_023) in
  58   let  F2_013 = (p01*p03+p01*p13+p03*p13)*sqrt(p2_013) in
  59   let  F3_012 = (p01*p02+p01*p12+p02*p12)*sqrt(p3_012) in
  60    &2*(F01_23+F02_13+F03_12+F0_123+F1_023+F2_013+F3_012)`;;
  61
  62 let selling_surface_nn = new_definition
  63   `selling_surface_nn p01 p02 p03 p12 p13 p23 =
  64   selling_surface_num p01 p02 p03 p12 p13 p23 - bcc_value`;;
  65
  66 let sqrt01 = new_definition
  67   `sqrt01 t = t* sqrt(#0.1) / (#0.1)`;;
  68 (* lower approx to sqrt on [0,0.1]. *)
  69
  70 (* replace sqrt with analytic sqrt01, near sqrt(0), to maintain analyticity *)
  71
  72 let  selling_surface_num2_013_approx = new_definition
  73   ` selling_surface_num2_013_approx p01 p02 p03 p12 p13 p23 =
  74   let  p0_123 = p01 + p02 + p03 in
  75   let  p1_023 = p01 + p12 + p13 in
  76   let  p2_013 = p02 + p12 + p23 in
  77   let  p3_012 = p03 + p13 + p23 in
  78   let  p01_23 = p02 + p03 + p12 + p13 in
  79   let  p02_13 = p01 + p03 + p12 + p23 in
  80   let  p03_12 = p01 + p02 + p13 + p23 in
  81   let  F01_23 = p01 * p23 * sqrt(p01_23) in
  82   let  F02_13 = p02 * p13 * sqrt(p02_13) in
  83   let  F03_12 = p03 * p12 * sqrt(p03_12) in
  84   let  F0_123 = (p12*p13+p12*p23+p13*p23)*sqrt(p0_123) in
  85   let  F1_023 = (p02*p03+p02*p23+p03*p23)*sqrt(p1_023) in
  86   let  F2_013 = (p01*p03+p01*p13+p03*p13)*sqrt01(p2_013) in
  87   let  F3_012 = (p01*p02+p01*p12+p02*p12)*sqrt(p3_012) in
  88    &2*(F01_23+F02_13+F03_12+F0_123+F1_023+F2_013+F3_012)`;;
  89
  90 let selling_surface_nn2_013 = new_definition
  91   `selling_surface_nn2_013 p01 p02 p03 p12 p13 p23 =
  92   selling_surface_num2_013_approx p01 p02 p03 p12 p13 p23 - bcc_value`;;
  93
  94
  95 let  selling_surface_num01_23_approx = new_definition
  96   ` selling_surface_num01_23_approx p01 p02 p03 p12 p13 p23 =
  97   let  p0_123 = p01 + p02 + p03 in
  98   let  p1_023 = p01 + p12 + p13 in
  99   let  p2_013 = p02 + p12 + p23 in
 100   let  p3_012 = p03 + p13 + p23 in
 101   let  p01_23 = p02 + p03 + p12 + p13 in
 102   let  p02_13 = p01 + p03 + p12 + p23 in
 103   let  p03_12 = p01 + p02 + p13 + p23 in
 104   let  F01_23 = p01 * p23 * sqrt01(p01_23) in
 105   let  F02_13 = p02 * p13 * sqrt(p02_13) in
 106   let  F03_12 = p03 * p12 * sqrt(p03_12) in
 107   let  F0_123 = (p12*p13+p12*p23+p13*p23)*sqrt(p0_123) in
 108   let  F1_023 = (p02*p03+p02*p23+p03*p23)*sqrt(p1_023) in
 109   let  F2_013 = (p01*p03+p01*p13+p03*p13)*sqrt(p2_013) in
 110   let  F3_012 = (p01*p02+p01*p12+p02*p12)*sqrt(p3_012) in
 111    &2*(F01_23+F02_13+F03_12+F0_123+F1_023+F2_013+F3_012)`;;
 112
 113 let selling_surface_nn01_23 = new_definition
 114   `selling_surface_nn01_23 p01 p02 p03 p12 p13 p23 =
 115   selling_surface_num01_23_approx p01 p02 p03 p12 p13 p23 - bcc_value`;;
 116
 117 let selling_homog = new_definition
 118  `selling_homog y1 y2 y3 y4 y5 y6 =
 119    (selling_surface_num y1 y2 y3 y4 y5 y6) -
 120       (bcc_value)*(root 6 (selling_volume2 y1 y2 y3 y4 y5 y6 pow 5))`;;
 121
 122 let fcc_ineq = new_definition
 123   `fcc_ineq y1 y2 y3 y4 y5 y6 =
 124   selling_surface_num y1 y2 y3 y4 y5 y6 -
 125     #12.0 * sqrt(&2) * sqrt(selling_volume2 y1 y2 y3 y4 y5 y6)`;;
 126
 127 let _ =
 128   let v = map Parse_ineq.prep_autogen
 129     [
 130       selling_volume2;
 131       bcc_value;
 132       selling_surface_num;
 133       selling_surface_nn;
 134       sqrt01;
 135       selling_surface_num2_013_approx;
 136       selling_surface_nn2_013;
 137       selling_surface_num01_23_approx;
 138       selling_surface_nn01_23;
 139       selling_homog;
 140       fcc_ineq;] in
 141     if (Lib.mem (hd v) !Parse_ineq.autogen) then () else
 142       (Parse_ineq.autogen := !Parse_ineq.autogen @ v);;
 143
 144
 145 let all_forall = Sphere.all_forall;;
 146
 147 Ineq.skip {
 148   idv="EIFIOKD";
 149   ineq = all_forall `ineq
 150    [
 151     (&0,y1,&81);
 152     (&0,y2,&81);
 153     (&0,y3,&81);
 154     (&0,y4,&81);
 155     (&0,y5,&81);
 156     (&0,y6,&81)
 157   ]
 158       ((selling_volume2 y1 y2 y3 y4 y5 y6 < &1) \/
 159    (selling_surface_num y1 y2 y3 y4 y5 y6 >= bcc_value)    )`;
 160   doc = "BCC cell main inequality.
 161      This is the reference inequality, which is partitioned below into pieces.
 162      We must give special treatment at non-analytic points sqrt(0).
 163      BCC cell occurs when all variables are (16)^(-1/3) approx 0.39685.
 164    y1 = p01, y2 = p02, y3 = p03, y4 = p12, y5 = p13, y6 = p23";
 165   tags = [Cfsqp;Strongdodec];
 166 };;
 167
 168 Ineq.add {
 169   idv="EIFIOKD-a";
 170   ineq = all_forall `ineq
 171    [
 172     (&5,y1,&81);
 173     (&0,y2,&81);
 174     (&0,y3,&81);
 175     (&0,y4,&81);
 176     (&0,y5,&81);
 177     (&0,y6,&81)
 178   ]
 179       ( (selling_surface_nn y1 y2 y3 y4 y5 y6 > &0) \/
 180           (  y2 + y3 + y4 + y5 < #0.1) \/
 181          (  y2 + y4 + y6  < #0.1) \/
 182          ( y3 + y5 + y6  < #0.1) \/
 183          (selling_volume2 y1 y2 y3 y4 y5 y6 < &1)
 184       )`;
 185   doc = "BCC cell main inequality.
 186     Region: some variable at least 5, moved to the first slot.
 187     p01_23, p2_013, p3_012 >= 0.1.
 188      We are bounded away from sqrt(0) on this domain.
 189      BCC cell occurs when all variables are (16)^(-1/3) approx 0.39685.
 190    y1 = p01, y2 = p02, y3 = p03, y4 = p12, y5 = p13, y6 = p23";
 191   tags = [Cfsqp;Strongdodec;Widthcutoff 0.05];
 192 };;
 193
 194
 195
 196
 197
 198 Ineq.add {
 199   idv="EIFIOKD-b";
 200   ineq = all_forall `ineq
 201    [
 202     (&5,y1,&81);
 203     (&0,y2,#0.1);
 204     (&0,y3,#0.1);
 205     (&0,y4,#0.1);
 206     (&0,y5,#0.1);
 207     (&0,y6,&81)
 208   ]
 209       (    (selling_surface_nn01_23 y1 y2 y3 y4 y5 y6 > &0) \/
 210              (  y2 + y3 + y4 + y5  > #0.1) \/
 211          (selling_volume2 y1 y2 y3 y4 y5 y6 < &1)
 212     )`;
 213   doc = "BCC cell main inequality.
 214     Region: some variable at least 5, moved to the first slot.
 215     p01_23 <= 0.1.
 216      We are bounded away from sqrt(0) on this domain, except for sqrt(p01_23).
 217      BCC cell occurs when all variables are (16)^(-1/3) approx 0.39685.
 218    y1 = p01, y2 = p02, y3 = p03, y4 = p12, y5 = p13, y6 = p23";
 219   tags = [Cfsqp;Strongdodec];
 220 };;
 221
 222
 223 Ineq.add {
 224   idv="EIFIOKD-c";
 225   ineq = all_forall `ineq
 226    [
 227     (&5,y1,&81);
 228     (&0,y2,#0.1);
 229     (&0,y3,&81);
 230     (&0,y4,#0.1);
 231     (&0,y5,&81);
 232     (&0,y6,#0.1)
 233   ]
 234       (  (selling_surface_nn2_013 y1 y2 y3 y4 y5 y6 > &0)  \/ (  y2 + y4 + y6  > #0.1) \/
 235          (selling_volume2 y1 y2 y3 y4 y5 y6 < &1)
 236      )`;
 237   doc = "BCC cell main inequality.
 238     Region: some variable at least 5, moved to the first slot.
 239     p2_013 <= 0.1 or p3_012 < 0.1, by symmetry assume p2_013 < 0.1.
 240      We are bounded away from sqrt(0) on this domain, except for sqrt(p01_23).
 241      BCC cell occurs when all variables are (16)^(-1/3) approx 0.39685.
 242    y1 = p01, y2 = p02, y3 = p03, y4 = p12, y5 = p13, y6 = p23";
 243   tags = [Cfsqp;Strongdodec];
 244 };;
 245
 246 Ineq.add {
 247   idv="EIFIOKD1";
 248   ineq = all_forall `ineq
 249    [
 250     (#0.4,y1,&5);
 251     (&0,y2,&5);
 252     (&0,y3,&5);
 253     (&0,y4,&5);
 254     (&0,y5,&5);
 255     (&0,y6,&5)
 256   ]
 257       ( (selling_surface_nn y1 y2 y3 y4 y5 y6 > &0) \/
 258    ((selling_volume2 y1 y2 y3 y4 y5 y6 < &1)     ))`;
 259   doc = "BCC cell main inequality.
 260      Region: all 0-5, some variable at least 0.4
 261      We are bounded away from sqrt(0) on this domain.
 262      BCC cell occurs when all variables are (16)^(-1/3) approx 0.39685.
 263    y1 = p01, y2 = p02, y3 = p03, y4 = p12, y5 = p13, y6 = p23";
 264   tags = [Cfsqp;Strongdodec];
 265 };;
 266
 267 Ineq.add {
 268   idv="EIFIOKD2";
 269   ineq = all_forall `ineq
 270    [
 271     (&0,y1,#0.39);
 272     (&0,y2,#0.4);
 273     (&0,y3,#0.4);
 274     (&0,y4,#0.4);
 275     (&0,y5,#0.4);
 276     (&0,y6,#0.4)
 277   ]
 278       ((selling_volume2 y1 y2 y3 y4 y5 y6 < &1) \/
 279    (selling_surface_nn y1 y2 y3 y4 y5 y6 > &0)    )`;
 280   doc = "BCC cell main inequality.
 281      Region: all 0-0.4, some variable at most 0.39,
 282      We are bounded away from sqrt(0) on this domain.
 283      BCC cell occurs when all variables are (16)^(-1/3) approx 0.39685.
 284    y1 = p01, y2 = p02, y3 = p03, y4 = p12, y5 = p13, y6 = p23";
 285   tags = [Cfsqp;Strongdodec];
 286 };;
 287
 288 Ineq.add {
 289   idv="EIFIOKD3";
 290   ineq = all_forall `ineq
 291    [
 292     (&1,y1,&1);
 293     (#0.975,y2,#0.999991);
 294     (#0.975,y3,&1);
 295     (#0.975,y4,&1);
 296     (#0.975,y5,&1);
 297     (#0.975,y6,&1)
 298   ]
 299    (selling_homog  y1 y2 y3 y4 y5 y6 >= &0)`;
 300   doc = "BCC cell main inequality.
 301      Local analysis near a critical point.
 302    Remove constraint by forming homogeneous (deg 15 expression in p)
 303     0.975 = 39/40  40/39 < 1.026.
 304     Homogeneity used to make the largest variable 1.  Symmetry shifted into first position.
 305     This is the case that an adjacent variable is at most 0.999991.
 306      We are bounded away from sqrt(0) on this domain.
 307      BCC cell occurs when all variables are (16)^(-1/3) approx 0.39685.
 308    y1 = p01, y2 = p02, y3 = p03, y4 = p12, y5 = p13, y6 = p23";
 309   tags = [Cfsqp;Strongdodec];
 310 };;
 311
 312 Ineq.add {
 313   idv="EIFIOKD4";
 314   ineq = all_forall `ineq
 315    [
 316     (&1,y1,&1);
 317     (#0.975,y2,&1);
 318     (#0.975,y3,&1);
 319     (#0.975,y4,&1);
 320     (#0.975,y5,&1);
 321     (#0.975,y6,#0.999991)
 322   ]
 323    (selling_homog  y1 y2 y3 y4 y5 y6 >= &0)`;
 324   doc = "BCC cell main inequality.
 325      Local analysis near a critical point.
 326    Remove constraint by forming homogeneous (deg 15 expression in p)
 327     0.975 = 39/40  40/39 < 1.026.
 328     Homogeneity used to make the largest variable 1.  Symmetry shifted into first position.
 329     This is the case that an adjacent variable is at most 0.999991.
 330      We are bounded away from sqrt(0) on this domain.
 331      BCC cell occurs when all variables are (16)^(-1/3) approx 0.39685.
 332    y1 = p01, y2 = p02, y3 = p03, y4 = p12, y5 = p13, y6 = p23";
 333   tags = [Cfsqp;Strongdodec];
 334 };;
 335
 336 Ineq.skip {
 337   idv="EIFIOKD5";
 338   ineq = all_forall `ineq
 339    [
 340     (&1,y1,&1);
 341     (#0.99999,y2,&1);
 342     (#0.99999,y3,&1);
 343     (#0.99999,y4,&1);
 344     (#0.99999,y5,&1);
 345     (#0.99999,y6,&1)
 346   ]
 347    (selling_homog  y1 y2 y3 y4 y5 y6 >= &0)`;
 348   doc = "BCC cell main inequality.
 349      Local analysis near a critical point.
 350    Remove constraint by forming homogeneous (deg 15 expression in p)
 351     Homogeneity used to make the largest variable 1.  Symmetry shifted into first position.
 352     This is that all variables are at least 0.99999.
 353    This has been done with interval arithmetic in Mathematica, by checking
 354     local optimality conditions in a nbd of (1,1,1,1,1,1).  (See code below.)
 355    Thus, we can skip this case in the C++ interval arithmetic verifications.
 356      We are bounded away from sqrt(0) on this domain.
 357      BCC cell occurs when all variables are (16)^(-1/3) approx 0.39685.
 358    y1 = p01, y2 = p02, y3 = p03, y4 = p12, y5 = p13, y6 = p23";
 359   tags = [Cfsqp;Strongdodec];
 360 };;
 361
 362
 363
 364 Ineq.skip {
 365   idv="FXZXPNS";
 366   ineq = all_forall `ineq
 367    [
 368     (#1.3,y1,&30);
 369     (&1,y2,&30);
 370     (&0,y3,&30);
 371     (&0,y4,&30);
 372     (&0,y5,&30);
 373     (&0,y6,&30)
 374   ]
 375       ( (fcc_ineq y1 y2 y3 y4 y5 y6 >= &0) \/
 376           (  y1 + y2 + y3  < &4) \/
 377           (  y1 + y4 + y5  < &4) \/
 378           (  y2 + y4 + y6  < &4) \/
 379           (  y3 + y5 + y6  < &4) \/
 380           (  y2 + y3 + y4 + y5 < &4) \/
 381           (  y1 + y3 + y4  + y6 < &4) \/
 382           (  y1 + y2 + y5 + y6  < &4)
 383       )`;
 384   doc = "FCC main inequality
 385     Region: assume p01 largest wlog, then p0X123 >= 4 ==> p01 >= 4/3.
 386     Assume wlog p02 >= p03, p12, p13, ==> p01X23 >=4 ==> p02 >= 1.
 387      We are bounded away from sqrt(0) on this domain.
 388      FCC cell occurs when {y1,..,y6} = {2,2,0,0,2,2};
 389     Local analysis at FCC is done in mathematica below.
 390    This allows us to assume that some constraint (1)..(6) goes out to 4.04.
 391    y1 = p01, y2 = p02, y3 = p03, y4 = p12, y5 = p13, y6 = p23";
 392   tags = [Cfsqp;Strongdodec];
 393 };;
 394
 395 Ineq.add {
 396   idv="FXZXPNS-1a";
 397   ineq = all_forall `ineq
 398    [
 399     (#1.3,y1,&30);
 400     (&1,y2,&30);
 401     (&1,y3,&30);
 402     (&0,y4,&30);
 403     (&0,y5,&30);
 404     (&0,y6,&30)
 405   ]
 406       (
 407           (  y1 + y2 + y3  < #4.04) \/
 408           (  y1 + y4 + y5  < &4) \/
 409           (  y2 + y4 + y6  < &4) \/
 410           (  y3 + y5 + y6  < &4) \/
 411           (  y2 + y3 + y4 + y5 < &4) \/
 412           (  y1 + y3 + y4  + y6 < &4) \/
 413           (  y1 + y2 + y5 + y6  < &4) \/ (fcc_ineq y1 y2 y3 y4 y5 y6 > &0)
 414       )`;
 415   doc = "FCC main inequality, case 1a.   a: 1<=y3
 416     Region: assume p01 largest wlog, then p0X123 >= 4 ==> p01 >= 4/3.
 417     Assume wlog p02 >= p03, p12, p13, ==> p01X23 >=4 ==> p02 >= 1.
 418      We are bounded away from sqrt(0) on this domain.
 419      FCC cell occurs when {y1,..,y6} = {2,2,0,0,2,2};
 420     Local analysis at FCC is done in mathematica below.
 421    This allows us to assume that some constraint (1)..(6) goes out to 4.04.
 422    y1 = p01, y2 = p02, y3 = p03, y4 = p12, y5 = p13, y6 = p23";
 423   tags = [Cfsqp;Strongdodec;Cfsqp_branch 7];
 424 };;
 425
 426
 427 Ineq.add {
 428   idv="FXZXPNS-1b";
 429   ineq = all_forall `ineq
 430    [
 431     (#1.3,y1,&30);
 432     (&1,y2,&30);
 433     (&0,y3,&1);
 434     (&0,y4,&30);
 435     (&1,y5,&30);
 436     (&0,y6,&30)
 437   ]
 438       (
 439           (  y1 + y2 + y3  < #4.04) \/
 440           (  y1 + y4 + y5  < &4) \/
 441           (  y2 + y4 + y6  < &4) \/
 442           (  y3 + y5 + y6  < &4) \/
 443           (  y2 + y3 + y4 + y5 < &4) \/
 444           (  y1 + y3 + y4  + y6 < &4) \/
 445           (  y1 + y2 + y5 + y6  < &4) \/ (fcc_ineq y1 y2 y3 y4 y5 y6 > &0)
 446       )`;
 447   doc = "FCC main inequality, case 1b.   b: 1<=y5
 448     Region: assume p01 largest wlog, then p0X123 >= 4 ==> p01 >= 4/3.
 449     Assume wlog p02 >= p03, p12, p13, ==> p01X23 >=4 ==> p02 >= 1.
 450      We are bounded away from sqrt(0) on this domain.
 451      FCC cell occurs when {y1,..,y6} = {2,2,0,0,2,2};
 452     Local analysis at FCC is done in mathematica below.
 453    This allows us to assume that some constraint (1)..(6) goes out to 4.04.
 454    y1 = p01, y2 = p02, y3 = p03, y4 = p12, y5 = p13, y6 = p23";
 455   tags = [Cfsqp;Strongdodec];
 456 };;
 457
 458
 459 Ineq.add {
 460   idv="FXZXPNS-1c";
 461   ineq = all_forall `ineq
 462    [
 463     (#1.3,y1,&30);
 464     (&1,y2,&30);
 465     (&0,y3,&1);
 466     (&0,y4,&30);
 467     (&0,y5,&1);
 468     (&1,y6,&30)
 469   ]
 470       (
 471           (  y1 + y2 + y3  < #4.04) \/
 472           (  y1 + y4 + y5  < &4) \/
 473           (  y2 + y4 + y6  < &4) \/
 474           (  y3 + y5 + y6  < &4) \/
 475           (  y2 + y3 + y4 + y5 < &4) \/
 476           (  y1 + y3 + y4  + y6 < &4) \/
 477           (  y1 + y2 + y5 + y6  < &4) \/ (fcc_ineq y1 y2 y3 y4 y5 y6 > &0)
 478       )`;
 479   doc = "FCC main inequality, case 1c.   c: 1<=y6.
 480    (Note: if y3,y5,y6<1, then y3+y5+y6<4 and we are done.)
 481     Region: assume p01 largest wlog, then p0X123 >= 4 ==> p01 >= 4/3.
 482     Assume wlog p02 >= p03, p12, p13, ==> p01X23 >=4 ==> p02 >= 1.
 483      We are bounded away from sqrt(0) on this domain.
 484      FCC cell occurs when {y1,..,y6} = {2,2,0,0,2,2};
 485     Local analysis at FCC is done in mathematica below.
 486    This allows us to assume that some constraint (1)..(6) goes out to 4.04.
 487    y1 = p01, y2 = p02, y3 = p03, y4 = p12, y5 = p13, y6 = p23";
 488   tags = [Cfsqp;Cfsqp_branch 7;Strongdodec];
 489 };;
 490
 491
 492
 493
 494 Ineq.add {
 495   idv="FXZXPNS-2";
 496   ineq = all_forall `ineq
 497    [
 498     (#1.3,y1,#4.04);
 499     (&1,y2,#4.04);
 500     (&0,y3,#4.04);
 501     (&0,y4,#4.04);
 502     (&0,y5,#4.04);
 503     (&0,y6,#4.04)
 504   ]
 505       (
 506           (  y1 + y2 + y3  < &4) \/
 507           (  y1 + y2 + y3  > #4.04) \/
 508           (  y1 + y4 + y5  < #4.04) \/
 509           (  y2 + y4 + y6  < &4) \/
 510           (  y3 + y5 + y6  < &4) \/
 511           (  y2 + y3 + y4 + y5 < &4) \/
 512           (  y1 + y3 + y4  + y6 < &4) \/
 513           (  y1 + y2 + y5 + y6  < &4) \/ (fcc_ineq y1 y2 y3 y4 y5 y6 > &0)
 514       )`;
 515   doc = "FCC main inequality, case 2.  second constraint.
 516     Region: assume p01 largest wlog, then p0X123 >= 4 ==> p01 >= 4/3.
 517     Assume wlog p02 >= p03, p12, p13, ==> p01X23 >=4 ==> p02 >= 1.
 518      We are bounded away from sqrt(0) on this domain.
 519      FCC cell occurs when {y1,..,y6} = {2,2,0,0,2,2};
 520     Local analysis at FCC is done in mathematica below.
 521    This allows us to assume that some constraint (1)..(6) goes out to 4.04.
 522    y1 = p01, y2 = p02, y3 = p03, y4 = p12, y5 = p13, y6 = p23";
 523   tags = [Cfsqp;Strongdodec];
 524 };;
 525
 526
 527 Ineq.add {
 528   idv="FXZXPNS-3";
 529   ineq = all_forall `ineq
 530    [
 531     (#1.3,y1,#4.04);
 532     (&1,y2,#4.04);
 533     (&0,y3,#4.04);
 534     (&0,y4,#4.04);
 535     (&0,y5,#4.04);
 536     (&0,y6,#4.04)
 537   ]
 538       (
 539           (  y1 + y2 + y3  < &4) \/
 540           (  y1 + y2 + y3  > #4.04) \/
 541           (  y1 + y4 + y5  < &4) \/
 542           (  y1 + y4 + y5  > #4.04) \/
 543           (  y2 + y4 + y6  < #4.04) \/
 544           (  y3 + y5 + y6  < &4) \/
 545           (  y2 + y3 + y4 + y5 < &4) \/
 546           (  y1 + y3 + y4  + y6 < &4) \/
 547           (  y1 + y2 + y5 + y6  < &4) \/ (fcc_ineq y1 y2 y3 y4 y5 y6 > &0)
 548       )`;
 549   doc = "FCC main inequality, case 3.  third constraint.
 550     Region: assume p01 largest wlog, then p0X123 >= 4 ==> p01 >= 4/3.
 551     Assume wlog p02 >= p03, p12, p13, ==> p01X23 >=4 ==> p02 >= 1.
 552      We are bounded away from sqrt(0) on this domain.
 553      FCC cell occurs when {y1,..,y6} = {2,2,0,0,2,2};
 554     Local analysis at FCC is done in mathematica below.
 555    This allows us to assume that some constraint (1)..(6) goes out to 4.04.
 556    y1 = p01, y2 = p02, y3 = p03, y4 = p12, y5 = p13, y6 = p23";
 557   tags = [Cfsqp;Strongdodec];
 558 };;
 559
 560
 561 Ineq.add {
 562   idv="FXZXPNS-4";
 563   ineq = all_forall `ineq
 564    [
 565     (#1.3,y1,#4.04);
 566     (&1,y2,#4.04);
 567     (&0,y3,#4.04);
 568     (&0,y4,#4.04);
 569     (&0,y5,#4.04);
 570     (&0,y6,#4.04)
 571   ]
 572       (
 573           (  y1 + y2 + y3  < &4) \/
 574           (  y1 + y2 + y3  > #4.04) \/
 575           (  y1 + y4 + y5  < &4) \/
 576           (  y1 + y4 + y5  > #4.04) \/
 577           (  y2 + y4 + y6  < &4) \/
 578           (  y2 + y4 + y6  > #4.04) \/
 579           (  y3 + y5 + y6  < #4.04) \/
 580           (  y2 + y3 + y4 + y5 < &4) \/
 581           (  y1 + y3 + y4  + y6 < &4) \/
 582           (  y1 + y2 + y5 + y6  < &4) \/ (fcc_ineq y1 y2 y3 y4 y5 y6 > &0)
 583       )`;
 584   doc = "FCC main inequality, case 4.  fourth constraint.
 585     Region: assume p01 largest wlog, then p0X123 >= 4 ==> p01 >= 4/3.
 586     Assume wlog p02 >= p03, p12, p13, ==> p01X23 >=4 ==> p02 >= 1.
 587      We are bounded away from sqrt(0) on this domain.
 588      FCC cell occurs when {y1,..,y6} = {2,2,0,0,2,2};
 589     Local analysis at FCC is done in mathematica below.
 590    This allows us to assume that some constraint (1)..(6) goes out to 4.04.
 591    y1 = p01, y2 = p02, y3 = p03, y4 = p12, y5 = p13, y6 = p23";
 592   tags = [Cfsqp;Strongdodec];
 593 };;
 594
 595
 596 Ineq.add {
 597   idv="FXZXPNS-5";
 598   ineq = all_forall `ineq
 599    [
 600     (#1.3,y1,#4.04);
 601     (&1,y2,#4.04);
 602     (&0,y3,#4.04);
 603     (&0,y4,#4.04);
 604     (&0,y5,#4.04);
 605     (&0,y6,#4.04)
 606   ]
 607       (
 608           (  y1 + y2 + y3  < &4) \/
 609           (  y1 + y2 + y3  > #4.04) \/
 610           (  y1 + y4 + y5  < &4) \/
 611           (  y1 + y4 + y5  > #4.04) \/
 612           (  y2 + y4 + y6  < &4) \/
 613           (  y2 + y4 + y6  > #4.04) \/
 614           (  y3 + y5 + y6  < &4) \/
 615           (  y3 + y5 + y6  > #4.04) \/
 616           (  y2 + y3 + y4 + y5 < #4.04) \/
 617           (  y1 + y3 + y4  + y6 < &4) \/
 618           (  y1 + y2 + y5 + y6  < &4) \/ (fcc_ineq y1 y2 y3 y4 y5 y6 > &0)
 619       )`;
 620   doc = "FCC main inequality, case 5.
 621     Region: assume p01 largest wlog, then p0X123 >= 4 ==> p01 >= 4/3.
 622     Assume wlog p02 >= p03, p12, p13, ==> p01X23 >=4 ==> p02 >= 1.
 623      We are bounded away from sqrt(0) on this domain.
 624      FCC cell occurs when {y1,..,y6} = {2,2,0,0,2,2};
 625     Local analysis at FCC is done in mathematica below.
 626    This allows us to assume that some constraint (1)..(6) goes out to 4.04.
 627    y1 = p01, y2 = p02, y3 = p03, y4 = p12, y5 = p13, y6 = p23";
 628   tags = [Cfsqp;Strongdodec];
 629 };;
 630
 631
 632 Ineq.add {
 633   idv="FXZXPNS-6";
 634   ineq = all_forall `ineq
 635    [
 636     (#1.3,y1,#4.04);
 637     (&1,y2,#4.04);
 638     (&0,y3,#4.04);
 639     (&0,y4,#4.04);
 640     (&0,y5,#4.04);
 641     (&0,y6,#4.04)
 642   ]
 643       (
 644           (  y1 + y2 + y3  < &4) \/
 645           (  y1 + y2 + y3  > #4.04) \/
 646           (  y1 + y4 + y5  < &4) \/
 647           (  y1 + y4 + y5  > #4.04) \/
 648           (  y2 + y4 + y6  < &4) \/
 649           (  y2 + y4 + y6  > #4.04) \/
 650           (  y3 + y5 + y6  < &4) \/
 651           (  y3 + y5 + y6  > #4.04) \/
 652           (  y2 + y3 + y4 + y5 < &4) \/
 653           (  y2 + y3 + y4 + y5 > #4.04) \/
 654           (  y1 + y3 + y4  + y6 < #4.04) \/
 655           (  y1 + y2 + y5 + y6  < &4) \/ (fcc_ineq y1 y2 y3 y4 y5 y6 > &0)
 656       )`;
 657   doc = "FCC main inequality, case 6. Final case.
 658     Region: assume p01 largest wlog, then p0X123 >= 4 ==> p01 >= 4/3.
 659     Assume wlog p02 >= p03, p12, p13, ==> p01X23 >=4 ==> p02 >= 1.
 660      We are bounded away from sqrt(0) on this domain.
 661      FCC cell occurs when {y1,..,y6} = {2,2,0,0,2,2};
 662     Local analysis at FCC is done in mathematica below.
 663    This allows us to assume that some constraint (1)..(6) goes out to 4.04.
 664    y1 = p01, y2 = p02, y3 = p03, y4 = p12, y5 = p13, y6 = p23";
 665   tags = [Cfsqp;Strongdodec];
 666 };;
 667
 668
 669
 670 (*
 671
 672 "Mathematica code used to prove local optimality in an explicit epsilon neighborhood
 673   of BCC";
 674
 675 "definition of objective function mu:";
 676
 677 SellingVolume2 = p01*p02*p03 + p01*p03*p12 + p02*p03*p12 +
 678       p01*p02*p13 + p02*p03*p13 + p01*p12*p13 + p02*p12*p13 +
 679       p03*p12*p13 + p01*p02*p23 + p01*p03*p23 + p01*p12*p23 +
 680       p02*p12*p23 + p03*p12*p23 + p01*p13*p23 + p02*p13*p23 + p03*p13*p23;
 681
 682 p0X123 = p01 + p02 + p03;
 683 p1X023 = p01 + p12 + p13;
 684 p2X013 = p02 + p12 + p23;
 685 p3X012 = p03 + p13 + p23;
 686 p01X23 = p02 + p03 + p12 + p13;
 687 p02X13 = p01 + p03 + p12 + p23;
 688 p03X12 = p01 + p02 + p13 + p23;
 689 F01X23 = p01*p23*Sqrt[p01X23];
 690 F02X13 = p02*p13*Sqrt[p02X13];
 691 F03X12 = p03*p12*Sqrt[p03X12];
 692 F0X123 = (p12*p13 + p12*p23 + p13*p23)*Sqrt[p0X123];
 693 F1X023 = (p02*p03 + p02*p23 + p03*p23)*Sqrt[p1X023];
 694 F2X013 = (p01*p03 + p01*p13 + p03*p13)*Sqrt[p2X013];
 695 F3X012 = (p01*p02 + p01*p12 + p02*p12)*Sqrt[p3X012];
 696
 697 SellingSurfaceNum = 2*(
 698     F01X23 + F02X13 + F03X12 + F0X123 + F1X023 + F2X013 + F3X012);
 699
 700 muSelling = ((SellingSurfaceNum/2)^3/(SellingVolume2)^(5/2));
 701
 702 "Interval Arithmetic Evaluation";
 703 eval[r_, ii_] := (r /. {p01 -> 1, p02 -> ii, p03 -> ii, p12 ->
 704     ii, p13 -> ii, p23 -> ii});
 705 ii = Interval[{0.99999, 1}];
 706
 707 "Gradient";
 708 dd1[i_] := D[muSelling1, {pv[[i]]}];
 709 dds = Table[dd1[i], {i, 1, 5}];
 710 {"Gradient", tab1i = Table[eval[dds[[i]], 1], {i, 1, 5}] // Simplify}
 711
 712 "Principal Minors";
 713 pv = {p02, p03, p12, p13, p23};
 714 dd[i_, j_] :=  D[muSelling1, {pv[[i]]}, {pv[[j]]}];
 715 ddij = Table[eval[dd[i, j], ii], {i, 1, 5}, {j, 1, 5}];
 716 minor[k_] := Det[Table[ddij[[i]][[j]], {i, 1, k}, {j, 1, k}]];
 717 {"Table of Principal Minors", Table[minor[k], {k, 1, 5}]  }
 718
 719 *)
 720
 721 (*
 722 (* Local Optimality of fcc lattice with respect to surface area of Voronoi
 723 cell among those lattices containing a unit ball packing. *)
 724 (* Corollary : Gauss on fcc lattice optimality.  *)
 725
 726 surfcc = SellingSurfaceNum/(SellingVolume2)^(1/2);
 727   (* at fcc, we get value 12 Sqrt[2], or about 16.9705627484771405856202646905 *)
 728
 729 epssub = {p01 -> 2 + eps01, p02 -> 2 + eps02, p03 -> 0 +
 730     eps03, p12 -> 0 + eps12, p13 -> 2 + eps13, p23 -> 2 + eps23};
 731 systemFcc = {t1 == p0X123 - 4, t2 == p1X023 - 4, t3 == p2X013 -
 732           4, t4 == p3X012 - 4,
 733         t5 == p01X23 - 4, t6 == p02X13 - 4} /. epssub;
 734 tsub = Solve[systemFcc, {eps01, eps02, eps03, eps12, eps13, eps23}];
 735 surt = (surfcc /. epssub) /. tsub;
 736 tvar = {t1, t2, t3, t4, t5, t6}; (* positive coordinates on constrained domain *)
 737 grad = Table[D[surt, tvar[[i]]], {i, 1, 6}];
 738 evalt[f_, a_] := f /. {t1 -> a, t2 -> a, t3 -> a, t4 -> a, t5 -> a, t6 -> a};
 739 {"gradient is pos", evalt[grad, Interval[{0, 0.04}]]}
 740 *)
 741
 742
 743 end;;