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