author  haftmann 
Fri, 14 Jun 2019 08:34:27 +0000  
changeset 70334  bba46912379e 
parent 69064  5840724b1d71 
permissions  rwrr 
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
adhoc 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
adhoc 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
adhoc 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
adhoc 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
adhoc 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 nonconstant 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 nonconstant 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 spaceseparated 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 floatingpoint, 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
adhoc 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
adhoc 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
adhoc 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
adhoc 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
adhoc 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 equationparametrized 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. 363374, 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
adhoc 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: nonsquare 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
adhoc 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 

70334  1009 
val sos_ss = @{context} 
1010 
> fold Simplifier.add_simp @{thms field_simps} 

1011 
> Simplifier.del_simp @{thm mult_nonneg_nonneg} 

1012 
> Simplifier.simpset_of; 

1013 

55508  1014 
fun elim_one_denom_tac ctxt = CSUBGOAL (fn (P, i) => 
1015 
(case get_denom false P of 

1016 
NONE => no_tac 

1017 
 SOME (d, ord) => 

1018 
let 

70334  1019 
val simp_ctxt = ctxt 
1020 
> Simplifier.put_simpset sos_ss 

55508  1021 
val th = 
60801  1022 
Thm.instantiate' [] [SOME d, SOME (Thm.dest_arg P)] 
67091  1023 
(if ord then @{lemma "(d=0 \<longrightarrow> P) \<and> (d>0 \<longrightarrow> P) \<and> (d<(0::real) \<longrightarrow> P) \<Longrightarrow> P" by auto} 
1024 
else @{lemma "(d=0 \<longrightarrow> P) \<and> (d \<noteq> (0::real) \<longrightarrow> P) \<Longrightarrow> P" by blast}) 

60752  1025 
in resolve_tac ctxt [th] i THEN Simplifier.asm_full_simp_tac simp_ctxt i end)); 
31131  1026 

1027 
fun elim_denom_tac ctxt i = REPEAT (elim_one_denom_tac ctxt i); 

1028 

32949
aa6c470a962a
eliminated slightly odd get/set operations in favour of Unsynchronized.ref;
wenzelm
parents:
32839
diff
changeset

1029 
fun sos_tac print_cert prover ctxt = 
70334  1030 
Object_Logic.full_atomize_tac ctxt THEN' 
1031 
elim_denom_tac ctxt THEN' 

1032 
core_sos_tac print_cert prover ctxt; 

31131  1033 

31512  1034 
end; 