src/HOL/Library/Float.thy
author haftmann
Fri Mar 22 19:18:08 2019 +0000 (3 months ago)
changeset 69946 494934c30f38
parent 69593 3dda49e08b9d
child 70347 e5cd5471c18a
permissions -rw-r--r--
improved code equations taken over from AFP
     1 (*  Title:      HOL/Library/Float.thy
     2     Author:     Johannes Hölzl, Fabian Immler
     3     Copyright   2012  TU München
     4 *)
     5 
     6 section \<open>Floating-Point Numbers\<close>
     7 
     8 theory Float
     9 imports Log_Nat Lattice_Algebras
    10 begin
    11 
    12 definition "float = {m * 2 powr e | (m :: int) (e :: int). True}"
    13 
    14 typedef float = float
    15   morphisms real_of_float float_of
    16   unfolding float_def by auto
    17 
    18 setup_lifting type_definition_float
    19 
    20 declare real_of_float [code_unfold]
    21 
    22 lemmas float_of_inject[simp]
    23 
    24 declare [[coercion "real_of_float :: float \<Rightarrow> real"]]
    25 
    26 lemma real_of_float_eq: "f1 = f2 \<longleftrightarrow> real_of_float f1 = real_of_float f2" for f1 f2 :: float
    27   unfolding real_of_float_inject ..
    28 
    29 declare real_of_float_inverse[simp] float_of_inverse [simp]
    30 declare real_of_float [simp]
    31 
    32 
    33 subsection \<open>Real operations preserving the representation as floating point number\<close>
    34 
    35 lemma floatI: "m * 2 powr e = x \<Longrightarrow> x \<in> float" for m e :: int
    36   by (auto simp: float_def)
    37 
    38 lemma zero_float[simp]: "0 \<in> float"
    39   by (auto simp: float_def)
    40 
    41 lemma one_float[simp]: "1 \<in> float"
    42   by (intro floatI[of 1 0]) simp
    43 
    44 lemma numeral_float[simp]: "numeral i \<in> float"
    45   by (intro floatI[of "numeral i" 0]) simp
    46 
    47 lemma neg_numeral_float[simp]: "- numeral i \<in> float"
    48   by (intro floatI[of "- numeral i" 0]) simp
    49 
    50 lemma real_of_int_float[simp]: "real_of_int x \<in> float" for x :: int
    51   by (intro floatI[of x 0]) simp
    52 
    53 lemma real_of_nat_float[simp]: "real x \<in> float" for x :: nat
    54   by (intro floatI[of x 0]) simp
    55 
    56 lemma two_powr_int_float[simp]: "2 powr (real_of_int i) \<in> float" for i :: int
    57   by (intro floatI[of 1 i]) simp
    58 
    59 lemma two_powr_nat_float[simp]: "2 powr (real i) \<in> float" for i :: nat
    60   by (intro floatI[of 1 i]) simp
    61 
    62 lemma two_powr_minus_int_float[simp]: "2 powr - (real_of_int i) \<in> float" for i :: int
    63   by (intro floatI[of 1 "-i"]) simp
    64 
    65 lemma two_powr_minus_nat_float[simp]: "2 powr - (real i) \<in> float" for i :: nat
    66   by (intro floatI[of 1 "-i"]) simp
    67 
    68 lemma two_powr_numeral_float[simp]: "2 powr numeral i \<in> float"
    69   by (intro floatI[of 1 "numeral i"]) simp
    70 
    71 lemma two_powr_neg_numeral_float[simp]: "2 powr - numeral i \<in> float"
    72   by (intro floatI[of 1 "- numeral i"]) simp
    73 
    74 lemma two_pow_float[simp]: "2 ^ n \<in> float"
    75   by (intro floatI[of 1 n]) (simp add: powr_realpow)
    76 
    77 
    78 lemma plus_float[simp]: "r \<in> float \<Longrightarrow> p \<in> float \<Longrightarrow> r + p \<in> float"
    79   unfolding float_def
    80 proof (safe, simp)
    81   have *: "\<exists>(m::int) (e::int). m1 * 2 powr e1 + m2 * 2 powr e2 = m * 2 powr e"
    82     if "e1 \<le> e2" for e1 m1 e2 m2 :: int
    83   proof -
    84     from that have "m1 * 2 powr e1 + m2 * 2 powr e2 = (m1 + m2 * 2 ^ nat (e2 - e1)) * 2 powr e1"
    85       by (simp add: powr_diff field_simps flip: powr_realpow)
    86     then show ?thesis
    87       by blast
    88   qed
    89   fix e1 m1 e2 m2 :: int
    90   consider "e2 \<le> e1" | "e1 \<le> e2" by (rule linorder_le_cases)
    91   then show "\<exists>(m::int) (e::int). m1 * 2 powr e1 + m2 * 2 powr e2 = m * 2 powr e"
    92   proof cases
    93     case 1
    94     from *[OF this, of m2 m1] show ?thesis
    95       by (simp add: ac_simps)
    96   next
    97     case 2
    98     then show ?thesis by (rule *)
    99   qed
   100 qed
   101 
   102 lemma uminus_float[simp]: "x \<in> float \<Longrightarrow> -x \<in> float"
   103   apply (auto simp: float_def)
   104   apply hypsubst_thin
   105   apply (rename_tac m e)
   106   apply (rule_tac x="-m" in exI)
   107   apply (rule_tac x="e" in exI)
   108   apply (simp add: field_simps)
   109   done
   110 
   111 lemma times_float[simp]: "x \<in> float \<Longrightarrow> y \<in> float \<Longrightarrow> x * y \<in> float"
   112   apply (auto simp: float_def)
   113   apply hypsubst_thin
   114   apply (rename_tac mx my ex ey)
   115   apply (rule_tac x="mx * my" in exI)
   116   apply (rule_tac x="ex + ey" in exI)
   117   apply (simp add: powr_add)
   118   done
   119 
   120 lemma minus_float[simp]: "x \<in> float \<Longrightarrow> y \<in> float \<Longrightarrow> x - y \<in> float"
   121   using plus_float [of x "- y"] by simp
   122 
   123 lemma abs_float[simp]: "x \<in> float \<Longrightarrow> \<bar>x\<bar> \<in> float"
   124   by (cases x rule: linorder_cases[of 0]) auto
   125 
   126 lemma sgn_of_float[simp]: "x \<in> float \<Longrightarrow> sgn x \<in> float"
   127   by (cases x rule: linorder_cases[of 0]) (auto intro!: uminus_float)
   128 
   129 lemma div_power_2_float[simp]: "x \<in> float \<Longrightarrow> x / 2^d \<in> float"
   130   apply (auto simp add: float_def)
   131   apply hypsubst_thin
   132   apply (rename_tac m e)
   133   apply (rule_tac x="m" in exI)
   134   apply (rule_tac x="e - d" in exI)
   135   apply (simp flip: powr_realpow powr_add add: field_simps)
   136   done
   137 
   138 lemma div_power_2_int_float[simp]: "x \<in> float \<Longrightarrow> x / (2::int)^d \<in> float"
   139   apply (auto simp add: float_def)
   140   apply hypsubst_thin
   141   apply (rename_tac m e)
   142   apply (rule_tac x="m" in exI)
   143   apply (rule_tac x="e - d" in exI)
   144   apply (simp flip: powr_realpow powr_add add: field_simps)
   145   done
   146 
   147 lemma div_numeral_Bit0_float[simp]:
   148   assumes "x / numeral n \<in> float"
   149   shows "x / (numeral (Num.Bit0 n)) \<in> float"
   150 proof -
   151   have "(x / numeral n) / 2^1 \<in> float"
   152     by (intro assms div_power_2_float)
   153   also have "(x / numeral n) / 2^1 = x / (numeral (Num.Bit0 n))"
   154     by (induct n) auto
   155   finally show ?thesis .
   156 qed
   157 
   158 lemma div_neg_numeral_Bit0_float[simp]:
   159   assumes "x / numeral n \<in> float"
   160   shows "x / (- numeral (Num.Bit0 n)) \<in> float"
   161 proof -
   162   have "- (x / numeral (Num.Bit0 n)) \<in> float"
   163     using assms by simp
   164   also have "- (x / numeral (Num.Bit0 n)) = x / - numeral (Num.Bit0 n)"
   165     by simp
   166   finally show ?thesis .
   167 qed
   168 
   169 lemma power_float[simp]:
   170   assumes "a \<in> float"
   171   shows "a ^ b \<in> float"
   172 proof -
   173   from assms obtain m e :: int where "a = m * 2 powr e"
   174     by (auto simp: float_def)
   175   then show ?thesis
   176     by (auto intro!: floatI[where m="m^b" and e = "e*b"]
   177       simp: power_mult_distrib powr_realpow[symmetric] powr_powr)
   178 qed
   179 
   180 lift_definition Float :: "int \<Rightarrow> int \<Rightarrow> float" is "\<lambda>(m::int) (e::int). m * 2 powr e"
   181   by simp
   182 declare Float.rep_eq[simp]
   183 
   184 code_datatype Float
   185 
   186 lemma compute_real_of_float[code]:
   187   "real_of_float (Float m e) = (if e \<ge> 0 then m * 2 ^ nat e else m / 2 ^ (nat (-e)))"
   188   by (simp add: powr_int)
   189 
   190 
   191 subsection \<open>Arithmetic operations on floating point numbers\<close>
   192 
   193 instantiation float :: "{ring_1,linorder,linordered_ring,linordered_idom,numeral,equal}"
   194 begin
   195 
   196 lift_definition zero_float :: float is 0 by simp
   197 declare zero_float.rep_eq[simp]
   198 
   199 lift_definition one_float :: float is 1 by simp
   200 declare one_float.rep_eq[simp]
   201 
   202 lift_definition plus_float :: "float \<Rightarrow> float \<Rightarrow> float" is "(+)" by simp
   203 declare plus_float.rep_eq[simp]
   204 
   205 lift_definition times_float :: "float \<Rightarrow> float \<Rightarrow> float" is "(*)" by simp
   206 declare times_float.rep_eq[simp]
   207 
   208 lift_definition minus_float :: "float \<Rightarrow> float \<Rightarrow> float" is "(-)" by simp
   209 declare minus_float.rep_eq[simp]
   210 
   211 lift_definition uminus_float :: "float \<Rightarrow> float" is "uminus" by simp
   212 declare uminus_float.rep_eq[simp]
   213 
   214 lift_definition abs_float :: "float \<Rightarrow> float" is abs by simp
   215 declare abs_float.rep_eq[simp]
   216 
   217 lift_definition sgn_float :: "float \<Rightarrow> float" is sgn by simp
   218 declare sgn_float.rep_eq[simp]
   219 
   220 lift_definition equal_float :: "float \<Rightarrow> float \<Rightarrow> bool" is "(=) :: real \<Rightarrow> real \<Rightarrow> bool" .
   221 
   222 lift_definition less_eq_float :: "float \<Rightarrow> float \<Rightarrow> bool" is "(\<le>)" .
   223 declare less_eq_float.rep_eq[simp]
   224 
   225 lift_definition less_float :: "float \<Rightarrow> float \<Rightarrow> bool" is "(<)" .
   226 declare less_float.rep_eq[simp]
   227 
   228 instance
   229   by standard (transfer; fastforce simp add: field_simps intro: mult_left_mono mult_right_mono)+
   230 
   231 end
   232 
   233 lemma real_of_float [simp]: "real_of_float (of_nat n) = of_nat n"
   234   by (induct n) simp_all
   235 
   236 lemma real_of_float_of_int_eq [simp]: "real_of_float (of_int z) = of_int z"
   237   by (cases z rule: int_diff_cases) (simp_all add: of_rat_diff)
   238 
   239 lemma Float_0_eq_0[simp]: "Float 0 e = 0"
   240   by transfer simp
   241 
   242 lemma real_of_float_power[simp]: "real_of_float (f^n) = real_of_float f^n" for f :: float
   243   by (induct n) simp_all
   244 
   245 lemma real_of_float_min: "real_of_float (min x y) = min (real_of_float x) (real_of_float y)"
   246   and real_of_float_max: "real_of_float (max x y) = max (real_of_float x) (real_of_float y)"
   247   for x y :: float
   248   by (simp_all add: min_def max_def)
   249 
   250 instance float :: unbounded_dense_linorder
   251 proof
   252   fix a b :: float
   253   show "\<exists>c. a < c"
   254     apply (intro exI[of _ "a + 1"])
   255     apply transfer
   256     apply simp
   257     done
   258   show "\<exists>c. c < a"
   259     apply (intro exI[of _ "a - 1"])
   260     apply transfer
   261     apply simp
   262     done
   263   show "\<exists>c. a < c \<and> c < b" if "a < b"
   264     apply (rule exI[of _ "(a + b) * Float 1 (- 1)"])
   265     using that
   266     apply transfer
   267     apply (simp add: powr_minus)
   268     done
   269 qed
   270 
   271 instantiation float :: lattice_ab_group_add
   272 begin
   273 
   274 definition inf_float :: "float \<Rightarrow> float \<Rightarrow> float"
   275   where "inf_float a b = min a b"
   276 
   277 definition sup_float :: "float \<Rightarrow> float \<Rightarrow> float"
   278   where "sup_float a b = max a b"
   279 
   280 instance
   281   by standard (transfer; simp add: inf_float_def sup_float_def real_of_float_min real_of_float_max)+
   282 
   283 end
   284 
   285 lemma float_numeral[simp]: "real_of_float (numeral x :: float) = numeral x"
   286   apply (induct x)
   287   apply simp
   288   apply (simp_all only: numeral_Bit0 numeral_Bit1 real_of_float_eq float_of_inverse
   289           plus_float.rep_eq one_float.rep_eq plus_float numeral_float one_float)
   290   done
   291 
   292 lemma transfer_numeral [transfer_rule]:
   293   "rel_fun (=) pcr_float (numeral :: _ \<Rightarrow> real) (numeral :: _ \<Rightarrow> float)"
   294   by (simp add: rel_fun_def float.pcr_cr_eq cr_float_def)
   295 
   296 lemma float_neg_numeral[simp]: "real_of_float (- numeral x :: float) = - numeral x"
   297   by simp
   298 
   299 lemma transfer_neg_numeral [transfer_rule]:
   300   "rel_fun (=) pcr_float (- numeral :: _ \<Rightarrow> real) (- numeral :: _ \<Rightarrow> float)"
   301   by (simp add: rel_fun_def float.pcr_cr_eq cr_float_def)
   302 
   303 lemma float_of_numeral: "numeral k = float_of (numeral k)"
   304   and float_of_neg_numeral: "- numeral k = float_of (- numeral k)"
   305   unfolding real_of_float_eq by simp_all
   306 
   307 
   308 subsection \<open>Quickcheck\<close>
   309 
   310 instantiation float :: exhaustive
   311 begin
   312 
   313 definition exhaustive_float where
   314   "exhaustive_float f d =
   315     Quickcheck_Exhaustive.exhaustive (\<lambda>x. Quickcheck_Exhaustive.exhaustive (\<lambda>y. f (Float x y)) d) d"
   316 
   317 instance ..
   318 
   319 end
   320 
   321 definition (in term_syntax) [code_unfold]:
   322   "valtermify_float x y = Code_Evaluation.valtermify Float {\<cdot>} x {\<cdot>} y"
   323 
   324 instantiation float :: full_exhaustive
   325 begin
   326 
   327 definition
   328   "full_exhaustive_float f d =
   329     Quickcheck_Exhaustive.full_exhaustive
   330       (\<lambda>x. Quickcheck_Exhaustive.full_exhaustive (\<lambda>y. f (valtermify_float x y)) d) d"
   331 
   332 instance ..
   333 
   334 end
   335 
   336 instantiation float :: random
   337 begin
   338 
   339 definition "Quickcheck_Random.random i =
   340   scomp (Quickcheck_Random.random (2 ^ nat_of_natural i))
   341     (\<lambda>man. scomp (Quickcheck_Random.random i) (\<lambda>exp. Pair (valtermify_float man exp)))"
   342 
   343 instance ..
   344 
   345 end
   346 
   347 
   348 subsection \<open>Represent floats as unique mantissa and exponent\<close>
   349 
   350 lemma int_induct_abs[case_names less]:
   351   fixes j :: int
   352   assumes H: "\<And>n. (\<And>i. \<bar>i\<bar> < \<bar>n\<bar> \<Longrightarrow> P i) \<Longrightarrow> P n"
   353   shows "P j"
   354 proof (induct "nat \<bar>j\<bar>" arbitrary: j rule: less_induct)
   355   case less
   356   show ?case by (rule H[OF less]) simp
   357 qed
   358 
   359 lemma int_cancel_factors:
   360   fixes n :: int
   361   assumes "1 < r"
   362   shows "n = 0 \<or> (\<exists>k i. n = k * r ^ i \<and> \<not> r dvd k)"
   363 proof (induct n rule: int_induct_abs)
   364   case (less n)
   365   have "\<exists>k i. n = k * r ^ Suc i \<and> \<not> r dvd k" if "n \<noteq> 0" "n = m * r" for m
   366   proof -
   367     from that have "\<bar>m \<bar> < \<bar>n\<bar>"
   368       using \<open>1 < r\<close> by (simp add: abs_mult)
   369     from less[OF this] that show ?thesis by auto
   370   qed
   371   then show ?case
   372     by (metis dvd_def monoid_mult_class.mult.right_neutral mult.commute power_0)
   373 qed
   374 
   375 lemma mult_powr_eq_mult_powr_iff_asym:
   376   fixes m1 m2 e1 e2 :: int
   377   assumes m1: "\<not> 2 dvd m1"
   378     and "e1 \<le> e2"
   379   shows "m1 * 2 powr e1 = m2 * 2 powr e2 \<longleftrightarrow> m1 = m2 \<and> e1 = e2"
   380   (is "?lhs \<longleftrightarrow> ?rhs")
   381 proof
   382   show ?rhs if eq: ?lhs
   383   proof -
   384     have "m1 \<noteq> 0"
   385       using m1 unfolding dvd_def by auto
   386     from \<open>e1 \<le> e2\<close> eq have "m1 = m2 * 2 powr nat (e2 - e1)"
   387       by (simp add: powr_diff field_simps)
   388     also have "\<dots> = m2 * 2^nat (e2 - e1)"
   389       by (simp add: powr_realpow)
   390     finally have m1_eq: "m1 = m2 * 2^nat (e2 - e1)"
   391       by linarith
   392     with m1 have "m1 = m2"
   393       by (cases "nat (e2 - e1)") (auto simp add: dvd_def)
   394     then show ?thesis
   395       using eq \<open>m1 \<noteq> 0\<close> by (simp add: powr_inj)
   396   qed
   397   show ?lhs if ?rhs
   398     using that by simp
   399 qed
   400 
   401 lemma mult_powr_eq_mult_powr_iff:
   402   "\<not> 2 dvd m1 \<Longrightarrow> \<not> 2 dvd m2 \<Longrightarrow> m1 * 2 powr e1 = m2 * 2 powr e2 \<longleftrightarrow> m1 = m2 \<and> e1 = e2"
   403   for m1 m2 e1 e2 :: int
   404   using mult_powr_eq_mult_powr_iff_asym[of m1 e1 e2 m2]
   405   using mult_powr_eq_mult_powr_iff_asym[of m2 e2 e1 m1]
   406   by (cases e1 e2 rule: linorder_le_cases) auto
   407 
   408 lemma floatE_normed:
   409   assumes x: "x \<in> float"
   410   obtains (zero) "x = 0"
   411    | (powr) m e :: int where "x = m * 2 powr e" "\<not> 2 dvd m" "x \<noteq> 0"
   412 proof -
   413   have "\<exists>(m::int) (e::int). x = m * 2 powr e \<and> \<not> (2::int) dvd m" if "x \<noteq> 0"
   414   proof -
   415     from x obtain m e :: int where x: "x = m * 2 powr e"
   416       by (auto simp: float_def)
   417     with \<open>x \<noteq> 0\<close> int_cancel_factors[of 2 m] obtain k i where "m = k * 2 ^ i" "\<not> 2 dvd k"
   418       by auto
   419     with \<open>\<not> 2 dvd k\<close> x show ?thesis
   420       apply (rule_tac exI[of _ "k"])
   421       apply (rule_tac exI[of _ "e + int i"])
   422       apply (simp add: powr_add powr_realpow)
   423       done
   424   qed
   425   with that show thesis by blast
   426 qed
   427 
   428 lemma float_normed_cases:
   429   fixes f :: float
   430   obtains (zero) "f = 0"
   431    | (powr) m e :: int where "real_of_float f = m * 2 powr e" "\<not> 2 dvd m" "f \<noteq> 0"
   432 proof (atomize_elim, induct f)
   433   case (float_of y)
   434   then show ?case
   435     by (cases rule: floatE_normed) (auto simp: zero_float_def)
   436 qed
   437 
   438 definition mantissa :: "float \<Rightarrow> int"
   439   where "mantissa f =
   440     fst (SOME p::int \<times> int. (f = 0 \<and> fst p = 0 \<and> snd p = 0) \<or>
   441       (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))"
   442 
   443 definition exponent :: "float \<Rightarrow> int"
   444   where "exponent f =
   445     snd (SOME p::int \<times> int. (f = 0 \<and> fst p = 0 \<and> snd p = 0) \<or>
   446       (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))"
   447 
   448 lemma exponent_0[simp]: "exponent 0 = 0" (is ?E)
   449   and mantissa_0[simp]: "mantissa 0 = 0" (is ?M)
   450 proof -
   451   have "\<And>p::int \<times> int. fst p = 0 \<and> snd p = 0 \<longleftrightarrow> p = (0, 0)"
   452     by auto
   453   then show ?E ?M
   454     by (auto simp add: mantissa_def exponent_def zero_float_def)
   455 qed
   456 
   457 lemma mantissa_exponent: "real_of_float f = mantissa f * 2 powr exponent f" (is ?E)
   458   and mantissa_not_dvd: "f \<noteq> 0 \<Longrightarrow> \<not> 2 dvd mantissa f" (is "_ \<Longrightarrow> ?D")
   459 proof cases
   460   assume [simp]: "f \<noteq> 0"
   461   have "f = mantissa f * 2 powr exponent f \<and> \<not> 2 dvd mantissa f"
   462   proof (cases f rule: float_normed_cases)
   463     case zero
   464     then show ?thesis by simp
   465   next
   466     case (powr m e)
   467     then have "\<exists>p::int \<times> int. (f = 0 \<and> fst p = 0 \<and> snd p = 0) \<or>
   468       (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)"
   469       by auto
   470     then show ?thesis
   471       unfolding exponent_def mantissa_def
   472       by (rule someI2_ex) simp
   473   qed
   474   then show ?E ?D by auto
   475 qed simp
   476 
   477 lemma mantissa_noteq_0: "f \<noteq> 0 \<Longrightarrow> mantissa f \<noteq> 0"
   478   using mantissa_not_dvd[of f] by auto
   479 
   480 lemma mantissa_eq_zero_iff: "mantissa x = 0 \<longleftrightarrow> x = 0"
   481   (is "?lhs \<longleftrightarrow> ?rhs")
   482 proof
   483   show ?rhs if ?lhs
   484   proof -
   485     from that have z: "0 = real_of_float x"
   486       using mantissa_exponent by simp
   487     show ?thesis
   488       by (simp add: zero_float_def z)
   489   qed
   490   show ?lhs if ?rhs
   491     using that by simp
   492 qed
   493 
   494 lemma mantissa_pos_iff: "0 < mantissa x \<longleftrightarrow> 0 < x"
   495   by (auto simp: mantissa_exponent sign_simps)
   496 
   497 lemma mantissa_nonneg_iff: "0 \<le> mantissa x \<longleftrightarrow> 0 \<le> x"
   498   by (auto simp: mantissa_exponent sign_simps zero_le_mult_iff)
   499 
   500 lemma mantissa_neg_iff: "0 > mantissa x \<longleftrightarrow> 0 > x"
   501   by (auto simp: mantissa_exponent sign_simps zero_le_mult_iff)
   502 
   503 lemma
   504   fixes m e :: int
   505   defines "f \<equiv> float_of (m * 2 powr e)"
   506   assumes dvd: "\<not> 2 dvd m"
   507   shows mantissa_float: "mantissa f = m" (is "?M")
   508     and exponent_float: "m \<noteq> 0 \<Longrightarrow> exponent f = e" (is "_ \<Longrightarrow> ?E")
   509 proof cases
   510   assume "m = 0"
   511   with dvd show "mantissa f = m" by auto
   512 next
   513   assume "m \<noteq> 0"
   514   then have f_not_0: "f \<noteq> 0" by (simp add: f_def zero_float_def)
   515   from mantissa_exponent[of f] have "m * 2 powr e = mantissa f * 2 powr exponent f"
   516     by (auto simp add: f_def)
   517   then show ?M ?E
   518     using mantissa_not_dvd[OF f_not_0] dvd
   519     by (auto simp: mult_powr_eq_mult_powr_iff)
   520 qed
   521 
   522 
   523 subsection \<open>Compute arithmetic operations\<close>
   524 
   525 lemma Float_mantissa_exponent: "Float (mantissa f) (exponent f) = f"
   526   unfolding real_of_float_eq mantissa_exponent[of f] by simp
   527 
   528 lemma Float_cases [cases type: float]:
   529   fixes f :: float
   530   obtains (Float) m e :: int where "f = Float m e"
   531   using Float_mantissa_exponent[symmetric]
   532   by (atomize_elim) auto
   533 
   534 lemma denormalize_shift:
   535   assumes f_def: "f = Float m e"
   536     and not_0: "f \<noteq> 0"
   537   obtains i where "m = mantissa f * 2 ^ i" "e = exponent f - i"
   538 proof
   539   from mantissa_exponent[of f] f_def
   540   have "m * 2 powr e = mantissa f * 2 powr exponent f"
   541     by simp
   542   then have eq: "m = mantissa f * 2 powr (exponent f - e)"
   543     by (simp add: powr_diff field_simps)
   544   moreover
   545   have "e \<le> exponent f"
   546   proof (rule ccontr)
   547     assume "\<not> e \<le> exponent f"
   548     then have pos: "exponent f < e" by simp
   549     then have "2 powr (exponent f - e) = 2 powr - real_of_int (e - exponent f)"
   550       by simp
   551     also have "\<dots> = 1 / 2^nat (e - exponent f)"
   552       using pos by (simp flip: powr_realpow add: powr_diff)
   553     finally have "m * 2^nat (e - exponent f) = real_of_int (mantissa f)"
   554       using eq by simp
   555     then have "mantissa f = m * 2^nat (e - exponent f)"
   556       by linarith
   557     with \<open>exponent f < e\<close> have "2 dvd mantissa f"
   558       apply (intro dvdI[where k="m * 2^(nat (e-exponent f)) div 2"])
   559       apply (cases "nat (e - exponent f)")
   560       apply auto
   561       done
   562     then show False using mantissa_not_dvd[OF not_0] by simp
   563   qed
   564   ultimately have "real_of_int m = mantissa f * 2^nat (exponent f - e)"
   565     by (simp flip: powr_realpow)
   566   with \<open>e \<le> exponent f\<close>
   567   show "m = mantissa f * 2 ^ nat (exponent f - e)"
   568     by linarith
   569   show "e = exponent f - nat (exponent f - e)"
   570     using \<open>e \<le> exponent f\<close> by auto
   571 qed
   572 
   573 context
   574 begin
   575 
   576 qualified lemma compute_float_zero[code_unfold, code]: "0 = Float 0 0"
   577   by transfer simp
   578 
   579 qualified lemma compute_float_one[code_unfold, code]: "1 = Float 1 0"
   580   by transfer simp
   581 
   582 lift_definition normfloat :: "float \<Rightarrow> float" is "\<lambda>x. x" .
   583 lemma normloat_id[simp]: "normfloat x = x" by transfer rule
   584 
   585 qualified lemma compute_normfloat[code]:
   586   "normfloat (Float m e) =
   587     (if m mod 2 = 0 \<and> m \<noteq> 0 then normfloat (Float (m div 2) (e + 1))
   588      else if m = 0 then 0 else Float m e)"
   589   by transfer (auto simp add: powr_add zmod_eq_0_iff)
   590 
   591 qualified lemma compute_float_numeral[code_abbrev]: "Float (numeral k) 0 = numeral k"
   592   by transfer simp
   593 
   594 qualified lemma compute_float_neg_numeral[code_abbrev]: "Float (- numeral k) 0 = - numeral k"
   595   by transfer simp
   596 
   597 qualified lemma compute_float_uminus[code]: "- Float m1 e1 = Float (- m1) e1"
   598   by transfer simp
   599 
   600 qualified lemma compute_float_times[code]: "Float m1 e1 * Float m2 e2 = Float (m1 * m2) (e1 + e2)"
   601   by transfer (simp add: field_simps powr_add)
   602 
   603 qualified lemma compute_float_plus[code]:
   604   "Float m1 e1 + Float m2 e2 =
   605     (if m1 = 0 then Float m2 e2
   606      else if m2 = 0 then Float m1 e1
   607      else if e1 \<le> e2 then Float (m1 + m2 * 2^nat (e2 - e1)) e1
   608      else Float (m2 + m1 * 2^nat (e1 - e2)) e2)"
   609   by transfer (simp add: field_simps powr_realpow[symmetric] powr_diff)
   610 
   611 qualified lemma compute_float_minus[code]: "f - g = f + (-g)" for f g :: float
   612   by simp
   613 
   614 qualified lemma compute_float_sgn[code]:
   615   "sgn (Float m1 e1) = (if 0 < m1 then 1 else if m1 < 0 then -1 else 0)"
   616   by transfer (simp add: sgn_mult)
   617 
   618 lift_definition is_float_pos :: "float \<Rightarrow> bool" is "(<) 0 :: real \<Rightarrow> bool" .
   619 
   620 qualified lemma compute_is_float_pos[code]: "is_float_pos (Float m e) \<longleftrightarrow> 0 < m"
   621   by transfer (auto simp add: zero_less_mult_iff not_le[symmetric, of _ 0])
   622 
   623 lift_definition is_float_nonneg :: "float \<Rightarrow> bool" is "(\<le>) 0 :: real \<Rightarrow> bool" .
   624 
   625 qualified lemma compute_is_float_nonneg[code]: "is_float_nonneg (Float m e) \<longleftrightarrow> 0 \<le> m"
   626   by transfer (auto simp add: zero_le_mult_iff not_less[symmetric, of _ 0])
   627 
   628 lift_definition is_float_zero :: "float \<Rightarrow> bool"  is "(=) 0 :: real \<Rightarrow> bool" .
   629 
   630 qualified lemma compute_is_float_zero[code]: "is_float_zero (Float m e) \<longleftrightarrow> 0 = m"
   631   by transfer (auto simp add: is_float_zero_def)
   632 
   633 qualified lemma compute_float_abs[code]: "\<bar>Float m e\<bar> = Float \<bar>m\<bar> e"
   634   by transfer (simp add: abs_mult)
   635 
   636 qualified lemma compute_float_eq[code]: "equal_class.equal f g = is_float_zero (f - g)"
   637   by transfer simp
   638 
   639 end
   640 
   641 
   642 subsection \<open>Lemmas for types \<^typ>\<open>real\<close>, \<^typ>\<open>nat\<close>, \<^typ>\<open>int\<close>\<close>
   643 
   644 lemmas real_of_ints =
   645   of_int_add
   646   of_int_minus
   647   of_int_diff
   648   of_int_mult
   649   of_int_power
   650   of_int_numeral of_int_neg_numeral
   651 
   652 lemmas int_of_reals = real_of_ints[symmetric]
   653 
   654 
   655 subsection \<open>Rounding Real Numbers\<close>
   656 
   657 definition round_down :: "int \<Rightarrow> real \<Rightarrow> real"
   658   where "round_down prec x = \<lfloor>x * 2 powr prec\<rfloor> * 2 powr -prec"
   659 
   660 definition round_up :: "int \<Rightarrow> real \<Rightarrow> real"
   661   where "round_up prec x = \<lceil>x * 2 powr prec\<rceil> * 2 powr -prec"
   662 
   663 lemma round_down_float[simp]: "round_down prec x \<in> float"
   664   unfolding round_down_def
   665   by (auto intro!: times_float simp flip: of_int_minus)
   666 
   667 lemma round_up_float[simp]: "round_up prec x \<in> float"
   668   unfolding round_up_def
   669   by (auto intro!: times_float simp flip: of_int_minus)
   670 
   671 lemma round_up: "x \<le> round_up prec x"
   672   by (simp add: powr_minus_divide le_divide_eq round_up_def ceiling_correct)
   673 
   674 lemma round_down: "round_down prec x \<le> x"
   675   by (simp add: powr_minus_divide divide_le_eq round_down_def)
   676 
   677 lemma round_up_0[simp]: "round_up p 0 = 0"
   678   unfolding round_up_def by simp
   679 
   680 lemma round_down_0[simp]: "round_down p 0 = 0"
   681   unfolding round_down_def by simp
   682 
   683 lemma round_up_diff_round_down: "round_up prec x - round_down prec x \<le> 2 powr -prec"
   684 proof -
   685   have "round_up prec x - round_down prec x = (\<lceil>x * 2 powr prec\<rceil> - \<lfloor>x * 2 powr prec\<rfloor>) * 2 powr -prec"
   686     by (simp add: round_up_def round_down_def field_simps)
   687   also have "\<dots> \<le> 1 * 2 powr -prec"
   688     by (rule mult_mono)
   689       (auto simp flip: of_int_diff simp: ceiling_diff_floor_le_1)
   690   finally show ?thesis by simp
   691 qed
   692 
   693 lemma round_down_shift: "round_down p (x * 2 powr k) = 2 powr k * round_down (p + k) x"
   694   unfolding round_down_def
   695   by (simp add: powr_add powr_mult field_simps powr_diff)
   696     (simp flip: powr_add)
   697 
   698 lemma round_up_shift: "round_up p (x * 2 powr k) = 2 powr k * round_up (p + k) x"
   699   unfolding round_up_def
   700   by (simp add: powr_add powr_mult field_simps powr_diff)
   701     (simp flip: powr_add)
   702 
   703 lemma round_up_uminus_eq: "round_up p (-x) = - round_down p x"
   704   and round_down_uminus_eq: "round_down p (-x) = - round_up p x"
   705   by (auto simp: round_up_def round_down_def ceiling_def)
   706 
   707 lemma round_up_mono: "x \<le> y \<Longrightarrow> round_up p x \<le> round_up p y"
   708   by (auto intro!: ceiling_mono simp: round_up_def)
   709 
   710 lemma round_up_le1:
   711   assumes "x \<le> 1" "prec \<ge> 0"
   712   shows "round_up prec x \<le> 1"
   713 proof -
   714   have "real_of_int \<lceil>x * 2 powr prec\<rceil> \<le> real_of_int \<lceil>2 powr real_of_int prec\<rceil>"
   715     using assms by (auto intro!: ceiling_mono)
   716   also have "\<dots> = 2 powr prec" using assms by (auto simp: powr_int intro!: exI[where x="2^nat prec"])
   717   finally show ?thesis
   718     by (simp add: round_up_def) (simp add: powr_minus inverse_eq_divide)
   719 qed
   720 
   721 lemma round_up_less1:
   722   assumes "x < 1 / 2" "p > 0"
   723   shows "round_up p x < 1"
   724 proof -
   725   have "x * 2 powr p < 1 / 2 * 2 powr p"
   726     using assms by simp
   727   also have "\<dots> \<le> 2 powr p - 1" using \<open>p > 0\<close>
   728     by (auto simp: powr_diff powr_int field_simps self_le_power)
   729   finally show ?thesis using \<open>p > 0\<close>
   730     by (simp add: round_up_def field_simps powr_minus powr_int ceiling_less_iff)
   731 qed
   732 
   733 lemma round_down_ge1:
   734   assumes x: "x \<ge> 1"
   735   assumes prec: "p \<ge> - log 2 x"
   736   shows "1 \<le> round_down p x"
   737 proof cases
   738   assume nonneg: "0 \<le> p"
   739   have "2 powr p = real_of_int \<lfloor>2 powr real_of_int p\<rfloor>"
   740     using nonneg by (auto simp: powr_int)
   741   also have "\<dots> \<le> real_of_int \<lfloor>x * 2 powr p\<rfloor>"
   742     using assms by (auto intro!: floor_mono)
   743   finally show ?thesis
   744     by (simp add: round_down_def) (simp add: powr_minus inverse_eq_divide)
   745 next
   746   assume neg: "\<not> 0 \<le> p"
   747   have "x = 2 powr (log 2 x)"
   748     using x by simp
   749   also have "2 powr (log 2 x) \<ge> 2 powr - p"
   750     using prec by auto
   751   finally have x_le: "x \<ge> 2 powr -p" .
   752 
   753   from neg have "2 powr real_of_int p \<le> 2 powr 0"
   754     by (intro powr_mono) auto
   755   also have "\<dots> \<le> \<lfloor>2 powr 0::real\<rfloor>" by simp
   756   also have "\<dots> \<le> \<lfloor>x * 2 powr (real_of_int p)\<rfloor>"
   757     unfolding of_int_le_iff
   758     using x x_le by (intro floor_mono) (simp add: powr_minus_divide field_simps)
   759   finally show ?thesis
   760     using prec x
   761     by (simp add: round_down_def powr_minus_divide pos_le_divide_eq)
   762 qed
   763 
   764 lemma round_up_le0: "x \<le> 0 \<Longrightarrow> round_up p x \<le> 0"
   765   unfolding round_up_def
   766   by (auto simp: field_simps mult_le_0_iff zero_le_mult_iff)
   767 
   768 
   769 subsection \<open>Rounding Floats\<close>
   770 
   771 definition div_twopow :: "int \<Rightarrow> nat \<Rightarrow> int"
   772   where [simp]: "div_twopow x n = x div (2 ^ n)"
   773 
   774 definition mod_twopow :: "int \<Rightarrow> nat \<Rightarrow> int"
   775   where [simp]: "mod_twopow x n = x mod (2 ^ n)"
   776 
   777 lemma compute_div_twopow[code]:
   778   "div_twopow x n = (if x = 0 \<or> x = -1 \<or> n = 0 then x else div_twopow (x div 2) (n - 1))"
   779   by (cases n) (auto simp: zdiv_zmult2_eq div_eq_minus1)
   780 
   781 lemma compute_mod_twopow[code]:
   782   "mod_twopow x n = (if n = 0 then 0 else x mod 2 + 2 * mod_twopow (x div 2) (n - 1))"
   783   by (cases n) (auto simp: zmod_zmult2_eq)
   784 
   785 lift_definition float_up :: "int \<Rightarrow> float \<Rightarrow> float" is round_up by simp
   786 declare float_up.rep_eq[simp]
   787 
   788 lemma round_up_correct: "round_up e f - f \<in> {0..2 powr -e}"
   789   unfolding atLeastAtMost_iff
   790 proof
   791   have "round_up e f - f \<le> round_up e f - round_down e f"
   792     using round_down by simp
   793   also have "\<dots> \<le> 2 powr -e"
   794     using round_up_diff_round_down by simp
   795   finally show "round_up e f - f \<le> 2 powr - (real_of_int e)"
   796     by simp
   797 qed (simp add: algebra_simps round_up)
   798 
   799 lemma float_up_correct: "real_of_float (float_up e f) - real_of_float f \<in> {0..2 powr -e}"
   800   by transfer (rule round_up_correct)
   801 
   802 lift_definition float_down :: "int \<Rightarrow> float \<Rightarrow> float" is round_down by simp
   803 declare float_down.rep_eq[simp]
   804 
   805 lemma round_down_correct: "f - (round_down e f) \<in> {0..2 powr -e}"
   806   unfolding atLeastAtMost_iff
   807 proof
   808   have "f - round_down e f \<le> round_up e f - round_down e f"
   809     using round_up by simp
   810   also have "\<dots> \<le> 2 powr -e"
   811     using round_up_diff_round_down by simp
   812   finally show "f - round_down e f \<le> 2 powr - (real_of_int e)"
   813     by simp
   814 qed (simp add: algebra_simps round_down)
   815 
   816 lemma float_down_correct: "real_of_float f - real_of_float (float_down e f) \<in> {0..2 powr -e}"
   817   by transfer (rule round_down_correct)
   818 
   819 context
   820 begin
   821 
   822 qualified lemma compute_float_down[code]:
   823   "float_down p (Float m e) =
   824     (if p + e < 0 then Float (div_twopow m (nat (-(p + e)))) (-p) else Float m e)"
   825 proof (cases "p + e < 0")
   826   case True
   827   then have "real_of_int ((2::int) ^ nat (-(p + e))) = 2 powr (-(p + e))"
   828     using powr_realpow[of 2 "nat (-(p + e))"] by simp
   829   also have "\<dots> = 1 / 2 powr p / 2 powr e"
   830     unfolding powr_minus_divide of_int_minus by (simp add: powr_add)
   831   finally show ?thesis
   832     using \<open>p + e < 0\<close>
   833     apply transfer
   834     apply (simp add: ac_simps round_down_def flip: floor_divide_of_int_eq)
   835     proof - (*FIXME*)
   836       fix pa :: int and ea :: int and ma :: int
   837       assume a1: "2 ^ nat (- pa - ea) = 1 / (2 powr real_of_int pa * 2 powr real_of_int ea)"
   838       assume "pa + ea < 0"
   839       have "\<lfloor>real_of_int ma / real_of_int (int 2 ^ nat (- (pa + ea)))\<rfloor> =
   840           \<lfloor>real_of_float (Float ma (pa + ea))\<rfloor>"
   841         using a1 by (simp add: powr_add)
   842       then show "\<lfloor>real_of_int ma * (2 powr real_of_int pa * 2 powr real_of_int ea)\<rfloor> =
   843           ma div 2 ^ nat (- pa - ea)"
   844         by (metis Float.rep_eq add_uminus_conv_diff floor_divide_of_int_eq
   845             minus_add_distrib of_int_simps(3) of_nat_numeral powr_add)
   846     qed
   847 next
   848   case False
   849   then have r: "real_of_int e + real_of_int p = real (nat (e + p))"
   850     by simp
   851   have r: "\<lfloor>(m * 2 powr e) * 2 powr real_of_int p\<rfloor> = (m * 2 powr e) * 2 powr real_of_int p"
   852     by (auto intro: exI[where x="m*2^nat (e+p)"]
   853         simp add: ac_simps powr_add[symmetric] r powr_realpow)
   854   with \<open>\<not> p + e < 0\<close> show ?thesis
   855     by transfer (auto simp add: round_down_def field_simps powr_add powr_minus)
   856 qed
   857 
   858 lemma abs_round_down_le: "\<bar>f - (round_down e f)\<bar> \<le> 2 powr -e"
   859   using round_down_correct[of f e] by simp
   860 
   861 lemma abs_round_up_le: "\<bar>f - (round_up e f)\<bar> \<le> 2 powr -e"
   862   using round_up_correct[of e f] by simp
   863 
   864 lemma round_down_nonneg: "0 \<le> s \<Longrightarrow> 0 \<le> round_down p s"
   865   by (auto simp: round_down_def)
   866 
   867 lemma ceil_divide_floor_conv:
   868   assumes "b \<noteq> 0"
   869   shows "\<lceil>real_of_int a / real_of_int b\<rceil> =
   870     (if b dvd a then a div b else \<lfloor>real_of_int a / real_of_int b\<rfloor> + 1)"
   871 proof (cases "b dvd a")
   872   case True
   873   then show ?thesis
   874     by (simp add: ceiling_def floor_divide_of_int_eq dvd_neg_div
   875              flip: of_int_minus divide_minus_left)
   876 next
   877   case False
   878   then have "a mod b \<noteq> 0"
   879     by auto
   880   then have ne: "real_of_int (a mod b) / real_of_int b \<noteq> 0"
   881     using \<open>b \<noteq> 0\<close> by auto
   882   have "\<lceil>real_of_int a / real_of_int b\<rceil> = \<lfloor>real_of_int a / real_of_int b\<rfloor> + 1"
   883     apply (rule ceiling_eq)
   884     apply (auto simp flip: floor_divide_of_int_eq)
   885   proof -
   886     have "real_of_int \<lfloor>real_of_int a / real_of_int b\<rfloor> \<le> real_of_int a / real_of_int b"
   887       by simp
   888     moreover have "real_of_int \<lfloor>real_of_int a / real_of_int b\<rfloor> \<noteq> real_of_int a / real_of_int b"
   889       apply (subst (2) real_of_int_div_aux)
   890       unfolding floor_divide_of_int_eq
   891       using ne \<open>b \<noteq> 0\<close> apply auto
   892       done
   893     ultimately show "real_of_int \<lfloor>real_of_int a / real_of_int b\<rfloor> < real_of_int a / real_of_int b" by arith
   894   qed
   895   then show ?thesis
   896     using \<open>\<not> b dvd a\<close> by simp
   897 qed
   898 
   899 qualified lemma compute_float_up[code]: "float_up p x = - float_down p (-x)"
   900   by transfer (simp add: round_down_uminus_eq)
   901 
   902 end
   903 
   904 
   905 lemma bitlen_Float:
   906   fixes m e
   907   defines [THEN meta_eq_to_obj_eq]: "f \<equiv> Float m e"
   908   shows "bitlen \<bar>mantissa f\<bar> + exponent f = (if m = 0 then 0 else bitlen \<bar>m\<bar> + e)"
   909 proof (cases "m = 0")
   910   case True
   911   then show ?thesis by (simp add: f_def bitlen_alt_def)
   912 next
   913   case False
   914   then have "f \<noteq> 0"
   915     unfolding real_of_float_eq by (simp add: f_def)
   916   then have "mantissa f \<noteq> 0"
   917     by (simp add: mantissa_eq_zero_iff)
   918   moreover
   919   obtain i where "m = mantissa f * 2 ^ i" "e = exponent f - int i"
   920     by (rule f_def[THEN denormalize_shift, OF \<open>f \<noteq> 0\<close>])
   921   ultimately show ?thesis by (simp add: abs_mult)
   922 qed
   923 
   924 lemma float_gt1_scale:
   925   assumes "1 \<le> Float m e"
   926   shows "0 \<le> e + (bitlen m - 1)"
   927 proof -
   928   have "0 < Float m e" using assms by auto
   929   then have "0 < m" using powr_gt_zero[of 2 e]
   930     by (auto simp: zero_less_mult_iff)
   931   then have "m \<noteq> 0" by auto
   932   show ?thesis
   933   proof (cases "0 \<le> e")
   934     case True
   935     then show ?thesis
   936       using \<open>0 < m\<close> by (simp add: bitlen_alt_def)
   937   next
   938     case False
   939     have "(1::int) < 2" by simp
   940     let ?S = "2^(nat (-e))"
   941     have "inverse (2 ^ nat (- e)) = 2 powr e"
   942       using assms False powr_realpow[of 2 "nat (-e)"]
   943       by (auto simp: powr_minus field_simps)
   944     then have "1 \<le> real_of_int m * inverse ?S"
   945       using assms False powr_realpow[of 2 "nat (-e)"]
   946       by (auto simp: powr_minus)
   947     then have "1 * ?S \<le> real_of_int m * inverse ?S * ?S"
   948       by (rule mult_right_mono) auto
   949     then have "?S \<le> real_of_int m"
   950       unfolding mult.assoc by auto
   951     then have "?S \<le> m"
   952       unfolding of_int_le_iff[symmetric] by auto
   953     from this bitlen_bounds[OF \<open>0 < m\<close>, THEN conjunct2]
   954     have "nat (-e) < (nat (bitlen m))"
   955       unfolding power_strict_increasing_iff[OF \<open>1 < 2\<close>, symmetric]
   956       by (rule order_le_less_trans)
   957     then have "-e < bitlen m"
   958       using False by auto
   959     then show ?thesis
   960       by auto
   961   qed
   962 qed
   963 
   964 
   965 subsection \<open>Truncating Real Numbers\<close>
   966 
   967 definition truncate_down::"nat \<Rightarrow> real \<Rightarrow> real"
   968   where "truncate_down prec x = round_down (prec - \<lfloor>log 2 \<bar>x\<bar>\<rfloor>) x"
   969 
   970 lemma truncate_down: "truncate_down prec x \<le> x"
   971   using round_down by (simp add: truncate_down_def)
   972 
   973 lemma truncate_down_le: "x \<le> y \<Longrightarrow> truncate_down prec x \<le> y"
   974   by (rule order_trans[OF truncate_down])
   975 
   976 lemma truncate_down_zero[simp]: "truncate_down prec 0 = 0"
   977   by (simp add: truncate_down_def)
   978 
   979 lemma truncate_down_float[simp]: "truncate_down p x \<in> float"
   980   by (auto simp: truncate_down_def)
   981 
   982 definition truncate_up::"nat \<Rightarrow> real \<Rightarrow> real"
   983   where "truncate_up prec x = round_up (prec - \<lfloor>log 2 \<bar>x\<bar>\<rfloor>) x"
   984 
   985 lemma truncate_up: "x \<le> truncate_up prec x"
   986   using round_up by (simp add: truncate_up_def)
   987 
   988 lemma truncate_up_le: "x \<le> y \<Longrightarrow> x \<le> truncate_up prec y"
   989   by (rule order_trans[OF _ truncate_up])
   990 
   991 lemma truncate_up_zero[simp]: "truncate_up prec 0 = 0"
   992   by (simp add: truncate_up_def)
   993 
   994 lemma truncate_up_uminus_eq: "truncate_up prec (-x) = - truncate_down prec x"
   995   and truncate_down_uminus_eq: "truncate_down prec (-x) = - truncate_up prec x"
   996   by (auto simp: truncate_up_def round_up_def truncate_down_def round_down_def ceiling_def)
   997 
   998 lemma truncate_up_float[simp]: "truncate_up p x \<in> float"
   999   by (auto simp: truncate_up_def)
  1000 
  1001 lemma mult_powr_eq: "0 < b \<Longrightarrow> b \<noteq> 1 \<Longrightarrow> 0 < x \<Longrightarrow> x * b powr y = b powr (y + log b x)"
  1002   by (simp_all add: powr_add)
  1003 
  1004 lemma truncate_down_pos:
  1005   assumes "x > 0"
  1006   shows "truncate_down p x > 0"
  1007 proof -
  1008   have "0 \<le> log 2 x - real_of_int \<lfloor>log 2 x\<rfloor>"
  1009     by (simp add: algebra_simps)
  1010   with assms
  1011   show ?thesis
  1012     apply (auto simp: truncate_down_def round_down_def mult_powr_eq
  1013       intro!: ge_one_powr_ge_zero mult_pos_pos)
  1014     by linarith
  1015 qed
  1016 
  1017 lemma truncate_down_nonneg: "0 \<le> y \<Longrightarrow> 0 \<le> truncate_down prec y"
  1018   by (auto simp: truncate_down_def round_down_def)
  1019 
  1020 lemma truncate_down_ge1: "1 \<le> x \<Longrightarrow> 1 \<le> truncate_down p x"
  1021   apply (auto simp: truncate_down_def algebra_simps intro!: round_down_ge1)
  1022   apply linarith
  1023   done
  1024 
  1025 lemma truncate_up_nonpos: "x \<le> 0 \<Longrightarrow> truncate_up prec x \<le> 0"
  1026   by (auto simp: truncate_up_def round_up_def intro!: mult_nonpos_nonneg)
  1027 
  1028 lemma truncate_up_le1:
  1029   assumes "x \<le> 1"
  1030   shows "truncate_up p x \<le> 1"
  1031 proof -
  1032   consider "x \<le> 0" | "x > 0"
  1033     by arith
  1034   then show ?thesis
  1035   proof cases
  1036     case 1
  1037     with truncate_up_nonpos[OF this, of p] show ?thesis
  1038       by simp
  1039   next
  1040     case 2
  1041     then have le: "\<lfloor>log 2 \<bar>x\<bar>\<rfloor> \<le> 0"
  1042       using assms by (auto simp: log_less_iff)
  1043     from assms have "0 \<le> int p" by simp
  1044     from add_mono[OF this le]
  1045     show ?thesis
  1046       using assms by (simp add: truncate_up_def round_up_le1 add_mono)
  1047   qed
  1048 qed
  1049 
  1050 lemma truncate_down_shift_int:
  1051   "truncate_down p (x * 2 powr real_of_int k) = truncate_down p x * 2 powr k"
  1052   by (cases "x = 0")
  1053     (simp_all add: algebra_simps abs_mult log_mult truncate_down_def
  1054       round_down_shift[of _ _ k, simplified])
  1055 
  1056 lemma truncate_down_shift_nat: "truncate_down p (x * 2 powr real k) = truncate_down p x * 2 powr k"
  1057   by (metis of_int_of_nat_eq truncate_down_shift_int)
  1058 
  1059 lemma truncate_up_shift_int: "truncate_up p (x * 2 powr real_of_int k) = truncate_up p x * 2 powr k"
  1060   by (cases "x = 0")
  1061     (simp_all add: algebra_simps abs_mult log_mult truncate_up_def
  1062       round_up_shift[of _ _ k, simplified])
  1063 
  1064 lemma truncate_up_shift_nat: "truncate_up p (x * 2 powr real k) = truncate_up p x * 2 powr k"
  1065   by (metis of_int_of_nat_eq truncate_up_shift_int)
  1066 
  1067 
  1068 subsection \<open>Truncating Floats\<close>
  1069 
  1070 lift_definition float_round_up :: "nat \<Rightarrow> float \<Rightarrow> float" is truncate_up
  1071   by (simp add: truncate_up_def)
  1072 
  1073 lemma float_round_up: "real_of_float x \<le> real_of_float (float_round_up prec x)"
  1074   using truncate_up by transfer simp
  1075 
  1076 lemma float_round_up_zero[simp]: "float_round_up prec 0 = 0"
  1077   by transfer simp
  1078 
  1079 lift_definition float_round_down :: "nat \<Rightarrow> float \<Rightarrow> float" is truncate_down
  1080   by (simp add: truncate_down_def)
  1081 
  1082 lemma float_round_down: "real_of_float (float_round_down prec x) \<le> real_of_float x"
  1083   using truncate_down by transfer simp
  1084 
  1085 lemma float_round_down_zero[simp]: "float_round_down prec 0 = 0"
  1086   by transfer simp
  1087 
  1088 lemmas float_round_up_le = order_trans[OF _ float_round_up]
  1089   and float_round_down_le = order_trans[OF float_round_down]
  1090 
  1091 lemma minus_float_round_up_eq: "- float_round_up prec x = float_round_down prec (- x)"
  1092   and minus_float_round_down_eq: "- float_round_down prec x = float_round_up prec (- x)"
  1093   by (transfer; simp add: truncate_down_uminus_eq truncate_up_uminus_eq)+
  1094 
  1095 context
  1096 begin
  1097 
  1098 qualified lemma compute_float_round_down[code]:
  1099   "float_round_down prec (Float m e) =
  1100     (let d = bitlen \<bar>m\<bar> - int prec - 1 in
  1101       if 0 < d then Float (div_twopow m (nat d)) (e + d)
  1102       else Float m e)"
  1103   using Float.compute_float_down[of "Suc prec - bitlen \<bar>m\<bar> - e" m e, symmetric]
  1104   by transfer
  1105     (simp add: field_simps abs_mult log_mult bitlen_alt_def truncate_down_def
  1106       cong del: if_weak_cong)
  1107 
  1108 qualified lemma compute_float_round_up[code]:
  1109   "float_round_up prec x = - float_round_down prec (-x)"
  1110   by transfer (simp add: truncate_down_uminus_eq)
  1111 
  1112 end
  1113 
  1114 
  1115 subsection \<open>Approximation of positive rationals\<close>
  1116 
  1117 lemma div_mult_twopow_eq: "a div ((2::nat) ^ n) div b = a div (b * 2 ^ n)" for a b :: nat
  1118   by (cases "b = 0") (simp_all add: div_mult2_eq[symmetric] ac_simps)
  1119 
  1120 lemma real_div_nat_eq_floor_of_divide: "a div b = real_of_int \<lfloor>a / b\<rfloor>" for a b :: nat
  1121   by (simp add: floor_divide_of_nat_eq [of a b])
  1122 
  1123 definition "rat_precision prec x y =
  1124   (let d = bitlen x - bitlen y
  1125    in int prec - d + (if Float (abs x) 0 < Float (abs y) d then 1 else 0))"
  1126 
  1127 lemma floor_log_divide_eq:
  1128   assumes "i > 0" "j > 0" "p > 1"
  1129   shows "\<lfloor>log p (i / j)\<rfloor> = floor (log p i) - floor (log p j) -
  1130     (if i \<ge> j * p powr (floor (log p i) - floor (log p j)) then 0 else 1)"
  1131 proof -
  1132   let ?l = "log p"
  1133   let ?fl = "\<lambda>x. floor (?l x)"
  1134   have "\<lfloor>?l (i / j)\<rfloor> = \<lfloor>?l i - ?l j\<rfloor>" using assms
  1135     by (auto simp: log_divide)
  1136   also have "\<dots> = floor (real_of_int (?fl i - ?fl j) + (?l i - ?fl i - (?l j - ?fl j)))"
  1137     (is "_ = floor (_ + ?r)")
  1138     by (simp add: algebra_simps)
  1139   also note floor_add2
  1140   also note \<open>p > 1\<close>
  1141   note powr = powr_le_cancel_iff[symmetric, OF \<open>1 < p\<close>, THEN iffD2]
  1142   note powr_strict = powr_less_cancel_iff[symmetric, OF \<open>1 < p\<close>, THEN iffD2]
  1143   have "floor ?r = (if i \<ge> j * p powr (?fl i - ?fl j) then 0 else -1)" (is "_ = ?if")
  1144     using assms
  1145     by (linarith |
  1146       auto
  1147         intro!: floor_eq2
  1148         intro: powr_strict powr
  1149         simp: powr_diff powr_add divide_simps algebra_simps)+
  1150   finally
  1151   show ?thesis by simp
  1152 qed
  1153 
  1154 lemma truncate_down_rat_precision:
  1155     "truncate_down prec (real x / real y) = round_down (rat_precision prec x y) (real x / real y)"
  1156   and truncate_up_rat_precision:
  1157     "truncate_up prec (real x / real y) = round_up (rat_precision prec x y) (real x / real y)"
  1158   unfolding truncate_down_def truncate_up_def rat_precision_def
  1159   by (cases x; cases y) (auto simp: floor_log_divide_eq algebra_simps bitlen_alt_def)
  1160 
  1161 lift_definition lapprox_posrat :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float"
  1162   is "\<lambda>prec (x::nat) (y::nat). truncate_down prec (x / y)"
  1163   by simp
  1164 
  1165 context
  1166 begin
  1167 
  1168 qualified lemma compute_lapprox_posrat[code]:
  1169   "lapprox_posrat prec x y =
  1170    (let
  1171       l = rat_precision prec x y;
  1172       d = if 0 \<le> l then x * 2^nat l div y else x div 2^nat (- l) div y
  1173     in normfloat (Float d (- l)))"
  1174     unfolding div_mult_twopow_eq
  1175     by transfer
  1176       (simp add: round_down_def powr_int real_div_nat_eq_floor_of_divide field_simps Let_def
  1177         truncate_down_rat_precision del: two_powr_minus_int_float)
  1178 
  1179 end
  1180 
  1181 lift_definition rapprox_posrat :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float"
  1182   is "\<lambda>prec (x::nat) (y::nat). truncate_up prec (x / y)"
  1183   by simp
  1184 
  1185 context
  1186 begin
  1187 
  1188 qualified lemma compute_rapprox_posrat[code]:
  1189   fixes prec x y
  1190   defines "l \<equiv> rat_precision prec x y"
  1191   shows "rapprox_posrat prec x y =
  1192    (let
  1193       l = l;
  1194       (r, s) = if 0 \<le> l then (x * 2^nat l, y) else (x, y * 2^nat(-l));
  1195       d = r div s;
  1196       m = r mod s
  1197     in normfloat (Float (d + (if m = 0 \<or> y = 0 then 0 else 1)) (- l)))"
  1198 proof (cases "y = 0")
  1199   assume "y = 0"
  1200   then show ?thesis by transfer simp
  1201 next
  1202   assume "y \<noteq> 0"
  1203   show ?thesis
  1204   proof (cases "0 \<le> l")
  1205     case True
  1206     define x' where "x' = x * 2 ^ nat l"
  1207     have "int x * 2 ^ nat l = x'"
  1208       by (simp add: x'_def)
  1209     moreover have "real x * 2 powr l = real x'"
  1210       by (simp flip: powr_realpow add: \<open>0 \<le> l\<close> x'_def)
  1211     ultimately show ?thesis
  1212       using ceil_divide_floor_conv[of y x'] powr_realpow[of 2 "nat l"] \<open>0 \<le> l\<close> \<open>y \<noteq> 0\<close>
  1213         l_def[symmetric, THEN meta_eq_to_obj_eq]
  1214       apply transfer
  1215       apply (auto simp add: round_up_def truncate_up_rat_precision)
  1216       apply (metis dvd_triv_left of_nat_dvd_iff)
  1217       apply (metis floor_divide_of_int_eq of_int_of_nat_eq)
  1218       done
  1219    next
  1220     case False
  1221     define y' where "y' = y * 2 ^ nat (- l)"
  1222     from \<open>y \<noteq> 0\<close> have "y' \<noteq> 0" by (simp add: y'_def)
  1223     have "int y * 2 ^ nat (- l) = y'"
  1224       by (simp add: y'_def)
  1225     moreover have "real x * real_of_int (2::int) powr real_of_int l / real y = x / real y'"
  1226       using \<open>\<not> 0 \<le> l\<close> by (simp flip: powr_realpow add: powr_minus y'_def field_simps)
  1227     ultimately show ?thesis
  1228       using ceil_divide_floor_conv[of y' x] \<open>\<not> 0 \<le> l\<close> \<open>y' \<noteq> 0\<close> \<open>y \<noteq> 0\<close>
  1229         l_def[symmetric, THEN meta_eq_to_obj_eq]
  1230       apply transfer
  1231       apply (auto simp add: round_up_def ceil_divide_floor_conv truncate_up_rat_precision)
  1232       apply (metis dvd_triv_left of_nat_dvd_iff)
  1233       apply (metis floor_divide_of_int_eq of_int_of_nat_eq)
  1234       done
  1235   qed
  1236 qed
  1237 
  1238 end
  1239 
  1240 lemma rat_precision_pos:
  1241   assumes "0 \<le> x"
  1242     and "0 < y"
  1243     and "2 * x < y"
  1244   shows "rat_precision n (int x) (int y) > 0"
  1245 proof -
  1246   have "0 < x \<Longrightarrow> log 2 x + 1 = log 2 (2 * x)"
  1247     by (simp add: log_mult)
  1248   then have "bitlen (int x) < bitlen (int y)"
  1249     using assms
  1250     by (simp add: bitlen_alt_def)
  1251       (auto intro!: floor_mono simp add: one_add_floor)
  1252   then show ?thesis
  1253     using assms
  1254     by (auto intro!: pos_add_strict simp add: field_simps rat_precision_def)
  1255 qed
  1256 
  1257 lemma rapprox_posrat_less1:
  1258   "0 \<le> x \<Longrightarrow> 0 < y \<Longrightarrow> 2 * x < y \<Longrightarrow> real_of_float (rapprox_posrat n x y) < 1"
  1259   by transfer (simp add: rat_precision_pos round_up_less1 truncate_up_rat_precision)
  1260 
  1261 lift_definition lapprox_rat :: "nat \<Rightarrow> int \<Rightarrow> int \<Rightarrow> float" is
  1262   "\<lambda>prec (x::int) (y::int). truncate_down prec (x / y)"
  1263   by simp
  1264 
  1265 context
  1266 begin
  1267 
  1268 qualified lemma compute_lapprox_rat[code]:
  1269   "lapprox_rat prec x y =
  1270    (if y = 0 then 0
  1271     else if 0 \<le> x then
  1272      (if 0 < y then lapprox_posrat prec (nat x) (nat y)
  1273       else - (rapprox_posrat prec (nat x) (nat (-y))))
  1274       else
  1275        (if 0 < y
  1276         then - (rapprox_posrat prec (nat (-x)) (nat y))
  1277         else lapprox_posrat prec (nat (-x)) (nat (-y))))"
  1278   by transfer (simp add: truncate_up_uminus_eq)
  1279 
  1280 lift_definition rapprox_rat :: "nat \<Rightarrow> int \<Rightarrow> int \<Rightarrow> float" is
  1281   "\<lambda>prec (x::int) (y::int). truncate_up prec (x / y)"
  1282   by simp
  1283 
  1284 lemma "rapprox_rat = rapprox_posrat"
  1285   by transfer auto
  1286 
  1287 lemma "lapprox_rat = lapprox_posrat"
  1288   by transfer auto
  1289 
  1290 qualified lemma compute_rapprox_rat[code]:
  1291   "rapprox_rat prec x y = - lapprox_rat prec (-x) y"
  1292   by transfer (simp add: truncate_down_uminus_eq)
  1293 
  1294 qualified lemma compute_truncate_down[code]:
  1295   "truncate_down p (Ratreal r) = (let (a, b) = quotient_of r in lapprox_rat p a b)"
  1296   by transfer (auto split: prod.split simp: of_rat_divide dest!: quotient_of_div)
  1297 
  1298 qualified lemma compute_truncate_up[code]:
  1299   "truncate_up p (Ratreal r) = (let (a, b) = quotient_of r in rapprox_rat p a b)"
  1300   by transfer (auto split: prod.split simp:  of_rat_divide dest!: quotient_of_div)
  1301 
  1302 end
  1303 
  1304 
  1305 subsection \<open>Division\<close>
  1306 
  1307 definition "real_divl prec a b = truncate_down prec (a / b)"
  1308 
  1309 definition "real_divr prec a b = truncate_up prec (a / b)"
  1310 
  1311 lift_definition float_divl :: "nat \<Rightarrow> float \<Rightarrow> float \<Rightarrow> float" is real_divl
  1312   by (simp add: real_divl_def)
  1313 
  1314 context
  1315 begin
  1316 
  1317 qualified lemma compute_float_divl[code]:
  1318   "float_divl prec (Float m1 s1) (Float m2 s2) = lapprox_rat prec m1 m2 * Float 1 (s1 - s2)"
  1319   apply transfer
  1320   unfolding real_divl_def of_int_1 mult_1 truncate_down_shift_int[symmetric]
  1321   apply (simp add: powr_diff powr_minus)
  1322   done
  1323 
  1324 lift_definition float_divr :: "nat \<Rightarrow> float \<Rightarrow> float \<Rightarrow> float" is real_divr
  1325   by (simp add: real_divr_def)
  1326 
  1327 qualified lemma compute_float_divr[code]:
  1328   "float_divr prec x y = - float_divl prec (-x) y"
  1329   by transfer (simp add: real_divr_def real_divl_def truncate_down_uminus_eq)
  1330 
  1331 end
  1332 
  1333 
  1334 subsection \<open>Approximate Power\<close>
  1335 
  1336 lemma div2_less_self[termination_simp]: "odd n \<Longrightarrow> n div 2 < n" for n :: nat
  1337   by (simp add: odd_pos)
  1338 
  1339 fun power_down :: "nat \<Rightarrow> real \<Rightarrow> nat \<Rightarrow> real"
  1340 where
  1341   "power_down p x 0 = 1"
  1342 | "power_down p x (Suc n) =
  1343     (if odd n then truncate_down (Suc p) ((power_down p x (Suc n div 2))\<^sup>2)
  1344      else truncate_down (Suc p) (x * power_down p x n))"
  1345 
  1346 fun power_up :: "nat \<Rightarrow> real \<Rightarrow> nat \<Rightarrow> real"
  1347 where
  1348   "power_up p x 0 = 1"
  1349 | "power_up p x (Suc n) =
  1350     (if odd n then truncate_up p ((power_up p x (Suc n div 2))\<^sup>2)
  1351      else truncate_up p (x * power_up p x n))"
  1352 
  1353 lift_definition power_up_fl :: "nat \<Rightarrow> float \<Rightarrow> nat \<Rightarrow> float" is power_up
  1354   by (induct_tac rule: power_up.induct) simp_all
  1355 
  1356 lift_definition power_down_fl :: "nat \<Rightarrow> float \<Rightarrow> nat \<Rightarrow> float" is power_down
  1357   by (induct_tac rule: power_down.induct) simp_all
  1358 
  1359 lemma power_float_transfer[transfer_rule]:
  1360   "(rel_fun pcr_float (rel_fun (=) pcr_float)) (^) (^)"
  1361   unfolding power_def
  1362   by transfer_prover
  1363 
  1364 lemma compute_power_up_fl[code]:
  1365     "power_up_fl p x 0 = 1"
  1366     "power_up_fl p x (Suc n) =
  1367       (if odd n then float_round_up p ((power_up_fl p x (Suc n div 2))\<^sup>2)
  1368        else float_round_up p (x * power_up_fl p x n))"
  1369   and compute_power_down_fl[code]:
  1370     "power_down_fl p x 0 = 1"
  1371     "power_down_fl p x (Suc n) =
  1372       (if odd n then float_round_down (Suc p) ((power_down_fl p x (Suc n div 2))\<^sup>2)
  1373        else float_round_down (Suc p) (x * power_down_fl p x n))"
  1374   unfolding atomize_conj by transfer simp
  1375 
  1376 lemma power_down_pos: "0 < x \<Longrightarrow> 0 < power_down p x n"
  1377   by (induct p x n rule: power_down.induct)
  1378     (auto simp del: odd_Suc_div_two intro!: truncate_down_pos)
  1379 
  1380 lemma power_down_nonneg: "0 \<le> x \<Longrightarrow> 0 \<le> power_down p x n"
  1381   by (induct p x n rule: power_down.induct)
  1382     (auto simp del: odd_Suc_div_two intro!: truncate_down_nonneg mult_nonneg_nonneg)
  1383 
  1384 lemma power_down: "0 \<le> x \<Longrightarrow> power_down p x n \<le> x ^ n"
  1385 proof (induct p x n rule: power_down.induct)
  1386   case (2 p x n)
  1387   have ?case if "odd n"
  1388   proof -
  1389     from that 2 have "(power_down p x (Suc n div 2)) ^ 2 \<le> (x ^ (Suc n div 2)) ^ 2"
  1390       by (auto intro: power_mono power_down_nonneg simp del: odd_Suc_div_two)
  1391     also have "\<dots> = x ^ (Suc n div 2 * 2)"
  1392       by (simp flip: power_mult)
  1393     also have "Suc n div 2 * 2 = Suc n"
  1394       using \<open>odd n\<close> by presburger
  1395     finally show ?thesis
  1396       using that by (auto intro!: truncate_down_le simp del: odd_Suc_div_two)
  1397   qed
  1398   then show ?case
  1399     by (auto intro!: truncate_down_le mult_left_mono 2 mult_nonneg_nonneg power_down_nonneg)
  1400 qed simp
  1401 
  1402 lemma power_up: "0 \<le> x \<Longrightarrow> x ^ n \<le> power_up p x n"
  1403 proof (induct p x n rule: power_up.induct)
  1404   case (2 p x n)
  1405   have ?case if "odd n"
  1406   proof -
  1407     from that even_Suc have "Suc n = Suc n div 2 * 2"
  1408       by presburger
  1409     then have "x ^ Suc n \<le> (x ^ (Suc n div 2))\<^sup>2"
  1410       by (simp flip: power_mult)
  1411     also from that 2 have "\<dots> \<le> (power_up p x (Suc n div 2))\<^sup>2"
  1412       by (auto intro: power_mono simp del: odd_Suc_div_two)
  1413     finally show ?thesis
  1414       using that by (auto intro!: truncate_up_le simp del: odd_Suc_div_two)
  1415   qed
  1416   then show ?case
  1417     by (auto intro!: truncate_up_le mult_left_mono 2)
  1418 qed simp
  1419 
  1420 lemmas power_up_le = order_trans[OF _ power_up]
  1421   and power_up_less = less_le_trans[OF _ power_up]
  1422   and power_down_le = order_trans[OF power_down]
  1423 
  1424 lemma power_down_fl: "0 \<le> x \<Longrightarrow> power_down_fl p x n \<le> x ^ n"
  1425   by transfer (rule power_down)
  1426 
  1427 lemma power_up_fl: "0 \<le> x \<Longrightarrow> x ^ n \<le> power_up_fl p x n"
  1428   by transfer (rule power_up)
  1429 
  1430 lemma real_power_up_fl: "real_of_float (power_up_fl p x n) = power_up p x n"
  1431   by transfer simp
  1432 
  1433 lemma real_power_down_fl: "real_of_float (power_down_fl p x n) = power_down p x n"
  1434   by transfer simp
  1435 
  1436 
  1437 subsection \<open>Approximate Addition\<close>
  1438 
  1439 definition "plus_down prec x y = truncate_down prec (x + y)"
  1440 
  1441 definition "plus_up prec x y = truncate_up prec (x + y)"
  1442 
  1443 lemma float_plus_down_float[intro, simp]: "x \<in> float \<Longrightarrow> y \<in> float \<Longrightarrow> plus_down p x y \<in> float"
  1444   by (simp add: plus_down_def)
  1445 
  1446 lemma float_plus_up_float[intro, simp]: "x \<in> float \<Longrightarrow> y \<in> float \<Longrightarrow> plus_up p x y \<in> float"
  1447   by (simp add: plus_up_def)
  1448 
  1449 lift_definition float_plus_down :: "nat \<Rightarrow> float \<Rightarrow> float \<Rightarrow> float" is plus_down ..
  1450 
  1451 lift_definition float_plus_up :: "nat \<Rightarrow> float \<Rightarrow> float \<Rightarrow> float" is plus_up ..
  1452 
  1453 lemma plus_down: "plus_down prec x y \<le> x + y"
  1454   and plus_up: "x + y \<le> plus_up prec x y"
  1455   by (auto simp: plus_down_def truncate_down plus_up_def truncate_up)
  1456 
  1457 lemma float_plus_down: "real_of_float (float_plus_down prec x y) \<le> x + y"
  1458   and float_plus_up: "x + y \<le> real_of_float (float_plus_up prec x y)"
  1459   by (transfer; rule plus_down plus_up)+
  1460 
  1461 lemmas plus_down_le = order_trans[OF plus_down]
  1462   and plus_up_le = order_trans[OF _ plus_up]
  1463   and float_plus_down_le = order_trans[OF float_plus_down]
  1464   and float_plus_up_le = order_trans[OF _ float_plus_up]
  1465 
  1466 lemma compute_plus_up[code]: "plus_up p x y = - plus_down p (-x) (-y)"
  1467   using truncate_down_uminus_eq[of p "x + y"]
  1468   by (auto simp: plus_down_def plus_up_def)
  1469 
  1470 lemma truncate_down_log2_eqI:
  1471   assumes "\<lfloor>log 2 \<bar>x\<bar>\<rfloor> = \<lfloor>log 2 \<bar>y\<bar>\<rfloor>"
  1472   assumes "\<lfloor>x * 2 powr (p - \<lfloor>log 2 \<bar>x\<bar>\<rfloor>)\<rfloor> = \<lfloor>y * 2 powr (p - \<lfloor>log 2 \<bar>x\<bar>\<rfloor>)\<rfloor>"
  1473   shows "truncate_down p x = truncate_down p y"
  1474   using assms by (auto simp: truncate_down_def round_down_def)
  1475 
  1476 lemma sum_neq_zeroI:
  1477   "\<bar>a\<bar> \<ge> k \<Longrightarrow> \<bar>b\<bar> < k \<Longrightarrow> a + b \<noteq> 0"
  1478   "\<bar>a\<bar> > k \<Longrightarrow> \<bar>b\<bar> \<le> k \<Longrightarrow> a + b \<noteq> 0"
  1479   for a k :: real
  1480   by auto
  1481 
  1482 lemma abs_real_le_2_powr_bitlen[simp]: "\<bar>real_of_int m2\<bar> < 2 powr real_of_int (bitlen \<bar>m2\<bar>)"
  1483 proof (cases "m2 = 0")
  1484   case True
  1485   then show ?thesis by simp
  1486 next
  1487   case False
  1488   then have "\<bar>m2\<bar> < 2 ^ nat (bitlen \<bar>m2\<bar>)"
  1489     using bitlen_bounds[of "\<bar>m2\<bar>"]
  1490     by (auto simp: powr_add bitlen_nonneg)
  1491   then show ?thesis
  1492     by (metis bitlen_nonneg powr_int of_int_abs of_int_less_numeral_power_cancel_iff
  1493         zero_less_numeral)
  1494 qed
  1495 
  1496 lemma floor_sum_times_2_powr_sgn_eq:
  1497   fixes ai p q :: int
  1498     and a b :: real
  1499   assumes "a * 2 powr p = ai"
  1500     and b_le_1: "\<bar>b * 2 powr (p + 1)\<bar> \<le> 1"
  1501     and leqp: "q \<le> p"
  1502   shows "\<lfloor>(a + b) * 2 powr q\<rfloor> = \<lfloor>(2 * ai + sgn b) * 2 powr (q - p - 1)\<rfloor>"
  1503 proof -
  1504   consider "b = 0" | "b > 0" | "b < 0" by arith
  1505   then show ?thesis
  1506   proof cases
  1507     case 1
  1508     then show ?thesis
  1509       by (simp flip: assms(1) powr_add add: algebra_simps powr_mult_base)
  1510   next
  1511     case 2
  1512     then have "b * 2 powr p < \<bar>b * 2 powr (p + 1)\<bar>"
  1513       by simp
  1514     also note b_le_1
  1515     finally have b_less_1: "b * 2 powr real_of_int p < 1" .
  1516 
  1517     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"
  1518       by (simp_all add: floor_eq_iff)
  1519 
  1520     have "\<lfloor>(a + b) * 2 powr q\<rfloor> = \<lfloor>(a + b) * 2 powr p * 2 powr (q - p)\<rfloor>"
  1521       by (simp add: algebra_simps flip: powr_realpow powr_add)
  1522     also have "\<dots> = \<lfloor>(ai + b * 2 powr p) * 2 powr (q - p)\<rfloor>"
  1523       by (simp add: assms algebra_simps)
  1524     also have "\<dots> = \<lfloor>(ai + b * 2 powr p) / real_of_int ((2::int) ^ nat (p - q))\<rfloor>"
  1525       using assms
  1526       by (simp add: algebra_simps divide_powr_uminus flip: powr_realpow powr_add)
  1527     also have "\<dots> = \<lfloor>ai / real_of_int ((2::int) ^ nat (p - q))\<rfloor>"
  1528       by (simp del: of_int_power add: floor_divide_real_eq_div floor_eq)
  1529     finally have "\<lfloor>(a + b) * 2 powr real_of_int q\<rfloor> = \<lfloor>real_of_int ai / real_of_int ((2::int) ^ nat (p - q))\<rfloor>" .
  1530     moreover
  1531     have "\<lfloor>(2 * ai + (sgn b)) * 2 powr (real_of_int (q - p) - 1)\<rfloor> =
  1532         \<lfloor>real_of_int ai / real_of_int ((2::int) ^ nat (p - q))\<rfloor>"
  1533     proof -
  1534       have "\<lfloor>(2 * ai + sgn b) * 2 powr (real_of_int (q - p) - 1)\<rfloor> = \<lfloor>(ai + sgn b / 2) * 2 powr (q - p)\<rfloor>"
  1535         by (subst powr_diff) (simp add: field_simps)
  1536       also have "\<dots> = \<lfloor>(ai + sgn b / 2) / real_of_int ((2::int) ^ nat (p - q))\<rfloor>"
  1537         using leqp by (simp flip: powr_realpow add: powr_diff)
  1538       also have "\<dots> = \<lfloor>ai / real_of_int ((2::int) ^ nat (p - q))\<rfloor>"
  1539         by (simp del: of_int_power add: floor_divide_real_eq_div floor_eq)
  1540       finally show ?thesis .
  1541     qed
  1542     ultimately show ?thesis by simp
  1543   next
  1544     case 3
  1545     then have floor_eq: "\<lfloor>b * 2 powr (real_of_int p + 1)\<rfloor> = -1"
  1546       using b_le_1
  1547       by (auto simp: floor_eq_iff algebra_simps pos_divide_le_eq[symmetric] abs_if divide_powr_uminus
  1548         intro!: mult_neg_pos split: if_split_asm)
  1549     have "\<lfloor>(a + b) * 2 powr q\<rfloor> = \<lfloor>(2*a + 2*b) * 2 powr p * 2 powr (q - p - 1)\<rfloor>"
  1550       by (simp add: algebra_simps powr_mult_base flip: powr_realpow powr_add)
  1551     also have "\<dots> = \<lfloor>(2 * (a * 2 powr p) + 2 * b * 2 powr p) * 2 powr (q - p - 1)\<rfloor>"
  1552       by (simp add: algebra_simps)
  1553     also have "\<dots> = \<lfloor>(2 * ai + b * 2 powr (p + 1)) / 2 powr (1 - q + p)\<rfloor>"
  1554       using assms by (simp add: algebra_simps powr_mult_base divide_powr_uminus)
  1555     also have "\<dots> = \<lfloor>(2 * ai + b * 2 powr (p + 1)) / real_of_int ((2::int) ^ nat (p - q + 1))\<rfloor>"
  1556       using assms by (simp add: algebra_simps flip: powr_realpow)
  1557     also have "\<dots> = \<lfloor>(2 * ai - 1) / real_of_int ((2::int) ^ nat (p - q + 1))\<rfloor>"
  1558       using \<open>b < 0\<close> assms
  1559       by (simp add: floor_divide_of_int_eq floor_eq floor_divide_real_eq_div
  1560         del: of_int_mult of_int_power of_int_diff)
  1561     also have "\<dots> = \<lfloor>(2 * ai - 1) * 2 powr (q - p - 1)\<rfloor>"
  1562       using assms by (simp add: algebra_simps divide_powr_uminus flip: powr_realpow)
  1563     finally show ?thesis
  1564       using \<open>b < 0\<close> by simp
  1565   qed
  1566 qed
  1567 
  1568 lemma log2_abs_int_add_less_half_sgn_eq:
  1569   fixes ai :: int
  1570     and b :: real
  1571   assumes "\<bar>b\<bar> \<le> 1/2"
  1572     and "ai \<noteq> 0"
  1573   shows "\<lfloor>log 2 \<bar>real_of_int ai + b\<bar>\<rfloor> = \<lfloor>log 2 \<bar>ai + sgn b / 2\<bar>\<rfloor>"
  1574 proof (cases "b = 0")
  1575   case True
  1576   then show ?thesis by simp
  1577 next
  1578   case False
  1579   define k where "k = \<lfloor>log 2 \<bar>ai\<bar>\<rfloor>"
  1580   then have "\<lfloor>log 2 \<bar>ai\<bar>\<rfloor> = k"
  1581     by simp
  1582   then have k: "2 powr k \<le> \<bar>ai\<bar>" "\<bar>ai\<bar> < 2 powr (k + 1)"
  1583     by (simp_all add: floor_log_eq_powr_iff \<open>ai \<noteq> 0\<close>)
  1584   have "k \<ge> 0"
  1585     using assms by (auto simp: k_def)
  1586   define r where "r = \<bar>ai\<bar> - 2 ^ nat k"
  1587   have r: "0 \<le> r" "r < 2 powr k"
  1588     using \<open>k \<ge> 0\<close> k
  1589     by (auto simp: r_def k_def algebra_simps powr_add abs_if powr_int)
  1590   then have "r \<le> (2::int) ^ nat k - 1"
  1591     using \<open>k \<ge> 0\<close> by (auto simp: powr_int)
  1592   from this[simplified of_int_le_iff[symmetric]] \<open>0 \<le> k\<close>
  1593   have r_le: "r \<le> 2 powr k - 1"
  1594     by (auto simp: algebra_simps powr_int)
  1595       (metis of_int_1 of_int_add of_int_le_numeral_power_cancel_iff)
  1596 
  1597   have "\<bar>ai\<bar> = 2 powr k + r"
  1598     using \<open>k \<ge> 0\<close> by (auto simp: k_def r_def simp flip: powr_realpow)
  1599 
  1600   have pos: "\<bar>b\<bar> < 1 \<Longrightarrow> 0 < 2 powr k + (r + b)" for b :: real
  1601     using \<open>0 \<le> k\<close> \<open>ai \<noteq> 0\<close>
  1602     by (auto simp add: r_def powr_realpow[symmetric] abs_if sgn_if algebra_simps
  1603       split: if_split_asm)
  1604   have less: "\<bar>sgn ai * b\<bar> < 1"
  1605     and less': "\<bar>sgn (sgn ai * b) / 2\<bar> < 1"
  1606     using \<open>\<bar>b\<bar> \<le> _\<close> by (auto simp: abs_if sgn_if split: if_split_asm)
  1607 
  1608   have floor_eq: "\<And>b::real. \<bar>b\<bar> \<le> 1 / 2 \<Longrightarrow>
  1609       \<lfloor>log 2 (1 + (r + b) / 2 powr k)\<rfloor> = (if r = 0 \<and> b < 0 then -1 else 0)"
  1610     using \<open>k \<ge> 0\<close> r r_le
  1611     by (auto simp: floor_log_eq_powr_iff powr_minus_divide field_simps sgn_if)
  1612 
  1613   from \<open>real_of_int \<bar>ai\<bar> = _\<close> have "\<bar>ai + b\<bar> = 2 powr k + (r + sgn ai * b)"
  1614     using \<open>\<bar>b\<bar> \<le> _\<close> \<open>0 \<le> k\<close> r
  1615     by (auto simp add: sgn_if abs_if)
  1616   also have "\<lfloor>log 2 \<dots>\<rfloor> = \<lfloor>log 2 (2 powr k + r + sgn (sgn ai * b) / 2)\<rfloor>"
  1617   proof -
  1618     have "2 powr k + (r + (sgn ai) * b) = 2 powr k * (1 + (r + sgn ai * b) / 2 powr k)"
  1619       by (simp add: field_simps)
  1620     also have "\<lfloor>log 2 \<dots>\<rfloor> = k + \<lfloor>log 2 (1 + (r + sgn ai * b) / 2 powr k)\<rfloor>"
  1621       using pos[OF less]
  1622       by (subst log_mult) (simp_all add: log_mult powr_mult field_simps)
  1623     also
  1624     let ?if = "if r = 0 \<and> sgn ai * b < 0 then -1 else 0"
  1625     have "\<lfloor>log 2 (1 + (r + sgn ai * b) / 2 powr k)\<rfloor> = ?if"
  1626       using \<open>\<bar>b\<bar> \<le> _\<close>
  1627       by (intro floor_eq) (auto simp: abs_mult sgn_if)
  1628     also
  1629     have "\<dots> = \<lfloor>log 2 (1 + (r + sgn (sgn ai * b) / 2) / 2 powr k)\<rfloor>"
  1630       by (subst floor_eq) (auto simp: sgn_if)
  1631     also have "k + \<dots> = \<lfloor>log 2 (2 powr k * (1 + (r + sgn (sgn ai * b) / 2) / 2 powr k))\<rfloor>"
  1632       unfolding int_add_floor
  1633       using pos[OF less'] \<open>\<bar>b\<bar> \<le> _\<close>
  1634       by (simp add: field_simps add_log_eq_powr del: floor_add2)
  1635     also have "2 powr k * (1 + (r + sgn (sgn ai * b) / 2) / 2 powr k) =
  1636         2 powr k + r + sgn (sgn ai * b) / 2"
  1637       by (simp add: sgn_if field_simps)
  1638     finally show ?thesis .
  1639   qed
  1640   also have "2 powr k + r + sgn (sgn ai * b) / 2 = \<bar>ai + sgn b / 2\<bar>"
  1641     unfolding \<open>real_of_int \<bar>ai\<bar> = _\<close>[symmetric] using \<open>ai \<noteq> 0\<close>
  1642     by (auto simp: abs_if sgn_if algebra_simps)
  1643   finally show ?thesis .
  1644 qed
  1645 
  1646 context
  1647 begin
  1648 
  1649 qualified lemma compute_far_float_plus_down:
  1650   fixes m1 e1 m2 e2 :: int
  1651     and p :: nat
  1652   defines "k1 \<equiv> Suc p - nat (bitlen \<bar>m1\<bar>)"
  1653   assumes H: "bitlen \<bar>m2\<bar> \<le> e1 - e2 - k1 - 2" "m1 \<noteq> 0" "m2 \<noteq> 0" "e1 \<ge> e2"
  1654   shows "float_plus_down p (Float m1 e1) (Float m2 e2) =
  1655     float_round_down p (Float (m1 * 2 ^ (Suc (Suc k1)) + sgn m2) (e1 - int k1 - 2))"
  1656 proof -
  1657   let ?a = "real_of_float (Float m1 e1)"
  1658   let ?b = "real_of_float (Float m2 e2)"
  1659   let ?sum = "?a + ?b"
  1660   let ?shift = "real_of_int e2 - real_of_int e1 + real k1 + 1"
  1661   let ?m1 = "m1 * 2 ^ Suc k1"
  1662   let ?m2 = "m2 * 2 powr ?shift"
  1663   let ?m2' = "sgn m2 / 2"
  1664   let ?e = "e1 - int k1 - 1"
  1665 
  1666   have sum_eq: "?sum = (?m1 + ?m2) * 2 powr ?e"
  1667     by (auto simp flip: powr_add powr_mult powr_realpow simp: powr_mult_base algebra_simps)
  1668 
  1669   have "\<bar>?m2\<bar> * 2 < 2 powr (bitlen \<bar>m2\<bar> + ?shift + 1)"
  1670     by (auto simp: field_simps powr_add powr_mult_base powr_diff abs_mult)
  1671   also have "\<dots> \<le> 2 powr 0"
  1672     using H by (intro powr_mono) auto
  1673   finally have abs_m2_less_half: "\<bar>?m2\<bar> < 1 / 2"
  1674     by simp
  1675 
  1676   then have "\<bar>real_of_int m2\<bar> < 2 powr -(?shift + 1)"
  1677     unfolding powr_minus_divide by (auto simp: bitlen_alt_def field_simps powr_mult_base abs_mult)
  1678   also have "\<dots> \<le> 2 powr real_of_int (e1 - e2 - 2)"
  1679     by simp
  1680   finally have b_less_quarter: "\<bar>?b\<bar> < 1/4 * 2 powr real_of_int e1"
  1681     by (simp add: powr_add field_simps powr_diff abs_mult)
  1682   also have "1/4 < \<bar>real_of_int m1\<bar> / 2" using \<open>m1 \<noteq> 0\<close> by simp
  1683   finally have b_less_half_a: "\<bar>?b\<bar> < 1/2 * \<bar>?a\<bar>"
  1684     by (simp add: algebra_simps powr_mult_base abs_mult)
  1685   then have a_half_less_sum: "\<bar>?a\<bar> / 2 < \<bar>?sum\<bar>"
  1686     by (auto simp: field_simps abs_if split: if_split_asm)
  1687 
  1688   from b_less_half_a have "\<bar>?b\<bar> < \<bar>?a\<bar>" "\<bar>?b\<bar> \<le> \<bar>?a\<bar>"
  1689     by simp_all
  1690 
  1691   have "\<bar>real_of_float (Float m1 e1)\<bar> \<ge> 1/4 * 2 powr real_of_int e1"
  1692     using \<open>m1 \<noteq> 0\<close>
  1693     by (auto simp: powr_add powr_int bitlen_nonneg divide_right_mono abs_mult)
  1694   then have "?sum \<noteq> 0" using b_less_quarter
  1695     by (rule sum_neq_zeroI)
  1696   then have "?m1 + ?m2 \<noteq> 0"
  1697     unfolding sum_eq by (simp add: abs_mult zero_less_mult_iff)
  1698 
  1699   have "\<bar>real_of_int ?m1\<bar> \<ge> 2 ^ Suc k1" "\<bar>?m2'\<bar> < 2 ^ Suc k1"
  1700     using \<open>m1 \<noteq> 0\<close> \<open>m2 \<noteq> 0\<close> by (auto simp: sgn_if less_1_mult abs_mult simp del: power.simps)
  1701   then have sum'_nz: "?m1 + ?m2' \<noteq> 0"
  1702     by (intro sum_neq_zeroI)
  1703 
  1704   have "\<lfloor>log 2 \<bar>real_of_float (Float m1 e1) + real_of_float (Float m2 e2)\<bar>\<rfloor> = \<lfloor>log 2 \<bar>?m1 + ?m2\<bar>\<rfloor> + ?e"
  1705     using \<open>?m1 + ?m2 \<noteq> 0\<close>
  1706     unfolding floor_add[symmetric] sum_eq
  1707     by (simp add: abs_mult log_mult) linarith
  1708   also have "\<lfloor>log 2 \<bar>?m1 + ?m2\<bar>\<rfloor> = \<lfloor>log 2 \<bar>?m1 + sgn (real_of_int m2 * 2 powr ?shift) / 2\<bar>\<rfloor>"
  1709     using abs_m2_less_half \<open>m1 \<noteq> 0\<close>
  1710     by (intro log2_abs_int_add_less_half_sgn_eq) (auto simp: abs_mult)
  1711   also have "sgn (real_of_int m2 * 2 powr ?shift) = sgn m2"
  1712     by (auto simp: sgn_if zero_less_mult_iff less_not_sym)
  1713   also
  1714   have "\<bar>?m1 + ?m2'\<bar> * 2 powr ?e = \<bar>?m1 * 2 + sgn m2\<bar> * 2 powr (?e - 1)"
  1715     by (auto simp: field_simps powr_minus[symmetric] powr_diff powr_mult_base)
  1716   then have "\<lfloor>log 2 \<bar>?m1 + ?m2'\<bar>\<rfloor> + ?e = \<lfloor>log 2 \<bar>real_of_float (Float (?m1 * 2 + sgn m2) (?e - 1))\<bar>\<rfloor>"
  1717     using \<open>?m1 + ?m2' \<noteq> 0\<close>
  1718     unfolding floor_add_int
  1719     by (simp add: log_add_eq_powr abs_mult_pos del: floor_add2)
  1720   finally
  1721   have "\<lfloor>log 2 \<bar>?sum\<bar>\<rfloor> = \<lfloor>log 2 \<bar>real_of_float (Float (?m1*2 + sgn m2) (?e - 1))\<bar>\<rfloor>" .
  1722   then have "plus_down p (Float m1 e1) (Float m2 e2) =
  1723       truncate_down p (Float (?m1*2 + sgn m2) (?e - 1))"
  1724     unfolding plus_down_def
  1725   proof (rule truncate_down_log2_eqI)
  1726     let ?f = "(int p - \<lfloor>log 2 \<bar>real_of_float (Float m1 e1) + real_of_float (Float m2 e2)\<bar>\<rfloor>)"
  1727     let ?ai = "m1 * 2 ^ (Suc k1)"
  1728     have "\<lfloor>(?a + ?b) * 2 powr real_of_int ?f\<rfloor> = \<lfloor>(real_of_int (2 * ?ai) + sgn ?b) * 2 powr real_of_int (?f - - ?e - 1)\<rfloor>"
  1729     proof (rule floor_sum_times_2_powr_sgn_eq)
  1730       show "?a * 2 powr real_of_int (-?e) = real_of_int ?ai"
  1731         by (simp add: powr_add powr_realpow[symmetric] powr_diff)
  1732       show "\<bar>?b * 2 powr real_of_int (-?e + 1)\<bar> \<le> 1"
  1733         using abs_m2_less_half
  1734         by (simp add: abs_mult powr_add[symmetric] algebra_simps powr_mult_base)
  1735     next
  1736       have "e1 + \<lfloor>log 2 \<bar>real_of_int m1\<bar>\<rfloor> - 1 = \<lfloor>log 2 \<bar>?a\<bar>\<rfloor> - 1"
  1737         using \<open>m1 \<noteq> 0\<close>
  1738         by (simp add: int_add_floor algebra_simps log_mult abs_mult del: floor_add2)
  1739       also have "\<dots> \<le> \<lfloor>log 2 \<bar>?a + ?b\<bar>\<rfloor>"
  1740         using a_half_less_sum \<open>m1 \<noteq> 0\<close> \<open>?sum \<noteq> 0\<close>
  1741         unfolding floor_diff_of_int[symmetric]
  1742         by (auto simp add: log_minus_eq_powr powr_minus_divide intro!: floor_mono)
  1743       finally
  1744       have "int p - \<lfloor>log 2 \<bar>?a + ?b\<bar>\<rfloor> \<le> p - (bitlen \<bar>m1\<bar>) - e1 + 2"
  1745         by (auto simp: algebra_simps bitlen_alt_def \<open>m1 \<noteq> 0\<close>)
  1746       also have "\<dots> \<le> - ?e"
  1747         using bitlen_nonneg[of "\<bar>m1\<bar>"] by (simp add: k1_def)
  1748       finally show "?f \<le> - ?e" by simp
  1749     qed
  1750     also have "sgn ?b = sgn m2"
  1751       using powr_gt_zero[of 2 e2]
  1752       by (auto simp add: sgn_if zero_less_mult_iff simp del: powr_gt_zero)
  1753     also have "\<lfloor>(real_of_int (2 * ?m1) + real_of_int (sgn m2)) * 2 powr real_of_int (?f - - ?e - 1)\<rfloor> =
  1754         \<lfloor>Float (?m1 * 2 + sgn m2) (?e - 1) * 2 powr ?f\<rfloor>"
  1755       by (simp flip: powr_add powr_realpow add: algebra_simps)
  1756     finally
  1757     show "\<lfloor>(?a + ?b) * 2 powr ?f\<rfloor> = \<lfloor>real_of_float (Float (?m1 * 2 + sgn m2) (?e - 1)) * 2 powr ?f\<rfloor>" .
  1758   qed
  1759   then show ?thesis
  1760     by transfer (simp add: plus_down_def ac_simps Let_def)
  1761 qed
  1762 
  1763 lemma compute_float_plus_down_naive[code]: "float_plus_down p x y = float_round_down p (x + y)"
  1764   by transfer (auto simp: plus_down_def)
  1765 
  1766 qualified lemma compute_float_plus_down[code]:
  1767   fixes p::nat and m1 e1 m2 e2::int
  1768   shows "float_plus_down p (Float m1 e1) (Float m2 e2) =
  1769     (if m1 = 0 then float_round_down p (Float m2 e2)
  1770     else if m2 = 0 then float_round_down p (Float m1 e1)
  1771     else
  1772       (if e1 \<ge> e2 then
  1773         (let k1 = Suc p - nat (bitlen \<bar>m1\<bar>) in
  1774           if bitlen \<bar>m2\<bar> > e1 - e2 - k1 - 2
  1775           then float_round_down p ((Float m1 e1) + (Float m2 e2))
  1776           else float_round_down p (Float (m1 * 2 ^ (Suc (Suc k1)) + sgn m2) (e1 - int k1 - 2)))
  1777     else float_plus_down p (Float m2 e2) (Float m1 e1)))"
  1778 proof -
  1779   {
  1780     assume "bitlen \<bar>m2\<bar> \<le> e1 - e2 - (Suc p - nat (bitlen \<bar>m1\<bar>)) - 2" "m1 \<noteq> 0" "m2 \<noteq> 0" "e1 \<ge> e2"
  1781     note compute_far_float_plus_down[OF this]
  1782   }
  1783   then show ?thesis
  1784     by transfer (simp add: Let_def plus_down_def ac_simps)
  1785 qed
  1786 
  1787 qualified lemma compute_float_plus_up[code]: "float_plus_up p x y = - float_plus_down p (-x) (-y)"
  1788   using truncate_down_uminus_eq[of p "x + y"]
  1789   by transfer (simp add: plus_down_def plus_up_def ac_simps)
  1790 
  1791 lemma mantissa_zero[simp]: "mantissa 0 = 0"
  1792   by (metis mantissa_0 zero_float.abs_eq)
  1793 
  1794 qualified lemma compute_float_less[code]: "a < b \<longleftrightarrow> is_float_pos (float_plus_down 0 b (- a))"
  1795   using truncate_down[of 0 "b - a"] truncate_down_pos[of "b - a" 0]
  1796   by transfer (auto simp: plus_down_def)
  1797 
  1798 qualified lemma compute_float_le[code]: "a \<le> b \<longleftrightarrow> is_float_nonneg (float_plus_down 0 b (- a))"
  1799   using truncate_down[of 0 "b - a"] truncate_down_nonneg[of "b - a" 0]
  1800   by transfer (auto simp: plus_down_def)
  1801 
  1802 end
  1803 
  1804 
  1805 subsection \<open>Lemmas needed by Approximate\<close>
  1806 
  1807 lemma Float_num[simp]:
  1808    "real_of_float (Float 1 0) = 1"
  1809    "real_of_float (Float 1 1) = 2"
  1810    "real_of_float (Float 1 2) = 4"
  1811    "real_of_float (Float 1 (- 1)) = 1/2"
  1812    "real_of_float (Float 1 (- 2)) = 1/4"
  1813    "real_of_float (Float 1 (- 3)) = 1/8"
  1814    "real_of_float (Float (- 1) 0) = -1"
  1815    "real_of_float (Float (numeral n) 0) = numeral n"
  1816    "real_of_float (Float (- numeral n) 0) = - numeral n"
  1817   using two_powr_int_float[of 2] two_powr_int_float[of "-1"] two_powr_int_float[of "-2"]
  1818     two_powr_int_float[of "-3"]
  1819   using powr_realpow[of 2 2] powr_realpow[of 2 3]
  1820   using powr_minus[of "2::real" 1] powr_minus[of "2::real" 2] powr_minus[of "2::real" 3]
  1821   by auto
  1822 
  1823 lemma real_of_Float_int[simp]: "real_of_float (Float n 0) = real n"
  1824   by simp
  1825 
  1826 lemma float_zero[simp]: "real_of_float (Float 0 e) = 0"
  1827   by simp
  1828 
  1829 lemma abs_div_2_less: "a \<noteq> 0 \<Longrightarrow> a \<noteq> -1 \<Longrightarrow> \<bar>(a::int) div 2\<bar> < \<bar>a\<bar>"
  1830   by arith
  1831 
  1832 lemma lapprox_rat: "real_of_float (lapprox_rat prec x y) \<le> real_of_int x / real_of_int y"
  1833   by (simp add: lapprox_rat.rep_eq truncate_down)
  1834 
  1835 lemma mult_div_le:
  1836   fixes a b :: int
  1837   assumes "b > 0"
  1838   shows "a \<ge> b * (a div b)"
  1839 proof -
  1840   from minus_div_mult_eq_mod [symmetric, of a b] have "a = b * (a div b) + a mod b"
  1841     by simp
  1842   also have "\<dots> \<ge> b * (a div b) + 0"
  1843     apply (rule add_left_mono)
  1844     apply (rule pos_mod_sign)
  1845     using assms
  1846     apply simp
  1847     done
  1848   finally show ?thesis
  1849     by simp
  1850 qed
  1851 
  1852 lemma lapprox_rat_nonneg:
  1853   assumes "0 \<le> x" and "0 \<le> y"
  1854   shows "0 \<le> real_of_float (lapprox_rat n x y)"
  1855   using assms
  1856   by transfer (simp add: truncate_down_nonneg)
  1857 
  1858 lemma rapprox_rat: "real_of_int x / real_of_int y \<le> real_of_float (rapprox_rat prec x y)"
  1859   by transfer (simp add: truncate_up)
  1860 
  1861 lemma rapprox_rat_le1:
  1862   assumes "0 \<le> x" "0 < y" "x \<le> y"
  1863   shows "real_of_float (rapprox_rat n x y) \<le> 1"
  1864   using assms
  1865   by transfer (simp add: truncate_up_le1)
  1866 
  1867 lemma rapprox_rat_nonneg_nonpos: "0 \<le> x \<Longrightarrow> y \<le> 0 \<Longrightarrow> real_of_float (rapprox_rat n x y) \<le> 0"
  1868   by transfer (simp add: truncate_up_nonpos divide_nonneg_nonpos)
  1869 
  1870 lemma rapprox_rat_nonpos_nonneg: "x \<le> 0 \<Longrightarrow> 0 \<le> y \<Longrightarrow> real_of_float (rapprox_rat n x y) \<le> 0"
  1871   by transfer (simp add: truncate_up_nonpos divide_nonpos_nonneg)
  1872 
  1873 lemma real_divl: "real_divl prec x y \<le> x / y"
  1874   by (simp add: real_divl_def truncate_down)
  1875 
  1876 lemma real_divr: "x / y \<le> real_divr prec x y"
  1877   by (simp add: real_divr_def truncate_up)
  1878 
  1879 lemma float_divl: "real_of_float (float_divl prec x y) \<le> x / y"
  1880   by transfer (rule real_divl)
  1881 
  1882 lemma real_divl_lower_bound: "0 \<le> x \<Longrightarrow> 0 \<le> y \<Longrightarrow> 0 \<le> real_divl prec x y"
  1883   by (simp add: real_divl_def truncate_down_nonneg)
  1884 
  1885 lemma float_divl_lower_bound: "0 \<le> x \<Longrightarrow> 0 \<le> y \<Longrightarrow> 0 \<le> real_of_float (float_divl prec x y)"
  1886   by transfer (rule real_divl_lower_bound)
  1887 
  1888 lemma exponent_1: "exponent 1 = 0"
  1889   using exponent_float[of 1 0] by (simp add: one_float_def)
  1890 
  1891 lemma mantissa_1: "mantissa 1 = 1"
  1892   using mantissa_float[of 1 0] by (simp add: one_float_def)
  1893 
  1894 lemma bitlen_1: "bitlen 1 = 1"
  1895   by (simp add: bitlen_alt_def)
  1896 
  1897 lemma float_upper_bound: "x \<le> 2 powr (bitlen \<bar>mantissa x\<bar> + exponent x)"
  1898 proof (cases "x = 0")
  1899   case True
  1900   then show ?thesis by simp
  1901 next
  1902   case False
  1903   then have "mantissa x \<noteq> 0"
  1904     using mantissa_eq_zero_iff by auto
  1905   have "x = mantissa x * 2 powr (exponent x)"
  1906     by (rule mantissa_exponent)
  1907   also have "mantissa x \<le> \<bar>mantissa x\<bar>"
  1908     by simp
  1909   also have "\<dots> \<le> 2 powr (bitlen \<bar>mantissa x\<bar>)"
  1910     using bitlen_bounds[of "\<bar>mantissa x\<bar>"] bitlen_nonneg \<open>mantissa x \<noteq> 0\<close>
  1911     by (auto simp del: of_int_abs simp add: powr_int)
  1912   finally show ?thesis by (simp add: powr_add)
  1913 qed
  1914 
  1915 lemma real_divl_pos_less1_bound:
  1916   assumes "0 < x" "x \<le> 1"
  1917   shows "1 \<le> real_divl prec 1 x"
  1918   using assms
  1919   by (auto intro!: truncate_down_ge1 simp: real_divl_def)
  1920 
  1921 lemma float_divl_pos_less1_bound:
  1922   "0 < real_of_float x \<Longrightarrow> real_of_float x \<le> 1 \<Longrightarrow> prec \<ge> 1 \<Longrightarrow>
  1923     1 \<le> real_of_float (float_divl prec 1 x)"
  1924   by transfer (rule real_divl_pos_less1_bound)
  1925 
  1926 lemma float_divr: "real_of_float x / real_of_float y \<le> real_of_float (float_divr prec x y)"
  1927   by transfer (rule real_divr)
  1928 
  1929 lemma real_divr_pos_less1_lower_bound:
  1930   assumes "0 < x"
  1931     and "x \<le> 1"
  1932   shows "1 \<le> real_divr prec 1 x"
  1933 proof -
  1934   have "1 \<le> 1 / x"
  1935     using \<open>0 < x\<close> and \<open>x \<le> 1\<close> by auto
  1936   also have "\<dots> \<le> real_divr prec 1 x"
  1937     using real_divr[where x = 1 and y = x] by auto
  1938   finally show ?thesis by auto
  1939 qed
  1940 
  1941 lemma float_divr_pos_less1_lower_bound: "0 < x \<Longrightarrow> x \<le> 1 \<Longrightarrow> 1 \<le> float_divr prec 1 x"
  1942   by transfer (rule real_divr_pos_less1_lower_bound)
  1943 
  1944 lemma real_divr_nonpos_pos_upper_bound: "x \<le> 0 \<Longrightarrow> 0 \<le> y \<Longrightarrow> real_divr prec x y \<le> 0"
  1945   by (simp add: real_divr_def truncate_up_nonpos divide_le_0_iff)
  1946 
  1947 lemma float_divr_nonpos_pos_upper_bound:
  1948   "real_of_float x \<le> 0 \<Longrightarrow> 0 \<le> real_of_float y \<Longrightarrow> real_of_float (float_divr prec x y) \<le> 0"
  1949   by transfer (rule real_divr_nonpos_pos_upper_bound)
  1950 
  1951 lemma real_divr_nonneg_neg_upper_bound: "0 \<le> x \<Longrightarrow> y \<le> 0 \<Longrightarrow> real_divr prec x y \<le> 0"
  1952   by (simp add: real_divr_def truncate_up_nonpos divide_le_0_iff)
  1953 
  1954 lemma float_divr_nonneg_neg_upper_bound:
  1955   "0 \<le> real_of_float x \<Longrightarrow> real_of_float y \<le> 0 \<Longrightarrow> real_of_float (float_divr prec x y) \<le> 0"
  1956   by transfer (rule real_divr_nonneg_neg_upper_bound)
  1957 
  1958 lemma truncate_up_nonneg_mono:
  1959   assumes "0 \<le> x" "x \<le> y"
  1960   shows "truncate_up prec x \<le> truncate_up prec y"
  1961 proof -
  1962   consider "\<lfloor>log 2 x\<rfloor> = \<lfloor>log 2 y\<rfloor>" | "\<lfloor>log 2 x\<rfloor> \<noteq> \<lfloor>log 2 y\<rfloor>" "0 < x" | "x \<le> 0"
  1963     by arith
  1964   then show ?thesis
  1965   proof cases
  1966     case 1
  1967     then show ?thesis
  1968       using assms
  1969       by (auto simp: truncate_up_def round_up_def intro!: ceiling_mono)
  1970   next
  1971     case 2
  1972     from assms \<open>0 < x\<close> have "log 2 x \<le> log 2 y"
  1973       by auto
  1974     with \<open>\<lfloor>log 2 x\<rfloor> \<noteq> \<lfloor>log 2 y\<rfloor>\<close>
  1975     have logless: "log 2 x < log 2 y" and flogless: "\<lfloor>log 2 x\<rfloor> < \<lfloor>log 2 y\<rfloor>"
  1976       by (metis floor_less_cancel linorder_cases not_le)+
  1977     have "truncate_up prec x =
  1978       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>)"
  1979       using assms by (simp add: truncate_up_def round_up_def)
  1980     also have "\<lceil>x * 2 powr real_of_int (int prec - \<lfloor>log 2 x\<rfloor>)\<rceil> \<le> (2 ^ (Suc prec))"
  1981     proof (simp only: ceiling_le_iff)
  1982       have "x * 2 powr real_of_int (int prec - \<lfloor>log 2 x\<rfloor>) \<le>
  1983         x * (2 powr real (Suc prec) / (2 powr log 2 x))"
  1984         using real_of_int_floor_add_one_ge[of "log 2 x"] assms
  1985         by (auto simp: algebra_simps simp flip: powr_diff intro!: mult_left_mono)
  1986       then show "x * 2 powr real_of_int (int prec - \<lfloor>log 2 x\<rfloor>) \<le> real_of_int ((2::int) ^ (Suc prec))"
  1987         using \<open>0 < x\<close> by (simp add: powr_realpow powr_add)
  1988     qed
  1989     then have "real_of_int \<lceil>x * 2 powr real_of_int (int prec - \<lfloor>log 2 x\<rfloor>)\<rceil> \<le> 2 powr int (Suc prec)"
  1990       by (auto simp: powr_realpow powr_add)
  1991         (metis power_Suc of_int_le_numeral_power_cancel_iff)
  1992     also
  1993     have "2 powr - real_of_int (int prec - \<lfloor>log 2 x\<rfloor>) \<le> 2 powr - real_of_int (int prec - \<lfloor>log 2 y\<rfloor> + 1)"
  1994       using logless flogless by (auto intro!: floor_mono)
  1995     also have "2 powr real_of_int (int (Suc prec)) \<le>
  1996         2 powr (log 2 y + real_of_int (int prec - \<lfloor>log 2 y\<rfloor> + 1))"
  1997       using assms \<open>0 < x\<close>
  1998       by (auto simp: algebra_simps)
  1999     finally have "truncate_up prec x \<le>
  2000         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)"
  2001       by simp
  2002     also have "\<dots> = 2 powr (log 2 y + real_of_int (int prec - \<lfloor>log 2 y\<rfloor>) - real_of_int (int prec - \<lfloor>log 2 y\<rfloor>))"
  2003       by (subst powr_add[symmetric]) simp
  2004     also have "\<dots> = y"
  2005       using \<open>0 < x\<close> assms
  2006       by (simp add: powr_add)
  2007     also have "\<dots> \<le> truncate_up prec y"
  2008       by (rule truncate_up)
  2009     finally show ?thesis .
  2010   next
  2011     case 3
  2012     then show ?thesis
  2013       using assms
  2014       by (auto intro!: truncate_up_le)
  2015   qed
  2016 qed
  2017 
  2018 lemma truncate_up_switch_sign_mono:
  2019   assumes "x \<le> 0" "0 \<le> y"
  2020   shows "truncate_up prec x \<le> truncate_up prec y"
  2021 proof -
  2022   note truncate_up_nonpos[OF \<open>x \<le> 0\<close>]
  2023   also note truncate_up_le[OF \<open>0 \<le> y\<close>]
  2024   finally show ?thesis .
  2025 qed
  2026 
  2027 lemma truncate_down_switch_sign_mono:
  2028   assumes "x \<le> 0"
  2029     and "0 \<le> y"
  2030     and "x \<le> y"
  2031   shows "truncate_down prec x \<le> truncate_down prec y"
  2032 proof -
  2033   note truncate_down_le[OF \<open>x \<le> 0\<close>]
  2034   also note truncate_down_nonneg[OF \<open>0 \<le> y\<close>]
  2035   finally show ?thesis .
  2036 qed
  2037 
  2038 lemma truncate_down_nonneg_mono:
  2039   assumes "0 \<le> x" "x \<le> y"
  2040   shows "truncate_down prec x \<le> truncate_down prec y"
  2041 proof -
  2042   consider "x \<le> 0" | "\<lfloor>log 2 \<bar>x\<bar>\<rfloor> = \<lfloor>log 2 \<bar>y\<bar>\<rfloor>" |
  2043     "0 < x" "\<lfloor>log 2 \<bar>x\<bar>\<rfloor> \<noteq> \<lfloor>log 2 \<bar>y\<bar>\<rfloor>"
  2044     by arith
  2045   then show ?thesis
  2046   proof cases
  2047     case 1
  2048     with assms have "x = 0" "0 \<le> y" by simp_all
  2049     then show ?thesis
  2050       by (auto intro!: truncate_down_nonneg)
  2051   next
  2052     case 2
  2053     then show ?thesis
  2054       using assms
  2055       by (auto simp: truncate_down_def round_down_def intro!: floor_mono)
  2056   next
  2057     case 3
  2058     from \<open>0 < x\<close> have "log 2 x \<le> log 2 y" "0 < y" "0 \<le> y"
  2059       using assms by auto
  2060     with \<open>\<lfloor>log 2 \<bar>x\<bar>\<rfloor> \<noteq> \<lfloor>log 2 \<bar>y\<bar>\<rfloor>\<close>
  2061     have logless: "log 2 x < log 2 y" and flogless: "\<lfloor>log 2 x\<rfloor> < \<lfloor>log 2 y\<rfloor>"
  2062       unfolding atomize_conj abs_of_pos[OF \<open>0 < x\<close>] abs_of_pos[OF \<open>0 < y\<close>]
  2063       by (metis floor_less_cancel linorder_cases not_le)
  2064     have "2 powr prec \<le> y * 2 powr real prec / (2 powr log 2 y)"
  2065       using \<open>0 < y\<close> by simp
  2066     also have "\<dots> \<le> y * 2 powr real (Suc prec) / (2 powr (real_of_int \<lfloor>log 2 y\<rfloor> + 1))"
  2067       using \<open>0 \<le> y\<close> \<open>0 \<le> x\<close> assms(2)
  2068       by (auto intro!: powr_mono divide_left_mono
  2069           simp: of_nat_diff powr_add powr_diff)
  2070     also have "\<dots> = y * 2 powr real (Suc prec) / (2 powr real_of_int \<lfloor>log 2 y\<rfloor> * 2)"
  2071       by (auto simp: powr_add)
  2072     finally have "(2 ^ prec) \<le> \<lfloor>y * 2 powr real_of_int (int (Suc prec) - \<lfloor>log 2 \<bar>y\<bar>\<rfloor> - 1)\<rfloor>"
  2073       using \<open>0 \<le> y\<close>
  2074       by (auto simp: powr_diff le_floor_iff powr_realpow powr_add)
  2075     then have "(2 ^ (prec)) * 2 powr - real_of_int (int prec - \<lfloor>log 2 \<bar>y\<bar>\<rfloor>) \<le> truncate_down prec y"
  2076       by (auto simp: truncate_down_def round_down_def)
  2077     moreover have "x \<le> (2 ^ prec) * 2 powr - real_of_int (int prec - \<lfloor>log 2 \<bar>y\<bar>\<rfloor>)"
  2078     proof -
  2079       have "x = 2 powr (log 2 \<bar>x\<bar>)" using \<open>0 < x\<close> by simp
  2080       also have "\<dots> \<le> (2 ^ (Suc prec )) * 2 powr - real_of_int (int prec - \<lfloor>log 2 \<bar>x\<bar>\<rfloor>)"
  2081         using real_of_int_floor_add_one_ge[of "log 2 \<bar>x\<bar>"] \<open>0 < x\<close>
  2082         by (auto simp flip: powr_realpow powr_add simp: algebra_simps powr_mult_base le_powr_iff)
  2083       also
  2084       have "2 powr - real_of_int (int prec - \<lfloor>log 2 \<bar>x\<bar>\<rfloor>) \<le> 2 powr - real_of_int (int prec - \<lfloor>log 2 \<bar>y\<bar>\<rfloor> + 1)"
  2085         using logless flogless \<open>x > 0\<close> \<open>y > 0\<close>
  2086         by (auto intro!: floor_mono)
  2087       finally show ?thesis
  2088         by (auto simp flip: powr_realpow simp: powr_diff assms of_nat_diff)
  2089     qed
  2090     ultimately show ?thesis
  2091       by (metis dual_order.trans truncate_down)
  2092   qed
  2093 qed
  2094 
  2095 lemma truncate_down_eq_truncate_up: "truncate_down p x = - truncate_up p (-x)"
  2096   and truncate_up_eq_truncate_down: "truncate_up p x = - truncate_down p (-x)"
  2097   by (auto simp: truncate_up_uminus_eq truncate_down_uminus_eq)
  2098 
  2099 lemma truncate_down_mono: "x \<le> y \<Longrightarrow> truncate_down p x \<le> truncate_down p y"
  2100   apply (cases "0 \<le> x")
  2101   apply (rule truncate_down_nonneg_mono, assumption+)
  2102   apply (simp add: truncate_down_eq_truncate_up)
  2103   apply (cases "0 \<le> y")
  2104   apply (auto intro: truncate_up_nonneg_mono truncate_up_switch_sign_mono)
  2105   done
  2106 
  2107 lemma truncate_up_mono: "x \<le> y \<Longrightarrow> truncate_up p x \<le> truncate_up p y"
  2108   by (simp add: truncate_up_eq_truncate_down truncate_down_mono)
  2109 
  2110 lemma Float_le_zero_iff: "Float a b \<le> 0 \<longleftrightarrow> a \<le> 0"
  2111   by (auto simp: zero_float_def mult_le_0_iff)
  2112 
  2113 lemma real_of_float_pprt[simp]:
  2114   fixes a :: float
  2115   shows "real_of_float (pprt a) = pprt (real_of_float a)"
  2116   unfolding pprt_def sup_float_def max_def sup_real_def by auto
  2117 
  2118 lemma real_of_float_nprt[simp]:
  2119   fixes a :: float
  2120   shows "real_of_float (nprt a) = nprt (real_of_float a)"
  2121   unfolding nprt_def inf_float_def min_def inf_real_def by auto
  2122 
  2123 context
  2124 begin
  2125 
  2126 lift_definition int_floor_fl :: "float \<Rightarrow> int" is floor .
  2127 
  2128 qualified lemma compute_int_floor_fl[code]:
  2129   "int_floor_fl (Float m e) = (if 0 \<le> e then m * 2 ^ nat e else m div (2 ^ (nat (-e))))"
  2130   apply transfer
  2131   apply (simp add: powr_int floor_divide_of_int_eq)
  2132   apply (metis (no_types, hide_lams)floor_divide_of_int_eq of_int_numeral of_int_power floor_of_int of_int_mult)
  2133   done
  2134 
  2135 lift_definition floor_fl :: "float \<Rightarrow> float" is "\<lambda>x. real_of_int \<lfloor>x\<rfloor>"
  2136   by simp
  2137 
  2138 qualified lemma compute_floor_fl[code]:
  2139   "floor_fl (Float m e) = (if 0 \<le> e then Float m e else Float (m div (2 ^ (nat (-e)))) 0)"
  2140   apply transfer
  2141   apply (simp add: powr_int floor_divide_of_int_eq)
  2142   apply (metis (no_types, hide_lams)floor_divide_of_int_eq of_int_numeral of_int_power of_int_mult)
  2143   done
  2144 
  2145 end
  2146 
  2147 lemma floor_fl: "real_of_float (floor_fl x) \<le> real_of_float x"
  2148   by transfer simp
  2149 
  2150 lemma int_floor_fl: "real_of_int (int_floor_fl x) \<le> real_of_float x"
  2151   by transfer simp
  2152 
  2153 lemma floor_pos_exp: "exponent (floor_fl x) \<ge> 0"
  2154 proof (cases "floor_fl x = 0")
  2155   case True
  2156   then show ?thesis
  2157     by (simp add: floor_fl_def)
  2158 next
  2159   case False
  2160   have eq: "floor_fl x = Float \<lfloor>real_of_float x\<rfloor> 0"
  2161     by transfer simp
  2162   obtain i where "\<lfloor>real_of_float x\<rfloor> = mantissa (floor_fl x) * 2 ^ i" "0 = exponent (floor_fl x) - int i"
  2163     by (rule denormalize_shift[OF eq False])
  2164   then show ?thesis
  2165     by simp
  2166 qed
  2167 
  2168 lemma compute_mantissa[code]:
  2169   "mantissa (Float m e) =
  2170     (if m = 0 then 0 else if 2 dvd m then mantissa (normfloat (Float m e)) else m)"
  2171   by (auto simp: mantissa_float Float.abs_eq simp flip: zero_float_def)
  2172 
  2173 lemma compute_exponent[code]:
  2174   "exponent (Float m e) =
  2175     (if m = 0 then 0 else if 2 dvd m then exponent (normfloat (Float m e)) else e)"
  2176   by (auto simp: exponent_float Float.abs_eq simp flip: zero_float_def)
  2177 
  2178 lifting_update Float.float.lifting
  2179 lifting_forget Float.float.lifting
  2180 
  2181 end