1 (* ========================================================================== *)
2 (* FLYSPECK - BOOK FORMALIZATION *)
4 (* Definitions: (General definitions file) *)
6 (* Author: Thomas C. Hales *)
8 (* ========================================================================== *)
11 (* needs "Multivariate/flyspeck.ml";; *)
14 module Sphere = struct
19 `(!c. ineq [] c <=> c)
20 /\ (!a x b xs c. ineq (CONS (a,x,b) xs) c <=> a <= x /\ x <= b ==> ineq xs c)`;;
23 let mk_forall = mk_binder "!" in
24 itlist (curry mk_forall) (sort (<) (frees bod)) bod;;
26 let periodic = new_definition `periodic (f:num->A) n = (!i. (f (i+n) = f (i:num)))`;;
29 (* symmetric difference *)
31 let SDIFF = new_definition `SDIFF X Y = (X DIFF Y) UNION (Y DIFF X)`;;
34 With the following definition, we should have
35 atn2(x,y) = pi - Arg(complex(-- x, y)),
36 including in the degenerate cases. In this sense, it is the "right" definition.
39 let atn2 = new_definition(`atn2(x,y) =
40 if ( abs y < x ) then atn(y / x) else
41 (if (&0 < y) then ((pi / &2) - atn(x / y)) else
42 (if (y < &0) then (-- (pi/ &2) - atn (x / y)) else ( pi )))`);;
44 (* ------------------------------------------------------------------ *)
47 f = \x. a x^2 + b x + c, extract a b c
50 let abc_of_quadratic =
51 new_definition `abc_of_quadratic f =
55 ((p + n)/(&2) - c, (p -n)/(&2), c)` ;;
58 let quadratic_root_plus =
59 new_definition `quadratic_root_plus (a, b, c) =
60 ( -- b + sqrt(b pow 2 - &4 * a * c))/ (&2 * a)`;;
64 (* ------------------------------------------------------------------ *)
66 let sqrt8 = new_definition (`sqrt8 = sqrt (&8) `);;
67 let sqrt2 = new_definition (`sqrt2 = sqrt (&2) `);;
68 let sqrt3 = new_definition `sqrt3 = sqrt(&3)`;;
71 let pi_rt18 = new_definition(`pi_rt18= pi/(sqrt (&18))`);;
74 (* ------------------------------------------------------------------ *)
75 (* This polynomial is essentially the Cayley-Menger determinant. *)
76 (* ------------------------------------------------------------------ *)
77 let delta_x = new_definition (`delta_x x1 x2 x3 x4 x5 x6 =
78 x1*x4*(--x1 + x2 + x3 -x4 + x5 + x6) +
79 x2*x5*(x1 - x2 + x3 + x4 -x5 + x6) +
80 x3*x6*(x1 + x2 - x3 + x4 + x5 - x6)
81 -x2*x3*x4 - x1*x3*x5 - x1*x2*x6 -x4*x5*x6`);;
83 let delta_y = new_definition `delta_y y1 y2 y3 y4 y5 y6 =
84 delta_x (y1*y1) (y2*y2) (y3*y3) (y4*y4) (y5*y5) (y6*y6)`;;
86 let edge_flat = new_definition`edge_flat y1 y2 y3 y5 y6 =
87 sqrt(quadratic_root_plus (abc_of_quadratic (
88 \x4. -- delta_x (y1*y1) (y2*y2) (y3*y3) x4 (y5*y5) (y6*y6))))`;;
90 let edge_flat2_x = new_definition `edge_flat2_x x1 x2 x3 x4 x5 x6 =
91 (edge_flat (sqrt x1) (sqrt x2) (sqrt x3) (sqrt x5) (sqrt x6)) pow 2`;; (* x4 dummy *)
93 let edge_flat_x = new_definition`edge_flat_x x1 x2 x3 (x4:real) x5 x6 =
94 edge_flat (sqrt x1) (sqrt x2) (sqrt x3) (sqrt x5) (sqrt x6)`;;
97 (* ------------------------------------------------------------------ *)
98 (* The partial derivative of delta_x with respect to x4. *)
99 (* ------------------------------------------------------------------ *)
101 let delta_x4= new_definition(`delta_x4 x1 x2 x3 x4 x5 x6
102 = -- x2* x3 - x1* x4 + x2* x5
103 + x3* x6 - x5* x6 + x1* (-- x1 + x2 + x3 - x4 + x5 + x6)`);;
105 let delta_x6 = new_definition(`delta_x6 x1 x2 x3 x4 x5 x6
106 = -- x1 * x2 - x3*x6 + x1 * x4
107 + x2* x5 - x4* x5 + x3*(-- x3 + x1 + x2 - x6 + x4 + x5)`);;
109 (* ------------------------------------------------------------------ *)
111 (* ------------------------------------------------------------------ *)
113 let ups_x = new_definition(`ups_x x1 x2 x6 =
114 --x1*x1 - x2*x2 - x6*x6
115 + &2 *x1*x6 + &2 *x1*x2 + &2 *x2*x6`);;
118 let eta_x = new_definition(`eta_x x1 x2 x3 =
119 (sqrt ((x1*x2*x3)/(ups_x x1 x2 x3)))
122 let eta_y = new_definition(`eta_y y1 y2 y3 =
128 let rho_x = new_definition(`rho_x x1 x2 x3 x4 x5 x6 =
129 --x1*x1*x4*x4 - x2*x2*x5*x5 - x3*x3*x6*x6 +
130 (&2)*x1*x2*x4*x5 + (&2)*x1*x3*x4*x6 + (&2)*x2*x3*x5*x6`);;
133 let chi_x = new_definition(`chi_x x1 x2 x3 x4 x5 x6
134 = -- (x1*x4*x4) + x1*x4*x5 + x2*x4*x5 - x2*x5*x5
135 + x1*x4*x6 + x3*x4*x6 +
136 x2*x5*x6 + x3*x5*x6 - (&2) * x4*x5*x6 - x3*x6*x6`);;
138 (* ------------------------------------------------------------------ *)
139 (* The formula for the dihedral angle of a simplex. *)
140 (* The variables xi are the squares of the lengths of the edges. *)
141 (* The angle is computed along the first edge (x1). *)
142 (* ------------------------------------------------------------------ *)
144 let dih_x = new_definition(`dih_x x1 x2 x3 x4 x5 x6 =
145 let d_x4 = delta_x4 x1 x2 x3 x4 x5 x6 in
146 let d = delta_x x1 x2 x3 x4 x5 x6 in
147 pi/ (&2) + atn2( (sqrt ((&4) * x1 * d)),-- d_x4)`);;
150 let dih_y = new_definition(`dih_y y1 y2 y3 y4 y5 y6 =
151 let (x1,x2,x3,x4,x5,x6)= (y1*y1,y2*y2,y3*y3,y4*y4,y5*y5,y6*y6) in
152 dih_x x1 x2 x3 x4 x5 x6`);;
154 let dih2_y = new_definition(`dih2_y y1 y2 y3 y4 y5 y6 =
155 dih_y y2 y1 y3 y5 y4 y6`);;
157 let dih3_y = new_definition(`dih3_y y1 y2 y3 y4 y5 y6 =
158 dih_y y3 y1 y2 y6 y4 y5`);;
160 let dih2_x = new_definition(`dih2_x x1 x2 x3 x4 x5 x6 =
161 dih_x x2 x1 x3 x5 x4 x6`);;
163 let dih3_x = new_definition(`dih3_x x1 x2 x3 x4 x5 x6 =
164 dih_x x3 x1 x2 x6 x4 x5`);;
167 (* ------------------------------------------------------------------ *)
168 (* Harriot-Euler formula for the area of a spherical triangle *)
169 (* in terms of the angles: area = alpha+beta+gamma - pi *)
170 (* ------------------------------------------------------------------ *)
172 let sol_x = new_definition(`sol_x x1 x2 x3 x4 x5 x6 =
173 (dih_x x1 x2 x3 x4 x5 x6) +
174 (dih_x x2 x3 x1 x5 x6 x4) + (dih_x x3 x1 x2 x6 x4 x5) - pi`);;
176 let sol_y = new_definition(`sol_y y1 y2 y3 y4 y5 y6 =
177 (dih_y y1 y2 y3 y4 y5 y6) +
178 (dih_y y2 y3 y1 y5 y6 y4) + (dih_y y3 y1 y2 y6 y4 y5) - pi`);;
181 (* ------------------------------------------------------------------ *)
182 (* squander functions -2009 version *)
184 (* ------------------------------------------------------------------ *)
186 let interp = new_definition `interp x1 y1 x2 y2 x = y1 + (x - x1) * (y2- y1)/(x2-x1)`;;
188 let const1 = new_definition `const1 = sol_y (&2) (&2) (&2) (&2) (&2) (&2) / pi`;;
190 let ly = new_definition `ly y = interp (&2) (&1) (#2.52) (&0) y`;;
192 let rho = new_definition `rho y = &1 + const1 - const1* ly y`;;
194 let h0 = new_definition `h0 = #1.26`;;
196 let sol0 = new_definition `sol0 = &3 * acs (&1 / &3) - pi`;;
198 let rho_fun = new_definition `rho_fun y = &1 + (inv (&2 * h0 - &2)) * (inv pi) * sol0 * (y - &2)`;;
200 let rhazim = new_definition `rhazim y1 y2 y3 y4 y5 y6 = rho y1 * dih_y y1 y2 y3 y4 y5 y6`;;
202 let lnazim = new_definition `lnazim y1 y2 y3 y4 y5 y6 = ly y1 * dih_y y1 y2 y3 y4 y5 y6`;;
204 let taum = new_definition
205 `taum y1 y2 y3 y4 y5 y6 = sol_y y1 y2 y3 y4 y5 y6 * (&1 + const1) -
206 const1 * (lnazim y1 y2 y3 y4 y5 y6 + lnazim y2 y3 y1 y5 y6 y4 + lnazim y3 y1 y2 y6 y4 y5)`;;
208 let tauV = new_definition
209 `tauV (v1:real^3) v2 v3 = taum (norm v1) (norm v2) (norm v3) (dist(v2,v3)) (dist(v1,v3)) (dist(v1,v2))`;;
211 let node2_y = new_definition `node2_y f y1 y2 y3 y4 y5 y6 = f y2 y3 y1 y5 y6 y4`;;
213 let node3_y = new_definition `node3_y f y1 y2 y3 y4 y5 y6 = f y3 y1 y2 y6 y4 y5`;;
215 let rhazim2 = new_definition `rhazim2 = node2_y rhazim`;;
217 let rhazim3 = new_definition `rhazim3 = node3_y rhazim`;;
219 let dih4_y = define `dih4_y y1 y2 y3 y4 y5 y6 = dih_y y4 y2 y6 y1 y5 y3`;;
221 let dih5_y = define `dih5_y y1 y2 y3 y4 y5 y6 = dih_y y5 y1 y6 y2 y4 y3`;;
223 let dih6_y = define `dih6_y y1 y2 y3 y4 y5 y6 = dih_y y6 y1 y5 y3 y4 y2`;;
225 let rhazim4 = define `rhazim4 y1 y2 y3 y4 y5 y6 = rho y4 * dih4_y y1 y2 y3 y4 y5 y6`;;
227 let rhazim5 = define `rhazim5 y1 y2 y3 y4 y5 y6 = rho y5 * dih5_y y1 y2 y3 y4 y5 y6 `;;
229 let rhazim6 = define `rhazim6 y1 y2 y3 y4 y5 y6 = rho y6 * dih6_y y1 y2 y3 y4 y5 y6`;;
231 (* two tetrahedra, shared edges y2 y3 y4 *)
233 let tauq = new_definition `tauq y1 y2 y3 y4 y5 y6 y7 y8 y9 = taum y1 y2 y3 y4 y5 y6 + taum y7 y2 y3 y4 y8 y9`;;
235 (* ------------------------------------------------------------------ *)
236 (* The Cayley-Menger formula for the volume of a simplex *)
237 (* The variables xi are the squares of the lengths of the edges. *)
238 (* ------------------------------------------------------------------ *)
240 let vol_x = new_definition(`vol_x x1 x2 x3 x4 x5 x6 =
241 (sqrt (delta_x x1 x2 x3 x4 x5 x6))/ (&12)`);;
243 (* ------------------------------------------------------------------ *)
244 (* Some lower dimensional funcions and Rogers simplices. *)
245 (* ------------------------------------------------------------------ *)
247 let arclength = new_definition(`arclength a b c =
248 pi/(&2) + (atn2( (sqrt (ups_x (a*a) (b*b) (c*c))),(c*c - a*a -b*b)))`);;
251 let volR = new_definition(`volR a b c =
252 (sqrt (a*a*(b*b-a*a)*(c*c-b*b)))/(&6)`);;
254 let solR = new_definition(`solR a b c =
255 (&2)*atn2( sqrt(((c+b)*(b+a))), sqrt ((c-b)*(b-a)))`);;
257 let dihR = new_definition(`dihR a b c =
258 atn2( sqrt(b*b-a*a),sqrt (c*c-b*b))`);;
260 let rad2_x = new_definition(`rad2_x x1 x2 x3 x4 x5 x6 =
261 (rho_x x1 x2 x3 x4 x5 x6)/((delta_x x1 x2 x3 x4 x5 x6)*(&4))`);;
263 (* aff is deprecated *)
264 let aff = new_definition `aff = ( hull ) affine`;;
266 let lin_combo = new_definition `lin_combo V f = vsum V (\v. f v % (v:real^N))`;;
268 let affsign = new_definition `affsign sgn s t (v:real^A) = (?f.
269 (v = lin_combo (s UNION t) f) /\ (!w. t w ==> sgn (f w)) /\ (sum (s UNION t) f = &1))`;;
271 let sgn_gt = new_definition `sgn_gt = (\t. (&0 < t))`;;
273 let sgn_ge = new_definition `sgn_ge = (\t. (&0 <= t))`;;
275 let sgn_lt = new_definition `sgn_lt = (\t. (t < &0))`;;
277 let sgn_le = new_definition `sgn_le = (\t. (t <= &0))`;;
279 let cone = new_definition `cone v S:real^A->bool = affsign sgn_ge {v} S`;;
281 let cone0 = new_definition `cone0 v S:real^A->bool = affsign sgn_gt {v} S`;;
283 let conv0 = new_definition `conv0 S:real^A->bool = affsign sgn_gt {} S`;;
285 let aff_gt_def = new_definition `aff_gt = affsign sgn_gt`;;
287 let aff_ge_def = new_definition `aff_ge = affsign sgn_ge`;;
289 let aff_lt_def = new_definition `aff_lt = affsign sgn_lt`;;
291 let aff_le_def = new_definition `aff_le = affsign sgn_le`;;
293 let voronoi_open = new_definition
294 `voronoi_open S (v:real^N) = { x | !w. ((S w) /\ ~(w=v)) ==> (dist( x, v) < dist( x, w)) }`;;
296 let voronoi_closed = new_definition
297 `voronoi_closed S (v:real^N) = { x | !w. S w ==> (dist( x, v) <= dist( x, w)) }`;;
299 let VORONOI_SET = new_definition
300 `voronoi_set V W = INTERS { voronoi_closed V (v:real^3) | v IN W}`;;
303 new_definition `voronoi_list V wl = voronoi_set (V:real^3 -> bool) (set_of_list wl)`;;
306 new_definition `voronoi_nondg V ul =
307 ((LENGTH ul < 5) /\ (set_of_list ul SUBSET V) /\
308 (aff_dim (voronoi_list V ul) + &(LENGTH ul) = &4))`;;
310 let INITIAL_SUBLIST = new_definition
311 `initial_sublist (xl:(A)list) zl = (?yl. zl = APPEND xl yl)`;;
313 let BARV = new_definition `barV V k ul = ((LENGTH ul = k+1) /\
314 (!vl. (initial_sublist vl ul /\ 0 < LENGTH vl) ==> voronoi_nondg V vl))`;;
316 let TRUNCATE_SIMPLEX = new_definition
317 `truncate_simplex j (ul:(A) list) =
318 @vl. ( LENGTH vl = j+1 /\ initial_sublist vl ul)`;;
320 let OMEGA_LIST_N = define `(omega_list_n V ul 0 = HD ul) /\
321 (omega_list_n V ul (SUC i) =
322 closest_point (voronoi_list V (truncate_simplex (SUC i) ul)) (omega_list_n V ul i))`;;
324 let OMEGA_LIST = new_definition
325 `!V ul. omega_list V ul = omega_list_n V ul ((LENGTH ul) - 1)`;;
327 let ROGERS = new_definition `!V ul. rogers V ul =
328 convex hull (IMAGE (omega_list_n V ul) { j | j < LENGTH ul })`;;
330 (* AFFINE GEOMETRY *)
332 let line = new_definition `line x = (?(v:real^N) w. ~(v =w) /\ (x = affine hull {v,w}))`;;
334 let plane = new_definition
335 `plane x = (?(u:real^N) v w. ~(collinear {u,v,w}) /\ (x = affine hull {u,v,w}))`;;
337 let closed_half_plane = new_definition
338 `closed_half_plane x = (?(u:real^N) v w. ~(collinear {u,v,w}) /\ (x = aff_ge {u,v} {w}))`;;
340 let open_half_plane = new_definition
341 `open_half_plane x = (?(u:real^N) v w. ~(collinear {u,v,w}) /\ (x = aff_gt {u,v} {w}))`;;
343 let closed_half_space = new_definition
344 `closed_half_space x = (?u v w w'. ~(coplanar {u,v,w,w'}) /\ (x = aff_ge {u,v,w} {w'}))`;;
346 let open_half_space = new_definition
347 `open_half_space x = (?u v w w'. ~(coplanar {u,v,w,w'}) /\ (x = aff_gt {u,v,w} {w'}))`;;
349 let bis = new_definition `bis (u:real^N) v = {x | dist(x,u) = dist(x,v)}`;;
351 let bis_le = new_definition `bis_le (u:real^N) v = {x | dist(x,u) <= dist(x,v) }`;;
353 let bis_lt = new_definition `bis_lt (u:real^N) v = {x | dist(x,u) < dist(x,v) }`;;
355 let BIS_SYM = prove(`!p (q:real^N). bis p q = bis q p`,
356 REWRITE_TAC[bis] THEN SET_TAC[]);;
358 let circumcenter = new_definition `circumcenter S = @v. ( (affine hull S) v /\ (?c. !w. (S w) ==> (c = dist(v,w))))`;;
360 let radV = new_definition `radV S = @c. !w. (S w) ==> (c = dist(circumcenter S,w))`;;
362 let orientation = new_definition `orientation S v sgn = affsign sgn (S DIFF {v}) {v} (circumcenter S)`;;
364 let arcV = new_definition `arcV u v w = acs (( (v - u) dot (w - u))/((norm (v-u)) * (norm (w-u))))`;;
366 let dihV = new_definition `dihV w0 w1 w2 w3 =
370 let vap = ( vc dot vc) % va - ( va dot vc) % vc in
371 let vbp = ( vc dot vc) % vb - ( vb dot vc) % vc in
372 arcV (vec 0) vap vbp`;;
374 (* conventional ordering on variables *)
376 let ylist = new_definition `ylist w0 w1 w2 w3 =
377 ((dist (w0, w1)),(dist( w0, w2)),(dist( w0, w3)),(dist( w2, w3)),(dist( w1, w3)),(dist( w1, w2)))`;;
379 let xlist = new_definition `xlist w0 w1 w2 w3 =
380 let (y1,y2,y3,y4,y5,y6) = ylist w0 w1 w2 w3 in
381 (y1 pow 2, y2 pow 2, y3 pow 2, y4 pow 2, y5 pow 2, y6 pow 2)`;;
383 let euler_p = new_definition `euler_p v0 v1 v2 v3 =
384 (let (y1,y2,y3,y4,y5,y6) = ylist v0 v1 v2 v3 in
388 y1*y2*y3 + y1*( w2 dot w3) + y2*( w3 dot w1) + y3*( w1 dot w2))`;;
390 let orthonormal = new_definition `orthonormal e1 e2 e3 =
391 (( e1 dot e1 = &1) /\ (e2 dot e2 = &1) /\ ( e3 dot e3 = &1) /\
392 ( e1 dot e2 = &0) /\ ( e1 dot e3 = &0) /\ ( e2 dot e3 = &0) /\
393 (&0 < (e1 cross e2) dot e3))`;;
395 let cyclic_set = new_definition `cyclic_set W v w =
396 (~(v=w) /\ (FINITE W) /\ (!p q h. W p /\ W q /\ (p = q + h % (v - w)) ==> (p=q)) /\
397 (W INTER (affine hull {v,w}) = EMPTY))`;;
399 (* projection to plane orth to e *)
400 let projection = new_definition `projection e x =
401 x - ((x dot e)/(e dot e)) % e`;;
403 let azim_cycle = new_definition `azim_cycle W v w p =
404 if (W SUBSET {p}) then p else
405 (@u. ~(u=p) /\ W u /\
406 (!q. ~(q=p) /\ W q ==>
407 (azim v w p u < azim v w p q) \/
408 ((azim v w p u = azim v w p q) /\
409 (norm (projection (w-v) (u-v)) <= norm (projection (w-v) (q-v))))))`;;
411 let packing = new_definition `packing (S:real^3 -> bool) = (!u v. S u /\ S v /\ ~(u = v) ==> (&2 <= dist( u, v)))`;;
413 let packing_lt = prove(`packing (V:real^3 -> bool) =
414 (!u:real^3 v:real^3. (u IN V) /\ (v IN V) /\ (dist( u, v) < &2) ==>
416 REWRITE_TAC[packing;IN;REAL_ARITH `x<y <=> ~(y<= x)`]
419 let saturated =new_definition `saturated S= (!x. ?(y:real^N). y IN S /\ dist (x,y)< &2)`;;
421 let c_cone = new_definition `c_cone (v,w:real^3, r:real)={x:real^3 | ((x-v) dot w = norm (x-v)* norm w* r)}`;;
423 let circular_cone =new_definition
424 `circular_cone (V:real^3-> bool)=
425 (? (v,w:real^3)(r:real). ~(w = vec 0) /\ V = c_cone (v,w,r))`;;
427 let null_equiv = new_definition `null_equiv (s,t :real^3->bool)=(? (B:real^3-> bool). NULLSET B /\
428 ((s DIFF t) UNION (t DIFF s)) SUBSET B)`;;
430 let radial = new_definition `radial r (x:real^A) C <=> (C SUBSET ball (x,r)) /\ (!u. (x+u) IN C ==> (!t.(t> &0) /\ (t* norm u < r)==>(x+ t % u) IN C))`;;
432 let eventually_radial = new_definition `eventually_radial x C <=> (?r. (r> &0) /\ radial r x (C INTER ball (x, r)))`;;
434 let rconesgn = new_definition `rconesgn sgn v w h = {x:real^A | sgn ((x-v) dot (w-v)) (dist(x,v)*dist(w,v)*h)}`;;
436 let rcone_ge = new_definition `rcone_ge = rconesgn ( >= )`;;
438 let rcone_gt = new_definition `rcone_gt = rconesgn ( > )`;;
440 let rcone_lt = new_definition `rcone_lt = rconesgn ( < )`;;
442 let rcone_eq = new_definition `rcone_eq = rconesgn ( = )`;;
444 let scale = new_definition `scale (t:real^3) (u:real^3) = vector[t$1 * u$1; t$2 * u$2; t$3 * u$3]`;;
446 let vol_solid_triangle = new_definition `vol_solid_triangle v0 v1 v2 v3 r =
447 let a123 = dihV v0 v1 v2 v3 in
448 let a231 = dihV v0 v2 v3 v1 in
449 let a312 = dihV v0 v3 v1 v2 in
450 (a123 + a231 + a312 - pi)*(r pow 3)/(&3)`;;
452 let vol_frustt_wedge = new_definition `vol_frustt_wedge v0 v1 v2 v3 h a =
453 (azim v0 v1 v2 v3)*(h pow 3)*(&1/(a*a) - &1)/(&6)`;;
455 (* volume of intersection of conic cap and wedge *)
456 let vol_conic_cap_wedge = new_definition `vol_conic_cap_wedge v0 v1 v2 v3 r c =
457 (azim v0 v1 v2 v3)*(&1 - c)*(r pow 3)/(&3)`;;
459 let vol_conv = new_definition `vol_conv v1 v2 v3 v4 =
460 let x12 = dist(v1,v2) pow 2 in
461 let x13 = dist(v1,v3) pow 2 in
462 let x14 = dist(v1,v4) pow 2 in
463 let x23 = dist(v2,v3) pow 2 in
464 let x24 = dist(v2,v4) pow 2 in
465 let x34 = dist(v3,v4) pow 2 in
466 sqrt(delta_x x12 x13 x14 x34 x24 x23)/(&12)`;;
468 let vol_rect = new_definition `vol_rect a b =
469 if (a$1 < b$1) /\ (a$2 < b$2) /\ (a$3 < b$3) then (b$3-a$3)*(b$2-a$2)*(b$1-a$1) else &0`;;
471 let vol_ball_wedge = new_definition `vol_ball_wedge v0 v1 v2 v3 r =
472 (azim v0 v1 v2 v3)*(&2)*(r pow 3)/(&3)`;;
474 let ortho0 = new_definition `ortho0 x v1 v2 v3 = conv0 {x,x+v1,x+v1+v2,x+v1+v2+v3}`;;
476 let make_point = new_definition `make_point v1 v2 v3 w r1 r2 r3 = @v. (aff_ge {v1,v2,v3} {w} (v:real^3)) /\ (r1 = dist(v1,v)) /\ (r2 = dist(v2,v)) /\ (r3 = dist(v3,v))`;;
478 let abc_param = new_definition `abc_param v0 v1 v2 c =
479 let a = (&1/(&2)) * dist(v0,v1) in
480 let b = radV {v0,v1,v2} in
483 let res = new_definition `!f:A->A s:A->bool x:A. res f s x = if x IN s then f x else x`;;
485 let regular_spherical_polygon_area = new_definition
486 `regular_spherical_polygon_area ca k =
487 &2 * pi - &2 * k * asn (ca * sin (pi /k))`;; (* corrected 2010-06-06. *)
489 let tau0 = new_definition `tau0 = &4 * pi - &20 * sol0`;;
491 let mm1 = new_definition `mm1 = sol0 * sqrt(&8)/tau0`;;
493 let mm2 = new_definition `mm2 = (&6 * sol0 - pi) * sqrt(&2) /(&6 * tau0)`;;
495 let hplus = new_definition `hplus = #1.3254`;;
497 let h0cut = new_definition `h0cut y = if (y <= &2 * h0) then &1 else &0`;;
499 let marchal_quartic = new_definition `marchal_quartic h =
500 (sqrt(&2)-h)*(h- hplus )*(&9*(h pow 2) - &17*h + &3)/
501 ((sqrt(&2) - &1)* &5 *(hplus - &1))`;;
503 let lmfun = new_definition`lmfun h = if (h<=h0) then (h0 - h)/(h0 - &1) else &0`;;
505 let lfun = new_definition `lfun h = (h0 - h)/(h0 - &1)`;;
507 let flat_term = new_definition `flat_term y = sol0 * (y - &2 * h0)/(&2 * h0 - &2)`;;
509 let hminus = new_definition `hminus = @x. #1.2 <= x /\ x < #1.3 /\ marchal_quartic x = lmfun x`;;
511 let cstab=new_definition ` cstab= #3.01`;;
518 let y_of_x = new_definition `y_of_x fx y1 y2 y3 y4 y5 y6 =
519 fx (y1*y1) (y2*y2) (y3*y3) (y4*y4) (y5*y5) (y6*y6)`;;
521 let rad2_y = new_definition `rad2_y = y_of_x rad2_x`;;
523 let delta4_y = new_definition `delta4_y = y_of_x delta_x4`;;
525 (* real and fake 4-cell volumes *)
527 let vol_y = new_definition `vol_y = y_of_x vol_x`;;
529 let vol4f = new_definition `vol4f y1 y2 y3 y4 y5 y6 f =
531 (sol_y y1 y2 y3 y4 y5 y6 +
532 sol_y y1 y5 y6 y4 y2 y3 +
533 sol_y y4 y5 y3 y1 y2 y6 +
534 sol_y y4 y2 y6 y1 y5 y3)
536 (f(y1/ &2)* dih_y y1 y2 y3 y4 y5 y6 +
537 f(y2/ &2)* dih_y y2 y3 y1 y5 y6 y4 +
538 f(y3/ &2)* dih_y y3 y1 y2 y6 y4 y5 +
539 f(y4/ &2)* dih_y y4 y3 y5 y1 y6 y2 +
540 f(y5/ &2)* dih_y y5 y1 y6 y2 y4 y3 +
541 f(y6/ &2)* dih_y y6 y1 y5 y3 y4 y2)`;;
543 let gamma4f = new_definition `gamma4f y1 y2 y3 y4 y5 y6 f =
544 vol_y y1 y2 y3 y4 y5 y6 - vol4f y1 y2 y3 y4 y5 y6 f`;;
546 let gamma4fgcy = new_definition `gamma4fgcy y1 y2 y3 y4 y5 y6 f =
547 gamma4f y1 y2 y3 y4 y5 y6 f`;;
549 (* real and fake 3-cell volumes *)
551 let vol3r = new_definition `vol3r y1 y2 y3 r = vol_y r r r y1 y2 y3`;;
553 let vol3f = new_definition `vol3f y1 y2 y3 r f = (&2 * mm1 / pi) *
554 (sol_y y1 y2 r r r y3 +
555 sol_y y2 y3 r r r y1 +
556 sol_y y3 y1 r r r y2)
558 (f(y1/ &2)* dih_y y1 y2 r r r y3 +
559 f(y2/ &2)* dih_y y2 y3 r r r y1 +
560 f(y3/ &2)* dih_y y3 y1 r r r y2)`;;
562 let gamma3f = new_definition `gamma3f y1 y2 y3 r f = vol3r y1 y2 y3 r - vol3f y1 y2 y3 r f`;;
565 (* error XX. vol2r missing factor of (y/2). *)
567 let vol2r = new_definition `vol2r y r = &2 * pi * (r*r - (y / (&2)) pow 2)/(&3)`;;
569 let vol2f = new_definition `vol2f y r f = (&2 * mm1 / pi) * &2 *pi* (&1- y/ (r * &2))
570 - (&8 * mm2/pi) * &2 * pi * f (y/ (&2)) `;;
574 (* For nonlinear inequalities *)
577 let norm2hh = new_definition `norm2hh y1 y2 y3 y4 y5 y6 =
578 (y1- hminus - hplus) pow 2 + (y2 - &2) pow 2 + (y3 - &2) pow 2 + (y4 - &2) pow 2
579 + (y5 - &2) pow 2 + (y6 - &2) pow 2`;;
581 let bump = new_definition `!h. bump h = #0.005*(&1 - ((h- h0) pow 2)/((hplus - h0) pow 2))`;;
583 let critical_edge_y = new_definition `critical_edge_y y = ((&2*hminus <= y) /\ (y <= &2 *hplus))`;;
585 let beta_bumpA_y = new_definition `beta_bumpA_y y1 y2 y3 y4 y5 y6 =
586 (if critical_edge_y y1 then &1 else &0) *
587 (if y2 < &2 * hminus then &1 else &0) *
588 (if y3 < &2 * hminus then &1 else &0) *
589 (if critical_edge_y y4 then &1 else &0) *
590 (if y5 < &2 * hminus then &1 else &0) *
591 (if y6 < &2 * hminus then &1 else &0) *
592 (bump (y1 / &2) - bump (y4 / &2))`;;
594 let beta_bump_force_y = new_definition `beta_bump_force_y y1 y2 y3 y4 y5 y6 =
595 (bump (y1/ &2) - bump (y4 / &2))`;;
597 let wtcount3_y = new_definition `wtcount3_y y1 y2 y3 =
598 (if critical_edge_y y1 then 1 else 0) +
599 (if critical_edge_y y2 then 1 else 0) +
600 (if critical_edge_y y3 then 1 else 0) `;;
602 let wtcount6_y = new_definition
603 `wtcount6_y y1 y2 y3 y4 y5 y6 = wtcount3_y y1 y2 y3 + wtcount3_y y4 y5 y6`;;
606 hack for numerical procedures. It is irrelevant for formal proofs. *)
607 let machine_eps = new_definition `machine_eps = &0`;;
610 (* nonlinear inequalities *)
612 let a_spine5 = new_definition `a_spine5 = #0.0560305`;;
614 let b_spine5 = new_definition `b_spine5 = -- #0.0445813`;;
616 let beta_bump_lb = new_definition `beta_bump_lb = -- #0.005`;;
618 let gamma23f = new_definition `gamma23f y1 y2 y3 y4 y5 y6 w1 w2 r f =
619 (gamma3f y1 y2 y6 r f / &w1 + gamma3f y1 y3 y5 r f / &w2
620 + (dih_y y1 y2 y3 y4 y5 y6 - dih_y y1 y2 r r r y6 - dih_y y1 y3 r r r y5) * (vol2r y1 r - vol2f y1 r f)/(&2 * pi)) `;;
622 let gamma23f_126_03 = new_definition `gamma23f_126_03 y1 y2 y3 y4 y5 y6 w1 r f =
623 (gamma3f y1 y2 y6 r f / &w1
624 + (dih_y y1 y2 y3 y4 y5 y6 - dih_y y1 y2 r r r y6 - #0.03) * (vol2r y1 r - vol2f y1 r f)/(&2 * pi)) `;;
626 let gamma23f_red_03 = new_definition `gamma23f_red_03 y1 y2 y3 y4 y5 y6 r f =
627 (dih_y y1 y2 y3 y4 y5 y6 - &2 * #0.03) * (vol2r y1 r - vol2f y1 r f)/(&2 * pi)`;;
629 let pathL = new_definition `pathL (a,b) = (a,(a+b)/ &2)`;;
631 let pathR = new_definition `pathR (a,b) = ((a+b)/ &2,b)`;;
634 let rotate2 = new_definition `rotate2 f x1 x2 x3 x4 x5 x6 =
635 f x2 x3 x1 x5 x6 x4`;;
637 let rotate3 = new_definition `rotate3 f x1 x2 x3 x4 x5 x6 =
638 f x3 x1 x2 x6 x4 x5`;;
640 let rotate4 = new_definition `rotate4 f x1 x2 x3 x4 x5 x6 =
641 f x4 x2 x6 x1 x5 x3`;;
643 let rotate5 = new_definition `rotate5 f x1 x2 x3 x4 x5 x6 =
644 f x5 x3 x4 x2 x6 x1`;;
646 let rotate6 = new_definition `rotate6 f x1 x2 x3 x4 x5 x6 =
647 f x6 x1 x5 x3 x4 x2`;;
652 let rad2_y = new_definition `rad2_y = y_of_x rad2_x`;;
654 let norm2hh_x = new_definition `norm2hh_x x1 x2 x3 x4 x5 x6 =
655 norm2hh (sqrt x1) (sqrt x2) (sqrt x3) (sqrt x4) (sqrt x5) (sqrt x6)`;;
657 let rhazim_x = new_definition `rhazim_x x1 x2 x3 x4 x5 x6 =
658 rhazim (sqrt x1) (sqrt x2) (sqrt x3) (sqrt x4) (sqrt x5) (sqrt x6)`;;
660 let rhazim2_x = new_definition `rhazim2_x x1 x2 x3 x4 x5 x6 =
661 rhazim2 (sqrt x1) (sqrt x2) (sqrt x3) (sqrt x4) (sqrt x5) (sqrt x6)`;;
663 let rhazim3_x = new_definition `rhazim3_x x1 x2 x3 x4 x5 x6 =
664 rhazim3 (sqrt x1) (sqrt x2) (sqrt x3) (sqrt x4) (sqrt x5) (sqrt x6)`;;
666 let dih4_x = new_definition `dih4_x x1 x2 x3 x4 x5 x6 =
667 dih4_y (sqrt x1) (sqrt x2) (sqrt x3) (sqrt x4) (sqrt x5) (sqrt x6)`;;
669 let dih5_x = new_definition `dih5_x x1 x2 x3 x4 x5 x6 =
670 dih5_y (sqrt x1) (sqrt x2) (sqrt x3) (sqrt x4) (sqrt x5) (sqrt x6)`;;
672 let dih6_x = new_definition `dih6_x x1 x2 x3 x4 x5 x6 =
673 dih6_y (sqrt x1) (sqrt x2) (sqrt x3) (sqrt x4) (sqrt x5) (sqrt x6)`;;
675 let gcy = new_definition `gcy y = &4 * mm1/pi - (&8*mm2/pi) * lmfun(y/ &2)`;;
677 let gchi = new_definition `gchi y = (&4*mm1/pi) - (&504*(mm2/pi))/(&13) + (&200*y*(mm2/pi))/(&13) `;;
680 let gchi1_x = new_definition `gchi1_x x1 x2 x3 x4 x5 x6 = gchi (sqrt x1) * dih_x x1 x2 x3 x4 x5 x6`;;
682 let gchi2_x = new_definition `gchi2_x x1 x2 x3 x4 x5 x6 = gchi (sqrt x2) * dih2_x x1 x2 x3 x4 x5 x6`;;
684 let gchi3_x = new_definition `gchi3_x x1 x2 x3 x4 x5 x6 = gchi (sqrt x3) * dih3_x x1 x2 x3 x4 x5 x6`;;
686 let gchi4_x = new_definition `gchi4_x x1 x2 x3 x4 x5 x6 = gchi (sqrt x4) * dih4_x x1 x2 x3 x4 x5 x6`;;
688 let gchi5_x = new_definition `gchi5_x x1 x2 x3 x4 x5 x6 = gchi (sqrt x5) * dih5_x x1 x2 x3 x4 x5 x6`;;
690 let gchi6_x = new_definition `gchi6_x x1 x2 x3 x4 x5 x6 = gchi (sqrt x6) * dih6_x x1 x2 x3 x4 x5 x6`;;
692 let ldih_x = new_definition `ldih_x x1 x2 x3 x4 x5 x6 =
693 lfun(sqrt(x1) / &2) * dih_x x1 x2 x3 x4 x5 x6`;;
695 let ldih2_x = new_definition `ldih2_x x1 x2 x3 x4 x5 x6 =
696 lfun(sqrt(x2) / &2) * dih2_x x1 x2 x3 x4 x5 x6`;;
698 let ldih3_x = new_definition `ldih3_x x1 x2 x3 x4 x5 x6 =
699 lfun(sqrt(x3) / &2) * dih3_x x1 x2 x3 x4 x5 x6`;;
701 let ldih6_x = new_definition `ldih6_x x1 x2 x3 x4 x5 x6 =
702 lfun(sqrt(x6) / &2) * dih6_x x1 x2 x3 x4 x5 x6`;;
705 (* modified arctangent, analytic continuation of ArcTan[Sqrt[x]]/Sqrt[x] *)
707 let matan = new_definition `matan x =
709 else if (x > &0) then atn (sqrt x) / (sqrt x)
710 else (log ((&1 + sqrt( -- x))/(&1 - sqrt( -- x)))) / (&2 * sqrt (-- x))`;;
712 (* compare "Euler_main_theorem.EULER_ANGLE_SUM_rescal" *)
714 let sol_euler_x = new_definition `sol_euler_x x1 x2 x3 x4 x5 x6 =
715 (let a = sqrt(x1*x2*x3) + sqrt( x1)*(x2 + x3 - x4)/ &2 +
716 sqrt(x2)*(x1 + x3 - x5)/ &2 + sqrt(x3)*(x1 + x2 - x6)/ &2 in
717 &2 * atn2( &2 * a, sqrt (delta_x x1 x2 x3 x4 x5 x6)))`;;
720 let taum_y1 = new_definition
721 `taum_y1 a b y1 (y2:real) (y3:real) (y4:real) (y5:real) (y6:real) =
722 taum (&2) (&2) (&2) a b y1`;;
724 let taum_y2 = new_definition
725 `taum_y2 a b (y1:real) (y2:real) (y3:real) (y4:real) (y5:real) (y6:real) =
726 taum (&2) (&2) (&2) a b y2`;;
728 let taum_y1_y2 = new_definition
729 `taum_y1_y2 a (y1:real) (y2:real) (y3:real) (y4:real) (y5:real) (y6:real) =
730 taum (&2) (&2) (&2) a y1 y2`;;
732 let taum_x1 = new_definition `taum_x1 a b x1 x2 x3 x4 x5 x6 =
733 taum_y1 a b (sqrt x1) (sqrt x2) (sqrt x3) (sqrt x4) (sqrt x5) (sqrt x6)`;;
735 let taum_x2 = new_definition `taum_x2 a b x1 x2 x3 x4 x5 x6 =
736 taum_y2 a b (sqrt x1) (sqrt x2) (sqrt x3) (sqrt x4) (sqrt x5) (sqrt x6)`;;
738 let taum_x1_x2 = new_definition `taum_x1_x2 a x1 x2 x3 x4 x5 x6 =
739 taum_y1_y2 a (sqrt x1) (sqrt x2) (sqrt x3) (sqrt x4) (sqrt x5) (sqrt x6)`;;
741 let arclength_y1 = new_definition
743 (y1:real) (y2:real) (y3:real) (y4:real) (y5:real) (y6:real) =
746 let arclength_y2 = new_definition
748 (y1:real) (y2:real) (y3:real) (y4:real) (y5:real) (y6:real) =
751 let arclength_x1 = new_definition
752 `arclength_x1 a b x1 x2 x3 x4 x5 x6 =
753 arclength_y1 a b (sqrt x1) (sqrt x2) (sqrt x3) (sqrt x4) (sqrt x5) (sqrt x6)`;;
755 let arclength_x2 = new_definition
756 `arclength_x2 a b x1 x2 x3 x4 x5 x6 =
757 arclength_y2 a b (sqrt x1) (sqrt x2) (sqrt x3) (sqrt x4) (sqrt x5) (sqrt x6)`;;
759 let arc_hhn = new_definition `arc_hhn =
760 arclength (&2 * h0) (&2 * h0) (&2)`;;
762 let asn797k = new_definition `asn797k k x2 x3 x4 x5 x6 =
763 k * asn (cos (#0.797) * sin (pi / k))`;;
765 let asnFnhk = new_definition `asnFnhk h k x3 x4 x5 x6 =
766 k * asn (( h * sqrt3 / #4.0 + sqrt(&1 - (h/ &2) pow 2)/ &2) * sin (pi/ k))`;;
768 let lfun_y1 = new_definition `lfun_y1 (y1:real) (y2:real) (y3:real)
769 (y4:real) (y5:real) (y6:real) = lfun y1`;;
771 let acs_sqrt_x1_d4 = new_definition `acs_sqrt_x1_d4 (x1:real) (x2:real) (x3:real) (x4:real) (x5:real) (x6:real) =
772 acs (sqrt(x1)/ &4)`;;
774 let acs_sqrt_x2_d4 = new_definition `acs_sqrt_x2_d4 (x1:real) (x2:real) (x3:real) (x4:real) (x5:real) (x6:real) =
775 acs (sqrt(x2) / &4)`;;
777 let arclength_x_123 = new_definition `arclength_x_123 (x1:real) (x2:real) (x3:real) (x4:real) (x5:real) (x6:real) = arclength (sqrt x1) (sqrt x2) (sqrt x3)`;;
782 let tame_table_d = new_definition
783 `tame_table_d r s = if (r + 2*s > 3)
784 then #0.103 * (&2 - &s) + #0.2759 * (&r + &2* &s - &4)
785 else #0.0`;; (* preferred term over d2_tame *)
788 let eta2_126 = new_definition `eta2_126 x1 (x2:real) (x3:real) (x4:real) (x5:real) x6 =
789 (eta_y (sqrt x1) (sqrt x2) (sqrt x6)) pow 2`;;
791 let eta2_135 = new_definition `eta2_135 x1 (x2:real) (x3:real) (x4:real) (x5:real) (x6:real) =
792 (eta_y (sqrt x1) (sqrt x3) (sqrt x5)) pow 2`;;
794 let eta2_456 = new_definition `eta2_456 (x1:real) (x2:real) (x3:real) (x4:real) (x5:real) (x6:real) =
795 (eta_y (sqrt x4) (sqrt x5) (sqrt x6)) pow 2`;;
797 let num1 = new_definition `num1 e1 e2 e3 a2 b2 c2 =
798 -- &4*((a2 pow 2) *e1 + &8*(b2 - c2)*(e2 - e3) -
799 a2*(&16*e1 + ( b2 - &8 )*e2 + (c2 - &8)*e3))`;;
801 let flat_term_x = new_definition `flat_term_x x = flat_term (sqrt x)`;;
803 let taum_x = new_definition `taum_x x1 x2 x3 x4 x5 x6 =
804 rhazim_x x1 x2 x3 x4 x5 x6 + rhazim2_x x1 x2 x3 x4 x5 x6 +
805 rhazim3_x x1 x2 x3 x4 x5 x6 - (&1 + const1)* pi`;;
807 let eulerA_x = new_definition `eulerA_x x1 x2 x3 x4 x5 x6 =
808 sqrt(x1) * sqrt(x2) * sqrt(x3) + sqrt(x1) * (x2 + x3 - x4) / &2 + sqrt(x2) * (x1 + x3 - x5) / &2 +
809 sqrt(x3) * (x1 + x2 - x6) / &2`;;
811 let delta4_squared_x = new_definition
812 `delta4_squared_x x1 x2 x3 x4 x5 x6 = (delta_x4 x1 x2 x3 x4 x5 x6) pow 2`;;
814 let delta4_squared_y = new_definition `delta4_squared_y = y_of_x delta4_squared_x`;;
816 let x1_delta_x = new_definition `x1_delta_x x1 x2 x3 x4 x5 x6 = x1 * delta_x x1 x2 x3 x4 x5 x6`;;
818 let x1_delta_y = new_definition `x1_delta_y = y_of_x x1_delta_x`;;
820 let delta_126_x = new_definition
821 `delta_126_x (x3s:real) (x4s:real) (x5s:real) (x1:real) (x2:real) (x3:real) (x4:real) (x5:real) (x6:real) =
822 delta_x x1 x2 x3s x4s x5s x6`;;
824 let delta_234_x = new_definition
825 `delta_234_x (x1s:real) (x5s:real) (x6s:real) (x1:real) (x2:real) (x3:real) (x4:real) (x5:real) (x6:real) =
826 delta_x x1s x2 x3 x4 x5s x6s`;;
828 let delta_135_x = new_definition
829 `delta_135_x (x2s:real) (x4s:real) (x6s:real) (x1:real) (x2:real) (x3:real) (x4:real) (x5:real) (x6:real) =
830 delta_x x1 x2s x3 x4s x5 x6s`;;
833 (* Added May 22, 2011 *)
835 let delta_sub1_x = new_definition
836 `delta_sub1_x (x1s:real) (x1:real) (x2:real) (x3:real) (x4:real) (x5:real) (x6:real) =
837 delta_x x1s x2 x3 x4 x5 x6`;;
839 let taum_sub1_x = new_definition `taum_sub1_x x1s (x1:real) x2 x3 x4 x5 x6 =
840 taum_x x1s x2 x3 x4 x5 x6`;;
842 let taum_sub246_x = new_definition
843 `taum_sub246_x x2s x4s x6s (x1:real) (x2:real) x3 (x4:real) x5 (x6:real) =
844 taum_x x1 x2s x3 x4s x5 x6s`;;
846 let taum_sub345_x = new_definition
847 `taum_sub345_x x3s x4s x5s (x1:real) (x2:real) (x3:real) (x4:real) (x5:real) (x6:real) =
848 taum_x x1 x2 x3s x4s x5s x6`;;