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