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