src/HOL/Library/Sum_of_Squares/positivstellensatz_tools.ML
author Christian Sternagel
Thu Dec 13 13:11:38 2012 +0100 (2012-12-13)
changeset 50516 ed6b40d15d1c
parent 43946 ba88bb44c192
child 55508 90c42b130652
permissions -rw-r--r--
renamed "emb" to "list_hembeq";
make "list_hembeq" reflexive independent of the base order;
renamed "sub" to "sublisteq";
dropped "transp_on" (state transitivity explicitly instead);
no need to hide "sub" after renaming;
replaced some ASCII symbols by proper Isabelle symbols;
NEWS
wenzelm@41474
     1
(*  Title:      HOL/Library/Sum_of_Squares/positivstellensatz_tools.ML
wenzelm@41474
     2
    Author:     Philipp Meyer, TU Muenchen
Philipp@32645
     3
wenzelm@41474
     4
Functions for generating a certificate from a positivstellensatz and vice
wenzelm@41474
     5
versa.
Philipp@32645
     6
*)
Philipp@32645
     7
Philipp@32645
     8
signature POSITIVSTELLENSATZ_TOOLS =
Philipp@32645
     9
sig
Philipp@32645
    10
  val pss_tree_to_cert : RealArith.pss_tree -> string
Philipp@32645
    11
Philipp@32645
    12
  val cert_to_pss_tree : Proof.context -> string -> RealArith.pss_tree
Philipp@32645
    13
end
Philipp@32645
    14
Philipp@32645
    15
Philipp@32645
    16
structure PositivstellensatzTools : POSITIVSTELLENSATZ_TOOLS =
Philipp@32645
    17
struct
Philipp@32645
    18
Philipp@32645
    19
(*** certificate generation ***)
Philipp@32645
    20
Philipp@32645
    21
fun string_of_rat r =
Philipp@32645
    22
  let
Philipp@32645
    23
    val (nom, den) = Rat.quotient_of_rat r
Philipp@32645
    24
  in
Philipp@32645
    25
    if den = 1 then string_of_int nom
Philipp@32645
    26
    else string_of_int nom ^ "/" ^ string_of_int den
Philipp@32645
    27
  end
Philipp@32645
    28
Philipp@32645
    29
(* map polynomials to strings *)
Philipp@32645
    30
Philipp@32645
    31
fun string_of_varpow x k =
Philipp@32645
    32
  let
Philipp@32645
    33
    val term = term_of x
Philipp@32645
    34
    val name = case term of
Philipp@32645
    35
      Free (n, _) => n
Philipp@32645
    36
    | _ => error "Term in monomial not free variable"
Philipp@32645
    37
  in
Philipp@32645
    38
    if k = 1 then name else name ^ "^" ^ string_of_int k 
Philipp@32645
    39
  end
Philipp@32645
    40
Philipp@32645
    41
fun string_of_monomial m = 
Philipp@32829
    42
 if FuncUtil.Ctermfunc.is_empty m then "1" 
Philipp@32645
    43
 else 
Philipp@32645
    44
  let 
Philipp@32828
    45
   val m' = FuncUtil.dest_monomial m
Philipp@32645
    46
   val vps = fold_rev (fn (x,k) => cons (string_of_varpow x k)) m' [] 
Philipp@32645
    47
  in foldr1 (fn (s, t) => s ^ "*" ^ t) vps
Philipp@32645
    48
  end
Philipp@32645
    49
Philipp@32645
    50
fun string_of_cmonomial (m,c) =
Philipp@32829
    51
  if FuncUtil.Ctermfunc.is_empty m then string_of_rat c
Philipp@32645
    52
  else if c = Rat.one then string_of_monomial m
Philipp@32645
    53
  else (string_of_rat c) ^ "*" ^ (string_of_monomial m);
Philipp@32645
    54
Philipp@32645
    55
fun string_of_poly p = 
Philipp@32829
    56
 if FuncUtil.Monomialfunc.is_empty p then "0" 
Philipp@32645
    57
 else
