Theory LinArith

(*  Author:     Tobias Nipkow, 2007  *)

header{* Linear real arithmetic *}

theory LinArith
imports QE "~~/src/HOL/Library/ListVector" Complex_Main
begin

declare iprod_assoc[simp]

subsection{*Basics*}

subsubsection{*Syntax and Semantics*}

datatype atom = Less real "real list" | Eq real "real list"

fun is_Less :: "atom => bool" where
"is_Less (Less r rs) = True" |
"is_Less f = False"

abbreviation "is_Eq ≡ Not o is_Less"

lemma is_Less_iff: "is_Less f = (∃r rs. f = Less r rs)"
by(induct f) auto

lemma is_Eq_iff: "(∀i j. a ≠ Less i j) = (∃i j. a = Eq i j)"
by(cases a) auto

fun neg⇣R :: "atom => atom fm" where
"neg⇣R (Less r t) = Or (Atom(Less (-r) (-t))) (Atom(Eq r t))" |
"neg⇣R (Eq r t) = Or (Atom(Less r t)) (Atom(Less (-r) (-t)))"

fun hd_coeff :: "atom => real" where
"hd_coeff (Less r cs) = (case cs of [] => 0 | c#_ => c)" |
"hd_coeff (Eq r cs) = (case cs of [] => 0 | c#_ => c)"

definition "depends⇣R a = (hd_coeff a ≠ 0)"

fun decr⇣R :: "atom => atom" where
"decr⇣R (Less r rs) = Less r (tl rs)" |
"decr⇣R (Eq r rs) = Eq r (tl rs)"

fun I⇣R :: "atom => real list => bool" where
"I⇣R (Less r cs) xs = (r < ⟨cs,xs⟩)" |
"I⇣R (Eq r cs) xs = (r = ⟨cs,xs⟩)"

definition "atoms⇣0 = ATOM.atoms⇣0 depends⇣R"
(* FIXME !!! (incl: display should hide params)*)

interpretation R!: ATOM neg⇣R "(λa. True)" I⇣R depends⇣R decr⇣R
  where "ATOM.atoms⇣0 depends⇣R = atoms⇣0"
proof -
  case goal1
  thus ?case
    apply(unfold_locales)
        apply(case_tac a)
         apply simp_all
     apply(case_tac a)
      apply simp_all
      apply arith
     apply arith
    apply(case_tac a)
    apply(simp_all add:depends⇣R_def split:list.splits)
    done
next
  case goal2
  thus ?case by(simp add:atoms⇣0_def)
qed

setup {* Sign.revert_abbrev "" @{const_abbrev R.I} *}
setup {* Sign.revert_abbrev "" @{const_abbrev R.lift_nnf_qe} *}


subsubsection{*Shared constructions*}

fun combine :: "(real * real list) => (real * real list) => atom" where
"combine (r⇣1,cs⇣1) (r⇣2,cs⇣2) = Less (r⇣1-r⇣2) (cs⇣2 - cs⇣1)"

(* More efficient combination than rhs of combine_conv below; unused
lemma combine_code:
  "combine (r1,c1,cs1) (r2,c2,cs2) =
  Less (c1*r2-c2*r1) (zipwith0 (%x y. c1*y-c2*x) cs1 cs2)"
apply(rule sym)
apply(simp)
apply(induct cs1 arbitrary:cs2)
 apply simp
apply(case_tac cs2)
 apply simp
apply simp
done*)

definition "lbounds as = [(r/c, (-1/c) *⇩s cs). Less r (c#cs) \<leftarrow> as, c>0]"
definition "ubounds as = [(r/c, (-1/c) *⇩s cs). Less r (c#cs) \<leftarrow> as, c<0]"
definition "ebounds as = [(r/c, (-1/c) *⇩s cs). Eq r (c#cs) \<leftarrow> as, c≠0]"

lemma set_lbounds:
 "set(lbounds as) = {(r/c, (-1/c) *⇩s cs)|r c cs. Less r (c#cs) : set as ∧ c>0}"
by(force simp: lbounds_def split:list.splits atom.splits if_splits)
lemma set_ubounds:
 "set(ubounds as) = {(r/c, (-1/c) *⇩s cs)|r c cs. Less r (c#cs) : set as ∧ c<0}"
by(force simp: ubounds_def split:list.splits atom.splits if_splits)
lemma set_ebounds:
  "set(ebounds as) = {(r/c, (-1/c) *⇩s cs)|r c cs. Eq r (c#cs) : set as ∧ c≠0}"
by(force simp: ebounds_def split:list.splits atom.splits if_splits)


