1.1 --- a/src/HOL/Library/sum_of_squares.ML Wed Aug 05 17:10:10 2009 +0200
1.2 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000
1.3 @@ -1,1754 +0,0 @@
1.4 -(* Title: sum_of_squares.ML
1.5 - Authors: Amine Chaieb, University of Cambridge
1.6 - Philipp Meyer, TU Muenchen
1.7 -
1.8 -A tactic for proving nonlinear inequalities
1.9 -*)
1.10 -
1.11 -signature SOS =
1.12 -sig
1.13 -
1.14 - val sos_tac : (string -> string) -> Proof.context -> int -> Tactical.tactic
1.15 -
1.16 -end
1.17 -
1.18 -structure Sos : SOS =
1.19 -struct
1.20 -
1.21 -
1.22 -val rat_0 = Rat.zero;
1.23 -val rat_1 = Rat.one;
1.24 -val rat_2 = Rat.two;
1.25 -val rat_10 = Rat.rat_of_int 10;
1.26 -val rat_1_2 = rat_1 // rat_2;
1.27 -val max = curry IntInf.max;
1.28 -val min = curry IntInf.min;
1.29 -
1.30 -val denominator_rat = Rat.quotient_of_rat #> snd #> Rat.rat_of_int;
1.31 -val numerator_rat = Rat.quotient_of_rat #> fst #> Rat.rat_of_int;
1.32 -fun int_of_rat a =
1.33 - case Rat.quotient_of_rat a of (i,1) => i | _ => error "int_of_rat: not an int";
1.34 -fun lcm_rat x y = Rat.rat_of_int (Integer.lcm (int_of_rat x) (int_of_rat y));
1.35 -
1.36 -fun rat_pow r i =
1.37 - let fun pow r i =
1.38 - if i = 0 then rat_1 else
1.39 - let val d = pow r (i div 2)
1.40 - in d */ d */ (if i mod 2 = 0 then rat_1 else r)
1.41 - end
1.42 - in if i < 0 then pow (Rat.inv r) (~ i) else pow r i end;
1.43 -
1.44 -fun round_rat r =
1.45 - let val (a,b) = Rat.quotient_of_rat (Rat.abs r)
1.46 - val d = a div b
1.47 - val s = if r </ rat_0 then (Rat.neg o Rat.rat_of_int) else Rat.rat_of_int
1.48 - val x2 = 2 * (a - (b * d))
1.49 - in s (if x2 >= b then d + 1 else d) end
1.50 -
1.51 -val abs_rat = Rat.abs;
1.52 -val pow2 = rat_pow rat_2;
1.53 -val pow10 = rat_pow rat_10;
1.54 -
1.55 -val debugging = ref false;
1.56 -
1.57 -exception Sanity;
1.58 -
1.59 -exception Unsolvable;
1.60 -
1.61 -(* Turn a rational into a decimal string with d sig digits. *)
1.62 -
1.63 -local
1.64 -fun normalize y =
1.65 - if abs_rat y </ (rat_1 // rat_10) then normalize (rat_10 */ y) - 1
1.66 - else if abs_rat y >=/ rat_1 then normalize (y // rat_10) + 1
1.67 - else 0
1.68 - in
1.69 -fun decimalize d x =
1.70 - if x =/ rat_0 then "0.0" else
1.71 - let
1.72 - val y = Rat.abs x
1.73 - val e = normalize y
1.74 - val z = pow10(~ e) */ y +/ rat_1
1.75 - val k = int_of_rat (round_rat(pow10 d */ z))
1.76 - in (if x </ rat_0 then "-0." else "0.") ^
1.77 - implode(tl(explode(string_of_int k))) ^
1.78 - (if e = 0 then "" else "e"^string_of_int e)
1.79 - end
1.80 -end;
1.81 -
1.82 -(* Iterations over numbers, and lists indexed by numbers. *)
1.83 -
1.84 -fun itern k l f a =
1.85 - case l of
1.86 - [] => a
1.87 - | h::t => itern (k + 1) t f (f h k a);
1.88 -
1.89 -fun iter (m,n) f a =
1.90 - if n < m then a
1.91 - else iter (m+1,n) f (f m a);
1.92 -
1.93 -(* The main types. *)
1.94 -
1.95 -fun strict_ord ord (x,y) = case ord (x,y) of LESS => LESS | _ => GREATER
1.96 -
1.97 -structure Intpairfunc = FuncFun(type key = int*int val ord = prod_ord int_ord int_ord);
1.98 -
1.99 -type vector = int* Rat.rat Intfunc.T;
1.100 -
1.101 -type matrix = (int*int)*(Rat.rat Intpairfunc.T);
1.102 -
1.103 -type monomial = int Ctermfunc.T;
1.104 -
1.105 -val cterm_ord = (fn (s,t) => TermOrd.fast_term_ord(term_of s, term_of t))
1.106 - fun monomial_ord (m1,m2) = list_ord (prod_ord cterm_ord int_ord) (Ctermfunc.graph m1, Ctermfunc.graph m2)
1.107 -structure Monomialfunc = FuncFun(type key = monomial val ord = monomial_ord)
1.108 -
1.109 -type poly = Rat.rat Monomialfunc.T;
1.110 -
1.111 - fun iszero (k,r) = r =/ rat_0;
1.112 -
1.113 -fun fold_rev2 f l1 l2 b =
1.114 - case (l1,l2) of
1.115 - ([],[]) => b
1.116 - | (h1::t1,h2::t2) => f h1 h2 (fold_rev2 f t1 t2 b)
1.117 - | _ => error "fold_rev2";
1.118 -
1.119 -(* Vectors. Conventionally indexed 1..n. *)
1.120 -
1.121 -fun vector_0 n = (n,Intfunc.undefined):vector;
1.122 -
1.123 -fun dim (v:vector) = fst v;
1.124 -
1.125 -fun vector_const c n =
1.126 - if c =/ rat_0 then vector_0 n
1.127 - else (n,fold_rev (fn k => Intfunc.update (k,c)) (1 upto n) Intfunc.undefined) :vector;
1.128 -
1.129 -val vector_1 = vector_const rat_1;
1.130 -
1.131 -fun vector_cmul c (v:vector) =
1.132 - let val n = dim v
1.133 - in if c =/ rat_0 then vector_0 n
1.134 - else (n,Intfunc.mapf (fn x => c */ x) (snd v))
1.135 - end;
1.136 -
1.137 -fun vector_neg (v:vector) = (fst v,Intfunc.mapf Rat.neg (snd v)) :vector;
1.138 -
1.139 -fun vector_add (v1:vector) (v2:vector) =
1.140 - let val m = dim v1
1.141 - val n = dim v2
1.142 - in if m <> n then error "vector_add: incompatible dimensions"
1.143 - else (n,Intfunc.combine (curry op +/) (fn x => x =/ rat_0) (snd v1) (snd v2)) :vector
1.144 - end;
1.145 -
1.146 -fun vector_sub v1 v2 = vector_add v1 (vector_neg v2);
1.147 -
1.148 -fun vector_dot (v1:vector) (v2:vector) =
1.149 - let val m = dim v1
1.150 - val n = dim v2
1.151 - in if m <> n then error "vector_dot: incompatible dimensions"
1.152 - else Intfunc.fold (fn (i,x) => fn a => x +/ a)
1.153 - (Intfunc.combine (curry op */) (fn x => x =/ rat_0) (snd v1) (snd v2)) rat_0
1.154 - end;
1.155 -
1.156 -fun vector_of_list l =
1.157 - let val n = length l
1.158 - in (n,fold_rev2 (curry Intfunc.update) (1 upto n) l Intfunc.undefined) :vector
1.159 - end;
1.160 -
1.161 -(* Matrices; again rows and columns indexed from 1. *)
1.162 -
1.163 -fun matrix_0 (m,n) = ((m,n),Intpairfunc.undefined):matrix;
1.164 -
1.165 -fun dimensions (m:matrix) = fst m;
1.166 -
1.167 -fun matrix_const c (mn as (m,n)) =
1.168 - if m <> n then error "matrix_const: needs to be square"
1.169 - else if c =/ rat_0 then matrix_0 mn
1.170 - else (mn,fold_rev (fn k => Intpairfunc.update ((k,k), c)) (1 upto n) Intpairfunc.undefined) :matrix;;
1.171 -
1.172 -val matrix_1 = matrix_const rat_1;
1.173 -
1.174 -fun matrix_cmul c (m:matrix) =
1.175 - let val (i,j) = dimensions m
1.176 - in if c =/ rat_0 then matrix_0 (i,j)
1.177 - else ((i,j),Intpairfunc.mapf (fn x => c */ x) (snd m))
1.178 - end;
1.179 -
1.180 -fun matrix_neg (m:matrix) =
1.181 - (dimensions m, Intpairfunc.mapf Rat.neg (snd m)) :matrix;
1.182 -
1.183 -fun matrix_add (m1:matrix) (m2:matrix) =
1.184 - let val d1 = dimensions m1
1.185 - val d2 = dimensions m2
1.186 - in if d1 <> d2
1.187 - then error "matrix_add: incompatible dimensions"
1.188 - else (d1,Intpairfunc.combine (curry op +/) (fn x => x =/ rat_0) (snd m1) (snd m2)) :matrix
1.189 - end;;
1.190 -
1.191 -fun matrix_sub m1 m2 = matrix_add m1 (matrix_neg m2);
1.192 -
1.193 -fun row k (m:matrix) =
1.194 - let val (i,j) = dimensions m
1.195 - in (j,
1.196 - Intpairfunc.fold (fn ((i,j), c) => fn a => if i = k then Intfunc.update (j,c) a else a) (snd m) Intfunc.undefined ) : vector
1.197 - end;
1.198 -
1.199 -fun column k (m:matrix) =
1.200 - let val (i,j) = dimensions m
1.201 - in (i,
1.202 - Intpairfunc.fold (fn ((i,j), c) => fn a => if j = k then Intfunc.update (i,c) a else a) (snd m) Intfunc.undefined)
1.203 - : vector
1.204 - end;
1.205 -
1.206 -fun transp (m:matrix) =
1.207 - let val (i,j) = dimensions m
1.208 - in
1.209 - ((j,i),Intpairfunc.fold (fn ((i,j), c) => fn a => Intpairfunc.update ((j,i), c) a) (snd m) Intpairfunc.undefined) :matrix
1.210 - end;
1.211 -
1.212 -fun diagonal (v:vector) =
1.213 - let val n = dim v
1.214 - in ((n,n),Intfunc.fold (fn (i, c) => fn a => Intpairfunc.update ((i,i), c) a) (snd v) Intpairfunc.undefined) : matrix
1.215 - end;
1.216 -
1.217 -fun matrix_of_list l =
1.218 - let val m = length l
1.219 - in if m = 0 then matrix_0 (0,0) else
1.220 - let val n = length (hd l)
1.221 - in ((m,n),itern 1 l (fn v => fn i => itern 1 v (fn c => fn j => Intpairfunc.update ((i,j), c))) Intpairfunc.undefined)
1.222 - end
1.223 - end;
1.224 -
1.225 -(* Monomials. *)
1.226 -
1.227 -fun monomial_eval assig (m:monomial) =
1.228 - Ctermfunc.fold (fn (x, k) => fn a => a */ rat_pow (Ctermfunc.apply assig x) k)
1.229 - m rat_1;
1.230 -val monomial_1 = (Ctermfunc.undefined:monomial);
1.231 -
1.232 -fun monomial_var x = Ctermfunc.onefunc (x, 1) :monomial;
1.233 -
1.234 -val (monomial_mul:monomial->monomial->monomial) =
1.235 - Ctermfunc.combine (curry op +) (K false);
1.236 -
1.237 -fun monomial_pow (m:monomial) k =
1.238 - if k = 0 then monomial_1
1.239 - else Ctermfunc.mapf (fn x => k * x) m;
1.240 -
1.241 -fun monomial_divides (m1:monomial) (m2:monomial) =
1.242 - Ctermfunc.fold (fn (x, k) => fn a => Ctermfunc.tryapplyd m2 x 0 >= k andalso a) m1 true;;
1.243 -
1.244 -fun monomial_div (m1:monomial) (m2:monomial) =
1.245 - let val m = Ctermfunc.combine (curry op +)
1.246 - (fn x => x = 0) m1 (Ctermfunc.mapf (fn x => ~ x) m2)
1.247 - in if Ctermfunc.fold (fn (x, k) => fn a => k >= 0 andalso a) m true then m
1.248 - else error "monomial_div: non-divisible"
1.249 - end;
1.250 -
1.251 -fun monomial_degree x (m:monomial) =
1.252 - Ctermfunc.tryapplyd m x 0;;
1.253 -
1.254 -fun monomial_lcm (m1:monomial) (m2:monomial) =
1.255 - fold_rev (fn x => Ctermfunc.update (x, max (monomial_degree x m1) (monomial_degree x m2)))
1.256 - (gen_union (is_equal o cterm_ord) (Ctermfunc.dom m1, Ctermfunc.dom m2)) (Ctermfunc.undefined :monomial);
1.257 -
1.258 -fun monomial_multidegree (m:monomial) =
1.259 - Ctermfunc.fold (fn (x, k) => fn a => k + a) m 0;;
1.260 -
1.261 -fun monomial_variables m = Ctermfunc.dom m;;
1.262 -
1.263 -(* Polynomials. *)
1.264 -
1.265 -fun eval assig (p:poly) =
1.266 - Monomialfunc.fold (fn (m, c) => fn a => a +/ c */ monomial_eval assig m) p rat_0;
1.267 -
1.268 -val poly_0 = (Monomialfunc.undefined:poly);
1.269 -
1.270 -fun poly_isconst (p:poly) =
1.271 - Monomialfunc.fold (fn (m, c) => fn a => Ctermfunc.is_undefined m andalso a) p true;
1.272 -
1.273 -fun poly_var x = Monomialfunc.onefunc (monomial_var x,rat_1) :poly;
1.274 -
1.275 -fun poly_const c =
1.276 - if c =/ rat_0 then poly_0 else Monomialfunc.onefunc(monomial_1, c);
1.277 -
1.278 -fun poly_cmul c (p:poly) =
1.279 - if c =/ rat_0 then poly_0
1.280 - else Monomialfunc.mapf (fn x => c */ x) p;
1.281 -
1.282 -fun poly_neg (p:poly) = (Monomialfunc.mapf Rat.neg p :poly);;
1.283 -
1.284 -fun poly_add (p1:poly) (p2:poly) =
1.285 - (Monomialfunc.combine (curry op +/) (fn x => x =/ rat_0) p1 p2 :poly);
1.286 -
1.287 -fun poly_sub p1 p2 = poly_add p1 (poly_neg p2);
1.288 -
1.289 -fun poly_cmmul (c,m) (p:poly) =
1.290 - if c =/ rat_0 then poly_0
1.291 - else if Ctermfunc.is_undefined m
1.292 - then Monomialfunc.mapf (fn d => c */ d) p
1.293 - else Monomialfunc.fold (fn (m', d) => fn a => (Monomialfunc.update (monomial_mul m m', c */ d) a)) p poly_0;
1.294 -
1.295 -fun poly_mul (p1:poly) (p2:poly) =
1.296 - Monomialfunc.fold (fn (m, c) => fn a => poly_add (poly_cmmul (c,m) p2) a) p1 poly_0;
1.297 -
1.298 -fun poly_div (p1:poly) (p2:poly) =
1.299 - if not(poly_isconst p2)
1.300 - then error "poly_div: non-constant" else
1.301 - let val c = eval Ctermfunc.undefined p2
1.302 - in if c =/ rat_0 then error "poly_div: division by zero"
1.303 - else poly_cmul (Rat.inv c) p1
1.304 - end;
1.305 -
1.306 -fun poly_square p = poly_mul p p;
1.307 -
1.308 -fun poly_pow p k =
1.309 - if k = 0 then poly_const rat_1
1.310 - else if k = 1 then p
1.311 - else let val q = poly_square(poly_pow p (k div 2)) in
1.312 - if k mod 2 = 1 then poly_mul p q else q end;
1.313 -
1.314 -fun poly_exp p1 p2 =
1.315 - if not(poly_isconst p2)
1.316 - then error "poly_exp: not a constant"
1.317 - else poly_pow p1 (int_of_rat (eval Ctermfunc.undefined p2));
1.318 -
1.319 -fun degree x (p:poly) =
1.320 - Monomialfunc.fold (fn (m,c) => fn a => max (monomial_degree x m) a) p 0;
1.321 -
1.322 -fun multidegree (p:poly) =
1.323 - Monomialfunc.fold (fn (m, c) => fn a => max (monomial_multidegree m) a) p 0;
1.324 -
1.325 -fun poly_variables (p:poly) =
1.326 - sort cterm_ord (Monomialfunc.fold_rev (fn (m, c) => curry (gen_union (is_equal o cterm_ord)) (monomial_variables m)) p []);;
1.327 -
1.328 -(* Order monomials for human presentation. *)
1.329 -
1.330 -fun cterm_ord (t,t') = TermOrd.fast_term_ord (term_of t, term_of t');
1.331 -
1.332 -val humanorder_varpow = prod_ord cterm_ord (rev_order o int_ord);
1.333 -
1.334 -local
1.335 - fun ord (l1,l2) = case (l1,l2) of
1.336 - (_,[]) => LESS
1.337 - | ([],_) => GREATER
1.338 - | (h1::t1,h2::t2) =>
1.339 - (case humanorder_varpow (h1, h2) of
1.340 - LESS => LESS
1.341 - | EQUAL => ord (t1,t2)
1.342 - | GREATER => GREATER)
1.343 -in fun humanorder_monomial m1 m2 =
1.344 - ord (sort humanorder_varpow (Ctermfunc.graph m1),
1.345 - sort humanorder_varpow (Ctermfunc.graph m2))
1.346 -end;
1.347 -
1.348 -fun fold1 f l = case l of
1.349 - [] => error "fold1"
1.350 - | [x] => x
1.351 - | (h::t) => f h (fold1 f t);
1.352 -
1.353 -(* Conversions to strings. *)
1.354 -
1.355 -fun string_of_vector min_size max_size (v:vector) =
1.356 - let val n_raw = dim v
1.357 - in if n_raw = 0 then "[]" else
1.358 - let
1.359 - val n = max min_size (min n_raw max_size)
1.360 - val xs = map (Rat.string_of_rat o (fn i => Intfunc.tryapplyd (snd v) i rat_0)) (1 upto n)
1.361 - in "[" ^ fold1 (fn s => fn t => s ^ ", " ^ t) xs ^
1.362 - (if n_raw > max_size then ", ...]" else "]")
1.363 - end
1.364 - end;
1.365 -
1.366 -fun string_of_matrix max_size (m:matrix) =
1.367 - let
1.368 - val (i_raw,j_raw) = dimensions m
1.369 - val i = min max_size i_raw
1.370 - val j = min max_size j_raw
1.371 - val rstr = map (fn k => string_of_vector j j (row k m)) (1 upto i)
1.372 - in "["^ fold1 (fn s => fn t => s^";\n "^t) rstr ^
1.373 - (if j > max_size then "\n ...]" else "]")
1.374 - end;
1.375 -
1.376 -fun string_of_term t =
1.377 - case t of
1.378 - a$b => "("^(string_of_term a)^" "^(string_of_term b)^")"
1.379 - | Abs x =>
1.380 - let val (xn, b) = Term.dest_abs x
1.381 - in "(\\"^xn^"."^(string_of_term b)^")"
1.382 - end
1.383 - | Const(s,_) => s
1.384 - | Free (s,_) => s
1.385 - | Var((s,_),_) => s
1.386 - | _ => error "string_of_term";
1.387 -
1.388 -val string_of_cterm = string_of_term o term_of;
1.389 -
1.390 -fun string_of_varpow x k =
1.391 - if k = 1 then string_of_cterm x
1.392 - else string_of_cterm x^"^"^string_of_int k;
1.393 -
1.394 -fun string_of_monomial m =
1.395 - if Ctermfunc.is_undefined m then "1" else
1.396 - let val vps = fold_rev (fn (x,k) => fn a => string_of_varpow x k :: a)
1.397 - (sort humanorder_varpow (Ctermfunc.graph m)) []
1.398 - in fold1 (fn s => fn t => s^"*"^t) vps
1.399 - end;
1.400 -
1.401 -fun string_of_cmonomial (c,m) =
1.402 - if Ctermfunc.is_undefined m then Rat.string_of_rat c
1.403 - else if c =/ rat_1 then string_of_monomial m
1.404 - else Rat.string_of_rat c ^ "*" ^ string_of_monomial m;;
1.405 -
1.406 -fun string_of_poly (p:poly) =
1.407 - if Monomialfunc.is_undefined p then "<<0>>" else
1.408 - let
1.409 - val cms = sort (fn ((m1,_),(m2,_)) => humanorder_monomial m1 m2) (Monomialfunc.graph p)
1.410 - val s = fold (fn (m,c) => fn a =>
1.411 - if c </ rat_0 then a ^ " - " ^ string_of_cmonomial(Rat.neg c,m)
1.412 - else a ^ " + " ^ string_of_cmonomial(c,m))
1.413 - cms ""
1.414 - val s1 = String.substring (s, 0, 3)
1.415 - val s2 = String.substring (s, 3, String.size s - 3)
1.416 - in "<<" ^(if s1 = " + " then s2 else "-"^s2)^">>"
1.417 - end;
1.418 -
1.419 -(* Conversion from HOL term. *)
1.420 -
1.421 -local
1.422 - val neg_tm = @{cterm "uminus :: real => _"}
1.423 - val add_tm = @{cterm "op + :: real => _"}
1.424 - val sub_tm = @{cterm "op - :: real => _"}
1.425 - val mul_tm = @{cterm "op * :: real => _"}
1.426 - val inv_tm = @{cterm "inverse :: real => _"}
1.427 - val div_tm = @{cterm "op / :: real => _"}
1.428 - val pow_tm = @{cterm "op ^ :: real => _"}
1.429 - val zero_tm = @{cterm "0:: real"}
1.430 - val is_numeral = can (HOLogic.dest_number o term_of)
1.431 - fun is_comb t = case t of _$_ => true | _ => false
1.432 - fun poly_of_term tm =
1.433 - if tm aconvc zero_tm then poly_0
1.434 - else if RealArith.is_ratconst tm
1.435 - then poly_const(RealArith.dest_ratconst tm)
1.436 - else
1.437 - (let val (lop,r) = Thm.dest_comb tm
1.438 - in if lop aconvc neg_tm then poly_neg(poly_of_term r)
1.439 - else if lop aconvc inv_tm then
1.440 - let val p = poly_of_term r
1.441 - in if poly_isconst p
1.442 - then poly_const(Rat.inv (eval Ctermfunc.undefined p))
1.443 - else error "poly_of_term: inverse of non-constant polyomial"
1.444 - end
1.445 - else (let val (opr,l) = Thm.dest_comb lop
1.446 - in
1.447 - if opr aconvc pow_tm andalso is_numeral r
1.448 - then poly_pow (poly_of_term l) ((snd o HOLogic.dest_number o term_of) r)
1.449 - else if opr aconvc add_tm
1.450 - then poly_add (poly_of_term l) (poly_of_term r)
1.451 - else if opr aconvc sub_tm
1.452 - then poly_sub (poly_of_term l) (poly_of_term r)
1.453 - else if opr aconvc mul_tm
1.454 - then poly_mul (poly_of_term l) (poly_of_term r)
1.455 - else if opr aconvc div_tm
1.456 - then let
1.457 - val p = poly_of_term l
1.458 - val q = poly_of_term r
1.459 - in if poly_isconst q then poly_cmul (Rat.inv (eval Ctermfunc.undefined q)) p
1.460 - else error "poly_of_term: division by non-constant polynomial"
1.461 - end
1.462 - else poly_var tm
1.463 -
1.464 - end
1.465 - handle CTERM ("dest_comb",_) => poly_var tm)
1.466 - end
1.467 - handle CTERM ("dest_comb",_) => poly_var tm)
1.468 -in
1.469 -val poly_of_term = fn tm =>
1.470 - if type_of (term_of tm) = @{typ real} then poly_of_term tm
1.471 - else error "poly_of_term: term does not have real type"
1.472 -end;
1.473 -
1.474 -(* String of vector (just a list of space-separated numbers). *)
1.475 -
1.476 -fun sdpa_of_vector (v:vector) =
1.477 - let
1.478 - val n = dim v
1.479 - val strs = map (decimalize 20 o (fn i => Intfunc.tryapplyd (snd v) i rat_0)) (1 upto n)
1.480 - in fold1 (fn x => fn y => x ^ " " ^ y) strs ^ "\n"
1.481 - end;
1.482 -
1.483 -fun increasing f ord (x,y) = ord (f x, f y);
1.484 -fun triple_int_ord ((a,b,c),(a',b',c')) =
1.485 - prod_ord int_ord (prod_ord int_ord int_ord)
1.486 - ((a,(b,c)),(a',(b',c')));
1.487 -structure Inttriplefunc = FuncFun(type key = int*int*int val ord = triple_int_ord);
1.488 -
1.489 -(* String for block diagonal matrix numbered k. *)
1.490 -
1.491 -fun sdpa_of_blockdiagonal k m =
1.492 - let
1.493 - val pfx = string_of_int k ^" "
1.494 - val ents =
1.495 - Inttriplefunc.fold (fn ((b,i,j), c) => fn a => if i > j then a else ((b,i,j),c)::a) m []
1.496 - val entss = sort (increasing fst triple_int_ord ) ents
1.497 - in fold_rev (fn ((b,i,j),c) => fn a =>
1.498 - pfx ^ string_of_int b ^ " " ^ string_of_int i ^ " " ^ string_of_int j ^
1.499 - " " ^ decimalize 20 c ^ "\n" ^ a) entss ""
1.500 - end;
1.501 -
1.502 -(* String for a matrix numbered k, in SDPA sparse format. *)
1.503 -
1.504 -fun sdpa_of_matrix k (m:matrix) =
1.505 - let
1.506 - val pfx = string_of_int k ^ " 1 "
1.507 - val ms = Intpairfunc.fold (fn ((i,j), c) => fn a => if i > j then a else ((i,j),c)::a) (snd m) []
1.508 - val mss = sort (increasing fst (prod_ord int_ord int_ord)) ms
1.509 - in fold_rev (fn ((i,j),c) => fn a =>
1.510 - pfx ^ string_of_int i ^ " " ^ string_of_int j ^
1.511 - " " ^ decimalize 20 c ^ "\n" ^ a) mss ""
1.512 - end;;
1.513 -
1.514 -(* ------------------------------------------------------------------------- *)
1.515 -(* String in SDPA sparse format for standard SDP problem: *)
1.516 -(* *)
1.517 -(* X = v_1 * [M_1] + ... + v_m * [M_m] - [M_0] must be PSD *)
1.518 -(* Minimize obj_1 * v_1 + ... obj_m * v_m *)
1.519 -(* ------------------------------------------------------------------------- *)
1.520 -
1.521 -fun sdpa_of_problem obj mats =
1.522 - let
1.523 - val m = length mats - 1
1.524 - val (n,_) = dimensions (hd mats)
1.525 - in
1.526 - string_of_int m ^ "\n" ^
1.527 - "1\n" ^
1.528 - string_of_int n ^ "\n" ^
1.529 - sdpa_of_vector obj ^
1.530 - fold_rev2 (fn k => fn m => fn a => sdpa_of_matrix (k - 1) m ^ a) (1 upto length mats) mats ""
1.531 - end;
1.532 -
1.533 -fun index_char str chr pos =
1.534 - if pos >= String.size str then ~1
1.535 - else if String.sub(str,pos) = chr then pos
1.536 - else index_char str chr (pos + 1);
1.537 -fun rat_of_quotient (a,b) = if b = 0 then rat_0 else Rat.rat_of_quotient (a,b);
1.538 -fun rat_of_string s =
1.539 - let val n = index_char s #"/" 0 in
1.540 - if n = ~1 then s |> IntInf.fromString |> valOf |> Rat.rat_of_int
1.541 - else
1.542 - let val SOME numer = IntInf.fromString(String.substring(s,0,n))
1.543 - val SOME den = IntInf.fromString (String.substring(s,n+1,String.size s - n - 1))
1.544 - in rat_of_quotient(numer, den)
1.545 - end
1.546 - end;
1.547 -
1.548 -fun isspace x = x = " " ;
1.549 -fun isnum x = x mem_string ["0","1","2","3","4","5","6","7","8","9"]
1.550 -
1.551 -(* More parser basics. *)
1.552 -
1.553 -local
1.554 - open Scan
1.555 -in
1.556 - val word = this_string
1.557 - fun token s =
1.558 - repeat ($$ " ") |-- word s --| repeat ($$ " ")
1.559 - val numeral = one isnum
1.560 - val decimalint = bulk numeral >> (rat_of_string o implode)
1.561 - val decimalfrac = bulk numeral
1.562 - >> (fn s => rat_of_string(implode s) // pow10 (length s))
1.563 - val decimalsig =
1.564 - decimalint -- option (Scan.$$ "." |-- decimalfrac)
1.565 - >> (fn (h,NONE) => h | (h,SOME x) => h +/ x)
1.566 - fun signed prs =
1.567 - $$ "-" |-- prs >> Rat.neg
1.568 - || $$ "+" |-- prs
1.569 - || prs;
1.570 -
1.571 -fun emptyin def xs = if null xs then (def,xs) else Scan.fail xs
1.572 -
1.573 - val exponent = ($$ "e" || $$ "E") |-- signed decimalint;
1.574 -
1.575 - val decimal = signed decimalsig -- (emptyin rat_0|| exponent)
1.576 - >> (fn (h, x) => h */ pow10 (int_of_rat x));
1.577 -end;
1.578 -
1.579 - fun mkparser p s =
1.580 - let val (x,rst) = p (explode s)
1.581 - in if null rst then x
1.582 - else error "mkparser: unparsed input"
1.583 - end;;
1.584 -val parse_decimal = mkparser decimal;
1.585 -
1.586 -fun fix err prs =
1.587 - prs || (fn x=> error err);
1.588 -
1.589 -fun listof prs sep err =
1.590 - prs -- Scan.bulk (sep |-- fix err prs) >> uncurry cons;
1.591 -
1.592 -(* Parse back a vector. *)
1.593 -
1.594 - val vector =
1.595 - token "{" |-- listof decimal (token ",") "decimal" --| token "}"
1.596 - >> vector_of_list
1.597 - val parse_vector = mkparser vector
1.598 - fun skipupto dscr prs inp =
1.599 - (dscr |-- prs
1.600 - || Scan.one (K true) |-- skipupto dscr prs) inp
1.601 - fun ignore inp = ((),[])
1.602 - fun sdpaoutput inp = skipupto (word "xVec" -- token "=")
1.603 - (vector --| ignore) inp
1.604 - fun csdpoutput inp = ((decimal -- Scan.bulk (Scan.$$ " " |-- Scan.option decimal) >> (fn (h,to) => map_filter I ((SOME h)::to))) --| ignore >> vector_of_list) inp
1.605 - val parse_sdpaoutput = mkparser sdpaoutput
1.606 - val parse_csdpoutput = mkparser csdpoutput
1.607 -
1.608 -(* Run prover on a problem in linear form. *)
1.609 -
1.610 -fun run_problem prover obj mats =
1.611 - parse_csdpoutput (prover (sdpa_of_problem obj mats))
1.612 -
1.613 -(*
1.614 -UNUSED
1.615 -
1.616 -(* Also parse the SDPA output to test success (CSDP yields a return code). *)
1.617 -
1.618 -local
1.619 - val prs =
1.620 - skipupto (word "phase.value" -- token "=")
1.621 - (Scan.option (Scan.$$ "p") -- Scan.option (Scan.$$ "d")
1.622 - -- (word "OPT" || word "FEAS"))
1.623 -in
1.624 - fun sdpa_run_succeeded s =
1.625 - (prs (explode s); true) handle _ => false
1.626 -end;
1.627 -
1.628 -(* The default parameters. Unfortunately this goes to a fixed file. *)
1.629 -
1.630 -val sdpa_default_parameters =
1.631 -"100 unsigned int maxIteration; \n1.0E-7 double 0.0 < epsilonStar;\n1.0E2 double 0.0 < lambdaStar;\n2.0 double 1.0 < omegaStar;\n-1.0E5 double lowerBound;\n1.0E5 double upperBound;\n0.1 double 0.0 <= betaStar < 1.0;\n0.2 double 0.0 <= betaBar < 1.0, betaStar <= betaBar;\n0.9 double 0.0 < gammaStar < 1.0;\n1.0E-7 double 0.0 < epsilonDash;\n";;
1.632 -
1.633 -(* These were suggested by Makoto Yamashita for problems where we are *)
1.634 -(* right at the edge of the semidefinite cone, as sometimes happens. *)
1.635 -
1.636 -val sdpa_alt_parameters =
1.637 -"1000 unsigned int maxIteration;\n1.0E-7 double 0.0 < epsilonStar;\n1.0E4 double 0.0 < lambdaStar;\n2.0 double 1.0 < omegaStar;\n-1.0E5 double lowerBound;\n1.0E5 double upperBound;\n0.1 double 0.0 <= betaStar < 1.0;\n0.2 double 0.0 <= betaBar < 1.0, betaStar <= betaBar;\n0.9 double 0.0 < gammaStar < 1.0;\n1.0E-7 double 0.0 < epsilonDash;\n";;
1.638 -
1.639 -val sdpa_params = sdpa_alt_parameters;;
1.640 -
1.641 -(* CSDP parameters; so far I'm sticking with the defaults. *)
1.642 -
1.643 -val csdp_default_parameters =
1.644 -"axtol=1.0e-8\natytol=1.0e-8\nobjtol=1.0e-8\npinftol=1.0e8\ndinftol=1.0e8\nmaxiter=100\nminstepfrac=0.9\nmaxstepfrac=0.97\nminstepp=1.0e-8\nminstepd=1.0e-8\nusexzgap=1\ntweakgap=0\naffine=0\nprintlevel=1\n";;
1.645 -
1.646 -val csdp_params = csdp_default_parameters;;
1.647 -
1.648 -fun tmp_file pre suf =
1.649 - let val i = string_of_int (round (random()))
1.650 - val name = Path.append (Path.variable "ISABELLE_TMP") (Path.explode (pre ^ i ^ suf))
1.651 - in
1.652 - if File.exists name then tmp_file pre suf
1.653 - else name
1.654 - end;
1.655 -
1.656 -(* Now call SDPA on a problem and parse back the output. *)
1.657 -
1.658 -fun run_sdpa dbg obj mats =
1.659 - let
1.660 - val input_file = tmp_file "sos" ".dat-s"
1.661 - val output_file = tmp_file "sos" ".out"
1.662 - val params_file = tmp_file "param" ".sdpa"
1.663 - val current_dir = File.pwd()
1.664 - val _ = File.write input_file
1.665 - (sdpa_of_problem "" obj mats)
1.666 - val _ = File.write params_file sdpa_params
1.667 - val _ = File.cd (Path.variable "ISABELLE_TMP")
1.668 - val _ = File.system_command ("sdpa "^ (Path.implode input_file) ^ " " ^
1.669 - (Path.implode output_file) ^
1.670 - (if dbg then "" else "> /dev/null"))
1.671 - val opr = File.read output_file
1.672 - in if not(sdpa_run_succeeded opr) then error "sdpa: call failed"
1.673 - else
1.674 - let val res = parse_sdpaoutput opr
1.675 - in ((if dbg then ()
1.676 - else (File.rm input_file; File.rm output_file ; File.cd current_dir));
1.677 - res)
1.678 - end
1.679 - end;
1.680 -
1.681 -fun sdpa obj mats = run_sdpa (!debugging) obj mats;
1.682 -
1.683 -(* The same thing with CSDP. *)
1.684 -
1.685 -fun run_csdp dbg obj mats =
1.686 - let
1.687 - val input_file = tmp_file "sos" ".dat-s"
1.688 - val output_file = tmp_file "sos" ".out"
1.689 - val params_file = tmp_file "param" ".csdp"
1.690 - val current_dir = File.pwd()
1.691 - val _ = File.write input_file (sdpa_of_problem "" obj mats)
1.692 - val _ = File.write params_file csdp_params
1.693 - val _ = File.cd (Path.variable "ISABELLE_TMP")
1.694 - val rv = system ("csdp "^(Path.implode input_file) ^ " "
1.695 - ^ (Path.implode output_file) ^
1.696 - (if dbg then "" else "> /dev/null"))
1.697 - val opr = File.read output_file
1.698 - val res = parse_csdpoutput opr
1.699 - in
1.700 - ((if dbg then ()
1.701 - else (File.rm input_file; File.rm output_file ; File.cd current_dir));
1.702 - (rv,res))
1.703 - end;
1.704 -
1.705 -fun csdp obj mats =
1.706 - let
1.707 - val (rv,res) = run_csdp (!debugging) obj mats
1.708 - in
1.709 - ((if rv = 1 orelse rv = 2 then error "csdp: Problem is infeasible"
1.710 - else if rv = 3 then writeln "csdp warning: Reduced accuracy"
1.711 - else if rv <> 0 then error ("csdp: error "^string_of_int rv)
1.712 - else ());
1.713 - res)
1.714 - end;
1.715 -
1.716 -*)
1.717 -
1.718 -(* Try some apparently sensible scaling first. Note that this is purely to *)
1.719 -(* get a cleaner translation to floating-point, and doesn't affect any of *)
1.720 -(* the results, in principle. In practice it seems a lot better when there *)
1.721 -(* are extreme numbers in the original problem. *)
1.722 -
1.723 - (* Version for (int*int) keys *)
1.724 -local
1.725 - fun max_rat x y = if x </ y then y else x
1.726 - fun common_denominator fld amat acc =
1.727 - fld (fn (m,c) => fn a => lcm_rat (denominator_rat c) a) amat acc
1.728 - fun maximal_element fld amat acc =
1.729 - fld (fn (m,c) => fn maxa => max_rat maxa (abs_rat c)) amat acc
1.730 -fun float_of_rat x = let val (a,b) = Rat.quotient_of_rat x
1.731 - in Real.fromLargeInt a / Real.fromLargeInt b end;
1.732 -in
1.733 -
1.734 -fun pi_scale_then solver (obj:vector) mats =
1.735 - let
1.736 - val cd1 = fold_rev (common_denominator Intpairfunc.fold) mats (rat_1)
1.737 - val cd2 = common_denominator Intfunc.fold (snd obj) (rat_1)
1.738 - val mats' = map (Intpairfunc.mapf (fn x => cd1 */ x)) mats
1.739 - val obj' = vector_cmul cd2 obj
1.740 - val max1 = fold_rev (maximal_element Intpairfunc.fold) mats' (rat_0)
1.741 - val max2 = maximal_element Intfunc.fold (snd obj') (rat_0)
1.742 - val scal1 = pow2 (20 - trunc(Math.ln (float_of_rat max1) / Math.ln 2.0))
1.743 - val scal2 = pow2 (20 - trunc(Math.ln (float_of_rat max2) / Math.ln 2.0))
1.744 - val mats'' = map (Intpairfunc.mapf (fn x => x */ scal1)) mats'
1.745 - val obj'' = vector_cmul scal2 obj'
1.746 - in solver obj'' mats''
1.747 - end
1.748 -end;
1.749 -
1.750 -(* Try some apparently sensible scaling first. Note that this is purely to *)
1.751 -(* get a cleaner translation to floating-point, and doesn't affect any of *)
1.752 -(* the results, in principle. In practice it seems a lot better when there *)
1.753 -(* are extreme numbers in the original problem. *)
1.754 -
1.755 - (* Version for (int*int*int) keys *)
1.756 -local
1.757 - fun max_rat x y = if x </ y then y else x
1.758 - fun common_denominator fld amat acc =
1.759 - fld (fn (m,c) => fn a => lcm_rat (denominator_rat c) a) amat acc
1.760 - fun maximal_element fld amat acc =
1.761 - fld (fn (m,c) => fn maxa => max_rat maxa (abs_rat c)) amat acc
1.762 -fun float_of_rat x = let val (a,b) = Rat.quotient_of_rat x
1.763 - in Real.fromLargeInt a / Real.fromLargeInt b end;
1.764 -fun int_of_float x = (trunc x handle Overflow => 0 | Domain => 0)
1.765 -in
1.766 -
1.767 -fun tri_scale_then solver (obj:vector) mats =
1.768 - let
1.769 - val cd1 = fold_rev (common_denominator Inttriplefunc.fold) mats (rat_1)
1.770 - val cd2 = common_denominator Intfunc.fold (snd obj) (rat_1)
1.771 - val mats' = map (Inttriplefunc.mapf (fn x => cd1 */ x)) mats
1.772 - val obj' = vector_cmul cd2 obj
1.773 - val max1 = fold_rev (maximal_element Inttriplefunc.fold) mats' (rat_0)
1.774 - val max2 = maximal_element Intfunc.fold (snd obj') (rat_0)
1.775 - val scal1 = pow2 (20 - int_of_float(Math.ln (float_of_rat max1) / Math.ln 2.0))
1.776 - val scal2 = pow2 (20 - int_of_float(Math.ln (float_of_rat max2) / Math.ln 2.0))
1.777 - val mats'' = map (Inttriplefunc.mapf (fn x => x */ scal1)) mats'
1.778 - val obj'' = vector_cmul scal2 obj'
1.779 - in solver obj'' mats''
1.780 - end
1.781 -end;
1.782 -
1.783 -(* Round a vector to "nice" rationals. *)
1.784 -
1.785 -fun nice_rational n x = round_rat (n */ x) // n;;
1.786 -fun nice_vector n ((d,v) : vector) =
1.787 - (d, Intfunc.fold (fn (i,c) => fn a =>
1.788 - let val y = nice_rational n c
1.789 - in if c =/ rat_0 then a
1.790 - else Intfunc.update (i,y) a end) v Intfunc.undefined):vector
1.791 -
1.792 -(*
1.793 -UNUSED
1.794 -
1.795 -(* Reduce linear program to SDP (diagonal matrices) and test with CSDP. This *)
1.796 -(* one tests A [-1;x1;..;xn] >= 0 (i.e. left column is negated constants). *)
1.797 -
1.798 -fun linear_program_basic a =
1.799 - let
1.800 - val (m,n) = dimensions a
1.801 - val mats = map (fn j => diagonal (column j a)) (1 upto n)
1.802 - val obj = vector_const rat_1 m
1.803 - val (rv,res) = run_csdp false obj mats
1.804 - in if rv = 1 orelse rv = 2 then false
1.805 - else if rv = 0 then true
1.806 - else error "linear_program: An error occurred in the SDP solver"
1.807 - end;
1.808 -
1.809 -(* Alternative interface testing A x >= b for matrix A, vector b. *)
1.810 -
1.811 -fun linear_program a b =
1.812 - let val (m,n) = dimensions a
1.813 - in if dim b <> m then error "linear_program: incompatible dimensions"
1.814 - else
1.815 - let
1.816 - val mats = diagonal b :: map (fn j => diagonal (column j a)) (1 upto n)
1.817 - val obj = vector_const rat_1 m
1.818 - val (rv,res) = run_csdp false obj mats
1.819 - in if rv = 1 orelse rv = 2 then false
1.820 - else if rv = 0 then true
1.821 - else error "linear_program: An error occurred in the SDP solver"
1.822 - end
1.823 - end;
1.824 -
1.825 -(* Test whether a point is in the convex hull of others. Rather than use *)
1.826 -(* computational geometry, express as linear inequalities and call CSDP. *)
1.827 -(* This is a bit lazy of me, but it's easy and not such a bottleneck so far. *)
1.828 -
1.829 -fun in_convex_hull pts pt =
1.830 - let
1.831 - val pts1 = (1::pt) :: map (fn x => 1::x) pts
1.832 - val pts2 = map (fn p => map (fn x => ~x) p @ p) pts1
1.833 - val n = length pts + 1
1.834 - val v = 2 * (length pt + 1)
1.835 - val m = v + n - 1
1.836 - val mat = ((m,n),
1.837 - itern 1 pts2 (fn pts => fn j => itern 1 pts
1.838 - (fn x => fn i => Intpairfunc.update ((i,j), Rat.rat_of_int x)))
1.839 - (iter (1,n) (fn i => Intpairfunc.update((v + i,i+1), rat_1))
1.840 - Intpairfunc.undefined))
1.841 - in linear_program_basic mat
1.842 - end;
1.843 -
1.844 -(* Filter down a set of points to a minimal set with the same convex hull. *)
1.845 -
1.846 -local
1.847 - fun augment1 (m::ms) = if in_convex_hull ms m then ms else ms@[m]
1.848 - fun augment m ms = funpow 3 augment1 (m::ms)
1.849 -in
1.850 -fun minimal_convex_hull mons =
1.851 - let val mons' = fold_rev augment (tl mons) [hd mons]
1.852 - in funpow (length mons') augment1 mons'
1.853 - end
1.854 -end;
1.855 -
1.856 -*)
1.857 -
1.858 -fun dest_ord f x = is_equal (f x);
1.859 -
1.860 -
1.861 -
1.862 -(* Stuff for "equations" ((int*int*int)->num functions). *)
1.863 -
1.864 -fun tri_equation_cmul c eq =
1.865 - if c =/ rat_0 then Inttriplefunc.undefined else Inttriplefunc.mapf (fn d => c */ d) eq;
1.866 -
1.867 -fun tri_equation_add eq1 eq2 = Inttriplefunc.combine (curry op +/) (fn x => x =/ rat_0) eq1 eq2;
1.868 -
1.869 -fun tri_equation_eval assig eq =
1.870 - let fun value v = Inttriplefunc.apply assig v
1.871 - in Inttriplefunc.fold (fn (v, c) => fn a => a +/ value v */ c) eq rat_0
1.872 - end;
1.873 -
1.874 -(* Eliminate among linear equations: return unconstrained variables and *)
1.875 -(* assignments for the others in terms of them. We give one pseudo-variable *)
1.876 -(* "one" that's used for a constant term. *)
1.877 -
1.878 -local
1.879 - fun extract_first p l = case l of (* FIXME : use find_first instead *)
1.880 - [] => error "extract_first"
1.881 - | h::t => if p h then (h,t) else
1.882 - let val (k,s) = extract_first p t in (k,h::s) end
1.883 -fun eliminate vars dun eqs = case vars of
1.884 - [] => if forall Inttriplefunc.is_undefined eqs then dun
1.885 - else raise Unsolvable
1.886 - | v::vs =>
1.887 - ((let
1.888 - val (eq,oeqs) = extract_first (fn e => Inttriplefunc.defined e v) eqs
1.889 - val a = Inttriplefunc.apply eq v
1.890 - val eq' = tri_equation_cmul ((Rat.neg rat_1) // a) (Inttriplefunc.undefine v eq)
1.891 - fun elim e =
1.892 - let val b = Inttriplefunc.tryapplyd e v rat_0
1.893 - in if b =/ rat_0 then e else
1.894 - tri_equation_add e (tri_equation_cmul (Rat.neg b // a) eq)
1.895 - end
1.896 - in eliminate vs (Inttriplefunc.update (v,eq') (Inttriplefunc.mapf elim dun)) (map elim oeqs)
1.897 - end)
1.898 - handle ERROR _ => eliminate vs dun eqs)
1.899 -in
1.900 -fun tri_eliminate_equations one vars eqs =
1.901 - let
1.902 - val assig = eliminate vars Inttriplefunc.undefined eqs
1.903 - val vs = Inttriplefunc.fold (fn (x, f) => fn a => remove (dest_ord triple_int_ord) one (Inttriplefunc.dom f) @ a) assig []
1.904 - in (distinct (dest_ord triple_int_ord) vs, assig)
1.905 - end
1.906 -end;
1.907 -
1.908 -(* Eliminate all variables, in an essentially arbitrary order. *)
1.909 -
1.910 -fun tri_eliminate_all_equations one =
1.911 - let
1.912 - fun choose_variable eq =
1.913 - let val (v,_) = Inttriplefunc.choose eq
1.914 - in if is_equal (triple_int_ord(v,one)) then
1.915 - let val eq' = Inttriplefunc.undefine v eq
1.916 - in if Inttriplefunc.is_undefined eq' then error "choose_variable"
1.917 - else fst (Inttriplefunc.choose eq')
1.918 - end
1.919 - else v
1.920 - end
1.921 - fun eliminate dun eqs = case eqs of
1.922 - [] => dun
1.923 - | eq::oeqs =>
1.924 - if Inttriplefunc.is_undefined eq then eliminate dun oeqs else
1.925 - let val v = choose_variable eq
1.926 - val a = Inttriplefunc.apply eq v
1.927 - val eq' = tri_equation_cmul ((Rat.rat_of_int ~1) // a)
1.928 - (Inttriplefunc.undefine v eq)
1.929 - fun elim e =
1.930 - let val b = Inttriplefunc.tryapplyd e v rat_0
1.931 - in if b =/ rat_0 then e
1.932 - else tri_equation_add e (tri_equation_cmul (Rat.neg b // a) eq)
1.933 - end
1.934 - in eliminate (Inttriplefunc.update(v, eq') (Inttriplefunc.mapf elim dun))
1.935 - (map elim oeqs)
1.936 - end
1.937 -in fn eqs =>
1.938 - let
1.939 - val assig = eliminate Inttriplefunc.undefined eqs
1.940 - val vs = Inttriplefunc.fold (fn (x, f) => fn a => remove (dest_ord triple_int_ord) one (Inttriplefunc.dom f) @ a) assig []
1.941 - in (distinct (dest_ord triple_int_ord) vs,assig)
1.942 - end
1.943 -end;
1.944 -
1.945 -(* Solve equations by assigning arbitrary numbers. *)
1.946 -
1.947 -fun tri_solve_equations one eqs =
1.948 - let
1.949 - val (vars,assigs) = tri_eliminate_all_equations one eqs
1.950 - val vfn = fold_rev (fn v => Inttriplefunc.update(v,rat_0)) vars
1.951 - (Inttriplefunc.onefunc(one, Rat.rat_of_int ~1))
1.952 - val ass =
1.953 - Inttriplefunc.combine (curry op +/) (K false)
1.954 - (Inttriplefunc.mapf (tri_equation_eval vfn) assigs) vfn
1.955 - in if forall (fn e => tri_equation_eval ass e =/ rat_0) eqs
1.956 - then Inttriplefunc.undefine one ass else raise Sanity
1.957 - end;
1.958 -
1.959 -(* Multiply equation-parametrized poly by regular poly and add accumulator. *)
1.960 -
1.961 -fun tri_epoly_pmul p q acc =
1.962 - Monomialfunc.fold (fn (m1, c) => fn a =>
1.963 - Monomialfunc.fold (fn (m2,e) => fn b =>
1.964 - let val m = monomial_mul m1 m2
1.965 - val es = Monomialfunc.tryapplyd b m Inttriplefunc.undefined
1.966 - in Monomialfunc.update (m,tri_equation_add (tri_equation_cmul c e) es) b
1.967 - end) q a) p acc ;
1.968 -
1.969 -(* Usual operations on equation-parametrized poly. *)
1.970 -
1.971 -fun tri_epoly_cmul c l =
1.972 - if c =/ rat_0 then Inttriplefunc.undefined else Inttriplefunc.mapf (tri_equation_cmul c) l;;
1.973 -
1.974 -val tri_epoly_neg = tri_epoly_cmul (Rat.rat_of_int ~1);
1.975 -
1.976 -val tri_epoly_add = Inttriplefunc.combine tri_equation_add Inttriplefunc.is_undefined;
1.977 -
1.978 -fun tri_epoly_sub p q = tri_epoly_add p (tri_epoly_neg q);;
1.979 -
1.980 -(* Stuff for "equations" ((int*int)->num functions). *)
1.981 -
1.982 -fun pi_equation_cmul c eq =
1.983 - if c =/ rat_0 then Inttriplefunc.undefined else Inttriplefunc.mapf (fn d => c */ d) eq;
1.984 -
1.985 -fun pi_equation_add eq1 eq2 = Inttriplefunc.combine (curry op +/) (fn x => x =/ rat_0) eq1 eq2;
1.986 -
1.987 -fun pi_equation_eval assig eq =
1.988 - let fun value v = Inttriplefunc.apply assig v
1.989 - in Inttriplefunc.fold (fn (v, c) => fn a => a +/ value v */ c) eq rat_0
1.990 - end;
1.991 -
1.992 -(* Eliminate among linear equations: return unconstrained variables and *)
1.993 -(* assignments for the others in terms of them. We give one pseudo-variable *)
1.994 -(* "one" that's used for a constant term. *)
1.995 -
1.996 -local
1.997 -fun extract_first p l = case l of
1.998 - [] => error "extract_first"
1.999 - | h::t => if p h then (h,t) else
1.1000 - let val (k,s) = extract_first p t in (k,h::s) end
1.1001 -fun eliminate vars dun eqs = case vars of
1.1002 - [] => if forall Inttriplefunc.is_undefined eqs then dun
1.1003 - else raise Unsolvable
1.1004 - | v::vs =>
1.1005 - let
1.1006 - val (eq,oeqs) = extract_first (fn e => Inttriplefunc.defined e v) eqs
1.1007 - val a = Inttriplefunc.apply eq v
1.1008 - val eq' = pi_equation_cmul ((Rat.neg rat_1) // a) (Inttriplefunc.undefine v eq)
1.1009 - fun elim e =
1.1010 - let val b = Inttriplefunc.tryapplyd e v rat_0
1.1011 - in if b =/ rat_0 then e else
1.1012 - pi_equation_add e (pi_equation_cmul (Rat.neg b // a) eq)
1.1013 - end
1.1014 - in eliminate vs (Inttriplefunc.update (v,eq') (Inttriplefunc.mapf elim dun)) (map elim oeqs)
1.1015 - end
1.1016 - handle ERROR _ => eliminate vs dun eqs
1.1017 -in
1.1018 -fun pi_eliminate_equations one vars eqs =
1.1019 - let
1.1020 - val assig = eliminate vars Inttriplefunc.undefined eqs
1.1021 - val vs = Inttriplefunc.fold (fn (x, f) => fn a => remove (dest_ord triple_int_ord) one (Inttriplefunc.dom f) @ a) assig []
1.1022 - in (distinct (dest_ord triple_int_ord) vs, assig)
1.1023 - end
1.1024 -end;
1.1025 -
1.1026 -(* Eliminate all variables, in an essentially arbitrary order. *)
1.1027 -
1.1028 -fun pi_eliminate_all_equations one =
1.1029 - let
1.1030 - fun choose_variable eq =
1.1031 - let val (v,_) = Inttriplefunc.choose eq
1.1032 - in if is_equal (triple_int_ord(v,one)) then
1.1033 - let val eq' = Inttriplefunc.undefine v eq
1.1034 - in if Inttriplefunc.is_undefined eq' then error "choose_variable"
1.1035 - else fst (Inttriplefunc.choose eq')
1.1036 - end
1.1037 - else v
1.1038 - end
1.1039 - fun eliminate dun eqs = case eqs of
1.1040 - [] => dun
1.1041 - | eq::oeqs =>
1.1042 - if Inttriplefunc.is_undefined eq then eliminate dun oeqs else
1.1043 - let val v = choose_variable eq
1.1044 - val a = Inttriplefunc.apply eq v
1.1045 - val eq' = pi_equation_cmul ((Rat.rat_of_int ~1) // a)
1.1046 - (Inttriplefunc.undefine v eq)
1.1047 - fun elim e =
1.1048 - let val b = Inttriplefunc.tryapplyd e v rat_0
1.1049 - in if b =/ rat_0 then e
1.1050 - else pi_equation_add e (pi_equation_cmul (Rat.neg b // a) eq)
1.1051 - end
1.1052 - in eliminate (Inttriplefunc.update(v, eq') (Inttriplefunc.mapf elim dun))
1.1053 - (map elim oeqs)
1.1054 - end
1.1055 -in fn eqs =>
1.1056 - let
1.1057 - val assig = eliminate Inttriplefunc.undefined eqs
1.1058 - val vs = Inttriplefunc.fold (fn (x, f) => fn a => remove (dest_ord triple_int_ord) one (Inttriplefunc.dom f) @ a) assig []
1.1059 - in (distinct (dest_ord triple_int_ord) vs,assig)
1.1060 - end
1.1061 -end;
1.1062 -
1.1063 -(* Solve equations by assigning arbitrary numbers. *)
1.1064 -
1.1065 -fun pi_solve_equations one eqs =
1.1066 - let
1.1067 - val (vars,assigs) = pi_eliminate_all_equations one eqs
1.1068 - val vfn = fold_rev (fn v => Inttriplefunc.update(v,rat_0)) vars
1.1069 - (Inttriplefunc.onefunc(one, Rat.rat_of_int ~1))
1.1070 - val ass =
1.1071 - Inttriplefunc.combine (curry op +/) (K false)
1.1072 - (Inttriplefunc.mapf (pi_equation_eval vfn) assigs) vfn
1.1073 - in if forall (fn e => pi_equation_eval ass e =/ rat_0) eqs
1.1074 - then Inttriplefunc.undefine one ass else raise Sanity
1.1075 - end;
1.1076 -
1.1077 -(* Multiply equation-parametrized poly by regular poly and add accumulator. *)
1.1078 -
1.1079 -fun pi_epoly_pmul p q acc =
1.1080 - Monomialfunc.fold (fn (m1, c) => fn a =>
1.1081 - Monomialfunc.fold (fn (m2,e) => fn b =>
1.1082 - let val m = monomial_mul m1 m2
1.1083 - val es = Monomialfunc.tryapplyd b m Inttriplefunc.undefined
1.1084 - in Monomialfunc.update (m,pi_equation_add (pi_equation_cmul c e) es) b
1.1085 - end) q a) p acc ;
1.1086 -
1.1087 -(* Usual operations on equation-parametrized poly. *)
1.1088 -
1.1089 -fun pi_epoly_cmul c l =
1.1090 - if c =/ rat_0 then Inttriplefunc.undefined else Inttriplefunc.mapf (pi_equation_cmul c) l;;
1.1091 -
1.1092 -val pi_epoly_neg = pi_epoly_cmul (Rat.rat_of_int ~1);
1.1093 -
1.1094 -val pi_epoly_add = Inttriplefunc.combine pi_equation_add Inttriplefunc.is_undefined;
1.1095 -
1.1096 -fun pi_epoly_sub p q = pi_epoly_add p (pi_epoly_neg q);;
1.1097 -
1.1098 -fun allpairs f l1 l2 = fold_rev (fn x => (curry (op @)) (map (f x) l2)) l1 [];
1.1099 -
1.1100 -(* Hence produce the "relevant" monomials: those whose squares lie in the *)
1.1101 -(* Newton polytope of the monomials in the input. (This is enough according *)
1.1102 -(* to Reznik: "Extremal PSD forms with few terms", Duke Math. Journal, *)
1.1103 -(* vol 45, pp. 363--374, 1978. *)
1.1104 -(* *)
1.1105 -(* These are ordered in sort of decreasing degree. In particular the *)
1.1106 -(* constant monomial is last; this gives an order in diagonalization of the *)
1.1107 -(* quadratic form that will tend to display constants. *)
1.1108 -
1.1109 -(*
1.1110 -UNUSED
1.1111 -
1.1112 -fun newton_polytope pol =
1.1113 - let
1.1114 - val vars = poly_variables pol
1.1115 - val mons = map (fn m => map (fn x => monomial_degree x m) vars)
1.1116 - (Monomialfunc.dom pol)
1.1117 - val ds = map (fn x => (degree x pol + 1) div 2) vars
1.1118 - val all = fold_rev (fn n => allpairs cons (0 upto n)) ds [[]]
1.1119 - val mons' = minimal_convex_hull mons
1.1120 - val all' =
1.1121 - filter (fn m => in_convex_hull mons' (map (fn x => 2 * x) m)) all
1.1122 - in map (fn m => fold_rev2 (fn v => fn i => fn a => if i = 0 then a else Ctermfunc.update (v,i) a)
1.1123 - vars m monomial_1) (rev all')
1.1124 - end;
1.1125 -
1.1126 -*)
1.1127 -
1.1128 -(* Diagonalize (Cholesky/LDU) the matrix corresponding to a quadratic form. *)
1.1129 -
1.1130 -local
1.1131 -fun diagonalize n i m =
1.1132 - if Intpairfunc.is_undefined (snd m) then []
1.1133 - else
1.1134 - let val a11 = Intpairfunc.tryapplyd (snd m) (i,i) rat_0
1.1135 - in if a11 </ rat_0 then error "diagonalize: not PSD"
1.1136 - else if a11 =/ rat_0 then
1.1137 - if Intfunc.is_undefined (snd (row i m)) then diagonalize n (i + 1) m
1.1138 - else error "diagonalize: not PSD ___ "
1.1139 - else
1.1140 - let
1.1141 - val v = row i m
1.1142 - val v' = (fst v, Intfunc.fold (fn (i, c) => fn a =>
1.1143 - let val y = c // a11
1.1144 - in if y = rat_0 then a else Intfunc.update (i,y) a
1.1145 - end) (snd v) Intfunc.undefined)
1.1146 - fun upt0 x y a = if y = rat_0 then a else Intpairfunc.update (x,y) a
1.1147 - val m' =
1.1148 - ((n,n),
1.1149 - iter (i+1,n) (fn j =>
1.1150 - iter (i+1,n) (fn k =>
1.1151 - (upt0 (j,k) (Intpairfunc.tryapplyd (snd m) (j,k) rat_0 -/ Intfunc.tryapplyd (snd v) j rat_0 */ Intfunc.tryapplyd (snd v') k rat_0))))
1.1152 - Intpairfunc.undefined)
1.1153 - in (a11,v')::diagonalize n (i + 1) m'
1.1154 - end
1.1155 - end
1.1156 -in
1.1157 -fun diag m =
1.1158 - let
1.1159 - val nn = dimensions m
1.1160 - val n = fst nn
1.1161 - in if snd nn <> n then error "diagonalize: non-square matrix"
1.1162 - else diagonalize n 1 m
1.1163 - end
1.1164 -end;
1.1165 -
1.1166 -fun gcd_rat a b = Rat.rat_of_int (Integer.gcd (int_of_rat a) (int_of_rat b));
1.1167 -
1.1168 -(* Adjust a diagonalization to collect rationals at the start. *)
1.1169 - (* FIXME : Potentially polymorphic keys, but here only: integers!! *)
1.1170 -local
1.1171 - fun upd0 x y a = if y =/ rat_0 then a else Intfunc.update(x,y) a;
1.1172 - fun mapa f (d,v) =
1.1173 - (d, Intfunc.fold (fn (i,c) => fn a => upd0 i (f c) a) v Intfunc.undefined)
1.1174 - fun adj (c,l) =
1.1175 - let val a =
1.1176 - Intfunc.fold (fn (i,c) => fn a => lcm_rat a (denominator_rat c))
1.1177 - (snd l) rat_1 //
1.1178 - Intfunc.fold (fn (i,c) => fn a => gcd_rat a (numerator_rat c))
1.1179 - (snd l) rat_0
1.1180 - in ((c // (a */ a)),mapa (fn x => a */ x) l)
1.1181 - end
1.1182 -in
1.1183 -fun deration d = if null d then (rat_0,d) else
1.1184 - let val d' = map adj d
1.1185 - val a = fold (lcm_rat o denominator_rat o fst) d' rat_1 //
1.1186 - fold (gcd_rat o numerator_rat o fst) d' rat_0
1.1187 - in ((rat_1 // a),map (fn (c,l) => (a */ c,l)) d')
1.1188 - end
1.1189 -end;
1.1190 -
1.1191 -(* Enumeration of monomials with given multidegree bound. *)
1.1192 -
1.1193 -fun enumerate_monomials d vars =
1.1194 - if d < 0 then []
1.1195 - else if d = 0 then [Ctermfunc.undefined]
1.1196 - else if null vars then [monomial_1] else
1.1197 - let val alts =
1.1198 - map (fn k => let val oths = enumerate_monomials (d - k) (tl vars)
1.1199 - in map (fn ks => if k = 0 then ks else Ctermfunc.update (hd vars, k) ks) oths end) (0 upto d)
1.1200 - in fold1 (curry op @) alts
1.1201 - end;
1.1202 -
1.1203 -(* Enumerate products of distinct input polys with degree <= d. *)
1.1204 -(* We ignore any constant input polynomials. *)
1.1205 -(* Give the output polynomial and a record of how it was derived. *)
1.1206 -
1.1207 -local
1.1208 - open RealArith
1.1209 -in
1.1210 -fun enumerate_products d pols =
1.1211 -if d = 0 then [(poly_const rat_1,Rational_lt rat_1)]
1.1212 -else if d < 0 then [] else
1.1213 -case pols of
1.1214 - [] => [(poly_const rat_1,Rational_lt rat_1)]
1.1215 - | (p,b)::ps =>
1.1216 - let val e = multidegree p
1.1217 - in if e = 0 then enumerate_products d ps else
1.1218 - enumerate_products d ps @
1.1219 - map (fn (q,c) => (poly_mul p q,Product(b,c)))
1.1220 - (enumerate_products (d - e) ps)
1.1221 - end
1.1222 -end;
1.1223 -
1.1224 -(* Convert regular polynomial. Note that we treat (0,0,0) as -1. *)
1.1225 -
1.1226 -fun epoly_of_poly p =
1.1227 - Monomialfunc.fold (fn (m,c) => fn a => Monomialfunc.update (m, Inttriplefunc.onefunc ((0,0,0), Rat.neg c)) a) p Monomialfunc.undefined;
1.1228 -
1.1229 -(* String for block diagonal matrix numbered k. *)
1.1230 -
1.1231 -fun sdpa_of_blockdiagonal k m =
1.1232 - let
1.1233 - val pfx = string_of_int k ^" "
1.1234 - val ents =
1.1235 - Inttriplefunc.fold
1.1236 - (fn ((b,i,j),c) => fn a => if i > j then a else ((b,i,j),c)::a)
1.1237 - m []
1.1238 - val entss = sort (increasing fst triple_int_ord) ents
1.1239 - in fold_rev (fn ((b,i,j),c) => fn a =>
1.1240 - pfx ^ string_of_int b ^ " " ^ string_of_int i ^ " " ^ string_of_int j ^
1.1241 - " " ^ decimalize 20 c ^ "\n" ^ a) entss ""
1.1242 - end;
1.1243 -
1.1244 -(* SDPA for problem using block diagonal (i.e. multiple SDPs) *)
1.1245 -
1.1246 -fun sdpa_of_blockproblem nblocks blocksizes obj mats =
1.1247 - let val m = length mats - 1
1.1248 - in
1.1249 - string_of_int m ^ "\n" ^
1.1250 - string_of_int nblocks ^ "\n" ^
1.1251 - (fold1 (fn s => fn t => s^" "^t) (map string_of_int blocksizes)) ^
1.1252 - "\n" ^
1.1253 - sdpa_of_vector obj ^
1.1254 - fold_rev2 (fn k => fn m => fn a => sdpa_of_blockdiagonal (k - 1) m ^ a)
1.1255 - (1 upto length mats) mats ""
1.1256 - end;
1.1257 -
1.1258 -(* Run prover on a problem in block diagonal form. *)
1.1259 -
1.1260 -fun run_blockproblem prover nblocks blocksizes obj mats=
1.1261 - parse_csdpoutput (prover (sdpa_of_blockproblem nblocks blocksizes obj mats))
1.1262 -
1.1263 -(*
1.1264 -UNUSED
1.1265 -
1.1266 -(* Hence run CSDP on a problem in block diagonal form. *)
1.1267 -
1.1268 -fun run_csdp dbg nblocks blocksizes obj mats =
1.1269 - let
1.1270 - val input_file = tmp_file "sos" ".dat-s"
1.1271 - val output_file = tmp_file "sos" ".out"
1.1272 - val params_file = tmp_file "param" ".csdp"
1.1273 - val _ = File.write input_file
1.1274 - (sdpa_of_blockproblem "" nblocks blocksizes obj mats)
1.1275 - val _ = File.write params_file csdp_params
1.1276 - val current_dir = File.pwd()
1.1277 - val _ = File.cd (Path.variable "ISABELLE_TMP")
1.1278 - val rv = system ("csdp "^(Path.implode input_file) ^ " "
1.1279 - ^ (Path.implode output_file) ^
1.1280 - (if dbg then "" else "> /dev/null"))
1.1281 - val opr = File.read output_file
1.1282 - val res = parse_csdpoutput opr
1.1283 - in
1.1284 - ((if dbg then ()
1.1285 - else (File.rm input_file ; File.rm output_file ; File.cd current_dir));
1.1286 - (rv,res))
1.1287 - end;
1.1288 -
1.1289 -fun csdp nblocks blocksizes obj mats =
1.1290 - let
1.1291 - val (rv,res) = run_csdp (!debugging) nblocks blocksizes obj mats
1.1292 - in ((if rv = 1 orelse rv = 2 then error "csdp: Problem is infeasible"
1.1293 - else if rv = 3 then writeln "csdp warning: Reduced accuracy"
1.1294 - else if rv <> 0 then error ("csdp: error "^string_of_int rv)
1.1295 - else ());
1.1296 - res)
1.1297 - end;
1.1298 -*)
1.1299 -
1.1300 -(* 3D versions of matrix operations to consider blocks separately. *)
1.1301 -
1.1302 -val bmatrix_add = Inttriplefunc.combine (curry op +/) (fn x => x =/ rat_0);
1.1303 -fun bmatrix_cmul c bm =
1.1304 - if c =/ rat_0 then Inttriplefunc.undefined
1.1305 - else Inttriplefunc.mapf (fn x => c */ x) bm;
1.1306 -
1.1307 -val bmatrix_neg = bmatrix_cmul (Rat.rat_of_int ~1);
1.1308 -fun bmatrix_sub m1 m2 = bmatrix_add m1 (bmatrix_neg m2);;
1.1309 -
1.1310 -(* Smash a block matrix into components. *)
1.1311 -
1.1312 -fun blocks blocksizes bm =
1.1313 - map (fn (bs,b0) =>
1.1314 - let val m = Inttriplefunc.fold
1.1315 - (fn ((b,i,j),c) => fn a => if b = b0 then Intpairfunc.update ((i,j),c) a else a) bm Intpairfunc.undefined
1.1316 - val d = Intpairfunc.fold (fn ((i,j),c) => fn a => max a (max i j)) m 0
1.1317 - in (((bs,bs),m):matrix) end)
1.1318 - (blocksizes ~~ (1 upto length blocksizes));;
1.1319 -
1.1320 -(* FIXME : Get rid of this !!!*)
1.1321 -local
1.1322 - fun tryfind_with msg f [] = error msg
1.1323 - | tryfind_with msg f (x::xs) = (f x handle ERROR s => tryfind_with s f xs);
1.1324 -in
1.1325 - fun tryfind f = tryfind_with "tryfind" f
1.1326 -end
1.1327 -
1.1328 -(*
1.1329 -fun tryfind f [] = error "tryfind"
1.1330 - | tryfind f (x::xs) = (f x handle ERROR _ => tryfind f xs);
1.1331 -*)
1.1332 -
1.1333 -(* Positiv- and Nullstellensatz. Flag "linf" forces a linear representation. *)
1.1334 -
1.1335 -
1.1336 -local
1.1337 - open RealArith
1.1338 -in
1.1339 -fun real_positivnullstellensatz_general prover linf d eqs leqs pol =
1.1340 -let
1.1341 - val vars = fold_rev (curry (gen_union (op aconvc)) o poly_variables)
1.1342 - (pol::eqs @ map fst leqs) []
1.1343 - val monoid = if linf then
1.1344 - (poly_const rat_1,Rational_lt rat_1)::
1.1345 - (filter (fn (p,c) => multidegree p <= d) leqs)
1.1346 - else enumerate_products d leqs
1.1347 - val nblocks = length monoid
1.1348 - fun mk_idmultiplier k p =
1.1349 - let
1.1350 - val e = d - multidegree p
1.1351 - val mons = enumerate_monomials e vars
1.1352 - val nons = mons ~~ (1 upto length mons)
1.1353 - in (mons,
1.1354 - fold_rev (fn (m,n) => Monomialfunc.update(m,Inttriplefunc.onefunc((~k,~n,n),rat_1))) nons Monomialfunc.undefined)
1.1355 - end
1.1356 -
1.1357 - fun mk_sqmultiplier k (p,c) =
1.1358 - let
1.1359 - val e = (d - multidegree p) div 2
1.1360 - val mons = enumerate_monomials e vars
1.1361 - val nons = mons ~~ (1 upto length mons)
1.1362 - in (mons,
1.1363 - fold_rev (fn (m1,n1) =>
1.1364 - fold_rev (fn (m2,n2) => fn a =>
1.1365 - let val m = monomial_mul m1 m2
1.1366 - in if n1 > n2 then a else
1.1367 - let val c = if n1 = n2 then rat_1 else rat_2
1.1368 - val e = Monomialfunc.tryapplyd a m Inttriplefunc.undefined
1.1369 - in Monomialfunc.update(m, tri_equation_add (Inttriplefunc.onefunc((k,n1,n2), c)) e) a
1.1370 - end
1.1371 - end) nons)
1.1372 - nons Monomialfunc.undefined)
1.1373 - end
1.1374 -
1.1375 - val (sqmonlist,sqs) = split_list (map2 mk_sqmultiplier (1 upto length monoid) monoid)
1.1376 - val (idmonlist,ids) = split_list(map2 mk_idmultiplier (1 upto length eqs) eqs)
1.1377 - val blocksizes = map length sqmonlist
1.1378 - val bigsum =
1.1379 - fold_rev2 (fn p => fn q => fn a => tri_epoly_pmul p q a) eqs ids
1.1380 - (fold_rev2 (fn (p,c) => fn s => fn a => tri_epoly_pmul p s a) monoid sqs
1.1381 - (epoly_of_poly(poly_neg pol)))
1.1382 - val eqns = Monomialfunc.fold (fn (m,e) => fn a => e::a) bigsum []
1.1383 - val (pvs,assig) = tri_eliminate_all_equations (0,0,0) eqns
1.1384 - val qvars = (0,0,0)::pvs
1.1385 - val allassig = fold_rev (fn v => Inttriplefunc.update(v,(Inttriplefunc.onefunc(v,rat_1)))) pvs assig
1.1386 - fun mk_matrix v =
1.1387 - Inttriplefunc.fold (fn ((b,i,j), ass) => fn m =>
1.1388 - if b < 0 then m else
1.1389 - let val c = Inttriplefunc.tryapplyd ass v rat_0
1.1390 - in if c = rat_0 then m else
1.1391 - Inttriplefunc.update ((b,j,i), c) (Inttriplefunc.update ((b,i,j), c) m)
1.1392 - end)
1.1393 - allassig Inttriplefunc.undefined
1.1394 - val diagents = Inttriplefunc.fold
1.1395 - (fn ((b,i,j), e) => fn a => if b > 0 andalso i = j then tri_equation_add e a else a)
1.1396 - allassig Inttriplefunc.undefined
1.1397 -
1.1398 - val mats = map mk_matrix qvars
1.1399 - val obj = (length pvs,
1.1400 - itern 1 pvs (fn v => fn i => Intfunc.updatep iszero (i,Inttriplefunc.tryapplyd diagents v rat_0))
1.1401 - Intfunc.undefined)
1.1402 - val raw_vec = if null pvs then vector_0 0
1.1403 - else tri_scale_then (run_blockproblem prover nblocks blocksizes) obj mats
1.1404 - fun int_element (d,v) i = Intfunc.tryapplyd v i rat_0
1.1405 - fun cterm_element (d,v) i = Ctermfunc.tryapplyd v i rat_0
1.1406 -
1.1407 - fun find_rounding d =
1.1408 - let
1.1409 - val _ = if !debugging
1.1410 - then writeln ("Trying rounding with limit "^Rat.string_of_rat d ^ "\n")
1.1411 - else ()
1.1412 - val vec = nice_vector d raw_vec
1.1413 - val blockmat = iter (1,dim vec)
1.1414 - (fn i => fn a => bmatrix_add (bmatrix_cmul (int_element vec i) (nth mats i)) a)
1.1415 - (bmatrix_neg (nth mats 0))
1.1416 - val allmats = blocks blocksizes blockmat
1.1417 - in (vec,map diag allmats)
1.1418 - end
1.1419 - val (vec,ratdias) =
1.1420 - if null pvs then find_rounding rat_1
1.1421 - else tryfind find_rounding (map Rat.rat_of_int (1 upto 31) @
1.1422 - map pow2 (5 upto 66))
1.1423 - val newassigs =
1.1424 - fold_rev (fn k => Inttriplefunc.update (nth pvs (k - 1), int_element vec k))
1.1425 - (1 upto dim vec) (Inttriplefunc.onefunc ((0,0,0), Rat.rat_of_int ~1))
1.1426 - val finalassigs =
1.1427 - Inttriplefunc.fold (fn (v,e) => fn a => Inttriplefunc.update(v, tri_equation_eval newassigs e) a) allassig newassigs
1.1428 - fun poly_of_epoly p =
1.1429 - Monomialfunc.fold (fn (v,e) => fn a => Monomialfunc.updatep iszero (v,tri_equation_eval finalassigs e) a)
1.1430 - p Monomialfunc.undefined
1.1431 - fun mk_sos mons =
1.1432 - let fun mk_sq (c,m) =
1.1433 - (c,fold_rev (fn k=> fn a => Monomialfunc.updatep iszero (nth mons (k - 1), int_element m k) a)
1.1434 - (1 upto length mons) Monomialfunc.undefined)
1.1435 - in map mk_sq
1.1436 - end
1.1437 - val sqs = map2 mk_sos sqmonlist ratdias
1.1438 - val cfs = map poly_of_epoly ids
1.1439 - val msq = filter (fn (a,b) => not (null b)) (map2 pair monoid sqs)
1.1440 - fun eval_sq sqs = fold_rev (fn (c,q) => poly_add (poly_cmul c (poly_mul q q))) sqs poly_0
1.1441 - val sanity =
1.1442 - fold_rev (fn ((p,c),s) => poly_add (poly_mul p (eval_sq s))) msq
1.1443 - (fold_rev2 (fn p => fn q => poly_add (poly_mul p q)) cfs eqs
1.1444 - (poly_neg pol))
1.1445 -
1.1446 -in if not(Monomialfunc.is_undefined sanity) then raise Sanity else
1.1447 - (cfs,map (fn (a,b) => (snd a,b)) msq)
1.1448 - end
1.1449 -
1.1450 -
1.1451 -end;
1.1452 -
1.1453 -(* Iterative deepening. *)
1.1454 -
1.1455 -fun deepen f n =
1.1456 - (writeln ("Searching with depth limit " ^ string_of_int n) ; (f n handle ERROR s => (writeln ("failed with message: " ^ s) ; deepen f (n+1))))
1.1457 -
1.1458 -(* The ordering so we can create canonical HOL polynomials. *)
1.1459 -
1.1460 -fun dest_monomial mon = sort (increasing fst cterm_ord) (Ctermfunc.graph mon);
1.1461 -
1.1462 -fun monomial_order (m1,m2) =
1.1463 - if Ctermfunc.is_undefined m2 then LESS
1.1464 - else if Ctermfunc.is_undefined m1 then GREATER
1.1465 - else
1.1466 - let val mon1 = dest_monomial m1
1.1467 - val mon2 = dest_monomial m2
1.1468 - val deg1 = fold (curry op + o snd) mon1 0
1.1469 - val deg2 = fold (curry op + o snd) mon2 0
1.1470 - in if deg1 < deg2 then GREATER else if deg1 > deg2 then LESS
1.1471 - else list_ord (prod_ord cterm_ord int_ord) (mon1,mon2)
1.1472 - end;
1.1473 -
1.1474 -fun dest_poly p =
1.1475 - map (fn (m,c) => (c,dest_monomial m))
1.1476 - (sort (prod_ord monomial_order (K EQUAL)) (Monomialfunc.graph p));
1.1477 -
1.1478 -(* Map back polynomials and their composites to HOL. *)
1.1479 -
1.1480 -local
1.1481 - open Thm Numeral RealArith
1.1482 -in
1.1483 -
1.1484 -fun cterm_of_varpow x k = if k = 1 then x else capply (capply @{cterm "op ^ :: real => _"} x)
1.1485 - (mk_cnumber @{ctyp nat} k)
1.1486 -
1.1487 -fun cterm_of_monomial m =
1.1488 - if Ctermfunc.is_undefined m then @{cterm "1::real"}
1.1489 - else
1.1490 - let
1.1491 - val m' = dest_monomial m
1.1492 - val vps = fold_rev (fn (x,k) => cons (cterm_of_varpow x k)) m' []
1.1493 - in fold1 (fn s => fn t => capply (capply @{cterm "op * :: real => _"} s) t) vps
1.1494 - end
1.1495 -
1.1496 -fun cterm_of_cmonomial (m,c) = if Ctermfunc.is_undefined m then cterm_of_rat c
1.1497 - else if c = Rat.one then cterm_of_monomial m
1.1498 - else capply (capply @{cterm "op *::real => _"} (cterm_of_rat c)) (cterm_of_monomial m);
1.1499 -
1.1500 -fun cterm_of_poly p =
1.1501 - if Monomialfunc.is_undefined p then @{cterm "0::real"}
1.1502 - else
1.1503 - let
1.1504 - val cms = map cterm_of_cmonomial
1.1505 - (sort (prod_ord monomial_order (K EQUAL)) (Monomialfunc.graph p))
1.1506 - in fold1 (fn t1 => fn t2 => capply(capply @{cterm "op + :: real => _"} t1) t2) cms
1.1507 - end;
1.1508 -
1.1509 -fun cterm_of_sqterm (c,p) = Product(Rational_lt c,Square(cterm_of_poly p));
1.1510 -
1.1511 -fun cterm_of_sos (pr,sqs) = if null sqs then pr
1.1512 - else Product(pr,fold1 (fn a => fn b => Sum(a,b)) (map cterm_of_sqterm sqs));
1.1513 -
1.1514 -end
1.1515 -
1.1516 -(* Interface to HOL. *)
1.1517 -local
1.1518 - open Thm Conv RealArith
1.1519 - val concl = dest_arg o cprop_of
1.1520 - fun simple_cterm_ord t u = TermOrd.fast_term_ord (term_of t, term_of u) = LESS
1.1521 -in
1.1522 - (* FIXME: Replace tryfind by get_first !! *)
1.1523 -fun real_nonlinear_prover prover ctxt =
1.1524 - let
1.1525 - val {add,mul,neg,pow,sub,main} = Normalizer.semiring_normalizers_ord_wrapper ctxt
1.1526 - (valOf (NormalizerData.match ctxt @{cterm "(0::real) + 1"}))
1.1527 - simple_cterm_ord
1.1528 - val (real_poly_add_conv,real_poly_mul_conv,real_poly_neg_conv,
1.1529 - real_poly_pow_conv,real_poly_sub_conv,real_poly_conv) = (add,mul,neg,pow,sub,main)
1.1530 - fun mainf translator (eqs,les,lts) =
1.1531 - let
1.1532 - val eq0 = map (poly_of_term o dest_arg1 o concl) eqs
1.1533 - val le0 = map (poly_of_term o dest_arg o concl) les
1.1534 - val lt0 = map (poly_of_term o dest_arg o concl) lts
1.1535 - val eqp0 = map (fn (t,i) => (t,Axiom_eq i)) (eq0 ~~ (0 upto (length eq0 - 1)))
1.1536 - val lep0 = map (fn (t,i) => (t,Axiom_le i)) (le0 ~~ (0 upto (length le0 - 1)))
1.1537 - val ltp0 = map (fn (t,i) => (t,Axiom_lt i)) (lt0 ~~ (0 upto (length lt0 - 1)))
1.1538 - val (keq,eq) = List.partition (fn (p,_) => multidegree p = 0) eqp0
1.1539 - val (klep,lep) = List.partition (fn (p,_) => multidegree p = 0) lep0
1.1540 - val (kltp,ltp) = List.partition (fn (p,_) => multidegree p = 0) ltp0
1.1541 - fun trivial_axiom (p,ax) =
1.1542 - case ax of
1.1543 - Axiom_eq n => if eval Ctermfunc.undefined p <>/ Rat.zero then nth eqs n
1.1544 - else error "trivial_axiom: Not a trivial axiom"
1.1545 - | Axiom_le n => if eval Ctermfunc.undefined p </ Rat.zero then nth les n
1.1546 - else error "trivial_axiom: Not a trivial axiom"
1.1547 - | Axiom_lt n => if eval Ctermfunc.undefined p <=/ Rat.zero then nth lts n
1.1548 - else error "trivial_axiom: Not a trivial axiom"
1.1549 - | _ => error "trivial_axiom: Not a trivial axiom"
1.1550 - in
1.1551 - ((let val th = tryfind trivial_axiom (keq @ klep @ kltp)
1.1552 - in fconv_rule (arg_conv (arg1_conv real_poly_conv) then_conv field_comp_conv) th end)
1.1553 - handle ERROR _ => (
1.1554 - let
1.1555 - val pol = fold_rev poly_mul (map fst ltp) (poly_const Rat.one)
1.1556 - val leq = lep @ ltp
1.1557 - fun tryall d =
1.1558 - let val e = multidegree pol
1.1559 - val k = if e = 0 then 0 else d div e
1.1560 - val eq' = map fst eq
1.1561 - in tryfind (fn i => (d,i,real_positivnullstellensatz_general prover false d eq' leq
1.1562 - (poly_neg(poly_pow pol i))))
1.1563 - (0 upto k)
1.1564 - end
1.1565 - val (d,i,(cert_ideal,cert_cone)) = deepen tryall 0
1.1566 - val proofs_ideal =
1.1567 - map2 (fn q => fn (p,ax) => Eqmul(cterm_of_poly q,ax)) cert_ideal eq
1.1568 - val proofs_cone = map cterm_of_sos cert_cone
1.1569 - val proof_ne = if null ltp then Rational_lt Rat.one else
1.1570 - let val p = fold1 (fn s => fn t => Product(s,t)) (map snd ltp)
1.1571 - in funpow i (fn q => Product(p,q)) (Rational_lt Rat.one)
1.1572 - end
1.1573 - val proof = fold1 (fn s => fn t => Sum(s,t))
1.1574 - (proof_ne :: proofs_ideal @ proofs_cone)
1.1575 - in writeln "Translating proof certificate to HOL";
1.1576 - translator (eqs,les,lts) proof
1.1577 - end))
1.1578 - end
1.1579 - in mainf end
1.1580 -end
1.1581 -
1.1582 -fun C f x y = f y x;
1.1583 - (* FIXME : This is very bad!!!*)
1.1584 -fun subst_conv eqs t =
1.1585 - let
1.1586 - val t' = fold (Thm.cabs o Thm.lhs_of) eqs t
1.1587 - in Conv.fconv_rule (Thm.beta_conversion true) (fold (C combination) eqs (reflexive t'))
1.1588 - end
1.1589 -
1.1590 -(* A wrapper that tries to substitute away variables first. *)
1.1591 -
1.1592 -local
1.1593 - open Thm Conv RealArith
1.1594 - fun simple_cterm_ord t u = TermOrd.fast_term_ord (term_of t, term_of u) = LESS
1.1595 - val concl = dest_arg o cprop_of
1.1596 - val shuffle1 =
1.1597 - fconv_rule (rewr_conv @{lemma "(a + x == y) == (x == y - (a::real))" by (atomize (full)) (simp add: ring_simps) })
1.1598 - val shuffle2 =
1.1599 - fconv_rule (rewr_conv @{lemma "(x + a == y) == (x == y - (a::real))" by (atomize (full)) (simp add: ring_simps)})
1.1600 - fun substitutable_monomial fvs tm = case term_of tm of
1.1601 - Free(_,@{typ real}) => if not (member (op aconvc) fvs tm) then (Rat.one,tm)
1.1602 - else error "substitutable_monomial"
1.1603 - | @{term "op * :: real => _"}$c$(t as Free _ ) =>
1.1604 - if is_ratconst (dest_arg1 tm) andalso not (member (op aconvc) fvs (dest_arg tm))
1.1605 - then (dest_ratconst (dest_arg1 tm),dest_arg tm) else error "substitutable_monomial"
1.1606 - | @{term "op + :: real => _"}$s$t =>
1.1607 - (substitutable_monomial (add_cterm_frees (dest_arg tm) fvs) (dest_arg1 tm)
1.1608 - handle ERROR _ => substitutable_monomial (add_cterm_frees (dest_arg1 tm) fvs) (dest_arg tm))
1.1609 - | _ => error "substitutable_monomial"
1.1610 -
1.1611 - fun isolate_variable v th =
1.1612 - let val w = dest_arg1 (cprop_of th)
1.1613 - in if v aconvc w then th
1.1614 - else case term_of w of
1.1615 - @{term "op + :: real => _"}$s$t =>
1.1616 - if dest_arg1 w aconvc v then shuffle2 th
1.1617 - else isolate_variable v (shuffle1 th)
1.1618 - | _ => error "isolate variable : This should not happen?"
1.1619 - end
1.1620 -in
1.1621 -
1.1622 -fun real_nonlinear_subst_prover prover ctxt =
1.1623 - let
1.1624 - val {add,mul,neg,pow,sub,main} = Normalizer.semiring_normalizers_ord_wrapper ctxt
1.1625 - (valOf (NormalizerData.match ctxt @{cterm "(0::real) + 1"}))
1.1626 - simple_cterm_ord
1.1627 -
1.1628 - val (real_poly_add_conv,real_poly_mul_conv,real_poly_neg_conv,
1.1629 - real_poly_pow_conv,real_poly_sub_conv,real_poly_conv) = (add,mul,neg,pow,sub,main)
1.1630 -
1.1631 - fun make_substitution th =
1.1632 - let
1.1633 - val (c,v) = substitutable_monomial [] (dest_arg1(concl th))
1.1634 - val th1 = Drule.arg_cong_rule (capply @{cterm "op * :: real => _"} (cterm_of_rat (Rat.inv c))) (mk_meta_eq th)
1.1635 - val th2 = fconv_rule (binop_conv real_poly_mul_conv) th1
1.1636 - in fconv_rule (arg_conv real_poly_conv) (isolate_variable v th2)
1.1637 - end
1.1638 - fun oprconv cv ct =
1.1639 - let val g = Thm.dest_fun2 ct
1.1640 - in if g aconvc @{cterm "op <= :: real => _"}
1.1641 - orelse g aconvc @{cterm "op < :: real => _"}
1.1642 - then arg_conv cv ct else arg1_conv cv ct
1.1643 - end
1.1644 - fun mainf translator =
1.1645 - let
1.1646 - fun substfirst(eqs,les,lts) =
1.1647 - ((let
1.1648 - val eth = tryfind make_substitution eqs
1.1649 - val modify = fconv_rule (arg_conv (oprconv(subst_conv [eth] then_conv real_poly_conv)))
1.1650 - in substfirst
1.1651 - (filter_out (fn t => (Thm.dest_arg1 o Thm.dest_arg o cprop_of) t
1.1652 - aconvc @{cterm "0::real"}) (map modify eqs),
1.1653 - map modify les,map modify lts)
1.1654 - end)
1.1655 - handle ERROR _ => real_nonlinear_prover prover ctxt translator (rev eqs, rev les, rev lts))
1.1656 - in substfirst
1.1657 - end
1.1658 -
1.1659 -
1.1660 - in mainf
1.1661 - end
1.1662 -
1.1663 -(* Overall function. *)
1.1664 -
1.1665 -fun real_sos prover ctxt t = gen_prover_real_arith ctxt (real_nonlinear_subst_prover prover ctxt) t;
1.1666 -end;
1.1667 -
1.1668 -(* A tactic *)
1.1669 -fun strip_all ct =
1.1670 - case term_of ct of
1.1671 - Const("all",_) $ Abs (xn,xT,p) =>
1.1672 - let val (a,(v,t')) = (apsnd (Thm.dest_abs (SOME xn)) o Thm.dest_comb) ct
1.1673 - in apfst (cons v) (strip_all t')
1.1674 - end
1.1675 -| _ => ([],ct)
1.1676 -
1.1677 -fun core_sos_conv prover ctxt t = Drule.arg_cong_rule @{cterm Trueprop} (real_sos prover ctxt (Thm.dest_arg t) RS @{thm Eq_TrueI})
1.1678 -
1.1679 -val known_sos_constants =
1.1680 - [@{term "op ==>"}, @{term "Trueprop"},
1.1681 - @{term "op -->"}, @{term "op &"}, @{term "op |"},
1.1682 - @{term "Not"}, @{term "op = :: bool => _"},
1.1683 - @{term "All :: (real => _) => _"}, @{term "Ex :: (real => _) => _"},
1.1684 - @{term "op = :: real => _"}, @{term "op < :: real => _"},
1.1685 - @{term "op <= :: real => _"},
1.1686 - @{term "op + :: real => _"}, @{term "op - :: real => _"},
1.1687 - @{term "op * :: real => _"}, @{term "uminus :: real => _"},
1.1688 - @{term "op / :: real => _"}, @{term "inverse :: real => _"},
1.1689 - @{term "op ^ :: real => _"}, @{term "abs :: real => _"},
1.1690 - @{term "min :: real => _"}, @{term "max :: real => _"},
1.1691 - @{term "0::real"}, @{term "1::real"}, @{term "number_of :: int => real"},
1.1692 - @{term "number_of :: int => nat"},
1.1693 - @{term "Int.Bit0"}, @{term "Int.Bit1"},
1.1694 - @{term "Int.Pls"}, @{term "Int.Min"}];
1.1695 -
1.1696 -fun check_sos kcts ct =
1.1697 - let
1.1698 - val t = term_of ct
1.1699 - val _ = if not (null (Term.add_tfrees t [])
1.1700 - andalso null (Term.add_tvars t []))
1.1701 - then error "SOS: not sos. Additional type varables" else ()
1.1702 - val fs = Term.add_frees t []
1.1703 - val _ = if exists (fn ((_,T)) => not (T = @{typ "real"})) fs
1.1704 - then error "SOS: not sos. Variables with type not real" else ()
1.1705 - val vs = Term.add_vars t []
1.1706 - val _ = if exists (fn ((_,T)) => not (T = @{typ "real"})) fs
1.1707 - then error "SOS: not sos. Variables with type not real" else ()
1.1708 - val ukcs = subtract (fn (t,p) => Const p aconv t) kcts (Term.add_consts t [])
1.1709 - val _ = if null ukcs then ()
1.1710 - else error ("SOSO: Unknown constants in Subgoal:" ^ commas (map fst ukcs))
1.1711 -in () end
1.1712 -
1.1713 -fun core_sos_tac prover ctxt = CSUBGOAL (fn (ct, i) =>
1.1714 - let val _ = check_sos known_sos_constants ct
1.1715 - val (avs, p) = strip_all ct
1.1716 - val th = standard (fold_rev forall_intr avs (real_sos prover ctxt (Thm.dest_arg p)))
1.1717 - in rtac th i end);
1.1718 -
1.1719 -fun default_SOME f NONE v = SOME v
1.1720 - | default_SOME f (SOME v) _ = SOME v;
1.1721 -
1.1722 -fun lift_SOME f NONE a = f a
1.1723 - | lift_SOME f (SOME a) _ = SOME a;
1.1724 -
1.1725 -
1.1726 -local
1.1727 - val is_numeral = can (HOLogic.dest_number o term_of)
1.1728 -in
1.1729 -fun get_denom b ct = case term_of ct of
1.1730 - @{term "op / :: real => _"} $ _ $ _ =>
1.1731 - if is_numeral (Thm.dest_arg ct) then get_denom b (Thm.dest_arg1 ct)
1.1732 - else default_SOME (get_denom b) (get_denom b (Thm.dest_arg ct)) (Thm.dest_arg ct, b)
1.1733 - | @{term "op < :: real => _"} $ _ $ _ => lift_SOME (get_denom true) (get_denom true (Thm.dest_arg ct)) (Thm.dest_arg1 ct)
1.1734 - | @{term "op <= :: real => _"} $ _ $ _ => lift_SOME (get_denom true) (get_denom true (Thm.dest_arg ct)) (Thm.dest_arg1 ct)
1.1735 - | _ $ _ => lift_SOME (get_denom b) (get_denom b (Thm.dest_fun ct)) (Thm.dest_arg ct)
1.1736 - | _ => NONE
1.1737 -end;
1.1738 -
1.1739 -fun elim_one_denom_tac ctxt =
1.1740 -CSUBGOAL (fn (P,i) =>
1.1741 - case get_denom false P of
1.1742 - NONE => no_tac
1.1743 - | SOME (d,ord) =>
1.1744 - let
1.1745 - val ss = simpset_of ctxt addsimps @{thms field_simps}
1.1746 - addsimps [@{thm nonzero_power_divide}, @{thm power_divide}]
1.1747 - val th = instantiate' [] [SOME d, SOME (Thm.dest_arg P)]
1.1748 - (if ord then @{lemma "(d=0 --> P) & (d>0 --> P) & (d<(0::real) --> P) ==> P" by auto}
1.1749 - else @{lemma "(d=0 --> P) & (d ~= (0::real) --> P) ==> P" by blast})
1.1750 - in (rtac th i THEN Simplifier.asm_full_simp_tac ss i) end);
1.1751 -
1.1752 -fun elim_denom_tac ctxt i = REPEAT (elim_one_denom_tac ctxt i);
1.1753 -
1.1754 -fun sos_tac prover ctxt = ObjectLogic.full_atomize_tac THEN' elim_denom_tac ctxt THEN' core_sos_tac prover ctxt
1.1755 -
1.1756 -
1.1757 -end;