src/HOL/Library/positivstellensatz.ML
changeset 63205 97b721666890
parent 63201 f151704c08e4
child 63211 0bec0d1d9998
equal deleted inserted replaced
63204:921a5be54132 63205:97b721666890
   336     in foldr1 (fn (s, t) => Thm.apply (Thm.apply @{cterm "op * :: real => _"} s) t) vps
   336     in foldr1 (fn (s, t) => Thm.apply (Thm.apply @{cterm "op * :: real => _"} s) t) vps
   337     end
   337     end
   338 
   338 
   339 fun cterm_of_cmonomial (m,c) =
   339 fun cterm_of_cmonomial (m,c) =
   340   if FuncUtil.Ctermfunc.is_empty m then cterm_of_rat c
   340   if FuncUtil.Ctermfunc.is_empty m then cterm_of_rat c
   341   else if c = Rat.one then cterm_of_monomial m
   341   else if c = @1 then cterm_of_monomial m
   342   else Thm.apply (Thm.apply @{cterm "op *::real => _"} (cterm_of_rat c)) (cterm_of_monomial m);
   342   else Thm.apply (Thm.apply @{cterm "op *::real => _"} (cterm_of_rat c)) (cterm_of_monomial m);
   343 
   343 
   344 fun cterm_of_poly p =
   344 fun cterm_of_poly p =
   345   if FuncUtil.Monomialfunc.is_empty p then @{cterm "0::real"}
   345   if FuncUtil.Monomialfunc.is_empty p then @{cterm "0::real"}
   346   else
   346   else
   583   in f
   583   in f
   584   end;
   584   end;
   585 
   585 
   586 (* A linear arithmetic prover *)
   586 (* A linear arithmetic prover *)
   587 local
   587 local
   588   val linear_add = FuncUtil.Ctermfunc.combine (curry op +) (fn z => z = Rat.zero)
   588   val linear_add = FuncUtil.Ctermfunc.combine (curry op +) (fn z => z = @0)
   589   fun linear_cmul c = FuncUtil.Ctermfunc.map (fn _ => fn x => c * x)
   589   fun linear_cmul c = FuncUtil.Ctermfunc.map (fn _ => fn x => c * x)
   590   val one_tm = @{cterm "1::real"}
   590   val one_tm = @{cterm "1::real"}
   591   fun contradictory p (e,_) = ((FuncUtil.Ctermfunc.is_empty e) andalso not(p Rat.zero)) orelse
   591   fun contradictory p (e,_) = ((FuncUtil.Ctermfunc.is_empty e) andalso not(p @0)) orelse
   592      ((eq_set (op aconvc) (FuncUtil.Ctermfunc.dom e, [one_tm])) andalso
   592      ((eq_set (op aconvc) (FuncUtil.Ctermfunc.dom e, [one_tm])) andalso
   593        not(p(FuncUtil.Ctermfunc.apply e one_tm)))
   593        not(p(FuncUtil.Ctermfunc.apply e one_tm)))
   594 
   594 
   595   fun linear_ineqs vars (les,lts) =
   595   fun linear_ineqs vars (les,lts) =
   596     case find_first (contradictory (fn x => x > Rat.zero)) lts of
   596     case find_first (contradictory (fn x => x > @0)) lts of
   597       SOME r => r
   597       SOME r => r
   598     | NONE =>
   598     | NONE =>
   599       (case find_first (contradictory (fn x => x > Rat.zero)) les of
   599       (case find_first (contradictory (fn x => x > @0)) les of
   600          SOME r => r
   600          SOME r => r
   601        | NONE =>
   601        | NONE =>
   602          if null vars then error "linear_ineqs: no contradiction" else
   602          if null vars then error "linear_ineqs: no contradiction" else
   603          let
   603          let
   604            val ineqs = les @ lts
   604            val ineqs = les @ lts
   605            fun blowup v =
   605            fun blowup v =
   606              length(filter (fn (e,_) => FuncUtil.Ctermfunc.tryapplyd e v Rat.zero = Rat.zero) ineqs) +
   606              length(filter (fn (e,_) => FuncUtil.Ctermfunc.tryapplyd e v @0 = @0) ineqs) +
   607              length(filter (fn (e,_) => FuncUtil.Ctermfunc.tryapplyd e v Rat.zero > Rat.zero) ineqs) *
   607              length(filter (fn (e,_) => FuncUtil.Ctermfunc.tryapplyd e v @0 > @0) ineqs) *
   608              length(filter (fn (e,_) => FuncUtil.Ctermfunc.tryapplyd e v Rat.zero < Rat.zero) ineqs)
   608              length(filter (fn (e,_) => FuncUtil.Ctermfunc.tryapplyd e v @0 < @0) ineqs)
   609            val v = fst(hd(sort (fn ((_,i),(_,j)) => int_ord (i,j))
   609            val v = fst(hd(sort (fn ((_,i),(_,j)) => int_ord (i,j))
   610              (map (fn v => (v,blowup v)) vars)))
   610              (map (fn v => (v,blowup v)) vars)))
   611            fun addup (e1,p1) (e2,p2) acc =
   611            fun addup (e1,p1) (e2,p2) acc =
   612              let
   612              let
   613                val c1 = FuncUtil.Ctermfunc.tryapplyd e1 v Rat.zero
   613                val c1 = FuncUtil.Ctermfunc.tryapplyd e1 v @0
   614                val c2 = FuncUtil.Ctermfunc.tryapplyd e2 v Rat.zero
   614                val c2 = FuncUtil.Ctermfunc.tryapplyd e2 v @0
   615              in
   615              in
   616                if c1 * c2 >= Rat.zero then acc else
   616                if c1 * c2 >= @0 then acc else
   617                let
   617                let
   618                  val e1' = linear_cmul (Rat.abs c2) e1
   618                  val e1' = linear_cmul (Rat.abs c2) e1
   619                  val e2' = linear_cmul (Rat.abs c1) e2
   619                  val e2' = linear_cmul (Rat.abs c1) e2
   620                  val p1' = Product(Rational_lt(Rat.abs c2),p1)
   620                  val p1' = Product(Rational_lt(Rat.abs c2),p1)
   621                  val p2' = Product(Rational_lt(Rat.abs c1),p2)
   621                  val p2' = Product(Rational_lt(Rat.abs c1),p2)
   622                in (linear_add e1' e2',Sum(p1',p2'))::acc
   622                in (linear_add e1' e2',Sum(p1',p2'))::acc
   623                end
   623                end
   624              end
   624              end
   625            val (les0,les1) =
   625            val (les0,les1) =
   626              List.partition (fn (e,_) => FuncUtil.Ctermfunc.tryapplyd e v Rat.zero = Rat.zero) les
   626              List.partition (fn (e,_) => FuncUtil.Ctermfunc.tryapplyd e v @0 = @0) les
   627            val (lts0,lts1) =
   627            val (lts0,lts1) =
   628              List.partition (fn (e,_) => FuncUtil.Ctermfunc.tryapplyd e v Rat.zero = Rat.zero) lts
   628              List.partition (fn (e,_) => FuncUtil.Ctermfunc.tryapplyd e v @0 = @0) lts
   629            val (lesp,lesn) =
   629            val (lesp,lesn) =
   630              List.partition (fn (e,_) => FuncUtil.Ctermfunc.tryapplyd e v Rat.zero > Rat.zero) les1
   630              List.partition (fn (e,_) => FuncUtil.Ctermfunc.tryapplyd e v @0 > @0) les1
   631            val (ltsp,ltsn) =
   631            val (ltsp,ltsn) =
   632              List.partition (fn (e,_) => FuncUtil.Ctermfunc.tryapplyd e v Rat.zero > Rat.zero) lts1
   632              List.partition (fn (e,_) => FuncUtil.Ctermfunc.tryapplyd e v @0 > @0) lts1
   633            val les' = fold_rev (fn ep1 => fold_rev (addup ep1) lesp) lesn les0
   633            val les' = fold_rev (fn ep1 => fold_rev (addup ep1) lesp) lesn les0
   634            val lts' = fold_rev (fn ep1 => fold_rev (addup ep1) (lesp@ltsp)) ltsn
   634            val lts' = fold_rev (fn ep1 => fold_rev (addup ep1) (lesp@ltsp)) ltsn
   635                       (fold_rev (fn ep1 => fold_rev (addup ep1) (lesn@ltsn)) ltsp lts0)
   635                       (fold_rev (fn ep1 => fold_rev (addup ep1) (lesn@ltsn)) ltsp lts0)
   636          in linear_ineqs (remove (op aconvc) v vars) (les',lts')
   636          in linear_ineqs (remove (op aconvc) v vars) (les',lts')
   637          end)
   637          end)
   638 
   638 
   639   fun linear_eqs(eqs,les,lts) =
   639   fun linear_eqs(eqs,les,lts) =
   640     case find_first (contradictory (fn x => x = Rat.zero)) eqs of
   640     case find_first (contradictory (fn x => x = @0)) eqs of
   641       SOME r => r
   641       SOME r => r
   642     | NONE =>
   642     | NONE =>
   643       (case eqs of
   643       (case eqs of
   644          [] =>
   644          [] =>
   645          let val vars = remove (op aconvc) one_tm
   645          let val vars = remove (op aconvc) one_tm
   648        | (e,p)::es =>
   648        | (e,p)::es =>
   649          if FuncUtil.Ctermfunc.is_empty e then linear_eqs (es,les,lts) else
   649          if FuncUtil.Ctermfunc.is_empty e then linear_eqs (es,les,lts) else
   650          let
   650          let
   651            val (x,c) = FuncUtil.Ctermfunc.choose (FuncUtil.Ctermfunc.delete_safe one_tm e)
   651            val (x,c) = FuncUtil.Ctermfunc.choose (FuncUtil.Ctermfunc.delete_safe one_tm e)
   652            fun xform (inp as (t,q)) =
   652            fun xform (inp as (t,q)) =
   653              let val d = FuncUtil.Ctermfunc.tryapplyd t x Rat.zero in
   653              let val d = FuncUtil.Ctermfunc.tryapplyd t x @0 in
   654                if d = Rat.zero then inp else
   654                if d = @0 then inp else
   655                let
   655                let
   656                  val k = (Rat.neg d) * Rat.abs c / c
   656                  val k = (Rat.neg d) * Rat.abs c / c
   657                  val e' = linear_cmul k e
   657                  val e' = linear_cmul k e
   658                  val t' = linear_cmul (Rat.abs c) t
   658                  val t' = linear_cmul (Rat.abs c) t
   659                  val p' = Eqmul(FuncUtil.Monomialfunc.onefunc (FuncUtil.Ctermfunc.empty, k),p)
   659                  val p' = Eqmul(FuncUtil.Monomialfunc.onefunc (FuncUtil.Ctermfunc.empty, k),p)
   672     in linear_eqs(eqs,les,lts)
   672     in linear_eqs(eqs,les,lts)
   673     end
   673     end
   674 
   674 
   675   fun lin_of_hol ct =
   675   fun lin_of_hol ct =
   676     if ct aconvc @{cterm "0::real"} then FuncUtil.Ctermfunc.empty
   676     if ct aconvc @{cterm "0::real"} then FuncUtil.Ctermfunc.empty
   677     else if not (is_comb ct) then FuncUtil.Ctermfunc.onefunc (ct, Rat.one)
   677     else if not (is_comb ct) then FuncUtil.Ctermfunc.onefunc (ct, @1)
   678     else if is_ratconst ct then FuncUtil.Ctermfunc.onefunc (one_tm, dest_ratconst ct)
   678     else if is_ratconst ct then FuncUtil.Ctermfunc.onefunc (one_tm, dest_ratconst ct)
   679     else
   679     else
   680       let val (lop,r) = Thm.dest_comb ct
   680       let val (lop,r) = Thm.dest_comb ct
   681       in
   681       in
   682         if not (is_comb lop) then FuncUtil.Ctermfunc.onefunc (ct, Rat.one)
   682         if not (is_comb lop) then FuncUtil.Ctermfunc.onefunc (ct, @1)
   683         else
   683         else
   684           let val (opr,l) = Thm.dest_comb lop
   684           let val (opr,l) = Thm.dest_comb lop
   685           in
   685           in
   686             if opr aconvc @{cterm "op + :: real =>_"}
   686             if opr aconvc @{cterm "op + :: real =>_"}
   687             then linear_add (lin_of_hol l) (lin_of_hol r)
   687             then linear_add (lin_of_hol l) (lin_of_hol r)
   688             else if opr aconvc @{cterm "op * :: real =>_"}
   688             else if opr aconvc @{cterm "op * :: real =>_"}
   689                     andalso is_ratconst l then FuncUtil.Ctermfunc.onefunc (r, dest_ratconst l)
   689                     andalso is_ratconst l then FuncUtil.Ctermfunc.onefunc (r, dest_ratconst l)
   690             else FuncUtil.Ctermfunc.onefunc (ct, Rat.one)
   690             else FuncUtil.Ctermfunc.onefunc (ct, @1)
   691           end
   691           end
   692       end
   692       end
   693 
   693 
   694   fun is_alien ct =
   694   fun is_alien ct =
   695     case Thm.term_of ct of
   695     case Thm.term_of ct of
   705     val le_pols = map rhs le
   705     val le_pols = map rhs le
   706     val lt_pols = map rhs lt
   706     val lt_pols = map rhs lt
   707     val aliens = filter is_alien
   707     val aliens = filter is_alien
   708       (fold_rev (union (op aconvc) o FuncUtil.Ctermfunc.dom)
   708       (fold_rev (union (op aconvc) o FuncUtil.Ctermfunc.dom)
   709                 (eq_pols @ le_pols @ lt_pols) [])
   709                 (eq_pols @ le_pols @ lt_pols) [])
   710     val le_pols' = le_pols @ map (fn v => FuncUtil.Ctermfunc.onefunc (v,Rat.one)) aliens
   710     val le_pols' = le_pols @ map (fn v => FuncUtil.Ctermfunc.onefunc (v,@1)) aliens
   711     val (_,proof) = linear_prover (eq_pols,le_pols',lt_pols)
   711     val (_,proof) = linear_prover (eq_pols,le_pols',lt_pols)
   712     val le' = le @ map (fn a => Thm.instantiate' [] [SOME (Thm.dest_arg a)] @{thm of_nat_0_le_iff}) aliens
   712     val le' = le @ map (fn a => Thm.instantiate' [] [SOME (Thm.dest_arg a)] @{thm of_nat_0_le_iff}) aliens
   713   in ((translator (eq,le',lt) proof), Trivial)
   713   in ((translator (eq,le',lt) proof), Trivial)
   714   end
   714   end
   715 end;
   715 end;