abbreviation EQ where
"EQ f xs ≡ {(r - ⟨cs,xs⟩)/c|r c cs. Eq r (c#cs) : set(R.atoms⇣0 f) ∧ c≠0}"
abbreviation LB where
"LB f xs ≡ {(r - ⟨cs,xs⟩)/c|r c cs. Less r (c#cs) : set(R.atoms⇣0 f) ∧ c>0}"
abbreviation UB where
"UB f xs ≡ {(r - ⟨cs,xs⟩)/c|r c cs. Less r (c#cs) : set(R.atoms⇣0 f) ∧ c<0}"


fun asubst :: "real * real list => atom => atom" where
"asubst (r,cs) (Less s (d#ds)) = Less (s - d*r) (d *⇩s cs + ds)" |
"asubst (r,cs) (Eq s (d#ds)) = Eq (s - d*r) (d *⇩s cs + ds)" |
"asubst (r,cs) (Less s []) = Less s []" |
"asubst (r,cs) (Eq s []) = Eq s []"

abbreviation "subst φ rcs ≡ map⇘fm⇙ (asubst rcs) φ"

definition eval :: "real * real list => real list => real" where
"eval rcs xs = fst rcs + ⟨snd rcs,xs⟩"

lemma I_asubst:
 "I⇣R (asubst t a) xs = I⇣R a (eval t xs # xs)"
proof(cases a)
  case (Less r cs)
  thus ?thesis by(cases t, cases cs,
    simp_all add:eval_def distrib_left iprod_left_add_distrib)
    arith
next
  case (Eq r cs)
  thus ?thesis
    by(cases t, cases cs, simp_all add:eval_def distrib_left iprod_left_add_distrib)
      arith
qed

lemma I_subst:
  "qfree φ ==> R.I (subst φ t) xs = R.I φ (eval t xs # xs)"
by(induct φ)(simp_all add:I_asubst)
lemma I_subst_pretty:
  "qfree φ ==> R.I (subst φ (r,cs)) xs = R.I φ ((r + ⟨cs,xs⟩) # xs)"
by(simp add:I_subst eval_def)


fun min_inf :: "atom fm => atom fm" ("inf⇣-") where
"inf⇣- (And φ⇣1 φ⇣2) = and (inf⇣- φ⇣1) (inf⇣- φ⇣2)" |
"inf⇣- (Or φ⇣1 φ⇣2) = or (inf⇣- φ⇣1) (inf⇣- φ⇣2)" |
"inf⇣- (Atom(Less r (c#cs))) =
  (if c<0 then TrueF else if c>0 then FalseF else Atom(Less r cs))" |
"inf⇣- (Atom(Eq r (c#cs))) = (if c=0 then Atom(Eq r cs) else FalseF)" |
"inf⇣- φ = φ"

fun plus_inf :: "atom fm => atom fm" ("inf⇣+") where
"inf⇣+ (And φ⇣1 φ⇣2) = and (inf⇣+ φ⇣1) (inf⇣+ φ⇣2)" |
"inf⇣+ (Or φ⇣1 φ⇣2) = or (inf⇣+ φ⇣1) (inf⇣+ φ⇣2)" |
"inf⇣+ (Atom(Less r (c#cs))) =
  (if c>0 then TrueF else if c<0 then FalseF else Atom(Less r cs))" |
"inf⇣+ (Atom(Eq r (c#cs))) = (if c=0 then Atom(Eq r cs) else FalseF)" |
"inf⇣+ φ = φ"

lemma qfree_min_inf: "qfree φ ==> qfree(inf⇣- φ)"
by(induct φ rule:min_inf.induct) simp_all

lemma qfree_plus_inf: "qfree φ ==> qfree(inf⇣+ φ)"
by(induct φ rule:plus_inf.induct) simp_all

lemma min_inf:
  "nqfree f ==> ∃x. ∀y≤x. R.I (inf⇣- f) xs = R.I f (y # xs)"
  (is "_ ==> ∃x. ?P f x")
proof(induct f)
  case (Atom a)
  show ?case
  proof (cases a)
    case (Less r cs)
    show ?thesis
    proof(cases cs)
      case Nil thus ?thesis using Less by simp
    next
      case (Cons c cs)
      { assume "c=0" hence ?thesis using Less Cons by simp }
      moreover
      { assume "c<0"
        hence "?P (Atom a) ((r - ⟨cs,xs⟩ + 1)/c)" using Less Cons
          by(auto simp add: field_simps)
        hence ?thesis .. }
      moreover
      { assume "c>0"
        hence "?P (Atom a) ((r - ⟨cs,xs⟩ - 1)/c)" using Less Cons
          by(auto simp add: field_simps)
        hence ?thesis .. }
      ultimately show ?thesis by force
    qed
  next
    case (Eq r cs)
    show ?thesis
    proof(cases cs)
      case Nil thus ?thesis using Eq by simp
    next
      case (Cons c cs)
      { assume "c=0" hence ?thesis using Eq Cons by simp }
      moreover
      { assume "c<0"
        hence "?P (Atom a) ((r - ⟨cs,xs⟩ + 1)/c)" using Eq Cons
          by(auto simp add: field_simps)
        hence ?thesis .. }
      moreover
      { assume "c>0"
        hence "?P (Atom a) ((r - ⟨cs,xs⟩ - 1)/c)" using Eq Cons
          by(auto simp add: field_simps)
        hence ?thesis .. }
      ultimately show ?thesis by force
    qed
  qed
