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