author  wenzelm 
Sat, 18 Jul 2015 20:47:08 +0200  
changeset 60752  b48830b670a1 
parent 60329  e85ed7a36b2f 
child 60801  7664e0916eec 
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 

2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

22 
val rat_0 = Rat.zero; 
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

23 
val rat_1 = Rat.one; 
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

24 
val rat_2 = Rat.two; 
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

25 
val rat_10 = Rat.rat_of_int 10; 
33029  26 
val max = Integer.max; 
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

27 

2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

28 
val denominator_rat = Rat.quotient_of_rat #> snd #> Rat.rat_of_int; 
55508  29 

32839  30 
fun int_of_rat a = 
55508  31 
(case Rat.quotient_of_rat a of 
32 
(i, 1) => i 

33 
 _ => error "int_of_rat: not an int"); 

34 

35 
fun lcm_rat x y = 

36 
Rat.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

37 

32839  38 
fun rat_pow r i = 
39 
let fun pow r i = 

40 
if i = 0 then rat_1 else 

31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

41 
let val d = pow r (i div 2) 
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

42 
in d */ d */ (if i mod 2 = 0 then rat_1 else r) 
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

43 
end 
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

44 
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

45 

32839  46 
fun round_rat r = 
55508  47 
let 
48 
val (a,b) = Rat.quotient_of_rat (Rat.abs r) 

49 
val d = a div b 

50 
val s = if r </ rat_0 then (Rat.neg o Rat.rat_of_int) else Rat.rat_of_int 

51 
val x2 = 2 * (a  (b * d)) 

52 
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

53 

2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

54 
val abs_rat = Rat.abs; 
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

55 
val pow2 = rat_pow rat_2; 
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

56 
val pow10 = rat_pow rat_10; 
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

57 

58631  58 

42616
92715b528e78
added Attrib.setup_config_XXX conveniences, with implicit setup of the background theory;
wenzelm
parents:
41490
diff
changeset

59 
val trace = Attrib.setup_config_bool @{binding sos_trace} (K false); 
58631  60 
val debug = Attrib.setup_config_bool @{binding sos_debug} (K false); 
61 

62 
fun trace_message ctxt msg = 

63 
if Config.get ctxt trace orelse Config.get ctxt debug then tracing (msg ()) else (); 

64 
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

65 

2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

66 
exception Sanity; 
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

67 

2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

68 
exception Unsolvable; 
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

69 

32332  70 
exception Failure of string; 
71 

32645
1cc5b24f5a01
sos method generates and uses proof certificates
Philipp Meyer
parents:
32332
diff
changeset

72 
datatype proof_method = 
1cc5b24f5a01
sos method generates and uses proof certificates
Philipp Meyer
parents:
32332
diff
changeset

73 
Certificate of RealArith.pss_tree 
1cc5b24f5a01
sos method generates and uses proof certificates
Philipp Meyer
parents:
32332
diff
changeset

74 
 Prover of (string > string) 
1cc5b24f5a01
sos method generates and uses proof certificates
Philipp Meyer
parents:
32332
diff
changeset

75 

31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

