src/HOL/Library/sum_of_squares.ML
author huffman
Thu Jun 11 09:03:24 2009 -0700 (2009-06-11)
changeset 31563 ded2364d14d4
parent 31512 27118561c2e0
child 32150 4ed2865f3a56
permissions -rw-r--r--
cleaned up some proofs
     1 structure Sos = 
     2 struct
     3 
     4 val rat_0 = Rat.zero;
     5 val rat_1 = Rat.one;
     6 val rat_2 = Rat.two;
     7 val rat_10 = Rat.rat_of_int 10;
     8 val rat_1_2 = rat_1 // rat_2;
     9 val max = curry IntInf.max;
    10 val min = curry IntInf.min;
    11 
    12 val denominator_rat = Rat.quotient_of_rat #> snd #> Rat.rat_of_int;
    13 val numerator_rat = Rat.quotient_of_rat #> fst #> Rat.rat_of_int;
    14 fun int_of_rat a = 
    15     case Rat.quotient_of_rat a of (i,1) => i | _ => error "int_of_rat: not an int";
    16 fun lcm_rat x y = Rat.rat_of_int (Integer.lcm (int_of_rat x) (int_of_rat y));
    17 
    18 fun rat_pow r i = 
    19  let fun pow r i = 
    20    if i = 0 then rat_1 else 
    21    let val d = pow r (i div 2)
    22    in d */ d */ (if i mod 2 = 0 then rat_1 else r)
    23    end
    24  in if i < 0 then pow (Rat.inv r) (~ i) else pow r i end;
    25 
    26 fun round_rat r = 
    27  let val (a,b) = Rat.quotient_of_rat (Rat.abs r)
    28      val d = a div b
    29      val s = if r </ rat_0 then (Rat.neg o Rat.rat_of_int) else Rat.rat_of_int
    30      val x2 = 2 * (a - (b * d))
    31  in s (if x2 >= b then d + 1 else d) end
    32 
    33 val abs_rat = Rat.abs;
    34 val pow2 = rat_pow rat_2;
    35 val pow10 = rat_pow rat_10;
    36 
    37 val debugging = ref false;
    38 
    39 exception Sanity;
    40 
    41 exception Unsolvable;
    42 
    43 (* Turn a rational into a decimal string with d sig digits.                  *)
    44 
    45 local
    46 fun normalize y =
    47   if abs_rat y </ (rat_1 // rat_10) then normalize (rat_10 */ y) - 1
    48   else if abs_rat y >=/ rat_1 then normalize (y // rat_10) + 1
    49   else 0 
    50  in
    51 fun decimalize d x =
    52   if x =/ rat_0 then "0.0" else
    53   let 
    54    val y = Rat.abs x
    55    val e = normalize y
    56    val z = pow10(~ e) */ y +/ rat_1
    57    val k = int_of_rat (round_rat(pow10 d */ z))
    58   in (if x </ rat_0 then "-0." else "0.") ^
    59      implode(tl(explode(string_of_int k))) ^
    60      (if e = 0 then "" else "e"^string_of_int e)
    61   end
    62 end;
    63 
    64 (* Iterations over numbers, and lists indexed by numbers.                    *)
    65 
    66 fun itern k l f a =
    67   case l of
    68     [] => a
    69   | h::t => itern (k + 1) t f (f h k a);
    70 
    71 fun iter (m,n) f a =
    72   if n < m then a
    73   else iter (m+1,n) f (f m a);
    74 
    75 (* The main types.                                                           *)
    76 
    77 fun strict_ord ord (x,y) = case ord (x,y) of LESS => LESS | _ => GREATER
    78 
    79 structure Intpairfunc = FuncFun(type key = int*int val ord = prod_ord int_ord int_ord);
    80 
    81 type vector = int* Rat.rat Intfunc.T;
    82 
    83 type matrix = (int*int)*(Rat.rat Intpairfunc.T);
    84 
    85 type monomial = int Ctermfunc.T;
    86 
    87 val cterm_ord = (fn (s,t) => TermOrd.fast_term_ord(term_of s, term_of t))
    88  fun monomial_ord (m1,m2) = list_ord (prod_ord cterm_ord int_ord) (Ctermfunc.graph m1, Ctermfunc.graph m2)
    89 structure Monomialfunc = FuncFun(type key = monomial val ord = monomial_ord)
    90 
    91 type poly = Rat.rat Monomialfunc.T;
    92 
    93  fun iszero (k,r) = r =/ rat_0;
    94 
    95 fun fold_rev2 f l1 l2 b =
    96   case (l1,l2) of
    97     ([],[]) => b
    98   | (h1::t1,h2::t2) => f h1 h2 (fold_rev2 f t1 t2 b)
    99   | _ => error "fold_rev2";
   100  
   101 (* Vectors. Conventionally indexed 1..n.                                     *)
   102 
   103 fun vector_0 n = (n,Intfunc.undefined):vector;
   104 
   105 fun dim (v:vector) = fst v;
   106 
   107 fun vector_const c n =
   108   if c =/ rat_0 then vector_0 n
   109   else (n,fold_rev (fn k => Intfunc.update (k,c)) (1 upto n) Intfunc.undefined) :vector;
   110 
   111 val vector_1 = vector_const rat_1;
   112 
   113 fun vector_cmul c (v:vector) =
   114  let val n = dim v 
   115  in if c =/ rat_0 then vector_0 n
   116     else (n,Intfunc.mapf (fn x => c */ x) (snd v))
   117  end;
   118 
   119 fun vector_neg (v:vector) = (fst v,Intfunc.mapf Rat.neg (snd v)) :vector;
   120 
   121 fun vector_add (v1:vector) (v2:vector) =
   122  let val m = dim v1  
   123      val n = dim v2 
   124  in if m <> n then error "vector_add: incompatible dimensions"
   125     else (n,Intfunc.combine (curry op +/) (fn x => x =/ rat_0) (snd v1) (snd v2)) :vector 
   126  end;
   127 
   128 fun vector_sub v1 v2 = vector_add v1 (vector_neg v2);
   129 
   130 fun vector_dot (v1:vector) (v2:vector) =
   131  let val m = dim v1 
   132      val n = dim v2 
   133  in if m <> n then error "vector_dot: incompatible dimensions" 
   134     else Intfunc.fold (fn (i,x) => fn a => x +/ a) 
   135         (Intfunc.combine (curry op */) (fn x => x =/ rat_0) (snd v1) (snd v2)) rat_0
   136  end;
   137 
   138 fun vector_of_list l =
   139  let val n = length l 
   140  in (n,fold_rev2 (curry Intfunc.update) (1 upto n) l Intfunc.undefined) :vector
   141  end;
   142 
   143 (* Matrices; again rows and columns indexed from 1.                          *)
   144 
   145 fun matrix_0 (m,n) = ((m,n),Intpairfunc.undefined):matrix;
   146 
   147 fun dimensions (m:matrix) = fst m;
   148 
   149 fun matrix_const c (mn as (m,n)) =
   150   if m <> n then error "matrix_const: needs to be square"
   151   else if c =/ rat_0 then matrix_0 mn
   152   else (mn,fold_rev (fn k => Intpairfunc.update ((k,k), c)) (1 upto n) Intpairfunc.undefined) :matrix;;
   153 
   154 val matrix_1 = matrix_const rat_1;
   155 
   156 fun matrix_cmul c (m:matrix) =
   157  let val (i,j) = dimensions m 
   158  in if c =/ rat_0 then matrix_0 (i,j)
   159     else ((i,j),Intpairfunc.mapf (fn x => c */ x) (snd m))
   160  end;
   161 
   162 fun matrix_neg (m:matrix) = 
   163   (dimensions m, Intpairfunc.mapf Rat.neg (snd m)) :matrix;
   164 
   165 fun matrix_add (m1:matrix) (m2:matrix) =
   166  let val d1 = dimensions m1 
   167      val d2 = dimensions m2 
   168  in if d1 <> d2 
   169      then error "matrix_add: incompatible dimensions"
   170     else (d1,Intpairfunc.combine (curry op +/) (fn x => x =/ rat_0) (snd m1) (snd m2)) :matrix
   171  end;;
   172 
   173 fun matrix_sub m1 m2 = matrix_add m1 (matrix_neg m2);
   174 
   175 fun row k (m:matrix) =
   176  let val (i,j) = dimensions m 
   177  in (j,
   178    Intpairfunc.fold (fn ((i,j), c) => fn a => if i = k then Intfunc.update (j,c) a else a) (snd m) Intfunc.undefined ) : vector
   179  end;
   180 
   181 fun column k (m:matrix) =
   182   let val (i,j) = dimensions m 
   183   in (i,
   184    Intpairfunc.fold (fn ((i,j), c) => fn a => if j = k then Intfunc.update (i,c) a else a) (snd m)  Intfunc.undefined)
   185    : vector
   186  end;
   187 
   188 fun transp (m:matrix) =
   189   let val (i,j) = dimensions m 
   190   in
   191   ((j,i),Intpairfunc.fold (fn ((i,j), c) => fn a => Intpairfunc.update ((j,i), c) a) (snd m) Intpairfunc.undefined) :matrix
   192  end;
   193 
   194 fun diagonal (v:vector) =
   195  let val n = dim v 
   196  in ((n,n),Intfunc.fold (fn (i, c) => fn a => Intpairfunc.update ((i,i), c) a) (snd v) Intpairfunc.undefined) : matrix
   197  end;
   198 
   199 fun matrix_of_list l =
   200  let val m = length l 
   201  in if m = 0 then matrix_0 (0,0) else
   202    let val n = length (hd l) 
   203    in ((m,n),itern 1 l (fn v => fn i => itern 1 v (fn c => fn j => Intpairfunc.update ((i,j), c))) Intpairfunc.undefined)
   204    end
   205  end;
   206 
   207 (* Monomials.                                                                *)
   208 
   209 fun monomial_eval assig (m:monomial) =
   210   Ctermfunc.fold (fn (x, k) => fn a => a */ rat_pow (Ctermfunc.apply assig x) k)
   211         m rat_1;
   212 val monomial_1 = (Ctermfunc.undefined:monomial);
   213 
   214 fun monomial_var x = Ctermfunc.onefunc (x, 1) :monomial;
   215 
   216 val (monomial_mul:monomial->monomial->monomial) =
   217   Ctermfunc.combine (curry op +) (K false);
   218 
   219 fun monomial_pow (m:monomial) k =
   220   if k = 0 then monomial_1
   221   else Ctermfunc.mapf (fn x => k * x) m;
   222 
   223 fun monomial_divides (m1:monomial) (m2:monomial) =
   224   Ctermfunc.fold (fn (x, k) => fn a => Ctermfunc.tryapplyd m2 x 0 >= k andalso a) m1 true;;
   225 
   226 fun monomial_div (m1:monomial) (m2:monomial) =
   227  let val m = Ctermfunc.combine (curry op +) 
   228    (fn x => x = 0) m1 (Ctermfunc.mapf (fn x => ~ x) m2) 
   229  in if Ctermfunc.fold (fn (x, k) => fn a => k >= 0 andalso a) m true then m
   230   else error "monomial_div: non-divisible"
   231  end;
   232 
   233 fun monomial_degree x (m:monomial) = 
   234   Ctermfunc.tryapplyd m x 0;;
   235 
   236 fun monomial_lcm (m1:monomial) (m2:monomial) =
   237   fold_rev (fn x => Ctermfunc.update (x, max (monomial_degree x m1) (monomial_degree x m2)))
   238           (gen_union (is_equal o  cterm_ord) (Ctermfunc.dom m1, Ctermfunc.dom m2)) (Ctermfunc.undefined :monomial);
   239 
   240 fun monomial_multidegree (m:monomial) = 
   241  Ctermfunc.fold (fn (x, k) => fn a => k + a) m 0;;
   242 
   243 fun monomial_variables m = Ctermfunc.dom m;;
   244 
   245 (* Polynomials.                                                              *)
   246 
   247 fun eval assig (p:poly) =
   248   Monomialfunc.fold (fn (m, c) => fn a => a +/ c */ monomial_eval assig m) p rat_0;
   249 
   250 val poly_0 = (Monomialfunc.undefined:poly);
   251 
   252 fun poly_isconst (p:poly) = 
   253   Monomialfunc.fold (fn (m, c) => fn a => Ctermfunc.is_undefined m andalso a) p true;
   254 
   255 fun poly_var x = Monomialfunc.onefunc (monomial_var x,rat_1) :poly;
   256 
   257 fun poly_const c =
   258   if c =/ rat_0 then poly_0 else Monomialfunc.onefunc(monomial_1, c);
   259 
   260 fun poly_cmul c (p:poly) =
   261   if c =/ rat_0 then poly_0
   262   else Monomialfunc.mapf (fn x => c */ x) p;
   263 
   264 fun poly_neg (p:poly) = (Monomialfunc.mapf Rat.neg p :poly);;
   265 
   266 fun poly_add (p1:poly) (p2:poly) =
   267   (Monomialfunc.combine (curry op +/) (fn x => x =/ rat_0) p1 p2 :poly);
   268 
   269 fun poly_sub p1 p2 = poly_add p1 (poly_neg p2);
   270 
   271 fun poly_cmmul (c,m) (p:poly) =
   272  if c =/ rat_0 then poly_0
   273  else if Ctermfunc.is_undefined m 
   274       then Monomialfunc.mapf (fn d => c */ d) p
   275       else Monomialfunc.fold (fn (m', d) => fn a => (Monomialfunc.update (monomial_mul m m', c */ d) a)) p poly_0;
   276 
   277 fun poly_mul (p1:poly) (p2:poly) =
   278   Monomialfunc.fold (fn (m, c) => fn a => poly_add (poly_cmmul (c,m) p2) a) p1 poly_0;
   279 
   280 fun poly_div (p1:poly) (p2:poly) =
   281  if not(poly_isconst p2) 
   282  then error "poly_div: non-constant" else
   283  let val c = eval Ctermfunc.undefined p2 
   284  in if c =/ rat_0 then error "poly_div: division by zero"
   285     else poly_cmul (Rat.inv c) p1
   286  end;
   287 
   288 fun poly_square p = poly_mul p p;
   289 
   290 fun poly_pow p k =
   291  if k = 0 then poly_const rat_1
   292  else if k = 1 then p
   293  else let val q = poly_square(poly_pow p (k div 2)) in
   294       if k mod 2 = 1 then poly_mul p q else q end;
   295 
   296 fun poly_exp p1 p2 =
   297   if not(poly_isconst p2) 
   298   then error "poly_exp: not a constant" 
   299   else poly_pow p1 (int_of_rat (eval Ctermfunc.undefined p2));
   300 
   301 fun degree x (p:poly) = 
   302  Monomialfunc.fold (fn (m,c) => fn a => max (monomial_degree x m) a) p 0;
   303 
   304 fun multidegree (p:poly) =
   305   Monomialfunc.fold (fn (m, c) => fn a => max (monomial_multidegree m) a) p 0;
   306 
   307 fun poly_variables (p:poly) =
   308   sort cterm_ord (Monomialfunc.fold_rev (fn (m, c) => curry (gen_union (is_equal o  cterm_ord)) (monomial_variables m)) p []);;
   309 
   310 (* Order monomials for human presentation.                                   *)
   311 
   312 fun cterm_ord (t,t') = TermOrd.fast_term_ord (term_of t, term_of t');
   313 
   314 val humanorder_varpow = prod_ord cterm_ord (rev_order o int_ord);
   315 
   316 local
   317  fun ord (l1,l2) = case (l1,l2) of
   318   (_,[]) => LESS 
   319  | ([],_) => GREATER
   320  | (h1::t1,h2::t2) => 
   321    (case humanorder_varpow (h1, h2) of 
   322      LESS => LESS
   323    | EQUAL => ord (t1,t2)
   324    | GREATER => GREATER)
   325 in fun humanorder_monomial m1 m2 = 
   326  ord (sort humanorder_varpow (Ctermfunc.graph m1),
   327   sort humanorder_varpow (Ctermfunc.graph m2))
   328 end;
   329 
   330 fun fold1 f l =  case l of
   331    []     => error "fold1"
   332  | [x]    => x
   333  | (h::t) => f h (fold1 f t);
   334 
   335 (* Conversions to strings.                                                   *)
   336 
   337 fun string_of_vector min_size max_size (v:vector) =
   338  let val n_raw = dim v 
   339  in if n_raw = 0 then "[]" else
   340   let 
   341    val n = max min_size (min n_raw max_size) 
   342    val xs = map (Rat.string_of_rat o (fn i => Intfunc.tryapplyd (snd v) i rat_0)) (1 upto n) 
   343   in "[" ^ fold1 (fn s => fn t => s ^ ", " ^ t) xs ^
   344   (if n_raw > max_size then ", ...]" else "]")
   345   end
   346  end;
   347 
   348 fun string_of_matrix max_size (m:matrix) =
   349  let 
   350   val (i_raw,j_raw) = dimensions m
   351   val i = min max_size i_raw 
   352   val j = min max_size j_raw
   353   val rstr = map (fn k => string_of_vector j j (row k m)) (1 upto i) 
   354  in "["^ fold1 (fn s => fn t => s^";\n "^t) rstr ^
   355   (if j > max_size then "\n ...]" else "]")
   356  end;
   357 
   358 fun string_of_term t = 
   359  case t of
   360    a$b => "("^(string_of_term a)^" "^(string_of_term b)^")"
   361  | Abs x => 
   362     let val (xn, b) = Term.dest_abs x
   363     in "(\\"^xn^"."^(string_of_term b)^")"
   364     end
   365  | Const(s,_) => s
   366  | Free (s,_) => s
   367  | Var((s,_),_) => s
   368  | _ => error "string_of_term";
   369 
   370 val string_of_cterm = string_of_term o term_of;
   371 
   372 fun string_of_varpow x k =
   373   if k = 1 then string_of_cterm x 
   374   else string_of_cterm x^"^"^string_of_int k;
   375 
   376 fun string_of_monomial m =
   377  if Ctermfunc.is_undefined m then "1" else
   378  let val vps = fold_rev (fn (x,k) => fn a => string_of_varpow x k :: a)
   379   (sort humanorder_varpow (Ctermfunc.graph m)) [] 
   380  in fold1 (fn s => fn t => s^"*"^t) vps
   381  end;
   382 
   383 fun string_of_cmonomial (c,m) =
   384  if Ctermfunc.is_undefined m then Rat.string_of_rat c
   385  else if c =/ rat_1 then string_of_monomial m
   386  else Rat.string_of_rat c ^ "*" ^ string_of_monomial m;;
   387 
   388 fun string_of_poly (p:poly) =
   389  if Monomialfunc.is_undefined p then "<<0>>" else
   390  let 
   391   val cms = sort (fn ((m1,_),(m2,_)) => humanorder_monomial m1  m2) (Monomialfunc.graph p)
   392   val s = fold (fn (m,c) => fn a =>
   393              if c </ rat_0 then a ^ " - " ^ string_of_cmonomial(Rat.neg c,m)
   394              else a ^ " + " ^ string_of_cmonomial(c,m))
   395           cms ""
   396   val s1 = String.substring (s, 0, 3)
   397   val s2 = String.substring (s, 3, String.size s - 3) 
   398  in "<<" ^(if s1 = " + " then s2 else "-"^s2)^">>"
   399  end;
   400 
   401 (* Conversion from HOL term.                                                 *)
   402 
   403 local
   404  val neg_tm = @{cterm "uminus :: real => _"}
   405  val add_tm = @{cterm "op + :: real => _"}
   406  val sub_tm = @{cterm "op - :: real => _"}
   407  val mul_tm = @{cterm "op * :: real => _"}
   408  val inv_tm = @{cterm "inverse :: real => _"}
   409  val div_tm = @{cterm "op / :: real => _"}
   410  val pow_tm = @{cterm "op ^ :: real => _"}
   411  val zero_tm = @{cterm "0:: real"}
   412  val is_numeral = can (HOLogic.dest_number o term_of)
   413  fun is_comb t = case t of _$_ => true | _ => false
   414  fun poly_of_term tm =
   415   if tm aconvc zero_tm then poly_0
   416   else if RealArith.is_ratconst tm 
   417        then poly_const(RealArith.dest_ratconst tm)
   418   else 
   419   (let val (lop,r) = Thm.dest_comb tm
   420    in if lop aconvc neg_tm then poly_neg(poly_of_term r)
   421       else if lop aconvc inv_tm then
   422        let val p = poly_of_term r 
   423        in if poly_isconst p 
   424           then poly_const(Rat.inv (eval Ctermfunc.undefined p))
   425           else error "poly_of_term: inverse of non-constant polyomial"
   426        end
   427    else (let val (opr,l) = Thm.dest_comb lop
   428          in 
   429           if opr aconvc pow_tm andalso is_numeral r 
   430           then poly_pow (poly_of_term l) ((snd o HOLogic.dest_number o term_of) r)
   431           else if opr aconvc add_tm 
   432            then poly_add (poly_of_term l) (poly_of_term r)
   433           else if opr aconvc sub_tm 
   434            then poly_sub (poly_of_term l) (poly_of_term r)
   435           else if opr aconvc mul_tm 
   436            then poly_mul (poly_of_term l) (poly_of_term r)
   437           else if opr aconvc div_tm 
   438            then let 
   439                   val p = poly_of_term l 
   440                   val q = poly_of_term r 
   441                 in if poly_isconst q then poly_cmul (Rat.inv (eval Ctermfunc.undefined q)) p
   442                    else error "poly_of_term: division by non-constant polynomial"
   443                 end
   444           else poly_var tm
   445  
   446          end
   447          handle CTERM ("dest_comb",_) => poly_var tm)
   448    end
   449    handle CTERM ("dest_comb",_) => poly_var tm)
   450 in
   451 val poly_of_term = fn tm =>
   452  if type_of (term_of tm) = @{typ real} then poly_of_term tm
   453  else error "poly_of_term: term does not have real type"
   454 end;
   455 
   456 (* String of vector (just a list of space-separated numbers).                *)
   457 
   458 fun sdpa_of_vector (v:vector) =
   459  let 
   460   val n = dim v
   461   val strs = map (decimalize 20 o (fn i => Intfunc.tryapplyd (snd v) i rat_0)) (1 upto n) 
   462  in fold1 (fn x => fn y => x ^ " " ^ y) strs ^ "\n"
   463  end;
   464 
   465 fun increasing f ord (x,y) = ord (f x, f y);
   466 fun triple_int_ord ((a,b,c),(a',b',c')) = 
   467  prod_ord int_ord (prod_ord int_ord int_ord) 
   468     ((a,(b,c)),(a',(b',c')));
   469 structure Inttriplefunc = FuncFun(type key = int*int*int val ord = triple_int_ord);
   470 
   471 (* String for block diagonal matrix numbered k.                              *)
   472 
   473 fun sdpa_of_blockdiagonal k m =
   474  let 
   475   val pfx = string_of_int k ^" "
   476   val ents =
   477     Inttriplefunc.fold (fn ((b,i,j), c) => fn a => if i > j then a else ((b,i,j),c)::a) m []
   478   val entss = sort (increasing fst triple_int_ord ) ents
   479  in  fold_rev (fn ((b,i,j),c) => fn a =>
   480      pfx ^ string_of_int b ^ " " ^ string_of_int i ^ " " ^ string_of_int j ^
   481      " " ^ decimalize 20 c ^ "\n" ^ a) entss ""
   482  end;
   483 
   484 (* String for a matrix numbered k, in SDPA sparse format.                    *)
   485 
   486 fun sdpa_of_matrix k (m:matrix) =
   487  let 
   488   val pfx = string_of_int k ^ " 1 "
   489   val ms = Intpairfunc.fold (fn ((i,j), c) => fn  a => if i > j then a else ((i,j),c)::a) (snd m) [] 
   490   val mss = sort (increasing fst (prod_ord int_ord int_ord)) ms 
   491  in fold_rev (fn ((i,j),c) => fn a =>
   492      pfx ^ string_of_int i ^ " " ^ string_of_int j ^
   493      " " ^ decimalize 20 c ^ "\n" ^ a) mss ""
   494  end;;
   495 
   496 (* ------------------------------------------------------------------------- *)
   497 (* String in SDPA sparse format for standard SDP problem:                    *)
   498 (*                                                                           *)
   499 (*    X = v_1 * [M_1] + ... + v_m * [M_m] - [M_0] must be PSD                *)
   500 (*    Minimize obj_1 * v_1 + ... obj_m * v_m                                 *)
   501 (* ------------------------------------------------------------------------- *)
   502 
   503 fun sdpa_of_problem comment obj mats =
   504  let 
   505   val m = length mats - 1
   506   val (n,_) = dimensions (hd mats) 
   507  in "\"" ^ comment ^ "\"\n" ^
   508   string_of_int m ^ "\n" ^
   509   "1\n" ^
   510   string_of_int n ^ "\n" ^
   511   sdpa_of_vector obj ^
   512   fold_rev2 (fn k => fn m => fn a => sdpa_of_matrix (k - 1) m ^ a) (1 upto length mats) mats ""
   513  end;
   514 
   515 fun index_char str chr pos =
   516   if pos >= String.size str then ~1
   517   else if String.sub(str,pos) = chr then pos
   518   else index_char str chr (pos + 1);
   519 fun rat_of_quotient (a,b) = if b = 0 then rat_0 else Rat.rat_of_quotient (a,b);
   520 fun rat_of_string s = 
   521  let val n = index_char s #"/" 0 in
   522   if n = ~1 then s |> IntInf.fromString |> valOf |> Rat.rat_of_int
   523   else 
   524    let val SOME numer = IntInf.fromString(String.substring(s,0,n))
   525        val SOME den = IntInf.fromString (String.substring(s,n+1,String.size s - n - 1))
   526    in rat_of_quotient(numer, den)
   527    end
   528  end;
   529 
   530 fun isspace x = x = " " ;
   531 fun isnum x = x mem_string ["0","1","2","3","4","5","6","7","8","9"]
   532 
   533 (* More parser basics.                                                       *)
   534 
   535 local
   536  open Scan
   537 in 
   538  val word = this_string
   539  fun token s =
   540   repeat ($$ " ") |-- word s --| repeat ($$ " ")
   541  val numeral = one isnum
   542  val decimalint = bulk numeral >> (rat_of_string o implode)
   543  val decimalfrac = bulk numeral
   544     >> (fn s => rat_of_string(implode s) // pow10 (length s))
   545  val decimalsig =
   546     decimalint -- option (Scan.$$ "." |-- decimalfrac)
   547     >> (fn (h,NONE) => h | (h,SOME x) => h +/ x)
   548  fun signed prs =
   549        $$ "-" |-- prs >> Rat.neg 
   550     || $$ "+" |-- prs
   551     || prs;
   552 
   553 fun emptyin def xs = if null xs then (def,xs) else Scan.fail xs
   554 
   555  val exponent = ($$ "e" || $$ "E") |-- signed decimalint;
   556 
   557  val decimal = signed decimalsig -- (emptyin rat_0|| exponent)
   558     >> (fn (h, x) => h */ pow10 (int_of_rat x));
   559 end;
   560 
   561  fun mkparser p s =
   562   let val (x,rst) = p (explode s) 
   563   in if null rst then x 
   564      else error "mkparser: unparsed input"
   565   end;;
   566 val parse_decimal = mkparser decimal;
   567 
   568 fun fix err prs = 
   569   prs || (fn x=> error err);
   570 
   571 fun listof prs sep err =
   572   prs -- Scan.bulk (sep |-- fix err prs) >> uncurry cons;
   573 
   574 (* Parse back a vector.                                                      *)
   575 
   576  val vector = 
   577     token "{" |-- listof decimal (token ",") "decimal" --| token "}"
   578                >> vector_of_list 
   579  val parse_vector = mkparser vector
   580  fun skipupto dscr prs inp =
   581    (dscr |-- prs 
   582     || Scan.one (K true) |-- skipupto dscr prs) inp 
   583  fun ignore inp = ((),[])
   584  fun sdpaoutput inp =  skipupto (word "xVec" -- token "=")
   585              (vector --| ignore) inp
   586  fun csdpoutput inp =  ((decimal -- Scan.bulk (Scan.$$ " " |-- Scan.option decimal) >> (fn (h,to) => map_filter I ((SOME h)::to))) --| ignore >> vector_of_list) inp
   587  val parse_sdpaoutput = mkparser sdpaoutput
   588  val parse_csdpoutput = mkparser csdpoutput
   589 
   590 (* Also parse the SDPA output to test success (CSDP yields a return code).   *)
   591 
   592 local
   593  val prs = 
   594   skipupto (word "phase.value" -- token "=")
   595    (Scan.option (Scan.$$ "p") -- Scan.option (Scan.$$ "d") 
   596     -- (word "OPT" || word "FEAS")) 
   597 in
   598  fun sdpa_run_succeeded s = 
   599   (prs (explode s); true) handle _ => false
   600 end;
   601 
   602 (* The default parameters. Unfortunately this goes to a fixed file.          *)
   603 
   604 val sdpa_default_parameters =
   605 "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";;
   606 
   607 (* These were suggested by Makoto Yamashita for problems where we are        *)
   608 (* right at the edge of the semidefinite cone, as sometimes happens.         *)
   609 
   610 val sdpa_alt_parameters =
   611 "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";;
   612 
   613 val sdpa_params = sdpa_alt_parameters;;
   614 
   615 (* CSDP parameters; so far I'm sticking with the defaults.                   *)
   616 
   617 val csdp_default_parameters =
   618 "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";;
   619 
   620 val csdp_params = csdp_default_parameters;;
   621 
   622 fun tmp_file pre suf =
   623  let val i = string_of_int (round (random()))
   624    val name = Path.append (Path.variable "ISABELLE_TMP") (Path.explode (pre ^ i ^ suf))
   625  in 
   626    if File.exists name then tmp_file pre suf 
   627    else name 
   628  end;
   629 
   630 (* Now call SDPA on a problem and parse back the output.                     *)
   631 
   632 fun run_sdpa dbg obj mats =
   633  let 
   634   val input_file = tmp_file "sos" ".dat-s"
   635   val output_file = tmp_file "sos" ".out"
   636   val params_file = tmp_file "param" ".sdpa" 
   637   val current_dir = File.pwd()
   638   val _ = File.write input_file 
   639                          (sdpa_of_problem "" obj mats)
   640   val _ = File.write params_file sdpa_params
   641   val _ = File.cd (Path.variable "ISABELLE_TMP")
   642   val _ = File.system_command ("sdpa "^ (Path.implode input_file) ^ " " ^ 
   643                                (Path.implode output_file) ^
   644                                (if dbg then "" else "> /dev/null"))
   645   val opr = File.read output_file 
   646  in if not(sdpa_run_succeeded opr) then error "sdpa: call failed" 
   647     else
   648       let val res = parse_sdpaoutput opr 
   649       in ((if dbg then ()
   650            else (File.rm input_file; File.rm output_file ; File.cd current_dir));
   651           res)
   652       end
   653  end;
   654 
   655 fun sdpa obj mats = run_sdpa (!debugging) obj mats;
   656 
   657 (* The same thing with CSDP.                                                 *)
   658 
   659 fun run_csdp dbg obj mats =
   660  let 
   661   val input_file = tmp_file "sos" ".dat-s"
   662   val output_file = tmp_file "sos" ".out"
   663   val params_file = tmp_file "param" ".csdp"
   664   val current_dir = File.pwd()
   665   val _ = File.write input_file (sdpa_of_problem "" obj mats)
   666   val _ = File.write params_file csdp_params
   667   val _ = File.cd (Path.variable "ISABELLE_TMP")
   668   val rv = system ("csdp "^(Path.implode input_file) ^ " " 
   669                    ^ (Path.implode output_file) ^
   670                    (if dbg then "" else "> /dev/null"))
   671   val  opr = File.read output_file 
   672   val res = parse_csdpoutput opr 
   673  in
   674     ((if dbg then ()
   675       else (File.rm input_file; File.rm output_file ; File.cd current_dir));
   676      (rv,res))
   677  end;
   678 
   679 fun csdp obj mats =
   680  let 
   681   val (rv,res) = run_csdp (!debugging) obj mats 
   682  in
   683    ((if rv = 1 orelse rv = 2 then error "csdp: Problem is infeasible"
   684     else if rv = 3 then writeln "csdp warning: Reduced accuracy"
   685     else if rv <> 0 then error ("csdp: error "^string_of_int rv)
   686     else ());
   687    res)
   688  end;
   689 
   690 (* Try some apparently sensible scaling first. Note that this is purely to   *)
   691 (* get a cleaner translation to floating-point, and doesn't affect any of    *)
   692 (* the results, in principle. In practice it seems a lot better when there   *)
   693 (* are extreme numbers in the original problem.                              *)
   694 
   695   (* Version for (int*int) keys *)
   696 local
   697   fun max_rat x y = if x </ y then y else x
   698   fun common_denominator fld amat acc =
   699       fld (fn (m,c) => fn a => lcm_rat (denominator_rat c) a) amat acc
   700   fun maximal_element fld amat acc =
   701     fld (fn (m,c) => fn maxa => max_rat maxa (abs_rat c)) amat acc 
   702 fun float_of_rat x = let val (a,b) = Rat.quotient_of_rat x
   703                      in Real.fromLargeInt a / Real.fromLargeInt b end;
   704 in
   705 
   706 fun pi_scale_then solver (obj:vector)  mats =
   707  let 
   708   val cd1 = fold_rev (common_denominator Intpairfunc.fold) mats (rat_1)
   709   val cd2 = common_denominator Intfunc.fold (snd obj)  (rat_1) 
   710   val mats' = map (Intpairfunc.mapf (fn x => cd1 */ x)) mats
   711   val obj' = vector_cmul cd2 obj
   712   val max1 = fold_rev (maximal_element Intpairfunc.fold) mats' (rat_0)
   713   val max2 = maximal_element Intfunc.fold (snd obj') (rat_0) 
   714   val scal1 = pow2 (20 - trunc(Math.ln (float_of_rat max1) / Math.ln 2.0))
   715   val scal2 = pow2 (20 - trunc(Math.ln (float_of_rat max2) / Math.ln 2.0)) 
   716   val mats'' = map (Intpairfunc.mapf (fn x => x */ scal1)) mats'
   717   val obj'' = vector_cmul scal2 obj' 
   718  in solver obj'' mats''
   719   end
   720 end;
   721 
   722 (* Try some apparently sensible scaling first. Note that this is purely to   *)
   723 (* get a cleaner translation to floating-point, and doesn't affect any of    *)
   724 (* the results, in principle. In practice it seems a lot better when there   *)
   725 (* are extreme numbers in the original problem.                              *)
   726 
   727   (* Version for (int*int*int) keys *)
   728 local
   729   fun max_rat x y = if x </ y then y else x
   730   fun common_denominator fld amat acc =
   731       fld (fn (m,c) => fn a => lcm_rat (denominator_rat c) a) amat acc
   732   fun maximal_element fld amat acc =
   733     fld (fn (m,c) => fn maxa => max_rat maxa (abs_rat c)) amat acc 
   734 fun float_of_rat x = let val (a,b) = Rat.quotient_of_rat x
   735                      in Real.fromLargeInt a / Real.fromLargeInt b end;
   736 fun int_of_float x = (trunc x handle Overflow => 0 | Domain => 0)
   737 in
   738 
   739 fun tri_scale_then solver (obj:vector)  mats =
   740  let 
   741   val cd1 = fold_rev (common_denominator Inttriplefunc.fold) mats (rat_1)
   742   val cd2 = common_denominator Intfunc.fold (snd obj)  (rat_1) 
   743   val mats' = map (Inttriplefunc.mapf (fn x => cd1 */ x)) mats
   744   val obj' = vector_cmul cd2 obj
   745   val max1 = fold_rev (maximal_element Inttriplefunc.fold) mats' (rat_0)
   746   val max2 = maximal_element Intfunc.fold (snd obj') (rat_0) 
   747   val scal1 = pow2 (20 - int_of_float(Math.ln (float_of_rat max1) / Math.ln 2.0))
   748   val scal2 = pow2 (20 - int_of_float(Math.ln (float_of_rat max2) / Math.ln 2.0)) 
   749   val mats'' = map (Inttriplefunc.mapf (fn x => x */ scal1)) mats'
   750   val obj'' = vector_cmul scal2 obj' 
   751  in solver obj'' mats''
   752   end
   753 end;
   754 
   755 (* Round a vector to "nice" rationals.                                       *)
   756 
   757 fun nice_rational n x = round_rat (n */ x) // n;;
   758 fun nice_vector n ((d,v) : vector) = 
   759  (d, Intfunc.fold (fn (i,c) => fn a => 
   760    let val y = nice_rational n c 
   761    in if c =/ rat_0 then a 
   762       else Intfunc.update (i,y) a end) v Intfunc.undefined):vector
   763 
   764 
   765 (* Reduce linear program to SDP (diagonal matrices) and test with CSDP. This *)
   766 (* one tests A [-1;x1;..;xn] >= 0 (i.e. left column is negated constants).   *)
   767 
   768 fun linear_program_basic a =
   769  let 
   770   val (m,n) = dimensions a
   771   val mats =  map (fn j => diagonal (column j a)) (1 upto n)
   772   val obj = vector_const rat_1 m 
   773   val (rv,res) = run_csdp false obj mats 
   774  in if rv = 1 orelse rv = 2 then false
   775     else if rv = 0 then true
   776     else error "linear_program: An error occurred in the SDP solver"
   777  end;
   778 
   779 (* Alternative interface testing A x >= b for matrix A, vector b.            *)
   780 
   781 fun linear_program a b =
   782  let val (m,n) = dimensions a 
   783  in if dim b <> m then error "linear_program: incompatible dimensions" 
   784     else
   785     let 
   786      val mats = diagonal b :: map (fn j => diagonal (column j a)) (1 upto n)
   787      val obj = vector_const rat_1 m 
   788      val (rv,res) = run_csdp false obj mats 
   789     in if rv = 1 orelse rv = 2 then false
   790        else if rv = 0 then true
   791        else error "linear_program: An error occurred in the SDP solver"
   792     end
   793  end;
   794 
   795 (* Test whether a point is in the convex hull of others. Rather than use     *)
   796 (* computational geometry, express as linear inequalities and call CSDP.     *)
   797 (* This is a bit lazy of me, but it's easy and not such a bottleneck so far. *)
   798 
   799 fun in_convex_hull pts pt =
   800  let 
   801   val pts1 = (1::pt) :: map (fn x => 1::x) pts 
   802   val pts2 = map (fn p => map (fn x => ~x) p @ p) pts1
   803   val n = length pts + 1
   804   val v = 2 * (length pt + 1)
   805   val m = v + n - 1 
   806   val mat = ((m,n),
   807   itern 1 pts2 (fn pts => fn j => itern 1 pts 
   808                (fn x => fn i => Intpairfunc.update ((i,j), Rat.rat_of_int x)))
   809   (iter (1,n) (fn i => Intpairfunc.update((v + i,i+1), rat_1)) 
   810       Intpairfunc.undefined))
   811  in linear_program_basic mat
   812  end;
   813 
   814 (* Filter down a set of points to a minimal set with the same convex hull.   *)
   815 
   816 local
   817  fun augment1 (m::ms) = if in_convex_hull ms m then ms else ms@[m]
   818  fun augment m ms = funpow 3 augment1 (m::ms)
   819 in
   820 fun minimal_convex_hull mons =
   821  let val mons' = fold_rev augment (tl mons) [hd mons] 
   822  in funpow (length mons') augment1 mons'
   823  end
   824 end;
   825 
   826 fun dest_ord f x = is_equal (f x);
   827 
   828 (* Stuff for "equations" ((int*int*int)->num functions).                         *)
   829 
   830 fun tri_equation_cmul c eq =
   831   if c =/ rat_0 then Inttriplefunc.undefined else Inttriplefunc.mapf (fn d => c */ d) eq;
   832 
   833 fun tri_equation_add eq1 eq2 = Inttriplefunc.combine (curry op +/) (fn x => x =/ rat_0) eq1 eq2;
   834 
   835 fun tri_equation_eval assig eq =
   836  let fun value v = Inttriplefunc.apply assig v 
   837  in Inttriplefunc.fold (fn (v, c) => fn a => a +/ value v */ c) eq rat_0
   838  end;
   839 
   840 (* Eliminate among linear equations: return unconstrained variables and      *)
   841 (* assignments for the others in terms of them. We give one pseudo-variable  *)
   842 (* "one" that's used for a constant term.                                    *)
   843 
   844 local
   845   fun extract_first p l = case l of  (* FIXME : use find_first instead *)
   846    [] => error "extract_first"
   847  | h::t => if p h then (h,t) else
   848           let val (k,s) = extract_first p t in (k,h::s) end
   849 fun eliminate vars dun eqs = case vars of 
   850   [] => if forall Inttriplefunc.is_undefined eqs then dun
   851         else raise Unsolvable
   852  | v::vs =>
   853   ((let 
   854     val (eq,oeqs) = extract_first (fn e => Inttriplefunc.defined e v) eqs 
   855     val a = Inttriplefunc.apply eq v
   856     val eq' = tri_equation_cmul ((Rat.neg rat_1) // a) (Inttriplefunc.undefine v eq)
   857     fun elim e =
   858      let val b = Inttriplefunc.tryapplyd e v rat_0 
   859      in if b =/ rat_0 then e else
   860         tri_equation_add e (tri_equation_cmul (Rat.neg b // a) eq)
   861      end
   862    in eliminate vs (Inttriplefunc.update (v,eq') (Inttriplefunc.mapf elim dun)) (map elim oeqs)
   863    end)
   864   handle ERROR _ => eliminate vs dun eqs)
   865 in
   866 fun tri_eliminate_equations one vars eqs =
   867  let 
   868   val assig = eliminate vars Inttriplefunc.undefined eqs
   869   val vs = Inttriplefunc.fold (fn (x, f) => fn a => remove (dest_ord triple_int_ord) one (Inttriplefunc.dom f) @ a) assig []
   870   in (distinct (dest_ord triple_int_ord) vs, assig)
   871   end
   872 end;
   873 
   874 (* Eliminate all variables, in an essentially arbitrary order.               *)
   875 
   876 fun tri_eliminate_all_equations one =
   877  let 
   878   fun choose_variable eq =
   879    let val (v,_) = Inttriplefunc.choose eq 
   880    in if is_equal (triple_int_ord(v,one)) then
   881       let val eq' = Inttriplefunc.undefine v eq 
   882       in if Inttriplefunc.is_undefined eq' then error "choose_variable" 
   883          else fst (Inttriplefunc.choose eq')
   884       end
   885     else v 
   886    end
   887   fun eliminate dun eqs = case eqs of 
   888     [] => dun
   889   | eq::oeqs =>
   890     if Inttriplefunc.is_undefined eq then eliminate dun oeqs else
   891     let val v = choose_variable eq
   892         val a = Inttriplefunc.apply eq v
   893         val eq' = tri_equation_cmul ((Rat.rat_of_int ~1) // a) 
   894                    (Inttriplefunc.undefine v eq)
   895         fun elim e =
   896          let val b = Inttriplefunc.tryapplyd e v rat_0 
   897          in if b =/ rat_0 then e 
   898             else tri_equation_add e (tri_equation_cmul (Rat.neg b // a) eq)
   899          end
   900     in eliminate (Inttriplefunc.update(v, eq') (Inttriplefunc.mapf elim dun)) 
   901                  (map elim oeqs) 
   902     end
   903 in fn eqs =>
   904  let 
   905   val assig = eliminate Inttriplefunc.undefined eqs
   906   val vs = Inttriplefunc.fold (fn (x, f) => fn a => remove (dest_ord triple_int_ord) one (Inttriplefunc.dom f) @ a) assig []
   907  in (distinct (dest_ord triple_int_ord) vs,assig)
   908  end
   909 end;
   910  
   911 (* Solve equations by assigning arbitrary numbers.                           *)
   912 
   913 fun tri_solve_equations one eqs =
   914  let 
   915   val (vars,assigs) = tri_eliminate_all_equations one eqs
   916   val vfn = fold_rev (fn v => Inttriplefunc.update(v,rat_0)) vars 
   917             (Inttriplefunc.onefunc(one, Rat.rat_of_int ~1))
   918   val ass =
   919     Inttriplefunc.combine (curry op +/) (K false) 
   920     (Inttriplefunc.mapf (tri_equation_eval vfn) assigs) vfn 
   921  in if forall (fn e => tri_equation_eval ass e =/ rat_0) eqs
   922     then Inttriplefunc.undefine one ass else raise Sanity
   923  end;
   924 
   925 (* Multiply equation-parametrized poly by regular poly and add accumulator.  *)
   926 
   927 fun tri_epoly_pmul p q acc =
   928  Monomialfunc.fold (fn (m1, c) => fn a =>
   929   Monomialfunc.fold (fn (m2,e) => fn b =>
   930    let val m =  monomial_mul m1 m2
   931        val es = Monomialfunc.tryapplyd b m Inttriplefunc.undefined 
   932    in Monomialfunc.update (m,tri_equation_add (tri_equation_cmul c e) es) b 
   933    end) q a) p acc ;
   934 
   935 (* Usual operations on equation-parametrized poly.                           *)
   936 
   937 fun tri_epoly_cmul c l =
   938   if c =/ rat_0 then Inttriplefunc.undefined else Inttriplefunc.mapf (tri_equation_cmul c) l;;
   939 
   940 val tri_epoly_neg = tri_epoly_cmul (Rat.rat_of_int ~1);
   941 
   942 val tri_epoly_add = Inttriplefunc.combine tri_equation_add Inttriplefunc.is_undefined;
   943 
   944 fun tri_epoly_sub p q = tri_epoly_add p (tri_epoly_neg q);;
   945 
   946 (* Stuff for "equations" ((int*int)->num functions).                         *)
   947 
   948 fun pi_equation_cmul c eq =
   949   if c =/ rat_0 then Inttriplefunc.undefined else Inttriplefunc.mapf (fn d => c */ d) eq;
   950 
   951 fun pi_equation_add eq1 eq2 = Inttriplefunc.combine (curry op +/) (fn x => x =/ rat_0) eq1 eq2;
   952 
   953 fun pi_equation_eval assig eq =
   954  let fun value v = Inttriplefunc.apply assig v 
   955  in Inttriplefunc.fold (fn (v, c) => fn a => a +/ value v */ c) eq rat_0
   956  end;
   957 
   958 (* Eliminate among linear equations: return unconstrained variables and      *)
   959 (* assignments for the others in terms of them. We give one pseudo-variable  *)
   960 (* "one" that's used for a constant term.                                    *)
   961 
   962 local
   963 fun extract_first p l = case l of 
   964    [] => error "extract_first"
   965  | h::t => if p h then (h,t) else
   966           let val (k,s) = extract_first p t in (k,h::s) end
   967 fun eliminate vars dun eqs = case vars of 
   968   [] => if forall Inttriplefunc.is_undefined eqs then dun
   969         else raise Unsolvable
   970  | v::vs =>
   971    let 
   972     val (eq,oeqs) = extract_first (fn e => Inttriplefunc.defined e v) eqs 
   973     val a = Inttriplefunc.apply eq v
   974     val eq' = pi_equation_cmul ((Rat.neg rat_1) // a) (Inttriplefunc.undefine v eq)
   975     fun elim e =
   976      let val b = Inttriplefunc.tryapplyd e v rat_0 
   977      in if b =/ rat_0 then e else
   978         pi_equation_add e (pi_equation_cmul (Rat.neg b // a) eq)
   979      end
   980    in eliminate vs (Inttriplefunc.update (v,eq') (Inttriplefunc.mapf elim dun)) (map elim oeqs)
   981    end
   982   handle ERROR _ => eliminate vs dun eqs
   983 in
   984 fun pi_eliminate_equations one vars eqs =
   985  let 
   986   val assig = eliminate vars Inttriplefunc.undefined eqs
   987   val vs = Inttriplefunc.fold (fn (x, f) => fn a => remove (dest_ord triple_int_ord) one (Inttriplefunc.dom f) @ a) assig []
   988   in (distinct (dest_ord triple_int_ord) vs, assig)
   989   end
   990 end;
   991 
   992 (* Eliminate all variables, in an essentially arbitrary order.               *)
   993 
   994 fun pi_eliminate_all_equations one =
   995  let 
   996   fun choose_variable eq =
   997    let val (v,_) = Inttriplefunc.choose eq 
   998    in if is_equal (triple_int_ord(v,one)) then
   999       let val eq' = Inttriplefunc.undefine v eq 
  1000       in if Inttriplefunc.is_undefined eq' then error "choose_variable" 
  1001          else fst (Inttriplefunc.choose eq')
  1002       end
  1003     else v 
  1004    end
  1005   fun eliminate dun eqs = case eqs of 
  1006     [] => dun
  1007   | eq::oeqs =>
  1008     if Inttriplefunc.is_undefined eq then eliminate dun oeqs else
  1009     let val v = choose_variable eq
  1010         val a = Inttriplefunc.apply eq v
  1011         val eq' = pi_equation_cmul ((Rat.rat_of_int ~1) // a) 
  1012                    (Inttriplefunc.undefine v eq)
  1013         fun elim e =
  1014          let val b = Inttriplefunc.tryapplyd e v rat_0 
  1015          in if b =/ rat_0 then e 
  1016             else pi_equation_add e (pi_equation_cmul (Rat.neg b // a) eq)
  1017          end
  1018     in eliminate (Inttriplefunc.update(v, eq') (Inttriplefunc.mapf elim dun)) 
  1019                  (map elim oeqs) 
  1020     end
  1021 in fn eqs =>
  1022  let 
  1023   val assig = eliminate Inttriplefunc.undefined eqs
  1024   val vs = Inttriplefunc.fold (fn (x, f) => fn a => remove (dest_ord triple_int_ord) one (Inttriplefunc.dom f) @ a) assig []
  1025  in (distinct (dest_ord triple_int_ord) vs,assig)
  1026  end
  1027 end;
  1028  
  1029 (* Solve equations by assigning arbitrary numbers.                           *)
  1030 
  1031 fun pi_solve_equations one eqs =
  1032  let 
  1033   val (vars,assigs) = pi_eliminate_all_equations one eqs
  1034   val vfn = fold_rev (fn v => Inttriplefunc.update(v,rat_0)) vars 
  1035             (Inttriplefunc.onefunc(one, Rat.rat_of_int ~1))
  1036   val ass =
  1037     Inttriplefunc.combine (curry op +/) (K false) 
  1038     (Inttriplefunc.mapf (pi_equation_eval vfn) assigs) vfn 
  1039  in if forall (fn e => pi_equation_eval ass e =/ rat_0) eqs
  1040     then Inttriplefunc.undefine one ass else raise Sanity
  1041  end;
  1042 
  1043 (* Multiply equation-parametrized poly by regular poly and add accumulator.  *)
  1044 
  1045 fun pi_epoly_pmul p q acc =
  1046  Monomialfunc.fold (fn (m1, c) => fn a =>
  1047   Monomialfunc.fold (fn (m2,e) => fn b =>
  1048    let val m =  monomial_mul m1 m2
  1049        val es = Monomialfunc.tryapplyd b m Inttriplefunc.undefined 
  1050    in Monomialfunc.update (m,pi_equation_add (pi_equation_cmul c e) es) b 
  1051    end) q a) p acc ;
  1052 
  1053 (* Usual operations on equation-parametrized poly.                           *)
  1054 
  1055 fun pi_epoly_cmul c l =
  1056   if c =/ rat_0 then Inttriplefunc.undefined else Inttriplefunc.mapf (pi_equation_cmul c) l;;
  1057 
  1058 val pi_epoly_neg = pi_epoly_cmul (Rat.rat_of_int ~1);
  1059 
  1060 val pi_epoly_add = Inttriplefunc.combine pi_equation_add Inttriplefunc.is_undefined;
  1061 
  1062 fun pi_epoly_sub p q = pi_epoly_add p (pi_epoly_neg q);;
  1063 
  1064 fun allpairs f l1 l2 =  fold_rev (fn x => (curry (op @)) (map (f x) l2)) l1 [];
  1065 
  1066 (* Hence produce the "relevant" monomials: those whose squares lie in the    *)
  1067 (* Newton polytope of the monomials in the input. (This is enough according  *)
  1068 (* to Reznik: "Extremal PSD forms with few terms", Duke Math. Journal,       *)
  1069 (* vol 45, pp. 363--374, 1978.                                               *)
  1070 (*                                                                           *)
  1071 (* These are ordered in sort of decreasing degree. In particular the         *)
  1072 (* constant monomial is last; this gives an order in diagonalization of the  *)
  1073 (* quadratic form that will tend to display constants.                       *)
  1074 
  1075 fun newton_polytope pol =
  1076  let 
  1077   val vars = poly_variables pol
  1078   val mons = map (fn m => map (fn x => monomial_degree x m) vars) 
  1079              (Monomialfunc.dom pol)
  1080   val ds = map (fn x => (degree x pol + 1) div 2) vars
  1081   val all = fold_rev (fn n => allpairs cons (0 upto n)) ds [[]]
  1082   val mons' = minimal_convex_hull mons
  1083   val all' =
  1084     filter (fn m => in_convex_hull mons' (map (fn x => 2 * x) m)) all 
  1085  in map (fn m => fold_rev2 (fn v => fn i => fn a => if i = 0 then a else Ctermfunc.update (v,i) a)
  1086                         vars m monomial_1) (rev all')
  1087  end;
  1088 
  1089 (* Diagonalize (Cholesky/LDU) the matrix corresponding to a quadratic form.  *)
  1090 
  1091 local
  1092 fun diagonalize n i m =
  1093  if Intpairfunc.is_undefined (snd m) then [] 
  1094  else
  1095   let val a11 = Intpairfunc.tryapplyd (snd m) (i,i) rat_0 
  1096   in if a11 </ rat_0 then error "diagonalize: not PSD"
  1097     else if a11 =/ rat_0 then
  1098           if Intfunc.is_undefined (snd (row i m)) then diagonalize n (i + 1) m
  1099           else error "diagonalize: not PSD ___ "
  1100     else
  1101      let 
  1102       val v = row i m
  1103       val v' = (fst v, Intfunc.fold (fn (i, c) => fn a => 
  1104        let val y = c // a11 
  1105        in if y = rat_0 then a else Intfunc.update (i,y) a 
  1106        end)  (snd v) Intfunc.undefined)
  1107       fun upt0 x y a = if y = rat_0 then a else Intpairfunc.update (x,y) a
  1108       val m' =
  1109       ((n,n),
  1110       iter (i+1,n) (fn j =>
  1111           iter (i+1,n) (fn k =>
  1112               (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))))
  1113           Intpairfunc.undefined)
  1114      in (a11,v')::diagonalize n (i + 1) m' 
  1115      end
  1116   end
  1117 in
  1118 fun diag m =
  1119  let 
  1120    val nn = dimensions m 
  1121    val n = fst nn 
  1122  in if snd nn <> n then error "diagonalize: non-square matrix" 
  1123     else diagonalize n 1 m
  1124  end
  1125 end;
  1126 
  1127 fun gcd_rat a b = Rat.rat_of_int (Integer.gcd (int_of_rat a) (int_of_rat b));
  1128 
  1129 (* Adjust a diagonalization to collect rationals at the start.               *)
  1130   (* FIXME : Potentially polymorphic keys, but here only: integers!! *)
  1131 local
  1132  fun upd0 x y a = if y =/ rat_0 then a else Intfunc.update(x,y) a;
  1133  fun mapa f (d,v) = 
  1134   (d, Intfunc.fold (fn (i,c) => fn a => upd0 i (f c) a) v Intfunc.undefined)
  1135  fun adj (c,l) =
  1136  let val a = 
  1137   Intfunc.fold (fn (i,c) => fn a => lcm_rat a (denominator_rat c)) 
  1138     (snd l) rat_1 //
  1139   Intfunc.fold (fn (i,c) => fn a => gcd_rat a (numerator_rat c)) 
  1140     (snd l) rat_0
  1141   in ((c // (a */ a)),mapa (fn x => a */ x) l)
  1142   end
  1143 in
  1144 fun deration d = if null d then (rat_0,d) else
  1145  let val d' = map adj d
  1146      val a = fold (lcm_rat o denominator_rat o fst) d' rat_1 //
  1147           fold (gcd_rat o numerator_rat o fst) d' rat_0 
  1148  in ((rat_1 // a),map (fn (c,l) => (a */ c,l)) d')
  1149  end
  1150 end;
  1151  
  1152 (* Enumeration of monomials with given multidegree bound.                    *)
  1153 
  1154 fun enumerate_monomials d vars = 
  1155  if d < 0 then []
  1156  else if d = 0 then [Ctermfunc.undefined]
  1157  else if null vars then [monomial_1] else
  1158  let val alts =
  1159   map (fn k => let val oths = enumerate_monomials (d - k) (tl vars) 
  1160                in map (fn ks => if k = 0 then ks else Ctermfunc.update (hd vars, k) ks) oths end) (0 upto d) 
  1161  in fold1 (curry op @) alts
  1162  end;
  1163 
  1164 (* Enumerate products of distinct input polys with degree <= d.              *)
  1165 (* We ignore any constant input polynomials.                                 *)
  1166 (* Give the output polynomial and a record of how it was derived.            *)
  1167 
  1168 local
  1169  open RealArith
  1170 in
  1171 fun enumerate_products d pols =
  1172 if d = 0 then [(poly_const rat_1,Rational_lt rat_1)] 
  1173 else if d < 0 then [] else
  1174 case pols of 
  1175    [] => [(poly_const rat_1,Rational_lt rat_1)]
  1176  | (p,b)::ps => 
  1177     let val e = multidegree p 
  1178     in if e = 0 then enumerate_products d ps else
  1179        enumerate_products d ps @
  1180        map (fn (q,c) => (poly_mul p q,Product(b,c)))
  1181          (enumerate_products (d - e) ps)
  1182     end
  1183 end;
  1184 
  1185 (* Convert regular polynomial. Note that we treat (0,0,0) as -1.             *)
  1186 
  1187 fun epoly_of_poly p =
  1188   Monomialfunc.fold (fn (m,c) => fn a => Monomialfunc.update (m, Inttriplefunc.onefunc ((0,0,0), Rat.neg c)) a) p Monomialfunc.undefined;
  1189 
  1190 (* String for block diagonal matrix numbered k.                              *)
  1191 
  1192 fun sdpa_of_blockdiagonal k m =
  1193  let 
  1194   val pfx = string_of_int k ^" "
  1195   val ents =
  1196     Inttriplefunc.fold 
  1197       (fn ((b,i,j),c) => fn a => if i > j then a else ((b,i,j),c)::a) 
  1198       m [] 
  1199   val entss = sort (increasing fst triple_int_ord) ents 
  1200  in fold_rev (fn ((b,i,j),c) => fn a =>
  1201      pfx ^ string_of_int b ^ " " ^ string_of_int i ^ " " ^ string_of_int j ^
  1202      " " ^ decimalize 20 c ^ "\n" ^ a) entss ""
  1203  end;
  1204 
  1205 (* SDPA for problem using block diagonal (i.e. multiple SDPs)                *)
  1206 
  1207 fun sdpa_of_blockproblem comment nblocks blocksizes obj mats =
  1208  let val m = length mats - 1 
  1209  in  "\"" ^ comment ^ "\"\n" ^
  1210   string_of_int m ^ "\n" ^
  1211   string_of_int nblocks ^ "\n" ^
  1212   (fold1 (fn s => fn t => s^" "^t) (map string_of_int blocksizes)) ^
  1213   "\n" ^
  1214   sdpa_of_vector obj ^
  1215   fold_rev2 (fn k => fn m => fn a => sdpa_of_blockdiagonal (k - 1) m ^ a)
  1216     (1 upto length mats) mats ""
  1217  end;
  1218 
  1219 (* Hence run CSDP on a problem in block diagonal form.                       *)
  1220 
  1221 fun run_csdp dbg nblocks blocksizes obj mats =
  1222  let 
  1223   val input_file = tmp_file "sos" ".dat-s" 
  1224   val output_file = tmp_file "sos" ".out"
  1225   val params_file = tmp_file "param" ".csdp" 
  1226   val _ = File.write input_file
  1227    (sdpa_of_blockproblem "" nblocks blocksizes obj mats)
  1228   val _ = File.write params_file csdp_params
  1229   val current_dir = File.pwd()
  1230   val _ = File.cd (Path.variable "ISABELLE_TMP")
  1231   val rv = system ("csdp "^(Path.implode input_file) ^ " " 
  1232                    ^ (Path.implode output_file) ^
  1233                    (if dbg then "" else "> /dev/null"))
  1234   val  opr = File.read output_file 
  1235   val res = parse_csdpoutput opr 
  1236  in
  1237    ((if dbg then ()
  1238      else (File.rm input_file ; File.rm output_file ; File.cd current_dir));
  1239     (rv,res))
  1240  end;
  1241 
  1242 fun csdp nblocks blocksizes obj mats =
  1243  let 
  1244   val (rv,res) = run_csdp (!debugging) nblocks blocksizes obj mats 
  1245  in ((if rv = 1 orelse rv = 2 then error "csdp: Problem is infeasible"
  1246      else if rv = 3 then writeln "csdp warning: Reduced accuracy"
  1247      else if rv <> 0 then error ("csdp: error "^string_of_int rv)
  1248      else ());
  1249      res)
  1250  end;
  1251 
  1252 (* 3D versions of matrix operations to consider blocks separately.           *)
  1253 
  1254 val bmatrix_add = Inttriplefunc.combine (curry op +/) (fn x => x =/ rat_0);
  1255 fun bmatrix_cmul c bm =
  1256   if c =/ rat_0 then Inttriplefunc.undefined
  1257   else Inttriplefunc.mapf (fn x => c */ x) bm;
  1258 
  1259 val bmatrix_neg = bmatrix_cmul (Rat.rat_of_int ~1);
  1260 fun bmatrix_sub m1 m2 = bmatrix_add m1 (bmatrix_neg m2);;
  1261 
  1262 (* Smash a block matrix into components.                                     *)
  1263 
  1264 fun blocks blocksizes bm =
  1265  map (fn (bs,b0) =>
  1266       let val m = Inttriplefunc.fold
  1267           (fn ((b,i,j),c) => fn a => if b = b0 then Intpairfunc.update ((i,j),c) a else a) bm Intpairfunc.undefined
  1268           val d = Intpairfunc.fold (fn ((i,j),c) => fn a => max a (max i j)) m 0 
  1269       in (((bs,bs),m):matrix) end)
  1270  (blocksizes ~~ (1 upto length blocksizes));;
  1271 
  1272 (* FIXME : Get rid of this !!!*)
  1273 fun tryfind f [] = error "tryfind"
  1274   | tryfind f (x::xs) = (f x handle ERROR _ => tryfind f xs);
  1275 
  1276 
  1277 (* Positiv- and Nullstellensatz. Flag "linf" forces a linear representation. *)
  1278 
  1279 
  1280 local
  1281  open RealArith
  1282 in
  1283 fun real_positivnullstellensatz_general linf d eqs leqs pol =
  1284 let 
  1285  val vars = fold_rev (curry (gen_union (op aconvc)) o poly_variables) 
  1286               (pol::eqs @ map fst leqs) []
  1287  val monoid = if linf then 
  1288       (poly_const rat_1,Rational_lt rat_1)::
  1289       (filter (fn (p,c) => multidegree p <= d) leqs)
  1290     else enumerate_products d leqs
  1291  val nblocks = length monoid
  1292  fun mk_idmultiplier k p =
  1293   let 
  1294    val e = d - multidegree p
  1295    val mons = enumerate_monomials e vars
  1296    val nons = mons ~~ (1 upto length mons) 
  1297   in (mons,
  1298       fold_rev (fn (m,n) => Monomialfunc.update(m,Inttriplefunc.onefunc((~k,~n,n),rat_1))) nons Monomialfunc.undefined)
  1299   end
  1300 
  1301  fun mk_sqmultiplier k (p,c) =
  1302   let 
  1303    val e = (d - multidegree p) div 2
  1304    val mons = enumerate_monomials e vars
  1305    val nons = mons ~~ (1 upto length mons) 
  1306   in (mons, 
  1307       fold_rev (fn (m1,n1) =>
  1308        fold_rev (fn (m2,n2) => fn  a =>
  1309         let val m = monomial_mul m1 m2 
  1310         in if n1 > n2 then a else
  1311           let val c = if n1 = n2 then rat_1 else rat_2
  1312               val e = Monomialfunc.tryapplyd a m Inttriplefunc.undefined 
  1313           in Monomialfunc.update(m, tri_equation_add (Inttriplefunc.onefunc((k,n1,n2), c)) e) a
  1314           end
  1315         end)  nons)
  1316        nons Monomialfunc.undefined)
  1317   end
  1318 
  1319   val (sqmonlist,sqs) = split_list (map2 mk_sqmultiplier (1 upto length monoid) monoid)
  1320   val (idmonlist,ids) =  split_list(map2 mk_idmultiplier (1 upto length eqs) eqs)
  1321   val blocksizes = map length sqmonlist
  1322   val bigsum =
  1323     fold_rev2 (fn p => fn q => fn a => tri_epoly_pmul p q a) eqs ids
  1324             (fold_rev2 (fn (p,c) => fn s => fn a => tri_epoly_pmul p s a) monoid sqs
  1325                      (epoly_of_poly(poly_neg pol)))
  1326   val eqns = Monomialfunc.fold (fn (m,e) => fn a => e::a) bigsum []
  1327   val (pvs,assig) = tri_eliminate_all_equations (0,0,0) eqns
  1328   val qvars = (0,0,0)::pvs
  1329   val allassig = fold_rev (fn v => Inttriplefunc.update(v,(Inttriplefunc.onefunc(v,rat_1)))) pvs assig
  1330   fun mk_matrix v =
  1331     Inttriplefunc.fold (fn ((b,i,j), ass) => fn m => 
  1332         if b < 0 then m else
  1333          let val c = Inttriplefunc.tryapplyd ass v rat_0
  1334          in if c = rat_0 then m else
  1335             Inttriplefunc.update ((b,j,i), c) (Inttriplefunc.update ((b,i,j), c) m)
  1336          end)
  1337           allassig Inttriplefunc.undefined
  1338   val diagents = Inttriplefunc.fold
  1339     (fn ((b,i,j), e) => fn a => if b > 0 andalso i = j then tri_equation_add e a else a)
  1340     allassig Inttriplefunc.undefined
  1341 
  1342   val mats = map mk_matrix qvars
  1343   val obj = (length pvs,
  1344             itern 1 pvs (fn v => fn i => Intfunc.updatep iszero (i,Inttriplefunc.tryapplyd diagents v rat_0))
  1345                         Intfunc.undefined)
  1346   val raw_vec = if null pvs then vector_0 0
  1347                 else tri_scale_then (csdp nblocks blocksizes) obj mats
  1348   fun int_element (d,v) i = Intfunc.tryapplyd v i rat_0
  1349   fun cterm_element (d,v) i = Ctermfunc.tryapplyd v i rat_0
  1350 
  1351   fun find_rounding d =
  1352    let 
  1353     val _ = if !debugging 
  1354           then writeln ("Trying rounding with limit "^Rat.string_of_rat d ^ "\n") 
  1355           else ()
  1356     val vec = nice_vector d raw_vec
  1357     val blockmat = iter (1,dim vec)
  1358      (fn i => fn a => bmatrix_add (bmatrix_cmul (int_element vec i) (nth mats i)) a)
  1359      (bmatrix_neg (nth mats 0))
  1360     val allmats = blocks blocksizes blockmat 
  1361    in (vec,map diag allmats)
  1362    end
  1363   val (vec,ratdias) =
  1364     if null pvs then find_rounding rat_1
  1365     else tryfind find_rounding (map Rat.rat_of_int (1 upto 31) @
  1366                                 map pow2 (5 upto 66))
  1367   val newassigs =
  1368     fold_rev (fn k => Inttriplefunc.update (nth pvs (k - 1), int_element vec k))
  1369            (1 upto dim vec) (Inttriplefunc.onefunc ((0,0,0), Rat.rat_of_int ~1))
  1370   val finalassigs =
  1371     Inttriplefunc.fold (fn (v,e) => fn a => Inttriplefunc.update(v, tri_equation_eval newassigs e) a) allassig newassigs
  1372   fun poly_of_epoly p =
  1373     Monomialfunc.fold (fn (v,e) => fn a => Monomialfunc.updatep iszero (v,tri_equation_eval finalassigs e) a)
  1374           p Monomialfunc.undefined
  1375   fun  mk_sos mons =
  1376    let fun mk_sq (c,m) =
  1377     (c,fold_rev (fn k=> fn a => Monomialfunc.updatep iszero (nth mons (k - 1), int_element m k) a)
  1378                  (1 upto length mons) Monomialfunc.undefined)
  1379    in map mk_sq
  1380    end
  1381   val sqs = map2 mk_sos sqmonlist ratdias
  1382   val cfs = map poly_of_epoly ids
  1383   val msq = filter (fn (a,b) => not (null b)) (map2 pair monoid sqs)
  1384   fun eval_sq sqs = fold_rev (fn (c,q) => poly_add (poly_cmul c (poly_mul q q))) sqs poly_0
  1385   val sanity =
  1386     fold_rev (fn ((p,c),s) => poly_add (poly_mul p (eval_sq s))) msq
  1387            (fold_rev2 (fn p => fn q => poly_add (poly_mul p q)) cfs eqs
  1388                     (poly_neg pol))
  1389 
  1390 in if not(Monomialfunc.is_undefined sanity) then raise Sanity else
  1391   (cfs,map (fn (a,b) => (snd a,b)) msq)
  1392  end
  1393 
  1394 
  1395 end;
  1396 
  1397 (* Iterative deepening.                                                      *)
  1398 
  1399 fun deepen f n = 
  1400   (writeln ("Searching with depth limit " ^ string_of_int n) ; (f n handle ERROR s => (writeln ("failed with message: " ^ s) ; deepen f (n+1))))
  1401 
  1402 (* The ordering so we can create canonical HOL polynomials.                  *)
  1403 
  1404 fun dest_monomial mon = sort (increasing fst cterm_ord) (Ctermfunc.graph mon);
  1405 
  1406 fun monomial_order (m1,m2) =
  1407  if Ctermfunc.is_undefined m2 then LESS 
  1408  else if Ctermfunc.is_undefined m1 then GREATER 
  1409  else
  1410   let val mon1 = dest_monomial m1 
  1411       val mon2 = dest_monomial m2
  1412       val deg1 = fold (curry op + o snd) mon1 0
  1413       val deg2 = fold (curry op + o snd) mon2 0 
  1414   in if deg1 < deg2 then GREATER else if deg1 > deg2 then LESS
  1415      else list_ord (prod_ord cterm_ord int_ord) (mon1,mon2)
  1416   end;
  1417 
  1418 fun dest_poly p =
  1419   map (fn (m,c) => (c,dest_monomial m))
  1420       (sort (prod_ord monomial_order (K EQUAL)) (Monomialfunc.graph p));
  1421 
  1422 (* Map back polynomials and their composites to HOL.                         *)
  1423 
  1424 local
  1425  open Thm Numeral RealArith
  1426 in
  1427 
  1428 fun cterm_of_varpow x k = if k = 1 then x else capply (capply @{cterm "op ^ :: real => _"} x) 
  1429   (mk_cnumber @{ctyp nat} k)
  1430 
  1431 fun cterm_of_monomial m = 
  1432  if Ctermfunc.is_undefined m then @{cterm "1::real"} 
  1433  else 
  1434   let 
  1435    val m' = dest_monomial m
  1436    val vps = fold_rev (fn (x,k) => cons (cterm_of_varpow x k)) m' [] 
  1437   in fold1 (fn s => fn t => capply (capply @{cterm "op * :: real => _"} s) t) vps
  1438   end
  1439 
  1440 fun cterm_of_cmonomial (m,c) = if Ctermfunc.is_undefined m then cterm_of_rat c
  1441     else if c = Rat.one then cterm_of_monomial m
  1442     else capply (capply @{cterm "op *::real => _"} (cterm_of_rat c)) (cterm_of_monomial m);
  1443 
  1444 fun cterm_of_poly p = 
  1445  if Monomialfunc.is_undefined p then @{cterm "0::real"} 
  1446  else
  1447   let 
  1448    val cms = map cterm_of_cmonomial
  1449      (sort (prod_ord monomial_order (K EQUAL)) (Monomialfunc.graph p))
  1450   in fold1 (fn t1 => fn t2 => capply(capply @{cterm "op + :: real => _"} t1) t2) cms
  1451   end;
  1452 
  1453 fun cterm_of_sqterm (c,p) = Product(Rational_lt c,Square(cterm_of_poly p));
  1454 
  1455 fun cterm_of_sos (pr,sqs) = if null sqs then pr
  1456   else Product(pr,fold1 (fn a => fn b => Sum(a,b)) (map cterm_of_sqterm sqs));
  1457 
  1458 end
  1459 
  1460 (* Interface to HOL.                                                         *)
  1461 local
  1462   open Thm Conv RealArith
  1463   val concl = dest_arg o cprop_of
  1464   fun simple_cterm_ord t u = TermOrd.fast_term_ord (term_of t, term_of u) = LESS
  1465 in
  1466   (* FIXME: Replace tryfind by get_first !! *)
  1467 fun real_nonlinear_prover ctxt =
  1468  let 
  1469   val {add,mul,neg,pow,sub,main} =  Normalizer.semiring_normalizers_ord_wrapper ctxt
  1470       (valOf (NormalizerData.match ctxt @{cterm "(0::real) + 1"})) 
  1471      simple_cterm_ord
  1472   val (real_poly_add_conv,real_poly_mul_conv,real_poly_neg_conv,
  1473        real_poly_pow_conv,real_poly_sub_conv,real_poly_conv) = (add,mul,neg,pow,sub,main)
  1474   fun mainf  translator (eqs,les,lts) = 
  1475   let 
  1476    val eq0 = map (poly_of_term o dest_arg1 o concl) eqs
  1477    val le0 = map (poly_of_term o dest_arg o concl) les
  1478    val lt0 = map (poly_of_term o dest_arg o concl) lts
  1479    val eqp0 = map (fn (t,i) => (t,Axiom_eq i)) (eq0 ~~ (0 upto (length eq0 - 1)))
  1480    val lep0 = map (fn (t,i) => (t,Axiom_le i)) (le0 ~~ (0 upto (length le0 - 1)))
  1481    val ltp0 = map (fn (t,i) => (t,Axiom_lt i)) (lt0 ~~ (0 upto (length lt0 - 1)))
  1482    val (keq,eq) = List.partition (fn (p,_) => multidegree p = 0) eqp0
  1483    val (klep,lep) = List.partition (fn (p,_) => multidegree p = 0) lep0
  1484    val (kltp,ltp) = List.partition (fn (p,_) => multidegree p = 0) ltp0
  1485    fun trivial_axiom (p,ax) =
  1486     case ax of
  1487        Axiom_eq n => if eval Ctermfunc.undefined p <>/ Rat.zero then nth eqs n 
  1488                      else error "trivial_axiom: Not a trivial axiom"
  1489      | Axiom_le n => if eval Ctermfunc.undefined p </ Rat.zero then nth les n 
  1490                      else error "trivial_axiom: Not a trivial axiom"
  1491      | Axiom_lt n => if eval Ctermfunc.undefined p <=/ Rat.zero then nth lts n 
  1492                      else error "trivial_axiom: Not a trivial axiom"
  1493      | _ => error "trivial_axiom: Not a trivial axiom"
  1494    in 
  1495   ((let val th = tryfind trivial_axiom (keq @ klep @ kltp)
  1496    in fconv_rule (arg_conv (arg1_conv real_poly_conv) then_conv field_comp_conv) th end)
  1497    handle ERROR _ => (
  1498     let 
  1499      val pol = fold_rev poly_mul (map fst ltp) (poly_const Rat.one)
  1500      val leq = lep @ ltp
  1501      fun tryall d =
  1502       let val e = multidegree pol
  1503           val k = if e = 0 then 0 else d div e
  1504           val eq' = map fst eq 
  1505       in tryfind (fn i => (d,i,real_positivnullstellensatz_general false d eq' leq
  1506                             (poly_neg(poly_pow pol i))))
  1507               (0 upto k)
  1508       end
  1509     val (d,i,(cert_ideal,cert_cone)) = deepen tryall 0
  1510     val proofs_ideal =
  1511       map2 (fn q => fn (p,ax) => Eqmul(cterm_of_poly q,ax)) cert_ideal eq
  1512     val proofs_cone = map cterm_of_sos cert_cone
  1513     val proof_ne = if null ltp then Rational_lt Rat.one else
  1514       let val p = fold1 (fn s => fn t => Product(s,t)) (map snd ltp) 
  1515       in  funpow i (fn q => Product(p,q)) (Rational_lt Rat.one)
  1516       end
  1517     val proof = fold1 (fn s => fn t => Sum(s,t))
  1518                            (proof_ne :: proofs_ideal @ proofs_cone) 
  1519     in writeln "Translating proof certificate to HOL";
  1520        translator (eqs,les,lts) proof
  1521     end))
  1522    end
  1523  in mainf end
  1524 end
  1525 
  1526 fun C f x y = f y x;
  1527   (* FIXME : This is very bad!!!*)
  1528 fun subst_conv eqs t = 
  1529  let 
  1530   val t' = fold (Thm.cabs o Thm.lhs_of) eqs t
  1531  in Conv.fconv_rule (Thm.beta_conversion true) (fold (C combination) eqs (reflexive t'))
  1532  end
  1533 
  1534 (* A wrapper that tries to substitute away variables first.                  *)
  1535 
  1536 local
  1537  open Thm Conv RealArith
  1538   fun simple_cterm_ord t u = TermOrd.fast_term_ord (term_of t, term_of u) = LESS
  1539  val concl = dest_arg o cprop_of
  1540  val shuffle1 = 
  1541    fconv_rule (rewr_conv @{lemma "(a + x == y) == (x == y - (a::real))" by (atomize (full)) (simp add: ring_simps) })
  1542  val shuffle2 =
  1543     fconv_rule (rewr_conv @{lemma "(x + a == y) ==  (x == y - (a::real))" by (atomize (full)) (simp add: ring_simps)})
  1544  fun substitutable_monomial fvs tm = case term_of tm of
  1545     Free(_,@{typ real}) => if not (member (op aconvc) fvs tm) then (Rat.one,tm) 
  1546                            else error "substitutable_monomial"
  1547   | @{term "op * :: real => _"}$c$(t as Free _ ) => 
  1548      if is_ratconst (dest_arg1 tm) andalso not (member (op aconvc) fvs (dest_arg tm))
  1549          then (dest_ratconst (dest_arg1 tm),dest_arg tm) else error "substitutable_monomial"
  1550   | @{term "op + :: real => _"}$s$t => 
  1551        (substitutable_monomial (add_cterm_frees (dest_arg tm) fvs) (dest_arg1 tm)
  1552         handle ERROR _ => substitutable_monomial (add_cterm_frees (dest_arg1 tm) fvs) (dest_arg tm))
  1553   | _ => error "substitutable_monomial"
  1554 
  1555   fun isolate_variable v th = 
  1556    let val w = dest_arg1 (cprop_of th)
  1557    in if v aconvc w then th
  1558       else case term_of w of
  1559            @{term "op + :: real => _"}$s$t => 
  1560               if dest_arg1 w aconvc v then shuffle2 th 
  1561               else isolate_variable v (shuffle1 th)
  1562           | _ => error "isolate variable : This should not happen?"
  1563    end 
  1564 in
  1565 
  1566 fun real_nonlinear_subst_prover ctxt =
  1567  let 
  1568   val {add,mul,neg,pow,sub,main} =  Normalizer.semiring_normalizers_ord_wrapper ctxt
  1569       (valOf (NormalizerData.match ctxt @{cterm "(0::real) + 1"})) 
  1570      simple_cterm_ord
  1571 
  1572   val (real_poly_add_conv,real_poly_mul_conv,real_poly_neg_conv,
  1573        real_poly_pow_conv,real_poly_sub_conv,real_poly_conv) = (add,mul,neg,pow,sub,main)
  1574 
  1575   fun make_substitution th =
  1576    let 
  1577     val (c,v) = substitutable_monomial [] (dest_arg1(concl th))
  1578     val th1 = Drule.arg_cong_rule (capply @{cterm "op * :: real => _"} (cterm_of_rat (Rat.inv c))) (mk_meta_eq th)
  1579     val th2 = fconv_rule (binop_conv real_poly_mul_conv) th1
  1580    in fconv_rule (arg_conv real_poly_conv) (isolate_variable v th2)
  1581    end
  1582    fun oprconv cv ct = 
  1583     let val g = Thm.dest_fun2 ct
  1584     in if g aconvc @{cterm "op <= :: real => _"} 
  1585          orelse g aconvc @{cterm "op < :: real => _"} 
  1586        then arg_conv cv ct else arg1_conv cv ct
  1587     end
  1588   fun mainf translator =
  1589    let 
  1590     fun substfirst(eqs,les,lts) =
  1591       ((let 
  1592            val eth = tryfind make_substitution eqs
  1593            val modify = fconv_rule (arg_conv (oprconv(subst_conv [eth] then_conv real_poly_conv)))
  1594        in  substfirst
  1595              (filter_out (fn t => (Thm.dest_arg1 o Thm.dest_arg o cprop_of) t 
  1596                                    aconvc @{cterm "0::real"}) (map modify eqs),
  1597                                    map modify les,map modify lts)
  1598        end)
  1599        handle ERROR  _ => real_nonlinear_prover ctxt translator (rev eqs, rev les, rev lts))
  1600     in substfirst
  1601    end
  1602 
  1603 
  1604  in mainf
  1605  end
  1606 
  1607 (* Overall function. *)
  1608 
  1609 fun real_sos ctxt t = gen_prover_real_arith ctxt (real_nonlinear_subst_prover ctxt) t;
  1610 end;
  1611 
  1612 (* A tactic *)
  1613 fun strip_all ct = 
  1614  case term_of ct of 
  1615   Const("all",_) $ Abs (xn,xT,p) => 
  1616    let val (a,(v,t')) = (apsnd (Thm.dest_abs (SOME xn)) o Thm.dest_comb) ct
  1617    in apfst (cons v) (strip_all t')
  1618    end
  1619 | _ => ([],ct)
  1620 
  1621 fun core_sos_conv ctxt t = Drule.arg_cong_rule @{cterm Trueprop} (real_sos ctxt (Thm.dest_arg t) RS @{thm Eq_TrueI})
  1622 
  1623 val known_sos_constants = 
  1624   [@{term "op ==>"}, @{term "Trueprop"}, 
  1625    @{term "op -->"}, @{term "op &"}, @{term "op |"}, 
  1626    @{term "Not"}, @{term "op = :: bool => _"}, 
  1627    @{term "All :: (real => _) => _"}, @{term "Ex :: (real => _) => _"}, 
  1628    @{term "op = :: real => _"}, @{term "op < :: real => _"}, 
  1629    @{term "op <= :: real => _"}, 
  1630    @{term "op + :: real => _"}, @{term "op - :: real => _"}, 
  1631    @{term "op * :: real => _"}, @{term "uminus :: real => _"}, 
  1632    @{term "op / :: real => _"}, @{term "inverse :: real => _"},
  1633    @{term "op ^ :: real => _"}, @{term "abs :: real => _"}, 
  1634    @{term "min :: real => _"}, @{term "max :: real => _"},
  1635    @{term "0::real"}, @{term "1::real"}, @{term "number_of :: int => real"},
  1636    @{term "number_of :: int => nat"},
  1637    @{term "Int.Bit0"}, @{term "Int.Bit1"}, 
  1638    @{term "Int.Pls"}, @{term "Int.Min"}];
  1639 
  1640 fun check_sos kcts ct = 
  1641  let
  1642   val t = term_of ct
  1643   val _ = if not (null (Term.add_tfrees t []) 
  1644                   andalso null (Term.add_tvars t [])) 
  1645           then error "SOS: not sos. Additional type varables" else ()
  1646   val fs = Term.add_frees t []
  1647   val _ = if exists (fn ((_,T)) => not (T = @{typ "real"})) fs 
  1648           then error "SOS: not sos. Variables with type not real" else ()
  1649   val vs = Term.add_vars t []
  1650   val _ = if exists (fn ((_,T)) => not (T = @{typ "real"})) fs 
  1651           then error "SOS: not sos. Variables with type not real" else ()
  1652   val ukcs = subtract (fn (t,p) => Const p aconv t) kcts (Term.add_consts t [])
  1653   val _ = if  null ukcs then () 
  1654               else error ("SOSO: Unknown constants in Subgoal:" ^ commas (map fst ukcs))
  1655 in () end
  1656 
  1657 fun core_sos_tac ctxt = CSUBGOAL (fn (ct, i) => 
  1658   let val _ = check_sos known_sos_constants ct
  1659       val (avs, p) = strip_all ct
  1660       val th = standard (fold_rev forall_intr avs (real_sos ctxt (Thm.dest_arg p)))
  1661   in rtac th i end);
  1662 
  1663 fun default_SOME f NONE v = SOME v
  1664   | default_SOME f (SOME v) _ = SOME v;
  1665 
  1666 fun lift_SOME f NONE a = f a
  1667   | lift_SOME f (SOME a) _ = SOME a;
  1668 
  1669 
  1670 local
  1671  val is_numeral = can (HOLogic.dest_number o term_of)
  1672 in
  1673 fun get_denom b ct = case term_of ct of
  1674   @{term "op / :: real => _"} $ _ $ _ => 
  1675      if is_numeral (Thm.dest_arg ct) then get_denom b (Thm.dest_arg1 ct)
  1676      else default_SOME (get_denom b) (get_denom b (Thm.dest_arg ct))   (Thm.dest_arg ct, b)
  1677  | @{term "op < :: real => _"} $ _ $ _ => lift_SOME (get_denom true) (get_denom true (Thm.dest_arg ct)) (Thm.dest_arg1 ct)
  1678  | @{term "op <= :: real => _"} $ _ $ _ => lift_SOME (get_denom true) (get_denom true (Thm.dest_arg ct)) (Thm.dest_arg1 ct)
  1679  | _ $ _ => lift_SOME (get_denom b) (get_denom b (Thm.dest_fun ct)) (Thm.dest_arg ct)
  1680  | _ => NONE
  1681 end;
  1682 
  1683 fun elim_one_denom_tac ctxt = 
  1684 CSUBGOAL (fn (P,i) => 
  1685  case get_denom false P of 
  1686    NONE => no_tac
  1687  | SOME (d,ord) => 
  1688      let 
  1689       val ss = simpset_of (ProofContext.theory_of ctxt) addsimps @{thms field_simps} 
  1690                addsimps [@{thm nonzero_power_divide}, @{thm power_divide}]
  1691       val th = instantiate' [] [SOME d, SOME (Thm.dest_arg P)] 
  1692          (if ord then @{lemma "(d=0 --> P) & (d>0 --> P) & (d<(0::real) --> P) ==> P" by auto}
  1693           else @{lemma "(d=0 --> P) & (d ~= (0::real) --> P) ==> P" by blast})
  1694      in (rtac th i THEN Simplifier.asm_full_simp_tac ss i) end);
  1695 
  1696 fun elim_denom_tac ctxt i = REPEAT (elim_one_denom_tac ctxt i);
  1697 
  1698 fun sos_tac ctxt = ObjectLogic.full_atomize_tac THEN' elim_denom_tac ctxt THEN' core_sos_tac ctxt
  1699 
  1700 
  1701 end;