author | nipkow |
Mon, 24 Sep 2018 14:30:09 +0200 | |
changeset 69064 | 5840724b1d71 |
parent 67562 | 2427d3e72b6e |
child 70334 | bba46912379e |
permissions | -rw-r--r-- |
41474 | 1 |
(* Title: HOL/Library/Sum_of_Squares/sum_of_squares.ML |
32949
aa6c470a962a
eliminated slightly odd get/set operations in favour of Unsynchronized.ref;
wenzelm
parents:
32839
diff
changeset
|
2 |
Author: Amine Chaieb, University of Cambridge |
aa6c470a962a
eliminated slightly odd get/set operations in favour of Unsynchronized.ref;
wenzelm
parents:
32839
diff
changeset
|
3 |
Author: Philipp Meyer, TU Muenchen |
32268
d50f0cb67578
Functionality for sum of squares to call a remote csdp prover
Philipp Meyer
parents:
32150
diff
changeset
|
4 |
|
32839 | 5 |
A tactic for proving nonlinear inequalities. |
32268
d50f0cb67578
Functionality for sum of squares to call a remote csdp prover
Philipp Meyer
parents:
32150
diff
changeset
|
6 |
*) |
d50f0cb67578
Functionality for sum of squares to call a remote csdp prover
Philipp Meyer
parents:
32150
diff
changeset
|
7 |
|
32949
aa6c470a962a
eliminated slightly odd get/set operations in favour of Unsynchronized.ref;
wenzelm
parents:
32839
diff
changeset
|
8 |
signature SUM_OF_SQUARES = |
32268
d50f0cb67578
Functionality for sum of squares to call a remote csdp prover
Philipp Meyer
parents:
32150
diff
changeset
|
9 |
sig |
32949
aa6c470a962a
eliminated slightly odd get/set operations in favour of Unsynchronized.ref;
wenzelm
parents:
32839
diff
changeset
|
10 |
datatype proof_method = Certificate of RealArith.pss_tree | Prover of string -> string |
aa6c470a962a
eliminated slightly odd get/set operations in favour of Unsynchronized.ref;
wenzelm
parents:
32839
diff
changeset
|
11 |
val sos_tac: (RealArith.pss_tree -> unit) -> proof_method -> Proof.context -> int -> tactic |
38805 | 12 |
val trace: bool Config.T |
58631 | 13 |
val debug: bool Config.T |
14 |
val trace_message: Proof.context -> (unit -> string) -> unit |
|
15 |
val debug_message: Proof.context -> (unit -> string) -> unit |
|
32332 | 16 |
exception Failure of string; |
32268
d50f0cb67578
Functionality for sum of squares to call a remote csdp prover
Philipp Meyer
parents:
32150
diff
changeset
|
17 |
end |
d50f0cb67578
Functionality for sum of squares to call a remote csdp prover
Philipp Meyer
parents:
32150
diff
changeset
|
18 |
|
41474 | 19 |
structure Sum_of_Squares: SUM_OF_SQUARES = |
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
20 |
struct |
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
21 |
|
33029 | 22 |
val max = Integer.max; |
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
23 |
|
63201 | 24 |
val denominator_rat = Rat.dest #> snd #> Rat.of_int; |
55508 | 25 |
|
32839 | 26 |
fun int_of_rat a = |
63201 | 27 |
(case Rat.dest a of |
55508 | 28 |
(i, 1) => i |
29 |
| _ => error "int_of_rat: not an int"); |
|
30 |
||
31 |
fun lcm_rat x y = |
|
63201 | 32 |
Rat.of_int (Integer.lcm (int_of_rat x) (int_of_rat y)); |
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
33 |
|
32839 | 34 |
fun rat_pow r i = |
35 |
let fun pow r i = |
|
63205 | 36 |
if i = 0 then @1 else |
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
37 |
let val d = pow r (i div 2) |
63205 | 38 |
in d * d * (if i mod 2 = 0 then @1 else r) |
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
39 |
end |
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
40 |
in if i < 0 then pow (Rat.inv r) (~ i) else pow r i end; |
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
41 |
|
32839 | 42 |
fun round_rat r = |
55508 | 43 |
let |
63211 | 44 |
val (a,b) = Rat.dest (abs r) |
55508 | 45 |
val d = a div b |
63211 | 46 |
val s = if r < @0 then ~ o Rat.of_int else Rat.of_int |
55508 | 47 |
val x2 = 2 * (a - (b * d)) |
48 |
in s (if x2 >= b then d + 1 else d) end |
|
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
49 |
|
58631 | 50 |
|
67271 | 51 |
val trace = Attrib.setup_config_bool \<^binding>\<open>sos_trace\<close> (K false); |
52 |
val debug = Attrib.setup_config_bool \<^binding>\<open>sos_debug\<close> (K false); |
|
58631 | 53 |
|
54 |
fun trace_message ctxt msg = |
|
55 |
if Config.get ctxt trace orelse Config.get ctxt debug then tracing (msg ()) else (); |
|
56 |
fun debug_message ctxt msg = if Config.get ctxt debug then tracing (msg ()) else (); |
|
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
57 |
|
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
58 |
exception Sanity; |
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
59 |
|
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
60 |
exception Unsolvable; |
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
61 |
|
32332 | 62 |
exception Failure of string; |
63 |
||
32645
1cc5b24f5a01
sos method generates and uses proof certificates
Philipp Meyer
parents:
32332
diff
changeset
|
64 |
datatype proof_method = |
1cc5b24f5a01
sos method generates and uses proof certificates
Philipp Meyer
parents:
32332
diff
changeset
|
65 |
Certificate of RealArith.pss_tree |
1cc5b24f5a01
sos method generates and uses proof certificates
Philipp Meyer
parents:
32332
diff
changeset
|
66 |
| Prover of (string -> string) |
1cc5b24f5a01
sos method generates and uses proof certificates
Philipp Meyer
parents:
32332
diff
changeset
|
67 |
|
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
68 |
(* Turn a rational into a decimal string with d sig digits. *) |
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
69 |
|
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
70 |
local |
55508 | 71 |
|
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
72 |
fun normalize y = |
63211 | 73 |
if abs y < @1/10 then normalize (@10 * y) - 1 |
74 |
else if abs y >= @1 then normalize (y / @10) + 1 |
|
32839 | 75 |
else 0 |
55508 | 76 |
|
77 |
in |
|
78 |
||
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
79 |
fun decimalize d x = |
63205 | 80 |
if x = @0 then "0.0" |
55508 | 81 |
else |
82 |
let |
|
63211 | 83 |
val y = abs x |
55508 | 84 |
val e = normalize y |
63211 | 85 |
val z = rat_pow @10 (~ e) * y + @1 |
86 |
val k = int_of_rat (round_rat (rat_pow @10 d * z)) |
|
55508 | 87 |
in |
63205 | 88 |
(if x < @0 then "-0." else "0.") ^ |
55508 | 89 |
implode (tl (raw_explode(string_of_int k))) ^ |
90 |
(if e = 0 then "" else "e" ^ string_of_int e) |
|
91 |
end |
|
92 |
||
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
93 |
end; |
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
94 |
|
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
95 |
(* Iterations over numbers, and lists indexed by numbers. *) |
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
96 |
|
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
97 |
fun itern k l f a = |
55508 | 98 |
(case l of |
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
99 |
[] => a |
55508 | 100 |
| h::t => itern (k + 1) t f (f h k a)); |
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
101 |
|
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
102 |
fun iter (m,n) f a = |
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
103 |
if n < m then a |
55508 | 104 |
else iter (m + 1, n) f (f m a); |
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
105 |
|
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
106 |
(* The main types. *) |
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
107 |
|
55508 | 108 |
type vector = int * Rat.rat FuncUtil.Intfunc.table; |
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
109 |
|
55508 | 110 |
type matrix = (int * int) * Rat.rat FuncUtil.Intpairfunc.table; |
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
111 |
|
63205 | 112 |
fun iszero (_, r) = r = @0; |
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
113 |
|
32839 | 114 |
|
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
115 |
(* Vectors. Conventionally indexed 1..n. *) |
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
116 |
|
55508 | 117 |
fun vector_0 n = (n, FuncUtil.Intfunc.empty): vector; |
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
118 |
|
55508 | 119 |
fun dim (v: vector) = fst v; |
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
120 |
|
55508 | 121 |
fun vector_cmul c (v: vector) = |
122 |
let val n = dim v in |
|
63205 | 123 |
if c = @0 then vector_0 n |
63198
c583ca33076a
ad-hoc overloading for standard operations on type Rat.rat;
wenzelm
parents:
60867
diff
changeset
|
124 |
else (n,FuncUtil.Intfunc.map (fn _ => fn x => c * x) (snd v)) |
55508 | 125 |
end; |
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
126 |
|
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
127 |
fun vector_of_list l = |
55508 | 128 |
let val n = length l in |
58635 | 129 |
(n, fold_rev FuncUtil.Intfunc.update (1 upto n ~~ l) FuncUtil.Intfunc.empty): vector |
55508 | 130 |
end; |
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
131 |
|
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
132 |
(* Matrices; again rows and columns indexed from 1. *) |
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
133 |
|
55508 | 134 |
fun dimensions (m: matrix) = fst m; |
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
135 |
|
55508 | 136 |
fun row k (m: matrix) : vector = |
137 |
let val (_, j) = dimensions m in |
|
138 |
(j, |
|
139 |
FuncUtil.Intpairfunc.fold (fn ((i, j), c) => fn a => |
|
140 |
if i = k then FuncUtil.Intfunc.update (j, c) a else a) (snd m) FuncUtil.Intfunc.empty) |
|
141 |
end; |
|
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
142 |
|
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
143 |
(* Monomials. *) |
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
144 |
|
32828 | 145 |
fun monomial_eval assig m = |
63198
c583ca33076a
ad-hoc overloading for standard operations on type Rat.rat;
wenzelm
parents:
60867
diff
changeset
|
146 |
FuncUtil.Ctermfunc.fold (fn (x, k) => fn a => a * rat_pow (FuncUtil.Ctermfunc.apply assig x) k) |
63205 | 147 |
m @1; |
55508 | 148 |
|
32829
671eb46eb0a3
tuned FuncFun and FuncUtil structure in positivstellensatz.ML
Philipp Meyer
parents:
32828
diff
changeset
|
149 |
val monomial_1 = FuncUtil.Ctermfunc.empty; |
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
150 |
|
32828 | 151 |
fun monomial_var x = FuncUtil.Ctermfunc.onefunc (x, 1); |
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
152 |
|
32828 | 153 |
val monomial_mul = |
33002 | 154 |
FuncUtil.Ctermfunc.combine Integer.add (K false); |
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
155 |
|
32839 | 156 |
fun monomial_multidegree m = |
55508 | 157 |
FuncUtil.Ctermfunc.fold (fn (_, k) => fn a => k + a) m 0; |
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
158 |
|
55508 | 159 |
fun monomial_variables m = FuncUtil.Ctermfunc.dom m; |
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
160 |
|
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
161 |
(* Polynomials. *) |
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
162 |
|
32828 | 163 |
fun eval assig p = |
63205 | 164 |
FuncUtil.Monomialfunc.fold (fn (m, c) => fn a => a + c * monomial_eval assig m) p @0; |
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
165 |
|
32829
671eb46eb0a3
tuned FuncFun and FuncUtil structure in positivstellensatz.ML
Philipp Meyer
parents:
32828
diff
changeset
|
166 |
val poly_0 = FuncUtil.Monomialfunc.empty; |
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
167 |
|
32839 | 168 |
fun poly_isconst p = |
55508 | 169 |
FuncUtil.Monomialfunc.fold (fn (m, _) => fn a => FuncUtil.Ctermfunc.is_empty m andalso a) |
170 |
p true; |
|
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
171 |
|
63205 | 172 |
fun poly_var x = FuncUtil.Monomialfunc.onefunc (monomial_var x, @1); |
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
173 |
|
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
174 |
fun poly_const c = |
63205 | 175 |
if c = @0 then poly_0 else FuncUtil.Monomialfunc.onefunc (monomial_1, c); |
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
176 |
|
32828 | 177 |
fun poly_cmul c p = |
63205 | 178 |
if c = @0 then poly_0 |
63198
c583ca33076a
ad-hoc overloading for standard operations on type Rat.rat;
wenzelm
parents:
60867
diff
changeset
|
179 |
else FuncUtil.Monomialfunc.map (fn _ => fn x => c * x) p; |
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
180 |
|
63211 | 181 |
fun poly_neg p = FuncUtil.Monomialfunc.map (K ~) p; |
55508 | 182 |
|
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
183 |
|
32828 | 184 |
fun poly_add p1 p2 = |
63205 | 185 |
FuncUtil.Monomialfunc.combine (curry op +) (fn x => x = @0) p1 p2; |
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
186 |
|
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
187 |
fun poly_sub p1 p2 = poly_add p1 (poly_neg p2); |
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
188 |
|
32828 | 189 |
fun poly_cmmul (c,m) p = |
63205 | 190 |
if c = @0 then poly_0 |
55508 | 191 |
else |
192 |
if FuncUtil.Ctermfunc.is_empty m |
|
63198
c583ca33076a
ad-hoc overloading for standard operations on type Rat.rat;
wenzelm
parents:
60867
diff
changeset
|
193 |
then FuncUtil.Monomialfunc.map (fn _ => fn d => c * d) p |
55508 | 194 |
else |
195 |
FuncUtil.Monomialfunc.fold (fn (m', d) => fn a => |
|
63198
c583ca33076a
ad-hoc overloading for standard operations on type Rat.rat;
wenzelm
parents:
60867
diff
changeset
|
196 |
(FuncUtil.Monomialfunc.update (monomial_mul m m', c * d) a)) p poly_0; |
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
197 |
|
32828 | 198 |
fun poly_mul p1 p2 = |
199 |
FuncUtil.Monomialfunc.fold (fn (m, c) => fn a => poly_add (poly_cmmul (c,m) p2) a) p1 poly_0; |
|
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
200 |
|
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
201 |
fun poly_square p = poly_mul p p; |
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
202 |
|
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
203 |
fun poly_pow p k = |
63205 | 204 |
if k = 0 then poly_const @1 |
55508 | 205 |
else if k = 1 then p |
206 |
else |
|
207 |
let val q = poly_square(poly_pow p (k div 2)) |
|
208 |
in if k mod 2 = 1 then poly_mul p q else q end; |
|
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
209 |
|
32828 | 210 |
fun multidegree p = |
44453 | 211 |
FuncUtil.Monomialfunc.fold (fn (m, _) => fn a => max (monomial_multidegree m) a) p 0; |
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
212 |
|
32828 | 213 |
fun poly_variables p = |
67559 | 214 |
sort Thm.fast_term_ord |
55508 | 215 |
(FuncUtil.Monomialfunc.fold_rev |
67559 | 216 |
(fn (m, _) => union (is_equal o Thm.fast_term_ord) (monomial_variables m)) p []); |
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
217 |
|
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
218 |
(* Conversion from HOL term. *) |
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
219 |
|
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
220 |
local |
67271 | 221 |
val neg_tm = \<^cterm>\<open>uminus :: real \<Rightarrow> _\<close> |
67399 | 222 |
val add_tm = \<^cterm>\<open>(+) :: real \<Rightarrow> _\<close> |
223 |
val sub_tm = \<^cterm>\<open>(-) :: real \<Rightarrow> _\<close> |
|
69064
5840724b1d71
Prefix form of infix with * on either side no longer needs special treatment
nipkow
parents:
67562
diff
changeset
|
224 |
val mul_tm = \<^cterm>\<open>(*) :: real \<Rightarrow> _\<close> |
67271 | 225 |
val inv_tm = \<^cterm>\<open>inverse :: real \<Rightarrow> _\<close> |
67399 | 226 |
val div_tm = \<^cterm>\<open>(/) :: real \<Rightarrow> _\<close> |
227 |
val pow_tm = \<^cterm>\<open>(^) :: real \<Rightarrow> _\<close> |
|
67271 | 228 |
val zero_tm = \<^cterm>\<open>0:: real\<close> |
59582 | 229 |
val is_numeral = can (HOLogic.dest_number o Thm.term_of) |
55508 | 230 |
fun poly_of_term tm = |
231 |
if tm aconvc zero_tm then poly_0 |
|
232 |
else |
|
233 |
if RealArith.is_ratconst tm |
|
234 |
then poly_const(RealArith.dest_ratconst tm) |
|
235 |
else |
|
236 |
(let |
|
237 |
val (lop, r) = Thm.dest_comb tm |
|
238 |
in |
|
239 |
if lop aconvc neg_tm then poly_neg(poly_of_term r) |
|
240 |
else if lop aconvc inv_tm then |
|
241 |
let val p = poly_of_term r in |
|
242 |
if poly_isconst p |
|
243 |
then poly_const(Rat.inv (eval FuncUtil.Ctermfunc.empty p)) |
|
244 |
else error "poly_of_term: inverse of non-constant polyomial" |
|
245 |
end |
|
246 |
else |
|
247 |
(let |
|
248 |
val (opr,l) = Thm.dest_comb lop |
|
249 |
in |
|
250 |
if opr aconvc pow_tm andalso is_numeral r |
|
59582 | 251 |
then poly_pow (poly_of_term l) ((snd o HOLogic.dest_number o Thm.term_of) r) |
55508 | 252 |
else if opr aconvc add_tm |
253 |
then poly_add (poly_of_term l) (poly_of_term r) |
|
254 |
else if opr aconvc sub_tm |
|
255 |
then poly_sub (poly_of_term l) (poly_of_term r) |
|
256 |
else if opr aconvc mul_tm |
|
257 |
then poly_mul (poly_of_term l) (poly_of_term r) |
|
258 |
else if opr aconvc div_tm |
|
259 |
then |
|
260 |
let |
|
32839 | 261 |
val p = poly_of_term l |
262 |
val q = poly_of_term r |
|
55508 | 263 |
in |
264 |
if poly_isconst q |
|
265 |
then poly_cmul (Rat.inv (eval FuncUtil.Ctermfunc.empty q)) p |
|
266 |
else error "poly_of_term: division by non-constant polynomial" |
|
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
267 |
end |
55508 | 268 |
else poly_var tm |
269 |
end handle CTERM ("dest_comb",_) => poly_var tm) |
|
270 |
end handle CTERM ("dest_comb",_) => poly_var tm) |
|
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
271 |
in |
55508 | 272 |
val poly_of_term = fn tm => |
67271 | 273 |
if type_of (Thm.term_of tm) = \<^typ>\<open>real\<close> |
55508 | 274 |
then poly_of_term tm |
275 |
else error "poly_of_term: term does not have real type" |
|
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
276 |
end; |
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
277 |
|
63522 | 278 |
|
279 |
(* String of vector (just a list of space-separated numbers). *) |
|
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
280 |
|
55508 | 281 |
fun sdpa_of_vector (v: vector) = |
282 |
let |
|
283 |
val n = dim v |
|
284 |
val strs = |
|
63205 | 285 |
map (decimalize 20 o (fn i => FuncUtil.Intfunc.tryapplyd (snd v) i @0)) (1 upto n) |
55508 | 286 |
in space_implode " " strs ^ "\n" end; |
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
287 |
|
55508 | 288 |
fun triple_int_ord ((a, b, c), (a', b', c')) = |
289 |
prod_ord int_ord (prod_ord int_ord int_ord) ((a, (b, c)), (a', (b', c'))); |
|
290 |
structure Inttriplefunc = FuncFun(type key = int * int * int val ord = triple_int_ord); |
|
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
291 |
|
63522 | 292 |
|
63523 | 293 |
(* Parse back csdp output. *) |
55508 | 294 |
|
63522 | 295 |
local |
296 |
||
297 |
val decimal_digits = Scan.many1 Symbol.is_ascii_digit |
|
298 |
val decimal_nat = decimal_digits >> (#1 o Library.read_int); |
|
299 |
val decimal_int = decimal_nat >> Rat.of_int; |
|
55508 | 300 |
|
63522 | 301 |
val decimal_sig = |
302 |
decimal_int -- Scan.option (Scan.$$ "." |-- decimal_digits) >> |
|
303 |
(fn (a, NONE) => a |
|
304 |
| (a, SOME bs) => a + Rat.of_int (#1 (Library.read_int bs)) / rat_pow @10 (length bs)); |
|
305 |
||
306 |
fun signed neg parse = $$ "-" |-- parse >> neg || $$ "+" |-- parse || parse; |
|
307 |
val exponent = ($$ "e" || $$ "E") |-- signed ~ decimal_nat; |
|
308 |
||
309 |
val decimal = |
|
310 |
signed ~ decimal_sig -- Scan.optional exponent 0 |
|
311 |
>> (fn (a, b) => a * rat_pow @10 b); |
|
312 |
||
63523 | 313 |
val csdp_output = |
314 |
decimal -- Scan.repeat (Scan.$$ " " |-- Scan.option decimal) --| Scan.many Symbol.not_eof |
|
315 |
>> (fn (a, bs) => vector_of_list (a :: map_filter I bs)); |
|
316 |
||
317 |
in |
|
318 |
||
319 |
fun parse_csdpoutput s = |
|
320 |
Symbol.scanner "Malformed CSDP output" csdp_output (raw_explode s); |
|
321 |
||
63522 | 322 |
end; |
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
323 |
|
53975
22ee3fb9d596
backout c6297fa1031a -- strange parsers are required to make this work;
wenzelm
parents:
53972
diff
changeset
|
324 |
|
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
325 |
(* Try some apparently sensible scaling first. Note that this is purely to *) |
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
326 |
(* get a cleaner translation to floating-point, and doesn't affect any of *) |
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
327 |
(* the results, in principle. In practice it seems a lot better when there *) |
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
328 |
(* are extreme numbers in the original problem. *) |
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
329 |
|
55508 | 330 |
(* Version for (int*int*int) keys *) |
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
331 |
local |
63198
c583ca33076a
ad-hoc overloading for standard operations on type Rat.rat;
wenzelm
parents:
60867
diff
changeset
|
332 |
fun max_rat x y = if x < y then y else x |
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
333 |
fun common_denominator fld amat acc = |
55508 | 334 |
fld (fn (_,c) => fn a => lcm_rat (denominator_rat c) a) amat acc |
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
335 |
fun maximal_element fld amat acc = |
63211 | 336 |
fld (fn (_,c) => fn maxa => max_rat maxa (abs c)) amat acc |
55508 | 337 |
fun float_of_rat x = |
63201 | 338 |
let val (a,b) = Rat.dest x |
55508 | 339 |
in Real.fromInt a / Real.fromInt b end; |
340 |
fun int_of_float x = (trunc x handle Overflow => 0 | Domain => 0) |
|
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
341 |
in |
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
342 |
|
55508 | 343 |
fun tri_scale_then solver (obj:vector) mats = |
344 |
let |
|
63205 | 345 |
val cd1 = fold_rev (common_denominator Inttriplefunc.fold) mats @1 |
346 |
val cd2 = common_denominator FuncUtil.Intfunc.fold (snd obj) @1 |
|
63198
c583ca33076a
ad-hoc overloading for standard operations on type Rat.rat;
wenzelm
parents:
60867
diff
changeset
|
347 |
val mats' = map (Inttriplefunc.map (fn _ => fn x => cd1 * x)) mats |
55508 | 348 |
val obj' = vector_cmul cd2 obj |
63205 | 349 |
val max1 = fold_rev (maximal_element Inttriplefunc.fold) mats' @0 |
350 |
val max2 = maximal_element FuncUtil.Intfunc.fold (snd obj') @0 |
|
63211 | 351 |
val scal1 = rat_pow @2 (20 - int_of_float(Math.ln (float_of_rat max1) / Math.ln 2.0)) |
352 |
val scal2 = rat_pow @2 (20 - int_of_float(Math.ln (float_of_rat max2) / Math.ln 2.0)) |
|
63198
c583ca33076a
ad-hoc overloading for standard operations on type Rat.rat;
wenzelm
parents:
60867
diff
changeset
|
353 |
val mats'' = map (Inttriplefunc.map (fn _ => fn x => x * scal1)) mats' |
55508 | 354 |
val obj'' = vector_cmul scal2 obj' |
355 |
in solver obj'' mats'' end |
|
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
356 |
end; |
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
357 |
|
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
358 |
(* Round a vector to "nice" rationals. *) |
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
359 |
|
63198
c583ca33076a
ad-hoc overloading for standard operations on type Rat.rat;
wenzelm
parents:
60867
diff
changeset
|
360 |
fun nice_rational n x = round_rat (n * x) / n; |
32839 | 361 |
fun nice_vector n ((d,v) : vector) = |
55508 | 362 |
(d, FuncUtil.Intfunc.fold (fn (i,c) => fn a => |
363 |
let val y = nice_rational n c in |
|
63205 | 364 |
if c = @0 then a |
55508 | 365 |
else FuncUtil.Intfunc.update (i,y) a |
366 |
end) v FuncUtil.Intfunc.empty): vector |
|
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
367 |
|
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
368 |
fun dest_ord f x = is_equal (f x); |
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
369 |
|
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
370 |
(* Stuff for "equations" ((int*int*int)->num functions). *) |
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
371 |
|
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
372 |
fun tri_equation_cmul c eq = |
63205 | 373 |
if c = @0 then Inttriplefunc.empty |
63198
c583ca33076a
ad-hoc overloading for standard operations on type Rat.rat;
wenzelm
parents:
60867
diff
changeset
|
374 |
else Inttriplefunc.map (fn _ => fn d => c * d) eq; |
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
375 |
|
55508 | 376 |
fun tri_equation_add eq1 eq2 = |
63205 | 377 |
Inttriplefunc.combine (curry op +) (fn x => x = @0) eq1 eq2; |
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
378 |
|
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
379 |
fun tri_equation_eval assig eq = |
55508 | 380 |
let |
381 |
fun value v = Inttriplefunc.apply assig v |
|
63205 | 382 |
in Inttriplefunc.fold (fn (v, c) => fn a => a + value v * c) eq @0 end; |
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
383 |
|
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
384 |
(* Eliminate all variables, in an essentially arbitrary order. *) |
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
385 |
|
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
386 |
fun tri_eliminate_all_equations one = |
55508 | 387 |
let |
388 |
fun choose_variable eq = |
|
389 |
let val (v,_) = Inttriplefunc.choose eq |
|
390 |
in |
|
391 |
if is_equal (triple_int_ord(v,one)) then |
|
392 |
let |
|
393 |
val eq' = Inttriplefunc.delete_safe v eq |
|
394 |
in |
|
395 |
if Inttriplefunc.is_empty eq' then error "choose_variable" |
|
396 |
else fst (Inttriplefunc.choose eq') |
|
397 |
end |
|
398 |
else v |
|
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
399 |
end |
55508 | 400 |
|
401 |
fun eliminate dun eqs = |
|
402 |
(case eqs of |
|
403 |
[] => dun |
|
404 |
| eq :: oeqs => |
|
405 |
if Inttriplefunc.is_empty eq then eliminate dun oeqs |
|
406 |
else |
|
407 |
let |
|
408 |
val v = choose_variable eq |
|
409 |
val a = Inttriplefunc.apply eq v |
|
410 |
val eq' = |
|
63201 | 411 |
tri_equation_cmul ((Rat.of_int ~1) / a) (Inttriplefunc.delete_safe v eq) |
55508 | 412 |
fun elim e = |
63205 | 413 |
let val b = Inttriplefunc.tryapplyd e v @0 in |
414 |
if b = @0 then e |
|
63211 | 415 |
else tri_equation_add e (tri_equation_cmul (~ b / a) eq) |
55508 | 416 |
end |
417 |
in |
|
418 |
eliminate (Inttriplefunc.update(v, eq') (Inttriplefunc.map (K elim) dun)) |
|
419 |
(map elim oeqs) |
|
420 |
end) |
|
421 |
in |
|
422 |
fn eqs => |
|
423 |
let |
|
424 |
val assig = eliminate Inttriplefunc.empty eqs |
|
425 |
val vs = Inttriplefunc.fold (fn (_, f) => fn a => |
|
426 |
remove (dest_ord triple_int_ord) one (Inttriplefunc.dom f) @ a) assig [] |
|
427 |
in (distinct (dest_ord triple_int_ord) vs,assig) end |
|
428 |
end; |
|
32839 | 429 |
|
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
430 |
(* Multiply equation-parametrized poly by regular poly and add accumulator. *) |
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
431 |
|
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
432 |
fun tri_epoly_pmul p q acc = |
55508 | 433 |
FuncUtil.Monomialfunc.fold (fn (m1, c) => fn a => |
434 |
FuncUtil.Monomialfunc.fold (fn (m2, e) => fn b => |
|
435 |
let |
|
436 |
val m = monomial_mul m1 m2 |
|
437 |
val es = FuncUtil.Monomialfunc.tryapplyd b m Inttriplefunc.empty |
|
438 |
in |
|
439 |
FuncUtil.Monomialfunc.update (m,tri_equation_add (tri_equation_cmul c e) es) b |
|
440 |
end) q a) p acc; |
|
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
441 |
|
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
442 |
(* Hence produce the "relevant" monomials: those whose squares lie in the *) |
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
443 |
(* Newton polytope of the monomials in the input. (This is enough according *) |
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
444 |
(* to Reznik: "Extremal PSD forms with few terms", Duke Math. Journal, *) |
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
445 |
(* vol 45, pp. 363--374, 1978. *) |
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
446 |
(* *) |
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
447 |
(* These are ordered in sort of decreasing degree. In particular the *) |
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
448 |
(* constant monomial is last; this gives an order in diagonalization of the *) |
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
449 |
(* quadratic form that will tend to display constants. *) |
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
450 |
|
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
451 |
(* Diagonalize (Cholesky/LDU) the matrix corresponding to a quadratic form. *) |
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
452 |
|
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
453 |
local |
55508 | 454 |
fun diagonalize n i m = |
455 |
if FuncUtil.Intpairfunc.is_empty (snd m) then [] |
|
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
456 |
else |
55508 | 457 |
let |
63205 | 458 |
val a11 = FuncUtil.Intpairfunc.tryapplyd (snd m) (i,i) @0 |
55508 | 459 |
in |
63205 | 460 |
if a11 < @0 then raise Failure "diagonalize: not PSD" |
461 |
else if a11 = @0 then |
|
55508 | 462 |
if FuncUtil.Intfunc.is_empty (snd (row i m)) |
463 |
then diagonalize n (i + 1) m |
|
464 |
else raise Failure "diagonalize: not PSD ___ " |
|
465 |
else |
|
466 |
let |
|
467 |
val v = row i m |
|
468 |
val v' = |
|
469 |
(fst v, FuncUtil.Intfunc.fold (fn (i, c) => fn a => |
|
63198
c583ca33076a
ad-hoc overloading for standard operations on type Rat.rat;
wenzelm
parents:
60867
diff
changeset
|
470 |
let val y = c / a11 |
63205 | 471 |
in if y = @0 then a else FuncUtil.Intfunc.update (i,y) a |
55508 | 472 |
end) (snd v) FuncUtil.Intfunc.empty) |
473 |
fun upt0 x y a = |
|
63205 | 474 |
if y = @0 then a |
55508 | 475 |
else FuncUtil.Intpairfunc.update (x,y) a |
476 |
val m' = |
|
477 |
((n, n), |
|
478 |
iter (i + 1, n) (fn j => |
|
479 |
iter (i + 1, n) (fn k => |
|
480 |
(upt0 (j, k) |
|
63205 | 481 |
(FuncUtil.Intpairfunc.tryapplyd (snd m) (j, k) @0 - |
482 |
FuncUtil.Intfunc.tryapplyd (snd v) j @0 * |
|
483 |
FuncUtil.Intfunc.tryapplyd (snd v') k @0)))) |
|
55508 | 484 |
FuncUtil.Intpairfunc.empty) |
485 |
in (a11, v') :: diagonalize n (i + 1) m' end |
|
486 |
end |
|
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
487 |
in |
55508 | 488 |
fun diag m = |
489 |
let |
|
490 |
val nn = dimensions m |
|
491 |
val n = fst nn |
|
492 |
in |
|
493 |
if snd nn <> n then error "diagonalize: non-square matrix" |
|
494 |
else diagonalize n 1 m |
|
495 |
end |
|
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
496 |
end; |
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
497 |
|
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
498 |
(* Enumeration of monomials with given multidegree bound. *) |
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
499 |
|
32839 | 500 |
fun enumerate_monomials d vars = |
55508 | 501 |
if d < 0 then [] |
502 |
else if d = 0 then [FuncUtil.Ctermfunc.empty] |
|
503 |
else if null vars then [monomial_1] |
|
504 |
else |
|
505 |
let val alts = |
|
506 |
map_range (fn k => |
|
507 |
let |
|
508 |
val oths = enumerate_monomials (d - k) (tl vars) |
|
509 |
in map (fn ks => if k = 0 then ks else FuncUtil.Ctermfunc.update (hd vars, k) ks) oths end) |
|
510 |
(d + 1) |
|
511 |
in flat alts end; |
|
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
512 |
|
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
513 |
(* Enumerate products of distinct input polys with degree <= d. *) |
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
514 |
(* We ignore any constant input polynomials. *) |
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
515 |
(* Give the output polynomial and a record of how it was derived. *) |
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
516 |
|
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
517 |
fun enumerate_products d pols = |
63205 | 518 |
if d = 0 then [(poly_const @1, RealArith.Rational_lt @1)] |
55508 | 519 |
else if d < 0 then [] |
520 |
else |
|
521 |
(case pols of |
|
63205 | 522 |
[] => [(poly_const @1, RealArith.Rational_lt @1)] |
55508 | 523 |
| (p, b) :: ps => |
524 |
let val e = multidegree p in |
|
525 |
if e = 0 then enumerate_products d ps |
|
526 |
else |
|
527 |
enumerate_products d ps @ |
|
528 |
map (fn (q, c) => (poly_mul p q, RealArith.Product (b, c))) |
|
529 |
(enumerate_products (d - e) ps) |
|
530 |
end) |
|
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
531 |
|
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
532 |
(* Convert regular polynomial. Note that we treat (0,0,0) as -1. *) |
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
533 |
|
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
534 |
fun epoly_of_poly p = |
55508 | 535 |
FuncUtil.Monomialfunc.fold (fn (m, c) => fn a => |
63211 | 536 |
FuncUtil.Monomialfunc.update (m, Inttriplefunc.onefunc ((0, 0, 0), ~ c)) a) |
55508 | 537 |
p FuncUtil.Monomialfunc.empty; |
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
538 |
|
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
539 |
(* String for block diagonal matrix numbered k. *) |
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
540 |
|
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
541 |
fun sdpa_of_blockdiagonal k m = |
55508 | 542 |
let |
543 |
val pfx = string_of_int k ^" " |
|
544 |
val ents = |
|
545 |
Inttriplefunc.fold |
|
546 |
(fn ((b, i, j), c) => fn a => if i > j then a else ((b, i, j), c) :: a) |
|
547 |
m [] |
|
59058
a78612c67ec0
renamed "pairself" to "apply2", in accordance to @{apply 2};
wenzelm
parents:
58635
diff
changeset
|
548 |
val entss = sort (triple_int_ord o apply2 fst) ents |
55508 | 549 |
in |
550 |
fold_rev (fn ((b,i,j),c) => fn a => |
|
551 |
pfx ^ string_of_int b ^ " " ^ string_of_int i ^ " " ^ string_of_int j ^ |
|
552 |
" " ^ decimalize 20 c ^ "\n" ^ a) entss "" |
|
553 |
end; |
|
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
554 |
|
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
555 |
(* SDPA for problem using block diagonal (i.e. multiple SDPs) *) |
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
556 |
|
32268
d50f0cb67578
Functionality for sum of squares to call a remote csdp prover
Philipp Meyer
parents:
32150
diff
changeset
|
557 |
fun sdpa_of_blockproblem nblocks blocksizes obj mats = |
55508 | 558 |
let val m = length mats - 1 |
559 |
in |
|
560 |
string_of_int m ^ "\n" ^ |
|
561 |
string_of_int nblocks ^ "\n" ^ |
|
562 |
(space_implode " " (map string_of_int blocksizes)) ^ |
|
563 |
"\n" ^ |
|
564 |
sdpa_of_vector obj ^ |
|
58635 | 565 |
fold_rev (fn (k, m) => fn a => sdpa_of_blockdiagonal (k - 1) m ^ a) |
566 |
(1 upto length mats ~~ mats) "" |
|
55508 | 567 |
end; |
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
568 |
|
32268
d50f0cb67578
Functionality for sum of squares to call a remote csdp prover
Philipp Meyer
parents:
32150
diff
changeset
|
569 |
(* Run prover on a problem in block diagonal form. *) |
d50f0cb67578
Functionality for sum of squares to call a remote csdp prover
Philipp Meyer
parents:
32150
diff
changeset
|
570 |
|
55508 | 571 |
fun run_blockproblem prover nblocks blocksizes obj mats = |
32268
d50f0cb67578
Functionality for sum of squares to call a remote csdp prover
Philipp Meyer
parents:
32150
diff
changeset
|
572 |
parse_csdpoutput (prover (sdpa_of_blockproblem nblocks blocksizes obj mats)) |
d50f0cb67578
Functionality for sum of squares to call a remote csdp prover
Philipp Meyer
parents:
32150
diff
changeset
|
573 |
|
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
574 |
(* 3D versions of matrix operations to consider blocks separately. *) |
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
575 |
|
63205 | 576 |
val bmatrix_add = Inttriplefunc.combine (curry op +) (fn x => x = @0); |
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
577 |
fun bmatrix_cmul c bm = |
63205 | 578 |
if c = @0 then Inttriplefunc.empty |
63198
c583ca33076a
ad-hoc overloading for standard operations on type Rat.rat;
wenzelm
parents:
60867
diff
changeset
|
579 |
else Inttriplefunc.map (fn _ => fn x => c * x) bm; |
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
580 |
|
63201 | 581 |
val bmatrix_neg = bmatrix_cmul (Rat.of_int ~1); |
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
582 |
|
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
583 |
(* Smash a block matrix into components. *) |
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
584 |
|
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
585 |
fun blocks blocksizes bm = |
55508 | 586 |
map (fn (bs, b0) => |
587 |
let |
|
588 |
val m = |
|
589 |
Inttriplefunc.fold |
|
590 |
(fn ((b, i, j), c) => fn a => |
|
591 |
if b = b0 then FuncUtil.Intpairfunc.update ((i, j), c) a else a) |
|
592 |
bm FuncUtil.Intpairfunc.empty |
|
593 |
val _ = FuncUtil.Intpairfunc.fold (fn ((i, j), _) => fn a => max a (max i j)) m 0 |
|
594 |
in (((bs, bs), m): matrix) end) |
|
595 |
(blocksizes ~~ (1 upto length blocksizes)); |
|
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
596 |
|
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
597 |
(* FIXME : Get rid of this !!!*) |
32268
d50f0cb67578
Functionality for sum of squares to call a remote csdp prover
Philipp Meyer
parents:
32150
diff
changeset
|
598 |
local |
44453 | 599 |
fun tryfind_with msg _ [] = raise Failure msg |
600 |
| tryfind_with _ f (x::xs) = (f x handle Failure s => tryfind_with s f xs); |
|
32839 | 601 |
in |
32268
d50f0cb67578
Functionality for sum of squares to call a remote csdp prover
Philipp Meyer
parents:
32150
diff
changeset
|
602 |
fun tryfind f = tryfind_with "tryfind" f |
d50f0cb67578
Functionality for sum of squares to call a remote csdp prover
Philipp Meyer
parents:
32150
diff
changeset
|
603 |
end |
d50f0cb67578
Functionality for sum of squares to call a remote csdp prover
Philipp Meyer
parents:
32150
diff
changeset
|
604 |
|
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
605 |
(* Positiv- and Nullstellensatz. Flag "linf" forces a linear representation. *) |
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
606 |
|
38805 | 607 |
fun real_positivnullstellensatz_general ctxt prover linf d eqs leqs pol = |
32839 | 608 |
let |
55508 | 609 |
val vars = |
610 |
fold_rev (union (op aconvc) o poly_variables) |
|
611 |
(pol :: eqs @ map fst leqs) [] |
|
612 |
val monoid = |
|
613 |
if linf then |
|
63205 | 614 |
(poly_const @1, RealArith.Rational_lt @1):: |
55508 | 615 |
(filter (fn (p,_) => multidegree p <= d) leqs) |
616 |
else enumerate_products d leqs |
|
617 |
val nblocks = length monoid |
|
618 |
fun mk_idmultiplier k p = |
|
619 |
let |
|
620 |
val e = d - multidegree p |
|
621 |
val mons = enumerate_monomials e vars |
|
622 |
val nons = mons ~~ (1 upto length mons) |
|
623 |
in |
|
624 |
(mons, |
|
625 |
fold_rev (fn (m, n) => |
|
63205 | 626 |
FuncUtil.Monomialfunc.update (m, Inttriplefunc.onefunc ((~k, ~n, n), @1))) |
55508 | 627 |
nons FuncUtil.Monomialfunc.empty) |
628 |
end |
|
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
629 |
|
55508 | 630 |
fun mk_sqmultiplier k (p,_) = |
631 |
let |
|
632 |
val e = (d - multidegree p) div 2 |
|
633 |
val mons = enumerate_monomials e vars |
|
634 |
val nons = mons ~~ (1 upto length mons) |
|
635 |
in |
|
636 |
(mons, |
|
637 |
fold_rev (fn (m1, n1) => |
|
638 |
fold_rev (fn (m2, n2) => fn a => |
|
639 |
let val m = monomial_mul m1 m2 in |
|
640 |
if n1 > n2 then a |
|
641 |
else |
|
642 |
let |
|
63205 | 643 |
val c = if n1 = n2 then @1 else @2 |
55508 | 644 |
val e = FuncUtil.Monomialfunc.tryapplyd a m Inttriplefunc.empty |
645 |
in |
|
646 |
FuncUtil.Monomialfunc.update |
|
647 |
(m, tri_equation_add (Inttriplefunc.onefunc ((k, n1, n2), c)) e) a |
|
648 |
end |
|
649 |
end) nons) nons FuncUtil.Monomialfunc.empty) |
|
650 |
end |
|
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
651 |
|
55508 | 652 |
val (sqmonlist,sqs) = split_list (map2 mk_sqmultiplier (1 upto length monoid) monoid) |
653 |
val (_(*idmonlist*),ids) = split_list (map2 mk_idmultiplier (1 upto length eqs) eqs) |
|
654 |
val blocksizes = map length sqmonlist |
|
655 |
val bigsum = |
|
58635 | 656 |
fold_rev (fn (p, q) => fn a => tri_epoly_pmul p q a) (eqs ~~ ids) |
657 |
(fold_rev (fn ((p, _), s) => fn a => tri_epoly_pmul p s a) (monoid ~~ sqs) |
|
658 |
(epoly_of_poly (poly_neg pol))) |
|
55508 | 659 |
val eqns = FuncUtil.Monomialfunc.fold (fn (_, e) => fn a => e :: a) bigsum [] |
660 |
val (pvs, assig) = tri_eliminate_all_equations (0, 0, 0) eqns |
|
661 |
val qvars = (0, 0, 0) :: pvs |
|
662 |
val allassig = |
|
63205 | 663 |
fold_rev (fn v => Inttriplefunc.update (v, (Inttriplefunc.onefunc (v, @1)))) pvs assig |
55508 | 664 |
fun mk_matrix v = |
665 |
Inttriplefunc.fold (fn ((b, i, j), ass) => fn m => |
|
666 |
if b < 0 then m |
|
667 |
else |
|
63205 | 668 |
let val c = Inttriplefunc.tryapplyd ass v @0 in |
669 |
if c = @0 then m |
|
55508 | 670 |
else Inttriplefunc.update ((b, j, i), c) (Inttriplefunc.update ((b, i, j), c) m) |
671 |
end) |
|
672 |
allassig Inttriplefunc.empty |
|
673 |
val diagents = |
|
674 |
Inttriplefunc.fold |
|
675 |
(fn ((b, i, j), e) => fn a => if b > 0 andalso i = j then tri_equation_add e a else a) |
|
676 |
allassig Inttriplefunc.empty |
|
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
677 |
|
55508 | 678 |
val mats = map mk_matrix qvars |
679 |
val obj = |
|
680 |
(length pvs, |
|
681 |
itern 1 pvs (fn v => fn i => |
|
63205 | 682 |
FuncUtil.Intfunc.updatep iszero (i,Inttriplefunc.tryapplyd diagents v @0)) |
55508 | 683 |
FuncUtil.Intfunc.empty) |
684 |
val raw_vec = |
|
685 |
if null pvs then vector_0 0 |
|
686 |
else tri_scale_then (run_blockproblem prover nblocks blocksizes) obj mats |
|
63205 | 687 |
fun int_element (_, v) i = FuncUtil.Intfunc.tryapplyd v i @0 |
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
688 |
|
55508 | 689 |
fun find_rounding d = |
690 |
let |
|
691 |
val _ = |
|
58631 | 692 |
debug_message ctxt (fn () => "Trying rounding with limit "^Rat.string_of_rat d ^ "\n") |
55508 | 693 |
val vec = nice_vector d raw_vec |
694 |
val blockmat = |
|
695 |
iter (1, dim vec) |
|
696 |
(fn i => fn a => bmatrix_add (bmatrix_cmul (int_element vec i) (nth mats i)) a) |
|
697 |
(bmatrix_neg (nth mats 0)) |
|
698 |
val allmats = blocks blocksizes blockmat |
|
699 |
in (vec, map diag allmats) end |
|
700 |
val (vec, ratdias) = |
|
63205 | 701 |
if null pvs then find_rounding @1 |
63211 | 702 |
else tryfind find_rounding (map Rat.of_int (1 upto 31) @ map (rat_pow @2) (5 upto 66)) |
55508 | 703 |
val newassigs = |
704 |
fold_rev (fn k => Inttriplefunc.update (nth pvs (k - 1), int_element vec k)) |
|
63201 | 705 |
(1 upto dim vec) (Inttriplefunc.onefunc ((0, 0, 0), Rat.of_int ~1)) |
55508 | 706 |
val finalassigs = |
707 |
Inttriplefunc.fold (fn (v, e) => fn a => |
|
708 |
Inttriplefunc.update (v, tri_equation_eval newassigs e) a) allassig newassigs |
|
709 |
fun poly_of_epoly p = |
|
710 |
FuncUtil.Monomialfunc.fold (fn (v, e) => fn a => |
|
711 |
FuncUtil.Monomialfunc.updatep iszero (v, tri_equation_eval finalassigs e) a) |
|
712 |
p FuncUtil.Monomialfunc.empty |
|
713 |
fun mk_sos mons = |
|
714 |
let |
|
715 |
fun mk_sq (c, m) = |
|
716 |
(c, fold_rev (fn k => fn a => |
|
717 |
FuncUtil.Monomialfunc.updatep iszero (nth mons (k - 1), int_element m k) a) |
|
718 |
(1 upto length mons) FuncUtil.Monomialfunc.empty) |
|
719 |
in map mk_sq end |
|
720 |
val sqs = map2 mk_sos sqmonlist ratdias |
|
721 |
val cfs = map poly_of_epoly ids |
|
722 |
val msq = filter (fn (_, b) => not (null b)) (map2 pair monoid sqs) |
|
723 |
fun eval_sq sqs = fold_rev (fn (c, q) => poly_add (poly_cmul c (poly_mul q q))) sqs poly_0 |
|
724 |
val sanity = |
|
725 |
fold_rev (fn ((p, _), s) => poly_add (poly_mul p (eval_sq s))) msq |
|
58635 | 726 |
(fold_rev (fn (p, q) => poly_add (poly_mul p q)) (cfs ~~ eqs) (poly_neg pol)) |
55508 | 727 |
in |
728 |
if not(FuncUtil.Monomialfunc.is_empty sanity) then raise Sanity |
|
729 |
else (cfs, map (fn (a, b) => (snd a, b)) msq) |
|
730 |
end |
|
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
731 |
|
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
732 |
|
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
733 |
(* Iterative deepening. *) |
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
734 |
|
58631 | 735 |
fun deepen ctxt f n = |
736 |
(trace_message ctxt (fn () => "Searching with depth limit " ^ string_of_int n); |
|
737 |
(f n handle Failure s => |
|
738 |
(trace_message ctxt (fn () => "failed with message: " ^ s); deepen ctxt f (n + 1)))); |
|
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
739 |
|
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
740 |
|
32645
1cc5b24f5a01
sos method generates and uses proof certificates
Philipp Meyer
parents:
32332
diff
changeset
|
741 |
(* Map back polynomials and their composites to a positivstellensatz. *) |
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
742 |
|
55508 | 743 |
fun cterm_of_sqterm (c, p) = RealArith.Product (RealArith.Rational_lt c, RealArith.Square p); |
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
744 |
|
55508 | 745 |
fun cterm_of_sos (pr,sqs) = |
746 |
if null sqs then pr |
|
747 |
else RealArith.Product (pr, foldr1 RealArith.Sum (map cterm_of_sqterm sqs)); |
|
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
748 |
|
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
749 |
(* Interface to HOL. *) |
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
750 |
local |
32828 | 751 |
open Conv |
59582 | 752 |
val concl = Thm.dest_arg o Thm.cprop_of |
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
753 |
in |
55508 | 754 |
(* FIXME: Replace tryfind by get_first !! *) |
32645
1cc5b24f5a01
sos method generates and uses proof certificates
Philipp Meyer
parents:
32332
diff
changeset
|
755 |
fun real_nonlinear_prover proof_method ctxt = |
32839 | 756 |
let |
55508 | 757 |
val {add = _, mul = _, neg = _, pow = _, sub = _, main = real_poly_conv} = |
758 |
Semiring_Normalizer.semiring_normalizers_ord_wrapper ctxt |
|
67271 | 759 |
(the (Semiring_Normalizer.match ctxt \<^cterm>\<open>(0::real) + 1\<close>)) |
67562 | 760 |
Thm.term_ord |
55508 | 761 |
fun mainf cert_choice translator (eqs, les, lts) = |
762 |
let |
|
763 |
val eq0 = map (poly_of_term o Thm.dest_arg1 o concl) eqs |
|
764 |
val le0 = map (poly_of_term o Thm.dest_arg o concl) les |
|
765 |
val lt0 = map (poly_of_term o Thm.dest_arg o concl) lts |
|
766 |
val eqp0 = map_index (fn (i, t) => (t,RealArith.Axiom_eq i)) eq0 |
|
767 |
val lep0 = map_index (fn (i, t) => (t,RealArith.Axiom_le i)) le0 |
|
768 |
val ltp0 = map_index (fn (i, t) => (t,RealArith.Axiom_lt i)) lt0 |
|
769 |
val (keq,eq) = List.partition (fn (p, _) => multidegree p = 0) eqp0 |
|
770 |
val (klep,lep) = List.partition (fn (p, _) => multidegree p = 0) lep0 |
|
771 |
val (kltp,ltp) = List.partition (fn (p, _) => multidegree p = 0) ltp0 |
|
772 |
fun trivial_axiom (p, ax) = |
|
773 |
(case ax of |
|
774 |
RealArith.Axiom_eq n => |
|
63205 | 775 |
if eval FuncUtil.Ctermfunc.empty p <> @0 then nth eqs n |
55508 | 776 |
else raise Failure "trivial_axiom: Not a trivial axiom" |
777 |
| RealArith.Axiom_le n => |
|
63205 | 778 |
if eval FuncUtil.Ctermfunc.empty p < @0 then nth les n |
55508 | 779 |
else raise Failure "trivial_axiom: Not a trivial axiom" |
780 |
| RealArith.Axiom_lt n => |
|
63205 | 781 |
if eval FuncUtil.Ctermfunc.empty p <= @0 then nth lts n |
55508 | 782 |
else raise Failure "trivial_axiom: Not a trivial axiom" |
783 |
| _ => error "trivial_axiom: Not a trivial axiom") |
|
784 |
in |
|
785 |
let val th = tryfind trivial_axiom (keq @ klep @ kltp) in |
|
786 |
(fconv_rule (arg_conv (arg1_conv (real_poly_conv ctxt)) |
|
787 |
then_conv Numeral_Simprocs.field_comp_conv ctxt) th, |
|
788 |
RealArith.Trivial) |
|
789 |
end handle Failure _ => |
|
790 |
let |
|
791 |
val proof = |
|
792 |
(case proof_method of |
|
793 |
Certificate certs => |
|
794 |
(* choose certificate *) |
|
795 |
let |
|
796 |
fun chose_cert [] (RealArith.Cert c) = c |
|
797 |
| chose_cert (RealArith.Left::s) (RealArith.Branch (l, _)) = chose_cert s l |
|
798 |
| chose_cert (RealArith.Right::s) (RealArith.Branch (_, r)) = chose_cert s r |
|
799 |
| chose_cert _ _ = error "certificate tree in invalid form" |
|
800 |
in |
|
801 |
chose_cert cert_choice certs |
|
802 |
end |
|
803 |
| Prover prover => |
|
804 |
(* call prover *) |
|
805 |
let |
|
63205 | 806 |
val pol = fold_rev poly_mul (map fst ltp) (poly_const @1) |
55508 | 807 |
val leq = lep @ ltp |
808 |
fun tryall d = |
|
809 |
let |
|
810 |
val e = multidegree pol |
|
811 |
val k = if e = 0 then 0 else d div e |
|
812 |
val eq' = map fst eq |
|
813 |
in |
|
814 |
tryfind (fn i => |
|
815 |
(d, i, real_positivnullstellensatz_general ctxt prover false d eq' leq |
|
816 |
(poly_neg(poly_pow pol i)))) |
|
817 |
(0 upto k) |
|
818 |
end |
|
58631 | 819 |
val (_,i,(cert_ideal,cert_cone)) = deepen ctxt tryall 0 |
55508 | 820 |
val proofs_ideal = |
821 |
map2 (fn q => fn (_,ax) => RealArith.Eqmul(q,ax)) cert_ideal eq |
|
822 |
val proofs_cone = map cterm_of_sos cert_cone |
|
823 |
val proof_ne = |
|
63205 | 824 |
if null ltp then RealArith.Rational_lt @1 |
55508 | 825 |
else |
826 |
let val p = foldr1 RealArith.Product (map snd ltp) in |
|
827 |
funpow i (fn q => RealArith.Product (p, q)) |
|
63205 | 828 |
(RealArith.Rational_lt @1) |
55508 | 829 |
end |
830 |
in |
|
831 |
foldr1 RealArith.Sum (proof_ne :: proofs_ideal @ proofs_cone) |
|
832 |
end) |
|
833 |
in |
|
834 |
(translator (eqs,les,lts) proof, RealArith.Cert proof) |
|
835 |
end |
|
836 |
end |
|
837 |
in mainf end |
|
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
838 |
end |
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
839 |
|
55508 | 840 |
(* FIXME : This is very bad!!!*) |
32839 | 841 |
fun subst_conv eqs t = |
55508 | 842 |
let |
843 |
val t' = fold (Thm.lambda o Thm.lhs_of) eqs t |
|
844 |
in |
|
60329
e85ed7a36b2f
eliminated odd C combinator -- Isabelle/ML usually has canonical argument order;
wenzelm
parents:
59582
diff
changeset
|
845 |
Conv.fconv_rule (Thm.beta_conversion true) |
e85ed7a36b2f
eliminated odd C combinator -- Isabelle/ML usually has canonical argument order;
wenzelm
parents:
59582
diff
changeset
|
846 |
(fold (fn a => fn b => Thm.combination b a) eqs (Thm.reflexive t')) |
55508 | 847 |
end |
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
848 |
|
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
849 |
(* A wrapper that tries to substitute away variables first. *) |
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
850 |
|
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
851 |
local |
55508 | 852 |
open Conv |
59582 | 853 |
val concl = Thm.dest_arg o Thm.cprop_of |
55508 | 854 |
val shuffle1 = |
855 |
fconv_rule (rewr_conv @{lemma "(a + x == y) == (x == y - (a::real))" |
|
856 |
by (atomize (full)) (simp add: field_simps)}) |
|
857 |
val shuffle2 = |
|
858 |
fconv_rule (rewr_conv @{lemma "(x + a == y) == (x == y - (a::real))" |
|
859 |
by (atomize (full)) (simp add: field_simps)}) |
|
860 |
fun substitutable_monomial fvs tm = |
|
59582 | 861 |
(case Thm.term_of tm of |
67271 | 862 |
Free (_, \<^typ>\<open>real\<close>) => |
63205 | 863 |
if not (member (op aconvc) fvs tm) then (@1, tm) |
55508 | 864 |
else raise Failure "substitutable_monomial" |
69064
5840724b1d71
Prefix form of infix with * on either side no longer needs special treatment
nipkow
parents:
67562
diff
changeset
|
865 |
| \<^term>\<open>(*) :: real \<Rightarrow> _\<close> $ _ $ (Free _) => |
55508 | 866 |
if RealArith.is_ratconst (Thm.dest_arg1 tm) andalso |
867 |
not (member (op aconvc) fvs (Thm.dest_arg tm)) |
|
868 |
then (RealArith.dest_ratconst (Thm.dest_arg1 tm), Thm.dest_arg tm) |
|
869 |
else raise Failure "substitutable_monomial" |
|
67399 | 870 |
| \<^term>\<open>(+) :: real \<Rightarrow> _\<close>$_$_ => |
60818 | 871 |
(substitutable_monomial (Drule.cterm_add_frees (Thm.dest_arg tm) fvs) (Thm.dest_arg1 tm) |
55508 | 872 |
handle Failure _ => |
60818 | 873 |
substitutable_monomial (Drule.cterm_add_frees (Thm.dest_arg1 tm) fvs) (Thm.dest_arg tm)) |
55508 | 874 |
| _ => raise Failure "substitutable_monomial") |
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
875 |
|
32839 | 876 |
fun isolate_variable v th = |
55508 | 877 |
let |
59582 | 878 |
val w = Thm.dest_arg1 (Thm.cprop_of th) |
55508 | 879 |
in |
880 |
if v aconvc w then th |
|
881 |
else |
|
59582 | 882 |
(case Thm.term_of w of |
67399 | 883 |
\<^term>\<open>(+) :: real \<Rightarrow> _\<close> $ _ $ _ => |
55508 | 884 |
if Thm.dest_arg1 w aconvc v then shuffle2 th |
885 |
else isolate_variable v (shuffle1 th) |
|
886 |
| _ => error "isolate variable : This should not happen?") |
|
32839 | 887 |
end |
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
888 |
in |
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
889 |
|
32268
d50f0cb67578
Functionality for sum of squares to call a remote csdp prover
Philipp Meyer
parents:
32150
diff
changeset
|
890 |
fun real_nonlinear_subst_prover prover ctxt = |
55508 | 891 |
let |
892 |
val {add = _, mul = real_poly_mul_conv, neg = _, pow = _, sub = _, main = real_poly_conv} = |
|
44453 | 893 |
Semiring_Normalizer.semiring_normalizers_ord_wrapper ctxt |
67271 | 894 |
(the (Semiring_Normalizer.match ctxt \<^cterm>\<open>(0::real) + 1\<close>)) |
67562 | 895 |
Thm.term_ord |
55508 | 896 |
|
897 |
fun make_substitution th = |
|
898 |
let |
|
899 |
val (c,v) = substitutable_monomial [] (Thm.dest_arg1(concl th)) |
|
900 |
val th1 = |
|
901 |
Drule.arg_cong_rule |
|
69064
5840724b1d71
Prefix form of infix with * on either side no longer needs special treatment
nipkow
parents:
67562
diff
changeset
|
902 |
(Thm.apply \<^cterm>\<open>(*) :: real \<Rightarrow> _\<close> (RealArith.cterm_of_rat (Rat.inv c))) |
55508 | 903 |
(mk_meta_eq th) |
904 |
val th2 = fconv_rule (binop_conv (real_poly_mul_conv ctxt)) th1 |
|
905 |
in fconv_rule (arg_conv (real_poly_conv ctxt)) (isolate_variable v th2) end |
|
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
906 |
|
55508 | 907 |
fun oprconv cv ct = |
908 |
let val g = Thm.dest_fun2 ct in |
|
67399 | 909 |
if g aconvc \<^cterm>\<open>(\<le>) :: real \<Rightarrow> _\<close> orelse g aconvc \<^cterm>\<open>(<) :: real \<Rightarrow> _\<close> |
55508 | 910 |
then arg_conv cv ct else arg1_conv cv ct |
911 |
end |
|
912 |
fun mainf cert_choice translator = |
|
913 |
let |
|
914 |
fun substfirst (eqs, les, lts) = |
|
915 |
(let |
|
916 |
val eth = tryfind make_substitution eqs |
|
917 |
val modify = |
|
918 |
fconv_rule (arg_conv (oprconv(subst_conv [eth] then_conv (real_poly_conv ctxt)))) |
|
919 |
in |
|
920 |
substfirst |
|
921 |
(filter_out |
|
67271 | 922 |
(fn t => (Thm.dest_arg1 o Thm.dest_arg o Thm.cprop_of) t aconvc \<^cterm>\<open>0::real\<close>) |
55508 | 923 |
(map modify eqs), |
924 |
map modify les, |
|
925 |
map modify lts) |
|
926 |
end handle Failure _ => |
|
927 |
real_nonlinear_prover prover ctxt cert_choice translator (rev eqs, rev les, rev lts)) |
|
928 |
in substfirst end |
|
929 |
in mainf end |
|
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
930 |
|
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
931 |
(* Overall function. *) |
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
932 |
|
32645
1cc5b24f5a01
sos method generates and uses proof certificates
Philipp Meyer
parents:
32332
diff
changeset
|
933 |
fun real_sos prover ctxt = |
32828 | 934 |
RealArith.gen_prover_real_arith ctxt (real_nonlinear_subst_prover prover ctxt) |
55508 | 935 |
|
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
936 |
end; |
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset
|
937 |
|
32839 | 938 |
val known_sos_constants = |
67399 | 939 |
[\<^term>\<open>(\<Longrightarrow>)\<close>, \<^term>\<open>Trueprop\<close>, |
67271 | 940 |
\<^term>\<open>HOL.False\<close>, \<^term>\<open>HOL.implies\<close>, \<^term>\<open>HOL.conj\<close>, \<^term>\<open>HOL.disj\<close>, |
67399 | 941 |
\<^term>\<open>Not\<close>, \<^term>\<open>(=) :: bool \<Rightarrow> _\<close>, |
67271 | 942 |
\<^term>\<open>All :: (real \<Rightarrow> _) \<Rightarrow> _\<close>, \<^term>\<open>Ex :: (real \<Rightarrow> _) \<Rightarrow> _\<close>, |
67399 | 943 |
\<^term>\<open>(=) :: real \<Rightarrow> _\<close>, \<^term>\<open>(<) :: real \<Rightarrow> _\<close>, |
944 |
\<^term>\<open>(\<le>) :: real \<Rightarrow> _\<close>, |
|
945 |
\<^term>\<open>(+) :: real \<Rightarrow> _\<close>, \<^term>\<open>(-) :: real \<Rightarrow> _\<close>, |
|
69064
5840724b1d71
Prefix form of infix with * on either side no longer needs special treatment
nipkow
parents:
67562
diff
changeset
|
946 |
\<^term>\<open>(*) :: real \<Rightarrow> _\<close>, \<^term>\<open>uminus :: real \<Rightarrow> _\<close>, |
67399 | 947 |
\<^term>\<open>(/) :: real \<Rightarrow> _\<close>, \<^term>\<open>inverse :: real \<Rightarrow> _\<close>, |
948 |
\<^term>\<open>(^) :: real \<Rightarrow> _\<close>, \<^term>\<open>abs :: real \<Rightarrow> _\<close>, |
|
67271 | 949 |
\<^term>\<open>min :: real \<Rightarrow> _\<close>, \<^term>\<open>max :: real \<Rightarrow> _\<close>, |
950 |
\<^term>\<open>0::real\<close>, \<^term>\<open>1::real\<close>, |
|
951 |
\<^term>\<open>numeral :: num \<Rightarrow> nat\<close>, |
|
952 |
\<^term>\<open>numeral :: num \<Rightarrow> real\<close>, |
|
953 |
\<^term>\<open>Num.Bit0\<close>, \<^term>\<open>Num.Bit1\<close>, \<^term>\<open>Num.One\<close>]; |
|
31512 | 954 |
|
32839 | 955 |
fun check_sos kcts ct = |
55508 | 956 |
let |
59582 | 957 |
val t = Thm.term_of ct |
55508 | 958 |
val _ = |
959 |
if not (null (Term.add_tfrees t []) andalso null (Term.add_tvars t [])) |
|
960 |
then error "SOS: not sos. Additional type varables" |
|
961 |
else () |
|
962 |
val fs = Term.add_frees t [] |
|
963 |
val _ = |
|
67271 | 964 |
if exists (fn ((_,T)) => not (T = \<^typ>\<open>real\<close>)) fs |
55508 | 965 |
then error "SOS: not sos. Variables with type not real" |
966 |
else () |
|
967 |
val vs = Term.add_vars t [] |
|
968 |
val _ = |
|
67271 | 969 |
if exists (fn ((_,T)) => not (T = \<^typ>\<open>real\<close>)) vs |
55508 | 970 |
then error "SOS: not sos. Variables with type not real" |
971 |
else () |
|
972 |
val ukcs = subtract (fn (t,p) => Const p aconv t) kcts (Term.add_consts t []) |
|
973 |
val _ = |
|
974 |
if null ukcs then () |
|
975 |
else error ("SOSO: Unknown constants in Subgoal:" ^ commas (map fst ukcs)) |
|
976 |
in () end |
|
31512 | 977 |
|
60752 | 978 |
fun core_sos_tac print_cert prover = SUBPROOF (fn {concl, context = ctxt, ...} => |
32831 | 979 |
let |
32838
d9dfd30af9ae
core_sos_tac: SUBPROOF body operates on subgoal 1;
wenzelm
parents:
32831
diff
changeset
|
980 |
val _ = check_sos known_sos_constants concl |
60752 | 981 |
val (th, certificates) = real_sos prover ctxt (Thm.dest_arg concl) |
32838
d9dfd30af9ae
core_sos_tac: SUBPROOF body operates on subgoal 1;
wenzelm
parents:
32831
diff
changeset
|
982 |
val _ = print_cert certificates |
60752 | 983 |
in resolve_tac ctxt [th] 1 end); |
31131 | 984 |
|
44453 | 985 |
fun default_SOME _ NONE v = SOME v |
986 |
| default_SOME _ (SOME v) _ = SOME v; |
|
31131 | 987 |
|
988 |
fun lift_SOME f NONE a = f a |
|
44453 | 989 |
| lift_SOME _ (SOME a) _ = SOME a; |
31131 | 990 |
|
991 |
||
992 |
local |
|
59582 | 993 |
val is_numeral = can (HOLogic.dest_number o Thm.term_of) |
31131 | 994 |
in |
55508 | 995 |
fun get_denom b ct = |
59582 | 996 |
(case Thm.term_of ct of |
67399 | 997 |
\<^term>\<open>(/) :: real \<Rightarrow> _\<close> $ _ $ _ => |
55508 | 998 |
if is_numeral (Thm.dest_arg ct) |
999 |
then get_denom b (Thm.dest_arg1 ct) |
|
1000 |
else default_SOME (get_denom b) (get_denom b (Thm.dest_arg ct)) (Thm.dest_arg ct, b) |
|
67399 | 1001 |
| \<^term>\<open>(<) :: real \<Rightarrow> _\<close> $ _ $ _ => |
55508 | 1002 |
lift_SOME (get_denom true) (get_denom true (Thm.dest_arg ct)) (Thm.dest_arg1 ct) |
67399 | 1003 |
| \<^term>\<open>(\<le>) :: real \<Rightarrow> _\<close> $ _ $ _ => |
55508 | 1004 |
lift_SOME (get_denom true) (get_denom true (Thm.dest_arg ct)) (Thm.dest_arg1 ct) |
1005 |
| _ $ _ => lift_SOME (get_denom b) (get_denom b (Thm.dest_fun ct)) (Thm.dest_arg ct) |
|
1006 |
| _ => NONE) |
|
31131 | 1007 |
end; |
1008 |
||
55508 | 1009 |
fun elim_one_denom_tac ctxt = CSUBGOAL (fn (P, i) => |
1010 |
(case get_denom false P of |
|
1011 |
NONE => no_tac |
|
1012 |
| SOME (d, ord) => |
|
1013 |
let |
|
1014 |
val simp_ctxt = |
|
1015 |
ctxt addsimps @{thms field_simps} |
|
60867 | 1016 |
addsimps [@{thm power_divide}] |
55508 | 1017 |
val th = |
60801 | 1018 |
Thm.instantiate' [] [SOME d, SOME (Thm.dest_arg P)] |
67091 | 1019 |
(if ord then @{lemma "(d=0 \<longrightarrow> P) \<and> (d>0 \<longrightarrow> P) \<and> (d<(0::real) \<longrightarrow> P) \<Longrightarrow> P" by auto} |
1020 |
else @{lemma "(d=0 \<longrightarrow> P) \<and> (d \<noteq> (0::real) \<longrightarrow> P) \<Longrightarrow> P" by blast}) |
|
60752 | 1021 |
in resolve_tac ctxt [th] i THEN Simplifier.asm_full_simp_tac simp_ctxt i end)); |
31131 | 1022 |
|
1023 |
fun elim_denom_tac ctxt i = REPEAT (elim_one_denom_tac ctxt i); |
|
1024 |
||
32949
aa6c470a962a
eliminated slightly odd get/set operations in favour of Unsynchronized.ref;
wenzelm
parents:
32839
diff
changeset
|
1025 |
fun sos_tac print_cert prover ctxt = |
56536 | 1026 |
(* The SOS prover breaks if mult_nonneg_nonneg is in the simpset *) |
57889 | 1027 |
let val ctxt' = Context_Position.set_visible false ctxt delsimps @{thms mult_nonneg_nonneg} |
56536 | 1028 |
in Object_Logic.full_atomize_tac ctxt' THEN' |
1029 |
elim_denom_tac ctxt' THEN' |
|
1030 |
core_sos_tac print_cert prover ctxt' |
|
1031 |
end; |
|
31131 | 1032 |
|
31512 | 1033 |
end; |