legacy/inequalities/definitions_keplerC.ml

   1 (* Start from the beginning of the text of Unabridged Proof\r
   2    of the Kepler Conjecture. Version Nov 26, 2003 *)\r
   3 \r
   4     (*  #use "load_def_kepler.ml";;   *)\r
   5 \r
   6 (*let mk_vec3 = kepler_def `mk_vec3 y1 y2 y3 = \r
   7   (\i. if (i=0) then y1 else if (i=1) then y2 else if (i=2) then y3\r
   8    else (&.0))`;;\r
   9 *)\r
  10 \r
  11 (*\r
  12 let mk_lattice = kepler_def `mk_lattice v1 v2 v3 =\r
  13   { z | ?m1 m2 m3. (z = (dest_int m1) *# v1 + (dest_int m2) *# v2 +\r
  14       (dest_int m3) *# v3 ) }`;;\r
  15 \r
  16 let fcc_packing = kepler_def `fcc_packing = \r
  17   let v1 = mk_vec3 (&.2) (&.0) (&.0) in\r
  18   let v2 = mk_vec3 (&.1) (sqrt (&.3)) (&.0) in\r
  19   let v3 = mk_vec3 (&.1) ((sqrt (&.3))/(&.3)) ((&.2)*sqrt(&.6)/(&.3)) in\r
  20   mk_lattice v1 v2 v3`;;\r
  21 \r
  22 \r
  23 let hcp_packing = kepler_def `hcp_packing = \r
  24   let v1 = mk_vec3 (&.2) (&.0) (&.0) in\r
  25   let v2 = mk_vec3 (&.1) (sqrt (&.3)) (&.0) in\r
  26   let v3 = mk_vec3 (&.1) ((sqrt (&.3))/(&.3)) ((&.2)*sqrt(&.6)/(&.3)) in\r
  27   let v4 = mk_vec3 (&.0) (&.0) ((&.4)*sqrt(&.6)/(&.3)) in\r
  28   let L = mk_lattice v1 v2 v4 in\r
  29   L UNION (set_translate v3 L)`;;\r
  30 *)\r
  31 \r
  32 (* B(x,r).  Changed from closed in text.  *)\r
  33 \r
  34 (* ------------------------------------------------------------------ *)\r
  35 (* The following definitions also appear in Jordan/misc_defs_and_lemmas.ml *)\r
  36 (* ------------------------------------------------------------------ *)\r
  37 \r
  38 (* mk_local_interface "kepler";; *)\r
  39 \r
  40 (* we are switching from real3 (based on num->real)  *)\r
  41 (* to real^3. *)\r
  42 (* The file definitions_keplerC.ml still depends on euclid.\r
  43    These definitions should be left in until the transition\r
  44    is complete. -tch 7/9/2008 *)\r
  45 \r
  46 overload_interface\r
  47  ("+", `euclid_plus:(num->real)->(num->real)->(num->real)`);;\r
  48 \r
  49 make_overloadable "*#" `:A -> B -> B`;;\r
  50 \r
  51 let euclid_scale = kepler_def\r
  52   `euclid_scale t f = \ (i:num). (t* (f i))`;;\r
  53 \r
  54 overload_interface ("*#",`euclid_scale`);; (* `kepler'euclid_scale`);; *)\r
  55 \r
  56 parse_as_infix("*#",(20,"right"));;\r
  57 \r
  58 let euclid_neg = kepler_def `euclid_neg (f:num->real) = \ (i:num). (-- (f i))`;;\r
  59 \r
  60 (* This is highly ambiguous: -- f x can be read as\r
  61    (-- f) x or as -- (f x).  *)\r
  62 overload_interface ("--",`euclid_neg`);; (* `kepler'euclid_neg`);; *)\r
  63 \r
  64 overload_interface\r
  65   ("-", `euclid_minus:(num->real)->(num->real)->(num->real)`);;\r
  66 \r
  67 let euclid_plus = kepler_def\r
  68   `euclid_plus (f:num->real) g = \ (i:num). (f i) + (g i)`;;\r
  69 \r
  70 let euclid_minus = kepler_def\r
  71   `euclid_minus (f:num->real) g = \(i:num). (f i) - (g i)`;;\r
  72 \r
  73 let euclid = kepler_def `euclid n v <=> !m. (n <= m) ==> (v (m:num) = &0)`;;\r
  74 \r
  75 \r
  76 let euclid0 = kepler_def `euclid0 = \(i:num). &0`;;\r
  77 \r
  78 let coord = kepler_def `coord i (f:num->real) = f i`;;\r
  79 \r
  80 let dot = kepler_def `dot f g =\r
  81   let (n = (min_num (\m. (euclid m f) /\ (euclid m g)))) in\r
  82   sum (0,n) (\i. (f i)*(g i))`;;\r
  83 \r
  84 let norm = kepler_def `norm f = sqrt(dot f f)`;;\r
  85 \r
  86 \r
  87 let d_euclid = kepler_def `d_euclid f g = norm (f - g)`;;\r
  88 \r
  89 \r
  90 \r
  91 let real3_exists = prove( `?f. (!n. (n> 2) ==> (f n = &0))`,\r
  92        EXISTS_TAC `(\j:num. &0)` THEN\r
  93        BETA_TAC THEN (REWRITE_TAC[])\r
  94                         );;\r
  95 \r
  96 let real3 = new_type_definition "real3" ("mk_real3","dest_real3") real3_exists;;\r
  97 \r
  98 overload_interface\r
  99  ("+", `real3_plus:real3->real3 ->real3`);;\r
 100 \r
 101 let real3_scale = new_definition\r
 102   `real3_scale t v = mk_real3 (t *# dest_real3 v)`;;\r
 103 \r
 104 overload_interface ("*#",`real3_scale`);;\r
 105 \r
 106 let real3_neg = new_definition `real3_neg v = mk_real3 (-- dest_real3 v)`;;\r
 107 \r
 108 (* This is highly ambiguous: -- f x can be read as\r
 109    (-- f) x or as -- (f x).  *)\r
 110 overload_interface ("--",`real3_neg`);; \r
 111 \r
 112 overload_interface\r
 113   ("-", `real3_minus:real3->real3->real3`);;\r
 114 \r
 115 let real3_plus =new_definition\r
 116   `real3_plus v w = mk_real3 (euclid_plus (dest_real3 v) (dest_real3 w))`;;\r
 117 \r
 118 let real3_minus = new_definition\r
 119   `real3_minus v w = mk_real3 (euclid_minus (dest_real3 v)  (dest_real3 w))`;;\r
 120 \r
 121 \r
 122 (* No need for this one.  v$i does the same thing. *)\r
 123 \r
 124 let coord3 = new_definition `coord3 i v = coord i (dest_real3 v)`;;\r
 125 \r
 126 let open_ball = new_definition\r
 127   `open_ball(X,d) (x:A) r = { y | (X x) /\ (X y) /\ (d x y <. r) }`;;\r
 128 \r
 129 let ball3 = kepler_def `ball3 x r = \r
 130       open_ball (euclid 3,d_euclid) x r`;;\r
 131 \r
 132 (* B(x,r,Lambda) *)\r
 133 let ball3_lambda = kepler_def \r
 134   `ball3_lambda x r Lambda =\r
 135    ( ball3 x r) INTER (UNIONS(IMAGE (\v. ball3 v (&.1)) Lambda))`;;\r
 136 \r
 137 (* delta(x,r,Lambda) - THIS DEFINITION IS NOT WORKING (VU KHAC KY)*)\r
 138 \r
 139 \r
 140 \r
 141 (*let delta_finite = kepler_def\r
 142   `delta_finite x r Lambda = \r
 143     (vol 3 (ball3_lambda x r Lambda))/(vol 3 (ball3 x r))`;;\r
 144 *)\r
 145 \r
 146 \r
 147 \r
 148 \r
 149 (* Lambda(x,r) *)\r
 150 let truncated_packing = kepler_def\r
 151   `truncated_packing x r Lambda = \r
 152    Lambda INTER (ball3_lambda x r Lambda)`;;\r
 153 \r
 154 (* Omega(v,Lambda) *)\r
 155 let closed_voronoi_cell = kepler_def\r
 156   `closed_voronoi_cell v Lambda = \r
 157     {x | euclid 3 x /\ \r
 158        (!w. (Lambda w) ==> (d_euclid x v <= d_euclid x w)) }`;;\r
 159 \r
 160 let open_voronoi_cell = kepler_def\r
 161   `open_voronoi_cell v Lambda = \r
 162     {x | euclid 3 x /\ \r
 163        (!w. (Lambda w) /\ ~(w = v) ==> (d_euclid x v <. d_euclid x w)) }`;;\r
 164 let fcc_voronoi_volume = \r
 165   `!v. (fcc_packing v) ==>\r
 166      (vol 3 (open_voronoi_cell v fcc_packing) = (sqrt(&.32)))`;;\r
 167 \r
 168 let hcp_voronoi_volume = \r
 169   `!v. (hcp_packing v) ==>\r
 170      (vol 3 (open_voronoi_cell v hcp_packing) = sqrt(&.32))`;;\r
 171 \r
 172 let negligible = kepler_def \r
 173   `negligible Lambda A = ?C. (!r x. (euclid 3 x) /\ (&.1 <=. r) ==>\r
 174      ITSET (\v acc. A(v) +. acc) (truncated_packing x r Lambda) (&.0)\r
 175      <=. (C * r * r))`;;\r
 176 \r
 177 \r
 178 (* fcc_compatible : THIS DEFINITION IS NOT WORKING (VUKHACKY) *)\r
 179 \r
 180 (*\r
 181 let fcc_compatible = kepler_def\r
 182   ` fcc_compatible Lambda (A:Lambda->real) = ( !v. (Lambda v) ==> \r
 183     (sqrt(&.32) <=. (( vol 3 (open_voronoi_cell v Lambda)) + (A v))))`;;\r
 184 *)\r
 185 \r
 186 let compatible_density = \r
 187   `!Lambda. \r
 188     (saturated_packing Lambda) /\\r
 189     (Lambda SUBSET (euclid 3)) /\\r
 190     (?A. (fcc_compatible Lambda A) /\ (negligible Lambda A)) ==> \r
 191     (?C. (! x r. (euclid 3 x) /\ (&.1 <=. r) ==>\r
 192         (delta_finite x r Lambda <=. pi/(sqrt (&.18)) + C/r)))`;;\r
 193 \r
 194 let kepler_conjecture = \r
 195   `!Lambda. ?C.  !x r. \r
 196      (Lambda SUBSET (euclid 3)) /\\r
 197      (saturated_packing Lambda) /\ \r
 198      (euclid 3 x) /\ (&.1 <=. r) ==>\r
 199       ( (delta_finite x r Lambda <=. pi/(sqrt (&. 18)) + C/r))`;;\r
 200 \r
 201 (* skipped some top-down material that doesn't make sense at this point *)\r
 202 \r
 203 (*\r
 204 \r
 205 \r
 206 let tetrahedron_vol = kepler_def \r
 207   `tetrahedron_vol = \r
 208   let v0 = mk_vec3 (&.0) (&.0) (&.0) in \r
 209   let v1 = mk_vec3 (&.2) (&.0) (&.0) in\r
 210   let v2 = mk_vec3 (&.1) (sqrt (&.3)) (&.0) in\r
 211   let v3 = mk_vec3 (&.1) ((sqrt (&.3))/(&.3)) ((&.2)*sqrt(&.6)/(&.3)) in\r
 212   vol 3 (convex_hull {v0,v1,v2,v3})`;;\r
 213 \r
 214 \r
 215 let tetrahedron_ball_vol = kepler_def \r
 216   `tetrahedron_ball_vol = \r
 217   let v0 = mk_vec3 (&.0) (&.0) (&.0) in \r
 218   let v1 = mk_vec3 (&.2) (&.0) (&.0) in\r
 219   let v2 = mk_vec3 (&.1) (sqrt (&.3)) (&.0) in\r
 220   let v3 = mk_vec3 (&.1) ((sqrt (&.3))/(&.3)) ((&.2)*sqrt(&.6)/(&.3)) in\r
 221   vol 3 (convex_hull {v0,v1,v2,v3} INTER \r
 222       (UNIONS (IMAGE (\v. ball3 v (&.1)) {v0,v1,v2,v3})))`;;\r
 223 \r
 224 let dtet_lemma = \r
 225   `dtet = (tetrahedron_ball_vol/tetrahedron_vol)`;;\r
 226 \r
 227 let octahedron_vol = kepler_def \r
 228   `octahedron_vol = \r
 229   let v0 = mk_vec3 (&.0) (&.0) (-- (sqrt(&.2))) in \r
 230   let v1 = mk_vec3 (&.1) (&.1) (&.0) in\r
 231   let v2 = mk_vec3 (&.1) (--(&.1)) (&.0) in\r
 232   let v3 = mk_vec3 (-- (&.1)) (&.1) (&.0) in\r
 233   let v4 = mk_vec3 (-- (&.1)) (-- (&.1)) (&.0) in\r
 234   let v5 = mk_vec3 (&.0) (&.0) ( (sqrt(&.2))) in \r
 235   let S = {v0,v1,v2,v3,v4,v5} in\r
 236   vol 3 (convex_hull S)`;;\r
 237 \r
 238 let octahedron_ball_vol = kepler_def \r
 239   `octahedron_ball_vol = \r
 240   let v0 = mk_vec3 (&.0) (&.0) (-- (sqrt(&.2))) in \r
 241   let v1 = mk_vec3 (&.1) (&.1) (&.0) in\r
 242   let v2 = mk_vec3 (&.1) (--(&.1)) (&.0) in\r
 243   let v3 = mk_vec3 (-- (&.1)) (&.1) (&.0) in\r
 244   let v4 = mk_vec3 (-- (&.1)) (-- (&.1)) (&.0) in\r
 245   let v5 = mk_vec3 (&.0) (&.0) ( (sqrt(&.2))) in \r
 246   let S = {v0,v1,v2,v3,v4,v5} in\r
 247   vol 3  (convex_hull S INTER \r
 248       (UNIONS (IMAGE (\v. ball3 v (&.1)) S)))`;;\r
 249 \r
 250 let doct_lemma = \r
 251   `doct = (octahedron_ball_vol/octahedron_vol)`;;\r
 252 \r
 253 let dtet_doct = \r
 254   `pi/(sqrt(&.18)) = dtet/(&.3) + (&.2)*doct/(&.3)`;;\r
 255 \r
 256 let pt_dtet = \r
 257   `pt = (sqrt(&.2))*dtet - pi/(&.3)`;;\r
 258 \r
 259 let pt_doct = \r
 260   `pt = (-- (&.2))*(sqrt(&.2)*doct - pi/(&.3))`;;\r
 261 *)\r
 262 \r
 263 \r
 264 \r
 265 \r
 266 \r
 267 \r
 268 \r
 269 \r
 270 (* Construction of the Q-system *) \r
 271 \r
 272 (* VU KHAC KY *)\r
 273 let packing = kepler_def \r
 274 `packing Lambda = (!v w. (((Lambda v)/\(Lambda w)/\(norm(v-w) < &2))==>(v=w)))`;;\r
 275 \r
 276 let two_t0 = kepler_def `two_t0 = #2.51 `;;\r
 277 (*THE NAME OF QUASI_REGULAR_TRIANGLE HAS BEEN CHANGED INTO QUASI_REGULAR_TRIG (VU KHAC KY)*)\r
 278 let quasi_regular_trig = kepler_def\r
 279   `quasi_regular_trig Lambda S = ((S HAS_SIZE 3) /\\r
 280    (S SUBSET Lambda) /\\r
 281    (!v w. (((S v ) /\ (S w)) ==> (d_euclid w v <= two_t0))))`;;\r
 282 (*THE NAME OF SIMPLEX HAS BEEN CHANGED INTO SIMPLX ( VU KHAC KY)*)\r
 283 let simplx = kepler_def `simplx Lambda S = ((S SUBSET Lambda) /\(S HAS_SIZE 4))`;;\r
 284 \r
 285 let quasi_regular_tet = kepler_def\r
 286   `quasi_regular_tet Lambda S = ((simplx Lambda S) /\ \r
 287     (!v w. ((S v) /\ (S w)) ==> (d_euclid w v <= two_t0)))`;;\r
 288 \r
 289 let two_to_2t0 = kepler_def `two_to_2t0 x =\r
 290         (((&2)<= x) /\ (x <= two_t0))`;;\r
 291 \r
 292 let twot0_to_sqrt8 = kepler_def `twot0_to_sqrt8 x = \r
 293         ((two_t0 <= x) /\ (x <= sqrt8))`;;\r
 294 \r
 295 let two_to_sqrt8 = kepler_def `two_to_sqrt8 x =\r
 296         (((&2)<= x) /\ (x <= sqrt8))`;;\r
 297 \r
 298 let strict_twot0_to_sqrt8 = kepler_def `strict_twot0_to_sqrt8 x = \r
 299         ((two_t0 < x) /\ (x < sqrt8))`;;\r
 300 \r
 301 let pre_quarter = kepler_def \r
 302 `pre_quarter Lambda S = ((simplx Lambda S) /\ (!v w. (((Lambda v)/\(Lambda w))==>(two_to_sqrt8 (d_euclid v w )))))`;;\r
 303 \r
 304 let quarter = kepler_def\r
 305   `quarter Lambda S = ((pre_quarter Lambda S) /\\r
 306     (?v w. (S v) /\ (S w) /\ (twot0_to_sqrt8 (d_euclid v w))/\ \r
 307        (!x y. (((S x) /\ (S y) /\ (two_t0 <= (d_euclid x y) )) ==>({x,y}={v,w}) ))))`;;\r
 308 \r
 309 let strict_quarter = kepler_def\r
 310   `strict_quarter Lambda S = ( (quarter Lambda S) /\\r
 311     (?v w. (S v) /\ (S w) /\ (strict_twot0_to_sqrt8 (d_euclid w v))))`;;\r
 312 \r
 313 let diagonal = kepler_def\r
 314   `diagonal S d = ((d SUBSET S) /\ \r
 315      (?d1 d2. (d = {d1,d2}) /\ \r
 316          (!u v. (S u) /\ (S v) ==>(d_euclid u v <=. d_euclid d1 d2))))`;;\r
 317 \r
 318 \r
 319 \r
 320 (* VU KHAC KY *)\r
 321 let six_point = new_definition  \r
 322   `six_point x1 x2 x3 x4 x5 x6 = \r
 323   (!x y. (((x IN {x1, x2 ,x3, x4, x5, x6})/\(y IN {x1, x2 ,x3, x4, x5, x6}))==>(norm(x-y) >(&0))))`;;\r
 324 \r
 325 let pre_quarter_oct = new_definition\r
 326   `pre_quarter_oct Lambda v w x1 x2 x3 x4 = \r
 327   let S1= {v,w,x1,x2} in\r
 328   let S2= {v,w,x2,x3} in\r
 329   let S3= {v,w,x3,x4} in\r
 330   let S4= {v,w,x4,x1} in\r
 331   (strict_quarter Lambda S1)/\\r
 332   (strict_quarter Lambda S2)/\\r
 333   (strict_quarter Lambda S3)/\\r
 334   (strict_quarter Lambda S4)/\\r
 335   (diagonal S1 {v,w})/\\r
 336   (diagonal S2 {v,w})/\\r
 337   (diagonal S3 {v,w})/\\r
 338   (diagonal S4 {v,w})/\\r
 339   ((convex hull (S1) INTER convex hull (S2))= {})/\\r
 340   ((convex hull (S1) INTER convex hull (S3))= {})/\\r
 341   ((convex hull (S1) INTER convex hull (S4))= {})/\\r
 342   ((convex hull (S2) INTER convex hull (S3))= {})/\\r
 343   ((convex hull (S2) INTER convex hull (S4))= {})/\\r
 344   ((convex hull (S3) INTER convex hull (S4))= {})`;;\r
 345 \r
 346 let quartered_oct = kepler_def \r
 347 `quartered_oct Lambda v w x1 x2 x3 x4 = \r
 348 ((six_point v w x1 x2 x3 x4 )/\(pre_quarter_oct Lambda v w x1 x2 x3 x4))`;;\r
 349 \r
 350 let adjacent_pair = kepler_def \r
 351   `adjacent_pair Lambda v v1 v3 v2 v4 =\r
 352    let Q  = {v, v1, v3, v2} in \r
 353    let Q1 = {v, v1,v3,v4}  in\r
 354   (strict_quarter Lambda Q)/\\r
 355   (strict_quarter Lambda Q1)/\\r
 356   (diagonal Q {v1,v3})/\\r
 357   (diagonal Q1 {v1,v3})/\\r
 358   (conv0 Q INTER conv Q1 = EMPTY )`;;\r
 359 \r
 360 let conflict_diagonal = new_definition \r
 361 `conflict_diagonal Lambda v v1 v3 v2 v4 = ((adjacent_pair Lambda v v1 v3 v2 v4)/\(adjacent_pair Lambda v v2 v4 v1 v3))`;;\r
 362 \r
 363 let inter_position = kepler_def \r
 364 ` inter_position Lambda v v1 v3 v2 v4 =\r
 365   ((conflict_diagonal Lambda v v1 v3 v2 v4) /\\r
 366 (conv {v1,v3} INTER conv0 {v,v2,v4}= EMPTY))`;;\r
 367 \r
 368 let isolated_pair = new_definition \r
 369 `isolated_pair Lambda  Q = \r
 370 ((quarter Lambda Q)/\\r
 371 (~(?v v1 v2 v3 v4. (Q v1)/\(Q v2)/\(Q v3)/\(Q v4)==>(adjacent_pair Lambda v v1 v3 v2 v4))))`;;\r
 372 \r
 373 \r
 374 let anchor = new_definition \r
 375 `anchor Lambda v v1 v2 = ( (Lambda v)/\\r
 376                           (twot0_to_sqrt8 (d_euclid v1 v2))/\\r
 377                           (~(v = v1))/\\r
 378                           (~(v = v2))/\\r
 379                           ( d_euclid v v1 <= two_t0)/\\r
 380                           ( d_euclid v v2 <= two_t0))`;; \r
 381 \r
 382 (* Definition of Q- System  ..to be continue.. *)\r
 383 (*\r
 384 let Q-system = new_definition \r
 385 `Q-system (Lambda:packing) Q = \r
 386  (!S. (Q S)==>\r
 387            (quasi_regular_tet Lambda S)\/\r
 388            (strict_quarter Lambda S) \r
 389 \r
 390 `;;\r
 391 *)\r
 392 \r
 393 \r
 394 \r
 395 \r
 396 \r
 397 \r
 398 \r
 399 \r
 400 \r
 401 \r
 402 \r
 403 \r
 404 \r
 405 \r
 406 let is_qrtet_y = kepler_def\r
 407  `is_qrtet_y y1 y2 y3 y4 y5 y6 =\r
 408         ((two_to_2t0 y1) /\\r
 409         (two_to_2t0 y2) /\\r
 410         (two_to_2t0 y3) /\\r
 411         (two_to_2t0 y4) /\\r
 412         (two_to_2t0 y5) /\\r
 413         (two_to_2t0 y6))`;;\r
 414 \r
 415 let s_to_8 = kepler_def `s_to_8 x =\r
 416         ( (#6.3001 <= x) /\ (x <= (&.8)) )`;;\r
 417 \r
 418 let four_to_s = kepler_def `four_to_s x =\r
 419         (((&.4)<= x) /\ (x <= #6.3001))`;;\r
 420 \r
 421 let is_qrtet_x = kepler_def(`is_qrtet_x x1 x2 x3 x4 x5 x6 =\r
 422         ((four_to_s x1) /\\r
 423         (four_to_s x2) /\\r
 424         (four_to_s x3) /\\r
 425         (four_to_s x4) /\\r
 426         (four_to_s x5) /\\r
 427         (four_to_s x6))`);;\r
 428 \r
 429 let is_upright_quarter_y = kepler_def\r
 430 (`is_upright_quarter_y y1 y2 y3 y4 y5 y6 = \r
 431         ((twot0_to_sqrt8 y1) /\\r
 432         (two_to_2t0 y2) /\\r
 433         (two_to_2t0 y3) /\\r
 434         (two_to_2t0 y4) /\\r
 435         (two_to_2t0 y5) /\\r
 436         (two_to_2t0 y6))`);;\r
 437 \r
 438 let is_upright_quarter_v = kepler_def\r
 439 (`is_upright_quarter_v v0 v1 v2 v3 =\r
 440         is_upright_quarter_y\r
 441         (d_euclid v0 v1) (d_euclid v0 v2) (d_euclid v0 v3)\r
 442         (d_euclid v2 v3) (d_euclid v1 v3) (d_euclid v2 v3)`);;\r
 443 \r
 444 let dih_v = kepler_def(`dih_v v0 v1 v2 v3 =\r
 445         dih_y (d_euclid v0 v1) (d_euclid v0 v2) (d_euclid v0 v3)\r
 446         (d_euclid v2 v3) (d_euclid v1 v3) (d_euclid v1 v2)`);;\r
 447 \r
 448 (* an equivalent definition of dih, except it works better in the\r
 449 degenerate case of delta = 0 in distinguishing between angle 0 and pi *)\r
 450 \r
 451 let dih_alt_x = kepler_def(`dih_alt_x x1 x2 x3 x4 x5 x6 =\r
 452         acs ((delta_x4 x1 x2 x3 x4 x5 x6)/\r
 453                 sqrt((ups_x x1 x2 x6)*(ups_x x1 x3 x5)))`);;\r
 454 \r
 455 let dih_alt_y = kepler_def(`dih_alt_y y1 y2 y3 y4 y5 y6 = \r
 456         let (x1,x2,x3,x4,x5,x6)=(y1*y1,y2*y2,y3*y3,y4*y4,y5*y5,y6*y6) in\r
 457         dih_alt_x x1 x2 x3 x4 x5 x6`);;\r
 458 \r
 459 let dih_alt_v = kepler_def(`dih_alt_v v0 v1 v2 v3 =\r
 460         dih_alt_y (d_euclid v0 v1) (d_euclid v0 v2) (d_euclid v0 v3)\r
 461         (d_euclid v2 v3) (d_euclid v1 v3) (d_euclid v1 v2)`);;\r
 462 \r
 463 (* oriented dihedral takes a value from 0 to < 2pi *)\r
 464 let or_dih_v = kepler_def(`or_dih_v v0 v1 v2 v3 =\r
 465         let w1 =  v1 - v0 in\r
 466         let w2 =  v2 - v0 in\r
 467         let w3 =  v3 - v0 in\r
 468         let triple = triple_product w1 w2 w3 in\r
 469         if (triple > (&0)) then (dih_v v0 v1 v2 v3)\r
 470         else if (triple < (&0)) then ((&2)*pi - (dih_v v0 v1 v2 v3))\r
 471         else (dih_alt_v v0 v1 v2 v3)`);;\r
 472 \r
 473 (* traverse the boundary v1 v2 v3 v4, with the region to the left\r
 474         as you circle around *)\r
 475 \r
 476 let solid_area_quad_v = kepler_def(`solid_area_quad_v v0 v1 v2 v3 v4 =\r
 477         (or_dih_v v0 v1 v2 v4) +\r
 478         (or_dih_v v0 v2 v3 v1) +\r
 479         (or_dih_v v0 v3 v4 v2) +\r
 480         (or_dih_v v0 v4 v1 v3) - (&2)*pi`);;\r
 481 \r
 482 let is_quad_cluster_v = kepler_def(`is_quad_cluster_v v0 v1 v2 v3 v4 =\r
 483         let y1 = d_euclid v0  v1 in\r
 484         let y2 = d_euclid v0  v2 in\r
 485         let y3 = d_euclid v0  v3 in\r
 486         let y4 = d_euclid v0  v4 in\r
 487         let e12 = d_euclid v1 v2 in\r
 488         let e23 = d_euclid v2 v3 in\r
 489         let e34 = d_euclid v3 v4 in\r
 490         let e41 = d_euclid v4 v1 in\r
 491         let d13 = d_euclid v1 v3 in\r
 492         let d24 = d_euclid v2 v4 in\r
 493         (two_to_2t0 y1) /\\r
 494         (two_to_2t0 y2) /\\r
 495         (two_to_2t0 y3) /\\r
 496         (two_to_2t0 y4) /\\r
 497         (two_to_2t0 e12) /\\r
 498         (two_to_2t0 e23) /\\r
 499         (two_to_2t0 e34) /\\r
 500         (two_to_2t0 e41) /\\r
 501         (two_t0 <= d13) /\\r
 502         (two_t0 <= d24) `);;\r
 503 \r
 504 let upright_oct_v = kepler_def(`upright_oct_v v0 w v1 v2 v3 v4 =\r
 505         ((is_upright_quarter_v v0 w v1 v2) /\\r
 506         (is_upright_quarter_v v0 w v2 v3) /\\r
 507         (is_upright_quarter_v v0 w v3 v4) /\\r
 508         (is_upright_quarter_v v0 w v4 v1))`);;\r
 509 \r
 510 let is_pairflat13 = kepler_def(`is_pairflat13 v0 v1 v2 v3 v4 =\r
 511         ((is_quad_cluster_v v0 v1  v2 v3 v4) /\\r
 512         ((d_euclid v1 v3) <= sqrt8))`);;\r
 513 \r
 514 let is_pairflat24 = kepler_def(`is_pairflat24 v0 v1 v2 v3 v4 =\r
 515         (is_quad_cluster_v v0 v1  v2 v3 v4) /\\r
 516         ((d_euclid v2 v4) <= sqrt8)`);;\r
 517 \r
 518 (* we make w lie in the open cone at v0 spanned by (v[i]-v0) *)\r
 519 let is_enclosed = kepler_def\r
 520   `is_enclosed w v0 v1 v2 v3 =\r
 521      (?a0 a1 a2 a3. ((w = a0 *# v0 + a1 *# v1 + a2 *# v2 + a3 *# v3)\r
 522         /\\r
 523                 ((&1) = a0 + a1 + a2 + a3) /\\r
 524                 (a1 > (&0)) /\\r
 525                 (a2 > (&0)) /\\r
 526                 (a3 > (&0))))`;;\r
 527 \r
 528 (* there are some geometry theorems that should be proved here to make\r
 529         sure that everything is as expected.  The edge constraints on\r
 530         a quad cluster constrain the polygon under v1 v2 v3 v4 on the\r
 531         unit sphere to be a simple polygon (fact).\r
 532         The edge constraints that the diagonals are at least sqrt8\r
 533         together with the constraint that the quad region has area\r
 534         at most 2Pi, constrains the region to be convex (fact).  Thus, we\r
 535         can get by without proving the Jordan curve theorem for now.\r
 536         An enclosed point will be one the lies over one of the two\r
 537         simplices formed by drawing either diagonal. *)\r
 538 \r
 539         (* a quad cluster with no flat quarters, and diagonals at least\r
 540                 sqrt8 *)\r
 541 let is_sqrt_quadable = kepler_def(`is_sqrt2_quadable v0 v1 v2 v3 v4 = \r
 542         (is_quad_cluster_v v0 v1 v2 v3 v4) /\\r
 543         (~(is_pairflat13 v0 v1 v2 v3 v4)) /\\r
 544         (~(is_pairflat24 v0 v1 v2 v3 v4))`);;\r
 545 \r
 546 \r
 547         (* define this the same way. The only difference will be in\r
 548         the scoring approximation.  These are the ones that have an\r
 549         upright diagonal of ht <= sqrt8, and for which the formulation\r
 550         bounds of vor0 and -1.04 pt apply *)\r
 551 let is_mixed_quadable = kepler_def(`is_mixed_quadable v0 v1 v2 v3 v4 =\r
 552         (is_quad_cluster_v v0 v1 v2 v3 v4) /\\r
 553         (~(is_pairflat13 v0 v1 v2 v3 v4))/\\r
 554         (~(is_pairflat24 v0 v1 v2 v3 v4))`);;\r
 555 \r
 556         (* define an approximation to sigma.  It will be\r
 557                 actual score if two flat quarters\r
 558                 highest score (w) if four upright quarters in an oct\r
 559                 vor(.,sqrt2) if sqrt_quadable\r
 560                 max(-1.04pt,vor0) if mixed_quadable. *)\r
 561 \r
 562 let eta_v = kepler_def(`\r
 563         eta_v v (i:num) j k = \r
 564               let v1 = (v i) in\r
 565               let v2 = (v j) in\r
 566               let v3 = (v k) in\r
 567               let y1 = d_euclid v2 v3 in\r
 568               let y2 = d_euclid v3 v1 in\r
 569               let y3 = d_euclid v1 v2 in\r
 570               eta_y y1 y2 y3`);;\r
 571 \r
 572 let edge_of_v = kepler_def(`edge_of_v v0 v1 v2 v3 =\r
 573         (d_euclid v0  v1,d_euclid v0  v2,d_euclid v0  v3,\r
 574         d_euclid v2 v3,d_euclid v1 v3,d_euclid v1 v2)`);;\r
 575 \r
 576 let mu_flat_v = kepler_def(`mu_flat_v v0 v1 v2 v3 =\r
 577         let (x1,x2,x3,x4,x5,x6) = edge_of_v v0 v1 v2 v3 in\r
 578         mu_flat_x x1 x2 x3 x4 x5 x6`);;\r
 579 \r
 580 let mu_upright_v = kepler_def(`mu_upright_v v0 v1 v2 v3 =\r
 581         let (x1,x2,x3,x4,x5,x6) = edge_of_v v0 v1 v2 v3 in\r
 582         mu_upright_x x1 x2 x3 x4 x5 x6`);;\r
 583 \r
 584 let sigma_upright_c21_x = kepler_def(`sigma_upright_c21_x \r
 585         x1 x2 x3 x4 x5 x6=\r
 586         mu_upright_x x1 x2 x3 x4 x5 x6`);;\r
 587 \r
 588 let sigma_upright_c40_x = kepler_def(`sigma_upright_c40_x \r
 589         x1 x2 x3 x4 x5 x6=\r
 590         ((&1)/(&2))*\r
 591         ((mu_upright_x x1 x2 x3 x4 x5 x6) +\r
 592          (mu_upright_x x1 x5 x6 x4 x2 x3))`);;\r
 593 \r
 594 let qy_v = kepler_def(`qy_v v0 v1 v2 =\r
 595         let y1 = d_euclid v0  v1 in\r
 596         let y2 = d_euclid v0  v2 in\r
 597         let y3 = d_euclid v1 v2 in\r
 598         qy y1 y2 y3`);;\r
 599 \r
 600 let full_phit = kepler_def(`full_phit h t =\r
 601         ((&1)- (h/t))*((phi h t)-(phi t t))`);;\r
 602 \r
 603 let vort_quad_v = kepler_def(`vort_quad_v v0 v1 v2 v3 v4 t=\r
 604         let sol = solid_area_quad_v v0 v1 v2 v3 v4 in\r
 605         let phit= phi t t in\r
 606         let qy12 = (qy_v v0 v1 v2 t) + (qy_v v0 v2 v1 t) in\r
 607         let qy23 = (qy_v v0 v2 v3 t) + (qy_v v0 v3 v2 t) in\r
 608         let qy34 = (qy_v v0 v3 v4 t) + (qy_v v0 v4 v3 t) in\r
 609         let qy41 = (qy_v v0 v4 v1 t) + (qy_v v0 v1 v4 t) in\r
 610         let d1 = or_dih_v v0 v1 v2 v4 in\r
 611         let d2 = or_dih_v v0 v2 v3 v1 in\r
 612         let d3 = or_dih_v v0 v3 v4 v2 in\r
 613         let d4 = or_dih_v v0 v4 v1 v3 in\r
 614         sol*phit - (&4)*doct*(qy12+qy23+qy34+qy41) +\r
 615                 (d1*(full_phit ((d_euclid v0  v1)/(&2)) t)) +\r
 616                 (d2*(full_phit ((d_euclid v0  v2)/(&2)) t)) +\r
 617                 (d3*(full_phit ((d_euclid v0  v3)/(&2)) t)) +\r
 618                 (d4*(full_phit ((d_euclid v0  v4)/(&2)) t))`);;\r
 619 \r
 620 \r
 621 \r
 622 (* score of an octahedron with an upright diagonal w *)\r
 623 let sigma_upoct_approx_w = kepler_def(`sigma_upoct_approx_w v0 w v1 v2 v3 v4=\r
 624         let t1 = \r
 625            (let (x1,x2,x3,x4,x5,x6) = edge_of_v v0 w v1 v2 in\r
 626            sigma_upright_c40_x x1 x2 x3 x4 x5 x6) in\r
 627         let t2 = \r
 628            (let (x1',x2',x3',x4',x5',x6') = edge_of_v v0 w v2 v3 in\r
 629            sigma_upright_c40_x x1' x2' x3' x4' x5' x6') in\r
 630         let t3 = \r
 631            (let (x1'',x2'',x3'',x4'',x5'',x6'') = edge_of_v v0 w v3 v4 in\r
 632            sigma_upright_c40_x x1'' x2'' x3'' x4'' x5'' x6'') in\r
 633         let t4 = \r
 634            (let (x1''',x2''',x3''',x4''',x5''',x6''') = edge_of_v v0 w v4 v1 in\r
 635            sigma_upright_c40_x x1''' x2''' x3''' x4''' x5''' x6''') in\r
 636         t1+t2+t3+t4`);;\r
 637 \r
 638 (* this is an upper bound on the score of a quad cluster *)\r
 639 let sigma_quad_approx1 = kepler_def(`sigma_quad_approx1 v0 v1 v2 v3 v4=\r
 640         let nonoctpart = (\r
 641         if ((is_pairflat13 v0 v1 v2 v3 v4)/\\r
 642                 (is_pairflat24 v0 v1 v2 v3 v4)) then\r
 643                 (max_real ((mu_flat_v v0 v2 v1 v3)+(mu_flat_v v0 v4 v1 v3))\r
 644                         ((mu_flat_v v0 v1 v2 v4)+(mu_flat_v v0 v3 v2 v4)))\r
 645         else if (is_pairflat13 v0 v1 v2 v3 v4)\r
 646                 then ((mu_flat_v v0 v2 v1 v3)+(mu_flat_v v0 v4 v1 v3))\r
 647         else if (is_pairflat24 v0 v1 v2 v3 v4)\r
 648                 then (((mu_flat_v v0 v1 v2 v4)+(mu_flat_v v0 v3 v2 v4)))\r
 649         else \r
 650                 (max_real\r
 651                         (min_real (--(&104)*pt/(&100)) \r
 652                                 (vort_quad_v v0 v1 v2 v3 v4 t0)) \r
 653                         (vort_quad_v v0 v1 v2 v3 v4 sqrt2))) in\r
 654         let octpart = (if (?w. (upright_oct_v v0 w v1 v2 v3 v4)) then \r
 655                 (sup (\x. ?w. ((upright_oct_v v0 w v1 v2 v3 v4) /\\r
 656                         (x = sigma_upoct_approx_w v0 w v1 v2 v3 v4))))\r
 657         else nonoctpart) in\r
 658                 (max_real octpart nonoctpart)`);;\r
 659 \r
 660 let sigma_quad_approx1_lambda = kepler_def(`sigma_quad_approx1_lambda\r
 661         v0 v1 v2 v3 v4 lambda =\r
 662                 (sigma_quad_approx1 v0 v1 v2 v3 v4) -\r
 663                 (solid_area_quad_v v0 v1 v2 v3 v4)*lambda*zeta*pt`);;\r
 664 \r
 665 \r
 666 \r
 667 (* VU KHAC KY *)\r
 668 \r
 669 \r
 670 \r
 671 \r
 672 \r
 673   \r