author  wenzelm 
Sun, 10 Aug 2014 20:45:48 +0200  
changeset 57889  049e13f616d4 
parent 56625  54505623a8ef 
child 58631  41333b45bff9 
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 
32332  13 
exception Failure of string; 
32268
d50f0cb67578
Functionality for sum of squares to call a remote csdp prover
Philipp Meyer
parents:
32150
diff
changeset

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

15 

41474  16 
structure Sum_of_Squares: SUM_OF_SQUARES = 
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

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

18 

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

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

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

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

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

24 

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

25 
val denominator_rat = Rat.quotient_of_rat #> snd #> Rat.rat_of_int; 
55508  26 

32839  27 
fun int_of_rat a = 
55508  28 
(case Rat.quotient_of_rat a of 
29 
(i, 1) => i 

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

31 

32 
fun lcm_rat x y = 

33 
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

34 

32839  35 
fun rat_pow r i = 
36 
let fun pow r i = 

37 
if i = 0 then rat_1 else 

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

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

39 
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

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

41 
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

42 

32839  43 
fun round_rat r = 
55508  44 
let 
45 
val (a,b) = Rat.quotient_of_rat (Rat.abs r) 

46 
val d = a div b 

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

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

49 
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

50 

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

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

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

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

54 

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

55 
val trace = Attrib.setup_config_bool @{binding sos_trace} (K false); 
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

56 

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

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

58 

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

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

60 

32332  61 
exception Failure of string; 
62 

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

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

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

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

66 

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

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

68 

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

69 
local 
55508  70 

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

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

72 
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

73 
else if abs_rat y >=/ rat_1 then normalize (y // rat_10) + 1 
32839  74 
else 0 
55508  75 

76 
in 

77 

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

78 
fun decimalize d x = 
55508  79 
if x =/ rat_0 then "0.0" 
80 
else 

81 
let 

82 
val y = Rat.abs x 

83 
val e = normalize y 

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

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

86 
in 

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

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

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

90 
end 

91 

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

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

93 

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

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

95 

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

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

98 
[] => a 
55508  99 
 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

100 

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

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

102 
if n < m then a 
55508  103 
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

104 

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

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

106 

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

108 

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

110 

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

112 

32839  113 

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

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

115 

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

117 

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

119 

55508  120 
fun vector_cmul c (v: vector) = 
121 
let val n = dim v in 

122 
if c =/ rat_0 then vector_0 n 

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

125 

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

126 
fun vector_of_list l = 
55508  127 
let val n = length l in 
128 
(n, fold_rev2 (curry FuncUtil.Intfunc.update) (1 upto n) l FuncUtil.Intfunc.empty): vector 

129 
end; 

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

130 

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

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

132 

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

134 

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

137 
(j, 

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

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

140 
end; 

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

141 

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

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

143 

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

55508  146 
m rat_1; 
147 

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

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

149 

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

151 

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

154 

32839  155 
fun monomial_multidegree m = 
55508  156 
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

157 

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

159 

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

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

161 

32828  162 
fun eval assig p = 
163 
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

164 

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

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

166 

32839  167 
fun poly_isconst p = 
55508  168 
FuncUtil.Monomialfunc.fold (fn (m, _) => fn a => FuncUtil.Ctermfunc.is_empty m andalso a) 
169 
p true; 

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

170 

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

172 

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

173 
fun poly_const c = 
55508  174 
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

175 

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

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

179 

55508  180 
fun poly_neg p = FuncUtil.Monomialfunc.map (K Rat.neg) p; 
181 

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

182 

32828  183 
fun poly_add p1 p2 = 
184 
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

185 

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

186 
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

187 

32828  188 
fun poly_cmmul (c,m) p = 
55508  189 
if c =/ rat_0 then poly_0 
190 
else 

191 
if FuncUtil.Ctermfunc.is_empty m 

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

193 
else 

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

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

196 

32828  197 
fun poly_mul p1 p2 = 
198 
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

199 

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

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

201 

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

202 
fun poly_pow p k = 
55508  203 
if k = 0 then poly_const rat_1 
204 
else if k = 1 then p 

205 
else 

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

207 
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

208 

32828  209 
fun multidegree p = 
44453  210 
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

211 

32828  212 
fun poly_variables p = 
55508  213 
sort FuncUtil.cterm_ord 
214 
(FuncUtil.Monomialfunc.fold_rev 

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

216 

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

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

218 

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

219 
local 
55508  220 
val neg_tm = @{cterm "uminus :: real => _"} 
221 
val add_tm = @{cterm "op + :: real => _"} 

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

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

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

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

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

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

228 
val is_numeral = can (HOLogic.dest_number o term_of) 

229 
fun poly_of_term tm = 

230 
if tm aconvc zero_tm then poly_0 

231 
else 

232 
if RealArith.is_ratconst tm 

233 
then poly_const(RealArith.dest_ratconst tm) 

234 
else 

235 
(let 

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

237 
in 

238 
if lop aconvc neg_tm then poly_neg(poly_of_term r) 

239 
else if lop aconvc inv_tm then 

240 
let val p = poly_of_term r in 

241 
if poly_isconst p 

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

243 
else error "poly_of_term: inverse of nonconstant polyomial" 

244 
end 

245 
else 

246 
(let 

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

248 
in 

249 
if opr aconvc pow_tm andalso is_numeral r 

250 
then poly_pow (poly_of_term l) ((snd o HOLogic.dest_number o term_of) r) 

251 
else if opr aconvc add_tm 

252 
then poly_add (poly_of_term l) (poly_of_term r) 

253 
else if opr aconvc sub_tm 

254 
then poly_sub (poly_of_term l) (poly_of_term r) 

255 
else if opr aconvc mul_tm 

256 
then poly_mul (poly_of_term l) (poly_of_term r) 

257 
else if opr aconvc div_tm 

258 
then 

259 
let 

32839  260 
val p = poly_of_term l 
261 
val q = poly_of_term r 

55508  262 
in 
263 
if poly_isconst q 

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

265 
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

266 
end 
55508  267 
else poly_var tm 
268 
end handle CTERM ("dest_comb",_) => poly_var tm) 

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

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

270 
in 
55508  271 
val poly_of_term = fn tm => 
272 
if type_of (term_of tm) = @{typ real} 

273 
then poly_of_term tm 

274 
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

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

276 

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

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

278 

55508  279 
fun sdpa_of_vector (v: vector) = 
280 
let 

281 
val n = dim v 

282 
val strs = 

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

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

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

285 

55508  286 
fun triple_int_ord ((a, b, c), (a', b', c')) = 
287 
prod_ord int_ord (prod_ord int_ord int_ord) ((a, (b, c)), (a', (b', c'))); 

288 
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

289 

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

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

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

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

293 
else index_char str chr (pos + 1); 
55508  294 

295 
fun rat_of_quotient (a,b) = 

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

297 

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

301 
else 

302 
let 

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

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

305 
in rat_of_quotient(numer, den) end 

306 
end; 

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

307 

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

308 

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

310 

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

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

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

313 

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

316 
val decimalfrac = Scan.repeat1 numeral 

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

318 
val decimalsig = 

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

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

321 
fun signed prs = 

322 
$$ ""  prs >> Rat.neg 

323 
 $$ "+"  prs 

324 
 prs; 

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

325 

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

327 

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

329 

55508  330 
val decimal = signed decimalsig  (emptyin rat_0 exponent) 
331 
>> (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

332 

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

334 
let val (x,rst) = p (raw_explode s) 
55508  335 
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

336 

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

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

339 

55508  340 
fun ignore _ = ((),[]) 
341 
fun csdpoutput inp = 

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

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

345 

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

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

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

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

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

350 

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

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

353 
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

354 
fun common_denominator fld amat acc = 
55508  355 
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

356 
fun maximal_element fld amat acc = 
44453  357 
fld (fn (_,c) => fn maxa => max_rat maxa (abs_rat c)) amat acc 
55508  358 
fun float_of_rat x = 
359 
let val (a,b) = Rat.quotient_of_rat x 

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

361 
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

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

363 

55508  364 
fun tri_scale_then solver (obj:vector) mats = 
365 
let 

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

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

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

369 
val obj' = vector_cmul cd2 obj 

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

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

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

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

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

375 
val obj'' = vector_cmul scal2 obj' 

376 
in solver obj'' mats'' end 

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

377 
end; 
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 
(* Round a vector to "nice" rationals. *) 
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

380 

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

385 
if c =/ rat_0 then a 

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

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

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

388 

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

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

390 

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

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

392 

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

393 
fun tri_equation_cmul c eq = 
55508  394 
if c =/ rat_0 then Inttriplefunc.empty 
395 
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

396 

55508  397 
fun tri_equation_add eq1 eq2 = 
398 
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

399 

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

400 
fun tri_equation_eval assig eq = 
55508  401 
let 
402 
fun value v = Inttriplefunc.apply assig v 

403 
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

404 

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

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

406 

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

407 
fun tri_eliminate_all_equations one = 
55508  408 
let 
409 
fun choose_variable eq = 

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

411 
in 

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

413 
let 

414 
val eq' = Inttriplefunc.delete_safe v eq 

415 
in 

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

417 
else fst (Inttriplefunc.choose eq') 

418 
end 

419 
else v 

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

420 
end 
55508  421 

422 
fun eliminate dun eqs = 

423 
(case eqs of 

424 
[] => dun 

425 
 eq :: oeqs => 

426 
if Inttriplefunc.is_empty eq then eliminate dun oeqs 

427 
else 

428 
let 

429 
val v = choose_variable eq 

430 
val a = Inttriplefunc.apply eq v 

431 
val eq' = 

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

433 
fun elim e = 

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

435 
if b =/ rat_0 then e 

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

437 
end 

438 
in 

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

440 
(map elim oeqs) 

441 
end) 

442 
in 

443 
fn eqs => 

444 
let 

445 
val assig = eliminate Inttriplefunc.empty eqs 

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

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

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

449 
end; 

32839  450 

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

451 
(* Multiply equationparametrized poly by regular poly and add accumulator. *) 
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 
fun tri_epoly_pmul p q acc = 
55508  454 
FuncUtil.Monomialfunc.fold (fn (m1, c) => fn a => 
455 
FuncUtil.Monomialfunc.fold (fn (m2, e) => fn b => 

456 
let 

457 
val m = monomial_mul m1 m2 

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

459 
in 

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

461 
end) q a) p acc; 

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

462 

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

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

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

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

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

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

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

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

470 
(* quadratic form that will tend to display constants. *) 
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 
(* 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

473 

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

474 
local 
55508  475 
fun diagonalize n i m = 
476 
if FuncUtil.Intpairfunc.is_empty (snd m) then [] 

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

477 
else 
55508  478 
let 
479 
val a11 = FuncUtil.Intpairfunc.tryapplyd (snd m) (i,i) rat_0 

480 
in 

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

482 
else if a11 =/ rat_0 then 

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

484 
then diagonalize n (i + 1) m 

485 
else raise Failure "diagonalize: not PSD ___ " 

486 
else 

487 
let 

488 
val v = row i m 

489 
val v' = 

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

491 
let val y = c // a11 

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

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

494 
fun upt0 x y a = 

495 
if y = rat_0 then a 

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

497 
val m' = 

498 
((n, n), 

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

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

501 
(upt0 (j, k) 

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

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

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

505 
FuncUtil.Intpairfunc.empty) 

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

507 
end 

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

508 
in 
55508  509 
fun diag m = 
510 
let 

511 
val nn = dimensions m 

512 
val n = fst nn 

513 
in 

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

515 
else diagonalize n 1 m 

516 
end 

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

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

518 

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

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

520 

32839  521 
fun enumerate_monomials d vars = 
55508  522 
if d < 0 then [] 
523 
else if d = 0 then [FuncUtil.Ctermfunc.empty] 

524 
else if null vars then [monomial_1] 

525 
else 

526 
let val alts = 

527 
map_range (fn k => 

528 
let 

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

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

531 
(d + 1) 

532 
in flat alts end; 

31119
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 
(* Enumerate products of distinct input polys with degree <= d. *) 
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

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

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

537 

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

538 
fun enumerate_products d pols = 
55508  539 
if d = 0 then [(poly_const rat_1,RealArith.Rational_lt rat_1)] 
540 
else if d < 0 then [] 

541 
else 

542 
(case pols of 

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

544 
 (p, b) :: ps => 

545 
let val e = multidegree p in 

546 
if e = 0 then enumerate_products d ps 

547 
else 

548 
enumerate_products d ps @ 

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

550 
(enumerate_products (d  e) ps) 

551 
end) 

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

552 

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

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

554 

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

555 
fun epoly_of_poly p = 
55508  556 
FuncUtil.Monomialfunc.fold (fn (m, c) => fn a => 
557 
FuncUtil.Monomialfunc.update (m, Inttriplefunc.onefunc ((0, 0, 0), Rat.neg c)) a) 

558 
p FuncUtil.Monomialfunc.empty; 

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

559 

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

560 
(* String for block diagonal matrix numbered k. *) 
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 
fun sdpa_of_blockdiagonal k m = 
55508  563 
let 
564 
val pfx = string_of_int k ^" " 

565 
val ents = 

566 
Inttriplefunc.fold 

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

568 
m [] 

569 
val entss = sort (triple_int_ord o pairself fst) ents 

570 
in 

571 
fold_rev (fn ((b,i,j),c) => fn a => 

572 
pfx ^ string_of_int b ^ " " ^ string_of_int i ^ " " ^ string_of_int j ^ 

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

574 
end; 

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

575 

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

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

577 

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

578 
fun sdpa_of_blockproblem nblocks blocksizes obj mats = 
55508  579 
let val m = length mats  1 
580 
in 

581 
string_of_int m ^ "\n" ^ 

582 
string_of_int nblocks ^ "\n" ^ 

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

584 
"\n" ^ 

585 
sdpa_of_vector obj ^ 

586 
fold_rev2 (fn k => fn m => fn a => sdpa_of_blockdiagonal (k  1) m ^ a) 

587 
(1 upto length mats) mats "" 

588 
end; 

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

589 

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

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

591 

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

593 
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

594 

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

595 
(* 3D versions of matrix operations to consider blocks separately. *) 
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 
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

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

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

601 

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

602 
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

603 

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

604 
(* Smash a block matrix into components. *) 
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 
fun blocks blocksizes bm = 
55508  607 
map (fn (bs, b0) => 
608 
let 

609 
val m = 

610 
Inttriplefunc.fold 

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

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

613 
bm FuncUtil.Intpairfunc.empty 

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

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

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

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

617 

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

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

619 
local 
44453  620 
fun tryfind_with msg _ [] = raise Failure msg 
621 
 tryfind_with _ f (x::xs) = (f x handle Failure s => tryfind_with s f xs); 

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

623 
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

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

625 

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

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

627 

38805  628 
fun real_positivnullstellensatz_general ctxt prover linf d eqs leqs pol = 
32839  629 
let 
55508  630 
val vars = 
631 
fold_rev (union (op aconvc) o poly_variables) 

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

633 
val monoid = 

634 
if linf then 

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

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

637 
else enumerate_products d leqs 

638 
val nblocks = length monoid 

639 
fun mk_idmultiplier k p = 

640 
let 

641 
val e = d  multidegree p 

642 
val mons = enumerate_monomials e vars 

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

644 
in 

645 
(mons, 

646 
fold_rev (fn (m, n) => 

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

648 
nons FuncUtil.Monomialfunc.empty) 

649 
end 

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

650 

55508  651 
fun mk_sqmultiplier k (p,_) = 
652 
let 

653 
val e = (d  multidegree p) div 2 

654 
val mons = enumerate_monomials e vars 

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

656 
in 

657 
(mons, 

658 
fold_rev (fn (m1, n1) => 

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

660 
let val m = monomial_mul m1 m2 in 

661 
if n1 > n2 then a 

662 
else 

663 
let 

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

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

666 
in 

667 
FuncUtil.Monomialfunc.update 

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

669 
end 

670 
end) nons) nons FuncUtil.Monomialfunc.empty) 

671 
end 

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

672 

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

675 
val blocksizes = map length sqmonlist 

676 
val bigsum = 

677 
fold_rev2 (fn p => fn q => fn a => tri_epoly_pmul p q a) eqs ids 

678 
(fold_rev2 (fn (p,_) => fn s => fn a => tri_epoly_pmul p s a) monoid sqs 

679 
(epoly_of_poly(poly_neg pol))) 

680 
val eqns = FuncUtil.Monomialfunc.fold (fn (_, e) => fn a => e :: a) bigsum [] 

681 
val (pvs, assig) = tri_eliminate_all_equations (0, 0, 0) eqns 

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

683 
val allassig = 

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

685 
fun mk_matrix v = 

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

687 
if b < 0 then m 

688 
else 

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

690 
if c = rat_0 then m 

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

692 
end) 

693 
allassig Inttriplefunc.empty 

694 
val diagents = 

695 
Inttriplefunc.fold 

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

697 
allassig Inttriplefunc.empty 

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

698 

55508  699 
val mats = map mk_matrix qvars 
700 
val obj = 

701 
(length pvs, 

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

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

704 
FuncUtil.Intfunc.empty) 

705 
val raw_vec = 

706 
if null pvs then vector_0 0 

707 
else tri_scale_then (run_blockproblem prover nblocks blocksizes) obj mats 

708 
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

709 

55508  710 
fun find_rounding d = 
711 
let 

712 
val _ = 

713 
if Config.get ctxt trace 

714 
then writeln ("Trying rounding with limit "^Rat.string_of_rat d ^ "\n") 

715 
else () 

716 
val vec = nice_vector d raw_vec 

717 
val blockmat = 

718 
iter (1, dim vec) 

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

720 
(bmatrix_neg (nth mats 0)) 

721 
val allmats = blocks blocksizes blockmat 

722 
in (vec, map diag allmats) end 

723 
val (vec, ratdias) = 

724 
if null pvs then find_rounding rat_1 

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

726 
val newassigs = 

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

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

729 
val finalassigs = 

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

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

732 
fun poly_of_epoly p = 

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

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

735 
p FuncUtil.Monomialfunc.empty 

736 
fun mk_sos mons = 

737 
let 

738 
fun mk_sq (c, m) = 

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

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

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

742 
in map mk_sq end 

743 
val sqs = map2 mk_sos sqmonlist ratdias 

744 
val cfs = map poly_of_epoly ids 

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

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

747 
val sanity = 

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

749 
(fold_rev2 (fn p => fn q => poly_add (poly_mul p q)) cfs eqs (poly_neg pol)) 

750 
in 

751 
if not(FuncUtil.Monomialfunc.is_empty sanity) then raise Sanity 

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

753 
end 

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

754 

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

755 

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

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

757 

32839  758 
fun deepen f n = 
32949
aa6c470a962a
eliminated slightly odd get/set operations in favour of Unsynchronized.ref;
wenzelm
parents:
32839
diff
changeset

759 
(writeln ("Searching with depth limit " ^ string_of_int n); 
aa6c470a962a
eliminated slightly odd get/set operations in favour of Unsynchronized.ref;
wenzelm
parents:
32839
diff
changeset

760 
(f n handle Failure s => (writeln ("failed with message: " ^ s); deepen f (n + 1)))); 
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 

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

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

764 

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

766 

55508  767 
fun cterm_of_sos (pr,sqs) = 
768 
if null sqs then pr 

769 
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

770 

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

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

772 
local 
32828  773 
open Conv 
774 
val concl = Thm.dest_arg o cprop_of 

35408  775 
fun simple_cterm_ord t u = Term_Ord.fast_term_ord (term_of t, term_of u) = LESS 
31119
2532bb2d65c7
A decision method for universal multivariate real arithmetic with add
chaieb
parents:
diff
changeset

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

778 
fun real_nonlinear_prover proof_method ctxt = 
32839  779 
let 
55508  780 
val {add = _, mul = _, neg = _, pow = _, sub = _, main = real_poly_conv} = 
781 
Semiring_Normalizer.semiring_normalizers_ord_wrapper ctxt 

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

783 
simple_cterm_ord 

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

785 
let 

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

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

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

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

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

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

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

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

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

795 
fun trivial_axiom (p, ax) = 

796 
(case ax of 

797 
RealArith.Axiom_eq n => 

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

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

800 
 RealArith.Axiom_le n => 

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

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

803 
 RealArith.Axiom_lt n => 

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

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

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

807 
in 

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

809 
(fconv_rule (arg_conv (arg1_conv (real_poly_conv ctxt)) 

810 
then_conv Numeral_Simprocs.field_comp_conv ctxt) th, 

811 
RealArith.Trivial) 

812 
end handle Failure _ => 

813 
let 

814 
val proof = 

815 
(case proof_method of 

816 
Certificate certs => 

817 
(* choose certificate *) 

818 
let 

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

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

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

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

823 
in 

824 
chose_cert cert_choice certs 

825 
end 

826 
 Prover prover => 

827 
(* call prover *) 

828 
let 

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

830 
val leq = lep @ ltp 

831 
fun tryall d = 

832 
let 

833 
val e = multidegree pol 

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

835 
val eq' = map fst eq 

836 
in 

837 
tryfind (fn i => 

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

839 
(poly_neg(poly_pow pol i)))) 

840 
(0 upto k) 

841 
end 

842 
val (_,i,(cert_ideal,cert_cone)) = deepen tryall 0 

843 
val proofs_ideal = 

844 
map2 (fn q => fn (_,ax) => RealArith.Eqmul(q,ax)) cert_ideal eq 

845 
val proofs_cone = map cterm_of_sos cert_cone 

846 
val proof_ne = 

847 
if null ltp then RealArith.Rational_lt Rat.one 

848 
else 

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

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

851 
(RealArith.Rational_lt Rat.one) 

852 
end 

853 
in 

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

855 
end) 

856 
in 

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

858 
end 

859 
end 

860 
in mainf end 

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

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

862 

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

863 
fun C f x y = f y x; 
55508  864 
(* FIXME : This is very bad!!!*) 
32839  865 
fun subst_conv eqs t = 
55508  866 
let 
867 
val t' = fold (Thm.lambda o Thm.lhs_of) eqs t 

868 
in 

869 
Conv.fconv_rule (Thm.beta_conversion true) (fold (C Thm.combination) eqs (Thm.reflexive t')) 

870 
end 

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

871 

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

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

873 

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

874 
local 
55508  875 
open Conv 
35408  876 
fun simple_cterm_ord t u = Term_Ord.fast_term_ord (term_of t, term_of u) = LESS 
55508  877 
val concl = Thm.dest_arg o cprop_of 
878 
val shuffle1 = 

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

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

881 
val shuffle2 = 

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

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

884 
fun substitutable_monomial fvs tm = 

885 
(case term_of tm of 

886 
Free (_, @{typ real}) => 

887 
if not (member (op aconvc) fvs tm) then (Rat.one, tm) 

888 
else raise Failure "substitutable_monomial" 

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

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

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

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

893 
else raise Failure "substitutable_monomial" 

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

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

896 
handle Failure _ => 

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

898 
 _ => raise Failure "substitutable_monomial") 

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

899 

32839  900 
fun isolate_variable v th = 
55508  901 
let 
902 
val w = Thm.dest_arg1 (cprop_of th) 

903 
in 

904 
if v aconvc w then th 

905 
else 

906 
(case term_of w of 

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

908 
if Thm.dest_arg1 w aconvc v then shuffle2 th 

909 
else isolate_variable v (shuffle1 th) 

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

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

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

913 

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

914 
fun real_nonlinear_subst_prover prover ctxt = 
55508  915 
let 
916 
val {add = _, mul = real_poly_mul_conv, neg = _, pow = _, sub = _, main = real_poly_conv} = 

44453  917 
Semiring_Normalizer.semiring_normalizers_ord_wrapper ctxt 
55508  918 
(the (Semiring_Normalizer.match ctxt @{cterm "(0::real) + 1"})) 
919 
simple_cterm_ord 

920 

921 
fun make_substitution th = 

922 
let 

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

924 
val th1 = 

925 
Drule.arg_cong_rule 

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

927 
(mk_meta_eq th) 

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

929 
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

930 

55508  931 
fun oprconv cv ct = 
932 
let val g = Thm.dest_fun2 ct in 

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

934 
then arg_conv cv ct else arg1_conv cv ct 

935 
end 

936 
fun mainf cert_choice translator = 

937 
let 

938 
fun substfirst (eqs, les, lts) = 

939 
(let 

940 
val eth = tryfind make_substitution eqs 

941 
val modify = 

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

943 
in 

944 
substfirst 

945 
(filter_out 

946 
(fn t => (Thm.dest_arg1 o Thm.dest_arg o cprop_of) t aconvc @{cterm "0::real"}) 

947 
(map modify eqs), 

948 
map modify les, 

949 
map modify lts) 

950 
end handle Failure _ => 

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

952 
in substfirst end 

953 
in mainf end 

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

954 

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

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

956 

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

957 
fun real_sos prover ctxt = 
32828  958 
RealArith.gen_prover_real_arith ctxt (real_nonlinear_subst_prover prover ctxt) 
55508  959 

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

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

961 

32839  962 
val known_sos_constants = 
963 
[@{term "op ==>"}, @{term "Trueprop"}, 

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

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

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

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

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

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

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

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

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

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

977 
@{term "Num.Bit0"}, @{term "Num.Bit1"}, @{term "Num.One"}]; 
31512  978 

32839  979 
fun check_sos kcts ct = 
55508  980 
let 
981 
val t = term_of ct 

982 
val _ = 

983 
if not (null (Term.add_tfrees t []) andalso null (Term.add_tvars t [])) 

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

985 
else () 

986 
val fs = Term.add_frees t [] 

987 
val _ = 

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

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

990 
else () 

991 
val vs = Term.add_vars t [] 

992 
val _ = 

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

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

995 
else () 

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

997 
val _ = 

998 
if null ukcs then () 

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

1000 
in () end 

31512  1001 

32838
d9dfd30af9ae
core_sos_tac: SUBPROOF body operates on subgoal 1;
wenzelm
parents:
32831
diff
changeset

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

1004 
val _ = check_sos known_sos_constants concl 
d9dfd30af9ae
core_sos_tac: SUBPROOF body operates on subgoal 1;
wenzelm
parents:
32831
diff
changeset

1005 
val (ths, certificates) = real_sos prover context (Thm.dest_arg concl) 
d9dfd30af9ae
core_sos_tac: SUBPROOF body operates on subgoal 1;
wenzelm
parents:
32831
diff
changeset

1006 
val _ = print_cert certificates 
55508  1007 
in rtac ths 1 end); 
31131  1008 

44453  1009 
fun default_SOME _ NONE v = SOME v 
1010 
 default_SOME _ (SOME v) _ = SOME v; 

31131  1011 

1012 
fun lift_SOME f NONE a = f a 

44453  1013 
 lift_SOME _ (SOME a) _ = SOME a; 
31131  1014 

1015 

1016 
local 

55508  1017 
val is_numeral = can (HOLogic.dest_number o term_of) 
31131  1018 
in 
55508  1019 
fun get_denom b ct = 
1020 
(case term_of ct of 

1021 
@{term "op / :: real => _"} $ _ $ _ => 

1022 
if is_numeral (Thm.dest_arg ct) 

1023 
then get_denom b (Thm.dest_arg1 ct) 

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

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

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

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

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

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

1030 
 _ => NONE) 

31131  1031 
end; 
1032 

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

1035 
NONE => no_tac 

1036 
 SOME (d, ord) => 

1037 
let 

1038 
val simp_ctxt = 

1039 
ctxt addsimps @{thms field_simps} 

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

1041 
val th = 

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

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

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

1045 
in rtac th i THEN Simplifier.asm_full_simp_tac simp_ctxt i end)); 

31131  1046 

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

1048 

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

1049 
fun sos_tac print_cert prover ctxt = 
56536  1050 
(* The SOS prover breaks if mult_nonneg_nonneg is in the simpset *) 
57889  1051 
let val ctxt' = Context_Position.set_visible false ctxt delsimps @{thms mult_nonneg_nonneg} 
56536  1052 
in Object_Logic.full_atomize_tac ctxt' THEN' 
1053 
elim_denom_tac ctxt' THEN' 

1054 
core_sos_tac print_cert prover ctxt' 

1055 
end; 

31131  1056 

31512  1057 
end; 