Theory LinArith

(*  Author:     Tobias Nipkow, 2007  *)

section‹Linear real arithmetic›

theory LinArith
imports QE "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 ∘ 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
  rewrites "ATOM.atoms⇩0 depends⇩R = atoms⇩0"
proof goal_cases
  case 1
  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 2
  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) ← as, c>0]"
definition "ubounds as = [(r/c, (-1/c) *⇩s cs). Less r (c#cs) ← as, c<0]"
definition "ebounds as = [(r/c, (-1/c) *⇩s cs). Eq r (c#cs) ← 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⇩f⇩m (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

context notes [[simp_depth_limit = 4]]
begin

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(rename_tac a)
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(rename_tac a)
apply(case_tac a)
apply(auto simp add: depends⇩R_def field_simps split:if_splits list.splits)
apply fastforce+
done

end

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