Quelle Float.thy Sprache: Isabelle

(*  Title:      HOL/Library/Float.thy
    Author:     Johannes Hölzl, Fabian Immler
    Copyright   2012  TU München
*)

section \<open>Floating-Point Numbers\<close>

theory Float
imports Log_Nat Lattice_Algebras
begin

definition "float = {m * 2 powr e | (m :: int) (e :: int). True}"

typedef float = float
  morphisms real_of_float float_of
  unfolding float_def by auto

setup_lifting type_definition_float

declare real_of_float [code_unfold]

lemmas float_of_inject[simp]

declare [[coercion "real_of_float :: float \ real"]]

lemma real_of_float_eq: "f1 = f2 \ real_of_float f1 = real_of_float f2" for f1 f2 :: float
  unfolding real_of_float_inject ..

declare real_of_float_inverse[simp] float_of_inverse [simp]
declare real_of_float [simp]

subsection \<open>Real operations preserving the representation as floating point number\<close>

lemma floatI: "m * 2 powr e = x \ x \ float" for m e :: int
  by (auto simp: float_def)

lemma zero_float[simp]: "0 \ float"
  by (auto simp: float_def)

lemma one_float[simp]: "1 \ float"
  by (intro floatI[of 1 0]) simp

lemma numeral_float[simp]: "numeral i \ float"
  by (intro floatI[of "numeral i" 0]) simp

lemma neg_numeral_float[simp]: "- numeral i \ float"
  by (intro floatI[of "- numeral i" 0]) simp

lemma real_of_int_float[simp]: "real_of_int x \ float" for x :: int
  by (intro floatI[of x 0]) simp

lemma real_of_nat_float[simp]: "real x \ float" for x :: nat
  by (intro floatI[of x 0]) simp

lemma two_powr_int_float[simp]: "2 powr (real_of_int i) \ float" for i :: int
  by (intro floatI[of 1 i]) simp

lemma two_powr_nat_float[simp]: "2 powr (real i) \ float" for i :: nat
  by (intro floatI[of 1 i]) simp

lemma two_powr_minus_int_float[simp]: "2 powr - (real_of_int i) \ float" for i :: int
  by (intro floatI[of 1 "-i"]) simp

lemma two_powr_minus_nat_float[simp]: "2 powr - (real i) \ float" for i :: nat
  by (intro floatI[of 1 "-i"]) simp

lemma two_powr_numeral_float[simp]: "2 powr numeral i \ float"
  by (intro floatI[of 1 "numeral i"]) simp

lemma two_powr_neg_numeral_float[simp]: "2 powr - numeral i \ float"
  by (intro floatI[of 1 "- numeral i"]) simp

lemma two_pow_float[simp]: "2 ^ n \ float"
  by (intro floatI[of 1 n]) (simp add: powr_realpow)

lemma plus_float[simp]: "r \ float \ p \ float \ r + p \ float"
  unfolding float_def
proof (safe, simp)
  have *: "\(m::int) (e::int). m1 * 2 powr e1 + m2 * 2 powr e2 = m * 2 powr e"
    if "e1 \ e2" for e1 m1 e2 m2 :: int
  proof -
    from that have "m1 * 2 powr e1 + m2 * 2 powr e2 = (m1 + m2 * 2 ^ nat (e2 - e1)) * 2 powr e1"
      by (simp add: powr_diff field_simps flip: powr_realpow)
    then show ?thesis
      by blast
  qed
  fix e1 m1 e2 m2 :: int
  consider "e2 \ e1" | "e1 \ e2" by (rule linorder_le_cases)
  then show "\(m::int) (e::int). m1 * 2 powr e1 + m2 * 2 powr e2 = m * 2 powr e"
  proof cases
    case 1
    from *[OF this, of m2 m1] show ?thesis
      by (simp add: ac_simps)
  next
    case 2
    then show ?thesis by (rule *)
  qed
qed

lemma uminus_float[simp]: "x \ float \ -x \ float"
  by (simp add: float_def) (metis mult_minus_left of_int_minus)

lemma times_float[simp]: "x \ float \ y \ float \ x * y \ float"
  apply (clarsimp simp: float_def mult_ac)
  by (metis mult.assoc of_int_mult of_int_add powr_add)

lemma minus_float[simp]: "x \ float \ y \ float \ x - y \ float"
  using plus_float [of x "- y"] by simp

lemma abs_float[simp]: "x \ float \ \x\ \ float"
  by (cases x rule: linorder_cases[of 0]) auto

lemma sgn_of_float[simp]: "x \ float \ sgn x \ float"
  by (simp add: sgn_real_def)

lemma div_power_2_float[simp]: "x \ float \ x / 2^d \ float"
  by (simp add: float_def) (metis of_int_diff of_int_of_nat_eq powr_diff powr_realpow zero_less_numeral times_divide_eq_right)

lemma div_power_2_int_float[simp]: "x \ float \ x / (2::int)^d \ float"
  by simp

lemma div_numeral_Bit0_float[simp]:
  assumes "x / numeral n \ float"
  shows "x / (numeral (Num.Bit0 n)) \ float"
proof -
  have "(x / numeral n) / 2^1 \ float"
    by (intro assms div_power_2_float)
  also have "(x / numeral n) / 2^1 = x / (numeral (Num.Bit0 n))"
    by (induct n) auto
  finally show ?thesis .
qed

lemma div_neg_numeral_Bit0_float[simp]:
  assumes "x / numeral n \ float"
  shows "x / (- numeral (Num.Bit0 n)) \ float"
  using assms by force

lemma power_float[simp]:
  assumes "a \ float"
  shows "a ^ b \ float"
proof -
  from assms obtain m e :: int where "a = m * 2 powr e"
    by (auto simp: float_def)
  then show ?thesis
    by (intro floatI[where m="m^b" and e = "e*b"])
       (auto simp: powr_powr power_mult_distrib simp flip: powr_realpow)
qed

lift_definition Float :: "int \ int \ float" is "\(m::int) (e::int). m * 2 powr e"
  by simp
declare Float.rep_eq[simp]

code_datatype Float

lemma compute_real_of_float[code]:
  "real_of_float (Float m e) = (if e \ 0 then m * 2 ^ nat e else m / 2 ^ (nat (-e)))"
  by (simp add: powr_int)

subsection \<open>Arithmetic operations on floating point numbers\<close>

instantiation float :: "{ring_1,linorder,linordered_ring,linordered_idom,numeral,equal}"
begin

lift_definition zero_float :: float is 0 by simp
declare zero_float.rep_eq[simp]

lift_definition one_float :: float is 1 by simp
declare one_float.rep_eq[simp]

lift_definition plus_float :: "float \ float \ float" is "(+)" by simp
declare plus_float.rep_eq[simp]

lift_definition times_float :: "float \ float \ float" is "(*)" by simp
declare times_float.rep_eq[simp]

lift_definition minus_float :: "float \ float \ float" is "(-)" by simp
declare minus_float.rep_eq[simp]

lift_definition uminus_float :: "float \ float" is "uminus" by simp
declare uminus_float.rep_eq[simp]

lift_definition abs_float :: "float \ float" is abs by simp
declare abs_float.rep_eq[simp]

lift_definition sgn_float :: "float \ float" is sgn by simp
declare sgn_float.rep_eq[simp]

lift_definition equal_float :: "float \ float \ bool" is "(=) :: real \ real \ bool" .

lift_definition less_eq_float :: "float \ float \ bool" is "(\)" .
declare less_eq_float.rep_eq[simp]

lift_definition less_float :: "float \ float \ bool" is "(<)" .
declare less_float.rep_eq[simp]

instance
  by standard (transfer; fastforce simp add: field_simps intro: mult_left_mono mult_right_mono)+

end

lemma real_of_float [simp]: "real_of_float (of_nat n) = of_nat n"
  by (induct n) simp_all

lemma real_of_float_of_int_eq [simp]: "real_of_float (of_int z) = of_int z"
  by (cases z rule: int_diff_cases) (simp_all add: of_rat_diff)

lemma Float_0_eq_0[simp]: "Float 0 e = 0"
  by transfer simp

lemma real_of_float_power[simp]: "real_of_float (f^n) = real_of_float f^n" for f :: float
  by (induct n) simp_all

lemma real_of_float_min: "real_of_float (min x y) = min (real_of_float x) (real_of_float y)"
  and real_of_float_max: "real_of_float (max x y) = max (real_of_float x) (real_of_float y)"
  for x y :: float
  by (simp_all add: min_def max_def)

instance float :: unbounded_dense_linorder
proof
  fix a b :: float
  show "\c. a < c"
    by (metis Float.real_of_float less_float.rep_eq reals_Archimedean2)
  show "\c. c < a"
    by (metis add_0 add_strict_right_mono neg_less_0_iff_less zero_less_one)
  show "\c. a < c \ c < b" if "a < b"
    apply (rule exI[of _ "(a + b) * Float 1 (- 1)"])
    using that
    apply transfer
    apply (simp add: powr_minus)
    done
qed

instantiation float :: lattice_ab_group_add
begin

definition inf_float :: "float \ float \ float"
  where "inf_float a b = min a b"

definition sup_float :: "float \ float \ float"
  where "sup_float a b = max a b"

instance
  by standard (transfer; simp add: inf_float_def sup_float_def real_of_float_min real_of_float_max)+

end

lemma float_numeral[simp]: "real_of_float (numeral x :: float) = numeral x"
  by (metis of_int_numeral real_of_float_of_int_eq)

