7 let atn2 = new_definition(`atn2(x,y) =
8 if ( abs y < x ) then atn(y / x) else
9 (if (&0 < y) then ((pi / &2) - atn(x / y)) else
10 (if (y < &0) then (-- (pi/ &2) - atn (x / y)) else ( pi )))`);;
12 (* ------------------------------------------------------------------ *)
14 let sqrt8 = new_definition (`sqrt8 = sqrt (&8) `);;
15 let sqrt2 = new_definition (`sqrt2 = sqrt (&2) `);;
17 let pi_rt18 = new_definition(`pi_rt18= pi/(sqrt (&18))`);;
20 (* ------------------------------------------------------------------ *)
21 (* This polynomial is essentially the Cayley-Menger determinant. *)
22 (* ------------------------------------------------------------------ *)
23 let delta_x = kepler_def (`delta_x x1 x2 x3 x4 x5 x6 =
24 x1*x4*(--x1 + x2 + x3 -x4 + x5 + x6) +
25 x2*x5*(x1 - x2 + x3 + x4 -x5 + x6) +
26 x3*x6*(x1 + x2 - x3 + x4 + x5 - x6)
27 -x2*x3*x4 - x1*x3*x5 - x1*x2*x6 -x4*x5*x6`);;
29 (* ------------------------------------------------------------------ *)
30 (* The partial derivative of delta_x with respect to x4. *)
31 (* ------------------------------------------------------------------ *)
33 let delta_x4 = kepler_def(`delta_x4 x1 x2 x3 x4 x5 x6
34 = -- x2* x3 - x1* x4 + x2* x5
35 + x3* x6 - x5* x6 + x1* (-- x1 + x2 + x3 - x4 + x5 + x6)`);;
37 let delta_x6 = kepler_def(`delta_x6 x1 x2 x3 x4 x5 x6
38 = -- x1 * x2 - x3*x6 + x1 * x4
39 + x2* x5 - x4* x5 + x3*(-- x3 + x1 + x2 - x6 + x4 + x5)`);;
41 (* ------------------------------------------------------------------ *)
43 (* ------------------------------------------------------------------ *)
47 `u_x x1 x2 x3 = (--(x1*x1+x2*x2+x3*x3)) +
48 (&2) * (x1*x2+x2*x3+x3*x1)`);;
51 let ups_x = kepler_def(`ups_x x1 x2 x6 =
52 --x1*x1 - x2*x2 - x6*x6
53 + &2 *x1*x6 + &2 *x1*x2 + &2 *x2*x6`);;
56 let eta_x = kepler_def(`eta_x x1 x2 x3 =
57 (sqrt ((x1*x2*x3)/(ups_x x1 x2 x3)))
60 let eta_y = kepler_def(`eta_y y1 y2 y3 =
66 let rho_x = kepler_def(`rho_x x1 x2 x3 x4 x5 x6 =
67 --x1*x1*x4*x4 - x2*x2*x5*x5 - x3*x3*x6*x6 +
68 (&2)*x1*x2*x4*x5 + (&2)*x1*x3*x4*x6 + (&2)*x2*x3*x5*x6`);;
70 let rad2_y = kepler_def(`rad2_y y1 y2 y3 y4 y5 y6 =
71 let (x1,x2,x3,x4,x5,x6)= (y1*y1,y2*y2,y3*y3,y4*y4,y5*y5,y6*y6) in
72 (rho_x x1 x2 x3 x4 x5 x6)/((delta_x x1 x2 x3 x4 x5 x6)*(&4))`);;
75 let chi_x = kepler_def(`chi_x x1 x2 x3 x4 x5 x6
76 = -- (x1*x4*x4) + x1*x4*x5 + x2*x4*x5 - x2*x5*x5
77 + x1*x4*x6 + x3*x4*x6 +
78 x2*x5*x6 + x3*x5*x6 - (&2) * x4*x5*x6 - x3*x6*x6`);;
82 (* ------------------------------------------------------------------ *)
83 (* The formula for the dihedral angle of a simplex. *)
84 (* The variables xi are the squares of the lengths of the edges. *)
85 (* The angle is computed along the first edge (x1). *)
86 (* ------------------------------------------------------------------ *)
88 let dih_x = kepler_def(`dih_x x1 x2 x3 x4 x5 x6 =
89 let d_x4 = delta_x4 x1 x2 x3 x4 x5 x6 in
90 let d = delta_x x1 x2 x3 x4 x5 x6 in
91 pi/ (&2) + atn2( (sqrt ((&4) * x1 * d)),-- d_x4)`);;
94 let dih_y = kepler_def(`dih_y y1 y2 y3 y4 y5 y6 =
95 let (x1,x2,x3,x4,x5,x6)= (y1*y1,y2*y2,y3*y3,y4*y4,y5*y5,y6*y6) in
96 dih_x x1 x2 x3 x4 x5 x6`);;
98 let dih2_y = kepler_def(`dih2_y y1 y2 y3 y4 y5 y6 =
99 dih_y y2 y1 y3 y5 y4 y6`);;
101 let dih3_y = kepler_def(`dih3_y y1 y2 y3 y4 y5 y6 =
102 dih_y y3 y1 y2 y6 y4 y5`);;
104 let dih2_x = kepler_def(`dih2_x x1 x2 x3 x4 x5 x6 =
105 dih_x x2 x1 x3 x5 x4 x6`);;
107 let dih3_x = kepler_def(`dih3_x x1 x2 x3 x4 x5 x6 =
108 dih_x x3 x1 x2 x6 x4 x5`);;
111 (* ------------------------------------------------------------------ *)
112 (* Harriot-Euler formula for the area of a spherical triangle *)
113 (* in terms of the angles: area = alpha+beta+gamma - pi *)
114 (* ------------------------------------------------------------------ *)
116 let sol_x = kepler_def(`sol_x x1 x2 x3 x4 x5 x6 =
117 (dih_x x1 x2 x3 x4 x5 x6) +
118 (dih_x x2 x3 x1 x5 x6 x4) + (dih_x x3 x1 x2 x6 x4 x5) - pi`);;
120 let sol_y = kepler_def(`sol_y y1 y2 y3 y4 y5 y6 =
121 (dih_y y1 y2 y3 y4 y5 y6) +
122 (dih_y y2 y3 y1 y5 y6 y4) + (dih_y y3 y1 y2 y6 y4 y5) - pi`);;
125 (* ------------------------------------------------------------------ *)
126 (* The Cayley-Menger formula for the volume of a simplex *)
127 (* The variables xi are the squares of the lengths of the edges. *)
128 (* ------------------------------------------------------------------ *)
130 let vol_x = kepler_def(`vol_x x1 x2 x3 x4 x5 x6 =
131 (sqrt (delta_x x1 x2 x3 x4 x5 x6))/ (&12)`);;
133 (* ------------------------------------------------------------------ *)
134 (* Some lower dimensional funcions and Rogers simplices. *)
135 (* ------------------------------------------------------------------ *)
137 let arclength = kepler_def(`arclength a b c =
138 pi/(&2) + (atn2( (sqrt (ups_x (a*a) (b*b) (c*c))),(c*c - a*a -b*b)))`);;
141 let volR = kepler_def(`volR a b c =
142 (sqrt (a*a*(b*b-a*a)*(c*c-b*b)))/(&6)`);;
144 let solR = kepler_def(`solR a b c =
145 (&2)*atn2( sqrt(((c+b)*(b+a))), sqrt ((c-b)*(b-a)))`);;
147 let dihR = kepler_def(`dihR a b c =
148 atn2( sqrt(b*b-a*a),sqrt (c*c-b*b))`);;
150 let vorR = kepler_def(`vorR a b c =
151 (&4)*(--doct*(volR a b c) + (solR a b c)/(&3))`);;
153 let rad2_x = kepler_def(`rad2_x x1 x2 x3 x4 x5 x6 =
154 (rho_x x1 x2 x3 x4 x5 x6)/((delta_x x1 x2 x3 x4 x5 x6)*(&4))`);;
158 let d3 = new_definition `d3 (v:real^3) w = dist(v,w)`;;
160 let mk_vec3 = new_definition `mk_vec3 a b c = vector[a; b; c]`;;
162 let real3_of_triple = new_definition `real3_of_triple (a,b,c) = (mk_vec3 a b c):real^3`;;
164 let triple_of_real3 = new_definition `triple_of_real3 (v:real^3) =
171 (* aff is deprecated *)
172 let aff = new_definition `aff = ( hull ) affine`;;
174 let lin_combo = new_definition `lin_combo V f = vsum V (\v. f v % (v:real^N))`;;
176 let affsign = new_definition `affsign sgn s t (v:real^A) = (?f.
177 (v = lin_combo (s UNION t) f) /\ (!w. t w ==> sgn (f w)) /\ (sum (s UNION t) f = &1))`;;
180 let sgn_gt = new_definition `sgn_gt = (\t. (&0 < t))`;;
181 let sgn_ge = new_definition `sgn_ge = (\t. (&0 <= t))`;;
182 let sgn_lt = new_definition `sgn_lt = (\t. (t < &0))`;;
183 let sgn_le = new_definition `sgn_le = (\t. (t <= &0))`;;
185 (* conv is deprecated. Use `convex hull S` instead *)
187 let conv = new_definition `conv S:real^A->bool = affsign sgn_ge {} S`;;
188 let conv0 = new_definition `conv0 S:real^A->bool = affsign sgn_gt {} S`;;
189 let cone = new_definition `cone v S:real^A->bool = affsign sgn_ge {v} S`;;
190 let cone0 = new_definition `cone0 v S:real^A->bool = affsign sgn_gt {v} S`;;
193 let aff_gt_def = new_definition `aff_gt = affsign sgn_gt`;;
194 let aff_ge_def = new_definition `aff_ge = affsign sgn_ge`;;
195 let aff_lt_def = new_definition `aff_lt = affsign sgn_lt`;;
196 let aff_le_def = new_definition `aff_le = affsign sgn_le`;;
198 let voronoi = new_definition `voronoi v S = { x | !w. ((S w) /\ ~(w=v)) ==> (dist( x, v) < dist( x, w)) }`;;
200 let voronoi_le = new_definition `voronoi_le v S = { x | !w. ((S w) /\ ~(w=v)) ==> (dist( x, v) <= dist( x, w)) }`;;
202 let line = new_definition `line x = (?v w. ~(v =w) /\ (x = affine hull {v,w}))`;;
204 let plane = new_definition `plane x = (?u v w. ~(collinear {u,v,w}) /\ (x = affine hull {u,v,w}))`;;
205 let closed_half_plane = new_definition `closed_half_plane x = (?u v w. ~(collinear {u,v,w}) /\ (x = aff_ge {u,v} {w}))`;;
206 let open_half_plane = new_definition `open_half_plane x = (?u v w. ~(collinear {u,v,w}) /\ (x = aff_gt {u,v} {w}))`;;
207 let coplanar = new_definition `coplanar S = (?x. plane x /\ S SUBSET x)`;;
208 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'}))`;;
209 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'}))`;;
212 let bis = new_definition `bis u v = {x | dist(x,u) = dist(x,v)}`;;
215 let bis_le = new_definition `bis_le u v = {x | dist(x,u) <= dist(x,v) }`;;
216 let bis_lt = new_definition `bis_lt u v = {x | dist(x,u) < dist(x,v) }`;;
219 let circumcenter = new_definition `circumcenter S = @v. ( (affine hull S) v /\ (?c. !w. (S w) ==> (c = dist(v,w))))`;;
223 let radV = new_definition `radV S = @c. !w. (S w) ==> (c = dist(circumcenter S,w))`;;
228 let orientation = new_definition `orientation S v sgn = affsign sgn (S DIFF {v}) {v} (circumcenter S)`;;
232 let arcV = new_definition `arcV u v w = acs (( (v - u) dot (w - u))/((norm (v-u)) * (norm (w-u))))`;;
234 let dihV = new_definition `dihV w0 w1 w2 w3 =
238 let vap = ( vc dot vc) % va - ( va dot vc) % vc in
239 let vbp = ( vc dot vc) % vb - ( vb dot vc) % vc in
240 arcV (vec 0) vap vbp`;;
242 (* polar coordinates *)
244 let radius = new_definition `radius x y = sqrt((x pow 2) + (y pow 2))`;;
245 let polar_angle = new_definition `polar_angle x y =
246 let theta = atn2(x,y) in
247 if (theta < &0) then (theta + (&2 * pi)) else theta`;;
248 let polar_c = new_definition `polar_c x y = (radius x y,polar_angle x y)`;;
250 let less_polar = new_definition `less_polar (x,y) (x',y') =
251 let (r,theta) = polar_c x y in
252 let (r',theta') = polar_c x' y' in
253 (theta < theta') \/ ((theta =theta') /\ (r < r'))`;;
255 let min_polar = new_definition `min_polar V = ( @ u. V u /\ (!w. V w /\ ~(u = w) ==> (less_polar u w)))`;;
257 let polar_cycle = new_definition `polar_cycle V v =
258 let W = {u | V u /\ less_polar v u} in
259 if (W = EMPTY) then min_polar V else min_polar W`;;
261 let iter_spec = prove(`?iter. !f u:A. (iter 0 f u = u) /\ (!n. (iter (SUC n) f u = f(iter n f u)))`,
262 (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
263 ([REWRITE_TAC[GSYM SKOLEM_THM;num_RECURSION_STD];REWRITE_TAC[]]) THEN
264 (EXISTS_TAC `\ (i:num) (f:A->A) (u:A). (g f u i):A`) THEN
266 (ASM_REWRITE_TAC[]));;
268 let iter = new_specification ["iter"] iter_spec;;
270 let orthonormal = new_definition `orthonormal e1 e2 e3 =
271 (( e1 dot e1 = &1) /\ (e2 dot e2 = &1) /\ ( e3 dot e3 = &1) /\
272 ( e1 dot e2 = &0) /\ ( e1 dot e3 = &0) /\ ( e2 dot e3 = &0) /\
273 (&0 < (cross e1 e2) dot e3))`;;
275 (* spherical coordinates *)
276 let azim_hyp_def = new_definition `azim_hyp = (!v w w1 w2. ?theta. !e1 e2 e3. ?psi h1 h2 r1 r2.
277 ~(collinear {v, w, w1}) /\ ~(collinear {v, w, w2}) /\
278 (orthonormal e1 e2 e3) /\ ((dist( w, v)) % e3 = (w - v)) ==>
279 ((&0 <= theta) /\ (theta < &2 * pi) /\ (&0 < r1) /\ (&0 < r2) /\
280 (w1 - v = (r1 * cos(psi)) % e1 + (r1 * sin(psi)) % e2 + h1 % (w-v)) /\
281 (w2 - v = (r2 * cos(psi + theta)) % e1 + (r2 * sin(psi + theta)) % e2 + h2 % (w-v))))`;;
283 let azim_spec = prove(`?theta. !v w w1 w2 e1 e2 e3. ?psi h1 h2 r1 r2.
285 ~(collinear {v, w, w1}) /\ ~(collinear {v, w, w2}) /\
286 (orthonormal e1 e2 e3) /\ ((dist( w, v)) % e3 = (w - v)) ==>
287 ((&0 <= theta v w w1 w2) /\ (theta v w w1 w2 < &2 * pi) /\ (&0 < r1) /\ (&0 < r2) /\
288 (w1 - v = (r1 * cos(psi)) % e1 + (r1 * sin(psi)) % e2 + h1 % (w-v)) /\
289 (w2 - v = (r2 * cos(psi + theta v w w1 w2)) % e1 + (r2 * sin(psi + theta v w w1 w2)) % e2 + h2 % (w-v)))`,
290 (REWRITE_TAC[GSYM SKOLEM_THM;GSYM RIGHT_IMP_EXISTS_THM;GSYM RIGHT_IMP_FORALL_THM]) THEN
291 (REWRITE_TAC[azim_hyp_def]) THEN
292 (REPEAT STRIP_TAC) THEN
293 (ASM_REWRITE_TAC[RIGHT_IMP_EXISTS_THM]));;
295 let azim_def = new_specification ["azim"] azim_spec;;
298 let cyclic_set = new_definition `cyclic_set W v w =
299 (~(v=w) /\ (FINITE W) /\ (!p q h. W p /\ W q /\ (p = q + h % (v - w)) ==> (p=q)) /\
300 (W INTER (affine hull {v,w}) = EMPTY))`;;
305 let azim_cycle_hyp_def = new_definition `azim_cycle_hyp =
306 (?sigma. !W proj v w e1 e2 e3 p.
308 (cyclic_set W v w) /\ ((dist( v ,w)) % e3 = (w-v)) /\
309 (orthonormal e1 e2 e3) /\
310 (!u x y. (proj u = (x,y)) <=> (?h. (u = v + x % e1 + y % e2 + h % e3))) ==>
311 (proj (sigma W v w p) = polar_cycle (IMAGE proj W) (proj p)))`;;
313 let azim_cycle_spec = prove(`?sigma. !W proj v w e1 e2 e3 p.
314 (azim_cycle_hyp) ==> ( (W p) /\
315 (cyclic_set W v w) /\ ((dist( v ,w)) % e3 = (w-v)) /\
316 (orthonormal e1 e2 e3) /\
317 (!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))`,
318 (REWRITE_TAC[GSYM RIGHT_IMP_EXISTS_THM;GSYM RIGHT_IMP_FORALL_THM]) THEN
319 (REWRITE_TAC[azim_cycle_hyp_def])
322 let azim_cycle_def = new_specification ["azim_cycle"] azim_cycle_spec;;
325 (* ------------------------------------------------------------------ *)
326 (* Definitions from the Collection in Elementary Geometry *)
327 (* ------------------------------------------------------------------ *)
331 let cayleyR = new_definition `cayleyR x12 x13 x14 x15 x23 x24 x25 x34 x35 x45 =
332 -- (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 +
333 &2 *x15*x15*x23*x24 - x13*x13*x24*x24 + &2 *x13*x15*x24*x24 - x15*x15*x24*x24 - &2 *x13*x14*x23*x25 +
334 &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 +
335 &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 -
336 &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 +
337 &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 -
338 &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 -
339 x12*x12*x34*x34 + &2 *x12*x15*x34*x34 - x15*x15*x34*x34 + &2 *x12*x25*x34*x34 + &2 *x15*x25*x34*x34 -
340 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 +
341 &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 +
342 &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 +
343 &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 +
344 &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 -
345 x12*x12*x35*x35 + &2 *x12*x14*x35*x35 - x14*x14*x35*x35 + &2 *x12*x24*x35*x35 + &2 *x14*x24*x35*x35 -
346 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 +
347 &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 -
348 &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 -
349 &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 -
350 &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 +
351 &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 -
352 &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 -
353 x13*x13*x45*x45 + &2 *x12*x23*x45*x45 + &2 *x13*x23*x45*x45 - x23*x23*x45*x45`;;
358 let packing = new_definition `packing S = (!u v. S u /\ S v /\ ~(u = v) ==> (&2 <= dist( u, v)))`;;
362 let wedge = new_definition (`wedge v1 v2 w1 w2 =
366 let n = cross z u1 in
368 if (aff_ge {v1,v2} {w1} w2) then {} else
369 if (aff_lt {v1,v2} {w1} w2) then aff_gt {v1,v2,w1} {n} else
370 if (d > &0) then aff_gt {v1,v2} {w1,w2} else
371 (:real^3) DIFF aff_ge {v1,v2} {w1,w2}`);;
373 (* ------------------------------------------------------------------ *)
374 (* Measure and Volume, following Nguyen Tat Thang *)
375 (* ------------------------------------------------------------------ *)
377 let sphere= new_definition`sphere x=(?(v:real^3)(r:real). (r> &0)/\ (x={w:real^3 | norm (w-v)= r}))`;;
379 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)}`;;
381 let circular_cone =new_definition `circular_cone (V:real^3-> bool)=
382 (? (v,w:real^3)(r:real). V= c_cone (v,w,r))`;;
384 let NULLSET_RULES,NULLSET_INDUCT,NULLSET_CASES =
385 new_inductive_definition
386 `(!P. ((plane P)\/ (sphere P) \/ (circular_cone P)) ==> NULLSET P) /\
387 !(s:real^3->bool) t. (NULLSET s /\ NULLSET t) ==> NULLSET (s UNION t)`;;
389 let null_equiv = new_definition `null_equiv (s,t :real^3->bool)=(? (B:real^3-> bool). NULLSET B /\
390 ((s DIFF t) UNION (t DIFF s)) SUBSET B)`;;
393 let normball = new_definition `normball x r = { y:real^A | dist(y,x) < r}`;;
395 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))`;;
397 let eventually_radial = new_definition `eventually_radial x C <=> (?r. (r> &0) /\ radial r x (C INTER normball x r))`;;
399 let solid_triangle = new_definition `solid_triangle v0 S r = normball v0 r INTER cone v0 S`;;
401 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)}`;;
405 let rcone_ge = new_definition `rcone_ge = rconesgn ( >= )`;;
406 let rcone_gt = new_definition `rcone_gt = rconesgn ( > )`;;
407 let rcone_lt = new_definition `rcone_lt = rconesgn ( < )`;;
408 let rcone_eq = new_definition `rcone_eq = rconesgn ( = )`;;
410 let scale = new_definition `scale (t:real^3) (u:real^3) = vector[t$1 * u$1; t$2 * u$2; t$3 * u$3]`;;
412 let ellipsoid = new_definition `ellipsoid t r = IMAGE (scale t) (normball(vec 0)r)`;;
414 let conic_cap = new_definition `conic_cap v0 v1 r a = normball v0 r INTER rcone_gt v0 v1 a`;;
416 let frustum = new_definition `frustum v0 v1 h1 h2 a = { y | rcone_gt v0 v1 a y /\
417 let d = (y - v0) dot (v1 - v0) in
418 let n = norm(v1 - v0) in
419 (h1*n < d /\ d < h2*n)}`;;
421 let frustt = new_definition `frustt v0 v1 h a = frustum v0 v1 (&0) h a`;;
423 let rect = new_definition `rect (a:real^3) (b:real^3) = {(v:real^3) | !i. ( a$i < v$i /\ v$i < b$i )}`;;
426 let is_tetrahedron = new_definition `is_tetrahedron S = ?v0 v1 v2 v3. (S = conv0 {v0,v1,v2,v3})`;;
429 let primitive = new_definition `primitive (C:real^3->bool) =
430 ((?v0 v1 v2 v3 r. (C = solid_triangle v0 {v1,v2,v3} r)) \/
431 (?v0 v1 v2 v3. (C = conv0 {v0,v1,v2,v3})) \/
432 (?v0 v1 v2 v3 h a. (C = frustt v0 v1 h a INTER wedge v0 v1 v2 v3)) \/
433 (?v0 v1 v2 v3 r c. (C = conic_cap v0 v1 r c INTER wedge v0 v1 v2 v3)) \/
434 (?a b. (C = rect a b)) \/
435 (?t r. (C = ellipsoid t r)) \/
436 (?v0 v1 v2 v3 r. (C = normball v0 r INTER wedge v0 v1 v2 v3)))`;;
438 let MEASURABLE_RULES,MEASURABLE_INDUCT,MEASURABLE_CASES =
439 new_inductive_definition
440 `(!C. primitive C ==> measurable C) /\
441 ( !Z. NULLSET Z ==> measurable Z) /\
442 (!X t. measurable X ==> (measurable (IMAGE (scale t) X))) /\
443 (!X v. measurable X ==> (measurable (IMAGE ((+) v) X))) /\
444 ( !(s:real^3->bool) t. (measurable s /\ measurable t) ==> measurable (s UNION t)) /\
445 ( !(s:real^3->bool) t. (measurable s /\ measurable t) ==> measurable (s INTER t)) /\
446 ( !(s:real^3->bool) t. (measurable s /\ measurable t) ==> measurable (s DIFF t))
450 let SDIFF = new_definition `SDIFF X Y = (X DIFF Y) UNION (Y DIFF X)`;;
453 let vol_solid_triangle = new_definition `vol_solid_triangle v0 v1 v2 v3 r =
454 let a123 = dihV v0 v1 v2 v3 in
455 let a231 = dihV v0 v2 v3 v1 in
456 let a312 = dihV v0 v3 v1 v2 in
457 (a123 + a231 + a312 - pi)*(r pow 3)/(&3)`;;
459 let vol_frustt_wedge = new_definition `vol_frustt_wedge v0 v1 v2 v3 h a =
460 (azim v0 v1 v2 v3)*(h pow 3)*(&1/(a*a) - &1)/(&6)`;;
462 (* volume of intersection of conic cap and wedge *)
463 let vol_conic_cap_wedge = new_definition `vol_conic_cap_wedge v0 v1 v2 v3 r c =
464 (azim v0 v1 v2 v3)*(&1 - c)*(r pow 3)/(&3)`;;
467 let vol_conv = new_definition `vol_conv v1 v2 v3 v4 =
468 let x12 = dist(v1,v2) pow 2 in
469 let x13 = dist(v1,v3) pow 2 in
470 let x14 = dist(v1,v4) pow 2 in
471 let x23 = dist(v2,v3) pow 2 in
472 let x24 = dist(v2,v4) pow 2 in
473 let x34 = dist(v3,v4) pow 2 in
474 sqrt(delta_x x12 x13 x14 x34 x24 x34)/(&12)`;;
476 let vol_rect = new_definition `vol_rect a b =
477 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`;;
479 let vol_ball_wedge = new_definition `vol_ball_wedge v0 v1 v2 v3 r =
480 (azim v0 v1 v2 v3)*(&2)*(r pow 3)/(&3)`;;
483 let volume_props = new_definition `volume_props (vol:(real^3->bool)->real) =
484 ( (!C. vol C >= &0) /\
485 (!Z. NULLSET Z ==> (vol Z = &0)) /\
486 (!X Y. measurable X /\ measurable Y /\ NULLSET (SDIFF X Y) ==> (vol X = vol Y)) /\
487 (!X t. (measurable X) /\ (measurable (IMAGE (scale t) X)) ==> (vol (IMAGE (scale t) X) = abs(t$1 * t$2 * t$3)*vol(X))) /\
488 (!X v. measurable X ==> (vol (IMAGE ((+) v) X) = vol X)) /\
489 (!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) /\
490 (!v0 v1 v2 v3. vol(conv0 {v0,v1,v2,v3}) = vol_conv v0 v1 v2 v3) /\
491 (!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) /\
492 (!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)) /\
493 (!(a:real^3) (b:real^3). vol(rect a b) = vol_rect a b) /\
494 (!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)))`;;
496 let vol_def = new_definition `vol = ( @ ) volume_props`;;
498 let solid = new_definition `solid x C = (@s. ?c. (c > &0) /\ (!r. (r > &0) /\ (r < c) ==>
499 (s = &3 * vol(C INTER normball x r)/(r pow 3)))) `;;
501 let sovo = new_definition `sovo x C (v,s) = v* vol(C) + s* solid x C`;;
503 let phivo = new_definition `phivo h t (v,s) = v*t*h*(t+h)/(&6) + s`;;
505 let avo = new_definition `avo h t l= (&1 - h/t)*(phivo h t l - phivo t t l)`;;
507 let ortho0 = new_definition `ortho0 x v1 v2 v3 = conv0 {x,x+v1,x+v1+v2,x+v1+v2+v3}`;;
509 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))`;;
511 let rogers = new_definition `rogers v0 v1 v2 v3 c =
512 let w = (&1/ (&2)) % (v0 + v1) in
513 let p = circumcenter {v0,v1,v2} in
514 let q = make_point v0 v1 v2 v3 c c c in
517 let rogers0 = new_definition `rogers0 v0 v1 v2 v3 c =
518 let w = (&1/ (&2)) % (v0 + v1) in
519 let p = circumcenter {v0,v1,v2} in
520 let q = make_point v0 v1 v2 v3 c c c in
523 let abc_param = new_definition `abc_param v0 v1 v2 c =
524 let a = (&1/(&2)) * dist(v0,v1) in
525 let b = radV {v0,v1,v2} in