Philipp@32645
    58
  let 
Philipp@32645
    59
   val cms = map string_of_cmonomial
Philipp@32829
    60
     (sort (prod_ord FuncUtil.monomial_order (K EQUAL)) (FuncUtil.Monomialfunc.dest p))
Philipp@32645
    61
  in foldr1 (fn (t1, t2) => t1 ^ " + " ^ t2) cms
Philipp@32645
    62
  end;
Philipp@32645
    63
Philipp@32828
    64
fun pss_to_cert (RealArith.Axiom_eq i) = "A=" ^ string_of_int i
Philipp@32828
    65
  | pss_to_cert (RealArith.Axiom_le i) = "A<=" ^ string_of_int i
Philipp@32828
    66
  | pss_to_cert (RealArith.Axiom_lt i) = "A<" ^ string_of_int i
Philipp@32828
    67
  | pss_to_cert (RealArith.Rational_eq r) = "R=" ^ string_of_rat r
Philipp@32828
    68
  | pss_to_cert (RealArith.Rational_le r) = "R<=" ^ string_of_rat r
Philipp@32828
    69
  | pss_to_cert (RealArith.Rational_lt r) = "R<" ^ string_of_rat r
Philipp@32828
    70
  | pss_to_cert (RealArith.Square p) = "[" ^ string_of_poly p ^ "]^2"
Philipp@32828
    71
  | pss_to_cert (RealArith.Eqmul (p, pss)) = "([" ^ string_of_poly p ^ "] * " ^ pss_to_cert pss ^ ")"
Philipp@32828
    72
  | pss_to_cert (RealArith.Sum (pss1, pss2)) = "(" ^ pss_to_cert pss1 ^ " + " ^ pss_to_cert pss2 ^ ")"
Philipp@32828
    73
  | pss_to_cert (RealArith.Product (pss1, pss2)) = "(" ^ pss_to_cert pss1 ^ " * " ^ pss_to_cert pss2 ^ ")"
Philipp@32645
    74
Philipp@32828
    75
fun pss_tree_to_cert RealArith.Trivial = "()"
Philipp@32828
    76
  | pss_tree_to_cert (RealArith.Cert pss) = "(" ^ pss_to_cert pss ^ ")"
Philipp@32828
    77
  | pss_tree_to_cert (RealArith.Branch (t1, t2)) = "(" ^ pss_tree_to_cert t1 ^ " & " ^ pss_tree_to_cert t2 ^ ")"
Philipp@32645
    78
Philipp@32645
    79
(*** certificate parsing ***)
Philipp@32645
    80
Philipp@32645
    81
(* basic parser *)
Philipp@32645
    82
Philipp@32646
    83
val str = Scan.this_string
Philipp@32645
    84
Philipp@32645
    85
val number = Scan.repeat1 (Scan.one Symbol.is_ascii_digit >>
Philipp@32645
    86
  (fn s => ord s - ord "0")) >>
Philipp@32645
    87
  foldl1 (fn (n, d) => n * 10 + d)
Philipp@32645
    88
Philipp@32645
    89
val nat = number
Philipp@32646
    90
val int = Scan.optional (str "~" >> K ~1) 1 -- nat >> op *;
Philipp@32646
    91
val rat = int --| str "/" -- int >> Rat.rat_of_quotient
Philipp@32645
    92
val rat_int = rat || int >> Rat.rat_of_int
Philipp@32645
    93
Philipp@32645
    94
(* polynomial parser *)
Philipp@32645
    95
Philipp@32646
    96
fun repeat_sep s f = f ::: Scan.repeat (str s |-- f)
Philipp@32645
    97
Philipp@32645
    98
val parse_id = Scan.one Symbol.is_letter ::: Scan.many Symbol.is_letdig >> implode
Philipp@32645
    99
Philipp@32646
   100
fun parse_varpow ctxt = parse_id -- Scan.optional (str "^" |-- nat) 1 >>
wenzelm@42361
   101
  (fn (x, k) => (cterm_of (Proof_Context.theory_of ctxt) (Free (x, @{typ real})), k)) 