76 
(* 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

77 

2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

78 
local 
55508  79 

31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

80 
fun normalize y = 
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

81 
if abs_rat y </ (rat_1 // rat_10) then normalize (rat_10 */ y)  1 
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

82 
else if abs_rat y >=/ rat_1 then normalize (y // rat_10) + 1 
32839  83 
else 0 
55508  84 

85 
in 

86 

31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

87 
fun decimalize d x = 
55508  88 
if x =/ rat_0 then "0.0" 
89 
else 

90 
let 

91 
val y = Rat.abs x 

92 
val e = normalize y 

93 
val z = pow10(~ e) */ y +/ rat_1 

94 
val k = int_of_rat (round_rat(pow10 d */ z)) 

95 
in 

96 
(if x </ rat_0 then "0." else "0.") ^ 

97 
implode (tl (raw_explode(string_of_int k))) ^ 

98 
(if e = 0 then "" else "e" ^ string_of_int e) 

99 
end 

100 

31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

101 
end; 
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

102 

2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

103 
(* Iterations over numbers, and lists indexed by numbers. *) 
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

104 

2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

105 
fun itern k l f a = 
55508  106 
(case l of 
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

107 
[] => a 
55508  108 
 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

109 

2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

110 
fun iter (m,n) f a = 
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

111 
if n < m then a 
55508  112 
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

113 

2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

114 
(* The main types. *) 
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

115 

55508  116 
type vector = int * Rat.rat FuncUtil.Intfunc.table; 
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

117 

55508  118 
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

119 

55508  120 
fun iszero (_, r) = r =/ rat_0; 
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

121 

32839  122 

31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

123 
(* Vectors. Conventionally indexed 1..n. *) 
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

124 

55508  125 
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

126 

55508  127 
fun dim (v: vector) = fst v; 
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

128 

55508  129 
fun vector_cmul c (v: vector) = 
130 
let val n = dim v in 

131 
if c =/ rat_0 then vector_0 n 

39027  132 
else (n,FuncUtil.Intfunc.map (fn _ => fn x => c */ x) (snd v)) 
55508  133 
end; 
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

134 

2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

135 
fun vector_of_list l = 
55508  136 
let val n = length l in 
58635  137 
(n, fold_rev FuncUtil.Intfunc.update (1 upto n ~~ l) FuncUtil.Intfunc.empty): vector 
55508  138 
end; 
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

139 

2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

140 
(* Matrices; again rows and columns indexed from 1. *) 
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

141 

55508  142 
fun dimensions (m: matrix) = fst m; 
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

143 

55508  144 
fun row k (m: matrix) : vector = 
145 
let val (_, j) = dimensions m in 

146 
(j, 

147 
FuncUtil.Intpairfunc.fold (fn ((i, j), c) => fn a => 

148 
if i = k then FuncUtil.Intfunc.update (j, c) a else a) (snd m) FuncUtil.Intfunc.empty) 

149 
end; 

31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

150 

2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

151 
(* Monomials. *) 
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

152 

32828  153 
fun monomial_eval assig m = 
154 
FuncUtil.Ctermfunc.fold (fn (x, k) => fn a => a */ rat_pow (FuncUtil.Ctermfunc.apply assig x) k) 

55508  155 
m rat_1; 
156 

32829
671eb46eb0a3
tuned FuncFun and FuncUtil structure in positivstellensatz.ML
Philipp Meyer
parents:
32828
diff
changeset

157 
val monomial_1 = FuncUtil.Ctermfunc.empty; 
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

158 

32828  159 
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

160 

32828  161 
val monomial_mul = 
33002  162 
FuncUtil.Ctermfunc.combine Integer.add (K false); 
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

163 

32839  164 
fun monomial_multidegree m = 
55508  165 
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

166 

55508  167 
fun monomial_variables m = FuncUtil.Ctermfunc.dom m; 
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

168 

2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

169 
(* Polynomials. *) 
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

170 

32828  171 
fun eval assig p = 
172 
FuncUtil.Monomialfunc.fold (fn (m, c) => fn a => a +/ c */ monomial_eval assig m) p rat_0; 

31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

173 

32829
671eb46eb0a3
tuned FuncFun and FuncUtil structure in positivstellensatz.ML
Philipp Meyer
parents:
32828
diff
changeset

174 
val poly_0 = FuncUtil.Monomialfunc.empty; 
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

175 

32839  176 
fun poly_isconst p = 
55508  177 
FuncUtil.Monomialfunc.fold (fn (m, _) => fn a => FuncUtil.Ctermfunc.is_empty m andalso a) 
178 
p true; 

31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

179 

55508  180 
fun poly_var x = FuncUtil.Monomialfunc.onefunc (monomial_var x, rat_1); 
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

181 

2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

182 
fun poly_const c = 
55508  183 
if c =/ rat_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

184 

32828  185 
fun poly_cmul c p = 
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

186 
if c =/ rat_0 then poly_0 
39027  187 
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

188 

55508  189 
fun poly_neg p = FuncUtil.Monomialfunc.map (K Rat.neg) p; 
190 

31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

191 

32828  192 
fun poly_add p1 p2 = 
193 
FuncUtil.Monomialfunc.combine (curry op +/) (fn x => x =/ rat_0) p1 p2; 

31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

194 

2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

195 
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

196 

32828  197 
fun poly_cmmul (c,m) p = 
55508  198 
if c =/ rat_0 then poly_0 
199 
else 

200 
if FuncUtil.Ctermfunc.is_empty m 

201 
then FuncUtil.Monomialfunc.map (fn _ => fn d => c */ d) p 

202 
else 

203 
FuncUtil.Monomialfunc.fold (fn (m', d) => fn a => 

204 
(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

205 

32828  206 
fun poly_mul p1 p2 = 
207 
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

208 

2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

209 
fun poly_square p = poly_mul p p; 
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

210 

2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

211 
fun poly_pow p k = 
55508  212 
if k = 0 then poly_const rat_1 
213 
else if k = 1 then p 

214 
else 

215 
let val q = poly_square(poly_pow p (k div 2)) 

216 
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

217 

32828  218 
fun multidegree p = 
44453  219 
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

220 

32828  221 
fun poly_variables p = 
55508  222 
sort FuncUtil.cterm_ord 
223 
(FuncUtil.Monomialfunc.fold_rev 

224 
(fn (m, _) => union (is_equal o FuncUtil.cterm_ord) (monomial_variables m)) p []); 

31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

225 

2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

226 
(* Conversion from HOL term. *) 
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

227 

2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

228 
local 
55508  229 
val neg_tm = @{cterm "uminus :: real => _"} 
230 
val add_tm = @{cterm "op + :: real => _"} 

231 
val sub_tm = @{cterm "op  :: real => _"} 

232 
val mul_tm = @{cterm "op * :: real => _"} 

233 
val inv_tm = @{cterm "inverse :: real => _"} 

234 
val div_tm = @{cterm "op / :: real => _"} 

235 
val pow_tm = @{cterm "op ^ :: real => _"} 

236 
val zero_tm = @{cterm "0:: real"} 

59582  237 
val is_numeral = can (HOLogic.dest_number o Thm.term_of) 
55508  238 
fun poly_of_term tm = 
239 
if tm aconvc zero_tm then poly_0 

240 
else 

241 
if RealArith.is_ratconst tm 

242 
then poly_const(RealArith.dest_ratconst tm) 

243 
else 

244 
(let 

245 
val (lop, r) = Thm.dest_comb tm 

246 
in 

247 
if lop aconvc neg_tm then poly_neg(poly_of_term r) 

248 
else if lop aconvc inv_tm then 

249 
let val p = poly_of_term r in 

250 
if poly_isconst p 

251 
then poly_const(Rat.inv (eval FuncUtil.Ctermfunc.empty p)) 

252 
else error "poly_of_term: inverse of nonconstant polyomial" 

253 
end 

254 
else 

255 
(let 

256 
val (opr,l) = Thm.dest_comb lop 

257 
in 

258 
if opr aconvc pow_tm andalso is_numeral r 

59582  259 
then poly_pow (poly_of_term l) ((snd o HOLogic.dest_number o Thm.term_of) r) 
55508  260 
else if opr aconvc add_tm 
261 
then poly_add (poly_of_term l) (poly_of_term r) 

262 
else if opr aconvc sub_tm 

263 
then poly_sub (poly_of_term l) (poly_of_term r) 

264 
else if opr aconvc mul_tm 

265 
then poly_mul (poly_of_term l) (poly_of_term r) 

266 
else if opr aconvc div_tm 

267 
then 

268 
let 

32839  269 
val p = poly_of_term l 
270 
val q = poly_of_term r 

55508  271 
in 
272 
if poly_isconst q 

273 
then poly_cmul (Rat.inv (eval FuncUtil.Ctermfunc.empty q)) p 

274 
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

275 
end 
55508  276 
else poly_var tm 
277 
end handle CTERM ("dest_comb",_) => poly_var tm) 

278 
end handle CTERM ("dest_comb",_) => poly_var tm) 

31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

279 
in 
55508  280 
val poly_of_term = fn tm => 
59582  281 
if type_of (Thm.term_of tm) = @{typ real} 
55508  282 
then poly_of_term tm 
283 
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

284 
end; 
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

285 

2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

286 
(* String of vector (just a list of spaceseparated numbers). *) 
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

287 

55508  288 
fun sdpa_of_vector (v: vector) = 
289 
let 

290 
val n = dim v 

291 
val strs = 

292 
map (decimalize 20 o (fn i => FuncUtil.Intfunc.tryapplyd (snd v) i rat_0)) (1 upto n) 

293 
in space_implode " " strs ^ "\n" end; 

31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

294 

55508  295 
fun triple_int_ord ((a, b, c), (a', b', c')) = 
296 
prod_ord int_ord (prod_ord int_ord int_ord) ((a, (b, c)), (a', (b', c'))); 

297 
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

298 

2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

299 
fun index_char str chr pos = 
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

300 
if pos >= String.size str then ~1 
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

301 
else if String.sub(str,pos) = chr then pos 
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

302 
else index_char str chr (pos + 1); 
55508  303 

304 
fun rat_of_quotient (a,b) = 

305 
if b = 0 then rat_0 else Rat.rat_of_quotient (a, b); 

306 

32839  307 
fun rat_of_string s = 
55508  308 
let val n = index_char s #"/" 0 in 
309 
if n = ~1 then s > Int.fromString > the > Rat.rat_of_int 

310 
else 

311 
let 

312 
val SOME numer = Int.fromString(String.substring(s,0,n)) 

313 
val SOME den = Int.fromString (String.substring(s,n+1,String.size s  n  1)) 

314 
in rat_of_quotient(numer, den) end 

315 
end; 

31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

316 

53975
22ee3fb9d596
backout c6297fa1031a  strange parsers are required to make this work;
wenzelm
parents:
53972
diff
changeset

317 

55508  318 
fun isnum x = member (op =) ["0", "1", "2", "3", "4", "5", "6", "7", "8", "9"] x; 
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

319 

53975
22ee3fb9d596
backout c6297fa1031a  strange parsers are required to make this work;
wenzelm
parents:
53972
diff
changeset

320 
(* More parser basics. *) 
22ee3fb9d596
backout c6297fa1031a  strange parsers are required to make this work;
wenzelm
parents:
53972
diff
changeset

321 
(* FIXME improper use of parser combinators ahead *) 
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

322 

55508  323 
val numeral = Scan.one isnum 
324 
val decimalint = Scan.repeat1 numeral >> (rat_of_string o implode) 

325 
val decimalfrac = Scan.repeat1 numeral 

326 
>> (fn s => rat_of_string(implode s) // pow10 (length s)) 

327 
val decimalsig = 

328 
decimalint  Scan.option (Scan.$$ "."  decimalfrac) 

329 
>> (fn (h,NONE) => h  (h,SOME x) => h +/ x) 

330 
fun signed prs = 

331 
$$ ""  prs >> Rat.neg 

332 
 $$ "+"  prs 

333 
 prs; 

31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

334 

55508  335 
fun emptyin def xs = if null xs then (def, xs) else Scan.fail xs 
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

336 

55508  337 
val exponent = ($$ "e"  $$ "E")  signed decimalint; 
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

338 

55508  339 
val decimal = signed decimalsig  (emptyin rat_0 exponent) 
340 
>> (fn (h, x) => h */ pow10 (int_of_rat x)); 

31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

341 

55508  342 
fun mkparser p s = 
53975
22ee3fb9d596
backout c6297fa1031a  strange parsers are required to make this work;
wenzelm
parents:
53972
diff
changeset

343 
let val (x,rst) = p (raw_explode s) 
55508  344 
in if null rst then x else error "mkparser: unparsed input" end; 
53975
22ee3fb9d596
backout c6297fa1031a  strange parsers are required to make this work;
wenzelm
parents:
53972
diff
changeset

345 

32332  346 
(* Parse back csdp output. *) 
53975
22ee3fb9d596
backout c6297fa1031a  strange parsers are required to make this work;
wenzelm
parents:
53972
diff
changeset

347 
(* FIXME improper use of parser combinators ahead *) 
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

348 

55508  349 
fun ignore _ = ((),[]) 
350 
fun csdpoutput inp = 

351 
((decimal  Scan.repeat (Scan.$$ " "  Scan.option decimal) >> 

33067
a70d5baa53ee
use plain Scan.repeat (NB: Scan.bulk is for cascading sources  mostly interna use);
wenzelm
parents:
33042
diff
changeset

352 
(fn (h,to) => map_filter I ((SOME h)::to)))  ignore >> vector_of_list) inp 
55508  353 
val parse_csdpoutput = mkparser csdpoutput 
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

354 

2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

355 
(* 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

356 
(* 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

357 
(* 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

358 
(* are extreme numbers in the original problem. *) 
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

359 

55508  360 
(* Version for (int*int*int) keys *) 
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

361 
local 
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

362 
fun max_rat x y = if x </ y then y else x 
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

363 
fun common_denominator fld amat acc = 
55508  364 
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

365 
fun maximal_element fld amat acc = 
44453  366 
fld (fn (_,c) => fn maxa => max_rat maxa (abs_rat c)) amat acc 
55508  367 
fun float_of_rat x = 
368 
let val (a,b) = Rat.quotient_of_rat x 

369 
in Real.fromInt a / Real.fromInt b end; 

370 
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

371 
in 
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

372 

55508  373 
fun tri_scale_then solver (obj:vector) mats = 
374 
let 

375 
val cd1 = fold_rev (common_denominator Inttriplefunc.fold) mats (rat_1) 

376 
val cd2 = common_denominator FuncUtil.Intfunc.fold (snd obj) (rat_1) 

377 
val mats' = map (Inttriplefunc.map (fn _ => fn x => cd1 */ x)) mats 

378 
val obj' = vector_cmul cd2 obj 

379 
val max1 = fold_rev (maximal_element Inttriplefunc.fold) mats' (rat_0) 

380 
val max2 = maximal_element FuncUtil.Intfunc.fold (snd obj') (rat_0) 

381 
val scal1 = pow2 (20  int_of_float(Math.ln (float_of_rat max1) / Math.ln 2.0)) 

382 
val scal2 = pow2 (20  int_of_float(Math.ln (float_of_rat max2) / Math.ln 2.0)) 

383 
val mats'' = map (Inttriplefunc.map (fn _ => fn x => x */ scal1)) mats' 

384 
val obj'' = vector_cmul scal2 obj' 

385 
in solver obj'' mats'' end 

31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

386 
end; 
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

387 

2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

388 
(* Round a vector to "nice" rationals. *) 
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

389 

55508  390 
fun nice_rational n x = round_rat (n */ x) // n; 
32839  391 
fun nice_vector n ((d,v) : vector) = 
55508  392 
(d, FuncUtil.Intfunc.fold (fn (i,c) => fn a => 
393 
let val y = nice_rational n c in 

394 
if c =/ rat_0 then a 

395 
else FuncUtil.Intfunc.update (i,y) a 

396 
end) v FuncUtil.Intfunc.empty): vector 

31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

397 

2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

398 
fun dest_ord f x = is_equal (f x); 
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

399 

2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

400 
(* Stuff for "equations" ((int*int*int)>num functions). *) 
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

401 

2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

402 
fun tri_equation_cmul c eq = 
55508  403 
if c =/ rat_0 then Inttriplefunc.empty 
404 
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

405 

55508  406 
fun tri_equation_add eq1 eq2 = 
407 
Inttriplefunc.combine (curry op +/) (fn x => x =/ rat_0) eq1 eq2; 

31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

408 

2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

409 
fun tri_equation_eval assig eq = 
55508  410 
let 
411 
fun value v = Inttriplefunc.apply assig v 

412 
in Inttriplefunc.fold (fn (v, c) => fn a => a +/ value v */ c) eq rat_0 end; 

31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

413 

2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

414 
(* Eliminate all variables, in an essentially arbitrary order. *) 
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

415 

2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

416 
fun tri_eliminate_all_equations one = 
55508  417 
let 
418 
fun choose_variable eq = 

419 
let val (v,_) = Inttriplefunc.choose eq 

420 
in 

421 
if is_equal (triple_int_ord(v,one)) then 

422 
let 

423 
val eq' = Inttriplefunc.delete_safe v eq 

424 
in 

425 
if Inttriplefunc.is_empty eq' then error "choose_variable" 

426 
else fst (Inttriplefunc.choose eq') 

427 
end 

428 
else v 

31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

429 
end 
55508  430 

431 
fun eliminate dun eqs = 

432 
(case eqs of 

433 
[] => dun 

434 
 eq :: oeqs => 

435 
if Inttriplefunc.is_empty eq then eliminate dun oeqs 

436 
else 

437 
let 

438 
val v = choose_variable eq 

439 
val a = Inttriplefunc.apply eq v 

440 
val eq' = 

441 
tri_equation_cmul ((Rat.rat_of_int ~1) // a) (Inttriplefunc.delete_safe v eq) 

442 
fun elim e = 

443 
let val b = Inttriplefunc.tryapplyd e v rat_0 in 

444 
if b =/ rat_0 then e 

445 
else tri_equation_add e (tri_equation_cmul (Rat.neg b // a) eq) 

446 
end 

447 
in 

448 
eliminate (Inttriplefunc.update(v, eq') (Inttriplefunc.map (K elim) dun)) 

449 
(map elim oeqs) 

450 
end) 

451 
in 

452 
fn eqs => 

453 
let 

454 
val assig = eliminate Inttriplefunc.empty eqs 

455 
val vs = Inttriplefunc.fold (fn (_, f) => fn a => 

456 
remove (dest_ord triple_int_ord) one (Inttriplefunc.dom f) @ a) assig [] 

457 
in (distinct (dest_ord triple_int_ord) vs,assig) end 

458 
end; 

32839  459 

31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

460 
(* Multiply equationparametrized poly by regular poly and add accumulator. *) 
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

461 

2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

462 
fun tri_epoly_pmul p q acc = 
55508  463 
FuncUtil.Monomialfunc.fold (fn (m1, c) => fn a => 
464 
FuncUtil.Monomialfunc.fold (fn (m2, e) => fn b => 

465 
let 

466 
val m = monomial_mul m1 m2 

467 
val es = FuncUtil.Monomialfunc.tryapplyd b m Inttriplefunc.empty 

468 
in 

469 
FuncUtil.Monomialfunc.update (m,tri_equation_add (tri_equation_cmul c e) es) b 

470 
end) q a) p acc; 

31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

471 

2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

472 
(* 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

473 
(* 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

474 
(* 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

475 
(* vol 45, pp. 363374, 1978. *) 
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

476 
(* *) 
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

477 
(* 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

478 
(* 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

479 
(* quadratic form that will tend to display constants. *) 
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

480 

2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

481 
(* 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

482 

2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

483 
local 
55508  484 
fun diagonalize n i m = 
485 
if FuncUtil.Intpairfunc.is_empty (snd m) then [] 

31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

486 
else 
55508  487 
let 
488 
val a11 = FuncUtil.Intpairfunc.tryapplyd (snd m) (i,i) rat_0 

489 
in 

490 
if a11 </ rat_0 then raise Failure "diagonalize: not PSD" 

491 
else if a11 =/ rat_0 then 

492 
if FuncUtil.Intfunc.is_empty (snd (row i m)) 

493 
then diagonalize n (i + 1) m 

494 
else raise Failure "diagonalize: not PSD ___ " 

495 
else 

496 
let 

497 
val v = row i m 

498 
val v' = 

499 
(fst v, FuncUtil.Intfunc.fold (fn (i, c) => fn a => 

500 
let val y = c // a11 

501 
in if y = rat_0 then a else FuncUtil.Intfunc.update (i,y) a 

502 
end) (snd v) FuncUtil.Intfunc.empty) 

503 
fun upt0 x y a = 

504 
if y = rat_0 then a 

505 
else FuncUtil.Intpairfunc.update (x,y) a 

506 
val m' = 

507 
((n, n), 

508 
iter (i + 1, n) (fn j => 

509 
iter (i + 1, n) (fn k => 

510 
(upt0 (j, k) 

511 
(FuncUtil.Intpairfunc.tryapplyd (snd m) (j, k) rat_0 / 

512 
FuncUtil.Intfunc.tryapplyd (snd v) j rat_0 */ 

513 
FuncUtil.Intfunc.tryapplyd (snd v') k rat_0)))) 

514 
FuncUtil.Intpairfunc.empty) 

515 
in (a11, v') :: diagonalize n (i + 1) m' end 

516 
end 

31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

517 
in 
55508  518 
fun diag m = 
519 
let 

520 
val nn = dimensions m 

521 
val n = fst nn 

522 
in 

523 
if snd nn <> n then error "diagonalize: nonsquare matrix" 

524 
else diagonalize n 1 m 

525 
end 

31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

526 
end; 
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

527 

2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

528 
(* Enumeration of monomials with given multidegree bound. *) 
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

529 

32839  530 
fun enumerate_monomials d vars = 
55508  531 
if d < 0 then [] 
532 
else if d = 0 then [FuncUtil.Ctermfunc.empty] 

533 
else if null vars then [monomial_1] 

534 
else 

535 
let val alts = 

536 
map_range (fn k => 

537 
let 

538 
val oths = enumerate_monomials (d  k) (tl vars) 

539 
in map (fn ks => if k = 0 then ks else FuncUtil.Ctermfunc.update (hd vars, k) ks) oths end) 

540 
(d + 1) 

541 
in flat alts end; 

31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

542 

2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

543 
(* Enumerate products of distinct input polys with degree <= d. *) 
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

544 
(* We ignore any constant input polynomials. *) 
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

545 
(* 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

546 

2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

547 
fun enumerate_products d pols = 
55508  548 
if d = 0 then [(poly_const rat_1,RealArith.Rational_lt rat_1)] 
549 
else if d < 0 then [] 

550 
else 

551 
(case pols of 

552 
[] => [(poly_const rat_1, RealArith.Rational_lt rat_1)] 

553 
 (p, b) :: ps => 

554 
let val e = multidegree p in 

555 
if e = 0 then enumerate_products d ps 

556 
else 

557 
enumerate_products d ps @ 

558 
map (fn (q, c) => (poly_mul p q, RealArith.Product (b, c))) 

559 
(enumerate_products (d  e) ps) 

560 
end) 

31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

561 

2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

562 
(* 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

563 

2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

564 
fun epoly_of_poly p = 
55508  565 
FuncUtil.Monomialfunc.fold (fn (m, c) => fn a => 
566 
FuncUtil.Monomialfunc.update (m, Inttriplefunc.onefunc ((0, 0, 0), Rat.neg c)) a) 

567 
p FuncUtil.Monomialfunc.empty; 

31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

568 

2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

569 
(* String for block diagonal matrix numbered k. *) 
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

570 

2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

571 
fun sdpa_of_blockdiagonal k m = 
55508  572 
let 
573 
val pfx = string_of_int k ^" " 

574 
val ents = 

575 
Inttriplefunc.fold 

576 
(fn ((b, i, j), c) => fn a => if i > j then a else ((b, i, j), c) :: a) 

577 
m [] 

59058
a78612c67ec0
renamed "pairself" to "apply2", in accordance to @{apply 2};
wenzelm
parents:
58635
diff
changeset

578 
val entss = sort (triple_int_ord o apply2 fst) ents 
55508  579 
in 
580 
fold_rev (fn ((b,i,j),c) => fn a => 

581 
pfx ^ string_of_int b ^ " " ^ string_of_int i ^ " " ^ string_of_int j ^ 

582 
" " ^ decimalize 20 c ^ "\n" ^ a) entss "" 

583 
end; 

31119
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 
(* 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

586 

32268
d50f0cb67578
Functionality for sum of squares to call a remote csdp prover
Philipp Meyer
parents:
32150
diff
changeset

587 
fun sdpa_of_blockproblem nblocks blocksizes obj mats = 
55508  588 
let val m = length mats  1 
589 
in 

590 
string_of_int m ^ "\n" ^ 

591 
string_of_int nblocks ^ "\n" ^ 

592 
(space_implode " " (map string_of_int blocksizes)) ^ 

593 
"\n" ^ 

594 
sdpa_of_vector obj ^ 

58635  595 
fold_rev (fn (k, m) => fn a => sdpa_of_blockdiagonal (k  1) m ^ a) 
596 
(1 upto length mats ~~ mats) "" 

55508  597 
end; 
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

598 

32268
d50f0cb67578
Functionality for sum of squares to call a remote csdp prover
Philipp Meyer
parents:
32150
diff
changeset

599 
(* 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

600 

55508  601 
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

602 
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

603 

31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

604 
(* 3D versions of matrix operations to consider blocks separately. *) 
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

605 

2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

606 
val bmatrix_add = Inttriplefunc.combine (curry op +/) (fn x => x =/ rat_0); 
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

607 
fun bmatrix_cmul c bm = 
32829
671eb46eb0a3
tuned FuncFun and FuncUtil structure in positivstellensatz.ML
Philipp Meyer
parents:
32828
diff
changeset

608 
if c =/ rat_0 then Inttriplefunc.empty 
39027  609 
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

610 

2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

611 
val bmatrix_neg = bmatrix_cmul (Rat.rat_of_int ~1); 
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

612 

2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

613 
(* Smash a block matrix into components. *) 
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

614 

2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

615 
fun blocks blocksizes bm = 
55508  616 
map (fn (bs, b0) => 
617 
let 

618 
val m = 

619 
Inttriplefunc.fold 

620 
(fn ((b, i, j), c) => fn a => 

621 
if b = b0 then FuncUtil.Intpairfunc.update ((i, j), c) a else a) 

622 
bm FuncUtil.Intpairfunc.empty 

623 
val _ = FuncUtil.Intpairfunc.fold (fn ((i, j), _) => fn a => max a (max i j)) m 0 

624 
in (((bs, bs), m): matrix) end) 

625 
(blocksizes ~~ (1 upto length blocksizes)); 

31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

626 

2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

627 
(* FIXME : Get rid of this !!!*) 
32268
d50f0cb67578
Functionality for sum of squares to call a remote csdp prover
Philipp Meyer
parents:
32150
diff
changeset

628 
local 
44453  629 
fun tryfind_with msg _ [] = raise Failure msg 
630 
 tryfind_with _ f (x::xs) = (f x handle Failure s => tryfind_with s f xs); 

32839  631 
in 
32268
d50f0cb67578
Functionality for sum of squares to call a remote csdp prover
Philipp Meyer
parents:
32150
diff
changeset

632 
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

633 
end 
d50f0cb67578
Functionality for sum of squares to call a remote csdp prover
Philipp Meyer
parents:
32150
diff
changeset

634 

31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

635 
(* Positiv and Nullstellensatz. Flag "linf" forces a linear representation. *) 
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

636 

38805  637 
fun real_positivnullstellensatz_general ctxt prover linf d eqs leqs pol = 
32839  638 
let 
55508  639 
val vars = 
640 
fold_rev (union (op aconvc) o poly_variables) 

641 
(pol :: eqs @ map fst leqs) [] 

642 
val monoid = 

643 
if linf then 

644 
(poly_const rat_1,RealArith.Rational_lt rat_1):: 

645 
(filter (fn (p,_) => multidegree p <= d) leqs) 

646 
else enumerate_products d leqs 

647 
val nblocks = length monoid 

648 
fun mk_idmultiplier k p = 

649 
let 

650 
val e = d  multidegree p 

651 
val mons = enumerate_monomials e vars 

652 
val nons = mons ~~ (1 upto length mons) 

653 
in 

654 
(mons, 

655 
fold_rev (fn (m, n) => 

656 
FuncUtil.Monomialfunc.update (m, Inttriplefunc.onefunc ((~k, ~n, n), rat_1))) 

657 
nons FuncUtil.Monomialfunc.empty) 

658 
end 

31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

659 

55508  660 
fun mk_sqmultiplier k (p,_) = 
661 
let 

662 
val e = (d  multidegree p) div 2 

663 
val mons = enumerate_monomials e vars 

664 
val nons = mons ~~ (1 upto length mons) 

665 
in 

666 
(mons, 

667 
fold_rev (fn (m1, n1) => 

668 
fold_rev (fn (m2, n2) => fn a => 

669 
let val m = monomial_mul m1 m2 in 

670 
if n1 > n2 then a 

671 
else 

672 
let 

673 
val c = if n1 = n2 then rat_1 else rat_2 

674 
val e = FuncUtil.Monomialfunc.tryapplyd a m Inttriplefunc.empty 

675 
in 

676 
FuncUtil.Monomialfunc.update 

677 
(m, tri_equation_add (Inttriplefunc.onefunc ((k, n1, n2), c)) e) a 

678 
end 

679 
end) nons) nons FuncUtil.Monomialfunc.empty) 

680 
end 

31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

681 

55508  682 
val (sqmonlist,sqs) = split_list (map2 mk_sqmultiplier (1 upto length monoid) monoid) 
683 
val (_(*idmonlist*),ids) = split_list (map2 mk_idmultiplier (1 upto length eqs) eqs) 

684 
val blocksizes = map length sqmonlist 

685 
val bigsum = 

58635  686 
fold_rev (fn (p, q) => fn a => tri_epoly_pmul p q a) (eqs ~~ ids) 
687 
(fold_rev (fn ((p, _), s) => fn a => tri_epoly_pmul p s a) (monoid ~~ sqs) 

688 
(epoly_of_poly (poly_neg pol))) 

55508  689 
val eqns = FuncUtil.Monomialfunc.fold (fn (_, e) => fn a => e :: a) bigsum [] 
690 
val (pvs, assig) = tri_eliminate_all_equations (0, 0, 0) eqns 

691 
val qvars = (0, 0, 0) :: pvs 

692 
val allassig = 

693 
fold_rev (fn v => Inttriplefunc.update (v, (Inttriplefunc.onefunc (v, rat_1)))) pvs assig 

694 
fun mk_matrix v = 

695 
Inttriplefunc.fold (fn ((b, i, j), ass) => fn m => 

696 
if b < 0 then m 

697 
else 

698 
let val c = Inttriplefunc.tryapplyd ass v rat_0 in 

699 
if c = rat_0 then m 

700 
else Inttriplefunc.update ((b, j, i), c) (Inttriplefunc.update ((b, i, j), c) m) 

701 
end) 

702 
allassig Inttriplefunc.empty 

703 
val diagents = 

704 
Inttriplefunc.fold 

705 
(fn ((b, i, j), e) => fn a => if b > 0 andalso i = j then tri_equation_add e a else a) 

706 
allassig Inttriplefunc.empty 

31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

707 

55508  708 
val mats = map mk_matrix qvars 
709 
val obj = 

710 
(length pvs, 

711 
itern 1 pvs (fn v => fn i => 

712 
FuncUtil.Intfunc.updatep iszero (i,Inttriplefunc.tryapplyd diagents v rat_0)) 

713 
FuncUtil.Intfunc.empty) 

714 
val raw_vec = 

715 
if null pvs then vector_0 0 

716 
else tri_scale_then (run_blockproblem prover nblocks blocksizes) obj mats 

717 
fun int_element (_, v) i = FuncUtil.Intfunc.tryapplyd v i rat_0 

31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

718 

55508  719 
fun find_rounding d = 
720 
let 

721 
val _ = 

58631  722 
debug_message ctxt (fn () => "Trying rounding with limit "^Rat.string_of_rat d ^ "\n") 
55508  723 
val vec = nice_vector d raw_vec 
724 
val blockmat = 

725 
iter (1, dim vec) 

726 
(fn i => fn a => bmatrix_add (bmatrix_cmul (int_element vec i) (nth mats i)) a) 

727 
(bmatrix_neg (nth mats 0)) 

728 
val allmats = blocks blocksizes blockmat 

729 
in (vec, map diag allmats) end 

730 
val (vec, ratdias) = 

731 
if null pvs then find_rounding rat_1 

732 
else tryfind find_rounding (map Rat.rat_of_int (1 upto 31) @ map pow2 (5 upto 66)) 

733 
val newassigs = 

734 
fold_rev (fn k => Inttriplefunc.update (nth pvs (k  1), int_element vec k)) 

735 
(1 upto dim vec) (Inttriplefunc.onefunc ((0, 0, 0), Rat.rat_of_int ~1)) 

736 
val finalassigs = 

737 
Inttriplefunc.fold (fn (v, e) => fn a => 

738 
Inttriplefunc.update (v, tri_equation_eval newassigs e) a) allassig newassigs 

739 
fun poly_of_epoly p = 

740 
FuncUtil.Monomialfunc.fold (fn (v, e) => fn a => 

741 
FuncUtil.Monomialfunc.updatep iszero (v, tri_equation_eval finalassigs e) a) 

742 
p FuncUtil.Monomialfunc.empty 

743 
fun mk_sos mons = 

744 
let 

745 
fun mk_sq (c, m) = 

746 
(c, fold_rev (fn k => fn a => 

747 
FuncUtil.Monomialfunc.updatep iszero (nth mons (k  1), int_element m k) a) 

748 
(1 upto length mons) FuncUtil.Monomialfunc.empty) 

749 
in map mk_sq end 

750 
val sqs = map2 mk_sos sqmonlist ratdias 

751 
val cfs = map poly_of_epoly ids 

752 
val msq = filter (fn (_, b) => not (null b)) (map2 pair monoid sqs) 

753 
fun eval_sq sqs = fold_rev (fn (c, q) => poly_add (poly_cmul c (poly_mul q q))) sqs poly_0 

754 
val sanity = 

755 
fold_rev (fn ((p, _), s) => poly_add (poly_mul p (eval_sq s))) msq 

58635  756 
(fold_rev (fn (p, q) => poly_add (poly_mul p q)) (cfs ~~ eqs) (poly_neg pol)) 
55508  757 
in 
758 
if not(FuncUtil.Monomialfunc.is_empty sanity) then raise Sanity 

759 
else (cfs, map (fn (a, b) => (snd a, b)) msq) 

760 
end 

31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

761 

2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

762 

2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

763 
(* Iterative deepening. *) 
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

764 

58631  765 
fun deepen ctxt f n = 
766 
(trace_message ctxt (fn () => "Searching with depth limit " ^ string_of_int n); 

767 
(f n handle Failure s => 

768 
(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

769 

2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

770 

32645
1cc5b24f5a01
sos method generates and uses proof certificates
Philipp Meyer
parents:
32332
diff
changeset

771 
(* 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

772 

55508  773 
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

774 

55508  775 
fun cterm_of_sos (pr,sqs) = 
776 
if null sqs then pr 

777 
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

778 

2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

779 
(* Interface to HOL. *) 
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

780 
local 
32828  781 
open Conv 
59582  782 
val concl = Thm.dest_arg o Thm.cprop_of 
783 
fun simple_cterm_ord t u = Term_Ord.fast_term_ord (Thm.term_of t, Thm.term_of u) = LESS 

31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

784 
in 
55508  785 
(* FIXME: Replace tryfind by get_first !! *) 
32645
1cc5b24f5a01
sos method generates and uses proof certificates
Philipp Meyer
parents:
32332
diff
changeset

786 
fun real_nonlinear_prover proof_method ctxt = 
32839  787 
let 
55508  788 
val {add = _, mul = _, neg = _, pow = _, sub = _, main = real_poly_conv} = 
789 
Semiring_Normalizer.semiring_normalizers_ord_wrapper ctxt 

790 
(the (Semiring_Normalizer.match ctxt @{cterm "(0::real) + 1"})) 

791 
simple_cterm_ord 

792 
fun mainf cert_choice translator (eqs, les, lts) = 

793 
let 

794 
val eq0 = map (poly_of_term o Thm.dest_arg1 o concl) eqs 

795 
val le0 = map (poly_of_term o Thm.dest_arg o concl) les 

796 
val lt0 = map (poly_of_term o Thm.dest_arg o concl) lts 

797 
val eqp0 = map_index (fn (i, t) => (t,RealArith.Axiom_eq i)) eq0 

798 
val lep0 = map_index (fn (i, t) => (t,RealArith.Axiom_le i)) le0 

799 
val ltp0 = map_index (fn (i, t) => (t,RealArith.Axiom_lt i)) lt0 

800 
val (keq,eq) = List.partition (fn (p, _) => multidegree p = 0) eqp0 

801 
val (klep,lep) = List.partition (fn (p, _) => multidegree p = 0) lep0 

802 
val (kltp,ltp) = List.partition (fn (p, _) => multidegree p = 0) ltp0 

803 
fun trivial_axiom (p, ax) = 

804 
(case ax of 

805 
RealArith.Axiom_eq n => 

806 
if eval FuncUtil.Ctermfunc.empty p <>/ Rat.zero then nth eqs n 

807 
else raise Failure "trivial_axiom: Not a trivial axiom" 

808 
 RealArith.Axiom_le n => 

809 
if eval FuncUtil.Ctermfunc.empty p </ Rat.zero then nth les n 

810 
else raise Failure "trivial_axiom: Not a trivial axiom" 

811 
 RealArith.Axiom_lt n => 

812 
if eval FuncUtil.Ctermfunc.empty p <=/ Rat.zero then nth lts n 

813 
else raise Failure "trivial_axiom: Not a trivial axiom" 

814 
 _ => error "trivial_axiom: Not a trivial axiom") 

815 
in 

816 
let val th = tryfind trivial_axiom (keq @ klep @ kltp) in 

817 
(fconv_rule (arg_conv (arg1_conv (real_poly_conv ctxt)) 

818 
then_conv Numeral_Simprocs.field_comp_conv ctxt) th, 

819 
RealArith.Trivial) 

820 
end handle Failure _ => 

821 
let 

822 
val proof = 

823 
(case proof_method of 

824 
Certificate certs => 

825 
(* choose certificate *) 

826 
let 

827 
fun chose_cert [] (RealArith.Cert c) = c 

828 
 chose_cert (RealArith.Left::s) (RealArith.Branch (l, _)) = chose_cert s l 

829 
 chose_cert (RealArith.Right::s) (RealArith.Branch (_, r)) = chose_cert s r 

830 
 chose_cert _ _ = error "certificate tree in invalid form" 

831 
in 

832 
chose_cert cert_choice certs 

833 
end 

834 
 Prover prover => 

835 
(* call prover *) 

836 
let 

837 
val pol = fold_rev poly_mul (map fst ltp) (poly_const Rat.one) 

838 
val leq = lep @ ltp 

839 
fun tryall d = 

840 
let 

841 
val e = multidegree pol 

842 
val k = if e = 0 then 0 else d div e 

843 
val eq' = map fst eq 

844 
in 

845 
tryfind (fn i => 

846 
(d, i, real_positivnullstellensatz_general ctxt prover false d eq' leq 

847 
(poly_neg(poly_pow pol i)))) 

848 
(0 upto k) 

849 
end 

58631  850 
val (_,i,(cert_ideal,cert_cone)) = deepen ctxt tryall 0 
55508  851 
val proofs_ideal = 
852 
map2 (fn q => fn (_,ax) => RealArith.Eqmul(q,ax)) cert_ideal eq 

853 
val proofs_cone = map cterm_of_sos cert_cone 

854 
val proof_ne = 

855 
if null ltp then RealArith.Rational_lt Rat.one 

856 
else 

857 
let val p = foldr1 RealArith.Product (map snd ltp) in 

858 
funpow i (fn q => RealArith.Product (p, q)) 

859 
(RealArith.Rational_lt Rat.one) 

860 
end 

861 
in 

862 
foldr1 RealArith.Sum (proof_ne :: proofs_ideal @ proofs_cone) 

863 
end) 

864 
in 

865 
(translator (eqs,les,lts) proof, RealArith.Cert proof) 

866 
end 

867 
end 

868 
in mainf end 

31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

869 
end 
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

870 

55508  871 
(* FIXME : This is very bad!!!*) 
32839  872 
fun subst_conv eqs t = 
55508  873 
let 
874 
val t' = fold (Thm.lambda o Thm.lhs_of) eqs t 

875 
in 

60329
e85ed7a36b2f
eliminated odd C combinator  Isabelle/ML usually has canonical argument order;
wenzelm
parents:
59582
diff
changeset

876 
Conv.fconv_rule (Thm.beta_conversion true) 
e85ed7a36b2f
eliminated odd C combinator  Isabelle/ML usually has canonical argument order;
wenzelm
parents:
59582
diff
changeset

877 
(fold (fn a => fn b => Thm.combination b a) eqs (Thm.reflexive t')) 
55508  878 
end 
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

879 

2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

880 
(* A wrapper that tries to substitute away variables first. *) 
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

881 

2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

882 
local 
55508  883 
open Conv 
59582  884 
fun simple_cterm_ord t u = Term_Ord.fast_term_ord (Thm.term_of t, Thm.term_of u) = LESS 
885 
val concl = Thm.dest_arg o Thm.cprop_of 

55508  886 
val shuffle1 = 
887 
fconv_rule (rewr_conv @{lemma "(a + x == y) == (x == y  (a::real))" 

888 
by (atomize (full)) (simp add: field_simps)}) 

889 
val shuffle2 = 

890 
fconv_rule (rewr_conv @{lemma "(x + a == y) == (x == y  (a::real))" 

891 
by (atomize (full)) (simp add: field_simps)}) 

892 
fun substitutable_monomial fvs tm = 

59582  893 
(case Thm.term_of tm of 
55508  894 
Free (_, @{typ real}) => 
895 
if not (member (op aconvc) fvs tm) then (Rat.one, tm) 

896 
else raise Failure "substitutable_monomial" 

897 
 @{term "op * :: real => _"} $ _ $ (Free _) => 

898 
if RealArith.is_ratconst (Thm.dest_arg1 tm) andalso 

899 
not (member (op aconvc) fvs (Thm.dest_arg tm)) 

900 
then (RealArith.dest_ratconst (Thm.dest_arg1 tm), Thm.dest_arg tm) 

901 
else raise Failure "substitutable_monomial" 

902 
 @{term "op + :: real => _"}$_$_ => 

903 
(substitutable_monomial (Thm.add_cterm_frees (Thm.dest_arg tm) fvs) (Thm.dest_arg1 tm) 

904 
handle Failure _ => 

905 
substitutable_monomial (Thm.add_cterm_frees (Thm.dest_arg1 tm) fvs) (Thm.dest_arg tm)) 

906 
 _ => raise Failure "substitutable_monomial") 

31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

907 

32839  908 
fun isolate_variable v th = 
55508  909 
let 
59582  910 
val w = Thm.dest_arg1 (Thm.cprop_of th) 
55508  911 
in 
912 
if v aconvc w then th 

913 
else 

59582  914 
(case Thm.term_of w of 
55508  915 
@{term "op + :: real => _"} $ _ $ _ => 
916 
if Thm.dest_arg1 w aconvc v then shuffle2 th 

917 
else isolate_variable v (shuffle1 th) 

918 
 _ => error "isolate variable : This should not happen?") 

32839  919 
end 
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

920 
in 
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

921 

32268
d50f0cb67578
Functionality for sum of squares to call a remote csdp prover
Philipp Meyer
parents:
32150
diff
changeset

922 
fun real_nonlinear_subst_prover prover ctxt = 
55508  923 
let 
924 
val {add = _, mul = real_poly_mul_conv, neg = _, pow = _, sub = _, main = real_poly_conv} = 

44453  925 
Semiring_Normalizer.semiring_normalizers_ord_wrapper ctxt 
55508  926 
(the (Semiring_Normalizer.match ctxt @{cterm "(0::real) + 1"})) 
927 
simple_cterm_ord 

928 

929 
fun make_substitution th = 

930 
let 

931 
val (c,v) = substitutable_monomial [] (Thm.dest_arg1(concl th)) 

932 
val th1 = 

933 
Drule.arg_cong_rule 

934 
(Thm.apply @{cterm "op * :: real => _"} (RealArith.cterm_of_rat (Rat.inv c))) 

935 
(mk_meta_eq th) 

936 
val th2 = fconv_rule (binop_conv (real_poly_mul_conv ctxt)) th1 

937 
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

938 

55508  939 
fun oprconv cv ct = 
940 
let val g = Thm.dest_fun2 ct in 

941 
if g aconvc @{cterm "op <= :: real => _"} orelse g aconvc @{cterm "op < :: real => _"} 

942 
then arg_conv cv ct else arg1_conv cv ct 

943 
end 

944 
fun mainf cert_choice translator = 

945 
let 

946 
fun substfirst (eqs, les, lts) = 

947 
(let 

948 
val eth = tryfind make_substitution eqs 

949 
val modify = 

950 
fconv_rule (arg_conv (oprconv(subst_conv [eth] then_conv (real_poly_conv ctxt)))) 

951 
in 

952 
substfirst 

953 
(filter_out 

59582  954 
(fn t => (Thm.dest_arg1 o Thm.dest_arg o Thm.cprop_of) t aconvc @{cterm "0::real"}) 
55508  955 
(map modify eqs), 
956 
map modify les, 

957 
map modify lts) 

958 
end handle Failure _ => 

959 
real_nonlinear_prover prover ctxt cert_choice translator (rev eqs, rev les, rev lts)) 

960 
in substfirst end 

961 
in mainf end 

31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

962 

2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

963 
(* Overall function. *) 
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

964 

32645
1cc5b24f5a01
sos method generates and uses proof certificates
Philipp Meyer
parents:
32332
diff
changeset

965 
fun real_sos prover ctxt = 
32828  966 
RealArith.gen_prover_real_arith ctxt (real_nonlinear_subst_prover prover ctxt) 
55508  967 

31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

968 
end; 
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

969 

32839  970 
val known_sos_constants = 
971 
[@{term "op ==>"}, @{term "Trueprop"}, 

56625
54505623a8ef
sos accepts False, returns apply command
paulson <lp15@cam.ac.uk>
parents:
56536
diff
changeset

972 
@{term HOL.False}, @{term HOL.implies}, @{term HOL.conj}, @{term HOL.disj}, 
32839  973 
@{term "Not"}, @{term "op = :: bool => _"}, 
974 
@{term "All :: (real => _) => _"}, @{term "Ex :: (real => _) => _"}, 

975 
@{term "op = :: real => _"}, @{term "op < :: real => _"}, 

976 
@{term "op <= :: real => _"}, 

977 
@{term "op + :: real => _"}, @{term "op  :: real => _"}, 

978 
@{term "op * :: real => _"}, @{term "uminus :: real => _"}, 

31512  979 
@{term "op / :: real => _"}, @{term "inverse :: real => _"}, 
32839  980 
@{term "op ^ :: real => _"}, @{term "abs :: real => _"}, 
31512  981 
@{term "min :: real => _"}, @{term "max :: real => _"}, 
47108
2a1953f0d20d
merged fork with new numeral representation (see NEWS)
huffman
parents:
46497
diff
changeset

982 
@{term "0::real"}, @{term "1::real"}, 
2a1953f0d20d
merged fork with new numeral representation (see NEWS)
huffman
parents:
46497
diff
changeset

983 
@{term "numeral :: num => nat"}, 
2a1953f0d20d
merged fork with new numeral representation (see NEWS)
huffman
parents:
46497
diff
changeset

984 
@{term "numeral :: num => real"}, 
2a1953f0d20d
merged fork with new numeral representation (see NEWS)
huffman
parents:
46497
diff
changeset

985 
@{term "Num.Bit0"}, @{term "Num.Bit1"}, @{term "Num.One"}]; 
31512  986 

32839  987 
fun check_sos kcts ct = 
55508  988 
let 
59582  989 
val t = Thm.term_of ct 
55508  990 
val _ = 
991 
if not (null (Term.add_tfrees t []) andalso null (Term.add_tvars t [])) 

992 
then error "SOS: not sos. Additional type varables" 

993 
else () 

994 
val fs = Term.add_frees t [] 

995 
val _ = 

996 
if exists (fn ((_,T)) => not (T = @{typ "real"})) fs 

997 
then error "SOS: not sos. Variables with type not real" 

998 
else () 

999 
val vs = Term.add_vars t [] 

1000 
val _ = 

1001 
if exists (fn ((_,T)) => not (T = @{typ "real"})) vs 

1002 
then error "SOS: not sos. Variables with type not real" 

1003 
else () 

1004 
val ukcs = subtract (fn (t,p) => Const p aconv t) kcts (Term.add_consts t []) 

1005 
val _ = 

1006 
if null ukcs then () 

1007 
else error ("SOSO: Unknown constants in Subgoal:" ^ commas (map fst ukcs)) 

1008 
in () end 

31512  1009 

60752  1010 
fun core_sos_tac print_cert prover = SUBPROOF (fn {concl, context = ctxt, ...} => 
32831  1011 
let 
32838
d9dfd30af9ae
core_sos_tac: SUBPROOF body operates on subgoal 1;
wenzelm
parents:
32831
diff
changeset

1012 
val _ = check_sos known_sos_constants concl 
60752  1013 
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

1014 
val _ = print_cert certificates 
60752  1015 
in resolve_tac ctxt [th] 1 end); 
31131  1016 

44453  1017 
fun default_SOME _ NONE v = SOME v 
1018 
 default_SOME _ (SOME v) _ = SOME v; 

31131  1019 

1020 
fun lift_SOME f NONE a = f a 

44453  1021 
 lift_SOME _ (SOME a) _ = SOME a; 
31131  1022 

1023 

1024 
local 

59582  1025 
val is_numeral = can (HOLogic.dest_number o Thm.term_of) 
31131  1026 
in 
55508  1027 
fun get_denom b ct = 
59582  1028 
(case Thm.term_of ct of 
55508  1029 
@{term "op / :: real => _"} $ _ $ _ => 
1030 
if is_numeral (Thm.dest_arg ct) 

1031 
then get_denom b (Thm.dest_arg1 ct) 

1032 
else default_SOME (get_denom b) (get_denom b (Thm.dest_arg ct)) (Thm.dest_arg ct, b) 

1033 
 @{term "op < :: real => _"} $ _ $ _ => 

1034 
lift_SOME (get_denom true) (get_denom true (Thm.dest_arg ct)) (Thm.dest_arg1 ct) 

1035 
 @{term "op <= :: real => _"} $ _ $ _ => 

1036 
lift_SOME (get_denom true) (get_denom true (Thm.dest_arg ct)) (Thm.dest_arg1 ct) 

1037 
 _ $ _ => lift_SOME (get_denom b) (get_denom b (Thm.dest_fun ct)) (Thm.dest_arg ct) 

1038 
 _ => NONE) 

31131  1039 
end; 
1040 

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

1043 
NONE => no_tac 

1044 
 SOME (d, ord) => 

1045 
let 

1046 
val simp_ctxt = 

1047 
ctxt addsimps @{thms field_simps} 

1048 
addsimps [@{thm nonzero_power_divide}, @{thm power_divide}] 

1049 
val th = 

1050 
instantiate' [] [SOME d, SOME (Thm.dest_arg P)] 

1051 
(if ord then @{lemma "(d=0 > P) & (d>0 > P) & (d<(0::real) > P) ==> P" by auto} 

1052 
else @{lemma "(d=0 > P) & (d ~= (0::real) > P) ==> P" by blast}) 

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

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

1056 

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

1057 
fun sos_tac print_cert prover ctxt = 
56536  1058 
(* The SOS prover breaks if mult_nonneg_nonneg is in the simpset *) 
57889  1059 
let val ctxt' = Context_Position.set_visible false ctxt delsimps @{thms mult_nonneg_nonneg} 
56536  1060 
in Object_Logic.full_atomize_tac ctxt' THEN' 
1061 
elim_denom_tac ctxt' THEN' 

1062 
core_sos_tac print_cert prover ctxt' 

1063 
end; 

31131  1064 

31512  1065 
end; 