author chaieb Tue May 12 17:32:49 2009 +0100 (2009-05-12) changeset 31119 2532bb2d65c7 parent 31118 541d43bee678 child 31120 fc654c95c29e
A decision method for universal multivariate real arithmetic with add
ition, multiplication and ordering using semidefinite programming
```     1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/src/HOL/Library/Sum_Of_Squares.thy	Tue May 12 17:32:49 2009 +0100
1.3 @@ -0,0 +1,114 @@
1.4 +(* Title:      Library/Sum_Of_Squares
1.5 +   Author:     Amine Chaieb, University of Cambridge
1.6 +*)
1.7 +
1.8 +header {* A decision method for universal multivariate real arithmetic with addition,
1.9 +          multiplication and ordering using semidefinite programming*}
1.10 +theory Sum_Of_Squares
1.11 +  imports Complex_Main (* "~~/src/HOL/Decision_Procs/Dense_Linear_Order" *)
1.12 +  uses "positivstellensatz.ML" "sum_of_squares.ML"
1.13 +  begin
1.14 +
1.15 +method_setup sos = {*
1.16 +let
1.17 + fun strip_all ct =
1.18 +  case term_of ct of
1.19 +   Const("all",_) \$ Abs (xn,xT,p) =>
1.20 +    let val (a,(v,t')) = (apsnd (Thm.dest_abs (SOME xn)) o Thm.dest_comb) ct
1.21 +    in apfst (cons v) (strip_all t')
1.22 +    end
1.23 + | _ => ([],ct)
1.24 +
1.25 + fun core_sos_conv ctxt t = Drule.arg_cong_rule @{cterm Trueprop} (Sos.real_sos ctxt (Thm.dest_arg t) RS @{thm Eq_TrueI})
1.26 + fun core_sos_tac ctxt = CSUBGOAL (fn (ct, i) =>
1.27 +   let val (avs, p) = strip_all ct
1.28 +       val th = standard (fold_rev forall_intr avs (Sos.real_sos ctxt (Thm.dest_arg p)))
1.29 +   in rtac th i end) (* CONVERSION o core_sos_conv *)
1.30 +in Scan.succeed (SIMPLE_METHOD' o core_sos_tac)
1.31 +end
1.32 +
1.33 +*} "Prove universal problems over the reals using sums of squares"
1.34 +
1.35 +text{* Tests -- commented since they work only when csdp is installed *}
1.36 +
1.37 +(*
1.38 +lemma "(3::real) * x + 7 * a < 4 & 3 < 2 * x --> a < 0" by sos
1.39 +
1.40 +lemma "a1 >= 0 & a2 >= 0 \<and> (a1 * a1 + a2 * a2 = b1 * b1 + b2 * b2 + 2) \<and> (a1 * b1 + a2 * b2 = 0) --> a1 * a2 - b1 * b2 >= (0::real)" by sos
1.41 +
1.42 +lemma "(3::real) * x + 7 * a < 4 & 3 < 2 * x --> a < 0" by sos
1.43 +
1.44 +lemma "(0::real) <= x & x <= 1 & 0 <= y & y <= 1  --> x^2 + y^2 < 1 |(x - 1)^2 + y^2 < 1 | x^2 + (y - 1)^2 < 1 | (x - 1)^2 + (y - 1)^2 < 1" by sos
1.45 +
1.46 +lemma "(0::real) <= x & 0 <= y & 0 <= z & x + y + z <= 3 --> x * y + x * z + y * z >= 3 * x * y * z" by sos
1.47 +
1.48 +lemma "((x::real)^2 + y^2 + z^2 = 1) --> (x + y + z)^2 <= 3" by sos
1.49 +
1.50 +lemma "(w^2 + x^2 + y^2 + z^2 = 1) --> (w + x + y + z)^2 <= (4::real)" by sos
1.51 +
1.52 +lemma "(x::real) >= 1 & y >= 1 --> x * y >= x + y - 1" by sos
1.53 +
1.54 +lemma "(x::real) > 1 & y > 1 --> x * y > x + y - 1" by sos;
1.55 +
1.56 +lemma "abs(x) <= 1 --> abs(64 * x^7 - 112 * x^5 + 56 * x^3 - 7 * x) <= (1::real)" by sos
1.57 +*)
1.58 +(* ------------------------------------------------------------------------- *)
1.59 +(* One component of denominator in dodecahedral example.                     *)
1.60 +(* ------------------------------------------------------------------------- *)
1.61 +(*
1.62 +lemma "2 <= x & x <= 125841 / 50000 & 2 <= y & y <= 125841 / 50000 & 2 <= z & z <= 125841 / 50000 --> 2 * (x * z + x * y + y * z) - (x * x + y * y + z * z) >= (0::real)" by sos;
1.63 +*)
1.64 +(* ------------------------------------------------------------------------- *)
1.65 +(* Over a larger but simpler interval.                                       *)
1.66 +(* ------------------------------------------------------------------------- *)
1.67 +(*
1.68 +lemma "(2::real) <= x & x <= 4 & 2 <= y & y <= 4 & 2 <= z & z <= 4 --> 0 <= 2 * (x * z + x * y + y * z) - (x * x + y * y + z * z)" by sos
1.69 +*)
1.70 +(* ------------------------------------------------------------------------- *)
1.71 +(* We can do 12. I think 12 is a sharp bound; see PP's certificate.          *)
1.72 +(* ------------------------------------------------------------------------- *)
1.73 +(*
1.74 +lemma "2 <= (x::real) & x <= 4 & 2 <= y & y <= 4 & 2 <= z & z <= 4 --> 12 <= 2 * (x * z + x * y + y * z) - (x * x + y * y + z * z)" by sos
1.75 +*)
1.76 +
1.77 +(* ------------------------------------------------------------------------- *)
1.78 +(* Inequality from sci.math (see "Leon-Sotelo, por favor").                  *)
1.79 +(* ------------------------------------------------------------------------- *)
1.80 +(*
1.81 +lemma "0 <= (x::real) & 0 <= y & (x * y = 1) --> x + y <= x^2 + y^2" by sos
1.82 +
1.83 +lemma "0 <= (x::real) & 0 <= y & (x * y = 1) --> x * y * (x + y) <= x^2 + y^2" by sos
1.84 +
1.85 +lemma "0 <= (x::real) & 0 <= y --> x * y * (x + y)^2 <= (x^2 + y^2)^2" by sos
1.86 +
1.87 +lemma "(0::real) <= a & 0 <= b & 0 <= c & c * (2 * a + b)^3/ 27 <= x \<longrightarrow> c * a^2 * b <= x" by sos
1.88 +
1.89 +lemma "(0::real) < x --> 0 < 1 + x + x^2" by sos
1.90 +
1.91 +lemma "(0::real) <= x --> 0 < 1 + x + x^2" by sos
1.92 +
1.93 +lemma "(0::real) < 1 + x^2" by sos
1.94 +
1.95 +lemma "(0::real) <= 1 + 2 * x + x^2" by sos
1.96 +
1.97 +lemma "(0::real) < 1 + abs x" by sos
1.98 +
1.99 +lemma "(0::real) < 1 + (1 + x)^2 * (abs x)" by sos
1.100 +
1.101 +
1.102 +
1.103 +lemma "abs ((1::real) + x^2) = (1::real) + x^2" by sos
1.104 +lemma "(3::real) * x + 7 * a < 4 \<and> 3 < 2 * x \<longrightarrow> a < 0" by sos
1.105 +
1.106 +lemma "(0::real) < x --> 1 < y --> y * x <= z --> x < z" by sos
1.107 +lemma "(1::real) < x --> x^2 < y --> 1 < y" by sos
1.108 +lemma "(b::real)^2 < 4 * a * c --> ~(a * x^2 + b * x + c = 0)" by sos
1.109 +lemma "(b::real)^2 < 4 * a * c --> ~(a * x^2 + b * x + c = 0)" by sos
1.110 +lemma "((a::real) * x^2 + b * x + c = 0) --> b^2 >= 4 * a * c" by sos
1.111 +lemma "(0::real) <= b & 0 <= c & 0 <= x & 0 <= y & (x^2 = c) & (y^2 = a^2 * c + b) --> a * c <= y * x" by sos
1.112 +lemma "abs(x - z) <= e & abs(y - z) <= e & 0 <= u & 0 <= v & (u + v = 1) --> abs((u * x + v * y) - z) <= (e::real)" by sos
1.113 +*)
1.114 +(*
1.115 +lemma "((x::real) - y - 2 * x^4 = 0) & 0 <= x & x <= 2 & 0 <= y & y <= 3 --> y^2 - 7 * y - 12 * x + 17 >= 0" by sos *) (* Too hard?*)
1.116 +
1.117 +end
```
```     2.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
2.2 +++ b/src/HOL/Library/sum_of_squares.ML	Tue May 12 17:32:49 2009 +0100
2.3 @@ -0,0 +1,1612 @@
2.4 +structure Sos =
2.5 +struct
2.6 +
2.7 +val rat_0 = Rat.zero;
2.8 +val rat_1 = Rat.one;
2.9 +val rat_2 = Rat.two;
2.10 +val rat_10 = Rat.rat_of_int 10;
2.11 +val rat_1_2 = rat_1 // rat_2;
2.12 +val max = curry IntInf.max;
2.13 +val min = curry IntInf.min;
2.14 +
2.15 +val denominator_rat = Rat.quotient_of_rat #> snd #> Rat.rat_of_int;
2.16 +val numerator_rat = Rat.quotient_of_rat #> fst #> Rat.rat_of_int;
2.17 +fun int_of_rat a =
2.18 +    case Rat.quotient_of_rat a of (i,1) => i | _ => error "int_of_rat: not an int";
2.19 +fun lcm_rat x y = Rat.rat_of_int (Integer.lcm (int_of_rat x) (int_of_rat y));
2.20 +
2.21 +fun rat_pow r i =
2.22 + let fun pow r i =
2.23 +   if i = 0 then rat_1 else
2.24 +   let val d = pow r (i div 2)
2.25 +   in d */ d */ (if i mod 2 = 0 then rat_1 else r)
2.26 +   end
2.27 + in if i < 0 then pow (Rat.inv r) (~ i) else pow r i end;
2.28 +
2.29 +fun round_rat r =
2.30 + let val (a,b) = Rat.quotient_of_rat (Rat.abs r)
2.31 +     val d = a div b
2.32 +     val s = if r </ rat_0 then (Rat.neg o Rat.rat_of_int) else Rat.rat_of_int
2.33 +     val x2 = 2 * (a - (b * d))
2.34 + in s (if x2 >= b then d + 1 else d) end
2.35 +
2.36 +val abs_rat = Rat.abs;
2.37 +val pow2 = rat_pow rat_2;
2.38 +val pow10 = rat_pow rat_10;
2.39 +
2.40 +val debugging = ref false;
2.41 +
2.42 +exception Sanity;
2.43 +
2.44 +exception Unsolvable;
2.45 +
2.46 +(* Turn a rational into a decimal string with d sig digits.                  *)
2.47 +
2.48 +local
2.49 +fun normalize y =
2.50 +  if abs_rat y </ (rat_1 // rat_10) then normalize (rat_10 */ y) - 1
2.51 +  else if abs_rat y >=/ rat_1 then normalize (y // rat_10) + 1
2.52 +  else 0
2.53 + in
2.54 +fun decimalize d x =
2.55 +  if x =/ rat_0 then "0.0" else
2.56 +  let
2.57 +   val y = Rat.abs x
2.58 +   val e = normalize y
2.59 +   val z = pow10(~ e) */ y +/ rat_1
2.60 +   val k = int_of_rat (round_rat(pow10 d */ z))
2.61 +  in (if x </ rat_0 then "-0." else "0.") ^
2.62 +     implode(tl(explode(string_of_int k))) ^
2.63 +     (if e = 0 then "" else "e"^string_of_int e)
2.64 +  end
2.65 +end;
2.66 +
2.67 +(* Iterations over numbers, and lists indexed by numbers.                    *)
2.68 +
2.69 +fun itern k l f a =
2.70 +  case l of
2.71 +    [] => a
2.72 +  | h::t => itern (k + 1) t f (f h k a);
2.73 +
2.74 +fun iter (m,n) f a =
2.75 +  if n < m then a
2.76 +  else iter (m+1,n) f (f m a);
2.77 +
2.78 +(* The main types.                                                           *)
2.79 +
2.80 +fun strict_ord ord (x,y) = case ord (x,y) of LESS => LESS | _ => GREATER
2.81 +
2.82 +structure Intpairfunc = FuncFun(type key = int*int val ord = prod_ord int_ord int_ord);
2.83 +
2.84 +type vector = int* Rat.rat Intfunc.T;
2.85 +
2.86 +type matrix = (int*int)*(Rat.rat Intpairfunc.T);
2.87 +
2.88 +type monomial = int Ctermfunc.T;
2.89 +
2.90 +val cterm_ord = (fn (s,t) => TermOrd.fast_term_ord(term_of s, term_of t))
2.91 + fun monomial_ord (m1,m2) = list_ord (prod_ord cterm_ord int_ord) (Ctermfunc.graph m1, Ctermfunc.graph m2)
2.92 +structure Monomialfunc = FuncFun(type key = monomial val ord = monomial_ord)
2.93 +
2.94 +type poly = Rat.rat Monomialfunc.T;
2.95 +
2.96 + fun iszero (k,r) = r =/ rat_0;
2.97 +
2.98 +fun fold_rev2 f l1 l2 b =
2.99 +  case (l1,l2) of
2.100 +    ([],[]) => b
2.101 +  | (h1::t1,h2::t2) => f h1 h2 (fold_rev2 f t1 t2 b)
2.102 +  | _ => error "fold_rev2";
2.103 +
2.104 +(* Vectors. Conventionally indexed 1..n.                                     *)
2.105 +
2.106 +fun vector_0 n = (n,Intfunc.undefined):vector;
2.107 +
2.108 +fun dim (v:vector) = fst v;
2.109 +
2.110 +fun vector_const c n =
2.111 +  if c =/ rat_0 then vector_0 n
2.112 +  else (n,fold_rev (fn k => Intfunc.update (k,c)) (1 upto n) Intfunc.undefined) :vector;
2.113 +
2.114 +val vector_1 = vector_const rat_1;
2.115 +
2.116 +fun vector_cmul c (v:vector) =
2.117 + let val n = dim v
2.118 + in if c =/ rat_0 then vector_0 n
2.119 +    else (n,Intfunc.mapf (fn x => c */ x) (snd v))
2.120 + end;
2.121 +
2.122 +fun vector_neg (v:vector) = (fst v,Intfunc.mapf Rat.neg (snd v)) :vector;
2.123 +
2.124 +fun vector_add (v1:vector) (v2:vector) =
2.125 + let val m = dim v1
2.126 +     val n = dim v2
2.127 + in if m <> n then error "vector_add: incompatible dimensions"
2.128 +    else (n,Intfunc.combine (curry op +/) (fn x => x =/ rat_0) (snd v1) (snd v2)) :vector
2.129 + end;
2.130 +
2.131 +fun vector_sub v1 v2 = vector_add v1 (vector_neg v2);
2.132 +
2.133 +fun vector_dot (v1:vector) (v2:vector) =
2.134 + let val m = dim v1
2.135 +     val n = dim v2
2.136 + in if m <> n then error "vector_dot: incompatible dimensions"
2.137 +    else Intfunc.fold (fn (i,x) => fn a => x +/ a)
2.138 +        (Intfunc.combine (curry op */) (fn x => x =/ rat_0) (snd v1) (snd v2)) rat_0
2.139 + end;
2.140 +
2.141 +fun vector_of_list l =
2.142 + let val n = length l
2.143 + in (n,fold_rev2 (curry Intfunc.update) (1 upto n) l Intfunc.undefined) :vector
2.144 + end;
2.145 +
2.146 +(* Matrices; again rows and columns indexed from 1.                          *)
2.147 +
2.148 +fun matrix_0 (m,n) = ((m,n),Intpairfunc.undefined):matrix;
2.149 +
2.150 +fun dimensions (m:matrix) = fst m;
2.151 +
2.152 +fun matrix_const c (mn as (m,n)) =
2.153 +  if m <> n then error "matrix_const: needs to be square"
2.154 +  else if c =/ rat_0 then matrix_0 mn
2.155 +  else (mn,fold_rev (fn k => Intpairfunc.update ((k,k), c)) (1 upto n) Intpairfunc.undefined) :matrix;;
2.156 +
2.157 +val matrix_1 = matrix_const rat_1;
2.158 +
2.159 +fun matrix_cmul c (m:matrix) =
2.160 + let val (i,j) = dimensions m
2.161 + in if c =/ rat_0 then matrix_0 (i,j)
2.162 +    else ((i,j),Intpairfunc.mapf (fn x => c */ x) (snd m))
2.163 + end;
2.164 +
2.165 +fun matrix_neg (m:matrix) =
2.166 +  (dimensions m, Intpairfunc.mapf Rat.neg (snd m)) :matrix;
2.167 +
2.168 +fun matrix_add (m1:matrix) (m2:matrix) =
2.169 + let val d1 = dimensions m1
2.170 +     val d2 = dimensions m2
2.171 + in if d1 <> d2
2.172 +     then error "matrix_add: incompatible dimensions"
2.173 +    else (d1,Intpairfunc.combine (curry op +/) (fn x => x =/ rat_0) (snd m1) (snd m2)) :matrix
2.174 + end;;
2.175 +
2.176 +fun matrix_sub m1 m2 = matrix_add m1 (matrix_neg m2);
2.177 +
2.178 +fun row k (m:matrix) =
2.179 + let val (i,j) = dimensions m
2.180 + in (j,
2.181 +   Intpairfunc.fold (fn ((i,j), c) => fn a => if i = k then Intfunc.update (j,c) a else a) (snd m) Intfunc.undefined ) : vector
2.182 + end;
2.183 +
2.184 +fun column k (m:matrix) =
2.185 +  let val (i,j) = dimensions m
2.186 +  in (i,
2.187 +   Intpairfunc.fold (fn ((i,j), c) => fn a => if j = k then Intfunc.update (i,c) a else a) (snd m)  Intfunc.undefined)
2.188 +   : vector
2.189 + end;
2.190 +
2.191 +fun transp (m:matrix) =
2.192 +  let val (i,j) = dimensions m
2.193 +  in
2.194 +  ((j,i),Intpairfunc.fold (fn ((i,j), c) => fn a => Intpairfunc.update ((j,i), c) a) (snd m) Intpairfunc.undefined) :matrix
2.195 + end;
2.196 +
2.197 +fun diagonal (v:vector) =
2.198 + let val n = dim v
2.199 + in ((n,n),Intfunc.fold (fn (i, c) => fn a => Intpairfunc.update ((i,i), c) a) (snd v) Intpairfunc.undefined) : matrix
2.200 + end;
2.201 +
2.202 +fun matrix_of_list l =
2.203 + let val m = length l
2.204 + in if m = 0 then matrix_0 (0,0) else
2.205 +   let val n = length (hd l)
2.206 +   in ((m,n),itern 1 l (fn v => fn i => itern 1 v (fn c => fn j => Intpairfunc.update ((i,j), c))) Intpairfunc.undefined)
2.207 +   end
2.208 + end;
2.209 +
2.210 +(* Monomials.                                                                *)
2.211 +
2.212 +fun monomial_eval assig (m:monomial) =
2.213 +  Ctermfunc.fold (fn (x, k) => fn a => a */ rat_pow (Ctermfunc.apply assig x) k)
2.214 +        m rat_1;
2.215 +val monomial_1 = (Ctermfunc.undefined:monomial);
2.216 +
2.217 +fun monomial_var x = Ctermfunc.onefunc (x, 1) :monomial;
2.218 +
2.219 +val (monomial_mul:monomial->monomial->monomial) =
2.220 +  Ctermfunc.combine (curry op +) (K false);
2.221 +
2.222 +fun monomial_pow (m:monomial) k =
2.223 +  if k = 0 then monomial_1
2.224 +  else Ctermfunc.mapf (fn x => k * x) m;
2.225 +
2.226 +fun monomial_divides (m1:monomial) (m2:monomial) =
2.227 +  Ctermfunc.fold (fn (x, k) => fn a => Ctermfunc.tryapplyd m2 x 0 >= k andalso a) m1 true;;
2.228 +
2.229 +fun monomial_div (m1:monomial) (m2:monomial) =
2.230 + let val m = Ctermfunc.combine (curry op +)
2.231 +   (fn x => x = 0) m1 (Ctermfunc.mapf (fn x => ~ x) m2)
2.232 + in if Ctermfunc.fold (fn (x, k) => fn a => k >= 0 andalso a) m true then m
2.233 +  else error "monomial_div: non-divisible"
2.234 + end;
2.235 +
2.236 +fun monomial_degree x (m:monomial) =
2.237 +  Ctermfunc.tryapplyd m x 0;;
2.238 +
2.239 +fun monomial_lcm (m1:monomial) (m2:monomial) =
2.240 +  fold_rev (fn x => Ctermfunc.update (x, max (monomial_degree x m1) (monomial_degree x m2)))
2.241 +          (gen_union (is_equal o  cterm_ord) (Ctermfunc.dom m1, Ctermfunc.dom m2)) (Ctermfunc.undefined :monomial);
2.242 +
2.243 +fun monomial_multidegree (m:monomial) =
2.244 + Ctermfunc.fold (fn (x, k) => fn a => k + a) m 0;;
2.245 +
2.246 +fun monomial_variables m = Ctermfunc.dom m;;
2.247 +
2.248 +(* Polynomials.                                                              *)
2.249 +
2.250 +fun eval assig (p:poly) =
2.251 +  Monomialfunc.fold (fn (m, c) => fn a => a +/ c */ monomial_eval assig m) p rat_0;
2.252 +
2.253 +val poly_0 = (Monomialfunc.undefined:poly);
2.254 +
2.255 +fun poly_isconst (p:poly) =
2.256 +  Monomialfunc.fold (fn (m, c) => fn a => Ctermfunc.is_undefined m andalso a) p true;
2.257 +
2.258 +fun poly_var x = Monomialfunc.onefunc (monomial_var x,rat_1) :poly;
2.259 +
2.260 +fun poly_const c =
2.261 +  if c =/ rat_0 then poly_0 else Monomialfunc.onefunc(monomial_1, c);
2.262 +
2.263 +fun poly_cmul c (p:poly) =
2.264 +  if c =/ rat_0 then poly_0
2.265 +  else Monomialfunc.mapf (fn x => c */ x) p;
2.266 +
2.267 +fun poly_neg (p:poly) = (Monomialfunc.mapf Rat.neg p :poly);;
2.268 +
2.269 +fun poly_add (p1:poly) (p2:poly) =
2.270 +  (Monomialfunc.combine (curry op +/) (fn x => x =/ rat_0) p1 p2 :poly);
2.271 +
2.272 +fun poly_sub p1 p2 = poly_add p1 (poly_neg p2);
2.273 +
2.274 +fun poly_cmmul (c,m) (p:poly) =
2.275 + if c =/ rat_0 then poly_0
2.276 + else if Ctermfunc.is_undefined m
2.277 +      then Monomialfunc.mapf (fn d => c */ d) p
2.278 +      else Monomialfunc.fold (fn (m', d) => fn a => (Monomialfunc.update (monomial_mul m m', c */ d) a)) p poly_0;
2.279 +
2.280 +fun poly_mul (p1:poly) (p2:poly) =
2.281 +  Monomialfunc.fold (fn (m, c) => fn a => poly_add (poly_cmmul (c,m) p2) a) p1 poly_0;
2.282 +
2.283 +fun poly_div (p1:poly) (p2:poly) =
2.284 + if not(poly_isconst p2)
2.285 + then error "poly_div: non-constant" else
2.286 + let val c = eval Ctermfunc.undefined p2
2.287 + in if c =/ rat_0 then error "poly_div: division by zero"
2.288 +    else poly_cmul (Rat.inv c) p1
2.289 + end;
2.290 +
2.291 +fun poly_square p = poly_mul p p;
2.292 +
2.293 +fun poly_pow p k =
2.294 + if k = 0 then poly_const rat_1
2.295 + else if k = 1 then p
2.296 + else let val q = poly_square(poly_pow p (k div 2)) in
2.297 +      if k mod 2 = 1 then poly_mul p q else q end;
2.298 +
2.299 +fun poly_exp p1 p2 =
2.300 +  if not(poly_isconst p2)
2.301 +  then error "poly_exp: not a constant"
2.302 +  else poly_pow p1 (int_of_rat (eval Ctermfunc.undefined p2));
2.303 +
2.304 +fun degree x (p:poly) =
2.305 + Monomialfunc.fold (fn (m,c) => fn a => max (monomial_degree x m) a) p 0;
2.306 +
2.307 +fun multidegree (p:poly) =
2.308 +  Monomialfunc.fold (fn (m, c) => fn a => max (monomial_multidegree m) a) p 0;
2.309 +
2.310 +fun poly_variables (p:poly) =
2.311 +  sort cterm_ord (Monomialfunc.fold_rev (fn (m, c) => curry (gen_union (is_equal o  cterm_ord)) (monomial_variables m)) p []);;
2.312 +
2.313 +(* Order monomials for human presentation.                                   *)
2.314 +
2.315 +fun cterm_ord (t,t') = TermOrd.fast_term_ord (term_of t, term_of t');
2.316 +
2.317 +val humanorder_varpow = prod_ord cterm_ord (rev_order o int_ord);
2.318 +
2.319 +local
2.320 + fun ord (l1,l2) = case (l1,l2) of
2.321 +  (_,[]) => LESS
2.322 + | ([],_) => GREATER
2.323 + | (h1::t1,h2::t2) =>
2.324 +   (case humanorder_varpow (h1, h2) of
2.325 +     LESS => LESS
2.326 +   | EQUAL => ord (t1,t2)
2.327 +   | GREATER => GREATER)
2.328 +in fun humanorder_monomial m1 m2 =
2.329 + ord (sort humanorder_varpow (Ctermfunc.graph m1),
2.330 +  sort humanorder_varpow (Ctermfunc.graph m2))
2.331 +end;
2.332 +
2.333 +fun fold1 f l =  case l of
2.334 +   []     => error "fold1"
2.335 + | [x]    => x
2.336 + | (h::t) => f h (fold1 f t);
2.337 +
2.338 +(* Conversions to strings.                                                   *)
2.339 +
2.340 +fun string_of_vector min_size max_size (v:vector) =
2.341 + let val n_raw = dim v
2.342 + in if n_raw = 0 then "[]" else
2.343 +  let
2.344 +   val n = max min_size (min n_raw max_size)
2.345 +   val xs = map (Rat.string_of_rat o (fn i => Intfunc.tryapplyd (snd v) i rat_0)) (1 upto n)
2.346 +  in "[" ^ fold1 (fn s => fn t => s ^ ", " ^ t) xs ^
2.347 +  (if n_raw > max_size then ", ...]" else "]")
2.348 +  end
2.349 + end;
2.350 +
2.351 +fun string_of_matrix max_size (m:matrix) =
2.352 + let
2.353 +  val (i_raw,j_raw) = dimensions m
2.354 +  val i = min max_size i_raw
2.355 +  val j = min max_size j_raw
2.356 +  val rstr = map (fn k => string_of_vector j j (row k m)) (1 upto i)
2.357 + in "["^ fold1 (fn s => fn t => s^";\n "^t) rstr ^
2.358 +  (if j > max_size then "\n ...]" else "]")
2.359 + end;
2.360 +
2.361 +fun string_of_term t =
2.362 + case t of
2.363 +   a\$b => "("^(string_of_term a)^" "^(string_of_term b)^")"
2.364 + | Abs x =>
2.365 +    let val (xn, b) = Term.dest_abs x
2.366 +    in "(\\"^xn^"."^(string_of_term b)^")"
2.367 +    end
2.368 + | Const(s,_) => s
2.369 + | Free (s,_) => s
2.370 + | Var((s,_),_) => s
2.371 + | _ => error "string_of_term";
2.372 +
2.373 +val string_of_cterm = string_of_term o term_of;
2.374 +
2.375 +fun string_of_varpow x k =
2.376 +  if k = 1 then string_of_cterm x
2.377 +  else string_of_cterm x^"^"^string_of_int k;
2.378 +
2.379 +fun string_of_monomial m =
2.380 + if Ctermfunc.is_undefined m then "1" else
2.381 + let val vps = fold_rev (fn (x,k) => fn a => string_of_varpow x k :: a)
2.382 +  (sort humanorder_varpow (Ctermfunc.graph m)) []
2.383 + in fold1 (fn s => fn t => s^"*"^t) vps
2.384 + end;
2.385 +
2.386 +fun string_of_cmonomial (c,m) =
2.387 + if Ctermfunc.is_undefined m then Rat.string_of_rat c
2.388 + else if c =/ rat_1 then string_of_monomial m
2.389 + else Rat.string_of_rat c ^ "*" ^ string_of_monomial m;;
2.390 +
2.391 +fun string_of_poly (p:poly) =
2.392 + if Monomialfunc.is_undefined p then "<<0>>" else
2.393 + let
2.394 +  val cms = sort (fn ((m1,_),(m2,_)) => humanorder_monomial m1  m2) (Monomialfunc.graph p)
2.395 +  val s = fold (fn (m,c) => fn a =>
2.396 +             if c </ rat_0 then a ^ " - " ^ string_of_cmonomial(Rat.neg c,m)
2.397 +             else a ^ " + " ^ string_of_cmonomial(c,m))
2.398 +          cms ""
2.399 +  val s1 = String.substring (s, 0, 3)
2.400 +  val s2 = String.substring (s, 3, String.size s - 3)
2.401 + in "<<" ^(if s1 = " + " then s2 else "-"^s2)^">>"
2.402 + end;
2.403 +
2.404 +(* Conversion from HOL term.                                                 *)
2.405 +
2.406 +local
2.407 + val neg_tm = @{cterm "uminus :: real => _"}
2.408 + val add_tm = @{cterm "op + :: real => _"}
2.409 + val sub_tm = @{cterm "op - :: real => _"}
2.410 + val mul_tm = @{cterm "op * :: real => _"}
2.411 + val inv_tm = @{cterm "inverse :: real => _"}
2.412 + val div_tm = @{cterm "op / :: real => _"}
2.413 + val pow_tm = @{cterm "op ^ :: real => _"}
2.414 + val zero_tm = @{cterm "0:: real"}
2.415 + val is_numeral = can (HOLogic.dest_number o term_of)
2.416 + fun is_comb t = case t of _\$_ => true | _ => false
2.417 + fun poly_of_term tm =
2.418 +  if tm aconvc zero_tm then poly_0
2.419 +  else if RealArith.is_ratconst tm
2.420 +       then poly_const(RealArith.dest_ratconst tm)
2.421 +  else
2.422 +  (let val (lop,r) = Thm.dest_comb tm
2.423 +   in if lop aconvc neg_tm then poly_neg(poly_of_term r)
2.424 +      else if lop aconvc inv_tm then
2.425 +       let val p = poly_of_term r
2.426 +       in if poly_isconst p
2.427 +          then poly_const(Rat.inv (eval Ctermfunc.undefined p))
2.428 +          else error "poly_of_term: inverse of non-constant polyomial"
2.429 +       end
2.430 +   else (let val (opr,l) = Thm.dest_comb lop
2.431 +         in
2.432 +          if opr aconvc pow_tm andalso is_numeral r
2.433 +          then poly_pow (poly_of_term l) ((snd o HOLogic.dest_number o term_of) r)
2.434 +          else if opr aconvc add_tm
2.435 +           then poly_add (poly_of_term l) (poly_of_term r)
2.436 +          else if opr aconvc sub_tm
2.437 +           then poly_sub (poly_of_term l) (poly_of_term r)
2.438 +          else if opr aconvc mul_tm
2.439 +           then poly_mul (poly_of_term l) (poly_of_term r)
2.440 +          else if opr aconvc div_tm
2.441 +           then let
2.442 +                  val p = poly_of_term l
2.443 +                  val q = poly_of_term r
2.444 +                in if poly_isconst q then poly_cmul (Rat.inv (eval Ctermfunc.undefined q)) p
2.445 +                   else error "poly_of_term: division by non-constant polynomial"
2.446 +                end
2.447 +          else poly_var tm
2.448 +
2.449 +         end
2.450 +         handle CTERM ("dest_comb",_) => poly_var tm)
2.451 +   end
2.452 +   handle CTERM ("dest_comb",_) => poly_var tm)
2.453 +in
2.454 +val poly_of_term = fn tm =>
2.455 + if type_of (term_of tm) = @{typ real} then poly_of_term tm
2.456 + else error "poly_of_term: term does not have real type"
2.457 +end;
2.458 +
2.459 +(* String of vector (just a list of space-separated numbers).                *)
2.460 +
2.461 +fun sdpa_of_vector (v:vector) =
2.462 + let
2.463 +  val n = dim v
2.464 +  val strs = map (decimalize 20 o (fn i => Intfunc.tryapplyd (snd v) i rat_0)) (1 upto n)
2.465 + in fold1 (fn x => fn y => x ^ " " ^ y) strs ^ "\n"
2.466 + end;
2.467 +
2.468 +fun increasing f ord (x,y) = ord (f x, f y);
2.469 +fun triple_int_ord ((a,b,c),(a',b',c')) =
2.470 + prod_ord int_ord (prod_ord int_ord int_ord)
2.471 +    ((a,(b,c)),(a',(b',c')));
2.472 +structure Inttriplefunc = FuncFun(type key = int*int*int val ord = triple_int_ord);
2.473 +
2.474 +(* String for block diagonal matrix numbered k.                              *)
2.475 +
2.476 +fun sdpa_of_blockdiagonal k m =
2.477 + let
2.478 +  val pfx = string_of_int k ^" "
2.479 +  val ents =
2.480 +    Inttriplefunc.fold (fn ((b,i,j), c) => fn a => if i > j then a else ((b,i,j),c)::a) m []
2.481 +  val entss = sort (increasing fst triple_int_ord ) ents
2.482 + in  fold_rev (fn ((b,i,j),c) => fn a =>
2.483 +     pfx ^ string_of_int b ^ " " ^ string_of_int i ^ " " ^ string_of_int j ^
2.484 +     " " ^ decimalize 20 c ^ "\n" ^ a) entss ""
2.485 + end;
2.486 +
2.487 +(* String for a matrix numbered k, in SDPA sparse format.                    *)
2.488 +
2.489 +fun sdpa_of_matrix k (m:matrix) =
2.490 + let
2.491 +  val pfx = string_of_int k ^ " 1 "
2.492 +  val ms = Intpairfunc.fold (fn ((i,j), c) => fn  a => if i > j then a else ((i,j),c)::a) (snd m) []
2.493 +  val mss = sort (increasing fst (prod_ord int_ord int_ord)) ms
2.494 + in fold_rev (fn ((i,j),c) => fn a =>
2.495 +     pfx ^ string_of_int i ^ " " ^ string_of_int j ^
2.496 +     " " ^ decimalize 20 c ^ "\n" ^ a) mss ""
2.497 + end;;
2.498 +
2.499 +(* ------------------------------------------------------------------------- *)
2.500 +(* String in SDPA sparse format for standard SDP problem:                    *)
2.501 +(*                                                                           *)
2.502 +(*    X = v_1 * [M_1] + ... + v_m * [M_m] - [M_0] must be PSD                *)
2.503 +(*    Minimize obj_1 * v_1 + ... obj_m * v_m                                 *)
2.504 +(* ------------------------------------------------------------------------- *)
2.505 +
2.506 +fun sdpa_of_problem comment obj mats =
2.507 + let
2.508 +  val m = length mats - 1
2.509 +  val (n,_) = dimensions (hd mats)
2.510 + in "\"" ^ comment ^ "\"\n" ^
2.511 +  string_of_int m ^ "\n" ^
2.512 +  "1\n" ^
2.513 +  string_of_int n ^ "\n" ^
2.514 +  sdpa_of_vector obj ^
2.515 +  fold_rev2 (fn k => fn m => fn a => sdpa_of_matrix (k - 1) m ^ a) (1 upto length mats) mats ""
2.516 + end;
2.517 +
2.518 +fun index_char str chr pos =
2.519 +  if pos >= String.size str then ~1
2.520 +  else if String.sub(str,pos) = chr then pos
2.521 +  else index_char str chr (pos + 1);
2.522 +fun rat_of_quotient (a,b) = if b = 0 then rat_0 else Rat.rat_of_quotient (a,b);
2.523 +fun rat_of_string s =
2.524 + let val n = index_char s #"/" 0 in
2.525 +  if n = ~1 then s |> IntInf.fromString |> valOf |> Rat.rat_of_int
2.526 +  else
2.527 +   let val SOME numer = IntInf.fromString(String.substring(s,0,n))
2.528 +       val SOME den = IntInf.fromString (String.substring(s,n+1,String.size s - n - 1))
2.529 +   in rat_of_quotient(numer, den)
2.530 +   end
2.531 + end;
2.532 +
2.533 +fun isspace x = x = " " ;
2.534 +fun isnum x = x mem_string ["0","1","2","3","4","5","6","7","8","9"]
2.535 +
2.536 +(* More parser basics.                                                       *)
2.537 +
2.538 +local
2.539 + open Scan
2.540 +in
2.541 + val word = this_string
2.542 + fun token s =
2.543 +  repeat (\$\$ " ") |-- word s --| repeat (\$\$ " ")
2.544 + val numeral = one isnum
2.545 + val decimalint = bulk numeral >> (rat_of_string o implode)
2.546 + val decimalfrac = bulk numeral
2.547 +    >> (fn s => rat_of_string(implode s) // pow10 (length s))
2.548 + val decimalsig =
2.549 +    decimalint -- option (Scan.\$\$ "." |-- decimalfrac)
2.550 +    >> (fn (h,NONE) => h | (h,SOME x) => h +/ x)
2.551 + fun signed prs =
2.552 +       \$\$ "-" |-- prs >> Rat.neg
2.553 +    || \$\$ "+" |-- prs
2.554 +    || prs;
2.555 +
2.556 +fun emptyin def xs = if null xs then (def,xs) else Scan.fail xs
2.557 +
2.558 + val exponent = (\$\$ "e" || \$\$ "E") |-- signed decimalint;
2.559 +
2.560 + val decimal = signed decimalsig -- (emptyin rat_0|| exponent)
2.561 +    >> (fn (h, x) => h */ pow10 (int_of_rat x));
2.562 +end;
2.563 +
2.564 + fun mkparser p s =
2.565 +  let val (x,rst) = p (explode s)
2.566 +  in if null rst then x
2.567 +     else error "mkparser: unparsed input"
2.568 +  end;;
2.569 +val parse_decimal = mkparser decimal;
2.570 +
2.571 +fun fix err prs =
2.572 +  prs || (fn x=> error err);
2.573 +
2.574 +fun listof prs sep err =
2.575 +  prs -- Scan.bulk (sep |-- fix err prs) >> uncurry cons;
2.576 +
2.577 +(* Parse back a vector.                                                      *)
2.578 +
2.579 + val vector =
2.580 +    token "{" |-- listof decimal (token ",") "decimal" --| token "}"
2.581 +               >> vector_of_list
2.582 + val parse_vector = mkparser vector
2.583 + fun skipupto dscr prs inp =
2.584 +   (dscr |-- prs
2.585 +    || Scan.one (K true) |-- skipupto dscr prs) inp
2.586 + fun ignore inp = ((),[])
2.587 + fun sdpaoutput inp =  skipupto (word "xVec" -- token "=")
2.588 +             (vector --| ignore) inp
2.589 + fun csdpoutput inp =  ((decimal -- Scan.bulk (Scan.\$\$ " " |-- Scan.option decimal) >> (fn (h,to) => map_filter I ((SOME h)::to))) --| ignore >> vector_of_list) inp
2.590 + val parse_sdpaoutput = mkparser sdpaoutput
2.591 + val parse_csdpoutput = mkparser csdpoutput
2.592 +
2.593 +(* Also parse the SDPA output to test success (CSDP yields a return code).   *)
2.594 +
2.595 +local
2.596 + val prs =
2.597 +  skipupto (word "phase.value" -- token "=")
2.598 +   (Scan.option (Scan.\$\$ "p") -- Scan.option (Scan.\$\$ "d")
2.599 +    -- (word "OPT" || word "FEAS"))
2.600 +in
2.601 + fun sdpa_run_succeeded s =
2.602 +  (prs (explode s); true) handle _ => false
2.603 +end;
2.604 +
2.605 +(* The default parameters. Unfortunately this goes to a fixed file.          *)
2.606 +
2.607 +val sdpa_default_parameters =
2.608 +"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";;
2.609 +
2.610 +(* These were suggested by Makoto Yamashita for problems where we are        *)
2.611 +(* right at the edge of the semidefinite cone, as sometimes happens.         *)
2.612 +
2.613 +val sdpa_alt_parameters =
2.614 +"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";;
2.615 +
2.616 +val sdpa_params = sdpa_alt_parameters;;
2.617 +
2.618 +(* CSDP parameters; so far I'm sticking with the defaults.                   *)
2.619 +
2.620 +val csdp_default_parameters =
2.621 +"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";;
2.622 +
2.623 +val csdp_params = csdp_default_parameters;;
2.624 +
2.625 +fun tmp_file pre suf =
2.626 + let val i = string_of_int (round (random()))
2.627 +   val name = Path.append (Path.variable "ISABELLE_TMP") (Path.explode (pre ^ i ^ suf))
2.628 + in
2.629 +   if File.exists name then tmp_file pre suf
2.630 +   else name
2.631 + end;
2.632 +
2.633 +(* Now call SDPA on a problem and parse back the output.                     *)
2.634 +
2.635 +fun run_sdpa dbg obj mats =
2.636 + let
2.637 +  val input_file = tmp_file "sos" ".dat-s"
2.638 +  val output_file = tmp_file "sos" ".out"
2.639 +  val params_file = tmp_file "param" ".sdpa"
2.640 +  val current_dir = File.pwd()
2.641 +  val _ = File.write input_file
2.642 +                         (sdpa_of_problem "" obj mats)
2.643 +  val _ = File.write params_file sdpa_params
2.644 +  val _ = File.cd (Path.variable "ISABELLE_TMP")
2.645 +  val _ = File.system_command ("sdpa "^ (Path.implode input_file) ^ " " ^
2.646 +                               (Path.implode output_file) ^
2.647 +                               (if dbg then "" else "> /dev/null"))
2.648 +  val opr = File.read output_file
2.649 + in if not(sdpa_run_succeeded opr) then error "sdpa: call failed"
2.650 +    else
2.651 +      let val res = parse_sdpaoutput opr
2.652 +      in ((if dbg then ()
2.653 +           else (File.rm input_file; File.rm output_file ; File.cd current_dir));
2.654 +          res)
2.655 +      end
2.656 + end;
2.657 +
2.658 +fun sdpa obj mats = run_sdpa (!debugging) obj mats;
2.659 +
2.660 +(* The same thing with CSDP.                                                 *)
2.661 +
2.662 +fun run_csdp dbg obj mats =
2.663 + let
2.664 +  val input_file = tmp_file "sos" ".dat-s"
2.665 +  val output_file = tmp_file "sos" ".out"
2.666 +  val params_file = tmp_file "param" ".csdp"
2.667 +  val current_dir = File.pwd()
2.668 +  val _ = File.write input_file (sdpa_of_problem "" obj mats)
2.669 +  val _ = File.write params_file csdp_params
2.670 +  val _ = File.cd (Path.variable "ISABELLE_TMP")
2.671 +  val rv = system ("csdp "^(Path.implode input_file) ^ " "
2.672 +                   ^ (Path.implode output_file) ^
2.673 +                   (if dbg then "" else "> /dev/null"))
2.674 +  val  opr = File.read output_file
2.675 +  val res = parse_csdpoutput opr
2.676 + in
2.677 +    ((if dbg then ()
2.678 +      else (File.rm input_file; File.rm output_file ; File.cd current_dir));
2.679 +     (rv,res))
2.680 + end;
2.681 +
2.682 +fun csdp obj mats =
2.683 + let
2.684 +  val (rv,res) = run_csdp (!debugging) obj mats
2.685 + in
2.686 +   ((if rv = 1 orelse rv = 2 then error "csdp: Problem is infeasible"
2.687 +    else if rv = 3 then writeln "csdp warning: Reduced accuracy"
2.688 +    else if rv <> 0 then error ("csdp: error "^string_of_int rv)
2.689 +    else ());
2.690 +   res)
2.691 + end;
2.692 +
2.693 +(* Try some apparently sensible scaling first. Note that this is purely to   *)
2.694 +(* get a cleaner translation to floating-point, and doesn't affect any of    *)
2.695 +(* the results, in principle. In practice it seems a lot better when there   *)
2.696 +(* are extreme numbers in the original problem.                              *)
2.697 +
2.698 +  (* Version for (int*int) keys *)
2.699 +local
2.700 +  fun max_rat x y = if x </ y then y else x
2.701 +  fun common_denominator fld amat acc =
2.702 +      fld (fn (m,c) => fn a => lcm_rat (denominator_rat c) a) amat acc
2.703 +  fun maximal_element fld amat acc =
2.704 +    fld (fn (m,c) => fn maxa => max_rat maxa (abs_rat c)) amat acc
2.705 +fun float_of_rat x = let val (a,b) = Rat.quotient_of_rat x
2.706 +                     in Real.fromLargeInt a / Real.fromLargeInt b end;
2.707 +in
2.708 +
2.709 +fun pi_scale_then solver (obj:vector)  mats =
2.710 + let
2.711 +  val cd1 = fold_rev (common_denominator Intpairfunc.fold) mats (rat_1)
2.712 +  val cd2 = common_denominator Intfunc.fold (snd obj)  (rat_1)
2.713 +  val mats' = map (Intpairfunc.mapf (fn x => cd1 */ x)) mats
2.714 +  val obj' = vector_cmul cd2 obj
2.715 +  val max1 = fold_rev (maximal_element Intpairfunc.fold) mats' (rat_0)
2.716 +  val max2 = maximal_element Intfunc.fold (snd obj') (rat_0)
2.717 +  val scal1 = pow2 (20 - trunc(Math.ln (float_of_rat max1) / Math.ln 2.0))
2.718 +  val scal2 = pow2 (20 - trunc(Math.ln (float_of_rat max2) / Math.ln 2.0))
2.719 +  val mats'' = map (Intpairfunc.mapf (fn x => x */ scal1)) mats'
2.720 +  val obj'' = vector_cmul scal2 obj'
2.721 + in solver obj'' mats''
2.722 +  end
2.723 +end;
2.724 +
2.725 +(* Try some apparently sensible scaling first. Note that this is purely to   *)
2.726 +(* get a cleaner translation to floating-point, and doesn't affect any of    *)
2.727 +(* the results, in principle. In practice it seems a lot better when there   *)
2.728 +(* are extreme numbers in the original problem.                              *)
2.729 +
2.730 +  (* Version for (int*int*int) keys *)
2.731 +local
2.732 +  fun max_rat x y = if x </ y then y else x
2.733 +  fun common_denominator fld amat acc =
2.734 +      fld (fn (m,c) => fn a => lcm_rat (denominator_rat c) a) amat acc
2.735 +  fun maximal_element fld amat acc =
2.736 +    fld (fn (m,c) => fn maxa => max_rat maxa (abs_rat c)) amat acc
2.737 +fun float_of_rat x = let val (a,b) = Rat.quotient_of_rat x
2.738 +                     in Real.fromLargeInt a / Real.fromLargeInt b end;
2.739 +fun int_of_float x = (trunc x handle Overflow => 0 | Domain => 0)
2.740 +in
2.741 +
2.742 +fun tri_scale_then solver (obj:vector)  mats =
2.743 + let
2.744 +  val cd1 = fold_rev (common_denominator Inttriplefunc.fold) mats (rat_1)
2.745 +  val cd2 = common_denominator Intfunc.fold (snd obj)  (rat_1)
2.746 +  val mats' = map (Inttriplefunc.mapf (fn x => cd1 */ x)) mats
2.747 +  val obj' = vector_cmul cd2 obj
2.748 +  val max1 = fold_rev (maximal_element Inttriplefunc.fold) mats' (rat_0)
2.749 +  val max2 = maximal_element Intfunc.fold (snd obj') (rat_0)
2.750 +  val scal1 = pow2 (20 - int_of_float(Math.ln (float_of_rat max1) / Math.ln 2.0))
2.751 +  val scal2 = pow2 (20 - int_of_float(Math.ln (float_of_rat max2) / Math.ln 2.0))
2.752 +  val mats'' = map (Inttriplefunc.mapf (fn x => x */ scal1)) mats'
2.753 +  val obj'' = vector_cmul scal2 obj'
2.754 + in solver obj'' mats''
2.755 +  end
2.756 +end;
2.757 +
2.758 +(* Round a vector to "nice" rationals.                                       *)
2.759 +
2.760 +fun nice_rational n x = round_rat (n */ x) // n;;
2.761 +fun nice_vector n ((d,v) : vector) =
2.762 + (d, Intfunc.fold (fn (i,c) => fn a =>
2.763 +   let val y = nice_rational n c
2.764 +   in if c =/ rat_0 then a
2.765 +      else Intfunc.update (i,y) a end) v Intfunc.undefined):vector
2.766 +
2.767 +
2.768 +(* Reduce linear program to SDP (diagonal matrices) and test with CSDP. This *)
2.769 +(* one tests A [-1;x1;..;xn] >= 0 (i.e. left column is negated constants).   *)
2.770 +
2.771 +fun linear_program_basic a =
2.772 + let
2.773 +  val (m,n) = dimensions a
2.774 +  val mats =  map (fn j => diagonal (column j a)) (1 upto n)
2.775 +  val obj = vector_const rat_1 m
2.776 +  val (rv,res) = run_csdp false obj mats
2.777 + in if rv = 1 orelse rv = 2 then false
2.778 +    else if rv = 0 then true
2.779 +    else error "linear_program: An error occurred in the SDP solver"
2.780 + end;
2.781 +
2.782 +(* Alternative interface testing A x >= b for matrix A, vector b.            *)
2.783 +
2.784 +fun linear_program a b =
2.785 + let val (m,n) = dimensions a
2.786 + in if dim b <> m then error "linear_program: incompatible dimensions"
2.787 +    else
2.788 +    let
2.789 +     val mats = diagonal b :: map (fn j => diagonal (column j a)) (1 upto n)
2.790 +     val obj = vector_const rat_1 m
2.791 +     val (rv,res) = run_csdp false obj mats
2.792 +    in if rv = 1 orelse rv = 2 then false
2.793 +       else if rv = 0 then true
2.794 +       else error "linear_program: An error occurred in the SDP solver"
2.795 +    end
2.796 + end;
2.797 +
2.798 +(* Test whether a point is in the convex hull of others. Rather than use     *)
2.799 +(* computational geometry, express as linear inequalities and call CSDP.     *)
2.800 +(* This is a bit lazy of me, but it's easy and not such a bottleneck so far. *)
2.801 +
2.802 +fun in_convex_hull pts pt =
2.803 + let
2.804 +  val pts1 = (1::pt) :: map (fn x => 1::x) pts
2.805 +  val pts2 = map (fn p => map (fn x => ~x) p @ p) pts1
2.806 +  val n = length pts + 1
2.807 +  val v = 2 * (length pt + 1)
2.808 +  val m = v + n - 1
2.809 +  val mat = ((m,n),
2.810 +  itern 1 pts2 (fn pts => fn j => itern 1 pts
2.811 +               (fn x => fn i => Intpairfunc.update ((i,j), Rat.rat_of_int x)))
2.812 +  (iter (1,n) (fn i => Intpairfunc.update((v + i,i+1), rat_1))
2.813 +      Intpairfunc.undefined))
2.814 + in linear_program_basic mat
2.815 + end;
2.816 +
2.817 +(* Filter down a set of points to a minimal set with the same convex hull.   *)
2.818 +
2.819 +local
2.820 + fun augment1 (m::ms) = if in_convex_hull ms m then ms else ms@[m]
2.821 + fun augment m ms = funpow 3 augment1 (m::ms)
2.822 +in
2.823 +fun minimal_convex_hull mons =
2.824 + let val mons' = fold_rev augment (tl mons) [hd mons]
2.825 + in funpow (length mons') augment1 mons'
2.826 + end
2.827 +end;
2.828 +
2.829 +fun dest_ord f x = is_equal (f x);
2.830 +
2.831 +(* Stuff for "equations" ((int*int*int)->num functions).                         *)
2.832 +
2.833 +fun tri_equation_cmul c eq =
2.834 +  if c =/ rat_0 then Inttriplefunc.undefined else Inttriplefunc.mapf (fn d => c */ d) eq;
2.835 +
2.836 +fun tri_equation_add eq1 eq2 = Inttriplefunc.combine (curry op +/) (fn x => x =/ rat_0) eq1 eq2;
2.837 +
2.838 +fun tri_equation_eval assig eq =
2.839 + let fun value v = Inttriplefunc.apply assig v
2.840 + in Inttriplefunc.fold (fn (v, c) => fn a => a +/ value v */ c) eq rat_0
2.841 + end;
2.842 +
2.843 +(* Eliminate among linear equations: return unconstrained variables and      *)
2.844 +(* assignments for the others in terms of them. We give one pseudo-variable  *)
2.845 +(* "one" that's used for a constant term.                                    *)
2.846 +
2.847 +local
2.848 +  fun extract_first p l = case l of  (* FIXME : use find_first instead *)
2.849 +   [] => error "extract_first"
2.850 + | h::t => if p h then (h,t) else
2.851 +          let val (k,s) = extract_first p t in (k,h::s) end
2.852 +fun eliminate vars dun eqs = case vars of
2.853 +  [] => if forall Inttriplefunc.is_undefined eqs then dun
2.854 +        else raise Unsolvable
2.855 + | v::vs =>
2.856 +  ((let
2.857 +    val (eq,oeqs) = extract_first (fn e => Inttriplefunc.defined e v) eqs
2.858 +    val a = Inttriplefunc.apply eq v
2.859 +    val eq' = tri_equation_cmul ((Rat.neg rat_1) // a) (Inttriplefunc.undefine v eq)
2.860 +    fun elim e =
2.861 +     let val b = Inttriplefunc.tryapplyd e v rat_0
2.862 +     in if b =/ rat_0 then e else
2.863 +        tri_equation_add e (tri_equation_cmul (Rat.neg b // a) eq)
2.864 +     end
2.865 +   in eliminate vs (Inttriplefunc.update (v,eq') (Inttriplefunc.mapf elim dun)) (map elim oeqs)
2.866 +   end)
2.867 +  handle ERROR _ => eliminate vs dun eqs)
2.868 +in
2.869 +fun tri_eliminate_equations one vars eqs =
2.870 + let
2.871 +  val assig = eliminate vars Inttriplefunc.undefined eqs
2.872 +  val vs = Inttriplefunc.fold (fn (x, f) => fn a => remove (dest_ord triple_int_ord) one (Inttriplefunc.dom f) @ a) assig []
2.873 +  in (distinct (dest_ord triple_int_ord) vs, assig)
2.874 +  end
2.875 +end;
2.876 +
2.877 +(* Eliminate all variables, in an essentially arbitrary order.               *)
2.878 +
2.879 +fun tri_eliminate_all_equations one =
2.880 + let
2.881 +  fun choose_variable eq =
2.882 +   let val (v,_) = Inttriplefunc.choose eq
2.883 +   in if is_equal (triple_int_ord(v,one)) then
2.884 +      let val eq' = Inttriplefunc.undefine v eq
2.885 +      in if Inttriplefunc.is_undefined eq' then error "choose_variable"
2.886 +         else fst (Inttriplefunc.choose eq')
2.887 +      end
2.888 +    else v
2.889 +   end
2.890 +  fun eliminate dun eqs = case eqs of
2.891 +    [] => dun
2.892 +  | eq::oeqs =>
2.893 +    if Inttriplefunc.is_undefined eq then eliminate dun oeqs else
2.894 +    let val v = choose_variable eq
2.895 +        val a = Inttriplefunc.apply eq v
2.896 +        val eq' = tri_equation_cmul ((Rat.rat_of_int ~1) // a)
2.897 +                   (Inttriplefunc.undefine v eq)
2.898 +        fun elim e =
2.899 +         let val b = Inttriplefunc.tryapplyd e v rat_0
2.900 +         in if b =/ rat_0 then e
2.901 +            else tri_equation_add e (tri_equation_cmul (Rat.neg b // a) eq)
2.902 +         end
2.903 +    in eliminate (Inttriplefunc.update(v, eq') (Inttriplefunc.mapf elim dun))
2.904 +                 (map elim oeqs)
2.905 +    end
2.906 +in fn eqs =>
2.907 + let
2.908 +  val assig = eliminate Inttriplefunc.undefined eqs
2.909 +  val vs = Inttriplefunc.fold (fn (x, f) => fn a => remove (dest_ord triple_int_ord) one (Inttriplefunc.dom f) @ a) assig []
2.910 + in (distinct (dest_ord triple_int_ord) vs,assig)
2.911 + end
2.912 +end;
2.913 +
2.914 +(* Solve equations by assigning arbitrary numbers.                           *)
2.915 +
2.916 +fun tri_solve_equations one eqs =
2.917 + let
2.918 +  val (vars,assigs) = tri_eliminate_all_equations one eqs
2.919 +  val vfn = fold_rev (fn v => Inttriplefunc.update(v,rat_0)) vars
2.920 +            (Inttriplefunc.onefunc(one, Rat.rat_of_int ~1))
2.921 +  val ass =
2.922 +    Inttriplefunc.combine (curry op +/) (K false)
2.923 +    (Inttriplefunc.mapf (tri_equation_eval vfn) assigs) vfn
2.924 + in if forall (fn e => tri_equation_eval ass e =/ rat_0) eqs
2.925 +    then Inttriplefunc.undefine one ass else raise Sanity
2.926 + end;
2.927 +
2.928 +(* Multiply equation-parametrized poly by regular poly and add accumulator.  *)
2.929 +
2.930 +fun tri_epoly_pmul p q acc =
2.931 + Monomialfunc.fold (fn (m1, c) => fn a =>
2.932 +  Monomialfunc.fold (fn (m2,e) => fn b =>
2.933 +   let val m =  monomial_mul m1 m2
2.934 +       val es = Monomialfunc.tryapplyd b m Inttriplefunc.undefined
2.935 +   in Monomialfunc.update (m,tri_equation_add (tri_equation_cmul c e) es) b
2.936 +   end) q a) p acc ;
2.937 +
2.938 +(* Usual operations on equation-parametrized poly.                           *)
2.939 +
2.940 +fun tri_epoly_cmul c l =
2.941 +  if c =/ rat_0 then Inttriplefunc.undefined else Inttriplefunc.mapf (tri_equation_cmul c) l;;
2.942 +
2.943 +val tri_epoly_neg = tri_epoly_cmul (Rat.rat_of_int ~1);
2.944 +
2.945 +val tri_epoly_add = Inttriplefunc.combine tri_equation_add Inttriplefunc.is_undefined;
2.946 +
2.947 +fun tri_epoly_sub p q = tri_epoly_add p (tri_epoly_neg q);;
2.948 +
2.949 +(* Stuff for "equations" ((int*int)->num functions).                         *)
2.950 +
2.951 +fun pi_equation_cmul c eq =
2.952 +  if c =/ rat_0 then Inttriplefunc.undefined else Inttriplefunc.mapf (fn d => c */ d) eq;
2.953 +
2.954 +fun pi_equation_add eq1 eq2 = Inttriplefunc.combine (curry op +/) (fn x => x =/ rat_0) eq1 eq2;
2.955 +
2.956 +fun pi_equation_eval assig eq =
2.957 + let fun value v = Inttriplefunc.apply assig v
2.958 + in Inttriplefunc.fold (fn (v, c) => fn a => a +/ value v */ c) eq rat_0
2.959 + end;
2.960 +
2.961 +(* Eliminate among linear equations: return unconstrained variables and      *)
2.962 +(* assignments for the others in terms of them. We give one pseudo-variable  *)
2.963 +(* "one" that's used for a constant term.                                    *)
2.964 +
2.965 +local
2.966 +fun extract_first p l = case l of
2.967 +   [] => error "extract_first"
2.968 + | h::t => if p h then (h,t) else
2.969 +          let val (k,s) = extract_first p t in (k,h::s) end
2.970 +fun eliminate vars dun eqs = case vars of
2.971 +  [] => if forall Inttriplefunc.is_undefined eqs then dun
2.972 +        else raise Unsolvable
2.973 + | v::vs =>
2.974 +   let
2.975 +    val (eq,oeqs) = extract_first (fn e => Inttriplefunc.defined e v) eqs
2.976 +    val a = Inttriplefunc.apply eq v
2.977 +    val eq' = pi_equation_cmul ((Rat.neg rat_1) // a) (Inttriplefunc.undefine v eq)
2.978 +    fun elim e =
2.979 +     let val b = Inttriplefunc.tryapplyd e v rat_0
2.980 +     in if b =/ rat_0 then e else
2.981 +        pi_equation_add e (pi_equation_cmul (Rat.neg b // a) eq)
2.982 +     end
2.983 +   in eliminate vs (Inttriplefunc.update (v,eq') (Inttriplefunc.mapf elim dun)) (map elim oeqs)
2.984 +   end
2.985 +  handle ERROR _ => eliminate vs dun eqs
2.986 +in
2.987 +fun pi_eliminate_equations one vars eqs =
2.988 + let
2.989 +  val assig = eliminate vars Inttriplefunc.undefined eqs
2.990 +  val vs = Inttriplefunc.fold (fn (x, f) => fn a => remove (dest_ord triple_int_ord) one (Inttriplefunc.dom f) @ a) assig []
2.991 +  in (distinct (dest_ord triple_int_ord) vs, assig)
2.992 +  end
2.993 +end;
2.994 +
2.995 +(* Eliminate all variables, in an essentially arbitrary order.               *)
2.996 +
2.997 +fun pi_eliminate_all_equations one =
2.998 + let
2.999 +  fun choose_variable eq =
2.1000 +   let val (v,_) = Inttriplefunc.choose eq
2.1001 +   in if is_equal (triple_int_ord(v,one)) then
2.1002 +      let val eq' = Inttriplefunc.undefine v eq
2.1003 +      in if Inttriplefunc.is_undefined eq' then error "choose_variable"
2.1004 +         else fst (Inttriplefunc.choose eq')
2.1005 +      end
2.1006 +    else v
2.1007 +   end
2.1008 +  fun eliminate dun eqs = case eqs of
2.1009 +    [] => dun
2.1010 +  | eq::oeqs =>
2.1011 +    if Inttriplefunc.is_undefined eq then eliminate dun oeqs else
2.1012 +    let val v = choose_variable eq
2.1013 +        val a = Inttriplefunc.apply eq v
2.1014 +        val eq' = pi_equation_cmul ((Rat.rat_of_int ~1) // a)
2.1015 +                   (Inttriplefunc.undefine v eq)
2.1016 +        fun elim e =
2.1017 +         let val b = Inttriplefunc.tryapplyd e v rat_0
2.1018 +         in if b =/ rat_0 then e
2.1019 +            else pi_equation_add e (pi_equation_cmul (Rat.neg b // a) eq)
2.1020 +         end
2.1021 +    in eliminate (Inttriplefunc.update(v, eq') (Inttriplefunc.mapf elim dun))
2.1022 +                 (map elim oeqs)
2.1023 +    end
2.1024 +in fn eqs =>
2.1025 + let
2.1026 +  val assig = eliminate Inttriplefunc.undefined eqs
2.1027 +  val vs = Inttriplefunc.fold (fn (x, f) => fn a => remove (dest_ord triple_int_ord) one (Inttriplefunc.dom f) @ a) assig []
2.1028 + in (distinct (dest_ord triple_int_ord) vs,assig)
2.1029 + end
2.1030 +end;
2.1031 +
2.1032 +(* Solve equations by assigning arbitrary numbers.                           *)
2.1033 +
2.1034 +fun pi_solve_equations one eqs =
2.1035 + let
2.1036 +  val (vars,assigs) = pi_eliminate_all_equations one eqs
2.1037 +  val vfn = fold_rev (fn v => Inttriplefunc.update(v,rat_0)) vars
2.1038 +            (Inttriplefunc.onefunc(one, Rat.rat_of_int ~1))
2.1039 +  val ass =
2.1040 +    Inttriplefunc.combine (curry op +/) (K false)
2.1041 +    (Inttriplefunc.mapf (pi_equation_eval vfn) assigs) vfn
2.1042 + in if forall (fn e => pi_equation_eval ass e =/ rat_0) eqs
2.1043 +    then Inttriplefunc.undefine one ass else raise Sanity
2.1044 + end;
2.1045 +
2.1046 +(* Multiply equation-parametrized poly by regular poly and add accumulator.  *)
2.1047 +
2.1048 +fun pi_epoly_pmul p q acc =
2.1049 + Monomialfunc.fold (fn (m1, c) => fn a =>
2.1050 +  Monomialfunc.fold (fn (m2,e) => fn b =>
2.1051 +   let val m =  monomial_mul m1 m2
2.1052 +       val es = Monomialfunc.tryapplyd b m Inttriplefunc.undefined
2.1053 +   in Monomialfunc.update (m,pi_equation_add (pi_equation_cmul c e) es) b
2.1054 +   end) q a) p acc ;
2.1055 +
2.1056 +(* Usual operations on equation-parametrized poly.                           *)
2.1057 +
2.1058 +fun pi_epoly_cmul c l =
2.1059 +  if c =/ rat_0 then Inttriplefunc.undefined else Inttriplefunc.mapf (pi_equation_cmul c) l;;
2.1060 +
2.1061 +val pi_epoly_neg = pi_epoly_cmul (Rat.rat_of_int ~1);
2.1062 +
2.1063 +val pi_epoly_add = Inttriplefunc.combine pi_equation_add Inttriplefunc.is_undefined;
2.1064 +
2.1065 +fun pi_epoly_sub p q = pi_epoly_add p (pi_epoly_neg q);;
2.1066 +
2.1067 +fun allpairs f l1 l2 =  fold_rev (fn x => (curry (op @)) (map (f x) l2)) l1 [];
2.1068 +
2.1069 +(* Hence produce the "relevant" monomials: those whose squares lie in the    *)
2.1070 +(* Newton polytope of the monomials in the input. (This is enough according  *)
2.1071 +(* to Reznik: "Extremal PSD forms with few terms", Duke Math. Journal,       *)
2.1072 +(* vol 45, pp. 363--374, 1978.                                               *)
2.1073 +(*                                                                           *)
2.1074 +(* These are ordered in sort of decreasing degree. In particular the         *)
2.1075 +(* constant monomial is last; this gives an order in diagonalization of the  *)
2.1076 +(* quadratic form that will tend to display constants.                       *)
2.1077 +
2.1078 +fun newton_polytope pol =
2.1079 + let
2.1080 +  val vars = poly_variables pol
2.1081 +  val mons = map (fn m => map (fn x => monomial_degree x m) vars)
2.1082 +             (Monomialfunc.dom pol)
2.1083 +  val ds = map (fn x => (degree x pol + 1) div 2) vars
2.1084 +  val all = fold_rev (fn n => allpairs cons (0 upto n)) ds [[]]
2.1085 +  val mons' = minimal_convex_hull mons
2.1086 +  val all' =
2.1087 +    filter (fn m => in_convex_hull mons' (map (fn x => 2 * x) m)) all
2.1088 + in map (fn m => fold_rev2 (fn v => fn i => fn a => if i = 0 then a else Ctermfunc.update (v,i) a)
2.1089 +                        vars m monomial_1) (rev all')
2.1090 + end;
2.1091 +
2.1092 +(* Diagonalize (Cholesky/LDU) the matrix corresponding to a quadratic form.  *)
2.1093 +
2.1094 +local
2.1095 +fun diagonalize n i m =
2.1096 + if Intpairfunc.is_undefined (snd m) then []
2.1097 + else
2.1098 +  let val a11 = Intpairfunc.tryapplyd (snd m) (i,i) rat_0
2.1099 +  in if a11 </ rat_0 then error "diagonalize: not PSD"
2.1100 +    else if a11 =/ rat_0 then
2.1101 +          if Intfunc.is_undefined (snd (row i m)) then diagonalize n (i + 1) m
2.1102 +          else error "diagonalize: not PSD ___ "
2.1103 +    else
2.1104 +     let
2.1105 +      val v = row i m
2.1106 +      val v' = (fst v, Intfunc.fold (fn (i, c) => fn a =>
2.1107 +       let val y = c // a11
2.1108 +       in if y = rat_0 then a else Intfunc.update (i,y) a
2.1109 +       end)  (snd v) Intfunc.undefined)
2.1110 +      fun upt0 x y a = if y = rat_0 then a else Intpairfunc.update (x,y) a
2.1111 +      val m' =
2.1112 +      ((n,n),
2.1113 +      iter (i+1,n) (fn j =>
2.1114 +          iter (i+1,n) (fn k =>
2.1115 +              (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))))
2.1116 +          Intpairfunc.undefined)
2.1117 +     in (a11,v')::diagonalize n (i + 1) m'
2.1118 +     end
2.1119 +  end
2.1120 +in
2.1121 +fun diag m =
2.1122 + let
2.1123 +   val nn = dimensions m
2.1124 +   val n = fst nn
2.1125 + in if snd nn <> n then error "diagonalize: non-square matrix"
2.1126 +    else diagonalize n 1 m
2.1127 + end
2.1128 +end;
2.1129 +
2.1130 +fun gcd_rat a b = Rat.rat_of_int (Integer.gcd (int_of_rat a) (int_of_rat b));
2.1131 +
2.1132 +(* Adjust a diagonalization to collect rationals at the start.               *)
2.1133 +  (* FIXME : Potentially polymorphic keys, but here only: integers!! *)
2.1134 +local
2.1135 + fun upd0 x y a = if y =/ rat_0 then a else Intfunc.update(x,y) a;
2.1136 + fun mapa f (d,v) =
2.1137 +  (d, Intfunc.fold (fn (i,c) => fn a => upd0 i (f c) a) v Intfunc.undefined)
2.1138 + fun adj (c,l) =
2.1139 + let val a =
2.1140 +  Intfunc.fold (fn (i,c) => fn a => lcm_rat a (denominator_rat c))
2.1141 +    (snd l) rat_1 //
2.1142 +  Intfunc.fold (fn (i,c) => fn a => gcd_rat a (numerator_rat c))
2.1143 +    (snd l) rat_0
2.1144 +  in ((c // (a */ a)),mapa (fn x => a */ x) l)
2.1145 +  end
2.1146 +in
2.1147 +fun deration d = if null d then (rat_0,d) else
2.1148 + let val d' = map adj d
2.1149 +     val a = fold (lcm_rat o denominator_rat o fst) d' rat_1 //
2.1150 +          fold (gcd_rat o numerator_rat o fst) d' rat_0
2.1151 + in ((rat_1 // a),map (fn (c,l) => (a */ c,l)) d')
2.1152 + end
2.1153 +end;
2.1154 +
2.1155 +(* Enumeration of monomials with given multidegree bound.                    *)
2.1156 +
2.1157 +fun enumerate_monomials d vars =
2.1158 + if d < 0 then []
2.1159 + else if d = 0 then [Ctermfunc.undefined]
2.1160 + else if null vars then [monomial_1] else
2.1161 + let val alts =
2.1162 +  map (fn k => let val oths = enumerate_monomials (d - k) (tl vars)
2.1163 +               in map (fn ks => if k = 0 then ks else Ctermfunc.update (hd vars, k) ks) oths end) (0 upto d)
2.1164 + in fold1 (curry op @) alts
2.1165 + end;
2.1166 +
2.1167 +(* Enumerate products of distinct input polys with degree <= d.              *)
2.1168 +(* We ignore any constant input polynomials.                                 *)
2.1169 +(* Give the output polynomial and a record of how it was derived.            *)
2.1170 +
2.1171 +local
2.1172 + open RealArith
2.1173 +in
2.1174 +fun enumerate_products d pols =
2.1175 +if d = 0 then [(poly_const rat_1,Rational_lt rat_1)]
2.1176 +else if d < 0 then [] else
2.1177 +case pols of
2.1178 +   [] => [(poly_const rat_1,Rational_lt rat_1)]
2.1179 + | (p,b)::ps =>
2.1180 +    let val e = multidegree p
2.1181 +    in if e = 0 then enumerate_products d ps else
2.1182 +       enumerate_products d ps @
2.1183 +       map (fn (q,c) => (poly_mul p q,Product(b,c)))
2.1184 +         (enumerate_products (d - e) ps)
2.1185 +    end
2.1186 +end;
2.1187 +
2.1188 +(* Convert regular polynomial. Note that we treat (0,0,0) as -1.             *)
2.1189 +
2.1190 +fun epoly_of_poly p =
2.1191 +  Monomialfunc.fold (fn (m,c) => fn a => Monomialfunc.update (m, Inttriplefunc.onefunc ((0,0,0), Rat.neg c)) a) p Monomialfunc.undefined;
2.1192 +
2.1193 +(* String for block diagonal matrix numbered k.                              *)
2.1194 +
2.1195 +fun sdpa_of_blockdiagonal k m =
2.1196 + let
2.1197 +  val pfx = string_of_int k ^" "
2.1198 +  val ents =
2.1199 +    Inttriplefunc.fold
2.1200 +      (fn ((b,i,j),c) => fn a => if i > j then a else ((b,i,j),c)::a)
2.1201 +      m []
2.1202 +  val entss = sort (increasing fst triple_int_ord) ents
2.1203 + in fold_rev (fn ((b,i,j),c) => fn a =>
2.1204 +     pfx ^ string_of_int b ^ " " ^ string_of_int i ^ " " ^ string_of_int j ^
2.1205 +     " " ^ decimalize 20 c ^ "\n" ^ a) entss ""
2.1206 + end;
2.1207 +
2.1208 +(* SDPA for problem using block diagonal (i.e. multiple SDPs)                *)
2.1209 +
2.1210 +fun sdpa_of_blockproblem comment nblocks blocksizes obj mats =
2.1211 + let val m = length mats - 1
2.1212 + in  "\"" ^ comment ^ "\"\n" ^
2.1213 +  string_of_int m ^ "\n" ^
2.1214 +  string_of_int nblocks ^ "\n" ^
2.1215 +  (fold1 (fn s => fn t => s^" "^t) (map string_of_int blocksizes)) ^
2.1216 +  "\n" ^
2.1217 +  sdpa_of_vector obj ^
2.1218 +  fold_rev2 (fn k => fn m => fn a => sdpa_of_blockdiagonal (k - 1) m ^ a)
2.1219 +    (1 upto length mats) mats ""
2.1220 + end;
2.1221 +
2.1222 +(* Hence run CSDP on a problem in block diagonal form.                       *)
2.1223 +
2.1224 +fun run_csdp dbg nblocks blocksizes obj mats =
2.1225 + let
2.1226 +  val input_file = tmp_file "sos" ".dat-s"
2.1227 +  val output_file = tmp_file "sos" ".out"
2.1228 +  val params_file = tmp_file "param" ".csdp"
2.1229 +  val _ = File.write input_file
2.1230 +   (sdpa_of_blockproblem "" nblocks blocksizes obj mats)
2.1231 +  val _ = File.write params_file csdp_params
2.1232 +  val current_dir = File.pwd()
2.1233 +  val _ = File.cd (Path.variable "ISABELLE_TMP")
2.1234 +  val rv = system ("csdp "^(Path.implode input_file) ^ " "
2.1235 +                   ^ (Path.implode output_file) ^
2.1236 +                   (if dbg then "" else "> /dev/null"))
2.1237 +  val  opr = File.read output_file
2.1238 +  val res = parse_csdpoutput opr
2.1239 + in
2.1240 +   ((if dbg then ()
2.1241 +     else (File.rm input_file ; File.rm output_file ; File.cd current_dir));
2.1242 +    (rv,res))
2.1243 + end;
2.1244 +
2.1245 +fun csdp nblocks blocksizes obj mats =
2.1246 + let
2.1247 +  val (rv,res) = run_csdp (!debugging) nblocks blocksizes obj mats
2.1248 + in ((if rv = 1 orelse rv = 2 then error "csdp: Problem is infeasible"
2.1249 +     else if rv = 3 then writeln "csdp warning: Reduced accuracy"
2.1250 +     else if rv <> 0 then error ("csdp: error "^string_of_int rv)
2.1251 +     else ());
2.1252 +     res)
2.1253 + end;
2.1254 +
2.1255 +(* 3D versions of matrix operations to consider blocks separately.           *)
2.1256 +
2.1257 +val bmatrix_add = Inttriplefunc.combine (curry op +/) (fn x => x =/ rat_0);
2.1258 +fun bmatrix_cmul c bm =
2.1259 +  if c =/ rat_0 then Inttriplefunc.undefined
2.1260 +  else Inttriplefunc.mapf (fn x => c */ x) bm;
2.1261 +
2.1262 +val bmatrix_neg = bmatrix_cmul (Rat.rat_of_int ~1);
2.1263 +fun bmatrix_sub m1 m2 = bmatrix_add m1 (bmatrix_neg m2);;
2.1264 +
2.1265 +(* Smash a block matrix into components.                                     *)
2.1266 +
2.1267 +fun blocks blocksizes bm =
2.1268 + map (fn (bs,b0) =>
2.1269 +      let val m = Inttriplefunc.fold
2.1270 +          (fn ((b,i,j),c) => fn a => if b = b0 then Intpairfunc.update ((i,j),c) a else a) bm Intpairfunc.undefined
2.1271 +          val d = Intpairfunc.fold (fn ((i,j),c) => fn a => max a (max i j)) m 0
2.1272 +      in (((bs,bs),m):matrix) end)
2.1273 + (blocksizes ~~ (1 upto length blocksizes));;
2.1274 +
2.1275 +(* FIXME : Get rid of this !!!*)
2.1276 +fun tryfind f [] = error "tryfind"
2.1277 +  | tryfind f (x::xs) = (f x handle ERROR _ => tryfind f xs);
2.1278 +
2.1279 +
2.1280 +(* Positiv- and Nullstellensatz. Flag "linf" forces a linear representation. *)
2.1281 +
2.1282 +
2.1283 +local
2.1284 + open RealArith
2.1285 +in
2.1286 +fun real_positivnullstellensatz_general linf d eqs leqs pol =
2.1287 +let
2.1288 + val vars = fold_rev (curry (gen_union (op aconvc)) o poly_variables)
2.1289 +              (pol::eqs @ map fst leqs) []
2.1290 + val monoid = if linf then
2.1291 +      (poly_const rat_1,Rational_lt rat_1)::
2.1292 +      (filter (fn (p,c) => multidegree p <= d) leqs)
2.1293 +    else enumerate_products d leqs
2.1294 + val nblocks = length monoid
2.1295 + fun mk_idmultiplier k p =
2.1296 +  let
2.1297 +   val e = d - multidegree p
2.1298 +   val mons = enumerate_monomials e vars
2.1299 +   val nons = mons ~~ (1 upto length mons)
2.1300 +  in (mons,
2.1301 +      fold_rev (fn (m,n) => Monomialfunc.update(m,Inttriplefunc.onefunc((~k,~n,n),rat_1))) nons Monomialfunc.undefined)
2.1302 +  end
2.1303 +
2.1304 + fun mk_sqmultiplier k (p,c) =
2.1305 +  let
2.1306 +   val e = (d - multidegree p) div 2
2.1307 +   val mons = enumerate_monomials e vars
2.1308 +   val nons = mons ~~ (1 upto length mons)
2.1309 +  in (mons,
2.1310 +      fold_rev (fn (m1,n1) =>
2.1311 +       fold_rev (fn (m2,n2) => fn  a =>
2.1312 +        let val m = monomial_mul m1 m2
2.1313 +        in if n1 > n2 then a else
2.1314 +          let val c = if n1 = n2 then rat_1 else rat_2
2.1315 +              val e = Monomialfunc.tryapplyd a m Inttriplefunc.undefined
2.1316 +          in Monomialfunc.update(m, tri_equation_add (Inttriplefunc.onefunc((k,n1,n2), c)) e) a
2.1317 +          end
2.1318 +        end)  nons)
2.1319 +       nons Monomialfunc.undefined)
2.1320 +  end
2.1321 +
2.1322 +  val (sqmonlist,sqs) = split_list (map2 mk_sqmultiplier (1 upto length monoid) monoid)
2.1323 +  val (idmonlist,ids) =  split_list(map2 mk_idmultiplier (1 upto length eqs) eqs)
2.1324 +  val blocksizes = map length sqmonlist
2.1325 +  val bigsum =
2.1326 +    fold_rev2 (fn p => fn q => fn a => tri_epoly_pmul p q a) eqs ids
2.1327 +            (fold_rev2 (fn (p,c) => fn s => fn a => tri_epoly_pmul p s a) monoid sqs
2.1328 +                     (epoly_of_poly(poly_neg pol)))
2.1329 +  val eqns = Monomialfunc.fold (fn (m,e) => fn a => e::a) bigsum []
2.1330 +  val (pvs,assig) = tri_eliminate_all_equations (0,0,0) eqns
2.1331 +  val qvars = (0,0,0)::pvs
2.1332 +  val allassig = fold_rev (fn v => Inttriplefunc.update(v,(Inttriplefunc.onefunc(v,rat_1)))) pvs assig
2.1333 +  fun mk_matrix v =
2.1334 +    Inttriplefunc.fold (fn ((b,i,j), ass) => fn m =>
2.1335 +        if b < 0 then m else
2.1336 +         let val c = Inttriplefunc.tryapplyd ass v rat_0
2.1337 +         in if c = rat_0 then m else
2.1338 +            Inttriplefunc.update ((b,j,i), c) (Inttriplefunc.update ((b,i,j), c) m)
2.1339 +         end)
2.1340 +          allassig Inttriplefunc.undefined
2.1341 +  val diagents = Inttriplefunc.fold
2.1342 +    (fn ((b,i,j), e) => fn a => if b > 0 andalso i = j then tri_equation_add e a else a)
2.1343 +    allassig Inttriplefunc.undefined
2.1344 +
2.1345 +  val mats = map mk_matrix qvars
2.1346 +  val obj = (length pvs,
2.1347 +            itern 1 pvs (fn v => fn i => Intfunc.updatep iszero (i,Inttriplefunc.tryapplyd diagents v rat_0))
2.1348 +                        Intfunc.undefined)
2.1349 +  val raw_vec = if null pvs then vector_0 0
2.1350 +                else tri_scale_then (csdp nblocks blocksizes) obj mats
2.1351 +  fun int_element (d,v) i = Intfunc.tryapplyd v i rat_0
2.1352 +  fun cterm_element (d,v) i = Ctermfunc.tryapplyd v i rat_0
2.1353 +
2.1354 +  fun find_rounding d =
2.1355 +   let
2.1356 +    val _ = if !debugging
2.1357 +          then writeln ("Trying rounding with limit "^Rat.string_of_rat d ^ "\n")
2.1358 +          else ()
2.1359 +    val vec = nice_vector d raw_vec
2.1360 +    val blockmat = iter (1,dim vec)
2.1361 +     (fn i => fn a => bmatrix_add (bmatrix_cmul (int_element vec i) (nth mats i)) a)
2.1362 +     (bmatrix_neg (nth mats 0))
2.1363 +    val allmats = blocks blocksizes blockmat
2.1364 +   in (vec,map diag allmats)
2.1365 +   end
2.1366 +  val (vec,ratdias) =
2.1367 +    if null pvs then find_rounding rat_1
2.1368 +    else tryfind find_rounding (map Rat.rat_of_int (1 upto 31) @
2.1369 +                                map pow2 (5 upto 66))
2.1370 +  val newassigs =
2.1371 +    fold_rev (fn k => Inttriplefunc.update (nth pvs (k - 1), int_element vec k))
2.1372 +           (1 upto dim vec) (Inttriplefunc.onefunc ((0,0,0), Rat.rat_of_int ~1))
2.1373 +  val finalassigs =
2.1374 +    Inttriplefunc.fold (fn (v,e) => fn a => Inttriplefunc.update(v, tri_equation_eval newassigs e) a) allassig newassigs
2.1375 +  fun poly_of_epoly p =
2.1376 +    Monomialfunc.fold (fn (v,e) => fn a => Monomialfunc.updatep iszero (v,tri_equation_eval finalassigs e) a)
2.1377 +          p Monomialfunc.undefined
2.1378 +  fun  mk_sos mons =
2.1379 +   let fun mk_sq (c,m) =
2.1380 +    (c,fold_rev (fn k=> fn a => Monomialfunc.updatep iszero (nth mons (k - 1), int_element m k) a)
2.1381 +                 (1 upto length mons) Monomialfunc.undefined)
2.1382 +   in map mk_sq
2.1383 +   end
2.1384 +  val sqs = map2 mk_sos sqmonlist ratdias
2.1385 +  val cfs = map poly_of_epoly ids
2.1386 +  val msq = filter (fn (a,b) => not (null b)) (map2 pair monoid sqs)
2.1387 +  fun eval_sq sqs = fold_rev (fn (c,q) => poly_add (poly_cmul c (poly_mul q q))) sqs poly_0
2.1388 +  val sanity =
2.1389 +    fold_rev (fn ((p,c),s) => poly_add (poly_mul p (eval_sq s))) msq
2.1390 +           (fold_rev2 (fn p => fn q => poly_add (poly_mul p q)) cfs eqs
2.1391 +                    (poly_neg pol))
2.1392 +
2.1393 +in if not(Monomialfunc.is_undefined sanity) then raise Sanity else
2.1394 +  (cfs,map (fn (a,b) => (snd a,b)) msq)
2.1395 + end
2.1396 +
2.1397 +
2.1398 +end;
2.1399 +
2.1400 +(* Iterative deepening.                                                      *)
2.1401 +
2.1402 +fun deepen f n =
2.1403 +  (writeln ("Searching with depth limit " ^ string_of_int n) ; (f n handle ERROR s => (writeln ("failed with message: " ^ s) ; deepen f (n+1))))
2.1404 +
2.1405 +(* The ordering so we can create canonical HOL polynomials.                  *)
2.1406 +
2.1407 +fun dest_monomial mon = sort (increasing fst cterm_ord) (Ctermfunc.graph mon);
2.1408 +
2.1409 +fun monomial_order (m1,m2) =
2.1410 + if Ctermfunc.is_undefined m2 then LESS
2.1411 + else if Ctermfunc.is_undefined m1 then GREATER
2.1412 + else
2.1413 +  let val mon1 = dest_monomial m1
2.1414 +      val mon2 = dest_monomial m2
2.1415 +      val deg1 = fold (curry op + o snd) mon1 0
2.1416 +      val deg2 = fold (curry op + o snd) mon2 0
2.1417 +  in if deg1 < deg2 then GREATER else if deg1 > deg2 then LESS
2.1418 +     else list_ord (prod_ord cterm_ord int_ord) (mon1,mon2)
2.1419 +  end;
2.1420 +
2.1421 +fun dest_poly p =
2.1422 +  map (fn (m,c) => (c,dest_monomial m))
2.1423 +      (sort (prod_ord monomial_order (K EQUAL)) (Monomialfunc.graph p));
2.1424 +
2.1425 +(* Map back polynomials and their composites to HOL.                         *)
2.1426 +
2.1427 +local
2.1428 + open Thm Numeral RealArith
2.1429 +in
2.1430 +
2.1431 +fun cterm_of_varpow x k = if k = 1 then x else capply (capply @{cterm "op ^ :: real => _"} x)
2.1432 +  (mk_cnumber @{ctyp nat} k)
2.1433 +
2.1434 +fun cterm_of_monomial m =
2.1435 + if Ctermfunc.is_undefined m then @{cterm "1::real"}
2.1436 + else
2.1437 +  let
2.1438 +   val m' = dest_monomial m
2.1439 +   val vps = fold_rev (fn (x,k) => cons (cterm_of_varpow x k)) m' []
2.1440 +  in fold1 (fn s => fn t => capply (capply @{cterm "op * :: real => _"} s) t) vps
2.1441 +  end
2.1442 +
2.1443 +fun cterm_of_cmonomial (m,c) = if Ctermfunc.is_undefined m then cterm_of_rat c
2.1444 +    else if c = Rat.one then cterm_of_monomial m
2.1445 +    else capply (capply @{cterm "op *::real => _"} (cterm_of_rat c)) (cterm_of_monomial m);
2.1446 +
2.1447 +fun cterm_of_poly p =
2.1448 + if Monomialfunc.is_undefined p then @{cterm "0::real"}
2.1449 + else
2.1450 +  let
2.1451 +   val cms = map cterm_of_cmonomial
2.1452 +     (sort (prod_ord monomial_order (K EQUAL)) (Monomialfunc.graph p))
2.1453 +  in fold1 (fn t1 => fn t2 => capply(capply @{cterm "op + :: real => _"} t1) t2) cms
2.1454 +  end;
2.1455 +
2.1456 +fun cterm_of_sqterm (c,p) = Product(Rational_lt c,Square(cterm_of_poly p));
2.1457 +
2.1458 +fun cterm_of_sos (pr,sqs) = if null sqs then pr
2.1459 +  else Product(pr,fold1 (fn a => fn b => Sum(a,b)) (map cterm_of_sqterm sqs));
2.1460 +
2.1461 +end
2.1462 +
2.1463 +(* Interface to HOL.                                                         *)
2.1464 +local
2.1465 +  open Thm Conv RealArith
2.1466 +  val concl = dest_arg o cprop_of
2.1467 +  fun simple_cterm_ord t u = TermOrd.fast_term_ord (term_of t, term_of u) = LESS
2.1468 +in
2.1469 +  (* FIXME: Replace tryfind by get_first !! *)
2.1470 +fun real_nonlinear_prover ctxt =
2.1471 + let
2.1472 +  val {add,mul,neg,pow,sub,main} =  Normalizer.semiring_normalizers_ord_wrapper ctxt
2.1473 +      (valOf (NormalizerData.match ctxt @{cterm "(0::real) + 1"}))
2.1474 +     simple_cterm_ord
2.1475 +  val (real_poly_add_conv,real_poly_mul_conv,real_poly_neg_conv,
2.1476 +       real_poly_pow_conv,real_poly_sub_conv,real_poly_conv) = (add,mul,neg,pow,sub,main)
2.1477 +  fun mainf  translator (eqs,les,lts) =
2.1478 +  let
2.1479 +   val eq0 = map (poly_of_term o dest_arg1 o concl) eqs
2.1480 +   val le0 = map (poly_of_term o dest_arg o concl) les
2.1481 +   val lt0 = map (poly_of_term o dest_arg o concl) lts
2.1482 +   val eqp0 = map (fn (t,i) => (t,Axiom_eq i)) (eq0 ~~ (0 upto (length eq0 - 1)))
2.1483 +   val lep0 = map (fn (t,i) => (t,Axiom_le i)) (le0 ~~ (0 upto (length le0 - 1)))
2.1484 +   val ltp0 = map (fn (t,i) => (t,Axiom_lt i)) (lt0 ~~ (0 upto (length lt0 - 1)))
2.1485 +   val (keq,eq) = List.partition (fn (p,_) => multidegree p = 0) eqp0
2.1486 +   val (klep,lep) = List.partition (fn (p,_) => multidegree p = 0) lep0
2.1487 +   val (kltp,ltp) = List.partition (fn (p,_) => multidegree p = 0) ltp0
2.1488 +   fun trivial_axiom (p,ax) =
2.1489 +    case ax of
2.1490 +       Axiom_eq n => if eval Ctermfunc.undefined p <>/ Rat.zero then nth eqs n
2.1491 +                     else error "trivial_axiom: Not a trivial axiom"
2.1492 +     | Axiom_le n => if eval Ctermfunc.undefined p </ Rat.zero then nth les n
2.1493 +                     else error "trivial_axiom: Not a trivial axiom"
2.1494 +     | Axiom_lt n => if eval Ctermfunc.undefined p <=/ Rat.zero then nth lts n
2.1495 +                     else error "trivial_axiom: Not a trivial axiom"
2.1496 +     | _ => error "trivial_axiom: Not a trivial axiom"
2.1497 +   in
2.1498 +  ((let val th = tryfind trivial_axiom (keq @ klep @ kltp)
2.1499 +   in fconv_rule (arg_conv (arg1_conv real_poly_conv) then_conv field_comp_conv) th end)
2.1500 +   handle ERROR _ => (
2.1501 +    let
2.1502 +     val pol = fold_rev poly_mul (map fst ltp) (poly_const Rat.one)
2.1503 +     val leq = lep @ ltp
2.1504 +     fun tryall d =
2.1505 +      let val e = multidegree pol
2.1506 +          val k = if e = 0 then 0 else d div e
2.1507 +          val eq' = map fst eq
2.1508 +      in tryfind (fn i => (d,i,real_positivnullstellensatz_general false d eq' leq
2.1509 +                            (poly_neg(poly_pow pol i))))
2.1510 +              (0 upto k)
2.1511 +      end
2.1512 +    val (d,i,(cert_ideal,cert_cone)) = deepen tryall 0
2.1513 +    val proofs_ideal =
2.1514 +      map2 (fn q => fn (p,ax) => Eqmul(cterm_of_poly q,ax)) cert_ideal eq
2.1515 +    val proofs_cone = map cterm_of_sos cert_cone
2.1516 +    val proof_ne = if null ltp then Rational_lt Rat.one else
2.1517 +      let val p = fold1 (fn s => fn t => Product(s,t)) (map snd ltp)
2.1518 +      in  funpow i (fn q => Product(p,q)) (Rational_lt Rat.one)
2.1519 +      end
2.1520 +    val proof = fold1 (fn s => fn t => Sum(s,t))
2.1521 +                           (proof_ne :: proofs_ideal @ proofs_cone)
2.1522 +    in writeln "Translating proof certificate to HOL";
2.1523 +       translator (eqs,les,lts) proof
2.1524 +    end))
2.1525 +   end
2.1526 + in mainf end
2.1527 +end
2.1528 +
2.1529 +fun C f x y = f y x;
2.1530 +  (* FIXME : This is very bad!!!*)
2.1531 +fun subst_conv eqs t =
2.1532 + let
2.1533 +  val t' = fold (Thm.cabs o Thm.lhs_of) eqs t
2.1534 + in Conv.fconv_rule (Thm.beta_conversion true) (fold (C combination) eqs (reflexive t'))
2.1535 + end
2.1536 +
2.1537 +(* A wrapper that tries to substitute away variables first.                  *)
2.1538 +
2.1539 +local
2.1540 + open Thm Conv RealArith
2.1541 +  fun simple_cterm_ord t u = TermOrd.fast_term_ord (term_of t, term_of u) = LESS
2.1542 + val concl = dest_arg o cprop_of
2.1543 + val shuffle1 =
2.1544 +   fconv_rule (rewr_conv @{lemma "(a + x == y) == (x == y - (a::real))" by (atomize (full)) (simp add: ring_simps) })
2.1545 + val shuffle2 =
2.1546 +    fconv_rule (rewr_conv @{lemma "(x + a == y) ==  (x == y - (a::real))" by (atomize (full)) (simp add: ring_simps)})
2.1547 + fun substitutable_monomial fvs tm = case term_of tm of
2.1548 +    Free(_,@{typ real}) => if not (member (op aconvc) fvs tm) then (Rat.one,tm)
2.1549 +                           else error "substitutable_monomial"
2.1550 +  | @{term "op * :: real => _"}\$c\$(t as Free _ ) =>
2.1551 +     if is_ratconst (dest_arg1 tm) andalso not (member (op aconvc) fvs (dest_arg tm))
2.1552 +         then (dest_ratconst (dest_arg1 tm),dest_arg tm) else error "substitutable_monomial"
2.1553 +  | @{term "op + :: real => _"}\$s\$t =>
2.1554 +       (substitutable_monomial (add_cterm_frees (dest_arg tm) fvs) (dest_arg1 tm)
2.1555 +        handle ERROR _ => substitutable_monomial (add_cterm_frees (dest_arg1 tm) fvs) (dest_arg tm))
2.1556 +  | _ => error "substitutable_monomial"
2.1557 +
2.1558 +  fun isolate_variable v th =
2.1559 +   let val w = dest_arg1 (cprop_of th)
2.1560 +   in if v aconvc w then th
2.1561 +      else case term_of w of
2.1562 +           @{term "op + :: real => _"}\$s\$t =>
2.1563 +              if dest_arg1 w aconvc v then shuffle2 th
2.1564 +              else isolate_variable v (shuffle1 th)
2.1565 +          | _ => error "isolate variable : This should not happen?"
2.1566 +   end
2.1567 +in
2.1568 +
2.1569 +fun real_nonlinear_subst_prover ctxt =
2.1570 + let
2.1571 +  val {add,mul,neg,pow,sub,main} =  Normalizer.semiring_normalizers_ord_wrapper ctxt
2.1572 +      (valOf (NormalizerData.match ctxt @{cterm "(0::real) + 1"}))
2.1573 +     simple_cterm_ord
2.1574 +
2.1575 +  val (real_poly_add_conv,real_poly_mul_conv,real_poly_neg_conv,
2.1576 +       real_poly_pow_conv,real_poly_sub_conv,real_poly_conv) = (add,mul,neg,pow,sub,main)
2.1577 +
2.1578 +  fun make_substitution th =
2.1579 +   let
2.1580 +    val (c,v) = substitutable_monomial [] (dest_arg1(concl th))
2.1581 +    val th1 = Drule.arg_cong_rule (capply @{cterm "op * :: real => _"} (cterm_of_rat (Rat.inv c))) (mk_meta_eq th)
2.1582 +    val th2 = fconv_rule (binop_conv real_poly_mul_conv) th1
2.1583 +   in fconv_rule (arg_conv real_poly_conv) (isolate_variable v th2)
2.1584 +   end
2.1585 +   fun oprconv cv ct =
2.1586 +    let val g = Thm.dest_fun2 ct
2.1587 +    in if g aconvc @{cterm "op <= :: real => _"}
2.1588 +         orelse g aconvc @{cterm "op < :: real => _"}
2.1589 +       then arg_conv cv ct else arg1_conv cv ct
2.1590 +    end
2.1591 +  fun mainf translator =
2.1592 +   let
2.1593 +    fun substfirst(eqs,les,lts) =
2.1594 +      ((let
2.1595 +           val eth = tryfind make_substitution eqs
2.1596 +           val modify = fconv_rule (arg_conv (oprconv(subst_conv [eth] then_conv real_poly_conv)))
2.1597 +       in  substfirst
2.1598 +             (filter_out (fn t => (Thm.dest_arg1 o Thm.dest_arg o cprop_of) t
2.1599 +                                   aconvc @{cterm "0::real"}) (map modify eqs),
2.1600 +                                   map modify les,map modify lts)
2.1601 +       end)
2.1602 +       handle ERROR  _ => real_nonlinear_prover ctxt translator (rev eqs, rev les, rev lts))
2.1603 +    in substfirst
2.1604 +   end
2.1605 +
2.1606 +
2.1607 + in mainf
2.1608 + end
2.1609 +
2.1610 +(* Overall function. *)
2.1611 +
2.1612 +fun real_sos ctxt t = gen_prover_real_arith ctxt (real_nonlinear_subst_prover ctxt) t;
2.1613 +end;
2.1614 +
2.1615 +end;
2.1616 \ No newline at end of file
```