Philipp@32645
   102
Philipp@32645
   103
fun parse_monomial ctxt = repeat_sep "*" (parse_varpow ctxt) >>
wenzelm@33339
   104
  (fn xs => fold FuncUtil.Ctermfunc.update xs FuncUtil.Ctermfunc.empty)
Philipp@32645
   105
Philipp@32645
   106
fun parse_cmonomial ctxt =
Philipp@32646
   107
  rat_int --| str "*" -- (parse_monomial ctxt) >> swap ||
Philipp@32645
   108
  (parse_monomial ctxt) >> (fn m => (m, Rat.one)) ||
Philipp@32829
   109
  rat_int >> (fn r => (FuncUtil.Ctermfunc.empty, r))
Philipp@32645
   110
Philipp@32645
   111
fun parse_poly ctxt = repeat_sep "+" (parse_cmonomial ctxt) >>
wenzelm@33339
   112
  (fn xs => fold FuncUtil.Monomialfunc.update xs FuncUtil.Monomialfunc.empty)
Philipp@32645
   113
Philipp@32645
   114
(* positivstellensatz parser *)
Philipp@32645
   115
Philipp@32645
   116
val parse_axiom =
Philipp@32828
   117
  (str "A=" |-- int >> RealArith.Axiom_eq) ||
Philipp@32828
   118
  (str "A<=" |-- int >> RealArith.Axiom_le) ||
Philipp@32828
   119
  (str "A<" |-- int >> RealArith.Axiom_lt)
Philipp@32645
   120
Philipp@32645
   121
val parse_rational =
Philipp@32828
   122
  (str "R=" |-- rat_int >> RealArith.Rational_eq) ||
Philipp@32828
   123
  (str "R<=" |-- rat_int >> RealArith.Rational_le) ||
Philipp@32828
   124
  (str "R<" |-- rat_int >> RealArith.Rational_lt)
Philipp@32645
   125
Philipp@32645
   126
fun parse_cert ctxt input =
Philipp@32645
   127
  let
Philipp@32645
   128
    val pc = parse_cert ctxt
Philipp@32645
   129
    val pp = parse_poly ctxt
Philipp@32645
   130
  in
Philipp@32645
   131
  (parse_axiom ||
Philipp@32645
   132
   parse_rational ||
Philipp@32828
   133
   str "[" |-- pp --| str "]^2" >> RealArith.Square ||
Philipp@32828
   134
   str "([" |-- pp --| str "]*" -- pc --| str ")" >> RealArith.Eqmul ||
Philipp@32828
   135
   str "(" |-- pc --| str "*" -- pc --| str ")" >> RealArith.Product ||
Philipp@32828
   136
   str "(" |-- pc --| str "+" -- pc --| str ")" >> RealArith.Sum) input
Philipp@32645
   137
  end
Philipp@32645
   138
Philipp@32645
   139
fun parse_cert_tree ctxt input =
Philipp@32645
   140
  let
Philipp@32645
   141
    val pc = parse_cert ctxt
Philipp@32645
   142
    val pt = parse_cert_tree ctxt
Philipp@32645
   143
  in
Philipp@32828
   144
  (str "()" >> K RealArith.Trivial ||
Philipp@32828
   145
   str "(" |-- pc --| str ")" >> RealArith.Cert ||
Philipp@32828
   146
   str "(" |-- pt --| str "&" -- pt --| str ")" >> RealArith.Branch) input
Philipp@32645
   147
  end
Philipp@32645
   148
Philipp@32645
   149
(* scanner *)
Philipp@32645
   150
wenzelm@43946
   151
fun cert_to_pss_tree ctxt input_str =
wenzelm@43946
   152
  Symbol.scanner "Bad certificate" (parse_cert_tree ctxt)
wenzelm@43946
   153
    (filter_out Symbol.is_blank (Symbol.explode input_str))
Philipp@32645
   154
Philipp@32645
   155
end
Philipp@32645
   156