lemma transfer_numeral [transfer_rule]:
  "rel_fun (=) pcr_float (numeral :: _ \ real) (numeral :: _ \ float)"
  by (simp add: rel_fun_def float.pcr_cr_eq cr_float_def)

lemma float_neg_numeral[simp]: "real_of_float (- numeral x :: float) = - numeral x"
  by simp

lemma transfer_neg_numeral [transfer_rule]:
  "rel_fun (=) pcr_float (- numeral :: _ \ real) (- numeral :: _ \ float)"
  by (simp add: rel_fun_def float.pcr_cr_eq cr_float_def)

lemma float_of_numeral: "numeral k = float_of (numeral k)"
  and float_of_neg_numeral: "- numeral k = float_of (- numeral k)"
  unfolding real_of_float_eq by simp_all

subsection \<open>Quickcheck\<close>

instantiation float :: exhaustive
begin

definition exhaustive_float where
  "exhaustive_float f d =
    Quickcheck_Exhaustive.exhaustive (\<lambda>x. Quickcheck_Exhaustive.exhaustive (\<lambda>y. f (Float x y)) d) d"

instance ..

end

context
  includes term_syntax
begin

definition [code_unfold]:
  "valtermify_float x y = Code_Evaluation.valtermify Float {\} x {\} y"

end

instantiation float :: full_exhaustive
begin

definition
  "full_exhaustive_float f d =
    Quickcheck_Exhaustive.full_exhaustive
      (\<lambda>x. Quickcheck_Exhaustive.full_exhaustive (\<lambda>y. f (valtermify_float x y)) d) d"

instance ..

end

instantiation float :: random
begin

definition "Quickcheck_Random.random i =
  scomp (Quickcheck_Random.random (2 ^ nat_of_natural i))
    (\<lambda>man. scomp (Quickcheck_Random.random i) (\<lambda>exp. Pair (valtermify_float man exp)))"

instance ..

end

subsection \<open>Represent floats as unique mantissa and exponent\<close>

lemma int_induct_abs[case_names less]:
  fixes j :: int
  assumes H: "\n. (\i. \i\ < \n\ \ P i) \ P n"
  shows "P j"