next
  case (And f1 f2)
  then obtain x1 x2 where "?P f1 x1" "?P f2 x2" by fastforce+
  hence "?P (And f1 f2) (min x1 x2)" by(force simp:and_def)
  thus ?case ..
next
  case (Or f1 f2)
  then obtain x1 x2 where "?P f1 x1" "?P f2 x2" by fastforce+
  hence "?P (Or f1 f2) (min x1 x2)" by(force simp:or_def)
  thus ?case ..
qed simp_all

lemma plus_inf:
  "nqfree f ==> ∃x. ∀y≥x. R.I (inf⇣+ f) xs = R.I f (y # xs)"
  (is "_ ==> ∃x. ?P f x")
proof(induct f)
  case (Atom a)
  show ?case
  proof (cases a)
    case (Less r cs)
    show ?thesis
    proof(cases cs)
      case Nil thus ?thesis using Less by simp
    next
      case (Cons c cs)
      { assume "c=0" hence ?thesis using Less Cons by simp }
      moreover
      { assume "c<0"
        hence "?P (Atom a) ((r - ⟨cs,xs⟩ - 1)/c)" using Less Cons
          by(auto simp add: field_simps)
        hence ?thesis .. }
      moreover
      { assume "c>0"
        hence "?P (Atom a) ((r - ⟨cs,xs⟩ + 1)/c)" using Less Cons
          by(auto simp add: field_simps)
        hence ?thesis .. }
      ultimately show ?thesis by force
    qed
  next
    case (Eq r cs)
    show ?thesis
    proof(cases cs)
      case Nil thus ?thesis using Eq by simp
    next
      case (Cons c cs)
      { assume "c=0" hence ?thesis using Eq Cons by simp }
      moreover
      { assume "c<0"
        hence "?P (Atom a) ((r - ⟨cs,xs⟩ - 1)/c)" using Eq Cons
          by(auto simp add: field_simps)
        hence ?thesis .. }
      moreover
      { assume "c>0"
        hence "?P (Atom a) ((r - ⟨cs,xs⟩ + 1)/c)" using Eq Cons
          by(auto simp add: field_simps)
        hence ?thesis .. }
      ultimately show ?thesis by force
    qed
  qed
next
  case (And f1 f2)
  then obtain x1 x2 where "?P f1 x1" "?P f2 x2" by fastforce+
  hence "?P (And f1 f2) (max x1 x2)" by(force simp:and_def)
  thus ?case ..
next
  case (Or f1 f2)
  then obtain x1 x2 where "?P f1 x1" "?P f2 x2" by fastforce+
  hence "?P (Or f1 f2) (max x1 x2)" by(force simp:or_def)
  thus ?case ..
qed simp_all


declare [[simp_depth_limit = 4]]

lemma LBex:
 "[| nqfree f; R.I f (x#xs); ¬R.I (inf⇣- f) xs; x ∉ EQ f xs |]
  ==> ∃l∈ LB f xs. l < x"
apply(induct f)
apply simp
apply simp
apply (case_tac a)
apply(auto simp add: depends⇣R_def field_simps split:if_splits list.splits)
apply fastforce+
done


lemma UBex:
 "[| nqfree f; R.I f (x#xs); ¬R.I (inf⇣+ f) xs; x ∉ EQ f xs |]
  ==> ∃u∈ UB f xs. x < u"
apply(induct f)
apply simp
apply simp
apply(case_tac a)
apply(auto simp add: depends⇣R_def field_simps split:if_splits list.splits)
apply fastforce+
done

declare [[simp_depth_limit = 50]]

lemma finite_LB: "finite(LB f xs)"
proof -
  have "LB f xs = (λ(r,cs). r + ⟨cs,xs⟩) ` set(lbounds(R.atoms⇣0 f))"
    by (force simp:set_lbounds image_def field_simps)
  thus ?thesis by simp
qed

lemma finite_UB: "finite(UB f xs)"
proof -
  have "UB f xs = (λ(r,cs). r + ⟨cs,xs⟩) ` set(ubounds(R.atoms⇣0 f))"
    by (force simp:set_ubounds image_def field_simps)
  thus ?thesis by simp
qed


end