| author | wenzelm | 
| Tue, 25 Mar 2025 23:05:15 +0100 | |
| changeset 82407 | fcc0f74ac086 | 
| parent 74623 | 9b1d33c7bbcc | 
| permissions | -rw-r--r-- | 
| 41474 | 1  | 
(* Title: HOL/Library/Sum_of_Squares/positivstellensatz_tools.ML  | 
2  | 
Author: Philipp Meyer, TU Muenchen  | 
|
| 32645 | 3  | 
|
| 41474 | 4  | 
Functions for generating a certificate from a positivstellensatz and vice  | 
5  | 
versa.  | 
|
| 32645 | 6  | 
*)  | 
7  | 
||
8  | 
signature POSITIVSTELLENSATZ_TOOLS =  | 
|
9  | 
sig  | 
|
| 58629 | 10  | 
val print_cert: RealArith.pss_tree -> string  | 
11  | 
val read_cert: Proof.context -> string -> RealArith.pss_tree  | 
|
| 32645 | 12  | 
end  | 
13  | 
||
| 58629 | 14  | 
structure Positivstellensatz_Tools : POSITIVSTELLENSATZ_TOOLS =  | 
| 32645 | 15  | 
struct  | 
16  | 
||
| 58629 | 17  | 
(** print certificate **)  | 
18  | 
||
19  | 
local  | 
|
| 32645 | 20  | 
|
21  | 
(* map polynomials to strings *)  | 
|
22  | 
||
23  | 
fun string_of_varpow x k =  | 
|
24  | 
let  | 
|
| 59582 | 25  | 
val term = Thm.term_of x  | 
| 55508 | 26  | 
val name =  | 
27  | 
(case term of  | 
|
28  | 
Free (n, _) => n  | 
|
29  | 
| _ => error "Term in monomial not free variable")  | 
|
| 32645 | 30  | 
in  | 
| 55508 | 31  | 
if k = 1 then name else name ^ "^" ^ string_of_int k  | 
| 32645 | 32  | 
end  | 
33  | 
||
| 55508 | 34  | 
fun string_of_monomial m =  | 
35  | 
if FuncUtil.Ctermfunc.is_empty m then "1"  | 
|
36  | 
else  | 
|
37  | 
let  | 
|
38  | 
val m' = FuncUtil.dest_monomial m  | 
|
39  | 
val vps = fold_rev (fn (x,k) => cons (string_of_varpow x k)) m' []  | 
|
40  | 
in foldr1 (fn (s, t) => s ^ "*" ^ t) vps end  | 
|
| 32645 | 41  | 
|
42  | 
fun string_of_cmonomial (m,c) =  | 
|
| 63208 | 43  | 
if FuncUtil.Ctermfunc.is_empty m then Rat.string_of_rat c  | 
| 63205 | 44  | 
else if c = @1 then string_of_monomial m  | 
| 63208 | 45  | 
else Rat.string_of_rat c ^ "*" ^ string_of_monomial m  | 
| 32645 | 46  | 
|
| 55508 | 47  | 
fun string_of_poly p =  | 
48  | 
if FuncUtil.Monomialfunc.is_empty p then "0"  | 
|
49  | 
else  | 
|
50  | 
let  | 
|
51  | 
val cms = map string_of_cmonomial  | 
|
52  | 
(sort (prod_ord FuncUtil.monomial_order (K EQUAL)) (FuncUtil.Monomialfunc.dest p))  | 
|
| 58629 | 53  | 
in foldr1 (fn (t1, t2) => t1 ^ " + " ^ t2) cms end  | 
54  | 
||
55  | 
||
56  | 
(* print cert *)  | 
|
| 32645 | 57  | 
|
| 32828 | 58  | 
fun pss_to_cert (RealArith.Axiom_eq i) = "A=" ^ string_of_int i  | 
59  | 
| pss_to_cert (RealArith.Axiom_le i) = "A<=" ^ string_of_int i  | 
|
60  | 
| pss_to_cert (RealArith.Axiom_lt i) = "A<" ^ string_of_int i  | 
|
| 63208 | 61  | 
| pss_to_cert (RealArith.Rational_eq r) = "R=" ^ Rat.string_of_rat r  | 
62  | 
| pss_to_cert (RealArith.Rational_le r) = "R<=" ^ Rat.string_of_rat r  | 
|
63  | 
| pss_to_cert (RealArith.Rational_lt r) = "R<" ^ Rat.string_of_rat r  | 
|
| 32828 | 64  | 
| pss_to_cert (RealArith.Square p) = "[" ^ string_of_poly p ^ "]^2"  | 
| 55508 | 65  | 
| pss_to_cert (RealArith.Eqmul (p, pss)) =  | 
66  | 
"([" ^ string_of_poly p ^ "] * " ^ pss_to_cert pss ^ ")"  | 
|
67  | 
| pss_to_cert (RealArith.Sum (pss1, pss2)) =  | 
|
68  | 
      "(" ^ pss_to_cert pss1 ^ " + " ^ pss_to_cert pss2 ^ ")"
 | 