proof (induct "nat \j\" arbitrary: j rule: less_induct)
  case less
  show ?case by (rule H[OF less]) simp
qed

lemma int_cancel_factors:
  fixes n :: int
  assumes "1 < r"
  shows "n = 0 \ (\k i. n = k * r ^ i \ \ r dvd k)"
proof (induct n rule: int_induct_abs)
  case (less n)
  have "\k i. n = k * r ^ Suc i \ \ r dvd k" if "n \ 0" "n = m * r" for m
  proof -
    from that have "\m \ < \n\"
      using \<open>1 < r\<close> by (simp add: abs_mult)
    from less[OF this] that show ?thesis by auto
  qed
  then show ?case
    by (metis dvd_def monoid_mult_class.mult.right_neutral mult.commute power_0)
qed

lemma mult_powr_eq_mult_powr_iff_asym:
  fixes m1 m2 e1 e2 :: int
  assumes m1: "\ 2 dvd m1"
    and "e1 \ e2"
  shows "m1 * 2 powr e1 = m2 * 2 powr e2 \ m1 = m2 \ e1 = e2"
  (is "?lhs \ ?rhs")
proof
  show ?rhs if eq: ?lhs
  proof -
    have "m1 \ 0"
      using m1 unfolding dvd_def by auto
    from \<open>e1 \<le> e2\<close> eq have "m1 = m2 * 2 powr nat (e2 - e1)"
      by (simp add: powr_diff field_simps)
    also have "\ = m2 * 2^nat (e2 - e1)"
      by (simp add: powr_realpow)
    finally have m1_eq: "m1 = m2 * 2^nat (e2 - e1)"
      by linarith
    with m1 have "m1 = m2"
      by (cases "nat (e2 - e1)") (auto simp add: dvd_def)
    then show ?thesis
      using eq \<open>m1 \<noteq> 0\<close> by (simp add: powr_inj)
  qed
  show ?lhs if ?rhs
    using that by simp
qed

lemma mult_powr_eq_mult_powr_iff:
  "\ 2 dvd m1 \ \ 2 dvd m2 \ m1 * 2 powr e1 = m2 * 2 powr e2 \ m1 = m2 \ e1 = e2"
  for m1 m2 e1 e2 :: int
  using mult_powr_eq_mult_powr_iff_asym[of m1 e1 e2 m2]
  using mult_powr_eq_mult_powr_iff_asym[of m2 e2 e1 m1]
  by (cases e1 e2 rule: linorder_le_cases) auto

lemma floatE_normed:
  assumes x: "x \ float"
  obtains (zero) "x = 0"
   | (powr) m e :: int where "x = m * 2 powr e" "\ 2 dvd m" "x \ 0"
proof -
  have "\(m::int) (e::int). x = m * 2 powr e \ \ (2::int) dvd m" if "x \ 0"
  proof -
    from x obtain m e :: int where x: "x = m * 2 powr e"
      by (auto simp: float_def)
    with \<open>x \<noteq> 0\<close> int_cancel_factors[of 2 m] obtain k i where "m = k * 2 ^ i" "\<not> 2 dvd k"
      by auto
    with \<open>\<not> 2 dvd k\<close> x have "x = real_of_int k * 2 powr real_of_int (e + int i) \<and> odd k"
      by (simp add: powr_add powr_realpow)
    then show ?thesis
      by blast
  qed
  with that show thesis by blast
qed

lemma float_normed_cases:
  fixes f :: float
  obtains (zero) "f = 0"
   | (powr) m e :: int where "real_of_float f = m * 2 powr e" "\ 2 dvd m" "f \ 0"
proof (atomize_elim, induct f)
  case (float_of y)
  then show ?case
    by (cases rule: floatE_normed) (auto simp: zero_float_def)
qed

definition mantissa :: "float \ int"
  where "mantissa f =
    fst (SOME p::int \<times> int. (f = 0 \<and> fst p = 0 \<and> snd p = 0) \<or>
      (f \<noteq> 0 \<and> real_of_float f = real_of_int (fst p) * 2 powr real_of_int (snd p) \<and> \<not> 2 dvd fst p))"

definition exponent :: "float \ int"
  where "exponent f =
    snd (SOME p::int \<times> int. (f = 0 \<and> fst p = 0 \<and> snd p = 0) \<or>
      (f \<noteq> 0 \<and> real_of_float f = real_of_int (fst p) * 2 powr real_of_int (snd p) \<and> \<not> 2 dvd fst p))"

lemma exponent_0[simp]: "exponent 0 = 0" (is ?E)
  and mantissa_0[simp]: "mantissa 0 = 0" (is ?M)
proof -
  have "\p::int \ int. fst p = 0 \ snd p = 0 \ p = (0, 0)"
    by auto
  then show ?E ?M
    by (auto simp add: mantissa_def exponent_def zero_float_def)
qed

lemma mantissa_exponent: "real_of_float f = mantissa f * 2 powr exponent f" (is ?E)
  and mantissa_not_dvd: "f \ 0 \ \ 2 dvd mantissa f" (is "_ \ ?D")
proof cases
  assume [simp]: "f \ 0"
  have "f = mantissa f * 2 powr exponent f \ \ 2 dvd mantissa f"
  proof (cases f rule: float_normed_cases)
    case zero
    then show ?thesis by simp
  next
    case (powr m e)
    then have "\p::int \ int. (f = 0 \ fst p = 0 \ snd p = 0) \
      (f \<noteq> 0 \<and> real_of_float f = real_of_int (fst p) * 2 powr real_of_int (snd p) \<and> \<not> 2 dvd fst p)"
      by auto
    then show ?thesis
      unfolding exponent_def mantissa_def
      by (rule someI2_ex) simp
  qed
  then show ?E ?D by auto
qed simp

lemma mantissa_noteq_0: "f \ 0 \ mantissa f \ 0"
  using mantissa_not_dvd[of f] by auto

lemma mantissa_eq_zero_iff: "mantissa x = 0 \ x = 0"
  (is "?lhs \ ?rhs")
proof
  show ?rhs if ?lhs
  proof -
    from that have z: "0 = real_of_float x"
      using mantissa_exponent by simp
    show ?thesis
      by (simp add: zero_float_def z)
  qed
  show ?lhs if ?rhs
    using that by simp
qed

lemma mantissa_pos_iff: "0 < mantissa x \ 0 < x"
  by (auto simp: mantissa_exponent algebra_split_simps)

lemma mantissa_nonneg_iff: "0 \ mantissa x \ 0 \ x"
  by (auto simp: mantissa_exponent algebra_split_simps)

lemma mantissa_neg_iff: "0 > mantissa x \ 0 > x"
  by (auto simp: mantissa_exponent algebra_split_simps)

lemma
  fixes m e :: int
  defines "f \ float_of (m * 2 powr e)"
  assumes dvd: "\ 2 dvd m"
  shows mantissa_float: "mantissa f = m" (is "?M")
    and exponent_float: "m \ 0 \ exponent f = e" (is "_ \ ?E")
proof cases
  assume "m = 0"
  with dvd show "mantissa f = m" by auto
next
  assume "m \ 0"
  then have f_not_0: "f \ 0" by (simp add: f_def zero_float_def)
  from mantissa_exponent[of f] have "m * 2 powr e = mantissa f * 2 powr exponent f"
    by (auto simp add: f_def)
  then show ?M ?E
    using mantissa_not_dvd[OF f_not_0] dvd
    by (auto simp: mult_powr_eq_mult_powr_iff)
qed

subsection \<open>Compute arithmetic operations\<close>

lemma Float_mantissa_exponent: "Float (mantissa f) (exponent f) = f"
  unfolding real_of_float_eq mantissa_exponent[of f] by simp

lemma Float_cases [cases type: float]:
  fixes f :: float
  obtains (Float) m e :: int where "f = Float m e"
  using Float_mantissa_exponent[symmetric]
  by (atomize_elim) auto

lemma denormalize_shift:
  assumes f_def: "f = Float m e"
    and not_0: "f \ 0"
  obtains i where "m = mantissa f * 2 ^ i" "e = exponent f - i"
proof
  from mantissa_exponent[of f] f_def
  have "m * 2 powr e = mantissa f * 2 powr exponent f"
    by simp
  then have eq: "m = mantissa f * 2 powr (exponent f - e)"
    by (simp add: powr_diff field_simps)
  moreover
  have "e \ exponent f"
  proof (rule ccontr)
    assume "\ e \ exponent f"
    then have pos: "exponent f < e" by simp
    then have "2 powr (exponent f - e) = 2 powr - real_of_int (e - exponent f)"
      by simp
    also have "\ = 1 / 2^nat (e - exponent f)"
      using pos by (simp flip: powr_realpow add: powr_diff)
    finally have "m * 2^nat (e - exponent f) = real_of_int (mantissa f)"
      using eq by simp
    then have "mantissa f = m * 2^nat (e - exponent f)"
      by linarith
    with \<open>exponent f < e\<close> have "2 dvd mantissa f"
      by (force intro: dvdI[where k="m * 2^(nat (e-exponent f)) div 2"])
    then show False using mantissa_not_dvd[OF not_0] by simp
  qed
  ultimately have "real_of_int m = mantissa f * 2^nat (exponent f - e)"
    by (simp flip: powr_realpow)
  with \<open>e \<le> exponent f\<close>
  show "m = mantissa f * 2 ^ nat (exponent f - e)"
    by linarith
  show "e = exponent f - nat (exponent f - e)"
    using \<open>e \<le> exponent f\<close> by auto
qed

context
begin

qualified lemma compute_float_zero[code_unfold, code]: "0 = Float 0 0"
  by transfer simp

qualified lemma compute_float_one[code_unfold, code]: "1 = Float 1 0"
  by transfer simp

lift_definition normfloat :: "float \ float" is "\x. x" .
lemma normloat_id[simp]: "normfloat x = x" by transfer rule

qualified lemma compute_normfloat[code]:
  "normfloat (Float m e) =
    (if m mod 2 = 0 \<and> m \<noteq> 0 then normfloat (Float (m div 2) (e + 1))
     else if m = 0 then 0 else Float m e)"
  by transfer (auto simp add: powr_add zmod_eq_0_iff)

qualified lemma compute_float_numeral[code_abbrev]: "Float (numeral k) 0 = numeral k"
  by transfer simp

qualified lemma compute_float_neg_numeral[code_abbrev]: "Float (- numeral k) 0 = - numeral k"
  by transfer simp

qualified lemma compute_float_uminus[code]: "- Float m1 e1 = Float (- m1) e1"
  by transfer simp

qualified lemma compute_float_times[code]: "Float m1 e1 * Float m2 e2 = Float (m1 * m2) (e1 + e2)"
  by transfer (simp add: field_simps powr_add)

qualified lemma compute_float_plus[code]:
  "Float m1 e1 + Float m2 e2 =
    (if m1 = 0 then Float m2 e2
     else if m2 = 0 then Float m1 e1
     else if e1 \<le> e2 then Float (m1 + m2 * 2^nat (e2 - e1)) e1
     else Float (m2 + m1 * 2^nat (e1 - e2)) e2)"
  by transfer (simp add: field_simps powr_realpow[symmetric] powr_diff)

qualified lemma compute_float_minus[code]: "f - g = f + (-g)" for f g :: float
  by simp

qualified lemma compute_float_sgn[code]:
  "sgn (Float m1 e1) = (if 0 < m1 then 1 else if m1 < 0 then -1 else 0)"
  by transfer (simp add: sgn_mult)

lift_definition is_float_pos :: "float \ bool" is "(<) 0 :: real \ bool" .

qualified lemma compute_is_float_pos[code]: "is_float_pos (Float m e) \ 0 < m"
  by transfer (auto simp add: zero_less_mult_iff not_le[symmetric, of _ 0])

lift_definition is_float_nonneg :: "float \ bool" is "(\) 0 :: real \ bool" .

qualified lemma compute_is_float_nonneg[code]: "is_float_nonneg (Float m e) \ 0 \ m"
  by transfer (auto simp add: zero_le_mult_iff not_less[symmetric, of _ 0])

lift_definition is_float_zero :: "float \ bool" is "(=) 0 :: real \ bool" .

qualified lemma compute_is_float_zero[code]: "is_float_zero (Float m e) \ 0 = m"
  by transfer (auto simp add: is_float_zero_def)

qualified lemma compute_float_abs[code]: "\Float m e\ = Float \m\ e"
  by transfer (simp add: abs_mult)

qualified lemma compute_float_eq[code]: "equal_class.equal f g = is_float_zero (f - g)"
  by transfer simp

end

subsection \<open>Lemmas for types \<^typ>\<open>real\<close>, \<^typ>\<open>nat\<close>, \<^typ>\<open>int\<close>\<close>

lemmas real_of_ints =
  of_int_add
  of_int_minus
  of_int_diff
  of_int_mult
  of_int_power
  of_int_numeral of_int_neg_numeral

lemmas int_of_reals = real_of_ints[symmetric]

subsection \<open>Rounding Real Numbers\<close>

definition round_down :: "int \ real \ real"
  where "round_down prec x = \x * 2 powr prec\ * 2 powr -prec"

definition round_up :: "int \ real \ real"
  where "round_up prec x = \x * 2 powr prec\ * 2 powr -prec"

lemma round_down_float[simp]: "round_down prec x \ float"
  unfolding round_down_def
  by (auto intro!: times_float simp flip: of_int_minus)

lemma round_up_float[simp]: "round_up prec x \ float"
  unfolding round_up_def
  by (auto intro!: times_float simp flip: of_int_minus)

lemma round_up: "x \ round_up prec x"
  by (simp add: powr_minus_divide le_divide_eq round_up_def ceiling_correct)

lemma round_down: "round_down prec x \ x"
  by (simp add: powr_minus_divide divide_le_eq round_down_def)

lemma round_up_0[simp]: "round_up p 0 = 0"
  unfolding round_up_def by simp

lemma round_down_0[simp]: "round_down p 0 = 0"
  unfolding round_down_def by simp

lemma round_up_diff_round_down: "round_up prec x - round_down prec x \ 2 powr -prec"
proof -
  have "round_up prec x - round_down prec x = (\x * 2 powr prec\ - \x * 2 powr prec\) * 2 powr -prec"
    by (simp add: round_up_def round_down_def field_simps)
  also have "\ \ 1 * 2 powr -prec"
    by (rule mult_mono)
      (auto simp flip: of_int_diff simp: ceiling_diff_floor_le_1)
  finally show ?thesis by simp
qed

lemma round_down_shift: "round_down p (x * 2 powr k) = 2 powr k * round_down (p + k) x"
  unfolding round_down_def
  by (simp add: powr_add powr_mult field_simps powr_diff)
    (simp flip: powr_add)

lemma round_up_shift: "round_up p (x * 2 powr k) = 2 powr k * round_up (p + k) x"
  unfolding round_up_def
  by (simp add: powr_add powr_mult field_simps powr_diff)
    (simp flip: powr_add)

lemma round_up_uminus_eq: "round_up p (-x) = - round_down p x"
  and round_down_uminus_eq: "round_down p (-x) = - round_up p x"
  by (auto simp: round_up_def round_down_def ceiling_def)

lemma round_up_mono: "x \ y \ round_up p x \ round_up p y"
  by (auto intro!: ceiling_mono simp: round_up_def)

lemma round_up_le1:
  assumes "x \ 1" "prec \ 0"
  shows "round_up prec x \ 1"
proof -
  have "real_of_int \x * 2 powr prec\ \ real_of_int \2 powr real_of_int prec\"
    using assms by (auto intro!: ceiling_mono)
  also have "\ = 2 powr prec" using assms by (auto simp: powr_int intro!: exI[where x="2^nat prec"])
  finally show ?thesis
    by (simp add: round_up_def) (simp add: powr_minus inverse_eq_divide)
qed

lemma round_up_less1:
  assumes "x < 1 / 2" "p > 0"
  shows "round_up p x < 1"
proof -
  have "x * 2 powr p < 1 / 2 * 2 powr p"
    using assms by simp
  also have "\ \ 2 powr p - 1" using \p > 0\
    by (auto simp: powr_diff powr_int field_simps self_le_power)
  finally show ?thesis using \<open>p > 0\<close>
    by (simp add: round_up_def field_simps powr_minus powr_int ceiling_less_iff)
qed

lemma round_down_ge1:
  assumes x: "x \ 1"
  assumes prec: "p \ - log 2 x"
  shows "1 \ round_down p x"
proof cases
  assume nonneg: "0 \ p"
  have "2 powr p = real_of_int \2 powr real_of_int p\"
    using nonneg by (auto simp: powr_int)
  also have "\ \ real_of_int \x * 2 powr p\"
    using assms by (auto intro!: floor_mono)
  finally show ?thesis
    by (simp add: round_down_def) (simp add: powr_minus inverse_eq_divide)
next
  assume neg: "\ 0 \ p"
  have "x = 2 powr (log 2 x)"
    using x by simp
  also have "2 powr (log 2 x) \ 2 powr - p"
    using prec by auto
  finally have x_le: "x \ 2 powr -p" .

  from neg have "2 powr real_of_int p \ 2 powr 0"
    by (intro powr_mono) auto
  also have "\ \ \2 powr 0::real\" by simp
  also have "\ \ \x * 2 powr (real_of_int p)\"
    unfolding of_int_le_iff
    using x x_le by (intro floor_mono) (simp add: powr_minus_divide field_simps)
  finally show ?thesis
    using prec x
    by (simp add: round_down_def powr_minus_divide pos_le_divide_eq)
qed

lemma round_up_le0: "x \ 0 \ round_up p x \ 0"
  unfolding round_up_def
  by (auto simp: field_simps mult_le_0_iff zero_le_mult_iff)

subsection \<open>Rounding Floats\<close>

definition div_twopow :: "int \ nat \ int"
  where [simp]: "div_twopow x n = x div (2 ^ n)"

definition mod_twopow :: "int \ nat \ int"
  where [simp]: "mod_twopow x n = x mod (2 ^ n)"

lemma compute_div_twopow[code]:
  "div_twopow x n = (if x = 0 \ x = -1 \ n = 0 then x else div_twopow (x div 2) (n - 1))"
  by (cases n) (auto simp: zdiv_zmult2_eq div_eq_minus1)

lemma compute_mod_twopow[code]:
  "mod_twopow x n = (if n = 0 then 0 else x mod 2 + 2 * mod_twopow (x div 2) (n - 1))"
  by (cases n) (auto simp: zmod_zmult2_eq)

lift_definition float_up :: "int \ float \ float" is round_up by simp
declare float_up.rep_eq[simp]

lemma round_up_correct: "round_up e f - f \ {0..2 powr -e}"
  unfolding atLeastAtMost_iff
proof
  have "round_up e f - f \ round_up e f - round_down e f"
    using round_down by simp
  also have "\ \ 2 powr -e"
    using round_up_diff_round_down by simp
  finally show "round_up e f - f \ 2 powr - (real_of_int e)"
    by simp
qed (simp add: algebra_simps round_up)

lemma float_up_correct: "real_of_float (float_up e f) - real_of_float f \ {0..2 powr -e}"
  by transfer (rule round_up_correct)

lift_definition float_down :: "int \ float \ float" is round_down by simp
declare float_down.rep_eq[simp]

lemma round_down_correct: "f - (round_down e f) \ {0..2 powr -e}"
  unfolding atLeastAtMost_iff
proof
  have "f - round_down e f \ round_up e f - round_down e f"
    using round_up by simp
  also have "\ \ 2 powr -e"
    using round_up_diff_round_down by simp
  finally show "f - round_down e f \ 2 powr - (real_of_int e)"
    by simp
qed (simp add: algebra_simps round_down)

lemma float_down_correct: "real_of_float f - real_of_float (float_down e f) \ {0..2 powr -e}"
  by transfer (rule round_down_correct)

context
begin

qualified lemma compute_float_down[code]:
  "float_down p (Float m e) =
    (if p + e < 0 then Float (div_twopow m (nat (-(p + e)))) (-p) else Float m e)"
proof (cases "p + e < 0")
  case True
  then have "real_of_int ((2::int) ^ nat (-(p + e))) = 2 powr (-(p + e))"
    using powr_realpow[of 2 "nat (-(p + e))"] by simp
  also have "\ = 1 / 2 powr p / 2 powr e"
    unfolding powr_minus_divide of_int_minus by (simp add: powr_add)
  finally show ?thesis
    using \<open>p + e < 0\<close>
    apply transfer
    apply (simp add: round_down_def field_simps flip: floor_divide_of_int_eq powr_add)
    apply (metis (no_types, opaque_lifting) Float.rep_eq
      add.inverse_inverse compute_real_of_float diff_minus_eq_add
      floor_divide_of_int_eq int_of_reals(1) linorder_not_le
      minus_add_distrib of_int_eq_numeral_power_cancel_iff)
    done
next
  case False
  then have r: "real_of_int e + real_of_int p = real (nat (e + p))"
    by simp
  have r: "\(m * 2 powr e) * 2 powr real_of_int p\ = (m * 2 powr e) * 2 powr real_of_int p"
    by (auto intro: exI[where x="m*2^nat (e+p)"]
        simp add: ac_simps powr_add[symmetric] r powr_realpow)
  with \<open>\<not> p + e < 0\<close> show ?thesis
    by transfer (auto simp add: round_down_def field_simps powr_add powr_minus)
qed

lemma abs_round_down_le: "\f - (round_down e f)\ \ 2 powr -e"
  using round_down_correct[of f e] by simp

lemma abs_round_up_le: "\f - (round_up e f)\ \ 2 powr -e"
  using round_up_correct[of e f] by simp

lemma round_down_nonneg: "0 \ s \ 0 \ round_down p s"
  by (auto simp: round_down_def)

lemma ceil_divide_floor_conv:
  assumes "b \ 0"
  shows "\real_of_int a / real_of_int b\ =
    (if b dvd a then a div b else \<lfloor>real_of_int a / real_of_int b\<rfloor> + 1)"
proof (cases "b dvd a")
  case True
  then show ?thesis
    by (simp add: ceiling_def floor_divide_of_int_eq dvd_neg_div
             flip: of_int_minus divide_minus_left)
next
  case False
  then have "a mod b \ 0"
    by auto
  then have ne: "real_of_int (a mod b) / real_of_int b \ 0"
    using \<open>b \<noteq> 0\<close> by auto
  have "\real_of_int a / real_of_int b\ = \real_of_int a / real_of_int b\ + 1"
    by (metis add_cancel_left_right ceiling_altdef floor_divide_of_int_eq ne of_int_div_aux)
  then show ?thesis
    using \<open>\<not> b dvd a\<close> by simp
qed

qualified lemma compute_float_up[code]: "float_up p x = - float_down p (-x)"
  by transfer (simp add: round_down_uminus_eq)

end

lemma bitlen_Float:
  fixes m e
  defines [THEN meta_eq_to_obj_eq]: "f \ Float m e"
  shows "bitlen \mantissa f\ + exponent f = (if m = 0 then 0 else bitlen \m\ + e)"
proof (cases "m = 0")
  case True
  then show ?thesis by (simp add: f_def bitlen_alt_def)
next
  case False
  then have "f \ 0"
    unfolding real_of_float_eq by (simp add: f_def)
  then have "mantissa f \ 0"
    by (simp add: mantissa_eq_zero_iff)
  moreover
  obtain i where "m = mantissa f * 2 ^ i" "e = exponent f - int i"
    by (rule f_def[THEN denormalize_shift, OF \<open>f \<noteq> 0\<close>])
  ultimately show ?thesis by (simp add: abs_mult)
qed

lemma float_gt1_scale:
  assumes "1 \ Float m e"
  shows "0 \ e + (bitlen m - 1)"
proof -
  have "0 < Float m e" using assms by auto
  then have "0 < m" using powr_gt_zero[of 2 e]
    by (auto simp: zero_less_mult_iff)
  then have "m \ 0" by auto
  show ?thesis
  proof (cases "0 \ e")
    case True
    then show ?thesis
      using \<open>0 < m\<close> by (simp add: bitlen_alt_def)
  next
    case False
    have "(1::int) < 2" by simp
    let ?S = "2^(nat (-e))"
    have "inverse (2 ^ nat (- e)) = 2 powr e"
      using assms False powr_realpow[of 2 "nat (-e)"]
      by (auto simp: powr_minus field_simps)
    then have "1 \ real_of_int m * inverse ?S"
      using assms False powr_realpow[of 2 "nat (-e)"]
      by (auto simp: powr_minus)
    then have "1 * ?S \ real_of_int m * inverse ?S * ?S"
      by (rule mult_right_mono) auto
    then have "?S \ real_of_int m"
      unfolding mult.assoc by auto
    then have "?S \ m"
      unfolding of_int_le_iff[symmetric] by auto
    from this bitlen_bounds[OF \<open>0 < m\<close>, THEN conjunct2]
    have "nat (-e) < (nat (bitlen m))"
      unfolding power_strict_increasing_iff[OF \<open>1 < 2\<close>, symmetric]
      by (rule order_le_less_trans)
    then have "-e < bitlen m"
      using False by auto
    then show ?thesis
      by auto
  qed
qed

subsection \<open>Truncating Real Numbers\<close>

definition truncate_down::"nat \ real \ real"
  where "truncate_down prec x = round_down (prec - \log 2 \x\\) x"

lemma truncate_down: "truncate_down prec x \ x"
  using round_down by (simp add: truncate_down_def)

lemma truncate_down_le: "x \ y \ truncate_down prec x \ y"
  by (rule order_trans[OF truncate_down])

lemma truncate_down_zero[simp]: "truncate_down prec 0 = 0"
  by (simp add: truncate_down_def)

lemma truncate_down_float[simp]: "truncate_down p x \ float"
  by (auto simp: truncate_down_def)

definition truncate_up::"nat \ real \ real"
  where "truncate_up prec x = round_up (prec - \log 2 \x\\) x"

lemma truncate_up: "x \ truncate_up prec x"
  using round_up by (simp add: truncate_up_def)

lemma truncate_up_le: "x \ y \ x \ truncate_up prec y"
  by (rule order_trans[OF _ truncate_up])

lemma truncate_up_zero[simp]: "truncate_up prec 0 = 0"
  by (simp add: truncate_up_def)

lemma truncate_up_uminus_eq: "truncate_up prec (-x) = - truncate_down prec x"
  and truncate_down_uminus_eq: "truncate_down prec (-x) = - truncate_up prec x"
  by (auto simp: truncate_up_def round_up_def truncate_down_def round_down_def ceiling_def)

lemma truncate_up_float[simp]: "truncate_up p x \ float"
  by (auto simp: truncate_up_def)

lemma mult_powr_eq: "0 < b \ b \ 1 \ 0 < x \ x * b powr y = b powr (y + log b x)"
  by (simp_all add: powr_add)

lemma truncate_down_pos:
  assumes "x > 0"
  shows "truncate_down p x > 0"
proof -
  have "0 \ log 2 x - real_of_int \log 2 x\"
    by (simp add: algebra_simps)
  moreover have "0 \ real p - real_of_int \log 2 x\ + log 2 x"
    by linarith
  ultimately show ?thesis
    using assms
    by (auto simp: truncate_down_def round_down_def mult_powr_eq
      intro!: ge_one_powr_ge_zero mult_pos_pos)
qed

lemma truncate_down_nonneg: "0 \ y \ 0 \ truncate_down prec y"
  by (auto simp: truncate_down_def round_down_def)

lemma truncate_down_ge1: "1 \ x \ 1 \ truncate_down p x"
  apply (auto simp: truncate_down_def algebra_simps intro!: round_down_ge1)
  apply linarith
  done

lemma truncate_up_nonpos: "x \ 0 \ truncate_up prec x \ 0"
  by (auto simp: truncate_up_def round_up_def intro!: mult_nonpos_nonneg)

lemma truncate_up_le1:
  assumes "x \ 1"
  shows "truncate_up p x \ 1"
proof -
  consider "x \ 0" | "x > 0"
    by arith
  then show ?thesis
  proof cases
    case 1
    with truncate_up_nonpos[OF this, of p] show ?thesis
      by simp
  next
    case 2
    then have le: "\log 2 \x\\ \ 0"
      using assms by (auto simp: log_less_iff)
    from assms have "0 \ int p" by simp
    from add_mono[OF this le]
    show ?thesis
      using assms by (simp add: truncate_up_def round_up_le1 add_mono)
  qed
qed

lemma truncate_down_shift_int:
  "truncate_down p (x * 2 powr real_of_int k) = truncate_down p x * 2 powr k"
  by (cases "x = 0")
    (simp_all add: algebra_simps abs_mult log_mult truncate_down_def
      round_down_shift[of _ _ k, simplified])

lemma truncate_down_shift_nat: "truncate_down p (x * 2 powr real k) = truncate_down p x * 2 powr k"
  by (metis of_int_of_nat_eq truncate_down_shift_int)

lemma truncate_up_shift_int: "truncate_up p (x * 2 powr real_of_int k) = truncate_up p x * 2 powr k"
  by (cases "x = 0")
    (simp_all add: algebra_simps abs_mult log_mult truncate_up_def
      round_up_shift[of _ _ k, simplified])

lemma truncate_up_shift_nat: "truncate_up p (x * 2 powr real k) = truncate_up p x * 2 powr k"
  by (metis of_int_of_nat_eq truncate_up_shift_int)

subsection \<open>Truncating Floats\<close>

lift_definition float_round_up :: "nat \ float \ float" is truncate_up
  by (simp add: truncate_up_def)

lemma float_round_up: "real_of_float x \ real_of_float (float_round_up prec x)"
  using truncate_up by transfer simp

lemma float_round_up_zero[simp]: "float_round_up prec 0 = 0"
  by transfer simp

lift_definition float_round_down :: "nat \ float \ float" is truncate_down
  by (simp add: truncate_down_def)

lemma float_round_down: "real_of_float (float_round_down prec x) \ real_of_float x"
  using truncate_down by transfer simp

lemma float_round_down_zero[simp]: "float_round_down prec 0 = 0"
  by transfer simp

lemmas float_round_up_le = order_trans[OF _ float_round_up]
  and float_round_down_le = order_trans[OF float_round_down]

lemma minus_float_round_up_eq: "- float_round_up prec x = float_round_down prec (- x)"
  and minus_float_round_down_eq: "- float_round_down prec x = float_round_up prec (- x)"
  by (transfer; simp add: truncate_down_uminus_eq truncate_up_uminus_eq)+

context
begin

qualified lemma compute_float_round_down[code]:
  "float_round_down prec (Float m e) =
    (let d = bitlen \<bar>m\<bar> - int prec - 1 in
      if 0 < d then Float (div_twopow m (nat d)) (e + d)
      else Float m e)"
  using Float.compute_float_down[of "Suc prec - bitlen \m\ - e" m e, symmetric]
  by transfer
    (simp add: field_simps abs_mult log_mult bitlen_alt_def truncate_down_def
      cong del: if_weak_cong)

qualified lemma compute_float_round_up[code]:
  "float_round_up prec x = - float_round_down prec (-x)"
  by transfer (simp add: truncate_down_uminus_eq)

end

lemma truncate_up_nonneg_mono:
  assumes "0 \ x" "x \ y"
  shows "truncate_up prec x \ truncate_up prec y"
proof -
  consider "\log 2 x\ = \log 2 y\" | "\log 2 x\ \ \log 2 y\" "0 < x" | "x \ 0"
    by arith
  then show ?thesis
  proof cases
    case 1
    then show ?thesis
      using assms
      by (auto simp: truncate_up_def round_up_def intro!: ceiling_mono)
  next
    case 2
    from assms \<open>0 < x\<close> have "log 2 x \<le> log 2 y"
      by auto
    with \<open>\<lfloor>log 2 x\<rfloor> \<noteq> \<lfloor>log 2 y\<rfloor>\<close>
    have logless: "log 2 x < log 2 y"
      by linarith
    have flogless: "\log 2 x\ < \log 2 y\"
      using \<open>\<lfloor>log 2 x\<rfloor> \<noteq> \<lfloor>log 2 y\<rfloor>\<close> \<open>log 2 x \<le> log 2 y\<close> by linarith
    have "truncate_up prec x =
      real_of_int \<lceil>x * 2 powr real_of_int (int prec - \<lfloor>log 2 x\<rfloor> )\<rceil> * 2 powr - real_of_int (int prec - \<lfloor>log 2 x\<rfloor>)"
      using assms by (simp add: truncate_up_def round_up_def)
    also have "\x * 2 powr real_of_int (int prec - \log 2 x\)\ \ (2 ^ (Suc prec))"
    proof (simp only: ceiling_le_iff)
      have "x * 2 powr real_of_int (int prec - \log 2 x\) \
        x * (2 powr real (Suc prec) / (2 powr log 2 x))"
        using real_of_int_floor_add_one_ge[of "log 2 x"] assms
        by (auto simp: algebra_simps simp flip: powr_diff intro!: mult_left_mono)
      then show "x * 2 powr real_of_int (int prec - \log 2 x\) \ real_of_int ((2::int) ^ (Suc prec))"
        using \<open>0 < x\<close> by (simp add: powr_realpow powr_add)
    qed
    then have "real_of_int \x * 2 powr real_of_int (int prec - \log 2 x\)\ \ 2 powr int (Suc prec)"
      by (auto simp: powr_realpow powr_add)
        (metis power_Suc of_int_le_numeral_power_cancel_iff)
    also
    have "2 powr - real_of_int (int prec - \log 2 x\) \ 2 powr - real_of_int (int prec - \log 2 y\ + 1)"
      using logless flogless by (auto intro!: floor_mono)
    also have "2 powr real_of_int (int (Suc prec)) \
        2 powr (log 2 y + real_of_int (int prec - \<lfloor>log 2 y\<rfloor> + 1))"
      using assms \<open>0 < x\<close>
      by (auto simp: algebra_simps)
    finally have "truncate_up prec x \
        2 powr (log 2 y + real_of_int (int prec - \<lfloor>log 2 y\<rfloor> + 1)) * 2 powr - real_of_int (int prec - \<lfloor>log 2 y\<rfloor> + 1)"
      by simp
    also have "\ = 2 powr (log 2 y + real_of_int (int prec - \log 2 y\) - real_of_int (int prec - \log 2 y\))"
      by (subst powr_add[symmetric]) simp
    also have "\ = y"
      using \<open>0 < x\<close> assms
      by (simp add: powr_add)
    also have "\ \ truncate_up prec y"
      by (rule truncate_up)
    finally show ?thesis .
  next
    case 3
    then show ?thesis
      using assms
      by (auto intro!: truncate_up_le)
  qed
qed

lemma truncate_up_switch_sign_mono:
  assumes "x \ 0" "0 \ y"
  shows "truncate_up prec x \ truncate_up prec y"
proof -
  note truncate_up_nonpos[OF \<open>x \<le> 0\<close>]
  also note truncate_up_le[OF \<open>0 \<le> y\<close>]
  finally show ?thesis .
qed

lemma truncate_down_switch_sign_mono:
  assumes "x \ 0"
    and "0 \ y"
    and "x \ y"
  shows "truncate_down prec x \ truncate_down prec y"
proof -
  note truncate_down_le[OF \<open>x \<le> 0\<close>]
  also note truncate_down_nonneg[OF \<open>0 \<le> y\<close>]
  finally show ?thesis .
qed

lemma truncate_down_nonneg_mono:
  assumes "0 \ x" "x \ y"
  shows "truncate_down prec x \ truncate_down prec y"
proof -
  consider "x \ 0" | "\log 2 \x\\ = \log 2 \y\\" |
    "0 < x" "\log 2 \x\\ \ \log 2 \y\\"
    by arith
  then show ?thesis
  proof cases
    case 1
    with assms have "x = 0" "0 \ y" by simp_all
    then show ?thesis
      by (auto intro!: truncate_down_nonneg)
  next
    case 2
    then show ?thesis
      using assms
      by (auto simp: truncate_down_def round_down_def intro!: floor_mono)
  next
    case 3
    from \<open>0 < x\<close> have "log 2 x \<le> log 2 y" "0 < y" "0 \<le> y"
      using assms by auto
    with \<open>\<lfloor>log 2 \<bar>x\<bar>\<rfloor> \<noteq> \<lfloor>log 2 \<bar>y\<bar>\<rfloor>\<close>
    have logless: "log 2 x < log 2 y" and flogless: "\log 2 x\ < \log 2 y\"
      unfolding atomize_conj abs_of_pos[OF \<open>0 < x\<close>] abs_of_pos[OF \<open>0 < y\<close>]
      by (metis floor_less_cancel linorder_cases not_le)
    have "2 powr prec \ y * 2 powr real prec / (2 powr log 2 y)"
      using \<open>0 < y\<close> by simp
    also have "\ \ y * 2 powr real (Suc prec) / (2 powr (real_of_int \log 2 y\ + 1))"
      using \<open>0 \<le> y\<close> \<open>0 \<le> x\<close> assms(2)
      by (auto intro!: powr_mono divide_left_mono
          simp: of_nat_diff powr_add powr_diff)
    also have "\ = y * 2 powr real (Suc prec) / (2 powr real_of_int \log 2 y\ * 2)"
      by (auto simp: powr_add)
    finally have "(2 ^ prec) \ \y * 2 powr real_of_int (int (Suc prec) - \log 2 \y\\ - 1)\"
      using \<open>0 \<le> y\<close>
      by (auto simp: powr_diff le_floor_iff powr_realpow powr_add)
    then have "(2 ^ (prec)) * 2 powr - real_of_int (int prec - \log 2 \y\\) \ truncate_down prec y"
      by (auto simp: truncate_down_def round_down_def)
    moreover have "x \ (2 ^ prec) * 2 powr - real_of_int (int prec - \log 2 \y\\)"
    proof -
      have "x = 2 powr (log 2 \x\)" using \0 < x\ by simp
      also have "\ \ (2 ^ (Suc prec )) * 2 powr - real_of_int (int prec - \log 2 \x\\)"
        using real_of_int_floor_add_one_ge[of "log 2 \x\"] \0 < x\
        by (auto simp flip: powr_realpow powr_add simp: algebra_simps powr_mult_base le_powr_iff)
      also
      have "2 powr - real_of_int (int prec - \log 2 \x\\) \ 2 powr - real_of_int (int prec - \log 2 \y\\ + 1)"
        using logless flogless \<open>x > 0\<close> \<open>y > 0\<close>
        by (auto intro!: floor_mono)
      finally show ?thesis
        by (auto simp flip: powr_realpow simp: powr_diff assms)
    qed
    ultimately show ?thesis
      by (metis dual_order.trans truncate_down)
  qed
qed

lemma truncate_down_eq_truncate_up: "truncate_down p x = - truncate_up p (-x)"
  and truncate_up_eq_truncate_down: "truncate_up p x = - truncate_down p (-x)"
  by (auto simp: truncate_up_uminus_eq truncate_down_uminus_eq)

lemma truncate_down_mono: "x \ y \ truncate_down p x \ truncate_down p y"
  by (smt (verit) truncate_down_nonneg_mono truncate_up_nonneg_mono truncate_up_uminus_eq)

lemma truncate_up_mono: "x \ y \ truncate_up p x \ truncate_up p y"
  by (simp add: truncate_up_eq_truncate_down truncate_down_mono)

lemma truncate_up_nonneg: "0 \ truncate_up p x" if "0 \ x"
  by (simp add: that truncate_up_le)

lemma truncate_up_pos: "0 < truncate_up p x" if "0 < x"
  by (meson less_le_trans that truncate_up)

lemma truncate_up_less_zero_iff[simp]: "truncate_up p x < 0 \ x < 0"
  by (smt (verit) truncate_down_pos truncate_down_uminus_eq truncate_up_nonneg)

lemma truncate_up_nonneg_iff[simp]: "truncate_up p x \ 0 \ x \ 0"
  using truncate_up_less_zero_iff[of p x] truncate_up_nonneg[of x]
  by linarith

lemma truncate_down_less_zero_iff[simp]: "truncate_down p x < 0 \ x < 0"
  by (metis le_less_trans not_less_iff_gr_or_eq truncate_down truncate_down_pos truncate_down_zero)

lemma truncate_down_nonneg_iff[simp]: "truncate_down p x \ 0 \ x \ 0"
  using truncate_down_less_zero_iff[of p x] truncate_down_nonneg[of x p]
  by linarith

lemma truncate_down_eq_zero_iff[simp]: "truncate_down prec x = 0 \ x = 0"
  by (metis not_less_iff_gr_or_eq truncate_down_less_zero_iff truncate_down_pos truncate_down_zero)

lemma truncate_up_eq_zero_iff[simp]: "truncate_up prec x = 0 \ x = 0"
  by (metis not_less_iff_gr_or_eq truncate_up_less_zero_iff truncate_up_pos truncate_up_zero)

subsection \<open>Approximation of positive rationals\<close>

lemma div_mult_twopow_eq: "a div ((2::nat) ^ n) div b = a div (b * 2 ^ n)" for a b :: nat
  by (cases "b = 0") (simp_all add: div_mult2_eq[symmetric] ac_simps)

lemma real_div_nat_eq_floor_of_divide: "a div b = real_of_int \a / b\" for a b :: nat
  by (simp add: floor_divide_of_nat_eq [of a b])

definition "rat_precision prec x y =
  (let d = bitlen x - bitlen y
   in int prec - d + (if Float (abs x) 0 < Float (abs y) d then 1 else 0))"

lemma floor_log_divide_eq:
  assumes "i > 0" "j > 0" "p > 1"
  shows "\log p (i / j)\ = floor (log p i) - floor (log p j) -
    (if i \<ge> j * p powr (floor (log p i) - floor (log p j)) then 0 else 1)"
proof -
  let ?l = "log p"
  let ?fl = "\x. floor (?l x)"
  have "\?l (i / j)\ = \?l i - ?l j\" using assms
    by (auto simp: log_divide)
  also have "\ = floor (real_of_int (?fl i - ?fl j) + (?l i - ?fl i - (?l j - ?fl j)))"
    (is "_ = floor (_ + ?r)")
    by (simp add: algebra_simps)
  also note floor_add2
  also note \<open>p > 1\<close>
  note powr = powr_le_cancel_iff[symmetric, OF \<open>1 < p\<close>, THEN iffD2]
  note powr_strict = powr_less_cancel_iff[symmetric, OF \<open>1 < p\<close>, THEN iffD2]
  have "floor ?r = (if i \ j * p powr (?fl i - ?fl j) then 0 else -1)" (is "_ = ?if")
    using assms
    apply simp
    by (smt (verit, ccfv_SIG) floor_less_iff floor_uminus_of_int le_log_iff mult_powr_eq
        of_int_1 real_of_int_floor_add_one_gt zero_le_floor)
  finally
  show ?thesis by simp
qed

lemma truncate_down_rat_precision:
    "truncate_down prec (real x / real y) = round_down (rat_precision prec x y) (real x / real y)"
  and truncate_up_rat_precision:
    "truncate_up prec (real x / real y) = round_up (rat_precision prec x y) (real x / real y)"
  unfolding truncate_down_def truncate_up_def rat_precision_def
  by (cases x; cases y) (auto simp: floor_log_divide_eq algebra_simps bitlen_alt_def)

lift_definition lapprox_posrat :: "nat \ nat \ nat \ float"
  is "\prec (x::nat) (y::nat). truncate_down prec (x / y)"
  by simp

context
begin

qualified lemma compute_lapprox_posrat[code]:
  "lapprox_posrat prec x y =
   (let
      l = rat_precision prec x y;
      d = if 0 \<le> l then x * 2^nat l div y else x div 2^nat (- l) div y
    in normfloat (Float d (- l)))"
    unfolding div_mult_twopow_eq
    by transfer
      (simp add: round_down_def powr_int real_div_nat_eq_floor_of_divide field_simps Let_def
        truncate_down_rat_precision del: two_powr_minus_int_float)

end

lift_definition rapprox_posrat :: "nat \ nat \ nat \ float"
  is "\prec (x::nat) (y::nat). truncate_up prec (x / y)"
  by simp

context
begin

qualified lemma compute_rapprox_posrat[code]:
  fixes prec x y
  defines "l \ rat_precision prec x y"
  shows "rapprox_posrat prec x y =
   (let
      l = l;
      (r, s) = if 0 \<le> l then (x * 2^nat l, y) else (x, y * 2^nat(-l));
      d = r div s;
      m = r mod s
    in normfloat (Float (d + (if m = 0 \<or> y = 0 then 0 else 1)) (- l)))"
proof (cases "y = 0")
  assume "y = 0"
  then show ?thesis by transfer simp
next
  assume "y \ 0"
  show ?thesis
  proof (cases "0 \ l")
    case True
    define x' where "x' = x * 2 ^ nat l"
    have "int x * 2 ^ nat l = x'"
      by (simp add: x'_def)
    moreover have "real x * 2 powr l = real x'"
      by (simp flip: powr_realpow add: \<open>0 \<le> l\<close> x'_def)
    ultimately show ?thesis
      using ceil_divide_floor_conv[of y x'] powr_realpow[of 2 "nat l"] \0 \ l\ \y \ 0\
        l_def[symmetric, THEN meta_eq_to_obj_eq]
      apply transfer
      apply (auto simp add: round_up_def truncate_up_rat_precision)
      apply (metis floor_divide_of_int_eq of_int_of_nat_eq)
      done
   next
    case False
    define y' where "y' = y * 2 ^ nat (- l)"
    from \<open>y \<noteq> 0\<close> have "y' \<noteq> 0" by (simp add: y'_def)
    have "int y * 2 ^ nat (- l) = y'"
      by (simp add: y'_def)
    moreover have "real x * real_of_int (2::int) powr real_of_int l / real y = x / real y'"
      using \<open>\<not> 0 \<le> l\<close> by (simp flip: powr_realpow add: powr_minus y'_def field_simps)
    ultimately show ?thesis
      using ceil_divide_floor_conv[of y' x] \\ 0 \ l\ \y' \ 0\ \y \ 0\
        l_def[symmetric, THEN meta_eq_to_obj_eq]
      apply transfer
      apply (auto simp add: round_up_def ceil_divide_floor_conv truncate_up_rat_precision)
      apply (metis floor_divide_of_int_eq of_int_of_nat_eq)
      done
  qed
qed

end

lemma rat_precision_pos:
  assumes "0 \ x"
    and "0 < y"
    and "2 * x < y"
  shows "rat_precision n (int x) (int y) > 0"
proof -
  have "0 < x \ log 2 x + 1 = log 2 (2 * x)"
    by (simp add: log_mult)
  then have "bitlen (int x) < bitlen (int y)"
    using assms
    by (simp add: bitlen_alt_def)
      (auto intro!: floor_mono simp add: one_add_floor)
  then show ?thesis
    using assms
    by (auto intro!: pos_add_strict simp add: field_simps rat_precision_def)
qed

lemma rapprox_posrat_less1:
  "0 \ x \ 0 < y \ 2 * x < y \ real_of_float (rapprox_posrat n x y) < 1"
  by transfer (simp add: rat_precision_pos round_up_less1 truncate_up_rat_precision)

lift_definition lapprox_rat :: "nat \ int \ int \ float" is
  "\prec (x::int) (y::int). truncate_down prec (x / y)"
  by simp

context
begin

qualified lemma compute_lapprox_rat[code]:
  "lapprox_rat prec x y =
   (if y = 0 then 0
    else if 0 \<le> x then
     (if 0 < y then lapprox_posrat prec (nat x) (nat y)
      else - (rapprox_posrat prec (nat x) (nat (-y))))
      else
       (if 0 < y
        then - (rapprox_posrat prec (nat (-x)) (nat y))
        else lapprox_posrat prec (nat (-x)) (nat (-y))))"
  by transfer (simp add: truncate_up_uminus_eq)

lift_definition rapprox_rat :: "nat \ int \ int \ float" is
  "\prec (x::int) (y::int). truncate_up prec (x / y)"
  by simp

lemma "rapprox_rat = rapprox_posrat"
  by transfer auto

lemma "lapprox_rat = lapprox_posrat"
  by transfer auto

qualified lemma compute_rapprox_rat[code]:
  "rapprox_rat prec x y = - lapprox_rat prec (-x) y"
  by transfer (simp add: truncate_down_uminus_eq)

qualified lemma compute_truncate_down[code]:
  "truncate_down p (Ratreal r) = (let (a, b) = quotient_of r in lapprox_rat p a b)"
  by transfer (auto split: prod.split simp: of_rat_divide dest!: quotient_of_div)

qualified lemma compute_truncate_up[code]:
  "truncate_up p (Ratreal r) = (let (a, b) = quotient_of r in rapprox_rat p a b)"
  by transfer (auto split: prod.split simp:  of_rat_divide dest!: quotient_of_div)

end

subsection \<open>Division\<close>

definition "real_divl prec a b = truncate_down prec (a / b)"

definition "real_divr prec a b = truncate_up prec (a / b)"

lift_definition float_divl :: "nat \ float \ float \ float" is real_divl
  by (simp add: real_divl_def)

context
begin

qualified lemma compute_float_divl[code]:
  "float_divl prec (Float m1 s1) (Float m2 s2) = lapprox_rat prec m1 m2 * Float 1 (s1 - s2)"
  apply transfer
  unfolding real_divl_def of_int_1 mult_1 truncate_down_shift_int[symmetric]
  apply (simp add: powr_diff powr_minus)
  done

lift_definition float_divr :: "nat \ float \ float \ float" is real_divr
  by (simp add: real_divr_def)

qualified lemma compute_float_divr[code]:
  "float_divr prec x y = - float_divl prec (-x) y"
  by transfer (simp add: real_divr_def real_divl_def truncate_down_uminus_eq)

end

subsection \<open>Approximate Addition\<close>

definition "plus_down prec x y = truncate_down prec (x + y)"

definition "plus_up prec x y = truncate_up prec (x + y)"

lemma float_plus_down_float[intro, simp]: "x \ float \ y \ float \ plus_down p x y \ float"
  by (simp add: plus_down_def)

lemma float_plus_up_float[intro, simp]: "x \ float \ y \ float \ plus_up p x y \ float"
  by (simp add: plus_up_def)

lift_definition float_plus_down :: "nat \ float \ float \ float" is plus_down ..

lift_definition float_plus_up :: "nat \ float \ float \ float" is plus_up ..

lemma plus_down: "plus_down prec x y \ x + y"
  and plus_up: "x + y \ plus_up prec x y"
  by (auto simp: plus_down_def truncate_down plus_up_def truncate_up)

lemma float_plus_down: "real_of_float (float_plus_down prec x y) \ x + y"
  and float_plus_up: "x + y \ real_of_float (float_plus_up prec x y)"
  by (transfer; rule plus_down plus_up)+

lemmas plus_down_le = order_trans[OF plus_down]
  and plus_up_le = order_trans[OF _ plus_up]
  and float_plus_down_le = order_trans[OF float_plus_down]
  and float_plus_up_le = order_trans[OF _ float_plus_up]

lemma compute_plus_up[code]: "plus_up p x y = - plus_down p (-x) (-y)"
  using truncate_down_uminus_eq[of p "x + y"]
  by (auto simp: plus_down_def plus_up_def)

lemma truncate_down_log2_eqI:
  assumes "\log 2 \x\\ = \log 2 \y\\"
  assumes "\x * 2 powr (p - \log 2 \x\\)\ = \y * 2 powr (p - \log 2 \x\\)\"
  shows "truncate_down p x = truncate_down p y"
  using assms by (auto simp: truncate_down_def round_down_def)

lemma sum_neq_zeroI:
  "\a\ \ k \ \b\ < k \ a + b \ 0"
  "\a\ > k \ \b\ \ k \ a + b \ 0"
  for a k :: real
  by auto

lemma abs_real_le_2_powr_bitlen[simp]: "\real_of_int m2\ < 2 powr real_of_int (bitlen \m2\)"
proof (cases "m2 = 0")
  case True
  then show ?thesis by simp
next
  case False
  then have "\m2\ < 2 ^ nat (bitlen \m2\)"
    using bitlen_bounds[of "\m2\"]
    by (auto simp: powr_add bitlen_nonneg)
  then show ?thesis
    by (metis bitlen_nonneg powr_int of_int_abs of_int_less_numeral_power_cancel_iff
        zero_less_numeral)
qed

lemma floor_sum_times_2_powr_sgn_eq:
  fixes ai p q :: int
    and a b :: real
  assumes "a * 2 powr p = ai"
    and b_le_1: "\b * 2 powr (p + 1)\ \ 1"
    and leqp: "q \ p"
  shows "\(a + b) * 2 powr q\ = \(2 * ai + sgn b) * 2 powr (q - p - 1)\"
proof -
  consider "b = 0" | "b > 0" | "b < 0" by arith
  then show ?thesis
  proof cases
    case 1
    then show ?thesis
      by (simp flip: assms(1) powr_add add: algebra_simps powr_mult_base)
  next
    case 2
    then have "b * 2 powr p < \b * 2 powr (p + 1)\"
      by simp
    also note b_le_1
    finally have b_less_1: "b * 2 powr real_of_int p < 1" .

    from b_less_1 \<open>b > 0\<close> have floor_eq: "\<lfloor>b * 2 powr real_of_int p\<rfloor> = 0" "\<lfloor>sgn b / 2\<rfloor> = 0"
      by (simp_all add: floor_eq_iff)

    have "\(a + b) * 2 powr q\ = \(a + b) * 2 powr p * 2 powr (q - p)\"
      by (simp add: algebra_simps flip: powr_realpow powr_add)
    also have "\ = \(ai + b * 2 powr p) * 2 powr (q - p)\"
      by (simp add: assms algebra_simps)
    also have "\ = \(ai + b * 2 powr p) / real_of_int ((2::int) ^ nat (p - q))\"
      using assms
      by (simp add: algebra_simps divide_powr_uminus flip: powr_realpow powr_add)
    also have "\ = \ai / real_of_int ((2::int) ^ nat (p - q))\"
      by (simp del: of_int_power add: floor_divide_real_eq_div floor_eq)
    finally have "\(a + b) * 2 powr real_of_int q\ = \real_of_int ai / real_of_int ((2::int) ^ nat (p - q))\" .
    moreover
    have "\(2 * ai + (sgn b)) * 2 powr (real_of_int (q - p) - 1)\ =
        \<lfloor>real_of_int ai / real_of_int ((2::int) ^ nat (p - q))\<rfloor>"
    proof -
      have "\(2 * ai + sgn b) * 2 powr (real_of_int (q - p) - 1)\ = \(ai + sgn b / 2) * 2 powr (q - p)\"
        by (subst powr_diff) (simp add: field_simps)
      also have "\ = \(ai + sgn b / 2) / real_of_int ((2::int) ^ nat (p - q))\"
        using leqp by (simp flip: powr_realpow add: powr_diff)
      also have "\ = \ai / real_of_int ((2::int) ^ nat (p - q))\"
        by (simp del: of_int_power add: floor_divide_real_eq_div floor_eq)
      finally show ?thesis .
    qed
    ultimately show ?thesis by simp
  next
    case 3
    then have floor_eq: "\b * 2 powr (real_of_int p + 1)\ = -1"
      using b_le_1
      by (auto simp: floor_eq_iff algebra_simps pos_divide_le_eq[symmetric] abs_if divide_powr_uminus
        intro!: mult_neg_pos split: if_split_asm)
    have "\(a + b) * 2 powr q\ = \(2*a + 2*b) * 2 powr p * 2 powr (q - p - 1)\"
      by (simp add: algebra_simps powr_mult_base flip: powr_realpow powr_add)
    also have "\ = \(2 * (a * 2 powr p) + 2 * b * 2 powr p) * 2 powr (q - p - 1)\"
      by (simp add: algebra_simps)
    also have "\ = \(2 * ai + b * 2 powr (p + 1)) / 2 powr (1 - q + p)\"
      using assms by (simp add: algebra_simps powr_mult_base divide_powr_uminus)
    also have "\ = \(2 * ai + b * 2 powr (p + 1)) / real_of_int ((2::int) ^ nat (p - q + 1))\"
      using assms by (simp add: algebra_simps flip: powr_realpow)
    also have "\ = \(2 * ai - 1) / real_of_int ((2::int) ^ nat (p - q + 1))\"
      using \<open>b < 0\<close> assms
      by (simp add: floor_divide_of_int_eq floor_eq floor_divide_real_eq_div
        del: of_int_mult of_int_power of_int_diff)
    also have "\ = \(2 * ai - 1) * 2 powr (q - p - 1)\"
      using assms by (simp add: algebra_simps divide_powr_uminus flip: powr_realpow)
    finally show ?thesis
      using \<open>b < 0\<close> by simp
  qed
qed

lemma log2_abs_int_add_less_half_sgn_eq:
  fixes ai :: int
    and b :: real
  assumes "\b\ \ 1/2"
    and "ai \ 0"
  shows "\log 2 \real_of_int ai + b\\ = \log 2 \ai + sgn b / 2\\"
proof (cases "b = 0")
  case True
  then show ?thesis by simp
next
  case False
  define k where "k = \log 2 \ai\\"
  then have "\log 2 \ai\\ = k"
    by simp
  then have k: "2 powr k \ \ai\" "\ai\ < 2 powr (k + 1)"
    by (simp_all add: floor_log_eq_powr_iff \<open>ai \<noteq> 0\<close>)
  have "k \ 0"
    using assms by (auto simp: k_def)
  define r where "r = \ai\ - 2 ^ nat k"
  have r: "0 \ r" "r < 2 powr k"
    using \<open>k \<ge> 0\<close> k
    by (auto simp: r_def k_def algebra_simps powr_add abs_if powr_int)
  then have "r \ (2::int) ^ nat k - 1"
    using \<open>k \<ge> 0\<close> by (auto simp: powr_int)
  from this[simplified of_int_le_iff[symmetric]] \<open>0 \<le> k\<close>
  have r_le: "r \ 2 powr k - 1"
    by (auto simp: algebra_simps powr_int)
      (metis of_int_1 of_int_add of_int_le_numeral_power_cancel_iff)

  have "\ai\ = 2 powr k + r"
    using \<open>k \<ge> 0\<close> by (auto simp: k_def r_def simp flip: powr_realpow)

  have pos: "\b\ < 1 \ 0 < 2 powr k + (r + b)" for b :: real
    using \<open>0 \<le> k\<close> \<open>ai \<noteq> 0\<close>
    by (auto simp add: r_def powr_realpow[symmetric] abs_if sgn_if algebra_simps
      split: if_split_asm)
  have less: "\sgn ai * b\ < 1"
    and less': "\sgn (sgn ai * b) / 2\ < 1"
    using \<open>\<bar>b\<bar> \<le> _\<close> by (auto simp: abs_if sgn_if split: if_split_asm)

  have floor_eq: "\b::real. \b\ \ 1 / 2 \
      \<lfloor>log 2 (1 + (r + b) / 2 powr k)\<rfloor> = (if r = 0 \<and> b < 0 then -1 else 0)"
    using \<open>k \<ge> 0\<close> r r_le
    by (auto simp: floor_log_eq_powr_iff powr_minus_divide field_simps sgn_if)

  from \<open>real_of_int \<bar>ai\<bar> = _\<close> have "\<bar>ai + b\<bar> = 2 powr k + (r + sgn ai * b)"
    using \<open>\<bar>b\<bar> \<le> _\<close> \<open>0 \<le> k\<close> r
    by (auto simp add: sgn_if abs_if)
  also have "\log 2 \\ = \log 2 (2 powr k + r + sgn (sgn ai * b) / 2)\"
  proof -
    have "2 powr k + (r + (sgn ai) * b) = 2 powr k * (1 + (r + sgn ai * b) / 2 powr k)"
      by (simp add: field_simps)
    also have "\log 2 \\ = k + \log 2 (1 + (r + sgn ai * b) / 2 powr k)\"
      using pos[OF less]
      by (subst log_mult) (simp_all add: log_mult powr_mult field_simps)
    also
    let ?if = "if r = 0 \ sgn ai * b < 0 then -1 else 0"
    have "\log 2 (1 + (r + sgn ai * b) / 2 powr k)\ = ?if"
      using \<open>\<bar>b\<bar> \<le> _\<close>
      by (intro floor_eq) (auto simp: abs_mult sgn_if)
    also
    have "\ = \log 2 (1 + (r + sgn (sgn ai * b) / 2) / 2 powr k)\"
      by (subst floor_eq) (auto simp: sgn_if)
    also have "k + \ = \log 2 (2 powr k * (1 + (r + sgn (sgn ai * b) / 2) / 2 powr k))\"
      unfolding int_add_floor
      using pos[OF less'] \\b\ \ _\
      by (simp add: field_simps add_log_eq_powr del: floor_add2)
    also have "2 powr k * (1 + (r + sgn (sgn ai * b) / 2) / 2 powr k) =
        2 powr k + r + sgn (sgn ai * b) / 2"
      by (simp add: sgn_if field_simps)
    finally show ?thesis .
  qed
  also have "2 powr k + r + sgn (sgn ai * b) / 2 = \ai + sgn b / 2\"
    unfolding \<open>real_of_int \<bar>ai\<bar> = _\<close>[symmetric] using \<open>ai \<noteq> 0\<close>
    by (auto simp: abs_if sgn_if algebra_simps)
  finally show ?thesis .
qed

context
begin

qualified lemma compute_far_float_plus_down:
  fixes m1 e1 m2 e2 :: int
    and p :: nat
  defines "k1 \ Suc p - nat (bitlen \m1\)"
  assumes H: "bitlen \m2\ \ e1 - e2 - k1 - 2" "m1 \ 0" "m2 \ 0" "e1 \ e2"
  shows "float_plus_down p (Float m1 e1) (Float m2 e2) =
    float_round_down p (Float (m1 * 2 ^ (Suc (Suc k1)) + sgn m2) (e1 - int k1 - 2))"
proof -
  let ?a = "real_of_float (Float m1 e1)"
  let ?b = "real_of_float (Float m2 e2)"
  let ?sum = "?a + ?b"
  let ?shift = "real_of_int e2 - real_of_int e1 + real k1 + 1"
  let ?m1 = "m1 * 2 ^ Suc k1"
  let ?m2 = "m2 * 2 powr ?shift"
  let ?m2' = "sgn m2 / 2"
  let ?e = "e1 - int k1 - 1"

  have sum_eq: "?sum = (?m1 + ?m2) * 2 powr ?e"
--> --------------------

--> maximum size reached

--> --------------------

quality99%

¤ Diese beiden folgenden Angebotsgruppen bietet das Unternehmen0.23Angebot Wie Sie bei der Firma Beratungs- und Dienstleistungen beauftragen können ¤

*Eine klare Vorstellung vom Zielzustand

Wurzel

Suchen

Beweissystem der NASA

Beweissystem Isabelle

NIST Cobol Testsuite

Cephes Mathematical Library

Wiener Entwicklungsmethode

Haftungshinweis

Die Informationen auf dieser Webseite wurden nach bestem Wissen sorgfältig zusammengestellt. Es wird jedoch weder Vollständigkeit, noch Richtigkeit, noch Qualität der bereit gestellten Informationen zugesichert.

Bemerkung:

Die farbliche Syntaxdarstellung ist noch experimentell.