renamed Thm.capply to Thm.apply, and Thm.cabs to Thm.lambda in conformance with similar operations in structure Term and Logic;
1 (* Title: HOL/Library/Sum_of_Squares/sum_of_squares.ML
2 Author: Amine Chaieb, University of Cambridge
3 Author: Philipp Meyer, TU Muenchen
5 A tactic for proving nonlinear inequalities.
8 signature SUM_OF_SQUARES =
10 datatype proof_method = Certificate of RealArith.pss_tree | Prover of string -> string
11 val sos_tac: (RealArith.pss_tree -> unit) -> proof_method -> Proof.context -> int -> tactic
12 val trace: bool Config.T
13 exception Failure of string;
16 structure Sum_of_Squares: SUM_OF_SQUARES =
22 val rat_10 = Rat.rat_of_int 10;
23 val max = Integer.max;
25 val denominator_rat = Rat.quotient_of_rat #> snd #> Rat.rat_of_int;
27 case Rat.quotient_of_rat a of (i,1) => i | _ => error "int_of_rat: not an int";
28 fun lcm_rat x y = Rat.rat_of_int (Integer.lcm (int_of_rat x) (int_of_rat y));
32 if i = 0 then rat_1 else
33 let val d = pow r (i div 2)
34 in d */ d */ (if i mod 2 = 0 then rat_1 else r)
36 in if i < 0 then pow (Rat.inv r) (~ i) else pow r i end;
39 let val (a,b) = Rat.quotient_of_rat (Rat.abs r)
41 val s = if r </ rat_0 then (Rat.neg o Rat.rat_of_int) else Rat.rat_of_int
42 val x2 = 2 * (a - (b * d))
43 in s (if x2 >= b then d + 1 else d) end
45 val abs_rat = Rat.abs;
46 val pow2 = rat_pow rat_2;
47 val pow10 = rat_pow rat_10;
49 val trace = Attrib.setup_config_bool @{binding sos_trace} (K false);
55 exception Failure of string;
57 datatype proof_method =
58 Certificate of RealArith.pss_tree
59 | Prover of (string -> string)
61 (* Turn a rational into a decimal string with d sig digits. *)
65 if abs_rat y </ (rat_1 // rat_10) then normalize (rat_10 */ y) - 1
66 else if abs_rat y >=/ rat_1 then normalize (y // rat_10) + 1
70 if x =/ rat_0 then "0.0" else
74 val z = pow10(~ e) */ y +/ rat_1
75 val k = int_of_rat (round_rat(pow10 d */ z))
76 in (if x </ rat_0 then "-0." else "0.") ^
77 implode(tl(raw_explode(string_of_int k))) ^
78 (if e = 0 then "" else "e"^string_of_int e)
82 (* Iterations over numbers, and lists indexed by numbers. *)
87 | h::t => itern (k + 1) t f (f h k a);
91 else iter (m+1,n) f (f m a);
95 type vector = int* Rat.rat FuncUtil.Intfunc.table;
97 type matrix = (int*int)*(Rat.rat FuncUtil.Intpairfunc.table);
99 fun iszero (_,r) = r =/ rat_0;
102 (* Vectors. Conventionally indexed 1..n. *)
104 fun vector_0 n = (n,FuncUtil.Intfunc.empty):vector;
106 fun dim (v:vector) = fst v;
108 fun vector_cmul c (v:vector) =
110 in if c =/ rat_0 then vector_0 n
111 else (n,FuncUtil.Intfunc.map (fn _ => fn x => c */ x) (snd v))
114 fun vector_of_list l =
116 in (n,fold_rev2 (curry FuncUtil.Intfunc.update) (1 upto n) l FuncUtil.Intfunc.empty) :vector
119 (* Matrices; again rows and columns indexed from 1. *)
121 fun dimensions (m:matrix) = fst m;
123 fun row k (m:matrix) =
124 let val (_,j) = dimensions m
126 FuncUtil.Intpairfunc.fold (fn ((i,j), c) => fn a => if i = k then FuncUtil.Intfunc.update (j,c) a else a) (snd m) FuncUtil.Intfunc.empty ) : vector
131 fun monomial_eval assig m =
132 FuncUtil.Ctermfunc.fold (fn (x, k) => fn a => a */ rat_pow (FuncUtil.Ctermfunc.apply assig x) k)
134 val monomial_1 = FuncUtil.Ctermfunc.empty;
136 fun monomial_var x = FuncUtil.Ctermfunc.onefunc (x, 1);
139 FuncUtil.Ctermfunc.combine Integer.add (K false);
141 fun monomial_multidegree m =
142 FuncUtil.Ctermfunc.fold (fn (_, k) => fn a => k + a) m 0;;
144 fun monomial_variables m = FuncUtil.Ctermfunc.dom m;;
149 FuncUtil.Monomialfunc.fold (fn (m, c) => fn a => a +/ c */ monomial_eval assig m) p rat_0;
151 val poly_0 = FuncUtil.Monomialfunc.empty;
154 FuncUtil.Monomialfunc.fold (fn (m, _) => fn a => FuncUtil.Ctermfunc.is_empty m andalso a) p true;
156 fun poly_var x = FuncUtil.Monomialfunc.onefunc (monomial_var x,rat_1);
159 if c =/ rat_0 then poly_0 else FuncUtil.Monomialfunc.onefunc(monomial_1, c);
162 if c =/ rat_0 then poly_0
163 else FuncUtil.Monomialfunc.map (fn _ => fn x => c */ x) p;
165 fun poly_neg p = FuncUtil.Monomialfunc.map (K Rat.neg) p;;
168 FuncUtil.Monomialfunc.combine (curry op +/) (fn x => x =/ rat_0) p1 p2;
170 fun poly_sub p1 p2 = poly_add p1 (poly_neg p2);
172 fun poly_cmmul (c,m) p =
173 if c =/ rat_0 then poly_0
174 else if FuncUtil.Ctermfunc.is_empty m
175 then FuncUtil.Monomialfunc.map (fn _ => fn d => c */ d) p
176 else FuncUtil.Monomialfunc.fold (fn (m', d) => fn a => (FuncUtil.Monomialfunc.update (monomial_mul m m', c */ d) a)) p poly_0;
179 FuncUtil.Monomialfunc.fold (fn (m, c) => fn a => poly_add (poly_cmmul (c,m) p2) a) p1 poly_0;
181 fun poly_square p = poly_mul p p;
184 if k = 0 then poly_const rat_1
186 else let val q = poly_square(poly_pow p (k div 2)) in
187 if k mod 2 = 1 then poly_mul p q else q end;
190 FuncUtil.Monomialfunc.fold (fn (m, _) => fn a => max (monomial_multidegree m) a) p 0;
192 fun poly_variables p =
193 sort FuncUtil.cterm_ord (FuncUtil.Monomialfunc.fold_rev (fn (m, _) => union (is_equal o FuncUtil.cterm_ord) (monomial_variables m)) p []);;
195 (* Conversion from HOL term. *)
198 val neg_tm = @{cterm "uminus :: real => _"}
199 val add_tm = @{cterm "op + :: real => _"}
200 val sub_tm = @{cterm "op - :: real => _"}
201 val mul_tm = @{cterm "op * :: real => _"}
202 val inv_tm = @{cterm "inverse :: real => _"}
203 val div_tm = @{cterm "op / :: real => _"}
204 val pow_tm = @{cterm "op ^ :: real => _"}
205 val zero_tm = @{cterm "0:: real"}
206 val is_numeral = can (HOLogic.dest_number o term_of)
207 fun poly_of_term tm =
208 if tm aconvc zero_tm then poly_0
209 else if RealArith.is_ratconst tm
210 then poly_const(RealArith.dest_ratconst tm)
212 (let val (lop,r) = Thm.dest_comb tm
213 in if lop aconvc neg_tm then poly_neg(poly_of_term r)
214 else if lop aconvc inv_tm then
215 let val p = poly_of_term r
217 then poly_const(Rat.inv (eval FuncUtil.Ctermfunc.empty p))
218 else error "poly_of_term: inverse of non-constant polyomial"
220 else (let val (opr,l) = Thm.dest_comb lop
222 if opr aconvc pow_tm andalso is_numeral r
223 then poly_pow (poly_of_term l) ((snd o HOLogic.dest_number o term_of) r)
224 else if opr aconvc add_tm
225 then poly_add (poly_of_term l) (poly_of_term r)
226 else if opr aconvc sub_tm
227 then poly_sub (poly_of_term l) (poly_of_term r)
228 else if opr aconvc mul_tm
229 then poly_mul (poly_of_term l) (poly_of_term r)
230 else if opr aconvc div_tm
232 val p = poly_of_term l
233 val q = poly_of_term r
234 in if poly_isconst q then poly_cmul (Rat.inv (eval FuncUtil.Ctermfunc.empty q)) p
235 else error "poly_of_term: division by non-constant polynomial"
240 handle CTERM ("dest_comb",_) => poly_var tm)
242 handle CTERM ("dest_comb",_) => poly_var tm)
244 val poly_of_term = fn tm =>
245 if type_of (term_of tm) = @{typ real} then poly_of_term tm
246 else error "poly_of_term: term does not have real type"
249 (* String of vector (just a list of space-separated numbers). *)
251 fun sdpa_of_vector (v:vector) =
254 val strs = map (decimalize 20 o (fn i => FuncUtil.Intfunc.tryapplyd (snd v) i rat_0)) (1 upto n)
255 in space_implode " " strs ^ "\n"
258 fun triple_int_ord ((a,b,c),(a',b',c')) =
259 prod_ord int_ord (prod_ord int_ord int_ord)
260 ((a,(b,c)),(a',(b',c')));
261 structure Inttriplefunc = FuncFun(type key = int*int*int val ord = triple_int_ord);
263 fun index_char str chr pos =
264 if pos >= String.size str then ~1
265 else if String.sub(str,pos) = chr then pos
266 else index_char str chr (pos + 1);
267 fun rat_of_quotient (a,b) = if b = 0 then rat_0 else Rat.rat_of_quotient (a,b);
268 fun rat_of_string s =
269 let val n = index_char s #"/" 0 in
270 if n = ~1 then s |> Int.fromString |> the |> Rat.rat_of_int
272 let val SOME numer = Int.fromString(String.substring(s,0,n))
273 val SOME den = Int.fromString (String.substring(s,n+1,String.size s - n - 1))
274 in rat_of_quotient(numer, den)
278 fun isnum x = member (op =) ["0","1","2","3","4","5","6","7","8","9"] x;
280 (* More parser basics. *)
282 val numeral = Scan.one isnum
283 val decimalint = Scan.repeat1 numeral >> (rat_of_string o implode)
284 val decimalfrac = Scan.repeat1 numeral
285 >> (fn s => rat_of_string(implode s) // pow10 (length s))
287 decimalint -- Scan.option (Scan.$$ "." |-- decimalfrac)
288 >> (fn (h,NONE) => h | (h,SOME x) => h +/ x)
290 $$ "-" |-- prs >> Rat.neg
294 fun emptyin def xs = if null xs then (def,xs) else Scan.fail xs
296 val exponent = ($$ "e" || $$ "E") |-- signed decimalint;
298 val decimal = signed decimalsig -- (emptyin rat_0|| exponent)
299 >> (fn (h, x) => h */ pow10 (int_of_rat x));
302 let val (x,rst) = p (raw_explode s)
303 in if null rst then x
304 else error "mkparser: unparsed input"
307 (* Parse back csdp output. *)
309 fun ignore _ = ((),[])
311 ((decimal -- Scan.repeat (Scan.$$ " " |-- Scan.option decimal) >>
312 (fn (h,to) => map_filter I ((SOME h)::to))) --| ignore >> vector_of_list) inp
313 val parse_csdpoutput = mkparser csdpoutput
315 (* Try some apparently sensible scaling first. Note that this is purely to *)
316 (* get a cleaner translation to floating-point, and doesn't affect any of *)
317 (* the results, in principle. In practice it seems a lot better when there *)
318 (* are extreme numbers in the original problem. *)
320 (* Version for (int*int*int) keys *)
322 fun max_rat x y = if x </ y then y else x
323 fun common_denominator fld amat acc =
324 fld (fn (_,c) => fn a => lcm_rat (denominator_rat c) a) amat acc
325 fun maximal_element fld amat acc =
326 fld (fn (_,c) => fn maxa => max_rat maxa (abs_rat c)) amat acc
327 fun float_of_rat x = let val (a,b) = Rat.quotient_of_rat x
328 in Real.fromInt a / Real.fromInt b end;
329 fun int_of_float x = (trunc x handle Overflow => 0 | Domain => 0)
332 fun tri_scale_then solver (obj:vector) mats =
334 val cd1 = fold_rev (common_denominator Inttriplefunc.fold) mats (rat_1)
335 val cd2 = common_denominator FuncUtil.Intfunc.fold (snd obj) (rat_1)
336 val mats' = map (Inttriplefunc.map (fn _ => fn x => cd1 */ x)) mats
337 val obj' = vector_cmul cd2 obj
338 val max1 = fold_rev (maximal_element Inttriplefunc.fold) mats' (rat_0)
339 val max2 = maximal_element FuncUtil.Intfunc.fold (snd obj') (rat_0)
340 val scal1 = pow2 (20 - int_of_float(Math.ln (float_of_rat max1) / Math.ln 2.0))
341 val scal2 = pow2 (20 - int_of_float(Math.ln (float_of_rat max2) / Math.ln 2.0))
342 val mats'' = map (Inttriplefunc.map (fn _ => fn x => x */ scal1)) mats'
343 val obj'' = vector_cmul scal2 obj'
344 in solver obj'' mats''
348 (* Round a vector to "nice" rationals. *)
350 fun nice_rational n x = round_rat (n */ x) // n;;
351 fun nice_vector n ((d,v) : vector) =
352 (d, FuncUtil.Intfunc.fold (fn (i,c) => fn a =>
353 let val y = nice_rational n c
354 in if c =/ rat_0 then a
355 else FuncUtil.Intfunc.update (i,y) a end) v FuncUtil.Intfunc.empty):vector
357 fun dest_ord f x = is_equal (f x);
359 (* Stuff for "equations" ((int*int*int)->num functions). *)
361 fun tri_equation_cmul c eq =
362 if c =/ rat_0 then Inttriplefunc.empty else Inttriplefunc.map (fn _ => fn d => c */ d) eq;
364 fun tri_equation_add eq1 eq2 = Inttriplefunc.combine (curry op +/) (fn x => x =/ rat_0) eq1 eq2;
366 fun tri_equation_eval assig eq =
367 let fun value v = Inttriplefunc.apply assig v
368 in Inttriplefunc.fold (fn (v, c) => fn a => a +/ value v */ c) eq rat_0
371 (* Eliminate all variables, in an essentially arbitrary order. *)
373 fun tri_eliminate_all_equations one =
375 fun choose_variable eq =
376 let val (v,_) = Inttriplefunc.choose eq
377 in if is_equal (triple_int_ord(v,one)) then
378 let val eq' = Inttriplefunc.delete_safe v eq
379 in if Inttriplefunc.is_empty eq' then error "choose_variable"
380 else fst (Inttriplefunc.choose eq')
384 fun eliminate dun eqs = case eqs of
387 if Inttriplefunc.is_empty eq then eliminate dun oeqs else
388 let val v = choose_variable eq
389 val a = Inttriplefunc.apply eq v
390 val eq' = tri_equation_cmul ((Rat.rat_of_int ~1) // a)
391 (Inttriplefunc.delete_safe v eq)
393 let val b = Inttriplefunc.tryapplyd e v rat_0
394 in if b =/ rat_0 then e
395 else tri_equation_add e (tri_equation_cmul (Rat.neg b // a) eq)
397 in eliminate (Inttriplefunc.update(v, eq') (Inttriplefunc.map (K elim) dun))
402 val assig = eliminate Inttriplefunc.empty eqs
403 val vs = Inttriplefunc.fold (fn (_, f) => fn a => remove (dest_ord triple_int_ord) one (Inttriplefunc.dom f) @ a) assig []
404 in (distinct (dest_ord triple_int_ord) vs,assig)
408 (* Multiply equation-parametrized poly by regular poly and add accumulator. *)
410 fun tri_epoly_pmul p q acc =
411 FuncUtil.Monomialfunc.fold (fn (m1, c) => fn a =>
412 FuncUtil.Monomialfunc.fold (fn (m2,e) => fn b =>
413 let val m = monomial_mul m1 m2
414 val es = FuncUtil.Monomialfunc.tryapplyd b m Inttriplefunc.empty
415 in FuncUtil.Monomialfunc.update (m,tri_equation_add (tri_equation_cmul c e) es) b
418 (* Hence produce the "relevant" monomials: those whose squares lie in the *)
419 (* Newton polytope of the monomials in the input. (This is enough according *)
420 (* to Reznik: "Extremal PSD forms with few terms", Duke Math. Journal, *)
421 (* vol 45, pp. 363--374, 1978. *)
423 (* These are ordered in sort of decreasing degree. In particular the *)
424 (* constant monomial is last; this gives an order in diagonalization of the *)
425 (* quadratic form that will tend to display constants. *)
427 (* Diagonalize (Cholesky/LDU) the matrix corresponding to a quadratic form. *)
430 fun diagonalize n i m =
431 if FuncUtil.Intpairfunc.is_empty (snd m) then []
433 let val a11 = FuncUtil.Intpairfunc.tryapplyd (snd m) (i,i) rat_0
434 in if a11 </ rat_0 then raise Failure "diagonalize: not PSD"
435 else if a11 =/ rat_0 then
436 if FuncUtil.Intfunc.is_empty (snd (row i m)) then diagonalize n (i + 1) m
437 else raise Failure "diagonalize: not PSD ___ "
441 val v' = (fst v, FuncUtil.Intfunc.fold (fn (i, c) => fn a =>
443 in if y = rat_0 then a else FuncUtil.Intfunc.update (i,y) a
444 end) (snd v) FuncUtil.Intfunc.empty)
445 fun upt0 x y a = if y = rat_0 then a else FuncUtil.Intpairfunc.update (x,y) a
448 iter (i+1,n) (fn j =>
449 iter (i+1,n) (fn k =>
450 (upt0 (j,k) (FuncUtil.Intpairfunc.tryapplyd (snd m) (j,k) rat_0 -/ FuncUtil.Intfunc.tryapplyd (snd v) j rat_0 */ FuncUtil.Intfunc.tryapplyd (snd v') k rat_0))))
451 FuncUtil.Intpairfunc.empty)
452 in (a11,v')::diagonalize n (i + 1) m'
458 val nn = dimensions m
460 in if snd nn <> n then error "diagonalize: non-square matrix"
461 else diagonalize n 1 m
465 (* Enumeration of monomials with given multidegree bound. *)
467 fun enumerate_monomials d vars =
469 else if d = 0 then [FuncUtil.Ctermfunc.empty]
470 else if null vars then [monomial_1] else
472 map_range (fn k => let val oths = enumerate_monomials (d - k) (tl vars)
473 in map (fn ks => if k = 0 then ks else FuncUtil.Ctermfunc.update (hd vars, k) ks) oths end) (d + 1)
477 (* Enumerate products of distinct input polys with degree <= d. *)
478 (* We ignore any constant input polynomials. *)
479 (* Give the output polynomial and a record of how it was derived. *)
481 fun enumerate_products d pols =
482 if d = 0 then [(poly_const rat_1,RealArith.Rational_lt rat_1)]
483 else if d < 0 then [] else
485 [] => [(poly_const rat_1,RealArith.Rational_lt rat_1)]
487 let val e = multidegree p
488 in if e = 0 then enumerate_products d ps else
489 enumerate_products d ps @
490 map (fn (q,c) => (poly_mul p q,RealArith.Product(b,c)))
491 (enumerate_products (d - e) ps)
494 (* Convert regular polynomial. Note that we treat (0,0,0) as -1. *)
496 fun epoly_of_poly p =
497 FuncUtil.Monomialfunc.fold (fn (m,c) => fn a => FuncUtil.Monomialfunc.update (m, Inttriplefunc.onefunc ((0,0,0), Rat.neg c)) a) p FuncUtil.Monomialfunc.empty;
499 (* String for block diagonal matrix numbered k. *)
501 fun sdpa_of_blockdiagonal k m =
503 val pfx = string_of_int k ^" "
506 (fn ((b,i,j),c) => fn a => if i > j then a else ((b,i,j),c)::a)
508 val entss = sort (triple_int_ord o pairself fst) ents
509 in fold_rev (fn ((b,i,j),c) => fn a =>
510 pfx ^ string_of_int b ^ " " ^ string_of_int i ^ " " ^ string_of_int j ^
511 " " ^ decimalize 20 c ^ "\n" ^ a) entss ""
514 (* SDPA for problem using block diagonal (i.e. multiple SDPs) *)
516 fun sdpa_of_blockproblem nblocks blocksizes obj mats =
517 let val m = length mats - 1
519 string_of_int m ^ "\n" ^
520 string_of_int nblocks ^ "\n" ^
521 (space_implode " " (map string_of_int blocksizes)) ^
524 fold_rev2 (fn k => fn m => fn a => sdpa_of_blockdiagonal (k - 1) m ^ a)
525 (1 upto length mats) mats ""
528 (* Run prover on a problem in block diagonal form. *)
530 fun run_blockproblem prover nblocks blocksizes obj mats=
531 parse_csdpoutput (prover (sdpa_of_blockproblem nblocks blocksizes obj mats))
533 (* 3D versions of matrix operations to consider blocks separately. *)
535 val bmatrix_add = Inttriplefunc.combine (curry op +/) (fn x => x =/ rat_0);
536 fun bmatrix_cmul c bm =
537 if c =/ rat_0 then Inttriplefunc.empty
538 else Inttriplefunc.map (fn _ => fn x => c */ x) bm;
540 val bmatrix_neg = bmatrix_cmul (Rat.rat_of_int ~1);
542 (* Smash a block matrix into components. *)
544 fun blocks blocksizes bm =
546 let val m = Inttriplefunc.fold
547 (fn ((b,i,j),c) => fn a => if b = b0 then FuncUtil.Intpairfunc.update ((i,j),c) a else a) bm FuncUtil.Intpairfunc.empty
548 val _ = FuncUtil.Intpairfunc.fold (fn ((i,j),_) => fn a => max a (max i j)) m 0
549 in (((bs,bs),m):matrix) end)
550 (blocksizes ~~ (1 upto length blocksizes));;
552 (* FIXME : Get rid of this !!!*)
554 fun tryfind_with msg _ [] = raise Failure msg
555 | tryfind_with _ f (x::xs) = (f x handle Failure s => tryfind_with s f xs);
557 fun tryfind f = tryfind_with "tryfind" f
560 (* Positiv- and Nullstellensatz. Flag "linf" forces a linear representation. *)
563 fun real_positivnullstellensatz_general ctxt prover linf d eqs leqs pol =
565 val vars = fold_rev (union (op aconvc) o poly_variables)
566 (pol :: eqs @ map fst leqs) []
567 val monoid = if linf then
568 (poly_const rat_1,RealArith.Rational_lt rat_1)::
569 (filter (fn (p,_) => multidegree p <= d) leqs)
570 else enumerate_products d leqs
571 val nblocks = length monoid
572 fun mk_idmultiplier k p =
574 val e = d - multidegree p
575 val mons = enumerate_monomials e vars
576 val nons = mons ~~ (1 upto length mons)
578 fold_rev (fn (m,n) => FuncUtil.Monomialfunc.update(m,Inttriplefunc.onefunc((~k,~n,n),rat_1))) nons FuncUtil.Monomialfunc.empty)
581 fun mk_sqmultiplier k (p,_) =
583 val e = (d - multidegree p) div 2
584 val mons = enumerate_monomials e vars
585 val nons = mons ~~ (1 upto length mons)
587 fold_rev (fn (m1,n1) =>
588 fold_rev (fn (m2,n2) => fn a =>
589 let val m = monomial_mul m1 m2
590 in if n1 > n2 then a else
591 let val c = if n1 = n2 then rat_1 else rat_2
592 val e = FuncUtil.Monomialfunc.tryapplyd a m Inttriplefunc.empty
593 in FuncUtil.Monomialfunc.update(m, tri_equation_add (Inttriplefunc.onefunc((k,n1,n2), c)) e) a
596 nons FuncUtil.Monomialfunc.empty)
599 val (sqmonlist,sqs) = split_list (map2 mk_sqmultiplier (1 upto length monoid) monoid)
600 val (_(*idmonlist*),ids) = split_list(map2 mk_idmultiplier (1 upto length eqs) eqs)
601 val blocksizes = map length sqmonlist
603 fold_rev2 (fn p => fn q => fn a => tri_epoly_pmul p q a) eqs ids
604 (fold_rev2 (fn (p,_) => fn s => fn a => tri_epoly_pmul p s a) monoid sqs
605 (epoly_of_poly(poly_neg pol)))
606 val eqns = FuncUtil.Monomialfunc.fold (fn (_,e) => fn a => e::a) bigsum []
607 val (pvs,assig) = tri_eliminate_all_equations (0,0,0) eqns
608 val qvars = (0,0,0)::pvs
609 val allassig = fold_rev (fn v => Inttriplefunc.update(v,(Inttriplefunc.onefunc(v,rat_1)))) pvs assig
611 Inttriplefunc.fold (fn ((b,i,j), ass) => fn m =>
613 let val c = Inttriplefunc.tryapplyd ass v rat_0
614 in if c = rat_0 then m else
615 Inttriplefunc.update ((b,j,i), c) (Inttriplefunc.update ((b,i,j), c) m)
617 allassig Inttriplefunc.empty
618 val diagents = Inttriplefunc.fold
619 (fn ((b,i,j), e) => fn a => if b > 0 andalso i = j then tri_equation_add e a else a)
620 allassig Inttriplefunc.empty
622 val mats = map mk_matrix qvars
623 val obj = (length pvs,
624 itern 1 pvs (fn v => fn i => FuncUtil.Intfunc.updatep iszero (i,Inttriplefunc.tryapplyd diagents v rat_0))
625 FuncUtil.Intfunc.empty)
626 val raw_vec = if null pvs then vector_0 0
627 else tri_scale_then (run_blockproblem prover nblocks blocksizes) obj mats
628 fun int_element (_,v) i = FuncUtil.Intfunc.tryapplyd v i rat_0
630 fun find_rounding d =
633 if Config.get ctxt trace
634 then writeln ("Trying rounding with limit "^Rat.string_of_rat d ^ "\n")
636 val vec = nice_vector d raw_vec
637 val blockmat = iter (1,dim vec)
638 (fn i => fn a => bmatrix_add (bmatrix_cmul (int_element vec i) (nth mats i)) a)
639 (bmatrix_neg (nth mats 0))
640 val allmats = blocks blocksizes blockmat
641 in (vec,map diag allmats)
644 if null pvs then find_rounding rat_1
645 else tryfind find_rounding (map Rat.rat_of_int (1 upto 31) @
646 map pow2 (5 upto 66))
648 fold_rev (fn k => Inttriplefunc.update (nth pvs (k - 1), int_element vec k))
649 (1 upto dim vec) (Inttriplefunc.onefunc ((0,0,0), Rat.rat_of_int ~1))
651 Inttriplefunc.fold (fn (v,e) => fn a => Inttriplefunc.update(v, tri_equation_eval newassigs e) a) allassig newassigs
652 fun poly_of_epoly p =
653 FuncUtil.Monomialfunc.fold (fn (v,e) => fn a => FuncUtil.Monomialfunc.updatep iszero (v,tri_equation_eval finalassigs e) a)
654 p FuncUtil.Monomialfunc.empty
656 let fun mk_sq (c,m) =
657 (c,fold_rev (fn k=> fn a => FuncUtil.Monomialfunc.updatep iszero (nth mons (k - 1), int_element m k) a)
658 (1 upto length mons) FuncUtil.Monomialfunc.empty)
661 val sqs = map2 mk_sos sqmonlist ratdias
662 val cfs = map poly_of_epoly ids
663 val msq = filter (fn (_,b) => not (null b)) (map2 pair monoid sqs)
664 fun eval_sq sqs = fold_rev (fn (c,q) => poly_add (poly_cmul c (poly_mul q q))) sqs poly_0
666 fold_rev (fn ((p,_),s) => poly_add (poly_mul p (eval_sq s))) msq
667 (fold_rev2 (fn p => fn q => poly_add (poly_mul p q)) cfs eqs
670 in if not(FuncUtil.Monomialfunc.is_empty sanity) then raise Sanity else
671 (cfs,map (fn (a,b) => (snd a,b)) msq)
675 (* Iterative deepening. *)
678 (writeln ("Searching with depth limit " ^ string_of_int n);
679 (f n handle Failure s => (writeln ("failed with message: " ^ s); deepen f (n + 1))));
682 (* Map back polynomials and their composites to a positivstellensatz. *)
684 fun cterm_of_sqterm (c,p) = RealArith.Product(RealArith.Rational_lt c,RealArith.Square p);
686 fun cterm_of_sos (pr,sqs) = if null sqs then pr
687 else RealArith.Product(pr,foldr1 RealArith.Sum (map cterm_of_sqterm sqs));
689 (* Interface to HOL. *)
692 val concl = Thm.dest_arg o cprop_of
693 fun simple_cterm_ord t u = Term_Ord.fast_term_ord (term_of t, term_of u) = LESS
695 (* FIXME: Replace tryfind by get_first !! *)
696 fun real_nonlinear_prover proof_method ctxt =
698 val {add = _, mul = _, neg = _, pow = _,
699 sub = _, main = real_poly_conv} =
700 Semiring_Normalizer.semiring_normalizers_ord_wrapper ctxt
701 (the (Semiring_Normalizer.match ctxt @{cterm "(0::real) + 1"}))
703 fun mainf cert_choice translator (eqs,les,lts) =
705 val eq0 = map (poly_of_term o Thm.dest_arg1 o concl) eqs
706 val le0 = map (poly_of_term o Thm.dest_arg o concl) les
707 val lt0 = map (poly_of_term o Thm.dest_arg o concl) lts
708 val eqp0 = map_index (fn (i, t) => (t,RealArith.Axiom_eq i)) eq0
709 val lep0 = map_index (fn (i, t) => (t,RealArith.Axiom_le i)) le0
710 val ltp0 = map_index (fn (i, t) => (t,RealArith.Axiom_lt i)) lt0
711 val (keq,eq) = List.partition (fn (p,_) => multidegree p = 0) eqp0
712 val (klep,lep) = List.partition (fn (p,_) => multidegree p = 0) lep0
713 val (kltp,ltp) = List.partition (fn (p,_) => multidegree p = 0) ltp0
714 fun trivial_axiom (p,ax) =
716 RealArith.Axiom_eq n => if eval FuncUtil.Ctermfunc.empty p <>/ Rat.zero then nth eqs n
717 else raise Failure "trivial_axiom: Not a trivial axiom"
718 | RealArith.Axiom_le n => if eval FuncUtil.Ctermfunc.empty p </ Rat.zero then nth les n
719 else raise Failure "trivial_axiom: Not a trivial axiom"
720 | RealArith.Axiom_lt n => if eval FuncUtil.Ctermfunc.empty p <=/ Rat.zero then nth lts n
721 else raise Failure "trivial_axiom: Not a trivial axiom"
722 | _ => error "trivial_axiom: Not a trivial axiom"
724 (let val th = tryfind trivial_axiom (keq @ klep @ kltp)
726 (fconv_rule (arg_conv (arg1_conv real_poly_conv) then_conv Numeral_Simprocs.field_comp_conv) th, RealArith.Trivial)
730 (case proof_method of Certificate certs =>
731 (* choose certificate *)
733 fun chose_cert [] (RealArith.Cert c) = c
734 | chose_cert (RealArith.Left::s) (RealArith.Branch (l, _)) = chose_cert s l
735 | chose_cert (RealArith.Right::s) (RealArith.Branch (_, r)) = chose_cert s r
736 | chose_cert _ _ = error "certificate tree in invalid form"
738 chose_cert cert_choice certs
743 val pol = fold_rev poly_mul (map fst ltp) (poly_const Rat.one)
746 let val e = multidegree pol
747 val k = if e = 0 then 0 else d div e
749 in tryfind (fn i => (d,i,real_positivnullstellensatz_general ctxt prover false d eq' leq
750 (poly_neg(poly_pow pol i))))
753 val (_,i,(cert_ideal,cert_cone)) = deepen tryall 0
755 map2 (fn q => fn (_,ax) => RealArith.Eqmul(q,ax)) cert_ideal eq
756 val proofs_cone = map cterm_of_sos cert_cone
757 val proof_ne = if null ltp then RealArith.Rational_lt Rat.one else
758 let val p = foldr1 RealArith.Product (map snd ltp)
759 in funpow i (fn q => RealArith.Product(p,q)) (RealArith.Rational_lt Rat.one)
762 foldr1 RealArith.Sum (proof_ne :: proofs_ideal @ proofs_cone)
765 (translator (eqs,les,lts) proof, RealArith.Cert proof)
772 (* FIXME : This is very bad!!!*)
773 fun subst_conv eqs t =
775 val t' = fold (Thm.lambda o Thm.lhs_of) eqs t
776 in Conv.fconv_rule (Thm.beta_conversion true) (fold (C Thm.combination) eqs (Thm.reflexive t'))
779 (* A wrapper that tries to substitute away variables first. *)
783 fun simple_cterm_ord t u = Term_Ord.fast_term_ord (term_of t, term_of u) = LESS
784 val concl = Thm.dest_arg o cprop_of
786 fconv_rule (rewr_conv @{lemma "(a + x == y) == (x == y - (a::real))" by (atomize (full)) (simp add: field_simps) })
788 fconv_rule (rewr_conv @{lemma "(x + a == y) == (x == y - (a::real))" by (atomize (full)) (simp add: field_simps)})
789 fun substitutable_monomial fvs tm = case term_of tm of
790 Free(_,@{typ real}) => if not (member (op aconvc) fvs tm) then (Rat.one,tm)
791 else raise Failure "substitutable_monomial"
792 | @{term "op * :: real => _"}$_$(Free _) =>
793 if RealArith.is_ratconst (Thm.dest_arg1 tm) andalso not (member (op aconvc) fvs (Thm.dest_arg tm))
794 then (RealArith.dest_ratconst (Thm.dest_arg1 tm),Thm.dest_arg tm) else raise Failure "substitutable_monomial"
795 | @{term "op + :: real => _"}$_$_ =>
796 (substitutable_monomial (Thm.add_cterm_frees (Thm.dest_arg tm) fvs) (Thm.dest_arg1 tm)
797 handle Failure _ => substitutable_monomial (Thm.add_cterm_frees (Thm.dest_arg1 tm) fvs) (Thm.dest_arg tm))
798 | _ => raise Failure "substitutable_monomial"
800 fun isolate_variable v th =
801 let val w = Thm.dest_arg1 (cprop_of th)
802 in if v aconvc w then th
803 else case term_of w of
804 @{term "op + :: real => _"}$_$_ =>
805 if Thm.dest_arg1 w aconvc v then shuffle2 th
806 else isolate_variable v (shuffle1 th)
807 | _ => error "isolate variable : This should not happen?"
811 fun real_nonlinear_subst_prover prover ctxt =
813 val {add = _, mul = real_poly_mul_conv, neg = _,
814 pow = _, sub = _, main = real_poly_conv} =
815 Semiring_Normalizer.semiring_normalizers_ord_wrapper ctxt
816 (the (Semiring_Normalizer.match ctxt @{cterm "(0::real) + 1"}))
819 fun make_substitution th =
821 val (c,v) = substitutable_monomial [] (Thm.dest_arg1(concl th))
822 val th1 = Drule.arg_cong_rule (Thm.apply @{cterm "op * :: real => _"} (RealArith.cterm_of_rat (Rat.inv c))) (mk_meta_eq th)
823 val th2 = fconv_rule (binop_conv real_poly_mul_conv) th1
824 in fconv_rule (arg_conv real_poly_conv) (isolate_variable v th2)
827 let val g = Thm.dest_fun2 ct
828 in if g aconvc @{cterm "op <= :: real => _"}
829 orelse g aconvc @{cterm "op < :: real => _"}
830 then arg_conv cv ct else arg1_conv cv ct
832 fun mainf cert_choice translator =
834 fun substfirst(eqs,les,lts) =
836 val eth = tryfind make_substitution eqs
837 val modify = fconv_rule (arg_conv (oprconv(subst_conv [eth] then_conv real_poly_conv)))
839 (filter_out (fn t => (Thm.dest_arg1 o Thm.dest_arg o cprop_of) t
840 aconvc @{cterm "0::real"}) (map modify eqs),
841 map modify les,map modify lts)
843 handle Failure _ => real_nonlinear_prover prover ctxt cert_choice translator (rev eqs, rev les, rev lts))
851 (* Overall function. *)
853 fun real_sos prover ctxt =
854 RealArith.gen_prover_real_arith ctxt (real_nonlinear_subst_prover prover ctxt)
857 val known_sos_constants =
858 [@{term "op ==>"}, @{term "Trueprop"},
859 @{term HOL.implies}, @{term HOL.conj}, @{term HOL.disj},
860 @{term "Not"}, @{term "op = :: bool => _"},
861 @{term "All :: (real => _) => _"}, @{term "Ex :: (real => _) => _"},
862 @{term "op = :: real => _"}, @{term "op < :: real => _"},
863 @{term "op <= :: real => _"},
864 @{term "op + :: real => _"}, @{term "op - :: real => _"},
865 @{term "op * :: real => _"}, @{term "uminus :: real => _"},
866 @{term "op / :: real => _"}, @{term "inverse :: real => _"},
867 @{term "op ^ :: real => _"}, @{term "abs :: real => _"},
868 @{term "min :: real => _"}, @{term "max :: real => _"},
869 @{term "0::real"}, @{term "1::real"}, @{term "number_of :: int => real"},
870 @{term "number_of :: int => nat"},
871 @{term "Int.Bit0"}, @{term "Int.Bit1"},
872 @{term "Int.Pls"}, @{term "Int.Min"}];
874 fun check_sos kcts ct =
877 val _ = if not (null (Term.add_tfrees t [])
878 andalso null (Term.add_tvars t []))
879 then error "SOS: not sos. Additional type varables" else ()
880 val fs = Term.add_frees t []
881 val _ = if exists (fn ((_,T)) => not (T = @{typ "real"})) fs
882 then error "SOS: not sos. Variables with type not real" else ()
883 val vs = Term.add_vars t []
884 val _ = if exists (fn ((_,T)) => not (T = @{typ "real"})) vs
885 then error "SOS: not sos. Variables with type not real" else ()
886 val ukcs = subtract (fn (t,p) => Const p aconv t) kcts (Term.add_consts t [])
887 val _ = if null ukcs then ()
888 else error ("SOSO: Unknown constants in Subgoal:" ^ commas (map fst ukcs))
891 fun core_sos_tac print_cert prover = SUBPROOF (fn {concl, context, ...} =>
893 val _ = check_sos known_sos_constants concl
894 val (ths, certificates) = real_sos prover context (Thm.dest_arg concl)
895 val _ = print_cert certificates
898 fun default_SOME _ NONE v = SOME v
899 | default_SOME _ (SOME v) _ = SOME v;
901 fun lift_SOME f NONE a = f a
902 | lift_SOME _ (SOME a) _ = SOME a;
906 val is_numeral = can (HOLogic.dest_number o term_of)
908 fun get_denom b ct = case term_of ct of
909 @{term "op / :: real => _"} $ _ $ _ =>
910 if is_numeral (Thm.dest_arg ct) then get_denom b (Thm.dest_arg1 ct)
911 else default_SOME (get_denom b) (get_denom b (Thm.dest_arg ct)) (Thm.dest_arg ct, b)
912 | @{term "op < :: real => _"} $ _ $ _ => lift_SOME (get_denom true) (get_denom true (Thm.dest_arg ct)) (Thm.dest_arg1 ct)
913 | @{term "op <= :: real => _"} $ _ $ _ => lift_SOME (get_denom true) (get_denom true (Thm.dest_arg ct)) (Thm.dest_arg1 ct)
914 | _ $ _ => lift_SOME (get_denom b) (get_denom b (Thm.dest_fun ct)) (Thm.dest_arg ct)
918 fun elim_one_denom_tac ctxt =
919 CSUBGOAL (fn (P,i) =>
920 case get_denom false P of
924 val ss = simpset_of ctxt addsimps @{thms field_simps}
925 addsimps [@{thm nonzero_power_divide}, @{thm power_divide}]
926 val th = instantiate' [] [SOME d, SOME (Thm.dest_arg P)]
927 (if ord then @{lemma "(d=0 --> P) & (d>0 --> P) & (d<(0::real) --> P) ==> P" by auto}
928 else @{lemma "(d=0 --> P) & (d ~= (0::real) --> P) ==> P" by blast})
929 in rtac th i THEN Simplifier.asm_full_simp_tac ss i end);
931 fun elim_denom_tac ctxt i = REPEAT (elim_one_denom_tac ctxt i);
933 fun sos_tac print_cert prover ctxt =
934 Object_Logic.full_atomize_tac THEN'
935 elim_denom_tac ctxt THEN'
936 core_sos_tac print_cert prover ctxt;