|
69  | 
| pss_to_cert (RealArith.Product (pss1, pss2)) =  | 
|
70  | 
      "(" ^ pss_to_cert pss1 ^ " * " ^ pss_to_cert pss2 ^ ")"
 | 
|
| 32645 | 71  | 
|
| 58629 | 72  | 
in  | 
73  | 
||
74  | 
fun print_cert RealArith.Trivial = "()"  | 
|
75  | 
  | print_cert (RealArith.Cert pss) = "(" ^ pss_to_cert pss ^ ")"
 | 
|
76  | 
| print_cert (RealArith.Branch (t1, t2)) =  | 
|
77  | 
      "(" ^ print_cert t1 ^ " & " ^ print_cert t2 ^ ")"
 | 
|
78  | 
||
79  | 
end  | 
|
| 55508 | 80  | 
|
| 32645 | 81  | 
|
| 58629 | 82  | 
|
83  | 
(** read certificate **)  | 
|
| 32645 | 84  | 
|
| 58629 | 85  | 
local  | 
86  | 
||
87  | 
(* basic parsers *)  | 
|
| 32645 | 88  | 
|
| 
32646
 
962b4354ed90
used standard fold function and type aliases
 
Philipp Meyer 
parents: 
32645 
diff
changeset
 | 
89  | 
val str = Scan.this_string  | 
| 32645 | 90  | 
|
| 55508 | 91  | 
val number =  | 
92  | 
Scan.repeat1 (Scan.one Symbol.is_ascii_digit >> (fn s => ord s - ord "0"))  | 
|
93  | 
>> foldl1 (fn (n, d) => n * 10 + d)  | 
|
| 32645 | 94  | 
|
95  | 
val nat = number  | 
|
| 58629 | 96  | 
val int = Scan.optional (str "~" >> K ~1) 1 -- nat >> op *  | 
| 63201 | 97  | 
val rat = int --| str "/" -- int >> Rat.make  | 
98  | 
val rat_int = rat || int >> Rat.of_int  | 
|
| 32645 | 99  | 
|
| 55508 | 100  | 
|
| 58629 | 101  | 
(* polynomial parsers *)  | 
| 32645 | 102  | 
|
| 
32646
 
962b4354ed90
used standard fold function and type aliases
 
Philipp Meyer 
parents: 
32645 
diff
changeset
 | 
103  | 
fun repeat_sep s f = f ::: Scan.repeat (str s |-- f)  | 
| 32645 | 104  | 
|
105  | 
val parse_id = Scan.one Symbol.is_letter ::: Scan.many Symbol.is_letdig >> implode  | 
|
106  | 
||
| 
32646
 
962b4354ed90
used standard fold function and type aliases
 
Philipp Meyer 
parents: 
32645 
diff
changeset
 | 
107  | 
fun parse_varpow ctxt = parse_id -- Scan.optional (str "^" |-- nat) 1 >>  | 
| 74623 | 108  | 
(fn (x, k) => (Thm.cterm_of ctxt (Free (x, \<^Type>\<open>real\<close>)), k))  | 
| 32645 | 109  | 
|
110  | 
fun parse_monomial ctxt = repeat_sep "*" (parse_varpow ctxt) >>  | 
|
| 33339 | 111  | 
(fn xs => fold FuncUtil.Ctermfunc.update xs FuncUtil.Ctermfunc.empty)  | 
| 32645 | 112  | 
|
113  | 
fun parse_cmonomial ctxt =  | 
|
| 
32646
 
