src/HOL/Library/positivstellensatz.ML
changeset 32645 1cc5b24f5a01
parent 32402 5731300da417
child 32646 962b4354ed90
     1.1 --- a/src/HOL/Library/positivstellensatz.ML	Tue Sep 22 20:25:31 2009 +0200
     1.2 +++ b/src/HOL/Library/positivstellensatz.ML	Mon Sep 21 15:05:26 2009 +0200
     1.3 @@ -1,10 +1,11 @@
     1.4 -(* Title:      Library/positivstellensatz
     1.5 +(* Title:      Library/Sum_Of_Squares/positivstellensatz
     1.6     Author:     Amine Chaieb, University of Cambridge
     1.7     Description: A generic arithmetic prover based on Positivstellensatz certificates --- 
     1.8      also implements Fourrier-Motzkin elimination as a special case Fourrier-Motzkin elimination.
     1.9  *)
    1.10  
    1.11  (* A functor for finite mappings based on Tables *)
    1.12 +
    1.13  signature FUNC = 
    1.14  sig
    1.15   type 'a T
    1.16 @@ -75,27 +76,54 @@
    1.17  end
    1.18  end;
    1.19  
    1.20 +(* Some standard functors and utility functions for them *)
    1.21 +
    1.22 +structure FuncUtil =
    1.23 +struct
    1.24 +
    1.25 +fun increasing f ord (x,y) = ord (f x, f y);
    1.26 +
    1.27  structure Intfunc = FuncFun(type key = int val ord = int_ord);
    1.28 +structure Ratfunc = FuncFun(type key = Rat.rat val ord = Rat.ord);
    1.29 +structure Intpairfunc = FuncFun(type key = int*int val ord = prod_ord int_ord int_ord);
    1.30  structure Symfunc = FuncFun(type key = string val ord = fast_string_ord);
    1.31  structure Termfunc = FuncFun(type key = term val ord = TermOrd.fast_term_ord);
    1.32 -structure Ctermfunc = FuncFun(type key = cterm val ord = (fn (s,t) => TermOrd.fast_term_ord(term_of s, term_of t)));
    1.33 +
    1.34 +val cterm_ord = (fn (s,t) => TermOrd.fast_term_ord(term_of s, term_of t))
    1.35 +
    1.36 +structure Ctermfunc = FuncFun(type key = cterm val ord = cterm_ord);
    1.37 +
    1.38 +type monomial = int Ctermfunc.T;
    1.39  
    1.40 -structure Ratfunc = FuncFun(type key = Rat.rat val ord = Rat.ord);
    1.41 +fun monomial_ord (m1,m2) = list_ord (prod_ord cterm_ord int_ord) (Ctermfunc.graph m1, Ctermfunc.graph m2)
    1.42 +
    1.43 +structure Monomialfunc = FuncFun(type key = monomial val ord = monomial_ord)
    1.44  
    1.45 +type poly = Rat.rat Monomialfunc.T;
    1.46 +
    1.47 +(* The ordering so we can create canonical HOL polynomials.                  *)
    1.48  
    1.49 -    (* Some useful derived rules *)
    1.50 -fun deduct_antisym_rule tha thb = 
    1.51 -    equal_intr (implies_intr (cprop_of thb) tha) 
    1.52 -     (implies_intr (cprop_of tha) thb);
    1.53 +fun dest_monomial mon = sort (increasing fst cterm_ord) (Ctermfunc.graph mon);
    1.54  
    1.55 -fun prove_hyp tha thb = 
    1.56 -  if exists (curry op aconv (concl_of tha)) (#hyps (rep_thm thb)) 
    1.57 -  then equal_elim (symmetric (deduct_antisym_rule tha thb)) tha else thb;
    1.58 +fun monomial_order (m1,m2) =
    1.59 + if Ctermfunc.is_undefined m2 then LESS 
    1.60 + else if Ctermfunc.is_undefined m1 then GREATER 
    1.61 + else
    1.62 +  let val mon1 = dest_monomial m1 
    1.63 +      val mon2 = dest_monomial m2
    1.64 +      val deg1 = fold (curry op + o snd) mon1 0
    1.65 +      val deg2 = fold (curry op + o snd) mon2 0 
    1.66 +  in if deg1 < deg2 then GREATER else if deg1 > deg2 then LESS
    1.67 +     else list_ord (prod_ord cterm_ord int_ord) (mon1,mon2)
    1.68 +  end;
    1.69  
    1.70 +end
    1.71  
    1.72 +(* positivstellensatz datatype and prover generation *)
    1.73  
    1.74  signature REAL_ARITH = 
    1.75  sig
    1.76 +  
    1.77    datatype positivstellensatz =
    1.78     Axiom_eq of int
    1.79   | Axiom_le of int
    1.80 @@ -103,34 +131,41 @@
    1.81   | Rational_eq of Rat.rat
    1.82   | Rational_le of Rat.rat
    1.83   | Rational_lt of Rat.rat
    1.84 - | Square of cterm
    1.85 - | Eqmul of cterm * positivstellensatz
    1.86 + | Square of FuncUtil.poly
    1.87 + | Eqmul of FuncUtil.poly * positivstellensatz
    1.88   | Sum of positivstellensatz * positivstellensatz
    1.89   | Product of positivstellensatz * positivstellensatz;
    1.90  
    1.91 +datatype pss_tree = Trivial | Cert of positivstellensatz | Branch of pss_tree * pss_tree
    1.92 +
    1.93 +datatype tree_choice = Left | Right
    1.94 +
    1.95 +type prover = tree_choice list -> 
    1.96 +  (thm list * thm list * thm list -> positivstellensatz -> thm) ->
    1.97 +  thm list * thm list * thm list -> thm * pss_tree
    1.98 +type cert_conv = cterm -> thm * pss_tree
    1.99 +
   1.100  val gen_gen_real_arith :
   1.101 -  Proof.context -> (Rat.rat -> Thm.cterm) * conv * conv * conv * 
   1.102 -   conv * conv * conv * conv * conv * conv * 
   1.103 -    ( (thm list * thm list * thm list -> positivstellensatz -> thm) ->
   1.104 -        thm list * thm list * thm list -> thm) -> conv
   1.105 -val real_linear_prover : 
   1.106 -  (thm list * thm list * thm list -> positivstellensatz -> thm) ->
   1.107 -   thm list * thm list * thm list -> thm
   1.108 +  Proof.context -> (Rat.rat -> Thm.cterm) * conv * conv * conv *
   1.109 +   conv * conv * conv * conv * conv * conv * prover -> cert_conv
   1.110 +val real_linear_prover : (thm list * thm list * thm list -> positivstellensatz -> thm) ->
   1.111 +  thm list * thm list * thm list -> thm * pss_tree
   1.112  
   1.113  val gen_real_arith : Proof.context ->
   1.114 -   (Rat.rat -> cterm) * conv * conv * conv * conv * conv * conv * conv *
   1.115 -   ( (thm list * thm list * thm list -> positivstellensatz -> thm) ->
   1.116 -       thm list * thm list * thm list -> thm) -> conv
   1.117 -val gen_prover_real_arith : Proof.context ->
   1.118 -   ((thm list * thm list * thm list -> positivstellensatz -> thm) ->
   1.119 -     thm list * thm list * thm list -> thm) -> conv
   1.120 -val real_arith : Proof.context -> conv
   1.121 +  (Rat.rat -> Thm.cterm) * conv * conv * conv * conv * conv * conv * conv * prover -> cert_conv
   1.122 +
   1.123 +val gen_prover_real_arith : Proof.context -> prover -> cert_conv
   1.124 +
   1.125 +val is_ratconst : Thm.cterm -> bool
   1.126 +val dest_ratconst : Thm.cterm -> Rat.rat
   1.127 +val cterm_of_rat : Rat.rat -> Thm.cterm
   1.128 +
   1.129  end
   1.130  
   1.131 -structure RealArith (* : REAL_ARITH *)=
   1.132 +structure RealArith : REAL_ARITH =
   1.133  struct
   1.134  
   1.135 - open Conv Thm;;
   1.136 + open Conv Thm FuncUtil;;
   1.137  (* ------------------------------------------------------------------------- *)
   1.138  (* Data structure for Positivstellensatz refutations.                        *)
   1.139  (* ------------------------------------------------------------------------- *)
   1.140 @@ -142,12 +177,18 @@
   1.141   | Rational_eq of Rat.rat
   1.142   | Rational_le of Rat.rat
   1.143   | Rational_lt of Rat.rat
   1.144 - | Square of cterm
   1.145 - | Eqmul of cterm * positivstellensatz
   1.146 + | Square of FuncUtil.poly
   1.147 + | Eqmul of FuncUtil.poly * positivstellensatz
   1.148   | Sum of positivstellensatz * positivstellensatz
   1.149   | Product of positivstellensatz * positivstellensatz;
   1.150           (* Theorems used in the procedure *)
   1.151  
   1.152 +datatype pss_tree = Trivial | Cert of positivstellensatz | Branch of pss_tree * pss_tree
   1.153 +datatype tree_choice = Left | Right
   1.154 +type prover = tree_choice list -> 
   1.155 +  (thm list * thm list * thm list -> positivstellensatz -> thm) ->
   1.156 +  thm list * thm list * thm list -> thm * pss_tree
   1.157 +type cert_conv = cterm -> thm * pss_tree
   1.158  
   1.159  val my_eqs = ref ([] : thm list);
   1.160  val my_les = ref ([] : thm list);
   1.161 @@ -164,6 +205,16 @@
   1.162  val my_poly_add_conv = ref no_conv;
   1.163  val my_poly_mul_conv = ref no_conv;
   1.164  
   1.165 +
   1.166 +    (* Some useful derived rules *)
   1.167 +fun deduct_antisym_rule tha thb = 
   1.168 +    equal_intr (implies_intr (cprop_of thb) tha) 
   1.169 +     (implies_intr (cprop_of tha) thb);
   1.170 +
   1.171 +fun prove_hyp tha thb = 
   1.172 +  if exists (curry op aconv (concl_of tha)) (#hyps (rep_thm thb)) 
   1.173 +  then equal_elim (symmetric (deduct_antisym_rule tha thb)) tha else thb;
   1.174 +
   1.175  fun conjunctions th = case try Conjunction.elim th of
   1.176     SOME (th1,th2) => (conjunctions th1) @ conjunctions th2
   1.177   | NONE => [th];
   1.178 @@ -292,7 +343,6 @@
   1.179   | Abs (_,_,t') => find_cterm p (Thm.dest_abs NONE t |> snd)
   1.180   | _ => raise CTERM ("find_cterm",[t]);
   1.181  
   1.182 -
   1.183      (* Some conversions-related stuff which has been forbidden entrance into Pure/conv.ML*)
   1.184  fun instantiate_cterm' ty tms = Drule.cterm_rule (Drule.instantiate' ty tms)
   1.185  fun is_comb t = case (term_of t) of _$_ => true | _ => false;
   1.186 @@ -300,6 +350,39 @@
   1.187  fun is_binop ct ct' = ct aconvc (Thm.dest_fun (Thm.dest_fun ct'))
   1.188    handle CTERM _ => false;
   1.189  
   1.190 +
   1.191 +(* Map back polynomials to HOL.                         *)
   1.192 +
   1.193 +local
   1.194 + open Thm Numeral
   1.195 +in
   1.196 +
   1.197 +fun cterm_of_varpow x k = if k = 1 then x else capply (capply @{cterm "op ^ :: real => _"} x) 
   1.198 +  (mk_cnumber @{ctyp nat} k)
   1.199 +
   1.200 +fun cterm_of_monomial m = 
   1.201 + if Ctermfunc.is_undefined m then @{cterm "1::real"} 
   1.202 + else 
   1.203 +  let 
   1.204 +   val m' = dest_monomial m
   1.205 +   val vps = fold_rev (fn (x,k) => cons (cterm_of_varpow x k)) m' [] 
   1.206 +  in foldr1 (fn (s, t) => capply (capply @{cterm "op * :: real => _"} s) t) vps
   1.207 +  end
   1.208 +
   1.209 +fun cterm_of_cmonomial (m,c) = if Ctermfunc.is_undefined m then cterm_of_rat c
   1.210 +    else if c = Rat.one then cterm_of_monomial m
   1.211 +    else capply (capply @{cterm "op *::real => _"} (cterm_of_rat c)) (cterm_of_monomial m);
   1.212 +
   1.213 +fun cterm_of_poly p = 
   1.214 + if Monomialfunc.is_undefined p then @{cterm "0::real"} 
   1.215 + else
   1.216 +  let 
   1.217 +   val cms = map cterm_of_cmonomial
   1.218 +     (sort (prod_ord monomial_order (K EQUAL)) (Monomialfunc.graph p))
   1.219 +  in foldr1 (fn (t1, t2) => capply(capply @{cterm "op + :: real => _"} t1) t2) cms
   1.220 +  end;
   1.221 +
   1.222 +end;
   1.223      (* A general real arithmetic prover *)
   1.224  
   1.225  fun gen_gen_real_arith ctxt (mk_numeric,
   1.226 @@ -368,8 +451,8 @@
   1.227        | Rational_lt x => eqT_elim(numeric_gt_conv(capply @{cterm Trueprop} 
   1.228                        (capply (capply @{cterm "op <::real => _"} @{cterm "0::real"})
   1.229                          (mk_numeric x))))
   1.230 -      | Square t => square_rule t
   1.231 -      | Eqmul(t,p) => emul_rule t (translate p)
   1.232 +      | Square pt => square_rule (cterm_of_poly pt)
   1.233 +      | Eqmul(pt,p) => emul_rule (cterm_of_poly pt) (translate p)
   1.234        | Sum(p1,p2) => add_rule (translate p1) (translate p2)
   1.235        | Product(p1,p2) => mul_rule (translate p1) (translate p2)
   1.236     in fconv_rule (first_conv [numeric_ge_conv, numeric_gt_conv, numeric_eq_conv, all_conv]) 
   1.237 @@ -394,13 +477,13 @@
   1.238         val _ = if c aconvc (concl th2) then () else error "disj_cases : conclusions not alpha convertible"
   1.239     in implies_elim (implies_elim (implies_elim (instantiate' [] (map SOME [p,q,c]) @{thm disjE}) th) (implies_intr (capply @{cterm Trueprop} p) th1)) (implies_intr (capply @{cterm Trueprop} q) th2)
   1.240     end
   1.241 - fun overall dun ths = case ths of
   1.242 + fun overall cert_choice dun ths = case ths of
   1.243    [] =>
   1.244     let 
   1.245      val (eq,ne) = List.partition (is_req o concl) dun
   1.246       val (le,nl) = List.partition (is_ge o concl) ne
   1.247       val lt = filter (is_gt o concl) nl 
   1.248 -    in prover hol_of_positivstellensatz (eq,le,lt) end
   1.249 +    in prover (rev cert_choice) hol_of_positivstellensatz (eq,le,lt) end
   1.250   | th::oths =>
   1.251     let 
   1.252      val ct = concl th 
   1.253 @@ -408,13 +491,13 @@
   1.254      if is_conj ct  then
   1.255       let 
   1.256        val (th1,th2) = conj_pair th in
   1.257 -      overall dun (th1::th2::oths) end
   1.258 +      overall cert_choice dun (th1::th2::oths) end
   1.259      else if is_disj ct then
   1.260        let 
   1.261 -       val th1 = overall dun (assume (capply @{cterm Trueprop} (dest_arg1 ct))::oths)
   1.262 -       val th2 = overall dun (assume (capply @{cterm Trueprop} (dest_arg ct))::oths)
   1.263 -      in disj_cases th th1 th2 end
   1.264 -   else overall (th::dun) oths
   1.265 +       val (th1, cert1) = overall (Left::cert_choice) dun (assume (capply @{cterm Trueprop} (dest_arg1 ct))::oths)
   1.266 +       val (th2, cert2) = overall (Right::cert_choice) dun (assume (capply @{cterm Trueprop} (dest_arg ct))::oths)
   1.267 +      in (disj_cases th th1 th2, Branch (cert1, cert2)) end
   1.268 +   else overall cert_choice (th::dun) oths
   1.269    end
   1.270    fun dest_binary b ct = if is_binop b ct then dest_binop ct 
   1.271                           else raise CTERM ("dest_binary",[b,ct])
   1.272 @@ -496,16 +579,16 @@
   1.273    val nct = capply @{cterm Trueprop} (capply @{cterm "Not"} ct)
   1.274    val th0 = (init_conv then_conv arg_conv nnf_norm_conv') nct
   1.275    val tm0 = dest_arg (rhs_of th0)
   1.276 -  val th = if tm0 aconvc @{cterm False} then equal_implies_1_rule th0 else
   1.277 +  val (th, certificates) = if tm0 aconvc @{cterm False} then (equal_implies_1_rule th0, Trivial) else
   1.278     let 
   1.279      val (evs,bod) = strip_exists tm0
   1.280      val (avs,ibod) = strip_forall bod
   1.281      val th1 = Drule.arg_cong_rule @{cterm Trueprop} (fold mk_forall avs (absremover ibod))
   1.282 -    val th2 = overall [] [specl avs (assume (rhs_of th1))]
   1.283 +    val (th2, certs) = overall [] [] [specl avs (assume (rhs_of th1))]
   1.284      val th3 = fold simple_choose evs (prove_hyp (equal_elim th1 (assume (capply @{cterm Trueprop} bod))) th2)
   1.285 -   in  Drule.implies_intr_hyps (prove_hyp (equal_elim th0 (assume nct)) th3)
   1.286 +   in (Drule.implies_intr_hyps (prove_hyp (equal_elim th0 (assume nct)) th3), certs)
   1.287     end
   1.288 -  in implies_elim (instantiate' [] [SOME ct] pth_final) th
   1.289 +  in (implies_elim (instantiate' [] [SOME ct] pth_final) th, certificates)
   1.290   end
   1.291  in f
   1.292  end;
   1.293 @@ -580,7 +663,7 @@
   1.294           val k = (Rat.neg d) */ Rat.abs c // c
   1.295           val e' = linear_cmul k e
   1.296           val t' = linear_cmul (Rat.abs c) t
   1.297 -         val p' = Eqmul(cterm_of_rat k,p)
   1.298 +         val p' = Eqmul(Monomialfunc.onefunc (Ctermfunc.undefined, k),p)
   1.299           val q' = Product(Rational_lt(Rat.abs c),q) 
   1.300          in (linear_add e' t',Sum(p',q')) 
   1.301          end 
   1.302 @@ -632,7 +715,7 @@
   1.303    val le_pols' = le_pols @ map (fn v => Ctermfunc.onefunc (v,Rat.one)) aliens
   1.304    val (_,proof) = linear_prover (eq_pols,le_pols',lt_pols)
   1.305    val le' = le @ map (fn a => instantiate' [] [SOME (dest_arg a)] @{thm real_of_nat_ge_zero}) aliens 
   1.306 - in (translator (eq,le',lt) proof) : thm
   1.307 + in ((translator (eq,le',lt) proof), Trivial)
   1.308   end
   1.309  end;
   1.310  
   1.311 @@ -698,5 +781,4 @@
   1.312      main,neg,add,mul, prover)
   1.313  end;
   1.314  
   1.315 -fun real_arith ctxt = gen_prover_real_arith ctxt real_linear_prover;
   1.316  end