A decision method for universal multivariate real arithmetic with add
authorchaieb
Tue May 12 17:32:49 2009 +0100 (2009-05-12)
changeset 311192532bb2d65c7
parent 31118 541d43bee678
child 31120 fc654c95c29e
A decision method for universal multivariate real arithmetic with add
ition, multiplication and ordering using semidefinite programming
src/HOL/Library/Sum_Of_Squares.thy
src/HOL/Library/sum_of_squares.ML
     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