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