1 (* This file now uses Harrison's R^n from hol-light version 2.20.
2 It should also be compatible with the Feb 2008 release of HOL-Light.
3 It is not compatible with pre-2.20 versions of HOL-LIGHT.
9 ["/Users/thomashales/Desktop/flyspeck_google/source/inequalities/"]
14 needs "Multivariate/vectors.ml";; (* Eventually should load entire *)
15 needs "Examples/analysis.ml";; (* multivariate-complex theory. *)
16 needs "Examples/transc.ml";; (* Then it won't need these three. *)
23 let kepler_def = (* local_definition "kepler";; *) new_definition;;
27 let max_real = new_definition(`max_real x y =
28 if (y < x) then x else y`);;
30 let min_real = new_definition(`min_real x y =
31 if (x < y) then x else y`);;
33 let min_num = new_definition
34 `min_num (X:num->bool) = @m. (m IN X) /\ (!n. (n IN X) ==> (m <= n))`;;
37 let deriv = new_definition(`deriv f x = @d. (f diffl d)(x)`);;
38 let deriv2 = new_definition(`deriv2 f = (deriv (deriv f))`);;
40 (* ------------------------------------------------------------------ *)
41 (* Extend atn to allow zero denominators. *)
42 (* ------------------------------------------------------------------ *)
44 (* new argument order 2/14/2008 to make it compatible with the
45 ANSI C arctan2 function. Also reworked for better numerical
46 stability in the regions that matter for us. The way things
47 are defined, it gives atn2(0,0) = pi. This is a bit strange,
48 but we never need its value at the origin anyway.
50 4/19/2008: changed final case to pi along the negative real axis.
53 let atn2 = new_definition(`atn2(x,y) =
54 if ( abs y < x ) then atn(y / x) else
55 (if (&0 < y) then ((pi / &2) - atn(x / y)) else
56 (if (y < &0) then (-- (pi/ &2) - atn (x / y)) else ( pi )))`);;
58 (* ------------------------------------------------------------------ *)
60 let sqrt8 = kepler_def (`sqrt8 = sqrt (&8) `);;
61 let sqrt2 = kepler_def (`sqrt2 = sqrt (&2) `);;
63 (* ------------------------------------------------------------------ *)
65 let t0 = kepler_def (`t0 = (#1.255)`);;
66 let two_t0 = kepler_def(`two_t0 = (#2.51)`);;
67 let square_2t0 = kepler_def(`square_2t0 = two_t0*two_t0`);;
68 let square_4t0 = kepler_def(`square_4t0 = (&4)*square_2t0`);;
69 let pt = kepler_def(`pt = (&4)*(atn (sqrt2/(&5))) - (pi/(&3))`);;
70 let square = kepler_def(`square x = x*x`);;
72 (* ------------------------------------------------------------------ *)
73 (* Standard constants. *)
74 (* ------------------------------------------------------------------ *)
76 let zeta = kepler_def(`zeta= (&1)/((&2) * atn (sqrt2/(&5)))`);;
77 let doct = kepler_def(`doct= (pi - (&2)/zeta)/sqrt8`);;
78 let dtet = kepler_def(`dtet = sqrt2/zeta`);;
79 let pi_rt18 = kepler_def(`pi_rt18= pi/(sqrt (&18))`);;
81 (* ------------------------------------------------------------------ *)
82 (* Technical constants. *)
83 (* ------------------------------------------------------------------ *)
85 let rogers_density=kepler_def(`rogers_density= sqrt2/zeta`);;
86 let compression_cut=kepler_def(`compression_cut=(#1.41)`);;
87 let squander_target=kepler_def(`squander_target =
88 ((&4)*pi*zeta - (&8))*pt`);;
89 let xiV=kepler_def(`xiV=(#0.003521)`);;
90 let xi_gamma=kepler_def(`xi_gamma=(#0.01561)`);;
91 let xi'_gamma=kepler_def(`xi'_gamma=(#0.00935)`);;
92 let xi_kappa=kepler_def(`xi_kappa= -- (#0.029)`);;
93 let xi_kappa_gamma=kepler_def(`xi_kappa_gamma=
95 let pi_max =kepler_def(`xi_max = (#0.006688)`);;
96 let t4 = kepler_def(`t4= (#0.1317)`);;
97 let t5 = kepler_def(`t5= (#0.27113)`);;
98 let t6 = kepler_def(`t6= (#0.41056)`);;
99 let t7 = kepler_def(`t7= (#0.54999)`);;
100 let t8 = kepler_def(`t8= (#0.6045)`);;
101 let t9 = kepler_def(`t9= (#0.6978)`);;
102 let t10= kepler_def(`t10=(#0.7891)`);;
103 let s5 = kepler_def(`s5= --(#0.05704)`);;
104 let s6 = kepler_def(`s6= --(#0.11408)`);;
105 let s7 = kepler_def(`s7= --(#0.17112)`);;
106 let s8 = kepler_def(`s8= --(#0.22816)`);;
107 let s9 = kepler_def(`s9= --(#0.1972)`);;
109 (* Note this is what s10 is in DCG p128, but for the blueprint
110 it should be made -eps0 so that the 8pt bound holds by margin eps0 *)
111 let s10= kepler_def(`s10= #0.0`);;
112 let eps0 = kepler_def(`eps0 = #0.000000000001`);; (* eps0 = 10^-12 *)
115 let Z31 = kepler_def(`Z31 = (#0.00005)`);;
116 let Z32 = kepler_def(`Z32 = -- (#0.05714)`);;
117 let Z33 = kepler_def(`Z33 = s6 - (#3.0)*Z31`);;
118 let Z41 = kepler_def(`Z41 = s5 - Z31`);;
119 let Z42 = kepler_def(`Z42 = s6 - (#2.0)*Z31`);;
121 let D31 = kepler_def(`D31 = t4 - (#0.06585)`);; (* = 0.06585 *)
122 let D32 = kepler_def(`D32 = (#0.13943)`);;
123 let D33 = kepler_def(`D33 = t6 - (#0.06585)*(#3.0)`);;
124 let D41 = kepler_def(`D41 = t5 - (#0.06585)`);;
125 let D42 = kepler_def(`D42 = t6 - (#2.0)*(#0.06585)`);;
126 let D51 = kepler_def(`D51 = t6 - (#0.06585)`);;
128 (* ------------------------------------------------------------------ *)
129 (* Ferguson's thesis constants from DCG-2006-Sec 17.4 *)
130 (* ------------------------------------------------------------------ *)
132 let pp_a1 = kepler_def(`pp_a1 = #0.38606588081240521`);;
133 let pp_a2 = kepler_def(`pp_a2 = #0.4198577862`);;
134 let pp_d0 = kepler_def(`pp_d0 = #1.4674`);;
135 let pp_m = kepler_def(`pp_m = #0.3621`);;
136 let pp_b = kepler_def(`pp_b = #0.49246`);;
137 let pp_a = kepler_def(`pp_a = #0.0739626`);;
138 let pp_bc = kepler_def(`pp_bc = #0.253095`);;
139 let pp_c = kepler_def(`pp_c = #0.1533667634670977`);;
140 let pp_d = kepler_def(`pp_d = #0.2265`);;
141 (* solt0 = Solid[2,2,2,2,2,Sqrt[8]] *)
142 let pp_solt0 = kepler_def(`pp_solt0 = &2 * atn2 (&1, sqrt8)`);;
144 (* ------------------------------------------------------------------ *)
145 (* This polynomial is essentially the Cayley-Menger determinant. *)
146 (* ------------------------------------------------------------------ *)
147 let delta_x = kepler_def (`delta_x x1 x2 x3 x4 x5 x6 =
148 x1*x4*(--x1 + x2 + x3 -x4 + x5 + x6) +
149 x2*x5*(x1 - x2 + x3 + x4 -x5 + x6) +
150 x3*x6*(x1 + x2 - x3 + x4 + x5 - x6)
151 -x2*x3*x4 - x1*x3*x5 - x1*x2*x6 -x4*x5*x6`);;
153 (* ------------------------------------------------------------------ *)
154 (* The partial derivative of delta_x with respect to x4. *)
155 (* ------------------------------------------------------------------ *)
157 let delta_x4 = kepler_def(`delta_x4 x1 x2 x3 x4 x5 x6
158 = -- x2* x3 - x1* x4 + x2* x5
159 + x3* x6 - x5* x6 + x1* (-- x1 + x2 + x3 - x4 + x5 + x6)`);;
161 let delta_x6 = kepler_def(`delta_x6 x1 x2 x3 x4 x5 x6
162 = -- x1 * x2 - x3*x6 + x1 * x4
163 + x2* x5 - x4* x5 + x3*(-- x3 + x1 + x2 - x6 + x4 + x5)`);;
165 (* ------------------------------------------------------------------ *)
167 (* ------------------------------------------------------------------ *)
170 let u_x = kepler_def(
171 `u_x x1 x2 x3 = (--(x1*x1+x2*x2+x3*x3)) +
172 (&2) * (x1*x2+x2*x3+x3*x1)`);;
175 let ups_x = kepler_def(`ups_x x1 x2 x6 =
176 --x1*x1 - x2*x2 - x6*x6
177 + &2 *x1*x6 + &2 *x1*x2 + &2 *x2*x6`);;
180 let eta_x = kepler_def(`eta_x x1 x2 x3 =
181 (sqrt ((x1*x2*x3)/(ups_x x1 x2 x3)))
184 let eta_y = kepler_def(`eta_y y1 y2 y3 =
190 let rho_x = kepler_def(`rho_x x1 x2 x3 x4 x5 x6 =
191 --x1*x1*x4*x4 - x2*x2*x5*x5 - x3*x3*x6*x6 +
192 (&2)*x1*x2*x4*x5 + (&2)*x1*x3*x4*x6 + (&2)*x2*x3*x5*x6`);;
194 let rad2_y = kepler_def(`rad2_y y1 y2 y3 y4 y5 y6 =
195 let (x1,x2,x3,x4,x5,x6)= (y1*y1,y2*y2,y3*y3,y4*y4,y5*y5,y6*y6) in
196 (rho_x x1 x2 x3 x4 x5 x6)/((delta_x x1 x2 x3 x4 x5 x6)*(&4))`);;
199 let chi_x = kepler_def(`chi_x x1 x2 x3 x4 x5 x6
200 = -- (x1*x4*x4) + x1*x4*x5 + x2*x4*x5 - x2*x5*x5
201 + x1*x4*x6 + x3*x4*x6 +
202 x2*x5*x6 + x3*x5*x6 - (&2) * x4*x5*x6 - x3*x6*x6`);;
206 (* ------------------------------------------------------------------ *)
207 (* The formula for the dihedral angle of a simplex. *)
208 (* The variables xi are the squares of the lengths of the edges. *)
209 (* The angle is computed along the first edge (x1). *)
210 (* ------------------------------------------------------------------ *)
212 let dih_x = kepler_def(`dih_x x1 x2 x3 x4 x5 x6 =
213 let d_x4 = delta_x4 x1 x2 x3 x4 x5 x6 in
214 let d = delta_x x1 x2 x3 x4 x5 x6 in
215 pi/ (&2) + atn2( (sqrt ((&4) * x1 * d)),-- d_x4)`);;
218 let dih_y = kepler_def(`dih_y y1 y2 y3 y4 y5 y6 =
219 let (x1,x2,x3,x4,x5,x6)= (y1*y1,y2*y2,y3*y3,y4*y4,y5*y5,y6*y6) in
220 dih_x x1 x2 x3 x4 x5 x6`);;
222 let dih2_y = kepler_def(`dih2_y y1 y2 y3 y4 y5 y6 =
223 dih_y y2 y1 y3 y5 y4 y6`);;
225 let dih3_y = kepler_def(`dih3_y y1 y2 y3 y4 y5 y6 =
226 dih_y y3 y1 y2 y6 y4 y5`);;
228 let dih2_x = kepler_def(`dih2_x x1 x2 x3 x4 x5 x6 =
229 dih_x x2 x1 x3 x5 x4 x6`);;
231 let dih3_x = kepler_def(`dih3_x x1 x2 x3 x4 x5 x6 =
232 dih_x x3 x1 x2 x6 x4 x5`);;
235 (* ------------------------------------------------------------------ *)
236 (* Harriot-Euler formula for the area of a spherical triangle *)
237 (* in terms of the angles: area = alpha+beta+gamma - pi *)
238 (* ------------------------------------------------------------------ *)
240 let sol_x = kepler_def(`sol_x x1 x2 x3 x4 x5 x6 =
241 (dih_x x1 x2 x3 x4 x5 x6) +
242 (dih_x x2 x3 x1 x5 x6 x4) + (dih_x x3 x1 x2 x6 x4 x5) - pi`);;
244 let sol_y = kepler_def(`sol_y y1 y2 y3 y4 y5 y6 =
245 (dih_y y1 y2 y3 y4 y5 y6) +
246 (dih_y y2 y3 y1 y5 y6 y4) + (dih_y y3 y1 y2 y6 y4 y5) - pi`);;
249 (* ------------------------------------------------------------------ *)
250 (* The Cayley-Menger formula for the volume of a simplex *)
251 (* The variables xi are the squares of the lengths of the edges. *)
252 (* ------------------------------------------------------------------ *)
254 let vol_x = kepler_def(`vol_x x1 x2 x3 x4 x5 x6 =
255 (sqrt (delta_x x1 x2 x3 x4 x5 x6))/ (&12)`);;
257 (* ------------------------------------------------------------------ *)
258 (* Some lower dimensional funcions and Rogers simplices. *)
259 (* ------------------------------------------------------------------ *)
261 let beta = kepler_def(`beta psi theta =
262 let arg = ((cos psi)*(cos psi) - (cos theta)*(cos theta))/
263 ((&1) - (cos theta)*(cos theta)) in
266 let arclength = kepler_def(`arclength a b c =
267 pi/(&2) + (atn2( (sqrt (ups_x (a*a) (b*b) (c*c))),(c*c - a*a -b*b)))`);;
270 let volR = kepler_def(`volR a b c =
271 (sqrt (a*a*(b*b-a*a)*(c*c-b*b)))/(&6)`);;
273 let solR = kepler_def(`solR a b c =
274 (&2)*atn2( sqrt(((c+b)*(b+a))), sqrt ((c-b)*(b-a)))`);;
276 let dihR = kepler_def(`dihR a b c =
277 atn2( sqrt(b*b-a*a),sqrt (c*c-b*b))`);;
279 let vorR = kepler_def(`vorR a b c =
280 (&4)*(--doct*(volR a b c) + (solR a b c)/(&3))`);;
282 let denR = kepler_def(`denR a b c =
283 (solR a b c)/((&3)*(volR a b c))`);;
285 let tauR = kepler_def(`tauR a b c =
286 --(volR a b c) + (solR a b c)*zeta*pt`);;
288 let quoin = kepler_def(`quoin a b c =
289 let u = sqrt ((c*c-b*b)/(b*b-a*a)) in
290 if ((a>=b) \/ (b>=c)) then (&0) else
291 (--(a*a + a*c-(&2)*c*c)*(c-a)*atn(u)/(&6) +
293 - ((&2)/(&3))*c*c*c*(atn2((b+c),(u*(b-a)))))`);;
295 let qy = kepler_def(`qy y1 y2 y3 t =
296 quoin (y1/(&2)) (eta_y y1 y2 y3) t`);;
298 let quo_x = kepler_def(`quo_x x y z = qy (sqrt x) (sqrt y) (sqrt z) t0`);;
300 let qn = kepler_def(`qn y1 y2 z t =
301 --(&4)*doct*((qy y1 y2 z t) +(qy y2 y1 z t))`);;
303 let phi = kepler_def(`phi h t =
304 (&2)*((&2) - doct*h*t*(h+t))/(&3)`);;
306 let phi0 = kepler_def(`phi0 =
309 let eta0 = kepler_def(`eta0 h =
310 eta_y ((&2)*h) (two_t0) (&2)`);;
312 let crown = kepler_def(`crown h =
314 (&2)*pi*((&1)- h/e)*(phi h e - phi0)`);;
316 let anc = kepler_def(`anc y1 y2 y6 =
319 let b = eta_y y1 y2 y6 in
321 if (b>c) then (&0) else
322 --(dihR h1 b c)*(crown h1)/((&2)*pi)
323 -(solR h1 b c)*phi0 + (vorR h1 b c)
324 -(dihR h2 b c)*((&1)-h2/t0)*((phi h2 t0)-phi0)
325 -(solR h2 b c)*(phi0) + (vorR h2 b c)`);;
327 let K0 = kepler_def(`K0 y1 y2 y6 =
328 (vorR (y1/(&2)) (eta_y y1 y2 y6) (sqrt2)) +
329 (vorR (y2/(&2)) (eta_y y1 y2 y6) (sqrt2)) -
330 (dihR (y1/(&2)) (eta_y y1 y2 y6) (sqrt2))*
331 (&1 - (y1/(sqrt8)))*(phi (y1/(&2)) sqrt2)`);;
333 let AH = kepler_def(`AH h t = (&1 - (h/t))*((phi h t) - (phi t t))`);;
335 let BHY = kepler_def(`BHY y = (AH (y/(&2)) t0) + phi0`);;
339 (* This definition still needs to be finished *)
340 let overlap_f = kepler_def(
343 let lam = (#1.945) in
344 let dih1 = dih_y y1 t0 y2 lam ell lam in
345 let dih2 = dih_y y2 t0 y1 lam ell lam in
346 let s = sol_y y2 t0 y1 lam ell lam in
347 let phi1 = phi (y1/(&2)) t0 in
348 let phi2 = phi (y2/(&2)) t0 in
349 (&2)*(zeta*pt - phi0)*s
350 + (&2)*dih1*((&1) - (y1/two_t0))*(phi0-phi1)
351 + (&2)*dih2*((&1) - (y2/two_t0))*(phi0-phi2)
352 + xxxxx-- need to insert tau terms ---xxxxx`);;
356 (* ------------------------------------------------------------------ *)
357 (* Analytic and truncated Voronoi function *)
358 (* ------------------------------------------------------------------ *)
360 let KY = kepler_def(`KY y1 y2 y3 y4 y5 y6 =
361 (K0 y1 y2 y6) + (K0 y1 y3 y5) +
362 (dih_y y1 y2 y3 y4 y5 y6)*
363 (&1 - (y1/(sqrt8)))*(phi (y1/(&2)) sqrt2)`);;
365 let KX = kepler_def(`KX x1 x2 x3 x4 x5 x6 =
366 KY (sqrt x1) (sqrt x2) (sqrt x3) (sqrt x4) (sqrt x5) (sqrt x6)`);;
368 let vor_analytic_x = kepler_def(`vor_analytic_x x1 x2 x3 x4 x5 x6 =
369 let del = sqrt (delta_x x1 x2 x3 x4 x5 x6) in
370 let u126 = ups_x x1 x2 x6 in
371 let u135 = ups_x x1 x3 x5 in
372 let u234 = ups_x x2 x3 x4 in
373 let vol = ((&1)/((&48)*del))*
374 ((x1*(x2+x6-x1)+x2*(x1+x6-x2))*(chi_x x4 x5 x3 x1 x2 x6)/u126
375 +(x2*(x3+x4-x2)+x3*(--x3+x4+x2))*(chi_x x6 x5 x1 x3 x2 x4)/u234
376 +(x1*(--x1+x3+x5)+x3*(x1-x3+x5))*(chi_x x4 x6 x2 x1 x3 x5)/u135)
378 (&4)*(--doct*vol + (sol_x x1 x2 x3 x4 x5 x6)/(&3))`);;
380 let vor_analytic_x_flipped = kepler_def(`vor_analytic_x_flipped x1 x2 x3 x4 x5 x6 =
381 vor_analytic_x x1 x5 x6 x4 x2 x3`);;
383 let octavor_analytic_x = kepler_def(`octavor_analytic_x
385 (#0.5)*((vor_analytic_x x1 x2 x3 x4 x5 x6) + (vor_analytic_x_flipped x1 x2 x3 x4 x5 x6))`);;
387 let tau_analytic_x = kepler_def(`tau_analytic_x x1 x2 x3 x4 x5 x6 =
388 (sol_x x1 x2 x3 x4 x5 x6)*zeta*pt -
389 (vor_analytic_x x1 x2 x3 x4 x5 x6)`);;
391 (* bug found 3/21/2008: had parenthesis misplaced. -tch *)
393 let kappa = kepler_def(`kappa y1 y2 y3 y4 y5 y6 =
394 (crown (y1/(&2)))*(dih_y y1 y2 y3 y4 y5 y6)/((&2)*pi)
395 + (anc y1 y2 y6) + (anc y1 y3 y5)`);;
397 let kappa_dih_y = kepler_def(`kappa_dih_y y1 y2 y3 y5 y6 d =
398 (crown (y1/(&2)))*d/((&2)*pi)
399 + (anc y1 y2 y6) + (anc y1 y3 y5)`);;
402 let level_at = kepler_def(`level_at h t = if (h<t) then h else t`);;
404 let vorstar = kepler_def(`vorstar h1 h2 h3 a1 a2 a3 t=
405 let phit = phi t t in
406 a1*((&1)-(level_at h1 t)/t)*(phi h1 t - phit)
407 +a2*((&1)-(level_at h2 t)/t)*(phi h2 t - phit)
408 +a3*((&1)-(level_at h3 t)/t)*(phi h3 t - phit)
409 +(a1+a2+a3-pi)*phit`);;
411 let vort_y = kepler_def(`vort_y y1 y2 y3 y4 y5 y6 t =
415 let a1 = dih_y y1 y2 y3 y4 y5 y6 in
416 let a2 = dih2_y y1 y2 y3 y4 y5 y6 in
417 let a3 = dih3_y y1 y2 y3 y4 y5 y6 in
418 (vorstar h1 h2 h3 a1 a2 a3 t)+
419 (qn y1 y2 y6 t)+(qn y2 y3 y4 t)+(qn y1 y3 y5 t)`);;
421 let vor_0_y = kepler_def(`vor_0_y y1 y2 y3 y4 y5 y6 =
422 vort_y y1 y2 y3 y4 y5 y6 t0`);;
424 let tau_0_y = kepler_def(`tau_0_y y1 y2 y3 y4 y5 y6 =
425 (sol_y y1 y2 y3 y4 y5 y6)*zeta*pt - (vor_0_y y1 y2 y3 y4 y5
428 let vor_0_x = kepler_def(`vor_0_x x1 x2 x3 x4 x5 x6=
429 let (y1,y2,y3,y4,y5,y6) = (sqrt x1,sqrt x2,sqrt x3,
430 sqrt x4,sqrt x5,sqrt x6) in
431 vor_0_y y1 y2 y3 y4 y5 y6`);;
433 let tau_0_x = kepler_def(`tau_0_x x1 x2 x3 x4 x5 x6 =
434 (sol_x x1 x2 x3 x4 x5 x6)*zeta*pt -
435 (vor_0_x x1 x2 x3 x4 x5 x6)`);;
437 let vort_x = kepler_def(`vort_x x1 x2 x3 x4 x5 x6 t =
438 let (y1,y2,y3,y4,y5,y6) = (sqrt x1,sqrt x2,sqrt x3,
439 sqrt x4,sqrt x5,sqrt x6) in
440 vort_y y1 y2 y3 y4 y5 y6 t`);;
442 let tauVt_x = kepler_def(`tauVt_x x1 x2 x3 x4 x5 x6 t =
443 (sol_x x1 x2 x3 x4 x5 x6)*zeta*pt -
444 (vort_x x1 x2 x3 x4 x5 x6 t)`);;
446 let vorA_x = kepler_def(`vorA_x x1 x2 x3 x4 x5 x6 =
447 if ((x5 <= (square (#2.77))) /\ (x6 <= (square (#2.77))) /\
448 (((eta_x x4 x5 x6) < sqrt2)))
449 then (vor_analytic_x x1 x2 x3 x4 x5 x6 )
450 else (vor_0_x x1 x2 x3 x4 x5 x6)`);;
452 let tauA_x = kepler_def(`tauA_x x1 x2 x3 x4 x5 x6 =
453 (sol_x x1 x2 x3 x4 x5 x6)*zeta*pt -
454 (vorA_x x1 x2 x3 x4 x5 x6)`);;
456 let vorC0_x = kepler_def(`vorC0_x x1 x2 x3 x4 x5 x6 =
457 if ((x4 <= (square (#2.77))) \/
458 ((eta_x x4 x5 x6 <= sqrt2) /\ (eta_x x2 x3 x4 <= sqrt2)))
459 then (vor_analytic_x x1 x2 x3 x4 x5 x6 )
460 else (vor_0_x x1 x2 x3 x4 x5 x6)`);;
462 let tauC0_x = kepler_def(`tauC0_x x1 x2 x3 x4 x5 x6 =
463 (sol_x x1 x2 x3 x4 x5 x6)*zeta*pt -
464 (vorC0_x x1 x2 x3 x4 x5 x6)`);;
466 let vorC_x = kepler_def(`vorC_x x1 x2 x3 x4 x5 x6 =
467 if ((x4 <= (square (#2.77))) \/
468 ((eta_x x4 x5 x6 <= sqrt2) /\ (eta_x x2 x3 x4 <= sqrt2)) \/
469 ((square (#2.45) <= x2) /\ ((square (#2.45) <= x6))))
470 then (vor_analytic_x x1 x2 x3 x4 x5 x6 )
471 else (vor_0_x x1 x2 x3 x4 x5 x6)`);;
473 let tauC_x = kepler_def(`tauC_x x1 x2 x3 x4 x5 x6 =
474 (sol_x x1 x2 x3 x4 x5 x6)*zeta*pt -
475 (vorC_x x1 x2 x3 x4 x5 x6)`);;
478 let v0x = kepler_def(`v0x x1 x2 x3 x4 x5 x6 =
479 let (y1,y2,y3) = (sqrt x1,sqrt x2,sqrt x3) in
480 (-- (BHY y1))*y1*(delta_x6 x1 x2 x3 x4 x5 x6) +
481 (BHY y2)* y2* (ups_x x1 x3 x5) +
482 (-- (BHY y3))*y3*(delta_x4 x1 x2 x3 x4 x5 x6)`);;
484 let v1x = kepler_def(`v1x x1 x2 x3 x4 x5 x6 =
485 let (y1,y2,y3) = (sqrt x1,sqrt x2,sqrt x3) in
486 (-- (BHY y1 - (zeta*pt)))*y1*(delta_x6 x1 x2 x3 x4 x5 x6) +
487 (BHY y2 - (zeta*pt))* y2* (ups_x x1 x3 x5) +
488 (-- (BHY y3 - (zeta*pt)))*y3*(delta_x4 x1 x2 x3 x4 x5 x6)`);;
491 (* ------------------------------------------------------------------ *)
492 (* The function gamma is called the "compression" in the proof *)
493 (* of the Kepler conjecture. It is interpreted as a linearized *)
494 (* density of a simplex. *)
495 (* ------------------------------------------------------------------ *)
497 let gamma_x = kepler_def(`gamma_x x1 x2 x3 x4 x5 x6 =
498 --doct*(vol_x x1 x2 x3 x4 x5 x6) + ((&2)/(&3))*
499 ((dih_x x1 x2 x3 x4 x5 x6) + (dih_x x2 x3 x1 x5 x6 x4)+
500 (dih_x x3 x1 x2 x6 x4 x5) + (dih_x x4 x5 x3 x1 x2 x6)+
501 (dih_x x5 x6 x1 x2 x3 x4) + (dih_x x6 x5 x1 x3 x2 x4))
504 let tau_gamma_x = kepler_def(`tau_gamma_x x1 x2 x3 x4 x5 x6 =
505 (sol_x x1 x2 x3 x4 x5 x6)*zeta*pt -
506 (gamma_x x1 x2 x3 x4 x5 x6)`);;
508 (* 1.41^2 = 1.9881 *)
509 let rad2_x = kepler_def(`rad2_x x1 x2 x3 x4 x5 x6 =
510 (rho_x x1 x2 x3 x4 x5 x6)/((delta_x x1 x2 x3 x4 x5 x6)*(&4))`);;
512 let sigma_qrtet_x = kepler_def(`sigma_qrtet_x x1 x2 x3 x4 x5 x6=
513 let r = rad2_x x1 x2 x3 x4 x5 x6 in
514 let r_cut = (#1.9881) in
516 vor_analytic_x x1 x2 x3 x4 x5 x6
517 else gamma_x x1 x2 x3 x4 x5 x6`);;
519 let sigma1_qrtet_x = kepler_def(`sigma1_qrtet_x x1 x2 x3 x4 x5 x6=
520 sigma_qrtet_x x1 x2 x3 x4 x5 x6 - (sol_x x1 x2 x3 x4 x5 x6)*zeta*pt`);;
522 let tau_sigma_x = kepler_def `tau_sigma_x x1 x2 x3 x4 x5 x6=
523 -- (sigma1_qrtet_x x1 x2 x3 x4 x5 x6)`;;
525 let sigma32_qrtet_x = kepler_def(`sigma32_qrtet_x x1 x2 x3 x4 x5 x6=
526 sigma_qrtet_x x1 x2 x3 x4 x5 x6 -
527 (#3.2)*(sol_x x1 x2 x3 x4 x5 x6)*zeta*pt`);;
529 let mu_flat_x = kepler_def(`mu_flat_x x1 x2 x3 x4 x5 x6 =
530 let r1 = eta_x x2 x3 x4 in
531 let r2 = eta_x x4 x5 x6 in
532 if ((sqrt2 <= r1)\/(sqrt2 <= r2))
533 then vor_analytic_x x1 x2 x3 x4 x5 x6
534 else gamma_x x1 x2 x3 x4 x5 x6`);;
536 let taumu_flat_x = kepler_def(`taumu_flat_x x1 x2 x3 x4 x5 x6 =
537 (sol_x x1 x2 x3 x4 x5 x6)*zeta*pt -
538 (mu_flat_x x1 x2 x3 x4 x5 x6)`);;
540 let mu_upright_x = kepler_def(`mu_upright_x x1 x2 x3 x4 x5 x6 =
541 let r1 = eta_x x1 x2 x6 in
542 let r2 = eta_x x1 x3 x5 in
543 if ((sqrt2 <= r1)\/(sqrt2 <= r2))
544 then vor_analytic_x x1 x2 x3 x4 x5 x6
545 else gamma_x x1 x2 x3 x4 x5 x6`);;
547 let mu_flipped_x = kepler_def(`mu_flipped_x x1 x2 x3 x4 x5 x6 =
548 mu_upright_x x1 x5 x6 x4 x2 x3`);;
552 let vor_0_x_flipped = kepler_def(`vor_0_x_flipped x1 x2 x3 x4 x5 x6=
553 vor_0_x x1 x5 x6 x4 x2 x3`);;
555 let octavor0_x = kepler_def(`octavor0_x x1 x2 x3 x4 x5 x6 =
556 (#0.5)* (vor_0_x x1 x2 x3 x4 x5 x6 + (vor_0_x_flipped x1 x2 x3 x4 x5 x6))`);;
558 (* STM changed to use mu_flipped_x and vor_0_x_flipped instead of the definition *)
559 (* let nu_x = kepler_def(`nu_x x1 x2 x3 x4 x5 x6 = *)
562 (* (mu_upright_x x1 x2 x3 x4 x5 x6)+ *)
563 (* (mu_upright_x x1 x5 x6 x4 x2 x3)+ *)
564 (* (vor_0_x x1 x2 x3 x4 x5 x6)- *)
565 (* (vor_0_x x1 x5 x6 x4 x2 x3))`);; *)
566 let nu_x = kepler_def(`nu_x x1 x2 x3 x4 x5 x6 =
569 (mu_upright_x x1 x2 x3 x4 x5 x6)+
570 (mu_flipped_x x1 x2 x3 x4 x5 x6)+
571 (vor_0_x x1 x2 x3 x4 x5 x6)-
572 (vor_0_x_flipped x1 x2 x3 x4 x5 x6))`);;
574 (* STM changed to use vor_0_x_flipped instead of the definition *)
575 (* let nu_gamma_x = kepler_def(`nu_gamma_x x1 x2 x3 x4 x5 x6 = *)
578 (* (&2 * (gamma_x x1 x2 x3 x4 x5 x6))+ *)
579 (* (vor_0_x x1 x2 x3 x4 x5 x6)- *)
580 (* (vor_0_x x1 x5 x6 x4 x2 x3))`);; *)
581 let nu_gamma_x = kepler_def(`nu_gamma_x x1 x2 x3 x4 x5 x6 =
584 (&2 * (gamma_x x1 x2 x3 x4 x5 x6))+
585 (vor_0_x x1 x2 x3 x4 x5 x6)-
586 (vor_0_x_flipped x1 x2 x3 x4 x5 x6))`);;
588 let taunu_x = kepler_def(`taunu_x x1 x2 x3 x4 x5 x6 =
589 (sol_x x1 x2 x3 x4 x5 x6)*zeta*pt - (nu_x x1 x2 x3 x4 x5 x6)`);;
593 (* score for upright quarters in a quasi-regular octahedron.
594 I don't think I had a name for this specifically. *)
595 (* let octa_x = kepler_def(`octa_x x1 x2 x3 x4 x5 x6 = *)
597 (* (mu_upright_x x1 x2 x3 x4 x5 x6)+ *)
598 (* (mu_upright_x x1 x5 x6 x4 x2 x3))`);; *)
600 let octa_x = kepler_def(`octa_x x1 x2 x3 x4 x5 x6 =
602 (mu_upright_x x1 x2 x3 x4 x5 x6)+
603 (mu_flipped_x x1 x2 x3 x4 x5 x6))`);;
605 let sigmahat_x = kepler_def(`sigmahat_x x1 x2 x3 x4 x5 x6 =
606 let r234 = eta_x x2 x3 x4 in
607 let r456 = eta_x x4 x5 x6 in
608 let v0 = (if (sqrt2 <= r456) then
609 (vor_0_x x1 x2 x3 x4 x5 x6)
610 else if (sqrt2 <= r234) then
611 (vor_analytic_x x1 x2 x3 x4 x5 x6)
613 (gamma_x x1 x2 x3 x4 x5 x6)) in
614 let v1 = (if ((r234 <= sqrt2) /\ (r456 <= sqrt2)) then
615 max_real v0 (gamma_x x1 x2 x3 x4 x5 x6)
617 let v2 = (if (sqrt2 <= r234) then
618 max_real v1 (vor_analytic_x x1 x2 x3 x4 x5 x6)
620 let v3 = (if ((square (#2.6)) <= x4) /\ ((square (#2.2)) <= x1)
621 then max_real v2 (vor_0_x x1 x2 x3 x4 x5 x6)
623 let v4 = (if ((square (#2.7) <= x4))
625 max_real v3 (vor_0_x x1 x2 x3 x4 x5 x6)
629 max_real v4 (vor_analytic_x x1 x2 x3 x4 x5 x6)
632 let sigmahat_x' = kepler_def(`sigmahat_x' x1 x2 x3 x4 x5 x6 =
633 let r234 = eta_x x2 x3 x4 in
634 let r456 = eta_x x4 x5 x6 in
635 let P1 = sqrt2 <= r456 in
636 let P2 = sqrt2 <= r234 in
637 let P3 = square (#2.2) <= x1 in
638 let P4 = square (#2.6) <= x4 in
639 let P5 = square (#2.7) <= x4 in
640 if ~P1 /\ P2 /\ ~P5 /\ (~P3 \/ ~P4) then
641 vor_analytic_x x1 x2 x3 x4 x5 x6
642 else if ~P1 /\ ~P2 /\ ~P5 /\ (~P3 \/ ~P4) then
643 gamma_x x1 x2 x3 x4 x5 x6
644 else if ~P1 /\ ~P2 /\ P4 /\ (P3 \/ P5) then
645 max_real (gamma_x x1 x2 x3 x4 x5 x6) (vor_0_x x1 x2 x3 x4 x5 x6)
647 max_real (vor_analytic_x x1 x2 x3 x4 x5 x6) (vor_0_x x1 x2 x3 x4 x5 x6)`);;
649 let sigmahatpi_x = kepler_def(`sigmahatpi_x x1 x2 x3 x4 x5 x6 =
650 let r234 = eta_x x2 x3 x4 in
651 let r456 = eta_x x4 x5 x6 in
652 let piF = (#2.0)*(xiV) + (xi_gamma) in
653 let v0 = (if (sqrt2 <= r456) then
654 (vor_0_x x1 x2 x3 x4 x5 x6)
655 else if (sqrt2 <= r234) then
656 (vor_analytic_x x1 x2 x3 x4 x5 x6)
658 (gamma_x x1 x2 x3 x4 x5 x6)) in
659 let v1 = (if ((r234 <= sqrt2) /\ (r456 <= sqrt2)) then
660 max_real v0 (gamma_x x1 x2 x3 x4 x5 x6)
662 let v2 = (if (sqrt2 <= r234) then
663 max_real v1 (vor_analytic_x x1 x2 x3 x4 x5 x6)
665 let v3 = (if ((square (#2.6)) <= x4) /\ ((square (#2.2)) <= x2)
666 then max_real v2 (piF + (vor_0_x x1 x2 x3 x4 x5 x6))
668 let v4 = (if ((square (#2.7) <= x4))
670 max_real v3 (piF + (vor_0_x x1 x2 x3 x4 x5 x6) )
674 max_real v4 (vor_analytic_x x1 x2 x3 x4 x5 x6)
678 let tauhat_x = kepler_def(`tauhat_x x1 x2 x3 x4 x5 x6 =
679 (sol_x x1 x2 x3 x4 x5 x6)*zeta*pt - (sigmahat_x x1 x2 x3 x4 x5 x6)`);;
681 let tauhatpi_x = kepler_def(`tauhatpi_x x1 x2 x3 x4 x5 x6 =
682 (sol_x x1 x2 x3 x4 x5 x6)*zeta*pt - (sigmahatpi_x x1 x2 x3 x4 x5 x6)`);;
684 let pi_prime_tau = kepler_def
685 `pi_prime_tau k0 k1 k2 =
687 else if (k0=1) /\ (k1=1) /\ (k2=1) then (#0.0254)
688 else (#0.04683 + (&k0 + &(2*k2) - #3.0)*((#0.008)/(#3.0))
689 + (&k2) * (#0.0066))`;;
691 let pi_prime_sigma = kepler_def
692 `pi_prime_sigma k0 (k1:num) k2 =
694 else if (k0=1) /\ (k2=1) then (#0.009)
695 else (& (k0 + 2*k2))*((#0.008)/(#3.0)) + (&k2)* (#0.009)`;;
699 (* ------------------------------------------------------------------ *)
701 (* ------------------------------------------------------------------ *)
703 (* We are swithing from real3 to real^3. *)
706 (* For now I'm keeping these definitons (at least the names, the definitions*)
707 (* themselves are changed radically), but it might be better to just get *)
708 (* rid of most of them. *)
711 (* let dot3 = new_definition `dot3 (v:real^3) w = v dot w`;; *)
714 (* let norm3 = new_definition `norm3 (v:real^3) = norm v`;; *)
716 (* deprecated *) let d3 = new_definition `d3 (v:real^3) w = dist(v,w)`;;
718 (* No need for this one. "basis" does something similar. *)
720 let dirac_delta = new_definition `dirac_delta (i:num) =
721 (\j. if (i=j) then (&1) else (&0))`;;
725 (* deprecated, use vector, instead *)
726 let mk_vec3 = new_definition `mk_vec3 a b c = vector[a; b; c]`;;
728 let real3_of_triple = new_definition `real3_of_triple (a,b,c) = (mk_vec3 a b c):real^3`;;
730 let triple_of_real3 = new_definition `triple_of_real3 (v:real^3) =
734 (* let orig3 = new_definition `orig3 = (vec 0):real^3`;; *)
736 (* ------------------------------------------------------------------ *)
737 (* Cross diagonal and Enclosed *)
738 (* ------------------------------------------------------------------ *)
742 (* find point in euclidean 3 space atdistance a b c
746 v3 = mk_vec3 v3_1 v3_2 0;
748 let findpoint = kepler_def `findpoint a b c y4 v3_1 v3_2 sgn =
749 let y5 = sqrt (v3_1*v3_1 + v3_2*v3_2) in
750 let w1 = (a*a + y4*y4 - b*b)/((&2)*y4) in
751 let w2 = (a*a + y5*y5 -c*c - (&2)*w1*v3_1)/((&2)*v3_2) in
752 let w3 = sgn* (sqrt(a*a - w1*w1 - w2*w2)) in
755 let enclosed = kepler_def `enclosed y1 y2 y3 y4 y5 y6 y1' y2' y3' =
756 let v1:real^3 = mk_vec3 (&0) (&0) (&0) in
757 let v2:real^3 = mk_vec3 y4 (&0) (&0) in
758 let a = ((y5*y5) + (y4*y4) - (y6*y6))/((&2)*y4) in
759 let b = sqrt((y5*y5) - (a*a)) in
760 let v3:real^3 = mk_vec3 a b (&0) in
761 let v4:real^3 = findpoint y3 y2 y1 y4 a b (#1.0) in
762 let v5 = findpoint y3' y2' y1' y4 a b (--(#1.0)) in
765 let cross_diag = kepler_def `cross_diag y1 y2 y3 y4 y5 y6 y7 y8 y9=
766 enclosed y1 y5 y6 y4 y2 y3 y7 y8 y9`;;
768 let cross_diag_x = kepler_def `cross_diag_x x1 x2 x3 x4 x5 x6 x7 x8 x9=
769 cross_diag (sqrt x1) (sqrt x2) (sqrt x3) (sqrt x4) (sqrt x5)
770 (sqrt x6) (sqrt x7) (sqrt x8) (sqrt x9)`;;
772 (* ------------------------------------------------------------------ *)
773 (* Definitions of Affine Geometry (from BLUEPRINT : Trigonometry ) *)
774 (* ------------------------------------------------------------------ *)
776 (* Notes on unique existence of definitions:
779 Exists uniquely on all nonempty subsets of N.
780 Can be uniquely extended to all subsets by defining min_num {} = 0
783 This is the derivative of a function of a real variable. Its domain is more difficult to describe.
786 Exists uniquely on all finite subsets of R3.
787 It will only be used on finite sets.
788 Can be uniquely extended to all subsets by defining aff S = S, when S is infinite
792 Exists uniquely on all nonempty finite sets of ordered pairs of real numbers
793 Can be uniquely extended to all sets of ordered pairs by setting min_polar X = ( &0, &0 ), when X is empty or infinite.
794 This definition actually holds for some infinite sets, but I never use it, except on finite sets, so you are free
795 to redefine it to be (&0,&0) on infinite sets.
798 iter is uniquely defined with no domain conditions, no preconditions
801 Exists uniquely when the stated non collinearity preconditions are met: ~(collinear {v, w, w1}) /\ ~(collinear {v, w, w2})
802 (The preconditions azim_hyp, orthonormality, and e3 normalization are not domain constraints, because they do not restrict v w w1 w2.)
803 Can be uniquely extended to all cases, by setting azim w1 w2 w3 w4 = &0, when the non-collinearity preconditions are not met.
807 azim_cycle W v w p exists uniquely under the preconditions (W p) /\ (cyclic_set W v w).
808 It is only used when the preconditions are met.
809 Can be uniquely extended to all cases, by setting azim_cycle W v w p = p, when these two preconditions are not met.
816 (* from convex.ml: *)
818 let affine = new_definition
819 `affine s <=> !x y u v. x IN s /\ y IN s /\ (u + v = &1)
820 ==> (u % x + v % y) IN s`;;
821 let convex = new_definition
823 !x y u v. x IN s /\ y IN s /\ &0 <= u /\ &0 <= v /\ (u + v = &1)
824 ==> (u % x + v % y) IN s`;;
826 (* aff is deprecated *)
827 let aff = new_definition `aff = ( hull ) affine`;;
829 let lin_combo = new_definition `lin_combo V f = vsum V (\v. f v % (v:real^N))`;;
831 (* Fix "sum" because Harrison's interface is too special in analysis.ml *)
833 reduce_interface("sum",`sum:(num->bool)->(num->real)->real`);;
834 reduce_interface("sum",`psum:(num#num)->(num->real)->real`);;
835 let remove_overload sym =
836 let overload_skeletons = filter ((<>)sym o fst) (!the_overload_skeletons) in
837 the_overload_skeletons := overload_skeletons;;
838 remove_overload "sum";;
839 make_overloadable "sum" `:A->(B->real)->real`;;
840 overload_interface("sum",`sum:(A->bool)->(A->real)->real`);;
841 overload_interface("sum",`psum:(num#num)->(num->real)->real`);;
843 let affsign = new_definition `affsign sgn s t (v:real^A) = (?f.
844 (v = lin_combo (s UNION t) f) /\ (!w. t w ==> sgn (f w)) /\ (sum (s UNION t) f = &1))`;;
847 let sgn_gt = new_definition `sgn_gt = (\t. (&0 < t))`;;
848 let sgn_ge = new_definition `sgn_ge = (\t. (&0 <= t))`;;
849 let sgn_lt = new_definition `sgn_lt = (\t. (t < &0))`;;
850 let sgn_le = new_definition `sgn_le = (\t. (t <= &0))`;;
852 (* conv is deprecated. Use `convex hull S` instead *)
854 let conv = new_definition `conv S:real^A->bool = affsign sgn_ge {} S`;;
855 let conv0 = new_definition `conv0 S:real^A->bool = affsign sgn_gt {} S`;;
856 let cone = new_definition `cone v S:real^A->bool = affsign sgn_ge {v} S`;;
857 let cone0 = new_definition `cone0 v S:real^A->bool = affsign sgn_gt {v} S`;;
861 let semiconvex = new_definition
862 `semiconvex sgn s t <=>
863 !x y z u v w. x IN (affine hull s) /\ y IN t /\ z IN t /\ sgn v /\ sgn w /\ (u + v + w = &1)
864 ==> (u % x + v % y + w % z) IN t`;;
868 let aff_gt_def = new_definition `aff_gt = affsign sgn_gt`;;
869 let aff_ge_def = new_definition `aff_ge = affsign sgn_ge`;;
870 let aff_lt_def = new_definition `aff_lt = affsign sgn_lt`;;
871 let aff_le_def = new_definition `aff_le = affsign sgn_le`;;
874 (* conv is deprecated. Use `convex hull S` instead *)
876 (* Vuong Quyen has pointed out that the definition of aff_insert
879 New definitions are based on Multivariate/convex.ml.
887 let voronoi = new_definition `voronoi v S = { x | !w. ((S w) /\ ~(w=v)) ==> (dist( x, v) < dist( x, w)) }`;;
889 let voronoi_le = new_definition `voronoi_le v S = { x | !w. ((S w) /\ ~(w=v)) ==> (dist( x, v) <= dist( x, w)) }`;;
895 let line = new_definition `line x = (?v w. ~(v =w) /\ (x = affine hull {v,w}))`;;
897 (* Done in Harrison's Multilinear/vectors.ml (Feb 2008 release only) : let collinear = new_definition `collinear S = (?x. line x /\ S SUBSET x)`;; *)
898 (* repeat of definition for 2.20 version *)
900 let collinear = new_definition
901 `collinear s <=> ?u. !x y. x IN s /\ y IN s ==> ?c. x - y = c % u`;;
907 let plane = new_definition `plane x = (?u v w. ~(collinear {u,v,w}) /\ (x = affine hull {u,v,w}))`;;
908 let closed_half_plane = new_definition `closed_half_plane x = (?u v w. ~(collinear {u,v,w}) /\ (x = aff_ge {u,v} {w}))`;;
909 let open_half_plane = new_definition `open_half_plane x = (?u v w. ~(collinear {u,v,w}) /\ (x = aff_gt {u,v} {w}))`;;
910 let coplanar = new_definition `coplanar S = (?x. plane x /\ S SUBSET x)`;;
911 let closed_half_space = new_definition `closed_half_space x = (?u v w w'. ~(coplanar {u,v,w,w'}) /\ (x = aff_ge {u,v,w} {w'}))`;;
912 let open_half_space = new_definition `open_half_space x = (?u v w w'. ~(coplanar {u,v,w,w'}) /\ (x = aff_gt {u,v,w} {w'}))`;;
915 let bis = new_definition `bis u v = {x | dist(x,u) = dist(x,v)}`;;
918 let bis_le = new_definition `bis_le u v = {x | dist(x,u) <= dist(x,v) }`;;
919 let bis_lt = new_definition `bis_lt u v = {x | dist(x,u) < dist(x,v) }`;;
922 let circumcenter = new_definition `circumcenter S = @v. ( (affine hull S) v /\ (?c. !w. (S w) ==> (c = dist(v,w))))`;;
926 let radV = new_definition `radV S = @c. !w. (S w) ==> (c = dist(circumcenter S,w))`;;
930 let orientation = new_definition `orientation S v sgn = affsign sgn (S DIFF {v}) {v} (circumcenter S)`;;
934 let arcV = new_definition `arcV u v w = acs (( (v - u) dot (w - u))/((norm (v-u)) * (norm (w-u))))`;;
936 let cross = new_definition `cross u v = let (x,y,z) = triple_of_real3 u in
937 let (x',y',z') = triple_of_real3 v in
938 (real3_of_triple (y*z' - y'*z, z*x' - z'*x, x*y' - x'*y))`;;
940 let dihV = new_definition `dihV w0 w1 w2 w3 =
944 let vap = ( vc dot vc) % va - ( va dot vc) % vc in
945 let vbp = ( vc dot vc) % vb - ( vb dot vc) % vc in
946 arcV (vec 0) vap vbp`;;
948 (* conventional ordering on variables *)
950 let ylist = new_definition `ylist w0 w1 w2 w3 =
951 ((dist (w0, w1)),(dist( w0, w2)),(dist( w0, w3)),(dist( w2, w3)),(dist( w1, w3)),(dist( w1, w2)))`;;
953 let xlist = new_definition `xlist w0 w1 w2 w3 =
954 let (y1,y2,y3,y4,y5,y6) = ylist w0 w1 w2 w3 in
955 (y1 pow 2, y2 pow 2, y3 pow 2, y4 pow 2, y5 pow 2, y6 pow 2)`;;
957 let euler_p = new_definition `euler_p v0 v1 v2 v3 =
958 (let (y1,y2,y3,y4,y5,y6) = ylist v0 v1 v2 v3 in
962 y1*y2*y3 + y1*( w2 dot w3) + y2*( w3 dot w1) + y3*( w1 dot w2))`;;
964 (* polar coordinates *)
966 let radius = new_definition `radius x y = sqrt((x pow 2) + (y pow 2))`;;
967 let polar_angle = new_definition `polar_angle x y =
968 let theta = atn2(x,y) in
969 if (theta < &0) then (theta + (&2 * pi)) else theta`;;
970 let polar_c = new_definition `polar_c x y = (radius x y,polar_angle x y)`;;
972 let less_polar = new_definition `less_polar (x,y) (x',y') =
973 let (r,theta) = polar_c x y in
974 let (r',theta') = polar_c x' y' in
975 (theta < theta') \/ ((theta =theta') /\ (r < r'))`;;
977 let min_polar = new_definition `min_polar V = ( @ u. V u /\ (!w. V w /\ ~(u = w) ==> (less_polar u w)))`;;
979 let polar_cycle = new_definition `polar_cycle V v =
980 let W = {u | V u /\ less_polar v u} in
981 if (W = EMPTY) then min_polar V else min_polar W`;;
983 (* iterates of a function must be done already, but I don't know where *)
985 let iter_spec = prove(`?iter. !f u:A. (iter 0 f u = u) /\ (!n. (iter (SUC n) f u = f(iter n f u)))`,
986 (SUBGOAL_THEN `?g. !f (u:A). (g f u 0 = u) /\ (!n. (g f u (SUC n) = f (g f u n)))` CHOOSE_TAC) THENL
987 ([REWRITE_TAC[GSYM SKOLEM_THM;num_RECURSION_STD];REWRITE_TAC[]]) THEN
988 (EXISTS_TAC `\ (i:num) (f:A->A) (u:A). (g f u i):A`) THEN
990 (ASM_REWRITE_TAC[]));;
992 let iter = new_specification ["iter"] iter_spec;;
995 let polar_power_spec = prove(`?fn. !V v. (fn V v 0 = v ) /\ (!n. (fn V v (SUC n) = polar_cycle V (fn V v n)))`,
996 (REWRITE_TAC[GSYM SKOLEM_THM;num_RECURSION_STD]));;
998 let polar_power = new_specification ["polar_power"] polar_power_spec;;
1001 let orthonormal = new_definition `orthonormal e1 e2 e3 =
1002 (( e1 dot e1 = &1) /\ (e2 dot e2 = &1) /\ ( e3 dot e3 = &1) /\
1003 ( e1 dot e2 = &0) /\ ( e1 dot e3 = &0) /\ ( e2 dot e3 = &0) /\
1004 (&0 < (cross e1 e2) dot e3))`;;
1006 (* spherical coordinates *)
1007 let azim_hyp_def = new_definition `azim_hyp = (!v w w1 w2. ?theta. !e1 e2 e3. ?psi h1 h2 r1 r2.
1008 ~(collinear {v, w, w1}) /\ ~(collinear {v, w, w2}) /\
1009 (orthonormal e1 e2 e3) /\ ((dist( w, v)) % e3 = (w - v)) ==>
1010 ((&0 <= theta) /\ (theta < &2 * pi) /\ (&0 < r1) /\ (&0 < r2) /\
1011 (w1 - v = (r1 * cos(psi)) % e1 + (r1 * sin(psi)) % e2 + h1 % (w-v)) /\
1012 (w2 - v = (r2 * cos(psi + theta)) % e1 + (r2 * sin(psi + theta)) % e2 + h2 % (w-v))))`;;
1014 let azim_spec = prove(`?theta. !v w w1 w2 e1 e2 e3. ?psi h1 h2 r1 r2.
1016 ~(collinear {v, w, w1}) /\ ~(collinear {v, w, w2}) /\
1017 (orthonormal e1 e2 e3) /\ ((dist( w, v)) % e3 = (w - v)) ==>
1018 ((&0 <= theta v w w1 w2) /\ (theta v w w1 w2 < &2 * pi) /\ (&0 < r1) /\ (&0 < r2) /\
1019 (w1 - v = (r1 * cos(psi)) % e1 + (r1 * sin(psi)) % e2 + h1 % (w-v)) /\
1020 (w2 - v = (r2 * cos(psi + theta v w w1 w2)) % e1 + (r2 * sin(psi + theta v w w1 w2)) % e2 + h2 % (w-v)))`,
1021 (REWRITE_TAC[GSYM SKOLEM_THM;GSYM RIGHT_IMP_EXISTS_THM;GSYM RIGHT_IMP_FORALL_THM]) THEN
1022 (REWRITE_TAC[azim_hyp_def]) THEN
1023 (REPEAT STRIP_TAC) THEN
1024 (ASM_REWRITE_TAC[RIGHT_IMP_EXISTS_THM]));;
1026 let azim_def = new_specification ["azim"] azim_spec;;
1029 let cyclic_set = new_definition `cyclic_set W v w =
1030 (~(v=w) /\ (FINITE W) /\ (!p q h. W p /\ W q /\ (p = q + h % (v - w)) ==> (p=q)) /\
1031 (W INTER (affine hull {v,w}) = EMPTY))`;;
1036 let azim_cycle_hyp_def = new_definition `azim_cycle_hyp =
1037 (?sigma. !W proj v w e1 e2 e3 p.
1039 (cyclic_set W v w) /\ ((dist( v ,w)) % e3 = (w-v)) /\
1040 (orthonormal e1 e2 e3) /\
1041 (!u x y. (proj u = (x,y)) <=> (?h. (u = v + x % e1 + y % e2 + h % e3))) ==>
1042 (proj (sigma W v w p) = polar_cycle (IMAGE proj W) (proj p)))`;;
1044 let azim_cycle_spec = prove(`?sigma. !W proj v w e1 e2 e3 p.
1045 (azim_cycle_hyp) ==> ( (W p) /\
1046 (cyclic_set W v w) /\ ((dist( v ,w)) % e3 = (w-v)) /\
1047 (orthonormal e1 e2 e3) /\
1048 (!u x y. (proj u = (x,y)) <=> (?h. (u = v + x % e1 + y % e2 + h % e3)))) ==> (proj (sigma W v w p) = polar_cycle (IMAGE proj W) (proj p))`,
1049 (REWRITE_TAC[GSYM RIGHT_IMP_EXISTS_THM;GSYM RIGHT_IMP_FORALL_THM]) THEN
1050 (REWRITE_TAC[azim_cycle_hyp_def])
1053 let azim_cycle_def = new_specification ["azim_cycle"] azim_cycle_spec;;
1056 (* ------------------------------------------------------------------ *)
1057 (* Definitions from the Collection in Elementary Geometry *)
1058 (* ------------------------------------------------------------------ *)
1062 let cayleyR = new_definition `cayleyR x12 x13 x14 x15 x23 x24 x25 x34 x35 x45 =
1063 -- (x14*x14*x23*x23) + &2 *x14*x15*x23*x23 - x15*x15*x23*x23 + &2 *x13*x14*x23*x24 - &2 *x13*x15*x23*x24 - &2 *x14*x15*x23*x24 +
1064 &2 *x15*x15*x23*x24 - x13*x13*x24*x24 + &2 *x13*x15*x24*x24 - x15*x15*x24*x24 - &2 *x13*x14*x23*x25 +
1065 &2 *x14*x14*x23*x25 + &2 *x13*x15*x23*x25 - &2 *x14*x15*x23*x25 + &2 *x13*x13*x24*x25 - &2 *x13*x14*x24*x25 - &2 *x13*x15*x24*x25 +
1066 &2 *x14*x15*x24*x25 - x13*x13*x25*x25 + &2 *x13*x14*x25*x25 - x14*x14*x25*x25 + &2 *x12*x14*x23*x34 - &2 *x12*x15*x23*x34 -
1067 &2 *x14*x15*x23*x34 + &2 *x15*x15*x23*x34 + &2 *x12*x13*x24*x34 - &2 *x12*x15*x24*x34 - &2 *x13*x15*x24*x34 + &2 *x15*x15*x24*x34 +
1068 &4 *x15*x23*x24*x34 - &2 *x12*x13*x25*x34 - &2 *x12*x14*x25*x34 + &4 *x13*x14*x25*x34 + &4 *x12*x15*x25*x34 - &2 *x13*x15*x25*x34 - &2 *x14*x15*x25*x34 -
1069 &2 *x14*x23*x25*x34 - &2 *x15*x23*x25*x34 - &2 *x13*x24*x25*x34 - &2 *x15*x24*x25*x34 + &2 *x13*x25*x25*x34 + &2 *x14*x25*x25*x34 -
1070 x12*x12*x34*x34 + &2 *x12*x15*x34*x34 - x15*x15*x34*x34 + &2 *x12*x25*x34*x34 + &2 *x15*x25*x34*x34 -
1071 x25*x25*x34*x34 - &2 *x12*x14*x23*x35 + &2 *x14*x14*x23*x35 + &2 *x12*x15*x23*x35 - &2 *x14*x15*x23*x35 - &2 *x12*x13*x24*x35 +
1072 &4 *x12*x14*x24*x35 - &2 *x13*x14*x24*x35 - &2 *x12*x15*x24*x35 + &4 *x13*x15*x24*x35 - &2 *x14*x15*x24*x35 - &2 *x14*x23*x24*x35 - &2 *x15*x23*x24*x35 +
1073 &2 *x13*x24*x24*x35 + &2 *x15*x24*x24*x35 + &2 *x12*x13*x25*x35 - &2 *x12*x14*x25*x35 - &2 *x13*x14*x25*x35 + &2 *x14*x14*x25*x35 +
1074 &4 *x14*x23*x25*x35 - &2 *x13*x24*x25*x35 - &2 *x14*x24*x25*x35 + &2 *x12*x12*x34*x35 - &2 *x12*x14*x34*x35 - &2 *x12*x15*x34*x35 +
1075 &2 *x14*x15*x34*x35 - &2 *x12*x24*x34*x35 - &2 *x15*x24*x34*x35 - &2 *x12*x25*x34*x35 - &2 *x14*x25*x34*x35 + &2 *x24*x25*x34*x35 -
1076 x12*x12*x35*x35 + &2 *x12*x14*x35*x35 - x14*x14*x35*x35 + &2 *x12*x24*x35*x35 + &2 *x14*x24*x35*x35 -
1077 x24*x24*x35*x35 + &4 *x12*x13*x23*x45 - &2 *x12*x14*x23*x45 - &2 *x13*x14*x23*x45 - &2 *x12*x15*x23*x45 - &2 *x13*x15*x23*x45 +
1078 &4 *x14*x15*x23*x45 + &2 *x14*x23*x23*x45 + &2 *x15*x23*x23*x45 - &2 *x12*x13*x24*x45 + &2 *x13*x13*x24*x45 + &2 *x12*x15*x24*x45 -
1079 &2 *x13*x15*x24*x45 - &2 *x13*x23*x24*x45 - &2 *x15*x23*x24*x45 - &2 *x12*x13*x25*x45 + &2 *x13*x13*x25*x45 + &2 *x12*x14*x25*x45 -
1080 &2 *x13*x14*x25*x45 - &2 *x13*x23*x25*x45 - &2 *x14*x23*x25*x45 + &4 *x13*x24*x25*x45 + &2 *x12*x12*x34*x45 - &2 *x12*x13*x34*x45 -
1081 &2 *x12*x15*x34*x45 + &2 *x13*x15*x34*x45 - &2 *x12*x23*x34*x45 - &2 *x15*x23*x34*x45 - &2 *x12*x25*x34*x45 - &2 *x13*x25*x34*x45 + &2 *x23*x25*x34*x45 +
1082 &2 *x12*x12*x35*x45 - &2 *x12*x13*x35*x45 - &2 *x12*x14*x35*x45 + &2 *x13*x14*x35*x45 - &2 *x12*x23*x35*x45 - &2 *x14*x23*x35*x45 -
1083 &2 *x12*x24*x35*x45 - &2 *x13*x24*x35*x45 + &2 *x23*x24*x35*x45 + &4 *x12*x34*x35*x45 - x12*x12*x45*x45 + &2 *x12*x13*x45*x45 -
1084 x13*x13*x45*x45 + &2 *x12*x23*x45*x45 + &2 *x13*x23*x45*x45 - x23*x23*x45*x45`;;
1089 let packing = new_definition `packing S = (!u v. S u /\ S v /\ ~(u = v) ==> (&2 <= dist( u, v)))`;;
1093 let wedge = new_definition (`wedge v1 v2 w1 w2 =
1097 let n = cross z u1 in
1099 if (aff_ge {v1,v2} {w1} w2) then {} else
1100 if (aff_lt {v1,v2} {w1} w2) then aff_gt {v1,v2,w1} {n} else
1101 if (d > &0) then aff_gt {v1,v2} {w1,w2} else
1102 (:real^3) DIFF aff_ge {v1,v2} {w1,w2}`);;
1104 (* ------------------------------------------------------------------ *)
1105 (* Measure and Volume, following Nguyen Tat Thang *)
1106 (* ------------------------------------------------------------------ *)
1108 let sphere= new_definition`sphere x=(?(v:real^3)(r:real). (r> &0)/\ (x={w:real^3 | norm (w-v)= r}))`;;
1110 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)}`;;
1112 let circular_cone =new_definition `circular_cone (V:real^3-> bool)=
1113 (? (v,w:real^3)(r:real). V= c_cone (v,w,r))`;;
1115 let NULLSET_RULES,NULLSET_INDUCT,NULLSET_CASES =
1116 new_inductive_definition
1117 `(!P. ((plane P)\/ (sphere P) \/ (circular_cone P)) ==> NULLSET P) /\
1118 !(s:real^3->bool) t. (NULLSET s /\ NULLSET t) ==> NULLSET (s UNION t)`;;
1120 let null_equiv = new_definition `null_equiv (s,t :real^3->bool)=(? (B:real^3-> bool). NULLSET B /\
1121 ((s DIFF t) UNION (t DIFF s)) SUBSET B)`;;
1124 let normball = new_definition `normball x r = { y:real^A | dist(y,x) < r}`;;
1126 let radial = new_definition `radial r x C <=> (C SUBSET normball x r) /\ (!u. (x+u) IN C ==> (!t.(t> &0) /\ (t* norm u < r)==>(x+ t % u) IN C))`;;
1128 let eventually_radial = new_definition `eventually_radial x C <=> (?r. (r> &0) /\ radial r x (C INTER normball x r))`;;
1130 let solid_triangle = new_definition `solid_triangle v0 S r = normball v0 r INTER cone v0 S`;;
1132 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)}`;;
1136 let rcone_ge = new_definition `rcone_ge = rconesgn ( >= )`;;
1137 let rcone_gt = new_definition `rcone_gt = rconesgn ( > )`;;
1138 let rcone_lt = new_definition `rcone_lt = rconesgn ( < )`;;
1139 let rcone_eq = new_definition `rcone_eq = rconesgn ( = )`;;
1141 let scale = new_definition `scale (t:real^3) (u:real^3) = vector[t$1 * u$1; t$2 * u$2; t$3 * u$3]`;;
1143 let ellipsoid = new_definition `ellipsoid t r = IMAGE (scale t) (normball(vec 0)r)`;;
1145 let conic_cap = new_definition `conic_cap v0 v1 r a = normball v0 r INTER rcone_gt v0 v1 a`;;
1147 let frustum = new_definition `frustum v0 v1 h1 h2 a = { y | rcone_gt v0 v1 a y /\
1148 let d = (y - v0) dot (v1 - v0) in
1149 let n = norm(v1 - v0) in
1150 (h1*n < d /\ d < h2*n)}`;;
1152 let frustt = new_definition `frustt v0 v1 h a = frustum v0 v1 (&0) h a`;;
1154 let rect = new_definition `rect (a:real^3) (b:real^3) = {(v:real^3) | !i. ( a$i < v$i /\ v$i < b$i )}`;;
1157 let is_tetrahedron = new_definition `is_tetrahedron S = ?v0 v1 v2 v3. (S = conv0 {v0,v1,v2,v3})`;;
1160 let primitive = new_definition `primitive (C:real^3->bool) =
1161 ((?v0 v1 v2 v3 r. (C = solid_triangle v0 {v1,v2,v3} r)) \/
1162 (?v0 v1 v2 v3. (C = conv0 {v0,v1,v2,v3})) \/
1163 (?v0 v1 v2 v3 h a. (C = frustt v0 v1 h a INTER wedge v0 v1 v2 v3)) \/
1164 (?v0 v1 v2 v3 r c. (C = conic_cap v0 v1 r c INTER wedge v0 v1 v2 v3)) \/
1165 (?a b. (C = rect a b)) \/
1166 (?t r. (C = ellipsoid t r)) \/
1167 (?v0 v1 v2 v3 r. (C = normball v0 r INTER wedge v0 v1 v2 v3)))`;;
1169 let MEASURABLE_RULES,MEASURABLE_INDUCT,MEASURABLE_CASES =
1170 new_inductive_definition
1171 `(!C. primitive C ==> measurable C) /\
1172 ( !Z. NULLSET Z ==> measurable Z) /\
1173 (!X t. measurable X ==> (measurable (IMAGE (scale t) X))) /\
1174 (!X v. measurable X ==> (measurable (IMAGE ((+) v) X))) /\
1175 ( !(s:real^3->bool) t. (measurable s /\ measurable t) ==> measurable (s UNION t)) /\
1176 ( !(s:real^3->bool) t. (measurable s /\ measurable t) ==> measurable (s INTER t)) /\
1177 ( !(s:real^3->bool) t. (measurable s /\ measurable t) ==> measurable (s DIFF t))
1181 let SDIFF = new_definition `SDIFF X Y = (X DIFF Y) UNION (Y DIFF X)`;;
1184 let vol_solid_triangle = new_definition `vol_solid_triangle v0 v1 v2 v3 r =
1185 let a123 = dihV v0 v1 v2 v3 in
1186 let a231 = dihV v0 v2 v3 v1 in
1187 let a312 = dihV v0 v3 v1 v2 in
1188 (a123 + a231 + a312 - pi)*(r pow 3)/(&3)`;;
1190 let vol_frustt_wedge = new_definition `vol_frustt_wedge v0 v1 v2 v3 h a =
1191 (azim v0 v1 v2 v3)*(h pow 3)*(&1/(a*a) - &1)/(&6)`;;
1193 (* volume of intersection of conic cap and wedge *)
1194 let vol_conic_cap_wedge = new_definition `vol_conic_cap_wedge v0 v1 v2 v3 r c =
1195 (azim v0 v1 v2 v3)*(&1 - c)*(r pow 3)/(&3)`;;
1198 let vol_conv = new_definition `vol_conv v1 v2 v3 v4 =
1199 let x12 = dist(v1,v2) pow 2 in
1200 let x13 = dist(v1,v3) pow 2 in
1201 let x14 = dist(v1,v4) pow 2 in
1202 let x23 = dist(v2,v3) pow 2 in
1203 let x24 = dist(v2,v4) pow 2 in
1204 let x34 = dist(v3,v4) pow 2 in
1205 sqrt(delta_x x12 x13 x14 x34 x24 x34)/(&12)`;;
1207 let vol_rect = new_definition `vol_rect a b =
1208 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`;;
1210 let vol_ball_wedge = new_definition `vol_ball_wedge v0 v1 v2 v3 r =
1211 (azim v0 v1 v2 v3)*(&2)*(r pow 3)/(&3)`;;
1214 let volume_props = new_definition `volume_props (vol:(real^3->bool)->real) =
1215 ( (!C. vol C >= &0) /\
1216 (!Z. NULLSET Z ==> (vol Z = &0)) /\
1217 (!X Y. measurable X /\ measurable Y /\ NULLSET (SDIFF X Y) ==> (vol X = vol Y)) /\
1218 (!X t. (measurable X) /\ (measurable (IMAGE (scale t) X)) ==> (vol (IMAGE (scale t) X) = abs(t$1 * t$2 * t$3)*vol(X))) /\
1219 (!X v. measurable X ==> (vol (IMAGE ((+) v) X) = vol X)) /\
1220 (!v0 v1 v2 v3 r. (r > &0) /\ (~(collinear {v0,v1,v2})) /\ ~(collinear {v0,v1,v3}) ==> vol (solid_triangle v0 {v1,v2,v3} r) = vol_solid_triangle v0 v1 v2 v3 r) /\
1221 (!v0 v1 v2 v3. vol(conv0 {v0,v1,v2,v3}) = vol_conv v0 v1 v2 v3) /\
1222 (!v0 v1 v2 v3 h a. ~(collinear {v0,v1,v2}) /\ ~(collinear {v0,v1,v3}) /\ (h >= &0) /\ (a > &0) /\ (a <= &1) ==> vol(frustt v0 v1 h a INTER wedge v0 v1 v2 v3) = vol_frustt_wedge v0 v1 v2 v3 h a) /\
1223 (!v0 v1 v2 v3 r c. ~(collinear {v0,v1,v2}) /\ ~(collinear {v0,v1,v3}) /\ (r >= &0) /\ (c >= -- (&1)) /\ (c <= &1) ==> (vol(conic_cap v0 v1 r c INTER wedge v0 v1 v2 v3) = vol_conic_cap_wedge v0 v1 v2 v3 r c)) /\
1224 (!(a:real^3) (b:real^3). vol(rect a b) = vol_rect a b) /\
1225 (!v0 v1 v2 v3 r. ~(collinear {v0,v1,v2}) /\ ~(collinear {v0,v1,v3}) /\ (r >= &0) ==> (vol(normball v0 r INTER wedge v0 v1 v2 v3) = vol_ball_wedge v0 v1 v2 v3 r)))`;;
1227 let vol_def = new_definition `vol = ( @ ) volume_props`;;
1229 let solid = new_definition `solid x C = (@s. ?c. (c > &0) /\ (!r. (r > &0) /\ (r < c) ==>
1230 (s = &3 * vol(C INTER normball x r)/(r pow 3)))) `;;
1232 let sovo = new_definition `sovo x C (v,s) = v* vol(C) + s* solid x C`;;
1234 let phivo = new_definition `phivo h t (v,s) = v*t*h*(t+h)/(&6) + s`;;
1236 let avo = new_definition `avo h t l= (&1 - h/t)*(phivo h t l - phivo t t l)`;;
1238 let ortho0 = new_definition `ortho0 x v1 v2 v3 = conv0 {x,x+v1,x+v1+v2,x+v1+v2+v3}`;;
1240 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))`;;
1242 let rogers = new_definition `rogers v0 v1 v2 v3 c =
1243 let w = (&1/ (&2)) % (v0 + v1) in
1244 let p = circumcenter {v0,v1,v2} in
1245 let q = make_point v0 v1 v2 v3 c c c in
1248 let rogers0 = new_definition `rogers0 v0 v1 v2 v3 c =
1249 let w = (&1/ (&2)) % (v0 + v1) in
1250 let p = circumcenter {v0,v1,v2} in
1251 let q = make_point v0 v1 v2 v3 c c c in
1254 let abc_param = new_definition `abc_param v0 v1 v2 c =
1255 let a = (&1/(&2)) * dist(v0,v1) in
1256 let b = radV {v0,v1,v2} in
1259 (* ------------------------------------------------------------------ *)
1260 (* Format of inequalities in the archive. *)
1261 (* ------------------------------------------------------------------ *)
1263 let ineq = kepler_def `ineq (bounds:(real#real#real) list) (A:bool) =
1264 ((!(a,x,b). ((MEM (a,x,b) bounds) ==> ((a <= x) /\ (x <= b)))) /\ A)`;;
1266 let all_forall bod =
1267 let mk_forall = mk_binder "!" in
1268 itlist (curry mk_forall) (frees bod) bod;;