1 (* ========================================================================== *)
2 (* FLYSPECK - BOOK FORMALIZATION *)
3 (* Section: Local Fan Main Estimate *)
4 (* Chapter: Local Fan *)
5 (* Author: Thomas C. Hales *)
7 (* Date: 2013-03, revised *)
8 (* ========================================================================== *)
11 Check completeness (informally)
12 of case anaysis in proof of Kepler conjecture, Main Estimate.
14 The purpose is to transform the list init_cs into terminal_cs
15 in a finite sequence of regulated moves.
17 This proves the arrow S_init ==> S_term described in the
18 appendix to the Local Fan chapter.
23 module Check_completeness = struct
27 We don't need the full flyspeck build, really just lib.ml.
29 #directory "/Users/thomashales/Desktop/googlecode/hol_light";;
33 (try Sys.getenv "FLYSPECK_DIR" with Not_found -> Sys.getcwd());;
35 loadt (flyspeck_dir^"/../glpk/sphere.ml");;
36 loadt (flyspeck_dir^"/general/lib.hl");;
42 ****************************************
44 ****************************************
50 let false_on_fail b x = try (b x) with _ -> false;;
52 let fail_on_false b x s = (b x or failwith "s");;
55 mapfilter (function |Some x -> x |None -> failwith "none" ) xs;;
57 let length = List.length;;
61 let sortuniq = setify;;
64 let sortuniq xs = uniq (Lib.sort (<) xs);;
67 let psort u = if fst u <= snd u then u else (snd u,fst u);;
72 | a::rest -> (map (fun x -> (a,x)) b) @ cart rest b;;
74 let cart2 a = cart a a;;
79 let unsplit d f = function
80 | (x::xs) -> List.fold_left (fun s t -> s^d^(f t)) (f x) xs
83 (* let join_comma = unsplit "," (fun x-> x);;*)
84 let join_semi = unsplit ";" (fun x-> x);;
85 let join_lines = unsplit "\n" (fun x-> x);;
86 (* let join_space = unsplit " " (fun x-> x);; *)
88 let string_of_list f xs =
90 "["^(join_semi ss)^"]";;
92 let string_of_pair (i,j) =
93 Printf.sprintf "(%d,%d)" i j;;
95 let string_of_triple (i,j,s) =
96 Printf.sprintf "(%d,%d,%3.2f)" i j s;;
99 let ks = filter (fun (i,j) -> i<j) (cart2 (0--(k-1))) in
100 let g (i,j) = (i,j,f i j) in
101 string_of_list string_of_triple (map g ks);;
103 let arc = Sphere_math.arclength;; (* in glpk directory *)
104 let sqrt = Pervasives.sqrt;;
105 let cos = Pervasives.cos;;
107 (* a few constants. We do tests of exact equality on floats,
108 but this should work because we never do arithmetic on floats
109 after initialization (except in one spot where a machine_eps
110 is thrown in), so there is no source of roundoff error
111 after initialization. *)
115 let sqrt8 = Pervasives.sqrt(8.0);;
117 (* 0.0 insert on djz on 2013-4-20 because when long edge is cstab get 0.11 + 0.1 (cstab-cstab) = 0.11. *)
118 let djz = 0.11 +. 0.0 *. 0.1 *.(cstab -. sqrt8);; (* ear-value on cstab edge *)
122 let sqrt1553 = sqrt(15.53);;
123 let arc1553 = arc two two sqrt1553;;
126 let fourh0 = 2.0 *. twoh0;;
127 let upperbd = 6.0;; (* 6.0 > 4 * h0, upper bound on diags in BB *)
130 if (i <= 3) then zero
131 else if (i=4) then 0.206
132 else if (i=5) then 0.4819
133 else if (i=6) then 0.712
136 let tame_table_d r s =
137 if (r + 2 * s <= 3) then 0.0
139 let r' = float_of_int r in
140 let s' = float_of_int s in
141 0.103 *. (2.0 -. s') +. 0.2759 *. (r' +. 2.0 *. s' -. 4.0);;
145 look only at global minumum points for tau*(s,v) (for s fixed)
146 and such that tau*(s,v)<=0.
147 Let index(v) = the number of edges st |vi-vj|=a(i,j).
148 and such that index(v) is minimized among global minimum points in Bs.
149 Call this set MM(s) (index-minimizers).
151 The semantics that guide us will be counterexample propagation;
152 if a stable constraint system contains a point, does the
153 output of the function contain a nonempty stable constraint system?
155 Formally, we will prove statements of the form for various
157 |- nonempty_M cs ==> (?cs'. mem cs'(f cs) /\ nonempty_M cs').
159 Termination of algorithms:
160 Every slice decreases k.
161 Every subdivision restricts the B-field through a finite set of choices
162 (retaining the M-field)
163 Every deformation sets a new M-field constraint.
165 As a matter of programming style.
166 Procedures that refer to particular cs's should be avoided,
167 except the list terminal list. Pass everything else in by argument.
168 Select from terminal list with filter_terminal
170 a_cs,b_cs,am_bs,bm_cs are symmetric functions for args i,j < k.
171 Their values for args k and large are not relevant.
174 type constraint_system =
178 a_cs : int->int ->float;
179 b_cs: int->int->float;
180 am_cs: int->int->float;
181 bm_cs: int->int->float;
182 js_cs : (int*int) list; (* psorted, both < k_cs *)
185 history_cs : string list;
189 ****************************************
191 ****************************************
195 let string_of_cs cs =
196 let s = string_of_list string_of_int in
197 let k = string_of_int cs.k_cs in
198 let d = string_of_float cs.d_cs in
199 let a = string_of_f cs.a_cs (cs.k_cs) in
200 let b = string_of_f cs.b_cs (cs.k_cs) in
201 let am = string_of_f cs.am_cs (cs.k_cs) in
202 let bm = string_of_f cs.bm_cs (cs.k_cs) in
203 let js = string_of_list string_of_pair cs.js_cs in
204 let lo = s cs.lo_cs in
205 let str = s cs.str_cs in
207 "{\n k=%s;\n d=%s;\n a=%s;\n b=%s;\n am=%s;\n bm=%s;\n js=%s;\n lo=%s;\n str=%s;\n}\n"
208 k d a b am bm js lo str;;
210 let pp_print_cs f cs = pp_print_string f (string_of_cs cs);;
213 (* report(string_of_cs quad_std_cs);; *)
214 (* #install_printer pp_print_cs;; *)
216 let ks_cs cs = (0-- (cs.k_cs -1));;
218 let ks_cart cs = cart2 (ks_cs cs);;
220 let extensional_equality cs cs' =
221 let bk = cs.k_cs = cs'.k_cs in
222 let bd = cs.d_cs = cs'.d_cs in
223 let m r r' = forall (fun (i,j) -> r i j = r' i j) (ks_cart cs) in
224 let ba = m cs.a_cs cs'.a_cs in
225 let bb = m cs.b_cs cs'.b_cs in
226 let bam = m cs.am_cs cs'.am_cs in
227 let bbm = m cs.bm_cs cs'.bm_cs in
228 let ps js = map psort js in
229 let bj = set_eq (ps cs.js_cs) (ps cs'.js_cs) in
230 let bstr = set_eq (cs.str_cs) cs'.str_cs in
231 let blo = set_eq cs.lo_cs cs'.lo_cs in
232 bk && bd && ba && bb && bam && bbm && bstr && blo && bj;;
236 let debug_cs = ref [];;
239 let catch_failure cs_term h =
244 if exists (extensional_equality cs_term) !debug_cs
245 then [cs_term] else [];;
247 let rec rev_assoc_e e a l =
249 (x,y)::t -> if e a y then x else rev_assoc_e e a t
250 | [] -> failwith "find";;
252 let rec build_path e dict cs cs' buf =
253 if (e cs cs') then cs::buf else
254 let v = rev_assoc_e e cs' dict in
255 build_path e dict cs v (cs'::buf);;
257 build_path (=) [(0,4);(4,9);(9,3);(9,7)] 0 7 [];;
259 let rec mk_assoc h cs_term dict cs_init =
260 let e = extensional_equality in
261 if can (Lib.find (e cs_term)) cs_init then dict
263 let out = map h cs_init in
264 let zs = zip cs_init out in
265 let f(z,css) = map (fun cs-> (z,cs)) css in
266 let zss = List.flatten (map f zs) in
267 mk_assoc h cs_term (zss @ dict) (List.flatten out);;
269 let debug_trace h cs_init cs_term =
270 let dict = mk_assoc h cs_term [] [cs_init] in
271 build_path extensional_equality dict cs_init cs_term [];;
274 let report_cs cs = report(string_of_cs cs);;
278 report_cs cs; debug_cs:= [cs]; failwith s;;
282 ****************************************
283 MORE BASIC OPERATIONS
284 ****************************************
290 There are B-fields (hard constraints) and M-fields (soft constraints).
291 The B-fields define the domain of the constraint system s.
292 The M-fields partition the minimizers in B_s.
294 k,a,b, are the primary B-fields
295 d,j are secondary B-fields that define the function tau.
297 lo,str,am,bm are M-fields.
299 js represented as pairs (i,j) with i<j and adjacent.
300 a,b are bounds defining BB_s.
302 am,bm do not affect the stable constraint system s.
303 The deformation lemmas are all partitioning the M-fields.
305 Subdivision partitions the B-fields.
306 Transfer resets the M-fields and changes the B-field.
311 fun ?k:(k=cs.k_cs) ?a:(a=cs.a_cs) ?b:(b=cs.b_cs)
312 ?am:(am=cs.am_cs) ?bm:(bm=cs.bm_cs)
314 ?js:(js=cs.js_cs) ?lo:(lo=cs.lo_cs)
315 ?str:(str=cs.str_cs) ?h:(history=cs.history_cs) () ->
326 history_cs = history;
329 let hist cs s = modify_cs cs ~h:(s::cs.history_cs) ();;
331 (* usage : modify_cs hex_std_cs ~k:4 ();; *)
333 let reset_to_unadorned cs =
334 modify_cs cs ~am:cs.a_cs ~bm:cs.b_cs ~lo:[] ~str:[] ();;
337 let ink k i = (i+1) mod k;;
339 let dek k i = (i+k-1) mod k;;
341 let inc cs i = ink cs.k_cs i;;
343 let dec cs i = dek cs.k_cs i;;
348 let mk_j cs i = psort (i,inc cs i);;
350 let memj cs (i,j) = mem (psort (i,j)) cs.js_cs;;
352 let compact u = length (sortuniq u) = length u;;
354 let ks_cs cs = (0-- (cs.k_cs -1));;
356 let ks_cart cs = cart2 (ks_cs cs);;
359 filter (fun (i,j) -> (i<j) && not (inc cs i = j) && not (inc cs j = i))
363 filter (fun (i,j) -> (i<j)) (ks_cart cs);;
367 let ed = map (fun i -> psort(i,inc cs i)) ks in
370 let all = (fun t -> true);;
372 let htmax cs p = if (mem p cs.lo_cs) then two else twoh0;;
376 let arcmin cs (i,j) = arc (htmax cs i) (htmax cs j) (cs.am_cs i j);;
378 let arcmax cs (i,j) = arc htmin htmin (cs.bm_cs i j);;
383 (cs.b_cs i j > twoh0 or cs.a_cs i j > two) in
384 List.length (filter m (ks_cs cs));;
388 ****************************************
389 BASIC OPS ON CONSTRAINT SYSTEMS
390 ****************************************
396 (i = j) or (2.0 <= cs.a_cs i j) in
398 let k2s = cart2 ks in
399 let bm = forall f k2s or failwith "is_st:bm" in
400 let bs = zero_diag cs.a_cs ks or failwith "is_st:bs" in
401 let g i = (cs.b_cs i (inc cs i) <= cstab) in
402 let bs' = forall g ks or failwith "is_st:bs'" in
403 let fj (i,j) = not(memj cs (i,j)) or
404 (sqrt8 = cs.a_cs i j && cs.b_cs i j = cstab) in
405 let bj = forall fj k2s or failwith "is_st:bj" in
406 bm && bs && bs' && bj;;
408 let is_tri_stable cs =
409 let bk = (3 = cs.k_cs) or failwith "ts:3" in
411 (i = j) or (2.0 <= cs.a_cs i j) in
413 let k2s = cart2 ks in
414 let bm = forall f k2s or failwith "ts:bm" in
415 let bs = zero_diag cs.a_cs ks or failwith "ts:bs" in
416 let g i = (cs.b_cs i (inc cs i) < four) in
417 let bs' = forall g ks or failwith "ts:bs'" in
418 let fj (i,j) = not(memj cs (i,j)) or
419 (sqrt8 = cs.a_cs i j && cs.b_cs i j = cstab) in
420 let bj = forall fj k2s or failwith "ts:bj" in
421 bk && bm && bs && bs' && bj;;
424 let is_unadorned cs =
425 let f (i,j) = (cs.a_cs i j = cs.am_cs i j) && (cs.b_cs i j = cs.bm_cs i j) in
426 let k2s = ks_cart cs in
427 let bm = forall f k2s (* or failwith "unadorned" *) in
428 bm && (cs.lo_cs=[]) && (cs.str_cs=[]);;
430 let zero_diag a xs = forall (fun i -> ((=) zero) (a i i)) xs;;
432 let stable_constraint_system_test1 cs =
434 let k2s = ks_cart cs in
436 (cs.a_cs i j = cs.a_cs j i) &&
437 (cs.b_cs i j = cs.b_cs j i) &&
438 (cs.am_cs i j = cs.am_cs j i) &&
439 (cs.bm_cs i j = cs.bm_cs j i) &&
440 (cs.a_cs i j <= cs.am_cs i j) &&
441 (cs.am_cs i j <= cs.bm_cs i j) &&
442 (cs.bm_cs i j <= cs.b_cs i j) in
443 let fj (i,j) = (i < j) && ((inc cs i = j) or (inc cs j = i)) in
444 let ls = cs.lo_cs @ cs.str_cs in
445 let _ = ((3 <= cs.k_cs) && (cs.k_cs <= 6) ) or failwith ( "test1:k_bound") in
446 let _ = (subset ls ks && subset cs.js_cs k2s) or failwith "test1:I_lo, I_str, J" in
447 let _ = (cs.d_cs < 0.9) or failwith "test1:d_bound" in
448 let _ = forall f k2s or failwith ("test1:ab_relations") in
449 let _ = (compact cs.js_cs) or failwith "test1:J set" in
450 let _ = (forall fj cs.js_cs) or failwith "test1:J edge" in
451 let _ = (length cs.js_cs + cs.k_cs <= 6) or failwith "test1:card_J" in
454 let stable_constraint_system_test2 cs =
456 let k2s = ks_cart cs in
458 (i = j) or (2.0 <= cs.a_cs i j) in
459 let _ = forall f2 k2s or failwith "test2:a packing" in
460 let _ = zero_diag cs.a_cs ks or failwith "test2:a diagonal" in
461 let b_edge i = cs.b_cs i (inc cs i) in
462 let b_bound i = if (cs.k_cs > 3) then (b_edge i <= cstab) else (b_edge i < four) in
463 let _ = forall b_bound ks or failwith "test2:b edge bounds" in
464 let j_bound (i,j) = not(memj cs (i,j)) or
465 (sqrt8 = cs.a_cs i j && cs.b_cs i j = cstab) in
466 let _ = forall j_bound k2s or failwith "test2:J bounds" in
467 let edge (i,j) = (j = inc cs i) or (i = inc cs j) in (* bug corrected 2/23/2013, was (i = inc c i) *)
468 let fm (i,j) = edge(i,j) && (i < j) && (two < cs.a_cs i j or twoh0< cs.b_cs i j) in
469 let m = length (filter fm k2s) in
470 let _ = (m + cs.k_cs <= 6) or failwith "test2:m" in
473 let assert_stable_cs cs =
475 stable_constraint_system_test1 cs;
476 stable_constraint_system_test2 cs
477 with Failure s -> failcs(cs,s);;
480 (* f is an edge preserving map *)
484 let _ = forall (fun i-> f(g i) = i) ks or failwith "map_cs:not inverse" in
485 let _ = forall (fun i-> g(f i) = i) ks or failwith "map_cs:not inverse" in
486 let a' a i j = a (f i) (f j) in
487 let g2 (i,j) = psort (g i , g j ) in
493 am_cs = a' (cs.am_cs);
494 bm_cs = a' (cs.bm_cs);
495 lo_cs = map g cs.lo_cs;
496 str_cs = map g cs.str_cs;
497 js_cs = map g2 cs.js_cs;
498 history_cs = cs.history_cs;
502 let f i = (cs.k_cs - (i+1)) in
505 let rotate_cs cs = map_cs (inc cs) (dec cs) cs;;
507 let rec rotatek_cs k cs =
508 funpow k rotate_cs cs;;
510 (* anadorned isomorphisms *)
512 let unadorned_iso_strict_cs cs cs' =
513 let cs = reset_to_unadorned cs in
514 let cs' = reset_to_unadorned cs' in
515 let bk = cs.k_cs = cs'.k_cs in
516 let bd = cs.d_cs = cs'.d_cs in
517 let m r r' = forall (fun (i,j) -> r i j = r' i j) (ks_cart cs) in
518 let ba = m cs.a_cs cs'.a_cs in
519 let bb = m cs.b_cs cs'.b_cs in
520 let ps js = map psort js in
521 let bj = set_eq (ps cs.js_cs) (ps cs'.js_cs) in
522 bk && bd && ba && bb && bj;;
524 let unadorned_iso_proper_cs cs cs' =
525 Lib.exists (fun i -> unadorned_iso_strict_cs (rotatek_cs i cs) cs') (ks_cs cs);;
527 let unadorned_iso_cs cs cs' =
528 (cs.k_cs = cs'.k_cs) &&
529 ( unadorned_iso_proper_cs cs cs' or unadorned_iso_proper_cs (opposite_cs cs) cs');;
534 ****************************************
536 ****************************************
539 let cs_adj adj diag k i j =
540 let s t = string_of_int t in
541 if (i<0) or (j<0) or (i>=k) or (j>=k)
542 then failwith ("adj out of range"^(s k)^(s i)^(s j))
543 else if (i=j) then zero
544 else if (j = ink k i) or (i = ink k j) then adj
547 let a_pro pro adj diag k i j =
548 if (i<0) or (j<0) or (i>=k) or (j>=k) then failwith "pro out of range"
549 else if (i=j) then zero
550 else if (i=0 && j=1) or (j=0 && i=1) then pro
551 else if (j = ink k i) or (i = ink k j) then adj
554 (* OLD BUG IF data not sorted:
555 let funlist data d i j =
557 else assocd (psort (i,j)) data d;;
561 let data' = map (fun ((i,j),d) -> (psort(i,j),d)) data in
564 else assocd (psort (i,j)) data' d;;
566 let override a (p,q,d) i j =
567 if psort (p,q) = psort (i,j) then d else a i j;;
569 let overrides a data =
570 let fd = funlist data in
577 ****************************************
578 TABLES OF AUGMENTED CONSTRAINT SYSTEMS
579 ****************************************
582 let mk_unadorned (k,d,a,b,h) = {
595 let hex_std_cs = mk_unadorned (6, d_tame 6,
597 cs_adj twoh0 upperbd 6,"hex std init");;
599 let pent_std_cs = mk_unadorned (5,d_tame 5,
601 cs_adj twoh0 upperbd 5,"pent std init");;
603 let quad_std_cs = mk_unadorned (4,d_tame 4,
605 cs_adj twoh0 upperbd 4,"quad std init");;
607 let tri_std_cs = mk_unadorned (3,d_tame 3,
609 cs_adj twoh0 upperbd 3,"tri std init");;
611 let pent_diag_cs = mk_unadorned (
615 cs_adj twoh0 upperbd 5,"pent diag init");;
617 let quad_diag_cs = mk_unadorned (
621 cs_adj twoh0 upperbd 4,"quad diag init");;
623 let pent_pro_cs = mk_unadorned (
626 a_pro twoh0 two twoh0 5,
627 a_pro sqrt8 twoh0 upperbd 5,"pent pro init");;
629 let quad_pro_cs = mk_unadorned (
632 a_pro twoh0 two sqrt8 4,
633 a_pro sqrt8 twoh0 upperbd 4,"quad pro init");;
646 map assert_stable_cs init_cs;;
647 (forall is_unadorned init_cs) or failwith "init_cs:unadorned";;
649 (* now for the terminal cases done by interval computer calculation *)
656 cs_adj two upperbd 6,"terminal hex");;
663 cs_adj two upperbd 5,"terminal pent");;
665 (* two cases for the 0.467 bound: all top edges 2 or both diags 3 *)
667 let terminal_adhoc_quad_5691615370 = (* doesn't exist: rhombus side 2, diags >= 3 *)
672 cs_adj two upperbd 4,"terminal");;
674 let terminal_adhoc_quad_9563139965B =
679 cs_adj twoh0 three 4,"terminal");;
681 let terminal_adhoc_quad_4680581274 = mk_unadorned(
684 funlist [(0,1),cstab ; (0,2),cstab ; (1,3),cstab ] two,
685 funlist [(0,1),cstab ; (0,2),upperbd ; (1,3),upperbd ] two,"terminal 4680581274");;
687 let terminal_adhoc_quad_7697147739 = mk_unadorned(
690 funlist [(0,1),sqrt8 ; (0,2),cstab ; (1,3),cstab ] two,
691 funlist [(0,1),sqrt8 ; (0,2),upperbd ; (1,3),upperbd ] two,"terminal 7697147739");;
693 (* special case of tri_492...
694 let terminal_tri_3456082115 = mk_unadorned(
697 funlist [(0,1), cstab; (0,2),twoh0; (1,2),two] two,
698 funlist [(0,1), 3.22; (0,2),twoh0; (1,2),two] two,"terminal 3456082115");;
701 let terminal_tri_7720405539 = mk_unadorned(
703 0.5518 /. 2.0 -. 0.2,
704 funlist [(0,1),cstab; (0,2),twoh0; (1,2),two] two,
705 funlist [(0,1),3.41; (0,2),twoh0; (1,2),two] two,"terminal 7720405539");;
707 let terminal_tri_2739661360 = mk_unadorned(
709 0.5518 /. 2.0 +. 0.2,
710 funlist [(0,1),cstab; (0,2),cstab; (1,2),two] two,
711 funlist [(0,1),3.41; (0,2),cstab; (1,2),two] two,"terminal 2739661360");;
713 (* range increased by combining with previous case *)
715 let terminal_tri_9269152105 = mk_unadorned(
718 funlist [(0,1),cstab; (0,2),cstab; (1,2),two] two,
719 funlist [(0,1),3.62; (0,2),cstab; (1,2),two] two,"terminal 9269152105");;
721 let terminal_tri_4922521904 = mk_unadorned(
724 funlist [(0,1),cstab; (0,2),twoh0; (1,2),two] two,
725 funlist [(0,1),3.339; (0,2),twoh0; (1,2),two] two,"terminal 4922521904");;
729 The interval ineq requires bounds on (0,2) of 3.41 3.634,
730 but we have Delta[3.634,2,2,3.01,2.52,3.01]<0, so the upper bound
732 In the range 3.01--3.41, we can slice and use
733 upper echelon inequalities _772 and _273.
734 I haven't implemented this, but it would be an easy extension of
735 the upper_echelon procedure.
737 OK. Implemented as upper_echelonC.
738 We no longer need quad_163.
739 We use 5512912661 deformation: 3.15/h0 instead
743 let terminal_quad_1637868761 = mk_unadorned(
746 funlist [(0,1),two; (1,2), cstab;
747 (2,3),twoh0; (0,3),two; (0,2),3.41; (1,3),cstab] two,
748 funlist [(0,1),two; (1,2), cstab;
749 (2,3),twoh0; (0,3),two; (0,2),3.634; (1,3),upperbd] two,"terminal 1637868761");;
753 let terminal_quad_1637868761 = mk_unadorned(
756 funlist [(0,1),two; (1,2), cstab;
757 (2,3),twoh0; (0,3),two; (0,2),cstab; (1,3),cstab] two,
758 funlist [(0,1),two; (1,2), cstab;
759 (2,3),twoh0; (0,3),two; (0,2),upperbd; (1,3),upperbd] two,"terminal 1637868761");;
767 a_cs = funlist [(0,1),sqrt8] two;
768 b_cs = funlist [(0,1),cstab] twoh0;
769 am_cs = funlist [(0,1),sqrt8] two;
770 bm_cs = funlist [(0,1),cstab] twoh0;
774 history_cs= ["ear 3603097872"];
777 (* two consequences of usual ear *)
778 let ear_jnull = modify_cs ear_cs ~js:[] ();;
780 let ear_stab = mk_unadorned(
783 funlist [(0,2),cstab] two,
784 funlist [(0,2),cstab] twoh0,
788 let terminal_ear_3603097872 = ear_cs;;
790 let terminal_tri_5405130650 =
793 d_cs = 0.477 -. 0.11;
794 a_cs = funlist [(0,1),sqrt8;(0,2),twoh0;(1,2),two] two;
795 b_cs = funlist [(0,1),cstab;(0,2),sqrt8;(1,2),twoh0] twoh0;
796 am_cs = funlist [(0,1),sqrt8;(0,2),twoh0;(1,2),two] two;
797 bm_cs = funlist [(0,1),cstab;(0,2),sqrt8;(1,2),twoh0] twoh0;
801 history_cs=["terminal 5405130650"];
804 let terminal_tri_5766053833 =
807 d_cs = 0.696 -. 2.0 *. 0.11;
808 a_cs = funlist [(0,1),sqrt8;(1,2),sqrt8] two;
809 b_cs = funlist [(0,1),cstab;(1,2),cstab] two;
810 am_cs = funlist [(0,1),sqrt8;(1,2),sqrt8] two;
811 bm_cs = funlist [(0,1),cstab;(1,2),cstab] two;
812 js_cs = [(0,1);(1,2)];
815 history_cs=["terminal 5766053833"];
818 let terminal_tri_5026777310a = mk_unadorned(
820 0.6548 -. 2.0 *. 0.11,
821 funlist [(0,1),sqrt8;(1,2),sqrt8] two,
822 funlist [(0,1),cstab;(1,2),cstab] twoh0,"terminal 5026777310a");;
824 let terminal_tri_7881254908 = mk_unadorned(
826 0.696 -. 2.0 *. 0.11,
827 funlist [(0,1),sqrt8;(1,2),sqrt8] twoh0,
828 funlist [(0,1),cstab;(1,2),cstab] twoh0,"terminal 7881254908");;
830 let terminal_std_tri_OMKYNLT_3336871894 = mk_unadorned(
834 funlist [] two,"terminal 3336871894");;
839 let terminal_std_tri_OMKYNLT_2_1 = mk_unadorned(
842 funlist [(0,1),twoh0] two,
843 funlist [(0,1),twoh0] two,"terminal 1107929058");;
845 let terminal_std_tri_7645170609 = mk_unadorned(
848 funlist [(0,1),sqrt8] two,
849 funlist [(0,1),sqrt8] two,"terminal 7645170609");;
852 let terminal_std_tri_OMKYNLT_1_2 = mk_unadorned(
855 funlist [(0,1),two] twoh0,
856 funlist [(0,1),two] twoh0,"terminal 1532755966");;
858 let terminal_std_tri_7097350062a = mk_unadorned(
861 funlist [(0,1),twoh0;(0,2),sqrt8] two,
862 funlist [(0,1),twoh0;(0,2),sqrt8] two,"terminal 7097350062a");;
864 let terminal_std_tri_2900061606 = mk_unadorned(
866 tame_table_d 1 2 +. (tame_table_d 2 1 -. 0.11),
867 funlist [(0,1),twoh0;(0,2),cstab] two,
868 funlist [(0,1),twoh0;(0,2),cstab] two,"terminal 2900061606");;
870 let terminal_std_tri_2200527225 = mk_unadorned(
872 tame_table_d 1 2 +. 2.0*. (tame_table_d 2 1 -. 0.11),
873 funlist [(0,1),two;] sqrt8,
874 funlist [(0,1),two;] sqrt8,"terminal 2200527225");;
876 let terminal_std_tri_3106201101 = mk_unadorned(
878 tame_table_d 1 2 +. 2.0*. (tame_table_d 2 1 -. 0.11),
879 funlist [(0,1),two;(0,2),cstab] sqrt8,
880 funlist [(0,1),two;(0,2),cstab] sqrt8,"terminal 3106201101");;
882 let terminal_std_tri_9816718044 = mk_unadorned(
884 tame_table_d 1 2 +. 2.0*. (tame_table_d 2 1 -. 0.11),
885 funlist [(0,1),two] cstab,
886 funlist [(0,1),two] cstab,"terminal 9816718044");;
888 let terminal_std_tri_1080462150 = mk_unadorned(
890 tame_table_d 0 3 +. 3.0 *.(tame_table_d 2 1 -. 0.11),
892 funlist [] twoh0,"terminal 1080462150");;
894 let terminal_std_tri_4143829594 = mk_unadorned(
896 tame_table_d 0 3 +. 3.0 *.(tame_table_d 2 1 -. 0.11),
897 funlist [(0,1),cstab] twoh0,
898 funlist [(0,1),cstab] twoh0,"terminal 4143829594");;
900 let terminal_std_tri_7459553847 = mk_unadorned(
902 tame_table_d 0 3 +. 3.0 *.(tame_table_d 2 1 -. 0.11),
903 funlist [(0,1),twoh0] cstab,
904 funlist [(0,1),twoh0] cstab,"terminal 7459553847");;
906 let terminal_std_tri_4528012043 = mk_unadorned(
908 tame_table_d 0 3 +. 3.0 *.(tame_table_d 2 1 -. 0.11),
910 funlist [] cstab,"terminal 4528012043");;
913 (* added 2013-06-04 *)
915 let terminal_tri_4010906068 = mk_unadorned(
919 funlist [] cstab,"terminal 4010906068");;
921 let terminal_tri_6833979866 = mk_unadorned(
924 funlist [(0,1),two] twoh0,
925 funlist [(0,1),twoh0] cstab,"terminal 6833979866");;
927 let terminal_tri_5541487347 = mk_unadorned(
930 funlist [(0,1),twoh0] two,
931 funlist [(0,1),sqrt8] twoh0,"terminal 5541487347");;
936 (* use unit_cs as a default terminal object
938 let unit_cs = terminal_std_tri_OMKYNLT_3336871894;;
945 terminal_adhoc_quad_5691615370;
946 terminal_adhoc_quad_9563139965B;
947 terminal_adhoc_quad_4680581274;
948 terminal_adhoc_quad_7697147739;
950 (* upper echelon related *)
951 terminal_tri_7720405539;
952 terminal_tri_2739661360;
953 terminal_tri_9269152105;
954 terminal_tri_4922521904;
957 terminal_ear_3603097872;
959 terminal_tri_5405130650;
962 terminal_tri_5026777310a;
963 terminal_tri_7881254908;
964 terminal_std_tri_OMKYNLT_3336871894;
966 (* added 2013-06-04 *)
967 terminal_tri_4010906068;
968 terminal_tri_6833979866;
969 terminal_tri_5541487347;
971 (* terminal_tri_3456082115; special case of tri_492... *)
972 (* terminal_quad_1637868761; *)
973 (* ear_stab; test 2013-06-03 *)
974 (* terminal_tri_5766053833; *) (* removed 2013-06-03 *)
977 (* removed 2013-06-04
978 terminal_std_tri_OMKYNLT_2_1 ; (* 1107929058 *)
979 terminal_std_tri_7645170609;
980 terminal_std_tri_OMKYNLT_1_2 ; (* 1532755966 *)
981 terminal_std_tri_7097350062a;
982 terminal_std_tri_2900061606;
983 terminal_std_tri_2200527225;
984 terminal_std_tri_3106201101;
985 terminal_std_tri_9816718044;
986 terminal_std_tri_1080462150;
987 terminal_std_tri_4143829594;
988 terminal_std_tri_7459553847;
989 terminal_std_tri_4528012043;
994 let filter_terminal f = filter f terminal_cs;;
996 map ( assert_stable_cs) (filter_terminal all);;
997 (forall ( is_unadorned) (filter_terminal all)) or failwith "terminal_cs:unadorned";;
1000 (cs.k_cs = 3) && (length cs.js_cs = 1) && (unadorned_iso_cs cs ear_cs);;
1004 ****************************************
1005 OPERATIONS AND TRANSFORMATIONS
1006 ****************************************
1009 (* transfer move an M-field to a larger
1010 B-fields resetting M-fields. b5 not used *)
1013 let machine_eps = 1e-08 in
1014 let b1 cs cs' = (cs.k_cs = cs'.k_cs) in
1015 let b2 cs cs' = (cs.d_cs <= cs'.d_cs +. machine_eps) in
1016 let c3 cs cs' (i,j) =
1017 cs.am_cs i j >= cs'.a_cs i j && cs.bm_cs i j <= cs'.b_cs i j in
1018 let b3 cs cs' = forall (c3 cs cs') (allpair cs) in
1020 if is_ear cs then is_ear cs'
1021 else subset (map psort cs'.js_cs) (map psort cs.js_cs) in
1023 subset cs'.lo_cs cs.lo_cs &&
1024 subset cs'.str_cs cs.str_cs in
1027 b1 cs cs' && b2 cs cs' && b3 cs cs' && b4 cs cs' ;;
1029 let proper_transfer_cs cs cs' =
1030 Lib.exists (fun i -> transfer_to (rotatek_cs i cs) cs') (ks_cs cs);;
1032 let equi_transfer_cs cs cs' =
1033 (cs.k_cs = cs'.k_cs) &&
1034 ( proper_transfer_cs cs cs' or proper_transfer_cs (opposite_cs cs) cs');;
1036 let equi_transfer_to_list terminals =
1037 let f1 = map (C equi_transfer_cs) terminals in
1038 fun cs -> exists (fun f -> f cs) f1;;
1040 let rec transfer_union a b =
1043 | a::aas -> if exists (equi_transfer_cs a) b
1044 then transfer_union aas b
1046 let b' = filter (not o (C equi_transfer_cs a)) b in
1047 transfer_union aas (a::b');;
1049 (* optimized version of equi_transfer_to_list *)
1051 let x_equi_transfer_to_list terms =
1052 let machine_eps = 1e-08 in
1053 let _ = not (can (Lib.find (not o is_unadorned)) terms) or failwith "eq:unadorned" in
1054 let has_ear = exists is_ear terms in
1055 let rotates cs = map (fun i -> rotatek_cs i cs) (ks_cs cs) in
1056 let orotates cs = map (fun i -> rotatek_cs i (opposite_cs cs)) (ks_cs cs) in
1057 let props = List.flatten (map rotates terms) in
1058 let oprops = List.flatten (map orotates terms) in
1059 let terms = (props @ oprops) in
1060 let list_le = forall2 (<=) in
1061 let eval_ls f cs = map (fun(i,j)-> f i j) (allpair cs) in
1062 let dump_e cs = (cs.k_cs,cs.d_cs+.machine_eps,
1063 eval_ls cs.a_cs cs,eval_ls cs.b_cs cs,
1065 let dump_m cs = (cs.k_cs,cs.d_cs,eval_ls cs.am_cs cs,eval_ls cs.bm_cs cs,
1066 sortuniq (map psort cs.js_cs)) in
1067 let dump' = setify (map (dump_e) terms) in
1069 if is_ear cs then has_ear
1071 let (k,d,am,bm,js) = dump_m cs in
1072 exists (fun (k',d',a',b',js') ->
1073 (k=k') && (d<=d') && (list_le a' am) &&
1074 (list_le bm b') && (subset (js') (map psort js))) dump';;
1076 (* changes B-field, M remains a minimizer on smaller domain. Type 1 restriction.
1077 The only places b gets modified is with restrict_type1_cs, restrict_type2_cs, subdivision, and
1078 initialization of pents.
1081 let restrict_type1_cs cs =
1082 (* preconditions added 2013-05-15 *)
1084 let jcount = List.length (filter (fun i -> memj cs (i,inc cs i)) (ks_cs cs)) in
1085 let _ = jcount = 0 or 3 < k or failwith "type1 cs 1" in
1086 (modify_cs cs ~b:cs.bm_cs ());;
1088 (* restrict shrinks to B-field down to the M-field,
1089 leaving M-field intact. *)
1091 let restrict_type2_cs cs =
1092 (* 2013-05-14, added preconditions *)
1094 let jcount = List.length (filter (fun i -> memj cs (i,inc cs i)) (ks_cs cs)) in
1095 let _ = jcount = 0 or failwith "type2 cs" in
1096 let _ = forall (fun i -> cs.am_cs i i = zero) (ks_cs cs) or failwith "type2 zero" in
1097 let cs' = modify_cs cs ~a:cs.am_cs ~b:cs.bm_cs () in
1098 let _ = m_count cs' + k <= 6 or failwith "type2 cs 2" in
1101 (* half_slice works on B-fields, resetting M-fields.
1102 This version does not create an ear.
1103 It does just one side. *)
1105 let half_slice cs p q dv =
1106 (* 2013-05-14, long diag preconditions added *)
1108 let _ = (k = 4) or (cs.bm_cs p q <= cstab) or failwith "half_slice long diag" in
1109 let _ = (cs.bm_cs p q < 4.0) or failwith "half_slice long diag" in
1112 let q' = if (q<p) then q+k else q in
1114 let av = cs.am_cs p q in
1115 let bv = cs.bm_cs p q in
1116 let a = override cs.a_cs (p,q,av) in
1117 let b = override cs.b_cs (p,q,bv) in
1118 let _ = (k' >2) or failwith "half_slice underflow" in
1119 let _ = (k' < k) or failwith "half_slice overflow" in
1120 let r i = (i + k - p) mod k in
1121 let s i' = (i' + p) mod k in
1122 let shift f i' j' = f (s i') (s j') in
1123 let cd1 = mk_unadorned (k',dv,shift a, shift b,"slice:"^(join_semi cs.history_cs)) in
1124 let js = map (fun (i,j) -> psort (r i,r j)) (intersect (cart2 (p--q)) cs.js_cs) in
1125 let cd2 = modify_cs cd1 ~js:js () in
1128 (* calculate the value of the new d from the old one,
1129 using standard assumptions about the desired terminal inequalities.
1130 This does cases where no "ear" is involved.
1131 This inequality is to be used when every edge has bounds in
1132 one of the three intervals [2,2h0], [2h0,sqrt8], [sqrt8,3.01] *)
1134 let calc_d_std_earless cs p q =
1139 let q' = if (q<p) then q+k else q in
1141 let _ = (k'=3) or raise Not_found in
1142 let av = cs.am_cs p q in
1143 let bv = cs.bm_cs p q in
1144 let a = override cs.a_cs (p,q,av) in
1145 let b = override cs.b_cs (p,q,bv) in
1146 let edge3(i,j) = (sqrt8 <= a i j && b i j <= cstab) in
1147 let edge2(i,j) = (twoh0 <= a i j && b i j <= sqrt8) && (not(edge3(i,j))) in
1148 let edge1(i,j) = (two <= a i j && b i j <= twoh0) && (not (edge2 (i,j))) in
1149 let edge3s(i,j) = (a i j = cstab && b i j = cstab) in
1150 let edgeL(i,j) = (twoh0 <= a i j && b i j <= cstab) in
1151 let edge (i, j) = edge1 (i ,j) or edge2 (i, j) or edge3 (i, j) or edgeL(i,j) in
1152 let (p1,p2,p3) = (p,inc cs p,funpow 2 (inc cs) p) in
1153 let triedge = [(p1,p2);(p2,p3);(p3,p1)] in
1154 let _ = forall edge triedge or failcs (cs, "slice_dstd") in
1155 let c_edge1 = length (filter (edge1) triedge) in
1156 let c_edge2 = length (filter (edge2) triedge) in
1157 let c_edge3= length (filter (edge3) triedge) in
1158 let c_edge3s = length (filter edge3s triedge) in
1159 let c_edgeL = length (filter edgeL triedge) in
1160 let eps = tame_table_d 2 1 -. 0.11 in
1161 let d12 = tame_table_d 1 2 in
1162 let m3 = eps *. float_of_int (c_edge3) in
1163 let flat1 = (c_edge1=2) && (c_edge2=1), tame_table_d 2 1 in
1164 let flat2s = (c_edge1=2) && (c_edge3s=1), djz in
1165 let flat2 = (c_edge1=2) && (c_edge3=1), 0.11 in
1166 let atype = (c_edge1=1) && (c_edge2+c_edge3 =2), d12 +. m3 in
1167 let btype = (c_edge1=0) && (c_edge2+c_edge3=3), tame_table_d 0 3 +. 3.0 *. eps in
1168 let typeL = (c_edge1=1) && (c_edgeL=2),0.2619 in
1169 let (_,dpq) = List.find (fst) [flat1;flat2s;flat2;atype;btype;typeL] in
1172 let sort_slice_order cs p q =
1173 let inc2 = funpow 2 (inc cs) in
1174 let edge1(i,j) = (two <= cs.a_cs i j && cs.b_cs i j <= twoh0) in
1177 (if (inc2 p = q) then (p,q,false) else (q,p,true))
1179 let p1 = inc cs p in
1180 let _ = (inc cs p1 = q) or failwith "sso:not a diag" in
1181 (if (edge1(p,p1) && edge1(p1,q)) then (p,q,false) else (q,p,true));;
1183 let calc_d_std cs p q =
1185 let swap (a,b) = (b,a) in
1186 let inc3 = funpow 3 (inc cs) in
1187 if (k=6) && (inc3 p = q)
1189 let d = cs.d_cs /. two in (d,d)
1191 let (p,q,swapped) = sort_slice_order cs p q in
1194 calc_d_std_earless cs p q
1195 with Not_found -> failwith "calc_d_std" in
1196 (if swapped then swap d else d);;
1198 let slice_cs cs p q dvpq dvqp mk_ear =
1199 let machine_eps = 1e-08 in (* will go away in formalization *)
1200 let _ = dvpq +. dvqp +. machine_eps >= cs.d_cs or failwith "slice_cs:bad d" in
1201 let cpq = half_slice cs p q dvpq in
1202 let cqp = half_slice cs q p dvqp in
1203 let addj cs = modify_cs cs ~js:((0,(cs.k_cs - 1))::cs.js_cs) () in
1204 let cpq = if (mk_ear) then addj cpq else cpq in
1205 let cqp = if (mk_ear) then addj cqp else cqp in
1206 let _ = not(mk_ear) or is_ear cpq or is_ear cqp or failwith "slice_cs:ear" in
1209 let slice_std cs p q =
1210 let (dvpq,dvqp) = calc_d_std cs p q in
1211 let mk_ear = false in
1212 let rv = slice_cs cs p q dvpq dvqp mk_ear in
1213 let _ = forall (fun cs' -> cs'.k_cs < cs.k_cs) rv in
1218 ****************************************
1220 ****************************************
1224 (* deformation acts on M-fields, keeping the domain fixed. *)
1226 let deform_ODXLSTC_cs p cs =
1227 let p0 = dec cs p in
1229 let p2 = inc cs p in
1230 let ks = ks_cs cs in
1231 let ksp = subtract ks [p] in
1232 let diag = subtract ks [p0;p1;p2] in
1233 let _ = mem p ks or failwith "odx:out of range" in
1234 let _ = not(memj cs (p0,p1)) or raise Unchanged in
1235 let _ = not(memj cs (p1,p2)) or raise Unchanged in
1236 let _ = not(mem p cs.lo_cs) or raise Unchanged in
1237 let m q = (cs.a_cs p1 q = cs.bm_cs p1 q) in
1238 let _ = forall (not o m) ksp or raise Unchanged in
1239 let n q = (fourh0 < cs.b_cs p1 q) in
1240 let _ = forall n diag or raise Unchanged in
1241 let cs1 = modify_cs cs ~lo:(sortuniq (p::cs.lo_cs)) () in
1243 if (cs.am_cs p q > cs.a_cs p q) then None
1244 else Some (modify_cs cs ~bm:(override cs.bm_cs (p,q,cs.a_cs p q)) ()) in
1245 cs1::(filter_some(map csp ksp));;
1247 let deform_IMJXPHR_cs p cs =
1248 let p0 = dec cs p in
1250 let p2 = inc cs p in
1251 let ks = ks_cs cs in
1252 let ksp = subtract ks [p] in
1253 let diag = subtract ks [p0;p1;p2] in
1254 let _ = mem p ks or failwith ("imj:out of range"^(string_of_int p)) in
1255 let _ = mem p cs.str_cs or raise Unchanged in
1256 let _ = not(memj cs (p0,p1)) or raise Unchanged in
1257 let _ = not(memj cs (p1,p2)) or raise Unchanged in
1258 let _ = not(mem p cs.lo_cs) or raise Unchanged in
1259 let m q = (cs.a_cs p1 q = cs.bm_cs p1 q) in
1260 let _ = forall (not o m) diag or raise Unchanged in
1261 let _ = (m p0 <> m p2) or raise Unchanged in (* boolean xor *)
1262 let p' = if (m p0) then p0 else p2 in
1263 let ksp' = subtract ksp [p'] in
1264 let n q = (fourh0 < cs.b_cs p1 q) in
1265 let _ = forall n diag or raise Unchanged in
1266 let cs1 = modify_cs cs ~lo:(sortuniq (p::cs.lo_cs)) () in
1268 if (cs.am_cs p q > cs.a_cs p q) then None
1269 else Some(modify_cs cs ~bm:(override cs.bm_cs (p,q,cs.a_cs p q)) ()) in
1270 cs1::(filter_some(map csp ksp'));;
1272 let deform_NUXCOEA_cs p cs =
1273 let p0 = dec cs p in
1275 let p2 = inc cs p in
1276 let ks = ks_cs cs in
1277 let ksp = subtract ks [p] in
1278 let diag = subtract ks [p0;p1;p2] in
1279 let _ = mem p ks or failwith ("nux:out of range"^(string_of_int p)) in
1280 let _ = mem p cs.str_cs or raise Unchanged in
1281 let _ = not(memj cs (p0,p1)) or raise Unchanged in
1282 let _ = not(memj cs (p1,p2)) or raise Unchanged in
1283 let m q = (cs.a_cs p1 q = cs.bm_cs p1 q) in
1284 let _ = forall (not o m) diag or raise Unchanged in
1285 let _ = (m p0 <> m p2) or raise Unchanged in (* boolean xor *)
1286 let p' = if (m p0) then p0 else p2 in
1287 let ksp' = subtract ksp [p'] in
1288 let n q = (fourh0 < cs.b_cs p1 q) in
1289 let _ = forall n diag or raise Unchanged in
1291 if (cs.am_cs p q > cs.a_cs p q) then None
1292 else Some(modify_cs cs ~bm:(override cs.bm_cs (p,q,cs.a_cs p q)) ()) in
1293 (filter_some(map csp ksp'));;
1296 apply M-field deformation at (p,p+1)=(p1,p2)
1297 This assumes not lunar.
1300 let deform_1834976373_A1_single p cs =
1302 let p2 = inc cs p in
1303 let ks = ks_cs cs in
1304 let alldiag' = alldiag cs in
1305 let _ = mem p ks or failwith "1834:out of range" in
1306 let _ = arcmax cs (p1, p2) < arc1553 or
1308 let _ = not(memj cs (p1,p2)) or raise Unchanged in
1309 let _ = not(mem p1 cs.str_cs) or raise Unchanged in
1310 let _ = not(mem p2 cs.str_cs) or raise Unchanged in
1311 let m (i,j) = (cs.a_cs i j = cs.bm_cs i j) in
1312 let _ = forall (not o m) alldiag' or raise Unchanged in
1313 let _ = not (m (p1,p2)) or raise Unchanged in
1314 let m' (i,j) = (cs.b_cs i j = cs.am_cs i j) in
1315 let _ = not (m'(p1,p2)) or raise Unchanged in
1316 let n (i,j) = (fourh0 < cs.b_cs i j) in
1317 let _ = forall n alldiag' or raise Unchanged in
1318 let cs1 = modify_cs cs ~str:(sortuniq (p1::cs.str_cs)) () in
1319 let cs2 = modify_cs cs ~str:(sortuniq (p2::cs.str_cs)) () in
1321 if (cs.am_cs p q > cs.a_cs p q) then None
1322 else Some (modify_cs cs ~bm:(override cs.bm_cs (p,q,cs.a_cs p q)) ()) in
1324 if (cs.bm_cs p q < cs.b_cs p q) then None
1325 else Some (modify_cs cs ~am:(override cs.am_cs (p,q,cs.b_cs p q)) ()) in
1327 (filter_some(cspq' (p1,p2) :: (cspq (p1,p2)) ::(map cspq alldiag')));;
1330 apply M-field deformation at straight p=p1, double edge (p0,p1) (p1,p2)
1331 This assumes not lunar.
1334 let deform_1834976373_A1_double p cs =
1335 let p0 = dec cs p in
1337 let p2 = inc cs p in
1338 let ks = ks_cs cs in
1339 let alldiag = alldiag cs in
1340 let _ = mem p ks or failwith "1834-double:out of range" in
1341 let _ = (arcmax cs (p1,p2) +. arcmax cs (p0,p1) < arc1553) or
1343 let _ = not(memj cs (p1,p2)) or raise Unchanged in
1344 let _ = not(memj cs (p0,p1)) or raise Unchanged in
1345 let _ = not(mem p0 cs.str_cs) or raise Unchanged in
1346 let _ = mem p1 cs.str_cs or raise Unchanged in
1347 let _ = not(mem p2 cs.str_cs) or raise Unchanged in
1348 let m (i,j) = (cs.a_cs i j = cs.bm_cs i j) in
1349 let _ = forall (not o m) alldiag or raise Unchanged in
1350 let _ = not (m (p1,p2) && m(p0,p1)) or raise Unchanged in
1351 let m' (i,j) = (cs.b_cs i j = cs.am_cs i j) in
1352 let _ = not (m'(p1,p2) && m'(p0,p1)) or raise Unchanged in
1353 let n (i,j) = (fourh0 < cs.b_cs i j) in
1354 let _ = forall n alldiag or raise Unchanged in
1355 let cs1 = modify_cs cs ~str:(sortuniq (p0::cs.str_cs)) () in
1356 let cs2 = modify_cs cs ~str:(sortuniq (p2::cs.str_cs)) () in
1358 if (cs.am_cs p q > cs.a_cs p q) then None
1359 else Some (modify_cs cs ~bm:(override cs.bm_cs (p,q,cs.a_cs p q)) ()) in
1361 if (cs.am_cs p1 p2 > cs.a_cs p1 p2 or cs.am_cs p0 p1 > cs.a_cs p0 p1)
1363 else Some (modify_cs cs ~bm:(overrides cs.bm_cs
1364 [(p0,p1),cs.a_cs p0 p1; (p1,p2),cs.a_cs p1 p2]) ()) in
1366 if (cs.bm_cs p1 p2 < cs.b_cs p1 p2 or cs.bm_cs p0 p1 < cs.b_cs p0 p1)
1368 else Some (modify_cs cs ~am:(overrides cs.am_cs
1369 [(p0,p1),cs.b_cs p0 p1; (p1,p2),cs.b_cs p1 p2]) ()) in
1370 cs1::cs2::(filter_some(csmin :: csmax:: (map cspq alldiag)));;
1374 apply M-field deformation at (p,p+1)=(p1,p2)
1378 correction -- originally I had constraints on the straightness
1379 of p0. By the beta lemma for nonreflexive local fans (or rather by the
1380 monotonocity of dih in the opposite edge), the
1381 azimuth angle is decreasing at p0 as p1 pivots in the direction of p2.
1385 let deform_4828966562 p0 p1 p2 cs =
1386 let ks = ks_cs cs in
1387 let diag = subtract ks [p0;p1;p2] in
1388 let _ = mem p1 ks or failwith "482:out of range" in
1389 let _ = (three <= cs.a_cs p0 p2) or raise Unchanged in
1390 let _ = (cs.b_cs p1 p2 <= twoh0) or raise Unchanged in
1391 let _ = (cs.b_cs p0 p1 <= cstab) or raise Unchanged in
1392 let _ = (cs.a_cs p1 p2 < cs.bm_cs p1 p2) or raise Unchanged in
1393 let _ = not(memj cs (p1,p2)) or raise Unchanged in
1394 (* let _ = not(mem p0 cs.str_cs) or raise Unchanged in *)
1395 let _ = not(mem p1 cs.str_cs) or raise Unchanged in
1396 let _ = not(mem p2 cs.str_cs) or raise Unchanged in
1397 let m q = (cs.a_cs p1 q = cs.bm_cs p1 q) in
1398 let _ = forall (not o m) diag or raise Unchanged in
1399 let n q = (fourh0 < cs.b_cs p1 q) in
1400 let _ = forall n diag or raise Unchanged in
1401 (* let cs0 = modify_cs cs ~str:(sortuniq (p0::cs.str_cs)) () in *)
1402 let cs1 = modify_cs cs ~str:(sortuniq (p1::cs.str_cs)) () in
1403 let cs2 = modify_cs cs ~str:(sortuniq (p2::cs.str_cs)) () in
1405 if (cs.am_cs p1 q > cs.a_cs p1 q) then None
1406 else Some (modify_cs cs ~bm:(override cs.bm_cs (p1,q,cs.a_cs p1 q)) ()) in
1407 (* cs0:: *) cs1::cs2::
1408 (filter_some( (cspq p2) ::(map cspq diag)));;
1410 let deform_4828966562A p cs =
1411 let p0 = dec cs p in
1413 let p2 = inc cs p in
1414 deform_4828966562 p0 p1 p2 cs;;
1416 let deform_4828966562B p cs =
1417 let p0 = dec cs p in
1419 let p2 = inc cs p in
1420 deform_4828966562 p2 p1 p0 cs;;
1424 Use obtuse angle at p1 to avoid node straigtenings.
1425 Uses "1117202051" and "4559601669" for obtuseness criteria.
1426 //This one could cause infinite looping because we are deleting
1427 //attributes str_cs at p0 p1.
1428 Only use as a last resort.
1430 Update: looping prevented.
1431 It is a bad choice of words to call these deformations.
1433 deforming, we are listing the constraints that prevent deformation.
1434 The deformation might be a diagonal at its lower bound,
1435 a straight node at p1, or an edge (p1,p2) at its lower bound.
1437 If there are other conditions in the stable constraint system, that is ok.
1438 If mem p2 cs.str_cs, then we don't need to eliminate it, because
1439 we never deform it. We are just saying that if an element of M_s
1440 has a straight node at p2, then it must also have some other constraint
1441 that prevents deformation.
1443 So I have modified the code so not to remove p0 p2 from str_cs.
1447 (* check for obtuse angle at p1, based on delta_x4 calcs in main_estimate_ineq.hl,
1448 117202051, 2449601669, 45596016696. *)
1450 let obtuse_crit cs p0 p1 p2 =
1451 ((cs.bm_cs p0 p1= two) && (mem p2 cs.lo_cs)) or
1452 (mem p0 cs.lo_cs && mem p2 cs.lo_cs) or
1453 ((cs.bm_cs p0 p1=two) && (mem p0 cs.lo_cs) );;
1456 Explanation of obtuse_sloc in the next deformation lemma.
1457 Preconditions: hexagon (k=6) and three straights: p0,p4 and either p3 or p5,
1458 making an effective triangle.
1459 Precondition: bm p1 p2 <= twoh0
1460 The spherical law of cosines arg has been replaced with ineq 8430954724,
1461 which does it as an ineq delta4_y < 0.
1463 Here is the old argument:
1464 new condition for obtuseness by spherical law of cosines.
1465 if p3 p4 are both straight, then the three segments
1466 give cos c <= cos(3.0 *. arc 2.52 2.52 2.) < -0.76 < -0.6.
1467 If there are a total of 3 straights, then we have an effective triangle arclengths a,b,c=opposite.
1468 assume p1 p2 are not straight and the side p1 p2 [2,2.52].
1469 Then (0.206,0.685)=(cos(arc 2. 2. 2.52),cos(arc 2.52 2.52 2.)).
1470 By the sloc, cos c - cos a cos b <= cos c + |cos a cos b|
1471 < cos c + 0.6 < 0. and we are obtuse.
1473 If we try to make p5 straight, we have three consec straights p4,p5,p0.
1474 making a chain of 4 edges. arc[2,2,2.52]*4 > pi, making a circular fan.
1475 Hence the preconditions imply that p3 (not p5) is straight.
1476 It is stated as it is, because we don't know if p0,p1,p2 are a inc sequence
1477 or dec sequence, so only p4 is well-defined.
1481 let deform_4828966562_obtuse p0 p1 p2 cs =
1482 let ks = ks_cs cs in
1483 let diag = subtract ks [p0;p1;p2] in
1484 let _ = mem p1 ks or failwith "482:out of range" in
1485 let _ = (cstab <= cs.a_cs p0 p2) or raise Unchanged in
1486 let _ = (cs.b_cs p1 p2 <= twoh0) or raise Unchanged in
1487 let _ = (cs.b_cs p0 p1 <= twoh0) or raise Unchanged in
1488 let obtuse = obtuse_crit cs p0 p1 p2 in
1490 let p4 = funpow 3 (inc cs) p1 in
1492 (length (sortuniq cs.str_cs) = 3) && (subset [p0;p4] cs.str_cs) &&
1493 (cs.bm_cs p1 p2 <= twoh0) && not (mem p1 cs.str_cs) &&
1494 not (mem p2 cs.str_cs) in
1495 let _ = obtuse or obtuse_sloc or raise Unchanged in
1496 let _ = (cs.a_cs p1 p2 < cs.bm_cs p1 p2) or raise Unchanged in
1497 let _ = not(memj cs (p1,p2)) or raise Unchanged in
1498 let _ = not(mem p1 cs.str_cs) or raise Unchanged in
1499 let m q = (cs.a_cs p1 q = cs.bm_cs p1 q) in
1500 let _ = forall (not o m) diag or raise Unchanged in
1501 let n q = (fourh0 < cs.b_cs p1 q) in
1502 let _ = forall n diag or raise Unchanged in
1503 let cs_str1 = modify_cs cs ~str:(sortuniq (p1::cs.str_cs)) () in
1505 if (cs.am_cs p1 q > cs.a_cs p1 q) then None
1506 else Some (modify_cs cs ~bm:(override cs.bm_cs (p1,q,cs.a_cs p1 q)) ()) in
1507 cs_str1:: filter_some (map cspq (p2::diag));;
1509 let deform_4828966562A_obtuse p cs =
1510 let p0 = dec cs p in
1512 let p2 = inc cs p in
1513 deform_4828966562_obtuse p0 p1 p2 cs;;
1515 let deform_4828966562B_obtuse p cs =
1516 let p0 = dec cs p in
1518 let p2 = inc cs p in
1519 deform_4828966562_obtuse p2 p1 p0 cs;;
1523 range calculations for 6843920790.
1524 2.0 *. arc 2.52 2.52 2.0 > arc 2.0 2.0 2.38;;
1525 2.0 *. arc 2.0 2.0 2.52 < arc1553;;
1526 As the documentation for this inequality indicates, it is
1527 applied to the triangle p4 p1 p2.
1529 This also applies to some quads and triangles,
1530 implemented as separate deformations.
1532 Edited: May 23, 2012.
1533 need to consider more diagonals. For instance (p1,p3) may bottom out.
1537 let deform_6843920790 p1 cs =
1538 let _ = (cs.k_cs = 5) or raise Unchanged in
1539 let ks = ks_cs cs in
1540 let p0 = dec cs p1 in
1541 let p2 = inc cs p1 in
1542 let p3 = inc cs p2 in
1543 let p4 = inc cs p3 in
1544 let _ = mem p1 ks or failwith "684:out of range" in
1545 let _ = (cs.am_cs p1 p2 = cstab && cstab = cs.bm_cs p1 p2)
1546 or raise Unchanged in
1547 let _ = (cstab <= cs.am_cs p1 p4) or raise Unchanged in
1548 let _ = (cstab <= cs.am_cs p4 p2) or raise Unchanged in
1549 let f (i,j) = (cs.bm_cs i j <= twoh0) in
1550 let _ = forall f [(p0,p1);(p2,p3);(p3,p4);(p4,p0)] or raise Unchanged in
1551 let _ = not(memj cs (p1,p2)) or raise Unchanged in
1552 let _ = not(mem p1 cs.str_cs) or raise Unchanged in
1553 let _ = not(mem p2 cs.str_cs) or raise Unchanged in
1554 let m (i,j) = (cs.a_cs i j < cs.bm_cs i j) in
1555 let _ = forall m (alldiag cs) or raise Unchanged in
1556 let _ = m (p1,p2) or raise Unchanged in
1557 let n (i,j) = (fourh0 < cs.b_cs i j) in
1558 let _ = forall n (alldiag cs) or raise Unchanged in
1559 let cs1 = modify_cs cs ~str:(sortuniq (p1::cs.str_cs)) () in
1560 let cs2 = modify_cs cs ~str:(sortuniq (p2::cs.str_cs)) () in
1561 (* diags at p4 can be excluded, their lengths are fixed under def. *)
1562 let non_p4diag = [(p3,p1);(p3,p0);(p0,p2)] in
1564 if (cs.am_cs i j > cs.a_cs i j) then None
1565 else Some (modify_cs cs ~bm:(override cs.bm_cs (i,j,cs.a_cs i j)) ()) in
1567 (filter_some(map cspq ((p1,p2)::non_p4diag)));;
1570 deformation 5512912661 functions in an almost identical manner
1571 to 6843920790 on quads. We take the union of the two deformations
1575 let deform_6843920790_quad p1 p2 p3 p4 cs =
1576 let _ = (cs.k_cs = 4) or raise Unchanged in
1577 let ks = ks_cs cs in
1578 let adj (i,j) = (inc cs i = j) or (inc cs j = i) in
1579 let _ = forall adj [(p1,p2);(p2,p3);(p3,p4);(p4,p1)] or failwith "684q" in
1580 let _ = subset [p1;p2;p3;p4] ks or failwith "684:range" in
1581 let arc238 = arc two two 2.38 in
1582 (* ineq 684 conditions. *)
1583 let a2 = (two <= cs.am_cs p1 p2 && cs.bm_cs p1 p2 <= cstab) in
1584 let b2 = (arc238 <= arcmin cs(p1,p4) && cs.bm_cs p1 p4 <= cstab) in
1585 let c2 = (cstab <= cs.am_cs p2 p4 && arcmax cs(p2,p3)+.arcmax cs(p3,p4) <= arc1553) in
1586 (* ineq 5512912661 conditions. *)
1587 let a2' = (arc238 <= arcmin cs (p1,p2) && cs.bm_cs p1 p2 <= cstab) in
1588 let b2' = (two <= cs.bm_cs p1 p4 && cs.bm_cs p1 p4 <= twoh0) in
1589 let arc315 = arc two two (3.15/.1.26) in
1590 let c2' = (arc315 <= arcmin cs ( p2, p4) &&
1591 ((arcmax cs(p2,p3)+.arcmax cs(p3,p4) <= arc1553) or
1592 (arcmax cs (p2,p1) +.arcmax cs (p1,p4) <= arc1553))) in
1593 let _ = (a2 && b2 && c2) or (a2' && b2' && c2') or raise Unchanged in
1594 let _ = not(memj cs (p1,p2)) or raise Unchanged in
1595 let _ = not(mem p1 cs.str_cs) or raise Unchanged in
1596 let _ = not(mem p2 cs.str_cs) or raise Unchanged in
1597 let m (i,j) = (cs.a_cs i j < cs.bm_cs i j) in
1598 let _ = forall m (alldiag cs) or raise Unchanged in
1599 let _ = m (p1,p2) or raise Unchanged in
1600 let n (i,j) = (fourh0 < cs.b_cs i j) in
1601 let _ = forall n (alldiag cs) or raise Unchanged in
1602 let cs1 = modify_cs cs ~str:(sortuniq (p1::cs.str_cs)) () in
1603 let cs2 = modify_cs cs ~str:(sortuniq (p2::cs.str_cs)) () in
1605 if (cs.am_cs i j > cs.a_cs i j) then None
1606 else Some (modify_cs cs ~bm:(override cs.bm_cs (i,j,cs.a_cs i j)) ()) in
1608 (filter_some( [cspq (p1,p2);cspq (p1,p3)] ));;
1610 let deform_6843920790_tri p1 cs =
1611 let p2 = inc cs p1 in
1612 let p3 = inc cs p2 in
1613 let _ = (cs.k_cs = 3) or raise Unchanged in
1614 let ks = ks_cs cs in
1615 let _ = mem p1 ks or failwith "684:tri" in
1616 let _ = subset [p1;p2;p3] ks or failwith "684:range" in
1617 let arc238 = arc two two 2.38 in
1618 let a2 = (two <= cs.am_cs p1 p2 && cs.bm_cs p1 p2 <= cstab) in
1619 let range(p,q) = (arc238 <= arcmin cs(p,q) && arcmax cs (p,q) <= arc1553) in
1620 let _ = (a2 && range(p3,p1) && range(p3,p2)) or raise Unchanged in
1621 let _ = not(memj cs (p1,p2)) or raise Unchanged in
1622 let _ = not(mem p1 cs.str_cs) or raise Unchanged in
1623 let _ = not(mem p2 cs.str_cs) or raise Unchanged in
1624 let m (i,j) = (cs.a_cs i j < cs.bm_cs i j) in
1625 let _ = m (p1,p2) or raise Unchanged in
1626 let cs1 = modify_cs cs ~str:(sortuniq (p1::cs.str_cs)) () in
1627 let cs2 = modify_cs cs ~str:(sortuniq (p2::cs.str_cs)) () in
1629 if (cs.am_cs i j > cs.a_cs i j) then None
1630 else Some (modify_cs cs ~bm:(override cs.bm_cs (i,j,cs.a_cs i j)) ()) in
1632 (filter_some( [cspq (p1,p2)]));;
1634 let deform_684_quadA p1 cs =
1635 let f j = funpow j (inc cs) in
1639 deform_6843920790_quad p1 p2 p3 p4 cs;;
1641 let deform_684_quadB p1 cs =
1642 let f j = funpow j (dec cs) in
1646 deform_6843920790_quad p1 p2 p3 p4 cs;;
1650 let deformations postfilter k =
1651 let r f p cs = filter postfilter (f p cs) in
1652 let m d = map (r d) (0--(k-1)) in
1653 let rtl d cs = map restrict_type1_cs (d cs) in
1654 let u = [deform_ODXLSTC_cs;
1657 deform_4828966562A_obtuse;
1658 deform_4828966562B_obtuse;
1659 deform_1834976373_A1_single;
1660 deform_1834976373_A1_double;
1666 deform_6843920790_tri] in
1667 map rtl (List.flatten (map m u));;
1669 let name_of_deformation k i =
1671 ["odx";"imj";"nux";"482ao";"482bo";"1834s";
1672 "1834d";"482a";"482b";"684";"684qA";"684qB";"684t"] in
1673 let offset = i mod k in
1675 (List.nth names_hex s) ^ "-" ^ (string_of_int offset);;
1679 ****************************************
1681 ****************************************
1684 (* division acts on B-fields, shrinks domain.
1685 The global minima Ms remain global minima on smaller domain.
1686 We can keep adornments.
1687 Avoid Unchanged. Pass through if c is out of range.
1688 Restricts the domain if possible. *)
1692 Five domain possibilities, starting with [a,am,bm,b]
1694 c <= a --> unchanged
1695 a < c <= am --> cs2:[c,am,bm,b]
1696 am < c < bm --> cs1:[a,am,c,c]; cs2:[c,c,bm,b]
1697 bm <= c < b --> cs1:[a,am,bm,c]
1698 b <= c --> unchanged.
1702 let subdivide_cs p q c cs =
1703 let _ = (cs.k_cs > 3) or failwith "subdiv: triangle" in
1704 let _ = (0 <= p && p< cs.k_cs) or failwith "p out of range divide_cs" in
1705 let _ = (0 <= q && q< cs.k_cs) or failwith "q out of range divide_cs" in
1706 let _ = not(p = q) or failwith "subdiv: p q must be unequal" in
1707 let _ = not( memj cs (p,q)) or failwith "subdiv: (p,q) in J" in
1708 let e = alledge cs in
1709 let (p,q) = psort(p,q) in
1710 let _ = not(mem (p,q) e) or (cs.a_cs p q > two) or (cs.b_cs p q > twoh0) or
1711 failwith "subdiv: cannot subdivide standard edge" in
1712 let _ = mem (p,q) (allpair cs) or failwith "subdivide range" in
1713 let _ = not (c = cs.am_cs p q) or not(mem(p,q) e) or failwith "subdiv: c is am on edge" in
1714 let a = cs.a_cs p q in
1715 let b = cs.b_cs p q in
1716 let am = cs.am_cs p q in
1717 let bm = cs.bm_cs p q in
1718 let cs1 = modify_cs cs ~b:(override cs.b_cs (p,q,c)) () in
1719 let cs2 = modify_cs cs ~a:(override cs.a_cs (p,q,c)) () in
1720 let _ = (a < c && c < b) or
1721 (report ("warning:unchanged subdivide "^(string_of_float c));
1722 debug_cs:= cs::!debug_cs; true) in
1723 if not (a < c && c < b) then [cs]
1724 else if bm <= c then [cs1]
1725 else if c <= am then [cs2]
1726 else (* am < c < bm *)
1727 let cs1' = modify_cs cs1 ~bm:(override cs.bm_cs (p,q,c)) () in
1728 let cs2' = modify_cs cs2 ~am:(override cs.am_cs (p,q,c)) () in
1731 let between_cs c cs (i,j) =
1732 (cs.a_cs i j < c && c < cs.b_cs i j );;
1734 let can_subdivide f_diag c cs =
1735 exists (between_cs c cs) (f_diag cs);;
1737 let find_subdivide_edge f_diag c cs =
1738 let diag = f_diag cs in
1739 let ind = index true (map (between_cs c cs) (diag)) in
1742 let subdivide_all_c_diag f_diag c =
1743 let p = partition (can_subdivide f_diag c) in
1744 let rec sub init term =
1748 let (i,j) = find_subdivide_edge f_diag c cs in
1749 let kss = subdivide_cs i j c cs in
1750 let (u,v) = p kss in
1751 sub (u @ css) (v @ term) in
1753 let (u,v) = p init in
1757 (* claims serve as documentation for
1758 how the calculations are progressing from init_cs to terminal_cs.
1760 To do: make this into a formal deductive system.
1761 With rules of inference given by deformation rules, etc.
1763 a = what we assert to prove at that point,
1764 b=what we leave for later
1766 reduce remaining if r transfers to a.
1767 remove from b those that transfer to terminal.
1768 put residual b back in remaining.
1772 let remaining = ref [];;
1773 remaining := init_cs;; (* arrow *)
1775 let claim_arrow(a,b) =
1776 let r = !remaining in
1777 let transfers_to_a cs = exists (equi_transfer_cs cs) a in
1778 let r' = filter (not o transfers_to_a) r in
1779 let transfers_from_b cs = exists (C equi_transfer_cs cs) b in
1780 let r'' = filter (not o transfers_from_b) r' in
1781 let _ = remaining := (b @ r'') in
1785 ****************************************
1787 ****************************************
1790 (* flow on hexagons.
1792 subdivide all diags at stab.
1794 in cyclic repetition.
1795 setting aside all cs st ?bm diag <= stab.
1796 checking they are assert_stable_cs.
1800 (* claim in this section goes as follows *)
1802 (* todo: construct from generic subdivide *)
1804 let hex_std_preslice_02,hex_std_preslice_03 =
1805 let c1 = subdivide_cs 0 2 cstab hex_std_cs in
1806 let c2 = subdivide_cs 0 3 cstab hex_std_cs in
1807 hist (hd c1) "hex_std_preslice 02" , hist (hd c2) "hex_std_preslice 03";;
1809 let hex_std_postslice = (* claim A *)
1810 let cs = hex_std_preslice_02 in
1812 let d = cs.d_cs -. d' in
1813 let vv = slice_cs cs 2 0 d' d false in
1814 let vv02 = map (C hist "slice_hex_02") vv in
1815 let cs = hex_std_preslice_03 in
1816 let d' = cs.d_cs /. 2.0 in
1817 let vv = slice_cs cs 0 3 d' d' false in
1818 let vv03 = transfer_union (map (C hist "slice_hex_03") vv) [] in
1821 let hex_assumptions = [hex_std_preslice_02;hex_std_preslice_03];;
1823 claim_arrow([hex_std_cs],hex_assumptions);; (* claim B *)
1824 claim_arrow(hex_assumptions,hex_std_postslice);; (* claim A *)
1830 let ok_for_more_hex =
1831 let is_hex cs = (cs.k_cs =6) in
1832 let terminals = filter_terminal is_hex in
1834 let _ = (cs.k_cs = 6) or failwith "ok:6" in
1835 let _ = assert_stable_cs cs in
1836 let bstr = 3 + length (cs.str_cs) <= cs.k_cs in
1839 let p1 = inc cs i in
1840 let p2 = inc cs p1 in
1841 let p3 = inc cs p2 in
1842 not (subset[p0;p1;p2;p3] cs.lo_cs && subset[p1;p2] cs.str_cs) in
1843 let bg = forall generic_at (ks_cs cs) in
1844 let bunfinished = not (exists (transfer_to cs) terminals) in
1845 bstr && bg && bunfinished;;
1847 let hex_deformations = deformations ok_for_more_hex 6;;
1849 let name_of_deformation_hex = name_of_deformation 6;;
1851 let filtered_subdivide postfilter init =
1852 let sub = subdivide_all_c_diag (alldiag) cstab init in
1853 filter postfilter sub;;
1855 let rec apply_first dl cs =
1857 [] -> failwith ( "all deformations fail")
1861 with Unchanged -> apply_first dls cs;;
1863 let preslice_ready cs =
1864 if (cs.k_cs=4) && (cs.d_cs=0.467)
1866 (cs.bm_cs 0 2 = three) && (cs.bm_cs 1 3 = three)
1868 exists (fun (i,j) -> cs.bm_cs i j <= cstab ) (alldiag cs);;
1870 let rec general_loop2 df postfilter c active stab_diags =
1871 if (c <= 0) or length stab_diags > 0 then (active,stab_diags)
1872 else match active with
1873 [] -> ([],stab_diags)
1876 let kss = apply_first df cs in
1877 let (u,v) = partition preslice_ready kss in
1878 let u' = filter postfilter u in
1879 general_loop2 df postfilter (c-1) (v @ css) (u' @ stab_diags)
1881 Failure s -> (active,stab_diags);;
1883 let execute_hexagons() = (* claim B *)
1884 let postfilter = not o (equi_transfer_to_list hex_assumptions) in
1885 let hex_ultra = filtered_subdivide postfilter [hex_std_cs] in
1886 let hex_loop = general_loop2 hex_deformations postfilter in
1887 (([],[]) =hex_loop 200000 hex_ultra []);;
1889 (* if true, it means that all hexagons have been reduced
1890 to the one terminal hexagon, together with two cases of hex_std_preslice
1893 (* time execute_hexagons ();; 117secs.
1894 working in svn:2819,2824m,2829m,2830m,2832m *)
1898 ****************************************
1900 ****************************************
1903 (* flow on pentagons is the same as for hexagons *)
1906 let ok_for_more_pent =
1907 let terminals = filter_terminal (fun cs -> cs.k_cs = 5) in
1909 let _ = (cs.k_cs = 5) or failwith "ok:5" in
1910 let _ = assert_stable_cs cs in
1911 let bstr = 3 + length (cs.str_cs) <= cs.k_cs in
1912 let sph_tri_ineq i =
1913 if (length cs.str_cs < 2) then true
1916 let p1 = inc cs p0 in
1917 let p2 = inc cs p1 in
1918 let p3 = inc cs p2 in
1919 let p4 = inc cs p3 in
1920 if not(subset [p1;p2] cs.str_cs) then true
1922 let e03min = arcmin cs (p0,p1) +.
1923 arcmin cs (p1,p2) +. arcmin cs (p2,p3) in
1924 let e03max = arcmax cs (p3,p4) +. arcmax cs (p4, p0) in
1926 let bunfinished = not (exists (transfer_to cs) terminals) in
1927 bstr && bunfinished && forall sph_tri_ineq (ks_cs cs);;
1930 let pent_deformations = deformations ok_for_more_pent 5;;
1932 let name_of_deformation_pent = name_of_deformation 5;;
1934 (* these are the cases handled in the pent section *)
1937 filter (fun cs ->(5=cs.k_cs)) (init_cs @ hex_std_postslice);;
1939 (* it doesn't matter how our assumptions are generated,
1940 as long as they are eventually discharged. *)
1942 let pent_composite_cs = mk_unadorned (
1945 a_pro two two cstab 5,
1946 a_pro cstab twoh0 upperbd 5,"pent composite") ;;
1948 let pent_assumptions =
1949 let pent_comp_rediag_cs (p,q) =
1950 let cs = pent_composite_cs in
1952 ~b:(override cs.b_cs (p,q,cstab))
1953 ~bm:(override cs.bm_cs (p,q,cstab)) () in
1954 let alld = alldiag pent_std_cs in
1955 let ffh ((p,q), cs) = subdivide_cs p q cstab cs in
1956 let preslices = List.flatten (map ffh (cart alld pent_init)) in
1957 let pent_comb = map pent_comp_rediag_cs alld in
1958 let cstab_preslices = filter preslice_ready (pent_comb @ preslices) in
1959 let union_cstab_preslices = transfer_union cstab_preslices [] in
1960 let pa = map (C hist "preslice") union_cstab_preslices in
1963 (* we introduce pent_composite_ultra_cs
1964 which combines all the pent_init cases,
1965 then we run execute_pentagons to get rid of it too.
1968 let pent_composite_ultra_cs = (* claim C *)
1969 let pl = filtered_subdivide
1970 (not o (equi_transfer_to_list pent_assumptions))
1971 (pent_composite_cs::pent_init) in
1972 transfer_union pl [];;
1975 claim_arrow(pent_init,pent_composite_ultra_cs @pent_assumptions);;
1977 (* pent section main claim *)
1979 let execute_pentagons() = (* claim D *)
1980 let pent_filter = not o (equi_transfer_to_list pent_assumptions) in
1981 let pent_loop = general_loop2 pent_deformations pent_filter in
1982 (([],[])=pent_loop 200000 (pent_composite_ultra_cs) []);;
1986 time execute_pentagons();;
1988 if true, it means that all pentagons have been reduced
1989 to the one terminal pentagon, together with cases of pent_assumptions
1990 worked 36sec. in svn:2821, May 20, 2012.
1994 claim_arrow(pent_composite_ultra_cs,pent_assumptions);; (* claim D *)
1998 ****************************************
2000 ****************************************
2003 (* These are the quadrilaterals that need subdivision in
2004 the 'ultra' range of diagonals >= cstab.
2005 There are two ways of doing this, either by picking the smaller
2006 diagonal (echelon B) or picking the "better" diagonal, even if
2009 This should be used as the last resort on a quadrilateral.
2010 It fails unless the slice reduces to terminal cases.
2013 let delta_am_diag2_neg cs d =
2014 let y01 = cs.am_cs 0 1 in
2015 let y03 = cs.am_cs 0 3 in
2016 let y21 = cs.am_cs 2 1 in
2017 let y23 = cs.am_cs 2 3 in
2018 Sphere_math.delta_y d y01 y03 d y23 y21 < 0.0;;
2020 let delta_am_diag_neg cs p1 d = (* d from p1 to p3, am from p2 to p4 *)
2021 let p2 = inc cs p1 in
2022 let p3 = inc cs p2 in
2023 let p4 = inc cs p3 in
2024 let y12 = cs.am_cs p1 p2 in
2025 let y14 = cs.am_cs p1 p4 in
2026 let y32 = cs.am_cs p3 p2 in
2027 let y34 = cs.am_cs p3 p4 in
2028 let am = cs.am_cs p2 p4 in
2029 Sphere_math.delta_y d y12 y14 am y34 y32 < 0.0;;
2031 let check_echelon_precondition cs =
2033 let _ = (cs.k_cs = 4) or failwith "echelon1" in
2034 let fixed (i,j) = (cs.am_cs i j = cs.bm_cs i j) in
2035 let _ = forall fixed (alledge cs) or failwith "echelon:edge" in
2036 let _ = forall (fun(i,j) ->cs.am_cs i j >= cstab) (alldiag cs) or
2037 failwith "ech:diag" in
2038 let _ = (cs.js_cs = []) or failwith "ech:j" in
2040 with Failure _ -> false;;
2042 let upper_echelonA =
2043 let assumptions = filter_terminal all in
2044 let transfer = x_equi_transfer_to_list assumptions in
2045 fun (p1,(dval,diag)) cs ->
2046 let _ = check_echelon_precondition cs or raise Unchanged in
2047 let _ = delta_am_diag_neg cs p1 diag or raise Unchanged in
2048 let p2 = inc cs p1 in
2049 let p3 = inc cs p2 in
2050 let dval' = cs.d_cs -. dval in
2052 (fun cs -> cs.bm_cs p1 p3 <= diag) (subdivide_cs p1 p3 diag cs) in
2053 let css' = List.flatten
2054 (map (fun cs -> slice_cs cs p1 p3 dval dval' false) css) in
2055 let css'' = filter (not o transfer) css' in
2056 let _ = (css''=[]) or raise Unchanged in
2059 let upper_echelonC =
2060 let assumptions = filter_terminal all in
2061 let transfer = x_equi_transfer_to_list assumptions in
2063 let diag341 = 3.41 in
2064 let _ = check_echelon_precondition cs or raise Unchanged in
2065 let p2 = inc cs p1 in
2066 let p3 = inc cs p2 in
2067 let p4 = inc cs p3 in
2068 let edge_is c (i,j) = (cs.am_cs i j = c) in
2069 let rhs = [(p1,p2);(p2,p3)] in
2070 let lhs = [(p3,p4);(p4,p1)] in
2071 let has_one c e = (length(filter(edge_is c) e) = 1) in
2072 let _ = (has_one two rhs) or raise Unchanged in
2073 let _ = (has_one two lhs) or raise Unchanged in
2074 let _ = (has_one cstab rhs) or raise Unchanged in
2075 let _ = (has_one twoh0 lhs) or (has_one cstab lhs) or raise Unchanged in
2076 let css = subdivide_cs p1 p3 diag341 cs in
2077 let (css1,css2) = partition (fun cs -> cs.bm_cs p1 p3 <= diag341) css in
2078 let _ = (length css1 = 1) or failwith "ech:C" in
2079 let cs1 = restrict_type2_cs (hd css1) in
2080 let dval = 0.4759 in
2081 let dval' = cs.d_cs -. dval in
2082 let css' = slice_cs cs1 p1 p3 dval dval' false in
2083 let css'' = filter (not o transfer) css' in
2084 let _ = (css'' =[]) or raise Unchanged in
2087 let upper_echelonB =
2088 let assumptions = filter_terminal all in
2089 let transfer = x_equi_transfer_to_list assumptions in
2091 let _ = check_echelon_precondition cs or raise Unchanged in
2092 let _ = delta_am_diag2_neg cs diag or raise Unchanged in
2094 (fun cs -> cs.bm_cs 0 2 <= diag or cs.bm_cs 1 3 <= diag)
2095 (subdivide_all_c_diag alldiag diag [cs]) in
2096 let css02,css13 = partition
2097 (fun cs -> cs.bm_cs 0 2 <= diag) css in
2098 let dcases = if (diag=3.41) then [0.0759;0.4759] else [0.2759] in
2099 let op(a,b) = (b,a) in
2100 let dv = map (fun d -> (d,cs.d_cs -. d)) dcases in
2101 let dv' = (dv @ map op dv) in
2102 let f (p,q) (d,d') cs =
2103 let css' = slice_cs cs p q d d' false in
2104 let css'' = filter (not o transfer) css' in
2105 let _ = (css''=[]) or raise Unchanged in
2108 ignore(map (apply_first (map (f(0,2)) dv')) css02);
2109 ignore(map (apply_first (map (f(1,3)) dv')) css13);
2111 with Failure _ -> raise Unchanged;;
2114 (* we get the case data from echelon data in terminal_cs. *)
2115 let dataA = cart (0--3)
2116 [(0.0759,3.41);(0.4759,3.41);(0.2759,3.339);(0.2759,3.62)] in
2117 let dataB = [3.41;3.339;3.62] in
2118 let dataC = (0--3) in
2119 let cases = map upper_echelonA dataA @ map upper_echelonB dataB @
2120 map upper_echelonC dataC in
2122 apply_first cases cs;;
2127 ****************************************
2129 ****************************************
2132 (* this handles quad_pro_cs modulo the return cs's *)
2134 let (quad_477_preslice_short,quad_477_preslice_long) = (* claim E *)
2135 let preslices = subdivide_all_c_diag alldiag cstab [quad_pro_cs] in
2136 let vv = (map (C hist "preslice pro") preslices) in
2137 let p = filter (fun cs -> cs.b_cs 0 2 = cstab)
2138 (subdivide_cs 0 2 cstab quad_pro_cs) in
2139 let ww = transfer_union (p @ vv) [] in
2140 partition preslice_ready ww;;
2142 (* This is the case where ears are required -- finally! *)
2144 let execute_quad_to_ear() = (* claim E' *)
2145 let transfer = x_equi_transfer_to_list (filter_terminal all) in
2146 let _ = (length quad_477_preslice_short = 1) or failwith "handle 477" in
2147 let cs = hd quad_477_preslice_short in
2148 let mk_ear = true in
2149 let inc2 = funpow 2 (inc cs) in
2150 let f i = slice_cs cs i (inc2 i) (0.11) (0.477 -. 0.11) mk_ear in
2151 let ind = filter (can f) (0--3) in
2152 let g i = forall (transfer) (f i) in
2156 (* let quad_remaining = !remaining;;
2157 remaining :=quad_remaining;;
2160 claim_arrow([quad_pro_cs],quad_477_preslice_long);; (* claim E *)
2162 let triquad_assumption =
2167 overrides (cs_adj two cstab 4) [((1,2),twoh0);((2,3),twoh0)],
2168 overrides (cs_adj twoh0 upperbd 4) [((1,2),cstab);((2,3),cstab)],
2169 "triquad assumption"
2174 overrides (cs_adj two cstab 4) [((1,2),twoh0);((0,3),twoh0)],
2175 overrides (cs_adj twoh0 upperbd 4) [((1,2),cstab);((0,3),cstab)],
2176 "triquad assumption"
2181 cs_adj twoh0 cstab 3,
2182 cs_adj cstab upperbd 3,
2189 funlist [(0,1),twoh0] two,
2190 funlist [(0,1),sqrt8] twoh0,
2196 cs_adj twoh0 cstab 3,
2201 claim_arrow([],triquad_assumption);; (* claim - no justification required *)
2203 (* 477 has been handled already.
2204 We make the case involving the ear part of the terminal
2205 set so that we don't need to deal with ears any further.
2206 can make it an assumption *)
2208 let terminalj_cs = quad_477_preslice_short @ (filter_terminal all);;
2211 triquad_assumption @ terminalj_cs;;
2213 (* see calc in main_estimate.hl. If opposite edges are too short
2214 the angle cannot be straight. *)
2215 let can_be_straight_2485876245b cs p1 =
2216 if (not (cs.k_cs = 4)) then true
2218 let p0 = dec cs p1 in
2219 let p2 = inc cs p1 in
2220 let p3 = inc cs p2 in
2221 not (cstab <= cs.am_cs p1 p3 &&
2222 cs.bm_cs p1 p2 <= cstab &&
2223 cs.bm_cs p1 p0 <= cstab &&
2224 cs.bm_cs p2 p3 <= twoh0 &&
2225 cs.bm_cs p3 p0 <= twoh0);;
2227 let ok_for_more_tri_quad assumptions =
2228 let terminal_transfer = x_equi_transfer_to_list assumptions in
2230 let _ = (mem cs.k_cs [3;4]) or failwith "ok:3,4" in
2231 let b467_2485876245a =
2232 if (cs.d_cs=0.467) && (cs.k_cs = 4) &&
2233 forall (fun (i,j)-> cs.bm_cs i j <= twoh0) [(0,1);(1,2);(2,3);(3,0)] &&
2234 forall (fun (i,j)-> cs.am_cs i j >= three) [(0,2);(1,3)]
2235 then ( cs.str_cs = []) else true in
2236 let deltay cs = Sphere_math.delta_y (* if neg then geometric impossibility *)
2237 (cs.bm_cs 0 1) (cs.bm_cs 0 3) (cs.am_cs 0 2)
2238 (cs.bm_cs 2 3) (cs.bm_cs 1 2) (cs.am_cs 1 3) >= 0.0 in
2239 let b4 = (cs.k_cs = 4) in
2240 let b4a = if b4 then deltay cs else true in
2241 let is_477 = b4 && (cs.d_cs=0.477) &&
2242 forall (fun (i,j)-> cs.bm_cs i j <= twoh0) [(0,1);(1,2);(2,3);(3,0)] &&
2243 forall (fun (i,j)-> cs.am_cs i j >= cstab) [(0,2);(1,3)] in
2247 let b477a = cs.str_cs = [] in
2249 not(exists (equi_transfer_cs cs) quad_477_preslice_short) in
2252 let sph_tri_ineq p1 =
2253 let p0 = dec cs p1 in
2254 let p2 = inc cs p1 in
2255 let p3 = inc cs p2 in
2256 if (not (mem p1 cs.str_cs)) then true (* in tri str=[] *)
2258 arcmin cs (p0,p1) +. arcmin cs (p1,p2) <
2259 arcmax cs (p0,p3) +. arcmax cs (p3,p2) in
2260 let sph_tri_ineq2 p1 = (* rather ad hoc to kill a case *)
2261 let p0 = dec cs p1 in
2262 let p2 = inc cs p1 in
2263 let p3 = inc cs p2 in
2264 if (not (mem p1 (intersect cs.str_cs cs.lo_cs))) then true
2266 cs.am_cs p0 p1 < cs.bm_cs p0 p3 or
2267 cs.am_cs p1 p2 < cs.bm_cs p2 p3 in
2268 let _ = assert_stable_cs cs in
2269 let bstr = 3 + length (cs.str_cs) <= cs.k_cs in
2270 let bunfinished = not (terminal_transfer cs) in
2271 bstr && b4a && bunfinished && b467_2485876245a &&
2272 b477 && forall sph_tri_ineq (cs.str_cs) &&
2273 forall sph_tri_ineq2 (intersect cs.str_cs cs.lo_cs) &&
2274 forall (can_be_straight_2485876245b cs) (cs.str_cs);;
2276 let special_quad_init =
2277 quad_diag_cs :: quad_477_preslice_long;;
2279 let execute_special_quads = (* claim F *)
2280 let quad_deformations =
2281 deformations (ok_for_more_tri_quad terminal_quad) 4 in
2282 let terminal_transfer = x_equi_transfer_to_list terminal_quad in
2284 let quad_loop = general_loop2 quad_deformations
2285 (not o terminal_transfer) in
2286 (([],[])=quad_loop 200000 special_quad_init []);;
2289 execute_special_quads();;
2292 claim_arrow(special_quad_init,[]);; (* claim F *)
2294 (* execute_special_quads is true, it means that
2295 all special quads have been reduced
2297 worked in svn:2821, May 20, 2012.
2302 (* now the general case remains
2303 policies and procedure:
2304 - no more use of js_cs and ears.
2305 - determine d by slice_dpq procedure.
2307 The order is extremely important!
2308 1- if the cs transfers to a terminal or is not "ready" for more,
2309 then were are done with it.
2310 2- if stab is between lower and upper bounds of a diag, then subdivide it
2311 3- if sqrt8 is between lower and upper bounds of a diag, then sub
2312 4- if there is a diag [2h0,sqrt8], then slice it.
2313 5- if there is a diag [sqrt8,cstab], then slice it.
2314 6- if a deformation applies, then apply it.
2315 7- do "upper echelon" treatment of quads (subdivisions on edges > 3.01)
2319 let ok_for_more assumptions =
2320 let ok3 = ok_for_more_tri_quad assumptions in
2322 ((cs.k_cs = 6) && (ok_for_more_hex cs)) or
2323 ((cs.k_cs=5) && (ok_for_more_pent cs)) or
2324 ((mem cs.k_cs [3;4]) && (ok3 cs));;
2327 let handle_general_case skip8 assumptions =
2328 let ok = ok_for_more assumptions in
2329 let dff k = deformations (fun cs -> true) k in
2330 let defs = [[];[];[];dff 3;dff 4;dff 5;dff 6] in
2331 let inrange cs a b (p,q) = (a <= cs.am_cs p q && cs.bm_cs p q <= b) in
2332 let allunderstable cs =
2333 if(cs.k_cs < 4) then []
2334 else filter (fun (i,j)-> (cs.b_cs i j <= cstab)) (allpair cs) in
2336 let alld = alldiag cs in
2338 if not(ok cs) then []
2339 (* else if exists (equi_transfer_cs cs) terminal_quad then [] *)
2341 else if not(skip8) && (can_subdivide allunderstable sqrt8 cs)
2342 then (subdivide_all_c_diag allunderstable sqrt8 [cs])
2343 else if (can_subdivide allunderstable twoh0 cs)
2344 then (subdivide_all_c_diag allunderstable twoh0 [cs])
2347 let (p,q) = Lib.find (inrange cs twoh0 sqrt8) alld in
2351 let (p,q) = Lib.find (inrange cs sqrt8 cstab) alld in
2355 if (can_subdivide alldiag cstab cs)
2356 then (subdivide_all_c_diag alldiag cstab [cs])
2359 let dl = List.nth defs k in
2364 with Failure _ -> failcs(cs, "no handler found");;
2367 let handle_loop skip8 assumptions =
2368 let handle_one = handle_general_case skip8 assumptions in
2369 let rec handle_loop_rec c ls =
2374 let v = handle_one cs in
2375 let (a,b) = partition (fun cs -> cs.k_cs > 4) (v @ css) in
2376 let (b',b'') = partition (fun cs -> cs.k_cs = 3) b in
2377 handle_loop_rec (c-1) ( a @ b' @ b'') in
2381 let r_init = !remaining;;
2383 let execute_triquad () = (* claim G*)
2384 let b1 = ([] = (handle_loop false terminal_quad 50000) r_init) in
2385 let b2 = ([] = (handle_loop true terminalj_cs 50000) triquad_assumption) in
2388 claim_arrow(r_init,terminal_quad);; (* claim G *)
2389 claim_arrow(triquad_assumption,terminalj_cs);; (* claim G *)
2390 claim_arrow(terminalj_cs,[]);; (* claim E *)
2392 let all_cases_done = (!remaining = []);;
2395 (time execute_hexagons() &&
2396 time execute_pentagons() &&
2397 time execute_quad_to_ear() &&
2398 time execute_special_quads() &&
2399 time execute_triquad() &&
2409 let handle = handle_general_case true terminalj_cs;;
2413 let cs1 = (hd !debug_cs);;
2415 let cs2 = upper_echelon cs1;;
2416 report_cs (hd cs2);;
2417 let cs3 = hd(handle (hd cs2));;
2421 let hl = time (handle_loop true terminalj_cs 50000) tri_cases_left;;
2422 let cs_term = hd (!debug_cs);;
2423 let cs_init = List.nth quad_cases_left 1;;
2425 let handle1 = handle_general_case true terminalj_cs;;
2426 let handle1' = catch_failure cs_term handle1;;
2427 let kl = debug_trace handle1' cs_init cs_term;;
2429 let cs' = List.nth kl 4;;
2432 slice_std cs' 1 3 ;;
2433 calc_d_std cs' 1 3;;
2435 let kl = handle1 cs1;;
2438 filter (not o (x_equi_transfer_to_list (filter_terminal all)) terminal_quad;;