962b4354ed90
used standard fold function and type aliases
 
Philipp Meyer 
parents: 
32645 
diff
changeset
 | 
114  | 
rat_int --| str "*" -- (parse_monomial ctxt) >> swap ||  | 
| 63205 | 115  | 
(parse_monomial ctxt) >> (fn m => (m, @1)) ||  | 
| 
32829
 
671eb46eb0a3
tuned FuncFun and FuncUtil structure in positivstellensatz.ML
 
Philipp Meyer 
parents: 
32828 
diff
changeset
 | 
116  | 
rat_int >> (fn r => (FuncUtil.Ctermfunc.empty, r))  | 
| 32645 | 117  | 
|
118  | 
fun parse_poly ctxt = repeat_sep "+" (parse_cmonomial ctxt) >>  | 
|
| 33339 | 119  | 
(fn xs => fold FuncUtil.Monomialfunc.update xs FuncUtil.Monomialfunc.empty)  | 
| 32645 | 120  | 
|
| 55508 | 121  | 
|
| 58629 | 122  | 
(* positivstellensatz parsers *)  | 
| 32645 | 123  | 
|
124  | 
val parse_axiom =  | 
|
| 32828 | 125  | 
(str "A=" |-- int >> RealArith.Axiom_eq) ||  | 
126  | 
(str "A<=" |-- int >> RealArith.Axiom_le) ||  | 
|
127  | 
(str "A<" |-- int >> RealArith.Axiom_lt)  | 
|
| 32645 | 128  | 
|
129  | 
val parse_rational =  | 
|
| 32828 | 130  | 
(str "R=" |-- rat_int >> RealArith.Rational_eq) ||  | 
131  | 
(str "R<=" |-- rat_int >> RealArith.Rational_le) ||  | 
|
132  | 
(str "R<" |-- rat_int >> RealArith.Rational_lt)  | 
|
| 32645 | 133  | 
|
134  | 
fun parse_cert ctxt input =  | 
|
135  | 
let  | 
|
136  | 
val pc = parse_cert ctxt  | 
|
137  | 
val pp = parse_poly ctxt  | 
|
138  | 
in  | 
|
| 55508 | 139  | 
(parse_axiom ||  | 
140  | 
parse_rational ||  | 
|
141  | 
str "[" |-- pp --| str "]^2" >> RealArith.Square ||  | 
|
142  | 
str "([" |-- pp --| str "]*" -- pc --| str ")" >> RealArith.Eqmul ||  | 
|
143  | 
     str "(" |-- pc --| str "*" -- pc --| str ")" >> RealArith.Product ||
 | 
|
144  | 
     str "(" |-- pc --| str "+" -- pc --| str ")" >> RealArith.Sum) input
 | 
|
| 32645 | 145  | 
end  | 
146  | 
||
147  | 
fun parse_cert_tree ctxt input =  | 
|
148  | 
let  | 
|
149  | 
val pc = parse_cert ctxt  | 
|
150  | 
val pt = parse_cert_tree ctxt  | 
|
151  | 
in  | 
|
| 55508 | 152  | 
(str "()" >> K RealArith.Trivial ||  | 
153  | 
     str "(" |-- pc --| str ")" >> RealArith.Cert ||
 | 
|
154  | 
     str "(" |-- pt --| str "&" -- pt --| str ")" >> RealArith.Branch) input
 | 
|
| 32645 | 155  | 
end  | 
156  | 
||
| 58629 | 157  | 
in  | 
| 55508 | 158  | 
|
| 58629 | 159  | 
fun read_cert ctxt input_str =  | 
| 43946 | 160  | 
Symbol.scanner "Bad certificate" (parse_cert_tree ctxt)  | 
161  | 
(filter_out Symbol.is_blank (Symbol.explode input_str))  | 
|
| 32645 | 162  | 
|
163  | 
end  | 
|
164  | 
||
| 58629 | 165  | 
end  | 
166  |