src/HOL/Library/Sum_of_Squares/sum_of_squares.ML
author Christian Sternagel
Thu Dec 13 13:11:38 2012 +0100 (2012-12-13)
changeset 50516 ed6b40d15d1c
parent 47108 2a1953f0d20d
child 51717 9e7d1c139569
permissions -rw-r--r--
renamed "emb" to "list_hembeq";
make "list_hembeq" reflexive independent of the base order;
renamed "sub" to "sublisteq";
dropped "transp_on" (state transitivity explicitly instead);
no need to hide "sub" after renaming;
replaced some ASCII symbols by proper Isabelle symbols;
NEWS
     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 trace: bool Config.T
    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 max = Integer.max;
    24 
    25 val denominator_rat = Rat.quotient_of_rat #> snd #> Rat.rat_of_int;
    26 fun int_of_rat a =
    27     case Rat.quotient_of_rat a of (i,1) => i | _ => error "int_of_rat: not an int";
    28 fun lcm_rat x y = Rat.rat_of_int (Integer.lcm (int_of_rat x) (int_of_rat y));
    29 
    30 fun rat_pow r i =
    31  let fun pow r i =
    32    if i = 0 then rat_1 else
    33    let val d = pow r (i div 2)
    34    in d */ d */ (if i mod 2 = 0 then rat_1 else r)
    35    end
    36  in if i < 0 then pow (Rat.inv r) (~ i) else pow r i end;
    37 
    38 fun round_rat r =
    39  let val (a,b) = Rat.quotient_of_rat (Rat.abs r)
    40      val d = a div b
    41      val s = if r </ rat_0 then (Rat.neg o Rat.rat_of_int) else Rat.rat_of_int
    42      val x2 = 2 * (a - (b * d))
    43  in s (if x2 >= b then d + 1 else d) end
    44 
    45 val abs_rat = Rat.abs;
    46 val pow2 = rat_pow rat_2;
    47 val pow10 = rat_pow rat_10;
    48 
    49 val trace = Attrib.setup_config_bool @{binding sos_trace} (K false);
    50 
    51 exception Sanity;
    52 
    53 exception Unsolvable;
    54 
    55 exception Failure of string;
    56 
    57 datatype proof_method =
    58     Certificate of RealArith.pss_tree
    59   | Prover of (string -> string)
    60 
    61 (* Turn a rational into a decimal string with d sig digits.                  *)
    62 
    63 local
    64 fun normalize y =
    65   if abs_rat y </ (rat_1 // rat_10) then normalize (rat_10 */ y) - 1
    66   else if abs_rat y >=/ rat_1 then normalize (y // rat_10) + 1
    67   else 0
    68  in
    69 fun decimalize d x =
    70   if x =/ rat_0 then "0.0" else
    71   let
    72    val y = Rat.abs x
    73    val e = normalize y
    74    val z = pow10(~ e) */ y +/ rat_1
    75    val k = int_of_rat (round_rat(pow10 d */ z))
    76   in (if x </ rat_0 then "-0." else "0.") ^
    77      implode(tl(raw_explode(string_of_int k))) ^
    78      (if e = 0 then "" else "e"^string_of_int e)
    79   end
    80 end;
    81 
    82 (* Iterations over numbers, and lists indexed by numbers.                    *)
    83 
    84 fun itern k l f a =
    85   case l of
    86     [] => a
    87   | h::t => itern (k + 1) t f (f h k a);
    88 
    89 fun iter (m,n) f a =
    90   if n < m then a
    91   else iter (m+1,n) f (f m a);
    92 
    93 (* The main types.                                                           *)
    94 
    95 type vector = int* Rat.rat FuncUtil.Intfunc.table;
    96 
    97 type matrix = (int*int)*(Rat.rat FuncUtil.Intpairfunc.table);
    98 
    99 fun iszero (_,r) = r =/ rat_0;
   100 
   101 
   102 (* Vectors. Conventionally indexed 1..n.                                     *)
   103 
   104 fun vector_0 n = (n,FuncUtil.Intfunc.empty):vector;
   105 
   106 fun dim (v:vector) = fst v;
   107 
   108 fun vector_cmul c (v:vector) =
   109  let val n = dim v
   110  in if c =/ rat_0 then vector_0 n
   111     else (n,FuncUtil.Intfunc.map (fn _ => fn x => c */ x) (snd v))
   112  end;
   113 
   114 fun vector_of_list l =
   115  let val n = length l
   116  in (n,fold_rev2 (curry FuncUtil.Intfunc.update) (1 upto n) l FuncUtil.Intfunc.empty) :vector
   117  end;
   118 
   119 (* Matrices; again rows and columns indexed from 1.                          *)
   120 
   121 fun dimensions (m:matrix) = fst m;
   122 
   123 fun row k (m:matrix) =
   124  let val (_,j) = dimensions m
   125  in (j,
   126    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
   127  end;
   128 
   129 (* Monomials.                                                                *)
   130 
   131 fun monomial_eval assig m =
   132   FuncUtil.Ctermfunc.fold (fn (x, k) => fn a => a */ rat_pow (FuncUtil.Ctermfunc.apply assig x) k)
   133         m rat_1;
   134 val monomial_1 = FuncUtil.Ctermfunc.empty;
   135 
   136 fun monomial_var x = FuncUtil.Ctermfunc.onefunc (x, 1);
   137 
   138 val monomial_mul =
   139   FuncUtil.Ctermfunc.combine Integer.add (K false);
   140 
   141 fun monomial_multidegree m =
   142  FuncUtil.Ctermfunc.fold (fn (_, k) => fn a => k + a) m 0;;
   143 
   144 fun monomial_variables m = FuncUtil.Ctermfunc.dom m;;
   145 
   146 (* Polynomials.                                                              *)
   147 
   148 fun eval assig p =
   149   FuncUtil.Monomialfunc.fold (fn (m, c) => fn a => a +/ c */ monomial_eval assig m) p rat_0;
   150 
   151 val poly_0 = FuncUtil.Monomialfunc.empty;
   152 
   153 fun poly_isconst p =
   154   FuncUtil.Monomialfunc.fold (fn (m, _) => fn a => FuncUtil.Ctermfunc.is_empty m andalso a) p true;
   155 
   156 fun poly_var x = FuncUtil.Monomialfunc.onefunc (monomial_var x,rat_1);
   157 
   158 fun poly_const c =
   159   if c =/ rat_0 then poly_0 else FuncUtil.Monomialfunc.onefunc(monomial_1, c);
   160 
   161 fun poly_cmul c p =
   162   if c =/ rat_0 then poly_0
   163   else FuncUtil.Monomialfunc.map (fn _ => fn x => c */ x) p;
   164 
   165 fun poly_neg p = FuncUtil.Monomialfunc.map (K Rat.neg) p;;
   166 
   167 fun poly_add p1 p2 =
   168   FuncUtil.Monomialfunc.combine (curry op +/) (fn x => x =/ rat_0) p1 p2;
   169 
   170 fun poly_sub p1 p2 = poly_add p1 (poly_neg p2);
   171 
   172 fun poly_cmmul (c,m) p =
   173  if c =/ rat_0 then poly_0
   174  else if FuncUtil.Ctermfunc.is_empty m
   175       then FuncUtil.Monomialfunc.map (fn _ => fn d => c */ d) p
   176       else FuncUtil.Monomialfunc.fold (fn (m', d) => fn a => (FuncUtil.Monomialfunc.update (monomial_mul m m', c */ d) a)) p poly_0;
   177 
   178 fun poly_mul p1 p2 =
   179   FuncUtil.Monomialfunc.fold (fn (m, c) => fn a => poly_add (poly_cmmul (c,m) p2) a) p1 poly_0;
   180 
   181 fun poly_square p = poly_mul p p;
   182 
   183 fun poly_pow p k =
   184  if k = 0 then poly_const rat_1
   185  else if k = 1 then p
   186  else let val q = poly_square(poly_pow p (k div 2)) in
   187       if k mod 2 = 1 then poly_mul p q else q end;
   188 
   189 fun multidegree p =
   190   FuncUtil.Monomialfunc.fold (fn (m, _) => fn a => max (monomial_multidegree m) a) p 0;
   191 
   192 fun poly_variables p =
   193   sort FuncUtil.cterm_ord (FuncUtil.Monomialfunc.fold_rev (fn (m, _) => union (is_equal o FuncUtil.cterm_ord) (monomial_variables m)) p []);;
   194 
   195 (* Conversion from HOL term.                                                 *)
   196 
   197 local
   198  val neg_tm = @{cterm "uminus :: real => _"}
   199  val add_tm = @{cterm "op + :: real => _"}
   200  val sub_tm = @{cterm "op - :: real => _"}
   201  val mul_tm = @{cterm "op * :: real => _"}
   202  val inv_tm = @{cterm "inverse :: real => _"}
   203  val div_tm = @{cterm "op / :: real => _"}
   204  val pow_tm = @{cterm "op ^ :: real => _"}
   205  val zero_tm = @{cterm "0:: real"}
   206  val is_numeral = can (HOLogic.dest_number o term_of)
   207  fun poly_of_term tm =
   208   if tm aconvc zero_tm then poly_0
   209   else if RealArith.is_ratconst tm
   210        then poly_const(RealArith.dest_ratconst tm)
   211   else
   212   (let val (lop,r) = Thm.dest_comb tm
   213    in if lop aconvc neg_tm then poly_neg(poly_of_term r)
   214       else if lop aconvc inv_tm then
   215        let val p = poly_of_term r
   216        in if poly_isconst p
   217           then poly_const(Rat.inv (eval FuncUtil.Ctermfunc.empty p))
   218           else error "poly_of_term: inverse of non-constant polyomial"
   219        end
   220    else (let val (opr,l) = Thm.dest_comb lop
   221          in
   222           if opr aconvc pow_tm andalso is_numeral r
   223           then poly_pow (poly_of_term l) ((snd o HOLogic.dest_number o term_of) r)
   224           else if opr aconvc add_tm
   225            then poly_add (poly_of_term l) (poly_of_term r)
   226           else if opr aconvc sub_tm
   227            then poly_sub (poly_of_term l) (poly_of_term r)
   228           else if opr aconvc mul_tm
   229            then poly_mul (poly_of_term l) (poly_of_term r)
   230           else if opr aconvc div_tm
   231            then let
   232                   val p = poly_of_term l
   233                   val q = poly_of_term r
   234                 in if poly_isconst q then poly_cmul (Rat.inv (eval FuncUtil.Ctermfunc.empty q)) p
   235                    else error "poly_of_term: division by non-constant polynomial"
   236                 end
   237           else poly_var tm
   238 
   239          end
   240          handle CTERM ("dest_comb",_) => poly_var tm)
   241    end
   242    handle CTERM ("dest_comb",_) => poly_var tm)
   243 in
   244 val poly_of_term = fn tm =>
   245  if type_of (term_of tm) = @{typ real} then poly_of_term tm
   246  else error "poly_of_term: term does not have real type"
   247 end;
   248 
   249 (* String of vector (just a list of space-separated numbers).                *)
   250 
   251 fun sdpa_of_vector (v:vector) =
   252  let
   253   val n = dim v
   254   val strs = map (decimalize 20 o (fn i => FuncUtil.Intfunc.tryapplyd (snd v) i rat_0)) (1 upto n)
   255  in space_implode " " strs ^ "\n"
   256  end;
   257 
   258 fun triple_int_ord ((a,b,c),(a',b',c')) =
   259  prod_ord int_ord (prod_ord int_ord int_ord)
   260     ((a,(b,c)),(a',(b',c')));
   261 structure Inttriplefunc = FuncFun(type key = int*int*int val ord = triple_int_ord);
   262 
   263 fun index_char str chr pos =
   264   if pos >= String.size str then ~1
   265   else if String.sub(str,pos) = chr then pos
   266   else index_char str chr (pos + 1);
   267 fun rat_of_quotient (a,b) = if b = 0 then rat_0 else Rat.rat_of_quotient (a,b);
   268 fun rat_of_string s =
   269  let val n = index_char s #"/" 0 in
   270   if n = ~1 then s |> Int.fromString |> the |> Rat.rat_of_int
   271   else
   272    let val SOME numer = Int.fromString(String.substring(s,0,n))
   273        val SOME den = Int.fromString (String.substring(s,n+1,String.size s - n - 1))
   274    in rat_of_quotient(numer, den)
   275    end
   276  end;
   277 
   278 fun isnum x = member (op =) ["0","1","2","3","4","5","6","7","8","9"] x;
   279 
   280 (* More parser basics.                                                       *)
   281 
   282  val numeral = Scan.one isnum
   283  val decimalint = Scan.repeat1 numeral >> (rat_of_string o implode)
   284  val decimalfrac = Scan.repeat1 numeral
   285     >> (fn s => rat_of_string(implode s) // pow10 (length s))
   286  val decimalsig =
   287     decimalint -- Scan.option (Scan.$$ "." |-- decimalfrac)
   288     >> (fn (h,NONE) => h | (h,SOME x) => h +/ x)
   289  fun signed prs =
   290        $$ "-" |-- prs >> Rat.neg
   291     || $$ "+" |-- prs
   292     || prs;
   293 
   294 fun emptyin def xs = if null xs then (def,xs) else Scan.fail xs
   295 
   296  val exponent = ($$ "e" || $$ "E") |-- signed decimalint;
   297 
   298  val decimal = signed decimalsig -- (emptyin rat_0|| exponent)
   299     >> (fn (h, x) => h */ pow10 (int_of_rat x));
   300 
   301  fun mkparser p s =
   302   let val (x,rst) = p (raw_explode s)
   303   in if null rst then x
   304      else error "mkparser: unparsed input"
   305   end;;
   306 
   307 (* Parse back csdp output.                                                      *)
   308 
   309  fun ignore _ = ((),[])
   310  fun csdpoutput inp =
   311    ((decimal -- Scan.repeat (Scan.$$ " " |-- Scan.option decimal) >>
   312     (fn (h,to) => map_filter I ((SOME h)::to))) --| ignore >> vector_of_list) inp
   313  val parse_csdpoutput = mkparser csdpoutput
   314 
   315 (* Try some apparently sensible scaling first. Note that this is purely to   *)
   316 (* get a cleaner translation to floating-point, and doesn't affect any of    *)
   317 (* the results, in principle. In practice it seems a lot better when there   *)
   318 (* are extreme numbers in the original problem.                              *)
   319 
   320   (* Version for (int*int*int) keys *)
   321 local
   322   fun max_rat x y = if x </ y then y else x
   323   fun common_denominator fld amat acc =
   324       fld (fn (_,c) => fn a => lcm_rat (denominator_rat c) a) amat acc
   325   fun maximal_element fld amat acc =
   326     fld (fn (_,c) => fn maxa => max_rat maxa (abs_rat c)) amat acc
   327 fun float_of_rat x = let val (a,b) = Rat.quotient_of_rat x
   328                      in Real.fromInt a / Real.fromInt b end;
   329 fun int_of_float x = (trunc x handle Overflow => 0 | Domain => 0)
   330 in
   331 
   332 fun tri_scale_then solver (obj:vector)  mats =
   333  let
   334   val cd1 = fold_rev (common_denominator Inttriplefunc.fold) mats (rat_1)
   335   val cd2 = common_denominator FuncUtil.Intfunc.fold (snd obj)  (rat_1)
   336   val mats' = map (Inttriplefunc.map (fn _ => fn x => cd1 */ x)) mats
   337   val obj' = vector_cmul cd2 obj
   338   val max1 = fold_rev (maximal_element Inttriplefunc.fold) mats' (rat_0)
   339   val max2 = maximal_element FuncUtil.Intfunc.fold (snd obj') (rat_0)
   340   val scal1 = pow2 (20 - int_of_float(Math.ln (float_of_rat max1) / Math.ln 2.0))
   341   val scal2 = pow2 (20 - int_of_float(Math.ln (float_of_rat max2) / Math.ln 2.0))
   342   val mats'' = map (Inttriplefunc.map (fn _ => fn x => x */ scal1)) mats'
   343   val obj'' = vector_cmul scal2 obj'
   344  in solver obj'' mats''
   345   end
   346 end;
   347 
   348 (* Round a vector to "nice" rationals.                                       *)
   349 
   350 fun nice_rational n x = round_rat (n */ x) // n;;
   351 fun nice_vector n ((d,v) : vector) =
   352  (d, FuncUtil.Intfunc.fold (fn (i,c) => fn a =>
   353    let val y = nice_rational n c
   354    in if c =/ rat_0 then a
   355       else FuncUtil.Intfunc.update (i,y) a end) v FuncUtil.Intfunc.empty):vector
   356 
   357 fun dest_ord f x = is_equal (f x);
   358 
   359 (* Stuff for "equations" ((int*int*int)->num functions).                         *)
   360 
   361 fun tri_equation_cmul c eq =
   362   if c =/ rat_0 then Inttriplefunc.empty else Inttriplefunc.map (fn _ => fn d => c */ d) eq;
   363 
   364 fun tri_equation_add eq1 eq2 = Inttriplefunc.combine (curry op +/) (fn x => x =/ rat_0) eq1 eq2;
   365 
   366 fun tri_equation_eval assig eq =
   367  let fun value v = Inttriplefunc.apply assig v
   368  in Inttriplefunc.fold (fn (v, c) => fn a => a +/ value v */ c) eq rat_0
   369  end;
   370 
   371 (* Eliminate all variables, in an essentially arbitrary order.               *)
   372 
   373 fun tri_eliminate_all_equations one =
   374  let
   375   fun choose_variable eq =
   376    let val (v,_) = Inttriplefunc.choose eq
   377    in if is_equal (triple_int_ord(v,one)) then
   378       let val eq' = Inttriplefunc.delete_safe v eq
   379       in if Inttriplefunc.is_empty eq' then error "choose_variable"
   380          else fst (Inttriplefunc.choose eq')
   381       end
   382     else v
   383    end
   384   fun eliminate dun eqs = case eqs of
   385     [] => dun
   386   | eq::oeqs =>
   387     if Inttriplefunc.is_empty eq then eliminate dun oeqs else
   388     let val v = choose_variable eq
   389         val a = Inttriplefunc.apply eq v
   390         val eq' = tri_equation_cmul ((Rat.rat_of_int ~1) // a)
   391                    (Inttriplefunc.delete_safe v eq)
   392         fun elim e =
   393          let val b = Inttriplefunc.tryapplyd e v rat_0
   394          in if b =/ rat_0 then e
   395             else tri_equation_add e (tri_equation_cmul (Rat.neg b // a) eq)
   396          end
   397     in eliminate (Inttriplefunc.update(v, eq') (Inttriplefunc.map (K elim) dun))
   398                  (map elim oeqs)
   399     end
   400 in fn eqs =>
   401  let
   402   val assig = eliminate Inttriplefunc.empty eqs
   403   val vs = Inttriplefunc.fold (fn (_, f) => fn a => remove (dest_ord triple_int_ord) one (Inttriplefunc.dom f) @ a) assig []
   404  in (distinct (dest_ord triple_int_ord) vs,assig)
   405  end
   406 end;
   407 
   408 (* Multiply equation-parametrized poly by regular poly and add accumulator.  *)
   409 
   410 fun tri_epoly_pmul p q acc =
   411  FuncUtil.Monomialfunc.fold (fn (m1, c) => fn a =>
   412   FuncUtil.Monomialfunc.fold (fn (m2,e) => fn b =>
   413    let val m =  monomial_mul m1 m2
   414        val es = FuncUtil.Monomialfunc.tryapplyd b m Inttriplefunc.empty
   415    in FuncUtil.Monomialfunc.update (m,tri_equation_add (tri_equation_cmul c e) es) b
   416    end) q a) p acc ;
   417 
   418 (* Hence produce the "relevant" monomials: those whose squares lie in the    *)
   419 (* Newton polytope of the monomials in the input. (This is enough according  *)
   420 (* to Reznik: "Extremal PSD forms with few terms", Duke Math. Journal,       *)
   421 (* vol 45, pp. 363--374, 1978.                                               *)
   422 (*                                                                           *)
   423 (* These are ordered in sort of decreasing degree. In particular the         *)
   424 (* constant monomial is last; this gives an order in diagonalization of the  *)
   425 (* quadratic form that will tend to display constants.                       *)
   426 
   427 (* Diagonalize (Cholesky/LDU) the matrix corresponding to a quadratic form.  *)
   428 
   429 local
   430 fun diagonalize n i m =
   431  if FuncUtil.Intpairfunc.is_empty (snd m) then []
   432  else
   433   let val a11 = FuncUtil.Intpairfunc.tryapplyd (snd m) (i,i) rat_0
   434   in if a11 </ rat_0 then raise Failure "diagonalize: not PSD"
   435     else if a11 =/ rat_0 then
   436           if FuncUtil.Intfunc.is_empty (snd (row i m)) then diagonalize n (i + 1) m
   437           else raise Failure "diagonalize: not PSD ___ "
   438     else
   439      let
   440       val v = row i m
   441       val v' = (fst v, FuncUtil.Intfunc.fold (fn (i, c) => fn a =>
   442        let val y = c // a11
   443        in if y = rat_0 then a else FuncUtil.Intfunc.update (i,y) a
   444        end)  (snd v) FuncUtil.Intfunc.empty)
   445       fun upt0 x y a = if y = rat_0 then a else FuncUtil.Intpairfunc.update (x,y) a
   446       val m' =
   447       ((n,n),
   448       iter (i+1,n) (fn j =>
   449           iter (i+1,n) (fn k =>
   450               (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))))
   451           FuncUtil.Intpairfunc.empty)
   452      in (a11,v')::diagonalize n (i + 1) m'
   453      end
   454   end
   455 in
   456 fun diag m =
   457  let
   458    val nn = dimensions m
   459    val n = fst nn
   460  in if snd nn <> n then error "diagonalize: non-square matrix"
   461     else diagonalize n 1 m
   462  end
   463 end;
   464 
   465 (* Enumeration of monomials with given multidegree bound.                    *)
   466 
   467 fun enumerate_monomials d vars =
   468  if d < 0 then []
   469  else if d = 0 then [FuncUtil.Ctermfunc.empty]
   470  else if null vars then [monomial_1] else
   471  let val alts =
   472   map_range (fn k => let val oths = enumerate_monomials (d - k) (tl vars)
   473                in map (fn ks => if k = 0 then ks else FuncUtil.Ctermfunc.update (hd vars, k) ks) oths end) (d + 1)
   474  in flat alts
   475  end;
   476 
   477 (* Enumerate products of distinct input polys with degree <= d.              *)
   478 (* We ignore any constant input polynomials.                                 *)
   479 (* Give the output polynomial and a record of how it was derived.            *)
   480 
   481 fun enumerate_products d pols =
   482 if d = 0 then [(poly_const rat_1,RealArith.Rational_lt rat_1)]
   483 else if d < 0 then [] else
   484 case pols of
   485    [] => [(poly_const rat_1,RealArith.Rational_lt rat_1)]
   486  | (p,b)::ps =>
   487     let val e = multidegree p
   488     in if e = 0 then enumerate_products d ps else
   489        enumerate_products d ps @
   490        map (fn (q,c) => (poly_mul p q,RealArith.Product(b,c)))
   491          (enumerate_products (d - e) ps)
   492     end
   493 
   494 (* Convert regular polynomial. Note that we treat (0,0,0) as -1.             *)
   495 
   496 fun epoly_of_poly p =
   497   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;
   498 
   499 (* String for block diagonal matrix numbered k.                              *)
   500 
   501 fun sdpa_of_blockdiagonal k m =
   502  let
   503   val pfx = string_of_int k ^" "
   504   val ents =
   505     Inttriplefunc.fold
   506       (fn ((b,i,j),c) => fn a => if i > j then a else ((b,i,j),c)::a)
   507       m []
   508   val entss = sort (triple_int_ord o pairself fst) ents
   509  in fold_rev (fn ((b,i,j),c) => fn a =>
   510      pfx ^ string_of_int b ^ " " ^ string_of_int i ^ " " ^ string_of_int j ^
   511      " " ^ decimalize 20 c ^ "\n" ^ a) entss ""
   512  end;
   513 
   514 (* SDPA for problem using block diagonal (i.e. multiple SDPs)                *)
   515 
   516 fun sdpa_of_blockproblem nblocks blocksizes obj mats =
   517  let val m = length mats - 1
   518  in
   519   string_of_int m ^ "\n" ^
   520   string_of_int nblocks ^ "\n" ^
   521   (space_implode " " (map string_of_int blocksizes)) ^
   522   "\n" ^
   523   sdpa_of_vector obj ^
   524   fold_rev2 (fn k => fn m => fn a => sdpa_of_blockdiagonal (k - 1) m ^ a)
   525     (1 upto length mats) mats ""
   526  end;
   527 
   528 (* Run prover on a problem in block diagonal form.                       *)
   529 
   530 fun run_blockproblem prover nblocks blocksizes obj mats=
   531   parse_csdpoutput (prover (sdpa_of_blockproblem nblocks blocksizes obj mats))
   532 
   533 (* 3D versions of matrix operations to consider blocks separately.           *)
   534 
   535 val bmatrix_add = Inttriplefunc.combine (curry op +/) (fn x => x =/ rat_0);
   536 fun bmatrix_cmul c bm =
   537   if c =/ rat_0 then Inttriplefunc.empty
   538   else Inttriplefunc.map (fn _ => fn x => c */ x) bm;
   539 
   540 val bmatrix_neg = bmatrix_cmul (Rat.rat_of_int ~1);
   541 
   542 (* Smash a block matrix into components.                                     *)
   543 
   544 fun blocks blocksizes bm =
   545  map (fn (bs,b0) =>
   546       let val m = Inttriplefunc.fold
   547           (fn ((b,i,j),c) => fn a => if b = b0 then FuncUtil.Intpairfunc.update ((i,j),c) a else a) bm FuncUtil.Intpairfunc.empty
   548           val _ = FuncUtil.Intpairfunc.fold (fn ((i,j),_) => fn a => max a (max i j)) m 0
   549       in (((bs,bs),m):matrix) end)
   550  (blocksizes ~~ (1 upto length blocksizes));;
   551 
   552 (* FIXME : Get rid of this !!!*)
   553 local
   554   fun tryfind_with msg _ [] = raise Failure msg
   555     | tryfind_with _ f (x::xs) = (f x handle Failure s => tryfind_with s f xs);
   556 in
   557   fun tryfind f = tryfind_with "tryfind" f
   558 end
   559 
   560 (* Positiv- and Nullstellensatz. Flag "linf" forces a linear representation. *)
   561 
   562 
   563 fun real_positivnullstellensatz_general ctxt prover linf d eqs leqs pol =
   564 let
   565  val vars = fold_rev (union (op aconvc) o poly_variables)
   566    (pol :: eqs @ map fst leqs) []
   567  val monoid = if linf then
   568       (poly_const rat_1,RealArith.Rational_lt rat_1)::
   569       (filter (fn (p,_) => multidegree p <= d) leqs)
   570     else enumerate_products d leqs
   571  val nblocks = length monoid
   572  fun mk_idmultiplier k p =
   573   let
   574    val e = d - multidegree p
   575    val mons = enumerate_monomials e vars
   576    val nons = mons ~~ (1 upto length mons)
   577   in (mons,
   578       fold_rev (fn (m,n) => FuncUtil.Monomialfunc.update(m,Inttriplefunc.onefunc((~k,~n,n),rat_1))) nons FuncUtil.Monomialfunc.empty)
   579   end
   580 
   581  fun mk_sqmultiplier k (p,_) =
   582   let
   583    val e = (d - multidegree p) div 2
   584    val mons = enumerate_monomials e vars
   585    val nons = mons ~~ (1 upto length mons)
   586   in (mons,
   587       fold_rev (fn (m1,n1) =>
   588        fold_rev (fn (m2,n2) => fn  a =>
   589         let val m = monomial_mul m1 m2
   590         in if n1 > n2 then a else
   591           let val c = if n1 = n2 then rat_1 else rat_2
   592               val e = FuncUtil.Monomialfunc.tryapplyd a m Inttriplefunc.empty
   593           in FuncUtil.Monomialfunc.update(m, tri_equation_add (Inttriplefunc.onefunc((k,n1,n2), c)) e) a
   594           end
   595         end)  nons)
   596        nons FuncUtil.Monomialfunc.empty)
   597   end
   598 
   599   val (sqmonlist,sqs) = split_list (map2 mk_sqmultiplier (1 upto length monoid) monoid)
   600   val (_(*idmonlist*),ids) =  split_list(map2 mk_idmultiplier (1 upto length eqs) eqs)
   601   val blocksizes = map length sqmonlist
   602   val bigsum =
   603     fold_rev2 (fn p => fn q => fn a => tri_epoly_pmul p q a) eqs ids
   604             (fold_rev2 (fn (p,_) => fn s => fn a => tri_epoly_pmul p s a) monoid sqs
   605                      (epoly_of_poly(poly_neg pol)))
   606   val eqns = FuncUtil.Monomialfunc.fold (fn (_,e) => fn a => e::a) bigsum []
   607   val (pvs,assig) = tri_eliminate_all_equations (0,0,0) eqns
   608   val qvars = (0,0,0)::pvs
   609   val allassig = fold_rev (fn v => Inttriplefunc.update(v,(Inttriplefunc.onefunc(v,rat_1)))) pvs assig
   610   fun mk_matrix v =
   611     Inttriplefunc.fold (fn ((b,i,j), ass) => fn m =>
   612         if b < 0 then m else
   613          let val c = Inttriplefunc.tryapplyd ass v rat_0
   614          in if c = rat_0 then m else
   615             Inttriplefunc.update ((b,j,i), c) (Inttriplefunc.update ((b,i,j), c) m)
   616          end)
   617           allassig Inttriplefunc.empty
   618   val diagents = Inttriplefunc.fold
   619     (fn ((b,i,j), e) => fn a => if b > 0 andalso i = j then tri_equation_add e a else a)
   620     allassig Inttriplefunc.empty
   621 
   622   val mats = map mk_matrix qvars
   623   val obj = (length pvs,
   624             itern 1 pvs (fn v => fn i => FuncUtil.Intfunc.updatep iszero (i,Inttriplefunc.tryapplyd diagents v rat_0))
   625                         FuncUtil.Intfunc.empty)
   626   val raw_vec = if null pvs then vector_0 0
   627                 else tri_scale_then (run_blockproblem prover nblocks blocksizes) obj mats
   628   fun int_element (_,v) i = FuncUtil.Intfunc.tryapplyd v i rat_0
   629 
   630   fun find_rounding d =
   631    let
   632     val _ =
   633       if Config.get ctxt trace
   634       then writeln ("Trying rounding with limit "^Rat.string_of_rat d ^ "\n")
   635       else ()
   636     val vec = nice_vector d raw_vec
   637     val blockmat = iter (1,dim vec)
   638      (fn i => fn a => bmatrix_add (bmatrix_cmul (int_element vec i) (nth mats i)) a)
   639      (bmatrix_neg (nth mats 0))
   640     val allmats = blocks blocksizes blockmat
   641    in (vec,map diag allmats)
   642    end
   643   val (vec,ratdias) =
   644     if null pvs then find_rounding rat_1
   645     else tryfind find_rounding (map Rat.rat_of_int (1 upto 31) @
   646                                 map pow2 (5 upto 66))
   647   val newassigs =
   648     fold_rev (fn k => Inttriplefunc.update (nth pvs (k - 1), int_element vec k))
   649            (1 upto dim vec) (Inttriplefunc.onefunc ((0,0,0), Rat.rat_of_int ~1))
   650   val finalassigs =
   651     Inttriplefunc.fold (fn (v,e) => fn a => Inttriplefunc.update(v, tri_equation_eval newassigs e) a) allassig newassigs
   652   fun poly_of_epoly p =
   653     FuncUtil.Monomialfunc.fold (fn (v,e) => fn a => FuncUtil.Monomialfunc.updatep iszero (v,tri_equation_eval finalassigs e) a)
   654           p FuncUtil.Monomialfunc.empty
   655   fun  mk_sos mons =
   656    let fun mk_sq (c,m) =
   657     (c,fold_rev (fn k=> fn a => FuncUtil.Monomialfunc.updatep iszero (nth mons (k - 1), int_element m k) a)
   658                  (1 upto length mons) FuncUtil.Monomialfunc.empty)
   659    in map mk_sq
   660    end
   661   val sqs = map2 mk_sos sqmonlist ratdias
   662   val cfs = map poly_of_epoly ids
   663   val msq = filter (fn (_,b) => not (null b)) (map2 pair monoid sqs)
   664   fun eval_sq sqs = fold_rev (fn (c,q) => poly_add (poly_cmul c (poly_mul q q))) sqs poly_0
   665   val sanity =
   666     fold_rev (fn ((p,_),s) => poly_add (poly_mul p (eval_sq s))) msq
   667            (fold_rev2 (fn p => fn q => poly_add (poly_mul p q)) cfs eqs
   668                     (poly_neg pol))
   669 
   670 in if not(FuncUtil.Monomialfunc.is_empty sanity) then raise Sanity else
   671   (cfs,map (fn (a,b) => (snd a,b)) msq)
   672  end
   673 
   674 
   675 (* Iterative deepening.                                                      *)
   676 
   677 fun deepen f n =
   678   (writeln ("Searching with depth limit " ^ string_of_int n);
   679     (f n handle Failure s => (writeln ("failed with message: " ^ s); deepen f (n + 1))));
   680 
   681 
   682 (* Map back polynomials and their composites to a positivstellensatz.        *)
   683 
   684 fun cterm_of_sqterm (c,p) = RealArith.Product(RealArith.Rational_lt c,RealArith.Square p);
   685 
   686 fun cterm_of_sos (pr,sqs) = if null sqs then pr
   687   else RealArith.Product(pr,foldr1 RealArith.Sum (map cterm_of_sqterm sqs));
   688 
   689 (* Interface to HOL.                                                         *)
   690 local
   691   open Conv
   692   val concl = Thm.dest_arg o cprop_of
   693   fun simple_cterm_ord t u = Term_Ord.fast_term_ord (term_of t, term_of u) = LESS
   694 in
   695   (* FIXME: Replace tryfind by get_first !! *)
   696 fun real_nonlinear_prover proof_method ctxt =
   697  let
   698   val {add = _, mul = _, neg = _, pow = _,
   699        sub = _, main = real_poly_conv} =
   700       Semiring_Normalizer.semiring_normalizers_ord_wrapper ctxt
   701       (the (Semiring_Normalizer.match ctxt @{cterm "(0::real) + 1"}))
   702      simple_cterm_ord
   703   fun mainf cert_choice translator (eqs,les,lts) =
   704   let
   705    val eq0 = map (poly_of_term o Thm.dest_arg1 o concl) eqs
   706    val le0 = map (poly_of_term o Thm.dest_arg o concl) les
   707    val lt0 = map (poly_of_term o Thm.dest_arg o concl) lts
   708    val eqp0 = map_index (fn (i, t) => (t,RealArith.Axiom_eq i)) eq0
   709    val lep0 = map_index (fn (i, t) => (t,RealArith.Axiom_le i)) le0
   710    val ltp0 = map_index (fn (i, t) => (t,RealArith.Axiom_lt i)) lt0
   711    val (keq,eq) = List.partition (fn (p,_) => multidegree p = 0) eqp0
   712    val (klep,lep) = List.partition (fn (p,_) => multidegree p = 0) lep0
   713    val (kltp,ltp) = List.partition (fn (p,_) => multidegree p = 0) ltp0
   714    fun trivial_axiom (p,ax) =
   715     case ax of
   716        RealArith.Axiom_eq n => if eval FuncUtil.Ctermfunc.empty p <>/ Rat.zero then nth eqs n
   717                      else raise Failure "trivial_axiom: Not a trivial axiom"
   718      | RealArith.Axiom_le n => if eval FuncUtil.Ctermfunc.empty p </ Rat.zero then nth les n
   719                      else raise Failure "trivial_axiom: Not a trivial axiom"
   720      | RealArith.Axiom_lt n => if eval FuncUtil.Ctermfunc.empty p <=/ Rat.zero then nth lts n
   721                      else raise Failure "trivial_axiom: Not a trivial axiom"
   722      | _ => error "trivial_axiom: Not a trivial axiom"
   723    in
   724   (let val th = tryfind trivial_axiom (keq @ klep @ kltp)
   725    in
   726     (fconv_rule (arg_conv (arg1_conv real_poly_conv) then_conv Numeral_Simprocs.field_comp_conv) th, RealArith.Trivial)
   727    end)
   728    handle Failure _ =>
   729      (let val proof =
   730        (case proof_method of Certificate certs =>
   731          (* choose certificate *)
   732          let
   733            fun chose_cert [] (RealArith.Cert c) = c
   734              | chose_cert (RealArith.Left::s) (RealArith.Branch (l, _)) = chose_cert s l
   735              | chose_cert (RealArith.Right::s) (RealArith.Branch (_, r)) = chose_cert s r
   736              | chose_cert _ _ = error "certificate tree in invalid form"
   737          in
   738            chose_cert cert_choice certs
   739          end
   740        | Prover prover =>
   741          (* call prover *)
   742          let
   743           val pol = fold_rev poly_mul (map fst ltp) (poly_const Rat.one)
   744           val leq = lep @ ltp
   745           fun tryall d =
   746            let val e = multidegree pol
   747                val k = if e = 0 then 0 else d div e
   748                val eq' = map fst eq
   749            in tryfind (fn i => (d,i,real_positivnullstellensatz_general ctxt prover false d eq' leq
   750                                  (poly_neg(poly_pow pol i))))
   751                    (0 upto k)
   752            end
   753          val (_,i,(cert_ideal,cert_cone)) = deepen tryall 0
   754          val proofs_ideal =
   755            map2 (fn q => fn (_,ax) => RealArith.Eqmul(q,ax)) cert_ideal eq
   756          val proofs_cone = map cterm_of_sos cert_cone
   757          val proof_ne = if null ltp then RealArith.Rational_lt Rat.one else
   758            let val p = foldr1 RealArith.Product (map snd ltp)
   759            in  funpow i (fn q => RealArith.Product(p,q)) (RealArith.Rational_lt Rat.one)
   760            end
   761          in
   762            foldr1 RealArith.Sum (proof_ne :: proofs_ideal @ proofs_cone)
   763          end)
   764      in
   765         (translator (eqs,les,lts) proof, RealArith.Cert proof)
   766      end)
   767    end
   768  in mainf end
   769 end
   770 
   771 fun C f x y = f y x;
   772   (* FIXME : This is very bad!!!*)
   773 fun subst_conv eqs t =
   774  let
   775   val t' = fold (Thm.lambda o Thm.lhs_of) eqs t
   776  in Conv.fconv_rule (Thm.beta_conversion true) (fold (C Thm.combination) eqs (Thm.reflexive t'))
   777  end
   778 
   779 (* A wrapper that tries to substitute away variables first.                  *)
   780 
   781 local
   782  open Conv
   783   fun simple_cterm_ord t u = Term_Ord.fast_term_ord (term_of t, term_of u) = LESS
   784  val concl = Thm.dest_arg o cprop_of
   785  val shuffle1 =
   786    fconv_rule (rewr_conv @{lemma "(a + x == y) == (x == y - (a::real))" by (atomize (full)) (simp add: field_simps) })
   787  val shuffle2 =
   788     fconv_rule (rewr_conv @{lemma "(x + a == y) ==  (x == y - (a::real))" by (atomize (full)) (simp add: field_simps)})
   789  fun substitutable_monomial fvs tm = case term_of tm of
   790     Free(_,@{typ real}) => if not (member (op aconvc) fvs tm) then (Rat.one,tm)
   791                            else raise Failure "substitutable_monomial"
   792   | @{term "op * :: real => _"}$_$(Free _) =>
   793      if RealArith.is_ratconst (Thm.dest_arg1 tm) andalso not (member (op aconvc) fvs (Thm.dest_arg tm))
   794          then (RealArith.dest_ratconst (Thm.dest_arg1 tm),Thm.dest_arg tm) else raise Failure "substitutable_monomial"
   795   | @{term "op + :: real => _"}$_$_ =>
   796        (substitutable_monomial (Thm.add_cterm_frees (Thm.dest_arg tm) fvs) (Thm.dest_arg1 tm)
   797         handle Failure _ => substitutable_monomial (Thm.add_cterm_frees (Thm.dest_arg1 tm) fvs) (Thm.dest_arg tm))
   798   | _ => raise Failure "substitutable_monomial"
   799 
   800   fun isolate_variable v th =
   801    let val w = Thm.dest_arg1 (cprop_of th)
   802    in if v aconvc w then th
   803       else case term_of w of
   804            @{term "op + :: real => _"}$_$_ =>
   805               if Thm.dest_arg1 w aconvc v then shuffle2 th
   806               else isolate_variable v (shuffle1 th)
   807           | _ => error "isolate variable : This should not happen?"
   808    end
   809 in
   810 
   811 fun real_nonlinear_subst_prover prover ctxt =
   812  let
   813   val {add = _, mul = real_poly_mul_conv, neg = _,
   814        pow = _, sub = _, main = real_poly_conv} =
   815       Semiring_Normalizer.semiring_normalizers_ord_wrapper ctxt
   816       (the (Semiring_Normalizer.match ctxt @{cterm "(0::real) + 1"}))
   817      simple_cterm_ord
   818 
   819   fun make_substitution th =
   820    let
   821     val (c,v) = substitutable_monomial [] (Thm.dest_arg1(concl th))
   822     val th1 = Drule.arg_cong_rule (Thm.apply @{cterm "op * :: real => _"} (RealArith.cterm_of_rat (Rat.inv c))) (mk_meta_eq th)
   823     val th2 = fconv_rule (binop_conv real_poly_mul_conv) th1
   824    in fconv_rule (arg_conv real_poly_conv) (isolate_variable v th2)
   825    end
   826    fun oprconv cv ct =
   827     let val g = Thm.dest_fun2 ct
   828     in if g aconvc @{cterm "op <= :: real => _"}
   829          orelse g aconvc @{cterm "op < :: real => _"}
   830        then arg_conv cv ct else arg1_conv cv ct
   831     end
   832   fun mainf cert_choice translator =
   833    let
   834     fun substfirst(eqs,les,lts) =
   835       ((let
   836            val eth = tryfind make_substitution eqs
   837            val modify = fconv_rule (arg_conv (oprconv(subst_conv [eth] then_conv real_poly_conv)))
   838        in  substfirst
   839              (filter_out (fn t => (Thm.dest_arg1 o Thm.dest_arg o cprop_of) t
   840                                    aconvc @{cterm "0::real"}) (map modify eqs),
   841                                    map modify les,map modify lts)
   842        end)
   843        handle Failure  _ => real_nonlinear_prover prover ctxt cert_choice translator (rev eqs, rev les, rev lts))
   844     in substfirst
   845    end
   846 
   847 
   848  in mainf
   849  end
   850 
   851 (* Overall function. *)
   852 
   853 fun real_sos prover ctxt =
   854   RealArith.gen_prover_real_arith ctxt (real_nonlinear_subst_prover prover ctxt)
   855 end;
   856 
   857 val known_sos_constants =
   858   [@{term "op ==>"}, @{term "Trueprop"},
   859    @{term HOL.implies}, @{term HOL.conj}, @{term HOL.disj},
   860    @{term "Not"}, @{term "op = :: bool => _"},
   861    @{term "All :: (real => _) => _"}, @{term "Ex :: (real => _) => _"},
   862    @{term "op = :: real => _"}, @{term "op < :: real => _"},
   863    @{term "op <= :: real => _"},
   864    @{term "op + :: real => _"}, @{term "op - :: real => _"},
   865    @{term "op * :: real => _"}, @{term "uminus :: real => _"},
   866    @{term "op / :: real => _"}, @{term "inverse :: real => _"},
   867    @{term "op ^ :: real => _"}, @{term "abs :: real => _"},
   868    @{term "min :: real => _"}, @{term "max :: real => _"},
   869    @{term "0::real"}, @{term "1::real"},
   870    @{term "numeral :: num => nat"},
   871    @{term "numeral :: num => real"},
   872    @{term "neg_numeral :: num => real"},
   873    @{term "Num.Bit0"}, @{term "Num.Bit1"}, @{term "Num.One"}];
   874 
   875 fun check_sos kcts ct =
   876  let
   877   val t = term_of ct
   878   val _ = if not (null (Term.add_tfrees t [])
   879                   andalso null (Term.add_tvars t []))
   880           then error "SOS: not sos. Additional type varables" else ()
   881   val fs = Term.add_frees t []
   882   val _ = if exists (fn ((_,T)) => not (T = @{typ "real"})) fs
   883           then error "SOS: not sos. Variables with type not real" else ()
   884   val vs = Term.add_vars t []
   885   val _ = if exists (fn ((_,T)) => not (T = @{typ "real"})) vs
   886           then error "SOS: not sos. Variables with type not real" else ()
   887   val ukcs = subtract (fn (t,p) => Const p aconv t) kcts (Term.add_consts t [])
   888   val _ = if  null ukcs then ()
   889               else error ("SOSO: Unknown constants in Subgoal:" ^ commas (map fst ukcs))
   890 in () end
   891 
   892 fun core_sos_tac print_cert prover = SUBPROOF (fn {concl, context, ...} =>
   893   let
   894     val _ = check_sos known_sos_constants concl
   895     val (ths, certificates) = real_sos prover context (Thm.dest_arg concl)
   896     val _ = print_cert certificates
   897   in rtac ths 1 end)
   898 
   899 fun default_SOME _ NONE v = SOME v
   900   | default_SOME _ (SOME v) _ = SOME v;
   901 
   902 fun lift_SOME f NONE a = f a
   903   | lift_SOME _ (SOME a) _ = SOME a;
   904 
   905 
   906 local
   907  val is_numeral = can (HOLogic.dest_number o term_of)
   908 in
   909 fun get_denom b ct = case term_of ct of
   910   @{term "op / :: real => _"} $ _ $ _ =>
   911      if is_numeral (Thm.dest_arg ct) then get_denom b (Thm.dest_arg1 ct)
   912      else default_SOME (get_denom b) (get_denom b (Thm.dest_arg ct))   (Thm.dest_arg ct, b)
   913  | @{term "op < :: real => _"} $ _ $ _ => lift_SOME (get_denom true) (get_denom true (Thm.dest_arg ct)) (Thm.dest_arg1 ct)
   914  | @{term "op <= :: real => _"} $ _ $ _ => lift_SOME (get_denom true) (get_denom true (Thm.dest_arg ct)) (Thm.dest_arg1 ct)
   915  | _ $ _ => lift_SOME (get_denom b) (get_denom b (Thm.dest_fun ct)) (Thm.dest_arg ct)
   916  | _ => NONE
   917 end;
   918 
   919 fun elim_one_denom_tac ctxt =
   920 CSUBGOAL (fn (P,i) =>
   921  case get_denom false P of
   922    NONE => no_tac
   923  | SOME (d,ord) =>
   924      let
   925       val ss = simpset_of ctxt addsimps @{thms field_simps}
   926                addsimps [@{thm nonzero_power_divide}, @{thm power_divide}]
   927       val th = instantiate' [] [SOME d, SOME (Thm.dest_arg P)]
   928          (if ord then @{lemma "(d=0 --> P) & (d>0 --> P) & (d<(0::real) --> P) ==> P" by auto}
   929           else @{lemma "(d=0 --> P) & (d ~= (0::real) --> P) ==> P" by blast})
   930      in rtac th i THEN Simplifier.asm_full_simp_tac ss i end);
   931 
   932 fun elim_denom_tac ctxt i = REPEAT (elim_one_denom_tac ctxt i);
   933 
   934 fun sos_tac print_cert prover ctxt =
   935   Object_Logic.full_atomize_tac THEN'
   936   elim_denom_tac ctxt THEN'
   937   core_sos_tac print_cert prover ctxt;
   938 
   939 end;