Update from HH
[Flyspeck/.git] / formal_lp / old / formal_interval / test.hl
1 needs "../formal_lp/formal_interval/verifier.hl";;
2
3
4 (************)
5 (* Examples *)
6 (************)
7
8 (* f(x) = x * x - atn x; f(x) < 0 when 0.01 <= x <= 0.8 *)
9 prove_ineq1d 5 `x * x - atn x` `#0.01` `#0.8`;;
10
11 (* f(x) = x * x * x - atn (x * x); f(x) < 0 when 0.01 <= x <= 0.8 *)
12 prove_ineq1d 5 `x * x * x - atn (x * x)` `#0.01` `#0.8`;;
13
14 (* f(x) = atn x - x; f(x) < 0 when 0.001 <= x <= sqrt(2) *)
15 prove_ineq1d 10 `atn x - x` `#0.001` `sqrt (&2) + pi`;;
16
17
18 (****************************************)
19 (* arclength #2.52 #2.52 (&2) >= #0.816 *)
20 (* arclength (&2) #2.52 (&2) - arclength #2.52 #2.52 (&2) >= #0.073 *)
21 (* arclength a b c = acs ((a * a + b * b - c * c) / (&2 * a * b)) *)
22
23 let arclength_tm = `acs ((a * a + b * b - c * c) / (&2 * a * b))`;;
24 let tm1 = subst[`#2.52`, `a:real`; `#2.52`, `b:real`; `&2`, `c:real`] arclength_tm and
25     tm2 = subst[`&2`, `a:real`; `#2.52`, `b:real`; `&2`, `c:real`] arclength_tm;;
26
27 let ineq1, ineq2 =
28   let (-) = mk_binop sub_op_real in
29     `#0.816` - tm1, `#0.073` - (tm2 - tm1);;
30
31 prove_ineq1d 6 ineq1 `&0` `&0`;;
32 prove_ineq1d 6 ineq2 `&0` `&0`;;
33
34
35 (********************************************)
36 (* !x. &2 <= x /\ x <= #2.52                                           *)
37 (*   ==> arclength (&2) x (&2) - arclength (#2.52) x (&2) >= #0.073`   *)
38
39 let ineq =
40   let tm1 = subst[`&2`, `a:real`; `x:real`, `b:real`; `&2`, `c:real`] arclength_tm and
41       tm2 = subst[`#2.52`, `a:real`; `x:real`, `b:real`; `&2`, `c:real`] arclength_tm and
42       (-) = mk_binop sub_op_real in
43     `#0.073` - (tm1 - tm2);;
44
45 prove_ineq1d 7 ineq `&2` `#2.52`;;
46
47
48 (***********************)
49 (* !x. &2 <= x /\ x <= #2.52  ==> &1 / &4 * --inv (sqrt (&1 - (x / &4) pow 2)) <= -- #0.28 *)
50
51 prove_ineq1d 4 `#0.28 + &1 / &4 * --inv (sqrt (&1 - (x * inv(&4)) * (x * inv(&4))))` `&2` `#2.52`;;
52
53 (* !x. &2 <= x /\ x <= #2.52 ==>
54         inv (sqrt (&1 - ((&3969 / &625 + x * x - &4) / (&2 * #2.52 * x)) pow 2)) <= &2 *)
55
56 let ineq = (rand o concl o REWRITE_CONV[REAL_POW_2])
57   `inv (sqrt (&1 - ((&3969 / &625 + x * x - &4) / (&2 * #2.52 * x)) pow 2)) - &2`;;
58
59 prove_ineq1d 3 ineq `&2` `#2.52`;;
60
61
62
63 (***********************************************)
64 (* |atn x - x / (1 + 0.28 * x * x)| < 0.005 when x IN [0, 1] *)
65 let abs_lemma = prove(`f - eps < &0 ==> -- eps - f < &0 ==> abs f < eps`, REAL_ARITH_TAC);;
66
67 let tm = `atn x - x / (&1 + #0.28 * x * x)`;;
68 let eps = `#0.005`;;
69 let hi_tm = mk_binop sub_op_real tm eps;;
70 let lo_tm = mk_binop sub_op_real (mk_comb (neg_op_real, eps)) tm;;
71 let hi_th = (UNDISCH_ALL o SPEC_ALL) (prove_ineq1d 5 hi_tm `&0` `&1`);;
72 let lo_th = (UNDISCH_ALL o SPEC_ALL) (prove_ineq1d 6 lo_tm `&0` `&1`);;
73 (GEN_ALL o DISCH_ALL) (MATCH_MP (MATCH_MP abs_lemma hi_th) lo_th);;
74
75
76
77 (*************************************)
78 (* A polynomial approximation of atn *)
79 let tm = `x * (&1 - (x * x) * (&11184811 / &33554432 - (x * x) * (&13421773 / &67108864))) - atn x`;;
80 let eps = `#0.0000001`;;
81 let hi_tm = mk_binop sub_op_real tm eps;;
82 let lo_tm = mk_binop sub_op_real (mk_comb (neg_op_real, eps)) tm;;
83 let hi_th = (UNDISCH_ALL o SPEC_ALL) (prove_ineq1d 8 hi_tm `&0` `&1 / &30`) and
84     lo_th = (UNDISCH_ALL o SPEC_ALL) (prove_ineq1d 8 lo_tm `&0` `&1 / &30`) and
85     hi_th' = (UNDISCH_ALL o SPEC_ALL) (prove_ineq1d 8 hi_tm `-- &1 / &30` `&0`) and
86     lo_th' = (UNDISCH_ALL o SPEC_ALL) (prove_ineq1d 8 lo_tm `-- &1 / &30` `&0`);;
87
88
89 (GEN_ALL o DISCH_ALL) (MATCH_MP (MATCH_MP abs_lemma hi_th) lo_th);;
90
91
92
93 (***********************************)
94 (* A rational approximation of atn *)
95
96
97 let atan1 = `let P00, P01, P02 = #4.26665376811246382, #3.291955407318148, #0.281939359003812 in
98   let Q00, Q01 = #4.26665376811354027, #4.714173328637605 in
99   let x1 = x * x in
100   let x2 = x1 * x1 in
101     x * (x1 * P01 + x2 * P02 + P00) / (x1 * Q01 + x2 + Q00) - atn x`;;
102 let atan1_tm = (rand o concl o REPEATC let_CONV) atan1;;
103
104 let eps = `#0.00001`;;
105 let hi_tm = mk_binop sub_op_real atan1_tm eps;;
106 let lo_tm = mk_binop sub_op_real (mk_comb (neg_op_real, eps)) atan1_tm;;
107 let hi_th = (UNDISCH_ALL o SPEC_ALL) (prove_ineq1d 10 hi_tm `&0` `#0.2`);;
108 let lo_th = (UNDISCH_ALL o SPEC_ALL) (prove_ineq1d 10 lo_tm `&0` `#0.2`);;
109 (GEN_ALL o DISCH_ALL) (MATCH_MP (MATCH_MP abs_lemma hi_th) lo_th);;
110
111
112 (*
113 (* A higher precision approximation *)
114 let eps = `#0.0000000001`;;
115 let hi_tm = mk_binop sub_op_real atan1_tm eps;;
116 let lo_tm = mk_binop sub_op_real (mk_comb (neg_op_real, eps)) atan1_tm;;
117 let hi_th = (UNDISCH_ALL o SPEC_ALL) (prove_ineq1d 15 hi_tm `&0` `#0.2`);;
118 let lo_th = (UNDISCH_ALL o SPEC_ALL) (prove_ineq1d 15 lo_tm `&0` `#0.2`);;
119 (GEN_ALL o DISCH_ALL) (MATCH_MP (MATCH_MP abs_lemma hi_th) lo_th);;
120 *)
121
122
123
124 (****************************************************************)
125 (* Test 1: f(x) = x * x - atn x; f(x) < 0 when 0.01 <= x <= 0.8 *)
126 (****************************************************************)
127
128 llet eps = `#0.00001`;;
129 let hi_tm = mk_binop sub_op_real atan1_tm eps;;
130 let lo_tm = mk_binop sub_op_real (mk_comb (neg_op_real, eps)) atan1_tm;;
131 let hi_th = (UNDISCH_ALL o SPEC_ALL) (prove_ineq1d 10 hi_tm `&0` `#0.2`);;
132 let lo_th = (UNDISCH_ALL o SPEC_ALL) (prove_ineq1d 10 lo_tm `&0` `#0.2`);;
133 (GEN_ALL o DISCH_ALL) (MATCH_MP (MATCH_MP abs_lemma hi_th) lo_th);;
134 let eps = `#0.00001`;;
135 let hi_tm = mk_binop sub_op_real atan1_tm eps;;
136 let lo_tm = mk_binop sub_op_real (mk_comb (neg_op_real, eps)) atan1_tm;;
137 let hi_th = (UNDISCH_ALL o SPEC_ALL) (prove_ineq1d 10 hi_tm `&0` `#0.2`);;
138 let lo_th = (UNDISCH_ALL o SPEC_ALL) (prove_ineq1d 10 lo_tm `&0` `#0.2`);;
139 (GEN_ALL o DISCH_ALL) (MATCH_MP (MATCH_MP abs_lemma hi_th) lo_th);;
140 let eps = `#0.00001`;;
141 let hi_tm = mk_binop sub_op_real atan1_tm eps;;
142 let lo_tm = mk_binop sub_op_real (mk_comb (neg_op_real, eps)) atan1_tm;;
143 let hi_th = (UNDISCH_ALL o SPEC_ALL) (prove_ineq1d 10 hi_tm `&0` `#0.2`);;
144 let lo_th = (UNDISCH_ALL o SPEC_ALL) (prove_ineq1d 10 lo_tm `&0` `#0.2`);;
145 (GEN_ALL o DISCH_ALL) (MATCH_MP (MATCH_MP abs_lemma hi_th) lo_th);;
146 et x_tm = mk_float 1 48;;
147 let z_tm = mk_float 8 49;;
148
149 let domain_th = mk_center_domain 5 x_tm z_tm;;
150 let _, y_tm, _, _ = dest_cell_domain (concl domain_th);;
151
152
153 let domain1_th = mk_center_domain 5 x_tm y_tm;;
154 let _, y1_tm, _, _ = dest_cell_domain (concl domain1_th);;
155 let domain11_th = mk_center_domain 5 x_tm y1_tm;;
156 let domain12_th = mk_center_domain 5 y1_tm y_tm;;
157
158 let domain2_th = mk_center_domain 5 y_tm z_tm;;
159
160
161 let eval_f pp domain_th =
162   let int_x = eval_taylor_x domain_th in
163   let int_atn = eval_taylor_atn pp domain_th in
164   let int_x2 = taylor_interval_mul pp int_x int_x in
165     taylor_interval_sub pp int_x2 int_atn;;
166
167
168 let int_th = eval_f 5 domain_th;;
169 let int_th11 = eval_f 5 domain11_th;;
170 let int_th12 = eval_f 5 domain12_th;;
171 let int_th2 = eval_f 5 domain2_th;;
172
173 taylor_cell_pass 5 int_th;;
174 let pass11 = taylor_cell_pass 5 int_th11;;
175 let pass12 = taylor_cell_pass 5 int_th12;;
176 let pass2 = taylor_cell_pass 5 int_th2;;
177
178 let pass1 = cell_pass_glue pass11 pass12;;
179 let pass = cell_pass_glue pass1 pass2;;
180
181
182
183 (*************************************************************************)
184 (* Test 2: f(x) = x * x * x - atn (x * x); f(x) < 0 when 0.1 <= x <= 0.8 *)
185 (*************************************************************************)
186
187 let x_tm = mk_float 1 49;;
188 let z_tm = mk_float 8 49;;
189
190 let domain_th = mk_center_domain 5 x_tm z_tm;;
191 let _, y_tm, _, _ = dest_cell_domain (concl domain_th);;
192
193 let domain1_th = mk_center_domain 5 x_tm y_tm;;
194 let _, y1_tm, _, _ = dest_cell_domain (concl domain1_th);;
195
196 let domain11_th = mk_center_domain 5 x_tm y1_tm;;
197 let _, y11_tm, _, _ = dest_cell_domain (concl domain11_th);;
198
199 let domain111_th = mk_center_domain 5 x_tm y11_tm;;
200 let domain112_th = mk_center_domain 5 y11_tm y1_tm;;
201
202 let domain12_th = mk_center_domain 5 y1_tm y_tm;;
203 let domain2_th = mk_center_domain 5 y_tm z_tm;;
204
205 let eval_sqr pp domain_th =
206   let x_th = eval_taylor_x domain_th in
207     taylor_interval_mul pp x_th x_th;;
208
209
210 let eval_f pp domain_th =
211   let int_x = eval_taylor_x domain_th in
212   let int_x2 = taylor_interval_mul pp int_x int_x in
213   let int_x3 = taylor_interval_mul pp int_x int_x2 in
214   let atn_x2 = taylor_interval_compose pp eval_taylor_atn eval_sqr domain_th in
215     taylor_interval_sub pp int_x3 atn_x2;;
216
217 let int_th = eval_f 5 domain_th;;
218 let int_th111 = eval_f 5 domain111_th;;
219 let int_th112 = eval_f 5 domain112_th;;
220 let int_th12 = eval_f 5 domain12_th;;
221 let int_th2 = eval_f 5 domain2_th;;
222
223 taylor_cell_pass 5 int_th;;
224 let pass111 = taylor_cell_pass 5 int_th111;;
225 let pass112 = taylor_cell_pass 5 int_th112;;
226 let pass12 = taylor_cell_pass 5 int_th12;;
227 let pass2 = taylor_cell_pass 5 int_th2;;
228
229 let pass11 = cell_pass_glue pass111 pass112;;
230 let pass1 = cell_pass_glue pass11 pass12;;
231 let pass = cell_pass_glue pass1 pass2;;
232
233 let th0 = REWRITE_RULE[cell_pass; interval_arith] pass;;
234 CONV_RULE (DEPTH_CONV FLOAT_TO_NUM_CONV) th0;;
235
236
237
238 (********************************************)
239 (* Test 3:                                                             *)
240 (* !x. &2 <= x /\ x <= #2.52                                           *)
241 (*   ==> arclength (&2) x (&2) - arclength (#2.52) x (&2) >= #0.073`   *)
242
243 let eval_arc1 pp domain_th =
244   let two = eval_taylor_const domain_th two_float and
245       ( * ) = taylor_interval_mul pp and
246       (+) = taylor_interval_add pp and
247       (-) = taylor_interval_sub pp in
248   let b = (fun pp domain_th ->
249              let x = eval_taylor_x domain_th in
250                two * two * x) in
251   let a = (fun pp domain_th -> 
252              let x = eval_taylor_x domain_th in
253              let rhs = taylor_interval_compose pp eval_taylor_inv b domain_th in
254                (two * two + x * x - two * two) * rhs) in
255     taylor_interval_compose pp eval_taylor_acs a domain_th;;
256
257
258 let h_float = mk_float 252 48;;
259 let eval_arc2 pp domain_th =
260   let two = eval_taylor_const domain_th two_float and
261       h = eval_taylor_const domain_th h_float and
262       ( * ) = taylor_interval_mul pp and
263       (+) = taylor_interval_add pp and
264       (-) = taylor_interval_sub pp in
265   let b = (fun pp domain_th ->
266              let x = eval_taylor_x domain_th in
267                two * h * x) in
268   let a = (fun pp domain_th -> 
269              let x = eval_taylor_x domain_th in
270              let rhs = taylor_interval_compose pp eval_taylor_inv b domain_th in
271                (h * h + x * x - two * two) * rhs) in
272     taylor_interval_compose pp eval_taylor_acs a domain_th;;
273
274
275 let t_float = mk_float 73 47;;
276 let eval_f pp domain_th =
277   let t = eval_taylor_const domain_th t_float and
278       (+) = taylor_interval_add pp and
279       (-) = taylor_interval_sub pp and
280       arc1 = eval_arc1 pp domain_th and
281       arc2 = eval_arc2 pp domain_th in
282     t + arc2 - arc1;;
283
284 let pp = 6;;  
285 let domain0 = mk_center_domain pp two_float h_float;;
286
287 let _, y1_tm, _, _ = dest_cell_domain (concl domain0);;
288 let domain1 = mk_center_domain pp two_float y1_tm and
289     domain2 = mk_center_domain pp y1_tm h_float;;
290
291 let _, y12_tm, _, _ = dest_cell_domain (concl domain2);;
292 let domain21 = mk_center_domain pp y1_tm y12_tm and
293     domain22 = mk_center_domain pp y12_tm h_float;;
294
295 let _, y122_tm, _, _ = dest_cell_domain (concl domain22);;
296 let domain221 = mk_center_domain pp y12_tm y122_tm and
297     domain222 = mk_center_domain pp y122_tm h_float;;
298
299 let int1 = eval_f pp domain1;;
300 let int21 = eval_f pp domain21;;
301 let int221 = eval_f pp domain221;;
302 let int222 = eval_f pp domain222;;
303
304 let pass1 = taylor_cell_pass pp int1;;
305 let pass21 = taylor_cell_pass pp int21;;
306
307 let pass221 = taylor_cell_pass pp int221;;
308 let pass222 = taylor_cell_pass pp int222;;
309
310 let pass22 = cell_pass_glue pass221 pass222;;
311 let pass2 = cell_pass_glue pass21 pass22;;
312
313 let pass = cell_pass_glue pass1 pass2;;
314
315 let th0 = REWRITE_RULE[cell_pass; interval_arith] pass;;
316 CONV_RULE (DEPTH_CONV FLOAT_TO_NUM_CONV) th0;;
317
318
319 let ineq_test pp =
320   let domain0 = mk_center_domain pp two_float h_float in
321   let _, y1_tm, _, _ = dest_cell_domain (concl domain0) in
322   let domain1 = mk_center_domain pp two_float y1_tm and
323       domain2 = mk_center_domain pp y1_tm h_float in
324   let _, y12_tm, _, _ = dest_cell_domain (concl domain2) in
325   let domain21 = mk_center_domain pp y1_tm y12_tm and
326       domain22 = mk_center_domain pp y12_tm h_float in
327   let _, y122_tm, _, _ = dest_cell_domain (concl domain22) in
328   let domain221 = mk_center_domain pp y12_tm y122_tm and
329       domain222 = mk_center_domain pp y122_tm h_float in
330
331   let int1 = eval_f pp domain1 and
332       int21 = eval_f pp domain21 and
333       int221 = eval_f pp domain221 and
334       int222 = eval_f pp domain222 in
335
336   let pass1 = taylor_cell_pass pp int1 in
337   let pass21 = taylor_cell_pass pp int21 in
338   let pass221 = taylor_cell_pass pp int221 in
339   let pass222 = taylor_cell_pass pp int222 in
340   let pass22 = cell_pass_glue pass221 pass222 in
341   let pass2 = cell_pass_glue pass21 pass22 in
342     cell_pass_glue pass1 pass2;;
343
344
345 ineq_test 5;;
346 ineq_test 6;;
347 ineq_test 7;;
348 ineq_test 10;;
349
350 (* 10: 9.768 *)
351 test 5 ineq_test 6;;
352 (* 10: 11.012 *)
353 test 5 ineq_test 7;;
354 (* 10: 16.381 *)
355 test 5 ineq_test 10;;
356
357
358
359
360 (***************************)
361
362 let x_tm = mk_float 20 48;;
363 let z_tm = mk_float 99 48;;
364
365 let domain_th = mk_center_domain 5 x_tm z_tm;;
366 test 1000 (mk_center_domain 5 x_tm) z_tm;;
367
368 let int_x = eval_taylor_x domain_th;;
369 let int_atn = eval_taylor_atn 5 domain_th;;
370 let int_x2 = taylor_interval_mul 5 int_x int_x;;
371 let int_x_atn = taylor_interval_mul 5 int_x int_atn;;
372
373 let int_th1 = taylor_interval_sub 5 int_x2 int_x;;
374 let int_th2 = taylor_interval_sub 5 int_x_atn int_x;;
375
376 eval_taylor_f_bounds 5 int_th1;;
377 eval_taylor_f_bounds 5 int_th2;;
378
379 let x_th =
380   let x_tm, y_tm, z_tm, w_tm = dest_cell_domain (concl domain_th) in
381     (mk_bounded_on_int o ASSUME) (mk_interval x_var_real (mk_pair (x_tm, z_tm)));;
382
383
384 let b_x2 = bounded_on_int_mul 5 x_th x_th;;
385 bounded_on_int_sub 5 b_x2 x_th;;
386
387
388
389
390
391
392 (*******************************************)
393 (* arclength #2.52 #2.52 (&2) >= #0.816 *)
394 (* arclength (&2) #2.52 (&2) - arclength #2.52 #2.52 (&2) >= #0.073 *)
395 (* arclength a b c = acs ((a * a + b * b - c * c) / (&2 * a * b)) *)
396
397 let two_tm = mk_float 2 50;;
398 let h_tm = mk_float 252 48;;
399 let pp = 6;;
400
401
402 let arclength_int pp a b c = 
403   let ( * ) = float_interval_mul pp and
404       (+) = float_interval_add pp and
405       (/) = float_interval_div pp and
406       (-) = float_interval_sub pp and
407       acs = float_interval_acs pp in
408     acs ((a * a + b * b - c * c) / (two_interval * a * b));;
409
410 let result =
411   let (!) = mk_const_interval and
412       (-) = float_interval_sub pp in
413     arclength_int pp !h_tm !h_tm !two_tm,
414   arclength_int pp !two_tm !h_tm !two_tm - arclength_int pp !h_tm !h_tm !two_tm;;
415   
416
417 (CONV_RULE (DEPTH_CONV FLOAT_TO_NUM_CONV)) (fst result);;
418 (CONV_RULE (DEPTH_CONV FLOAT_TO_NUM_CONV)) (snd result);;