src/HOL/Library/Float.thy
author haftmann
Mon Jun 05 15:59:41 2017 +0200 (2017-06-05)
changeset 66010 2f7d39285a1a
parent 65583 8d53b3bebab4
child 66912 a99a7cbf0fb5
permissions -rw-r--r--
executable domain membership checks
     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_realpow[symmetric] powr_diff field_simps)
    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 add: powr_realpow[symmetric] field_simps powr_add[symmetric])
   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 add: powr_realpow[symmetric] field_simps powr_add[symmetric])
   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 "op +" by simp
   203 declare plus_float.rep_eq[simp]
   204 
   205 lift_definition times_float :: "float \<Rightarrow> float \<Rightarrow> float" is "op *" by simp
   206 declare times_float.rep_eq[simp]
   207 
   208 lift_definition minus_float :: "float \<Rightarrow> float \<Rightarrow> float" is "op -" 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 "op = :: real \<Rightarrow> real \<Rightarrow> bool" .
   221 
   222 lift_definition less_eq_float :: "float \<Rightarrow> float \<Rightarrow> bool" is "op \<le>" .
   223 declare less_eq_float.rep_eq[simp]
   224 
   225 lift_definition less_float :: "float \<Rightarrow> float \<Rightarrow> bool" is "op <" .
   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 (op =) 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 (op =) 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[simp]: "numeral k = float_of (numeral k)"
   304   and float_of_neg_numeral[simp]: "- 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 (float_of 0) = 0" (is ?E)
   449   and mantissa_0[simp]: "mantissa (float_of 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> (float_of 0) \<Longrightarrow> \<not> 2 dvd mantissa f" (is "_ \<Longrightarrow> ?D")
   459 proof cases
   460   assume [simp]: "f \<noteq> float_of 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 add: zero_float_def)
   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 add: zero_float_def)
   473   qed
   474   then show ?E ?D by auto
   475 qed simp
   476 
   477 lemma mantissa_noteq_0: "f \<noteq> float_of 0 \<Longrightarrow> mantissa f \<noteq> 0"
   478   using mantissa_not_dvd[of f] by auto
   479 
   480 lemma
   481   fixes m e :: int
   482   defines "f \<equiv> float_of (m * 2 powr e)"
   483   assumes dvd: "\<not> 2 dvd m"
   484   shows mantissa_float: "mantissa f = m" (is "?M")
   485     and exponent_float: "m \<noteq> 0 \<Longrightarrow> exponent f = e" (is "_ \<Longrightarrow> ?E")
   486 proof cases
   487   assume "m = 0"
   488   with dvd show "mantissa f = m" by auto
   489 next
   490   assume "m \<noteq> 0"
   491   then have f_not_0: "f \<noteq> float_of 0" by (simp add: f_def)
   492   from mantissa_exponent[of f] have "m * 2 powr e = mantissa f * 2 powr exponent f"
   493     by (auto simp add: f_def)
   494   then show ?M ?E
   495     using mantissa_not_dvd[OF f_not_0] dvd
   496     by (auto simp: mult_powr_eq_mult_powr_iff)
   497 qed
   498 
   499 
   500 subsection \<open>Compute arithmetic operations\<close>
   501 
   502 lemma Float_mantissa_exponent: "Float (mantissa f) (exponent f) = f"
   503   unfolding real_of_float_eq mantissa_exponent[of f] by simp
   504 
   505 lemma Float_cases [cases type: float]:
   506   fixes f :: float
   507   obtains (Float) m e :: int where "f = Float m e"
   508   using Float_mantissa_exponent[symmetric]
   509   by (atomize_elim) auto
   510 
   511 lemma denormalize_shift:
   512   assumes f_def: "f \<equiv> Float m e"
   513     and not_0: "f \<noteq> float_of 0"
   514   obtains i where "m = mantissa f * 2 ^ i" "e = exponent f - i"
   515 proof
   516   from mantissa_exponent[of f] f_def
   517   have "m * 2 powr e = mantissa f * 2 powr exponent f"
   518     by simp
   519   then have eq: "m = mantissa f * 2 powr (exponent f - e)"
   520     by (simp add: powr_diff field_simps)
   521   moreover
   522   have "e \<le> exponent f"
   523   proof (rule ccontr)
   524     assume "\<not> e \<le> exponent f"
   525     then have pos: "exponent f < e" by simp
   526     then have "2 powr (exponent f - e) = 2 powr - real_of_int (e - exponent f)"
   527       by simp
   528     also have "\<dots> = 1 / 2^nat (e - exponent f)"
   529       using pos by (simp add: powr_realpow[symmetric] powr_diff)
   530     finally have "m * 2^nat (e - exponent f) = real_of_int (mantissa f)"
   531       using eq by simp
   532     then have "mantissa f = m * 2^nat (e - exponent f)"
   533       by linarith
   534     with \<open>exponent f < e\<close> have "2 dvd mantissa f"
   535       apply (intro dvdI[where k="m * 2^(nat (e-exponent f)) div 2"])
   536       apply (cases "nat (e - exponent f)")
   537       apply auto
   538       done
   539     then show False using mantissa_not_dvd[OF not_0] by simp
   540   qed
   541   ultimately have "real_of_int m = mantissa f * 2^nat (exponent f - e)"
   542     by (simp add: powr_realpow[symmetric])
   543   with \<open>e \<le> exponent f\<close>
   544   show "m = mantissa f * 2 ^ nat (exponent f - e)"
   545     by linarith
   546   show "e = exponent f - nat (exponent f - e)"
   547     using \<open>e \<le> exponent f\<close> by auto
   548 qed
   549 
   550 context
   551 begin
   552 
   553 qualified lemma compute_float_zero[code_unfold, code]: "0 = Float 0 0"
   554   by transfer simp
   555 
   556 qualified lemma compute_float_one[code_unfold, code]: "1 = Float 1 0"
   557   by transfer simp
   558 
   559 lift_definition normfloat :: "float \<Rightarrow> float" is "\<lambda>x. x" .
   560 lemma normloat_id[simp]: "normfloat x = x" by transfer rule
   561 
   562 qualified lemma compute_normfloat[code]:
   563   "normfloat (Float m e) =
   564     (if m mod 2 = 0 \<and> m \<noteq> 0 then normfloat (Float (m div 2) (e + 1))
   565      else if m = 0 then 0 else Float m e)"
   566   by transfer (auto simp add: powr_add zmod_eq_0_iff)
   567 
   568 qualified lemma compute_float_numeral[code_abbrev]: "Float (numeral k) 0 = numeral k"
   569   by transfer simp
   570 
   571 qualified lemma compute_float_neg_numeral[code_abbrev]: "Float (- numeral k) 0 = - numeral k"
   572   by transfer simp
   573 
   574 qualified lemma compute_float_uminus[code]: "- Float m1 e1 = Float (- m1) e1"
   575   by transfer simp
   576 
   577 qualified lemma compute_float_times[code]: "Float m1 e1 * Float m2 e2 = Float (m1 * m2) (e1 + e2)"
   578   by transfer (simp add: field_simps powr_add)
   579 
   580 qualified lemma compute_float_plus[code]:
   581   "Float m1 e1 + Float m2 e2 =
   582     (if m1 = 0 then Float m2 e2
   583      else if m2 = 0 then Float m1 e1
   584      else if e1 \<le> e2 then Float (m1 + m2 * 2^nat (e2 - e1)) e1
   585      else Float (m2 + m1 * 2^nat (e1 - e2)) e2)"
   586   by transfer (simp add: field_simps powr_realpow[symmetric] powr_diff)
   587 
   588 qualified lemma compute_float_minus[code]: "f - g = f + (-g)" for f g :: float
   589   by simp
   590 
   591 qualified lemma compute_float_sgn[code]:
   592   "sgn (Float m1 e1) = (if 0 < m1 then 1 else if m1 < 0 then -1 else 0)"
   593   by transfer (simp add: sgn_mult)
   594 
   595 lift_definition is_float_pos :: "float \<Rightarrow> bool" is "op < 0 :: real \<Rightarrow> bool" .
   596 
   597 qualified lemma compute_is_float_pos[code]: "is_float_pos (Float m e) \<longleftrightarrow> 0 < m"
   598   by transfer (auto simp add: zero_less_mult_iff not_le[symmetric, of _ 0])
   599 
   600 lift_definition is_float_nonneg :: "float \<Rightarrow> bool" is "op \<le> 0 :: real \<Rightarrow> bool" .
   601 
   602 qualified lemma compute_is_float_nonneg[code]: "is_float_nonneg (Float m e) \<longleftrightarrow> 0 \<le> m"
   603   by transfer (auto simp add: zero_le_mult_iff not_less[symmetric, of _ 0])
   604 
   605 lift_definition is_float_zero :: "float \<Rightarrow> bool"  is "op = 0 :: real \<Rightarrow> bool" .
   606 
   607 qualified lemma compute_is_float_zero[code]: "is_float_zero (Float m e) \<longleftrightarrow> 0 = m"
   608   by transfer (auto simp add: is_float_zero_def)
   609 
   610 qualified lemma compute_float_abs[code]: "\<bar>Float m e\<bar> = Float \<bar>m\<bar> e"
   611   by transfer (simp add: abs_mult)
   612 
   613 qualified lemma compute_float_eq[code]: "equal_class.equal f g = is_float_zero (f - g)"
   614   by transfer simp
   615 
   616 end
   617 
   618 
   619 subsection \<open>Lemmas for types @{typ real}, @{typ nat}, @{typ int}\<close>
   620 
   621 lemmas real_of_ints =
   622   of_int_add
   623   of_int_minus
   624   of_int_diff
   625   of_int_mult
   626   of_int_power
   627   of_int_numeral of_int_neg_numeral
   628 
   629 lemmas int_of_reals = real_of_ints[symmetric]
   630 
   631 
   632 subsection \<open>Rounding Real Numbers\<close>
   633 
   634 definition round_down :: "int \<Rightarrow> real \<Rightarrow> real"
   635   where "round_down prec x = \<lfloor>x * 2 powr prec\<rfloor> * 2 powr -prec"
   636 
   637 definition round_up :: "int \<Rightarrow> real \<Rightarrow> real"
   638   where "round_up prec x = \<lceil>x * 2 powr prec\<rceil> * 2 powr -prec"
   639 
   640 lemma round_down_float[simp]: "round_down prec x \<in> float"
   641   unfolding round_down_def
   642   by (auto intro!: times_float simp: of_int_minus[symmetric] simp del: of_int_minus)
   643 
   644 lemma round_up_float[simp]: "round_up prec x \<in> float"
   645   unfolding round_up_def
   646   by (auto intro!: times_float simp: of_int_minus[symmetric] simp del: of_int_minus)
   647 
   648 lemma round_up: "x \<le> round_up prec x"
   649   by (simp add: powr_minus_divide le_divide_eq round_up_def ceiling_correct)
   650 
   651 lemma round_down: "round_down prec x \<le> x"
   652   by (simp add: powr_minus_divide divide_le_eq round_down_def)
   653 
   654 lemma round_up_0[simp]: "round_up p 0 = 0"
   655   unfolding round_up_def by simp
   656 
   657 lemma round_down_0[simp]: "round_down p 0 = 0"
   658   unfolding round_down_def by simp
   659 
   660 lemma round_up_diff_round_down: "round_up prec x - round_down prec x \<le> 2 powr -prec"
   661 proof -
   662   have "round_up prec x - round_down prec x = (\<lceil>x * 2 powr prec\<rceil> - \<lfloor>x * 2 powr prec\<rfloor>) * 2 powr -prec"
   663     by (simp add: round_up_def round_down_def field_simps)
   664   also have "\<dots> \<le> 1 * 2 powr -prec"
   665     by (rule mult_mono)
   666       (auto simp del: of_int_diff simp: of_int_diff[symmetric] ceiling_diff_floor_le_1)
   667   finally show ?thesis by simp
   668 qed
   669 
   670 lemma round_down_shift: "round_down p (x * 2 powr k) = 2 powr k * round_down (p + k) x"
   671   unfolding round_down_def
   672   by (simp add: powr_add powr_mult field_simps powr_diff)
   673     (simp add: powr_add[symmetric])
   674 
   675 lemma round_up_shift: "round_up p (x * 2 powr k) = 2 powr k * round_up (p + k) x"
   676   unfolding round_up_def
   677   by (simp add: powr_add powr_mult field_simps powr_diff)
   678     (simp add: powr_add[symmetric])
   679 
   680 lemma round_up_uminus_eq: "round_up p (-x) = - round_down p x"
   681   and round_down_uminus_eq: "round_down p (-x) = - round_up p x"
   682   by (auto simp: round_up_def round_down_def ceiling_def)
   683 
   684 lemma round_up_mono: "x \<le> y \<Longrightarrow> round_up p x \<le> round_up p y"
   685   by (auto intro!: ceiling_mono simp: round_up_def)
   686 
   687 lemma round_up_le1:
   688   assumes "x \<le> 1" "prec \<ge> 0"
   689   shows "round_up prec x \<le> 1"
   690 proof -
   691   have "real_of_int \<lceil>x * 2 powr prec\<rceil> \<le> real_of_int \<lceil>2 powr real_of_int prec\<rceil>"
   692     using assms by (auto intro!: ceiling_mono)
   693   also have "\<dots> = 2 powr prec" using assms by (auto simp: powr_int intro!: exI[where x="2^nat prec"])
   694   finally show ?thesis
   695     by (simp add: round_up_def) (simp add: powr_minus inverse_eq_divide)
   696 qed
   697 
   698 lemma round_up_less1:
   699   assumes "x < 1 / 2" "p > 0"
   700   shows "round_up p x < 1"
   701 proof -
   702   have "x * 2 powr p < 1 / 2 * 2 powr p"
   703     using assms by simp
   704   also have "\<dots> \<le> 2 powr p - 1" using \<open>p > 0\<close>
   705     by (auto simp: powr_diff powr_int field_simps self_le_power)
   706   finally show ?thesis using \<open>p > 0\<close>
   707     by (simp add: round_up_def field_simps powr_minus powr_int ceiling_less_iff)
   708 qed
   709 
   710 lemma round_down_ge1:
   711   assumes x: "x \<ge> 1"
   712   assumes prec: "p \<ge> - log 2 x"
   713   shows "1 \<le> round_down p x"
   714 proof cases
   715   assume nonneg: "0 \<le> p"
   716   have "2 powr p = real_of_int \<lfloor>2 powr real_of_int p\<rfloor>"
   717     using nonneg by (auto simp: powr_int)
   718   also have "\<dots> \<le> real_of_int \<lfloor>x * 2 powr p\<rfloor>"
   719     using assms by (auto intro!: floor_mono)
   720   finally show ?thesis
   721     by (simp add: round_down_def) (simp add: powr_minus inverse_eq_divide)
   722 next
   723   assume neg: "\<not> 0 \<le> p"
   724   have "x = 2 powr (log 2 x)"
   725     using x by simp
   726   also have "2 powr (log 2 x) \<ge> 2 powr - p"
   727     using prec by auto
   728   finally have x_le: "x \<ge> 2 powr -p" .
   729 
   730   from neg have "2 powr real_of_int p \<le> 2 powr 0"
   731     by (intro powr_mono) auto
   732   also have "\<dots> \<le> \<lfloor>2 powr 0::real\<rfloor>" by simp
   733   also have "\<dots> \<le> \<lfloor>x * 2 powr (real_of_int p)\<rfloor>"
   734     unfolding of_int_le_iff
   735     using x x_le by (intro floor_mono) (simp add: powr_minus_divide field_simps)
   736   finally show ?thesis
   737     using prec x
   738     by (simp add: round_down_def powr_minus_divide pos_le_divide_eq)
   739 qed
   740 
   741 lemma round_up_le0: "x \<le> 0 \<Longrightarrow> round_up p x \<le> 0"
   742   unfolding round_up_def
   743   by (auto simp: field_simps mult_le_0_iff zero_le_mult_iff)
   744 
   745 
   746 subsection \<open>Rounding Floats\<close>
   747 
   748 definition div_twopow :: "int \<Rightarrow> nat \<Rightarrow> int"
   749   where [simp]: "div_twopow x n = x div (2 ^ n)"
   750 
   751 definition mod_twopow :: "int \<Rightarrow> nat \<Rightarrow> int"
   752   where [simp]: "mod_twopow x n = x mod (2 ^ n)"
   753 
   754 lemma compute_div_twopow[code]:
   755   "div_twopow x n = (if x = 0 \<or> x = -1 \<or> n = 0 then x else div_twopow (x div 2) (n - 1))"
   756   by (cases n) (auto simp: zdiv_zmult2_eq div_eq_minus1)
   757 
   758 lemma compute_mod_twopow[code]:
   759   "mod_twopow x n = (if n = 0 then 0 else x mod 2 + 2 * mod_twopow (x div 2) (n - 1))"
   760   by (cases n) (auto simp: zmod_zmult2_eq)
   761 
   762 lift_definition float_up :: "int \<Rightarrow> float \<Rightarrow> float" is round_up by simp
   763 declare float_up.rep_eq[simp]
   764 
   765 lemma round_up_correct: "round_up e f - f \<in> {0..2 powr -e}"
   766   unfolding atLeastAtMost_iff
   767 proof
   768   have "round_up e f - f \<le> round_up e f - round_down e f"
   769     using round_down by simp
   770   also have "\<dots> \<le> 2 powr -e"
   771     using round_up_diff_round_down by simp
   772   finally show "round_up e f - f \<le> 2 powr - (real_of_int e)"
   773     by simp
   774 qed (simp add: algebra_simps round_up)
   775 
   776 lemma float_up_correct: "real_of_float (float_up e f) - real_of_float f \<in> {0..2 powr -e}"
   777   by transfer (rule round_up_correct)
   778 
   779 lift_definition float_down :: "int \<Rightarrow> float \<Rightarrow> float" is round_down by simp
   780 declare float_down.rep_eq[simp]
   781 
   782 lemma round_down_correct: "f - (round_down e f) \<in> {0..2 powr -e}"
   783   unfolding atLeastAtMost_iff
   784 proof
   785   have "f - round_down e f \<le> round_up e f - round_down e f"
   786     using round_up by simp
   787   also have "\<dots> \<le> 2 powr -e"
   788     using round_up_diff_round_down by simp
   789   finally show "f - round_down e f \<le> 2 powr - (real_of_int e)"
   790     by simp
   791 qed (simp add: algebra_simps round_down)
   792 
   793 lemma float_down_correct: "real_of_float f - real_of_float (float_down e f) \<in> {0..2 powr -e}"
   794   by transfer (rule round_down_correct)
   795 
   796 context
   797 begin
   798 
   799 qualified lemma compute_float_down[code]:
   800   "float_down p (Float m e) =
   801     (if p + e < 0 then Float (div_twopow m (nat (-(p + e)))) (-p) else Float m e)"
   802 proof (cases "p + e < 0")
   803   case True
   804   then have "real_of_int ((2::int) ^ nat (-(p + e))) = 2 powr (-(p + e))"
   805     using powr_realpow[of 2 "nat (-(p + e))"] by simp
   806   also have "\<dots> = 1 / 2 powr p / 2 powr e"
   807     unfolding powr_minus_divide of_int_minus by (simp add: powr_add)
   808   finally show ?thesis
   809     using \<open>p + e < 0\<close>
   810     apply transfer
   811     apply (simp add: ac_simps round_down_def floor_divide_of_int_eq[symmetric])
   812     proof - (*FIXME*)
   813       fix pa :: int and ea :: int and ma :: int
   814       assume a1: "2 ^ nat (- pa - ea) = 1 / (2 powr real_of_int pa * 2 powr real_of_int ea)"
   815       assume "pa + ea < 0"
   816       have "\<lfloor>real_of_int ma / real_of_int (int 2 ^ nat (- (pa + ea)))\<rfloor> =
   817           \<lfloor>real_of_float (Float ma (pa + ea))\<rfloor>"
   818         using a1 by (simp add: powr_add)
   819       then show "\<lfloor>real_of_int ma * (2 powr real_of_int pa * 2 powr real_of_int ea)\<rfloor> =
   820           ma div 2 ^ nat (- pa - ea)"
   821         by (metis Float.rep_eq add_uminus_conv_diff floor_divide_of_int_eq
   822             minus_add_distrib of_int_simps(3) of_nat_numeral powr_add)
   823     qed
   824 next
   825   case False
   826   then have r: "real_of_int e + real_of_int p = real (nat (e + p))"
   827     by simp
   828   have r: "\<lfloor>(m * 2 powr e) * 2 powr real_of_int p\<rfloor> = (m * 2 powr e) * 2 powr real_of_int p"
   829     by (auto intro: exI[where x="m*2^nat (e+p)"]
   830         simp add: ac_simps powr_add[symmetric] r powr_realpow)
   831   with \<open>\<not> p + e < 0\<close> show ?thesis
   832     by transfer (auto simp add: round_down_def field_simps powr_add powr_minus)
   833 qed
   834 
   835 lemma abs_round_down_le: "\<bar>f - (round_down e f)\<bar> \<le> 2 powr -e"
   836   using round_down_correct[of f e] by simp
   837 
   838 lemma abs_round_up_le: "\<bar>f - (round_up e f)\<bar> \<le> 2 powr -e"
   839   using round_up_correct[of e f] by simp
   840 
   841 lemma round_down_nonneg: "0 \<le> s \<Longrightarrow> 0 \<le> round_down p s"
   842   by (auto simp: round_down_def)
   843 
   844 lemma ceil_divide_floor_conv:
   845   assumes "b \<noteq> 0"
   846   shows "\<lceil>real_of_int a / real_of_int b\<rceil> =
   847     (if b dvd a then a div b else \<lfloor>real_of_int a / real_of_int b\<rfloor> + 1)"
   848 proof (cases "b dvd a")
   849   case True
   850   then show ?thesis
   851     by (simp add: ceiling_def of_int_minus[symmetric] divide_minus_left[symmetric]
   852       floor_divide_of_int_eq dvd_neg_div del: divide_minus_left of_int_minus)
   853 next
   854   case False
   855   then have "a mod b \<noteq> 0"
   856     by auto
   857   then have ne: "real_of_int (a mod b) / real_of_int b \<noteq> 0"
   858     using \<open>b \<noteq> 0\<close> by auto
   859   have "\<lceil>real_of_int a / real_of_int b\<rceil> = \<lfloor>real_of_int a / real_of_int b\<rfloor> + 1"
   860     apply (rule ceiling_eq)
   861     apply (auto simp: floor_divide_of_int_eq[symmetric])
   862   proof -
   863     have "real_of_int \<lfloor>real_of_int a / real_of_int b\<rfloor> \<le> real_of_int a / real_of_int b"
   864       by simp
   865     moreover have "real_of_int \<lfloor>real_of_int a / real_of_int b\<rfloor> \<noteq> real_of_int a / real_of_int b"
   866       apply (subst (2) real_of_int_div_aux)
   867       unfolding floor_divide_of_int_eq
   868       using ne \<open>b \<noteq> 0\<close> apply auto
   869       done
   870     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
   871   qed
   872   then show ?thesis
   873     using \<open>\<not> b dvd a\<close> by simp
   874 qed
   875 
   876 qualified lemma compute_float_up[code]: "float_up p x = - float_down p (-x)"
   877   by transfer (simp add: round_down_uminus_eq)
   878 
   879 end
   880 
   881 
   882 lemma bitlen_Float:
   883   fixes m e
   884   defines "f \<equiv> Float m e"
   885   shows "bitlen \<bar>mantissa f\<bar> + exponent f = (if m = 0 then 0 else bitlen \<bar>m\<bar> + e)"
   886 proof (cases "m = 0")
   887   case True
   888   then show ?thesis by (simp add: f_def bitlen_alt_def Float_def)
   889 next
   890   case False
   891   then have "f \<noteq> float_of 0"
   892     unfolding real_of_float_eq by (simp add: f_def)
   893   then have "mantissa f \<noteq> 0"
   894     by (simp add: mantissa_noteq_0)
   895   moreover
   896   obtain i where "m = mantissa f * 2 ^ i" "e = exponent f - int i"
   897     by (rule f_def[THEN denormalize_shift, OF \<open>f \<noteq> float_of 0\<close>])
   898   ultimately show ?thesis by (simp add: abs_mult)
   899 qed
   900 
   901 lemma float_gt1_scale:
   902   assumes "1 \<le> Float m e"
   903   shows "0 \<le> e + (bitlen m - 1)"
   904 proof -
   905   have "0 < Float m e" using assms by auto
   906   then have "0 < m" using powr_gt_zero[of 2 e]
   907     apply (auto simp: zero_less_mult_iff)
   908     using not_le powr_ge_pzero
   909     apply blast
   910     done
   911   then have "m \<noteq> 0" by auto
   912   show ?thesis
   913   proof (cases "0 \<le> e")
   914     case True
   915     then show ?thesis
   916       using \<open>0 < m\<close> by (simp add: bitlen_alt_def)
   917   next
   918     case False
   919     have "(1::int) < 2" by simp
   920     let ?S = "2^(nat (-e))"
   921     have "inverse (2 ^ nat (- e)) = 2 powr e"
   922       using assms False powr_realpow[of 2 "nat (-e)"]
   923       by (auto simp: powr_minus field_simps)
   924     then have "1 \<le> real_of_int m * inverse ?S"
   925       using assms False powr_realpow[of 2 "nat (-e)"]
   926       by (auto simp: powr_minus)
   927     then have "1 * ?S \<le> real_of_int m * inverse ?S * ?S"
   928       by (rule mult_right_mono) auto
   929     then have "?S \<le> real_of_int m"
   930       unfolding mult.assoc by auto
   931     then have "?S \<le> m"
   932       unfolding of_int_le_iff[symmetric] by auto
   933     from this bitlen_bounds[OF \<open>0 < m\<close>, THEN conjunct2]
   934     have "nat (-e) < (nat (bitlen m))"
   935       unfolding power_strict_increasing_iff[OF \<open>1 < 2\<close>, symmetric]
   936       by (rule order_le_less_trans)
   937     then have "-e < bitlen m"
   938       using False by auto
   939     then show ?thesis
   940       by auto
   941   qed
   942 qed
   943 
   944 
   945 subsection \<open>Truncating Real Numbers\<close>
   946 
   947 definition truncate_down::"nat \<Rightarrow> real \<Rightarrow> real"
   948   where "truncate_down prec x = round_down (prec - \<lfloor>log 2 \<bar>x\<bar>\<rfloor>) x"
   949 
   950 lemma truncate_down: "truncate_down prec x \<le> x"
   951   using round_down by (simp add: truncate_down_def)
   952 
   953 lemma truncate_down_le: "x \<le> y \<Longrightarrow> truncate_down prec x \<le> y"
   954   by (rule order_trans[OF truncate_down])
   955 
   956 lemma truncate_down_zero[simp]: "truncate_down prec 0 = 0"
   957   by (simp add: truncate_down_def)
   958 
   959 lemma truncate_down_float[simp]: "truncate_down p x \<in> float"
   960   by (auto simp: truncate_down_def)
   961 
   962 definition truncate_up::"nat \<Rightarrow> real \<Rightarrow> real"
   963   where "truncate_up prec x = round_up (prec - \<lfloor>log 2 \<bar>x\<bar>\<rfloor>) x"
   964 
   965 lemma truncate_up: "x \<le> truncate_up prec x"
   966   using round_up by (simp add: truncate_up_def)
   967 
   968 lemma truncate_up_le: "x \<le> y \<Longrightarrow> x \<le> truncate_up prec y"
   969   by (rule order_trans[OF _ truncate_up])
   970 
   971 lemma truncate_up_zero[simp]: "truncate_up prec 0 = 0"
   972   by (simp add: truncate_up_def)
   973 
   974 lemma truncate_up_uminus_eq: "truncate_up prec (-x) = - truncate_down prec x"
   975   and truncate_down_uminus_eq: "truncate_down prec (-x) = - truncate_up prec x"
   976   by (auto simp: truncate_up_def round_up_def truncate_down_def round_down_def ceiling_def)
   977 
   978 lemma truncate_up_float[simp]: "truncate_up p x \<in> float"
   979   by (auto simp: truncate_up_def)
   980 
   981 lemma mult_powr_eq: "0 < b \<Longrightarrow> b \<noteq> 1 \<Longrightarrow> 0 < x \<Longrightarrow> x * b powr y = b powr (y + log b x)"
   982   by (simp_all add: powr_add)
   983 
   984 lemma truncate_down_pos:
   985   assumes "x > 0"
   986   shows "truncate_down p x > 0"
   987 proof -
   988   have "0 \<le> log 2 x - real_of_int \<lfloor>log 2 x\<rfloor>"
   989     by (simp add: algebra_simps)
   990   with assms
   991   show ?thesis
   992     apply (auto simp: truncate_down_def round_down_def mult_powr_eq
   993       intro!: ge_one_powr_ge_zero mult_pos_pos)
   994     by linarith
   995 qed
   996 
   997 lemma truncate_down_nonneg: "0 \<le> y \<Longrightarrow> 0 \<le> truncate_down prec y"
   998   by (auto simp: truncate_down_def round_down_def)
   999 
  1000 lemma truncate_down_ge1: "1 \<le> x \<Longrightarrow> 1 \<le> truncate_down p x"
  1001   apply (auto simp: truncate_down_def algebra_simps intro!: round_down_ge1)
  1002   apply linarith
  1003   done
  1004 
  1005 lemma truncate_up_nonpos: "x \<le> 0 \<Longrightarrow> truncate_up prec x \<le> 0"
  1006   by (auto simp: truncate_up_def round_up_def intro!: mult_nonpos_nonneg)
  1007 
  1008 lemma truncate_up_le1:
  1009   assumes "x \<le> 1"
  1010   shows "truncate_up p x \<le> 1"
  1011 proof -
  1012   consider "x \<le> 0" | "x > 0"
  1013     by arith
  1014   then show ?thesis
  1015   proof cases
  1016     case 1
  1017     with truncate_up_nonpos[OF this, of p] show ?thesis
  1018       by simp
  1019   next
  1020     case 2
  1021     then have le: "\<lfloor>log 2 \<bar>x\<bar>\<rfloor> \<le> 0"
  1022       using assms by (auto simp: log_less_iff)
  1023     from assms have "0 \<le> int p" by simp
  1024     from add_mono[OF this le]
  1025     show ?thesis
  1026       using assms by (simp add: truncate_up_def round_up_le1 add_mono)
  1027   qed
  1028 qed
  1029 
  1030 lemma truncate_down_shift_int:
  1031   "truncate_down p (x * 2 powr real_of_int k) = truncate_down p x * 2 powr k"
  1032   by (cases "x = 0")
  1033     (simp_all add: algebra_simps abs_mult log_mult truncate_down_def
  1034       round_down_shift[of _ _ k, simplified])
  1035 
  1036 lemma truncate_down_shift_nat: "truncate_down p (x * 2 powr real k) = truncate_down p x * 2 powr k"
  1037   by (metis of_int_of_nat_eq truncate_down_shift_int)
  1038 
  1039 lemma truncate_up_shift_int: "truncate_up p (x * 2 powr real_of_int k) = truncate_up p x * 2 powr k"
  1040   by (cases "x = 0")
  1041     (simp_all add: algebra_simps abs_mult log_mult truncate_up_def
  1042       round_up_shift[of _ _ k, simplified])
  1043 
  1044 lemma truncate_up_shift_nat: "truncate_up p (x * 2 powr real k) = truncate_up p x * 2 powr k"
  1045   by (metis of_int_of_nat_eq truncate_up_shift_int)
  1046 
  1047 
  1048 subsection \<open>Truncating Floats\<close>
  1049 
  1050 lift_definition float_round_up :: "nat \<Rightarrow> float \<Rightarrow> float" is truncate_up
  1051   by (simp add: truncate_up_def)
  1052 
  1053 lemma float_round_up: "real_of_float x \<le> real_of_float (float_round_up prec x)"
  1054   using truncate_up by transfer simp
  1055 
  1056 lemma float_round_up_zero[simp]: "float_round_up prec 0 = 0"
  1057   by transfer simp
  1058 
  1059 lift_definition float_round_down :: "nat \<Rightarrow> float \<Rightarrow> float" is truncate_down
  1060   by (simp add: truncate_down_def)
  1061 
  1062 lemma float_round_down: "real_of_float (float_round_down prec x) \<le> real_of_float x"
  1063   using truncate_down by transfer simp
  1064 
  1065 lemma float_round_down_zero[simp]: "float_round_down prec 0 = 0"
  1066   by transfer simp
  1067 
  1068 lemmas float_round_up_le = order_trans[OF _ float_round_up]
  1069   and float_round_down_le = order_trans[OF float_round_down]
  1070 
  1071 lemma minus_float_round_up_eq: "- float_round_up prec x = float_round_down prec (- x)"
  1072   and minus_float_round_down_eq: "- float_round_down prec x = float_round_up prec (- x)"
  1073   by (transfer; simp add: truncate_down_uminus_eq truncate_up_uminus_eq)+
  1074 
  1075 context
  1076 begin
  1077 
  1078 qualified lemma compute_float_round_down[code]:
  1079   "float_round_down prec (Float m e) =
  1080     (let d = bitlen \<bar>m\<bar> - int prec - 1 in
  1081       if 0 < d then Float (div_twopow m (nat d)) (e + d)
  1082       else Float m e)"
  1083   using Float.compute_float_down[of "Suc prec - bitlen \<bar>m\<bar> - e" m e, symmetric]
  1084   by transfer
  1085     (simp add: field_simps abs_mult log_mult bitlen_alt_def truncate_down_def
  1086       cong del: if_weak_cong)
  1087 
  1088 qualified lemma compute_float_round_up[code]:
  1089   "float_round_up prec x = - float_round_down prec (-x)"
  1090   by transfer (simp add: truncate_down_uminus_eq)
  1091 
  1092 end
  1093 
  1094 
  1095 subsection \<open>Approximation of positive rationals\<close>
  1096 
  1097 lemma div_mult_twopow_eq: "a div ((2::nat) ^ n) div b = a div (b * 2 ^ n)" for a b :: nat
  1098   by (cases "b = 0") (simp_all add: div_mult2_eq[symmetric] ac_simps)
  1099 
  1100 lemma real_div_nat_eq_floor_of_divide: "a div b = real_of_int \<lfloor>a / b\<rfloor>" for a b :: nat
  1101   by (simp add: floor_divide_of_nat_eq [of a b])
  1102 
  1103 definition "rat_precision prec x y =
  1104   (let d = bitlen x - bitlen y
  1105    in int prec - d + (if Float (abs x) 0 < Float (abs y) d then 1 else 0))"
  1106 
  1107 lemma floor_log_divide_eq:
  1108   assumes "i > 0" "j > 0" "p > 1"
  1109   shows "\<lfloor>log p (i / j)\<rfloor> = floor (log p i) - floor (log p j) -
  1110     (if i \<ge> j * p powr (floor (log p i) - floor (log p j)) then 0 else 1)"
  1111 proof -
  1112   let ?l = "log p"
  1113   let ?fl = "\<lambda>x. floor (?l x)"
  1114   have "\<lfloor>?l (i / j)\<rfloor> = \<lfloor>?l i - ?l j\<rfloor>" using assms
  1115     by (auto simp: log_divide)
  1116   also have "\<dots> = floor (real_of_int (?fl i - ?fl j) + (?l i - ?fl i - (?l j - ?fl j)))"
  1117     (is "_ = floor (_ + ?r)")
  1118     by (simp add: algebra_simps)
  1119   also note floor_add2
  1120   also note \<open>p > 1\<close>
  1121   note powr = powr_le_cancel_iff[symmetric, OF \<open>1 < p\<close>, THEN iffD2]
  1122   note powr_strict = powr_less_cancel_iff[symmetric, OF \<open>1 < p\<close>, THEN iffD2]
  1123   have "floor ?r = (if i \<ge> j * p powr (?fl i - ?fl j) then 0 else -1)" (is "_ = ?if")
  1124     using assms
  1125     by (linarith |
  1126       auto
  1127         intro!: floor_eq2
  1128         intro: powr_strict powr
  1129         simp: powr_diff powr_add divide_simps algebra_simps)+
  1130   finally
  1131   show ?thesis by simp
  1132 qed
  1133 
  1134 lemma truncate_down_rat_precision:
  1135     "truncate_down prec (real x / real y) = round_down (rat_precision prec x y) (real x / real y)"
  1136   and truncate_up_rat_precision:
  1137     "truncate_up prec (real x / real y) = round_up (rat_precision prec x y) (real x / real y)"
  1138   unfolding truncate_down_def truncate_up_def rat_precision_def
  1139   by (cases x; cases y) (auto simp: floor_log_divide_eq algebra_simps bitlen_alt_def)
  1140 
  1141 lift_definition lapprox_posrat :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float"
  1142   is "\<lambda>prec (x::nat) (y::nat). truncate_down prec (x / y)"
  1143   by simp
  1144 
  1145 context
  1146 begin
  1147 
  1148 qualified lemma compute_lapprox_posrat[code]:
  1149   "lapprox_posrat prec x y =
  1150    (let
  1151       l = rat_precision prec x y;
  1152       d = if 0 \<le> l then x * 2^nat l div y else x div 2^nat (- l) div y
  1153     in normfloat (Float d (- l)))"
  1154     unfolding div_mult_twopow_eq
  1155     by transfer
  1156       (simp add: round_down_def powr_int real_div_nat_eq_floor_of_divide field_simps Let_def
  1157         truncate_down_rat_precision del: two_powr_minus_int_float)
  1158 
  1159 end
  1160 
  1161 lift_definition rapprox_posrat :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float"
  1162   is "\<lambda>prec (x::nat) (y::nat). truncate_up prec (x / y)"
  1163   by simp
  1164 
  1165 context
  1166 begin
  1167 
  1168 qualified lemma compute_rapprox_posrat[code]:
  1169   fixes prec x y
  1170   defines "l \<equiv> rat_precision prec x y"
  1171   shows "rapprox_posrat prec x y =
  1172    (let
  1173       l = l;
  1174       (r, s) = if 0 \<le> l then (x * 2^nat l, y) else (x, y * 2^nat(-l));
  1175       d = r div s;
  1176       m = r mod s
  1177     in normfloat (Float (d + (if m = 0 \<or> y = 0 then 0 else 1)) (- l)))"
  1178 proof (cases "y = 0")
  1179   assume "y = 0"
  1180   then show ?thesis by transfer simp
  1181 next
  1182   assume "y \<noteq> 0"
  1183   show ?thesis
  1184   proof (cases "0 \<le> l")
  1185     case True
  1186     define x' where "x' = x * 2 ^ nat l"
  1187     have "int x * 2 ^ nat l = x'"
  1188       by (simp add: x'_def)
  1189     moreover have "real x * 2 powr l = real x'"
  1190       by (simp add: powr_realpow[symmetric] \<open>0 \<le> l\<close> x'_def)
  1191     ultimately show ?thesis
  1192       using ceil_divide_floor_conv[of y x'] powr_realpow[of 2 "nat l"] \<open>0 \<le> l\<close> \<open>y \<noteq> 0\<close>
  1193         l_def[symmetric, THEN meta_eq_to_obj_eq]
  1194       apply transfer
  1195       apply (auto simp add: round_up_def truncate_up_rat_precision)
  1196       apply (metis floor_divide_of_int_eq of_int_of_nat_eq)
  1197       done
  1198    next
  1199     case False
  1200     define y' where "y' = y * 2 ^ nat (- l)"
  1201     from \<open>y \<noteq> 0\<close> have "y' \<noteq> 0" by (simp add: y'_def)
  1202     have "int y * 2 ^ nat (- l) = y'"
  1203       by (simp add: y'_def)
  1204     moreover have "real x * real_of_int (2::int) powr real_of_int l / real y = x / real y'"
  1205       using \<open>\<not> 0 \<le> l\<close> by (simp add: powr_realpow[symmetric] powr_minus y'_def field_simps)
  1206     ultimately show ?thesis
  1207       using ceil_divide_floor_conv[of y' x] \<open>\<not> 0 \<le> l\<close> \<open>y' \<noteq> 0\<close> \<open>y \<noteq> 0\<close>
  1208         l_def[symmetric, THEN meta_eq_to_obj_eq]
  1209       apply transfer
  1210       apply (auto simp add: round_up_def ceil_divide_floor_conv truncate_up_rat_precision)
  1211       apply (metis floor_divide_of_int_eq of_int_of_nat_eq)
  1212       done
  1213   qed
  1214 qed
  1215 
  1216 end
  1217 
  1218 lemma rat_precision_pos:
  1219   assumes "0 \<le> x"
  1220     and "0 < y"
  1221     and "2 * x < y"
  1222   shows "rat_precision n (int x) (int y) > 0"
  1223 proof -
  1224   have "0 < x \<Longrightarrow> log 2 x + 1 = log 2 (2 * x)"
  1225     by (simp add: log_mult)
  1226   then have "bitlen (int x) < bitlen (int y)"
  1227     using assms
  1228     by (simp add: bitlen_alt_def)
  1229       (auto intro!: floor_mono simp add: one_add_floor)
  1230   then show ?thesis
  1231     using assms
  1232     by (auto intro!: pos_add_strict simp add: field_simps rat_precision_def)
  1233 qed
  1234 
  1235 lemma rapprox_posrat_less1:
  1236   "0 \<le> x \<Longrightarrow> 0 < y \<Longrightarrow> 2 * x < y \<Longrightarrow> real_of_float (rapprox_posrat n x y) < 1"
  1237   by transfer (simp add: rat_precision_pos round_up_less1 truncate_up_rat_precision)
  1238 
  1239 lift_definition lapprox_rat :: "nat \<Rightarrow> int \<Rightarrow> int \<Rightarrow> float" is
  1240   "\<lambda>prec (x::int) (y::int). truncate_down prec (x / y)"
  1241   by simp
  1242 
  1243 context
  1244 begin
  1245 
  1246 qualified lemma compute_lapprox_rat[code]:
  1247   "lapprox_rat prec x y =
  1248    (if y = 0 then 0
  1249     else if 0 \<le> x then
  1250      (if 0 < y then lapprox_posrat prec (nat x) (nat y)
  1251       else - (rapprox_posrat prec (nat x) (nat (-y))))
  1252       else
  1253        (if 0 < y
  1254         then - (rapprox_posrat prec (nat (-x)) (nat y))
  1255         else lapprox_posrat prec (nat (-x)) (nat (-y))))"
  1256   by transfer (simp add: truncate_up_uminus_eq)
  1257 
  1258 lift_definition rapprox_rat :: "nat \<Rightarrow> int \<Rightarrow> int \<Rightarrow> float" is
  1259   "\<lambda>prec (x::int) (y::int). truncate_up prec (x / y)"
  1260   by simp
  1261 
  1262 lemma "rapprox_rat = rapprox_posrat"
  1263   by transfer auto
  1264 
  1265 lemma "lapprox_rat = lapprox_posrat"
  1266   by transfer auto
  1267 
  1268 qualified lemma compute_rapprox_rat[code]:
  1269   "rapprox_rat prec x y = - lapprox_rat prec (-x) y"
  1270   by transfer (simp add: truncate_down_uminus_eq)
  1271 
  1272 qualified lemma compute_truncate_down[code]:
  1273   "truncate_down p (Ratreal r) = (let (a, b) = quotient_of r in lapprox_rat p a b)"
  1274   by transfer (auto split: prod.split simp: of_rat_divide dest!: quotient_of_div)
  1275 
  1276 qualified lemma compute_truncate_up[code]:
  1277   "truncate_up p (Ratreal r) = (let (a, b) = quotient_of r in rapprox_rat p a b)"
  1278   by transfer (auto split: prod.split simp:  of_rat_divide dest!: quotient_of_div)
  1279 
  1280 end
  1281 
  1282 
  1283 subsection \<open>Division\<close>
  1284 
  1285 definition "real_divl prec a b = truncate_down prec (a / b)"
  1286 
  1287 definition "real_divr prec a b = truncate_up prec (a / b)"
  1288 
  1289 lift_definition float_divl :: "nat \<Rightarrow> float \<Rightarrow> float \<Rightarrow> float" is real_divl
  1290   by (simp add: real_divl_def)
  1291 
  1292 context
  1293 begin
  1294 
  1295 qualified lemma compute_float_divl[code]:
  1296   "float_divl prec (Float m1 s1) (Float m2 s2) = lapprox_rat prec m1 m2 * Float 1 (s1 - s2)"
  1297   apply transfer
  1298   unfolding real_divl_def of_int_1 mult_1 truncate_down_shift_int[symmetric]
  1299   apply (simp add: powr_diff powr_minus)
  1300   done
  1301 
  1302 lift_definition float_divr :: "nat \<Rightarrow> float \<Rightarrow> float \<Rightarrow> float" is real_divr
  1303   by (simp add: real_divr_def)
  1304 
  1305 qualified lemma compute_float_divr[code]:
  1306   "float_divr prec x y = - float_divl prec (-x) y"
  1307   by transfer (simp add: real_divr_def real_divl_def truncate_down_uminus_eq)
  1308 
  1309 end
  1310 
  1311 
  1312 subsection \<open>Approximate Power\<close>
  1313 
  1314 lemma div2_less_self[termination_simp]: "odd n \<Longrightarrow> n div 2 < n" for n :: nat
  1315   by (simp add: odd_pos)
  1316 
  1317 fun power_down :: "nat \<Rightarrow> real \<Rightarrow> nat \<Rightarrow> real"
  1318 where
  1319   "power_down p x 0 = 1"
  1320 | "power_down p x (Suc n) =
  1321     (if odd n then truncate_down (Suc p) ((power_down p x (Suc n div 2))\<^sup>2)
  1322      else truncate_down (Suc p) (x * power_down p x n))"
  1323 
  1324 fun power_up :: "nat \<Rightarrow> real \<Rightarrow> nat \<Rightarrow> real"
  1325 where
  1326   "power_up p x 0 = 1"
  1327 | "power_up p x (Suc n) =
  1328     (if odd n then truncate_up p ((power_up p x (Suc n div 2))\<^sup>2)
  1329      else truncate_up p (x * power_up p x n))"
  1330 
  1331 lift_definition power_up_fl :: "nat \<Rightarrow> float \<Rightarrow> nat \<Rightarrow> float" is power_up
  1332   by (induct_tac rule: power_up.induct) simp_all
  1333 
  1334 lift_definition power_down_fl :: "nat \<Rightarrow> float \<Rightarrow> nat \<Rightarrow> float" is power_down
  1335   by (induct_tac rule: power_down.induct) simp_all
  1336 
  1337 lemma power_float_transfer[transfer_rule]:
  1338   "(rel_fun pcr_float (rel_fun op = pcr_float)) op ^ op ^"
  1339   unfolding power_def
  1340   by transfer_prover
  1341 
  1342 lemma compute_power_up_fl[code]:
  1343     "power_up_fl p x 0 = 1"
  1344     "power_up_fl p x (Suc n) =
  1345       (if odd n then float_round_up p ((power_up_fl p x (Suc n div 2))\<^sup>2)
  1346        else float_round_up p (x * power_up_fl p x n))"
  1347   and compute_power_down_fl[code]:
  1348     "power_down_fl p x 0 = 1"
  1349     "power_down_fl p x (Suc n) =
  1350       (if odd n then float_round_down (Suc p) ((power_down_fl p x (Suc n div 2))\<^sup>2)
  1351        else float_round_down (Suc p) (x * power_down_fl p x n))"
  1352   unfolding atomize_conj by transfer simp
  1353 
  1354 lemma power_down_pos: "0 < x \<Longrightarrow> 0 < power_down p x n"
  1355   by (induct p x n rule: power_down.induct)
  1356     (auto simp del: odd_Suc_div_two intro!: truncate_down_pos)
  1357 
  1358 lemma power_down_nonneg: "0 \<le> x \<Longrightarrow> 0 \<le> power_down p x n"
  1359   by (induct p x n rule: power_down.induct)
  1360     (auto simp del: odd_Suc_div_two intro!: truncate_down_nonneg mult_nonneg_nonneg)
  1361 
  1362 lemma power_down: "0 \<le> x \<Longrightarrow> power_down p x n \<le> x ^ n"
  1363 proof (induct p x n rule: power_down.induct)
  1364   case (2 p x n)
  1365   have ?case if "odd n"
  1366   proof -
  1367     from that 2 have "(power_down p x (Suc n div 2)) ^ 2 \<le> (x ^ (Suc n div 2)) ^ 2"
  1368       by (auto intro: power_mono power_down_nonneg simp del: odd_Suc_div_two)
  1369     also have "\<dots> = x ^ (Suc n div 2 * 2)"
  1370       by (simp add: power_mult[symmetric])
  1371     also have "Suc n div 2 * 2 = Suc n"
  1372       using \<open>odd n\<close> by presburger
  1373     finally show ?thesis
  1374       using that by (auto intro!: truncate_down_le simp del: odd_Suc_div_two)
  1375   qed
  1376   then show ?case
  1377     by (auto intro!: truncate_down_le mult_left_mono 2 mult_nonneg_nonneg power_down_nonneg)
  1378 qed simp
  1379 
  1380 lemma power_up: "0 \<le> x \<Longrightarrow> x ^ n \<le> power_up p x n"
  1381 proof (induct p x n rule: power_up.induct)
  1382   case (2 p x n)
  1383   have ?case if "odd n"
  1384   proof -
  1385     from that even_Suc have "Suc n = Suc n div 2 * 2"
  1386       by presburger
  1387     then have "x ^ Suc n \<le> (x ^ (Suc n div 2))\<^sup>2"
  1388       by (simp add: power_mult[symmetric])
  1389     also from that 2 have "\<dots> \<le> (power_up p x (Suc n div 2))\<^sup>2"
  1390       by (auto intro: power_mono simp del: odd_Suc_div_two)
  1391     finally show ?thesis
  1392       using that by (auto intro!: truncate_up_le simp del: odd_Suc_div_two)
  1393   qed
  1394   then show ?case
  1395     by (auto intro!: truncate_up_le mult_left_mono 2)
  1396 qed simp
  1397 
  1398 lemmas power_up_le = order_trans[OF _ power_up]
  1399   and power_up_less = less_le_trans[OF _ power_up]
  1400   and power_down_le = order_trans[OF power_down]
  1401 
  1402 lemma power_down_fl: "0 \<le> x \<Longrightarrow> power_down_fl p x n \<le> x ^ n"
  1403   by transfer (rule power_down)
  1404 
  1405 lemma power_up_fl: "0 \<le> x \<Longrightarrow> x ^ n \<le> power_up_fl p x n"
  1406   by transfer (rule power_up)
  1407 
  1408 lemma real_power_up_fl: "real_of_float (power_up_fl p x n) = power_up p x n"
  1409   by transfer simp
  1410 
  1411 lemma real_power_down_fl: "real_of_float (power_down_fl p x n) = power_down p x n"
  1412   by transfer simp
  1413 
  1414 
  1415 subsection \<open>Approximate Addition\<close>
  1416 
  1417 definition "plus_down prec x y = truncate_down prec (x + y)"
  1418 
  1419 definition "plus_up prec x y = truncate_up prec (x + y)"
  1420 
  1421 lemma float_plus_down_float[intro, simp]: "x \<in> float \<Longrightarrow> y \<in> float \<Longrightarrow> plus_down p x y \<in> float"
  1422   by (simp add: plus_down_def)
  1423 
  1424 lemma float_plus_up_float[intro, simp]: "x \<in> float \<Longrightarrow> y \<in> float \<Longrightarrow> plus_up p x y \<in> float"
  1425   by (simp add: plus_up_def)
  1426 
  1427 lift_definition float_plus_down :: "nat \<Rightarrow> float \<Rightarrow> float \<Rightarrow> float" is plus_down ..
  1428 
  1429 lift_definition float_plus_up :: "nat \<Rightarrow> float \<Rightarrow> float \<Rightarrow> float" is plus_up ..
  1430 
  1431 lemma plus_down: "plus_down prec x y \<le> x + y"
  1432   and plus_up: "x + y \<le> plus_up prec x y"
  1433   by (auto simp: plus_down_def truncate_down plus_up_def truncate_up)
  1434 
  1435 lemma float_plus_down: "real_of_float (float_plus_down prec x y) \<le> x + y"
  1436   and float_plus_up: "x + y \<le> real_of_float (float_plus_up prec x y)"
  1437   by (transfer; rule plus_down plus_up)+
  1438 
  1439 lemmas plus_down_le = order_trans[OF plus_down]
  1440   and plus_up_le = order_trans[OF _ plus_up]
  1441   and float_plus_down_le = order_trans[OF float_plus_down]
  1442   and float_plus_up_le = order_trans[OF _ float_plus_up]
  1443 
  1444 lemma compute_plus_up[code]: "plus_up p x y = - plus_down p (-x) (-y)"
  1445   using truncate_down_uminus_eq[of p "x + y"]
  1446   by (auto simp: plus_down_def plus_up_def)
  1447 
  1448 lemma truncate_down_log2_eqI:
  1449   assumes "\<lfloor>log 2 \<bar>x\<bar>\<rfloor> = \<lfloor>log 2 \<bar>y\<bar>\<rfloor>"
  1450   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>"
  1451   shows "truncate_down p x = truncate_down p y"
  1452   using assms by (auto simp: truncate_down_def round_down_def)
  1453 
  1454 lemma sum_neq_zeroI:
  1455   "\<bar>a\<bar> \<ge> k \<Longrightarrow> \<bar>b\<bar> < k \<Longrightarrow> a + b \<noteq> 0"
  1456   "\<bar>a\<bar> > k \<Longrightarrow> \<bar>b\<bar> \<le> k \<Longrightarrow> a + b \<noteq> 0"
  1457   for a k :: real
  1458   by auto
  1459 
  1460 lemma abs_real_le_2_powr_bitlen[simp]: "\<bar>real_of_int m2\<bar> < 2 powr real_of_int (bitlen \<bar>m2\<bar>)"
  1461 proof (cases "m2 = 0")
  1462   case True
  1463   then show ?thesis by simp
  1464 next
  1465   case False
  1466   then have "\<bar>m2\<bar> < 2 ^ nat (bitlen \<bar>m2\<bar>)"
  1467     using bitlen_bounds[of "\<bar>m2\<bar>"]
  1468     by (auto simp: powr_add bitlen_nonneg)
  1469   then show ?thesis
  1470     by (metis bitlen_nonneg powr_int of_int_abs real_of_int_less_numeral_power_cancel_iff
  1471       zero_less_numeral)
  1472 qed
  1473 
  1474 lemma floor_sum_times_2_powr_sgn_eq:
  1475   fixes ai p q :: int
  1476     and a b :: real
  1477   assumes "a * 2 powr p = ai"
  1478     and b_le_1: "\<bar>b * 2 powr (p + 1)\<bar> \<le> 1"
  1479     and leqp: "q \<le> p"
  1480   shows "\<lfloor>(a + b) * 2 powr q\<rfloor> = \<lfloor>(2 * ai + sgn b) * 2 powr (q - p - 1)\<rfloor>"
  1481 proof -
  1482   consider "b = 0" | "b > 0" | "b < 0" by arith
  1483   then show ?thesis
  1484   proof cases
  1485     case 1
  1486     then show ?thesis
  1487       by (simp add: assms(1)[symmetric] powr_add[symmetric] algebra_simps powr_mult_base)
  1488   next
  1489     case 2
  1490     then have "b * 2 powr p < \<bar>b * 2 powr (p + 1)\<bar>"
  1491       by simp
  1492     also note b_le_1
  1493     finally have b_less_1: "b * 2 powr real_of_int p < 1" .
  1494 
  1495     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"
  1496       by (simp_all add: floor_eq_iff)
  1497 
  1498     have "\<lfloor>(a + b) * 2 powr q\<rfloor> = \<lfloor>(a + b) * 2 powr p * 2 powr (q - p)\<rfloor>"
  1499       by (simp add: algebra_simps powr_realpow[symmetric] powr_add[symmetric])
  1500     also have "\<dots> = \<lfloor>(ai + b * 2 powr p) * 2 powr (q - p)\<rfloor>"
  1501       by (simp add: assms algebra_simps)
  1502     also have "\<dots> = \<lfloor>(ai + b * 2 powr p) / real_of_int ((2::int) ^ nat (p - q))\<rfloor>"
  1503       using assms
  1504       by (simp add: algebra_simps powr_realpow[symmetric] divide_powr_uminus powr_add[symmetric])
  1505     also have "\<dots> = \<lfloor>ai / real_of_int ((2::int) ^ nat (p - q))\<rfloor>"
  1506       by (simp del: of_int_power add: floor_divide_real_eq_div floor_eq)
  1507     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>" .
  1508     moreover
  1509     have "\<lfloor>(2 * ai + (sgn b)) * 2 powr (real_of_int (q - p) - 1)\<rfloor> =
  1510         \<lfloor>real_of_int ai / real_of_int ((2::int) ^ nat (p - q))\<rfloor>"
  1511     proof -
  1512       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>"
  1513         by (subst powr_diff) (simp add: field_simps)
  1514       also have "\<dots> = \<lfloor>(ai + sgn b / 2) / real_of_int ((2::int) ^ nat (p - q))\<rfloor>"
  1515         using leqp by (simp add: powr_realpow[symmetric] powr_diff)
  1516       also have "\<dots> = \<lfloor>ai / real_of_int ((2::int) ^ nat (p - q))\<rfloor>"
  1517         by (simp del: of_int_power add: floor_divide_real_eq_div floor_eq)
  1518       finally show ?thesis .
  1519     qed
  1520     ultimately show ?thesis by simp
  1521   next
  1522     case 3
  1523     then have floor_eq: "\<lfloor>b * 2 powr (real_of_int p + 1)\<rfloor> = -1"
  1524       using b_le_1
  1525       by (auto simp: floor_eq_iff algebra_simps pos_divide_le_eq[symmetric] abs_if divide_powr_uminus
  1526         intro!: mult_neg_pos split: if_split_asm)
  1527     have "\<lfloor>(a + b) * 2 powr q\<rfloor> = \<lfloor>(2*a + 2*b) * 2 powr p * 2 powr (q - p - 1)\<rfloor>"
  1528       by (simp add: algebra_simps powr_realpow[symmetric] powr_add[symmetric] powr_mult_base)
  1529     also have "\<dots> = \<lfloor>(2 * (a * 2 powr p) + 2 * b * 2 powr p) * 2 powr (q - p - 1)\<rfloor>"
  1530       by (simp add: algebra_simps)
  1531     also have "\<dots> = \<lfloor>(2 * ai + b * 2 powr (p + 1)) / 2 powr (1 - q + p)\<rfloor>"
  1532       using assms by (simp add: algebra_simps powr_mult_base divide_powr_uminus)
  1533     also have "\<dots> = \<lfloor>(2 * ai + b * 2 powr (p + 1)) / real_of_int ((2::int) ^ nat (p - q + 1))\<rfloor>"
  1534       using assms by (simp add: algebra_simps powr_realpow[symmetric])
  1535     also have "\<dots> = \<lfloor>(2 * ai - 1) / real_of_int ((2::int) ^ nat (p - q + 1))\<rfloor>"
  1536       using \<open>b < 0\<close> assms
  1537       by (simp add: floor_divide_of_int_eq floor_eq floor_divide_real_eq_div
  1538         del: of_int_mult of_int_power of_int_diff)
  1539     also have "\<dots> = \<lfloor>(2 * ai - 1) * 2 powr (q - p - 1)\<rfloor>"
  1540       using assms by (simp add: algebra_simps divide_powr_uminus powr_realpow[symmetric])
  1541     finally show ?thesis
  1542       using \<open>b < 0\<close> by simp
  1543   qed
  1544 qed
  1545 
  1546 lemma log2_abs_int_add_less_half_sgn_eq:
  1547   fixes ai :: int
  1548     and b :: real
  1549   assumes "\<bar>b\<bar> \<le> 1/2"
  1550     and "ai \<noteq> 0"
  1551   shows "\<lfloor>log 2 \<bar>real_of_int ai + b\<bar>\<rfloor> = \<lfloor>log 2 \<bar>ai + sgn b / 2\<bar>\<rfloor>"
  1552 proof (cases "b = 0")
  1553   case True
  1554   then show ?thesis by simp
  1555 next
  1556   case False
  1557   define k where "k = \<lfloor>log 2 \<bar>ai\<bar>\<rfloor>"
  1558   then have "\<lfloor>log 2 \<bar>ai\<bar>\<rfloor> = k"
  1559     by simp
  1560   then have k: "2 powr k \<le> \<bar>ai\<bar>" "\<bar>ai\<bar> < 2 powr (k + 1)"
  1561     by (simp_all add: floor_log_eq_powr_iff \<open>ai \<noteq> 0\<close>)
  1562   have "k \<ge> 0"
  1563     using assms by (auto simp: k_def)
  1564   define r where "r = \<bar>ai\<bar> - 2 ^ nat k"
  1565   have r: "0 \<le> r" "r < 2 powr k"
  1566     using \<open>k \<ge> 0\<close> k
  1567     by (auto simp: r_def k_def algebra_simps powr_add abs_if powr_int)
  1568   then have "r \<le> (2::int) ^ nat k - 1"
  1569     using \<open>k \<ge> 0\<close> by (auto simp: powr_int)
  1570   from this[simplified of_int_le_iff[symmetric]] \<open>0 \<le> k\<close>
  1571   have r_le: "r \<le> 2 powr k - 1"
  1572     by (auto simp: algebra_simps powr_int)
  1573       (metis of_int_1 of_int_add real_of_int_le_numeral_power_cancel_iff)
  1574 
  1575   have "\<bar>ai\<bar> = 2 powr k + r"
  1576     using \<open>k \<ge> 0\<close> by (auto simp: k_def r_def powr_realpow[symmetric])
  1577 
  1578   have pos: "\<bar>b\<bar> < 1 \<Longrightarrow> 0 < 2 powr k + (r + b)" for b :: real
  1579     using \<open>0 \<le> k\<close> \<open>ai \<noteq> 0\<close>
  1580     by (auto simp add: r_def powr_realpow[symmetric] abs_if sgn_if algebra_simps
  1581       split: if_split_asm)
  1582   have less: "\<bar>sgn ai * b\<bar> < 1"
  1583     and less': "\<bar>sgn (sgn ai * b) / 2\<bar> < 1"
  1584     using \<open>\<bar>b\<bar> \<le> _\<close> by (auto simp: abs_if sgn_if split: if_split_asm)
  1585 
  1586   have floor_eq: "\<And>b::real. \<bar>b\<bar> \<le> 1 / 2 \<Longrightarrow>
  1587       \<lfloor>log 2 (1 + (r + b) / 2 powr k)\<rfloor> = (if r = 0 \<and> b < 0 then -1 else 0)"
  1588     using \<open>k \<ge> 0\<close> r r_le
  1589     by (auto simp: floor_log_eq_powr_iff powr_minus_divide field_simps sgn_if)
  1590 
  1591   from \<open>real_of_int \<bar>ai\<bar> = _\<close> have "\<bar>ai + b\<bar> = 2 powr k + (r + sgn ai * b)"
  1592     using \<open>\<bar>b\<bar> \<le> _\<close> \<open>0 \<le> k\<close> r
  1593     by (auto simp add: sgn_if abs_if)
  1594   also have "\<lfloor>log 2 \<dots>\<rfloor> = \<lfloor>log 2 (2 powr k + r + sgn (sgn ai * b) / 2)\<rfloor>"
  1595   proof -
  1596     have "2 powr k + (r + (sgn ai) * b) = 2 powr k * (1 + (r + sgn ai * b) / 2 powr k)"
  1597       by (simp add: field_simps)
  1598     also have "\<lfloor>log 2 \<dots>\<rfloor> = k + \<lfloor>log 2 (1 + (r + sgn ai * b) / 2 powr k)\<rfloor>"
  1599       using pos[OF less]
  1600       by (subst log_mult) (simp_all add: log_mult powr_mult field_simps)
  1601     also
  1602     let ?if = "if r = 0 \<and> sgn ai * b < 0 then -1 else 0"
  1603     have "\<lfloor>log 2 (1 + (r + sgn ai * b) / 2 powr k)\<rfloor> = ?if"
  1604       using \<open>\<bar>b\<bar> \<le> _\<close>
  1605       by (intro floor_eq) (auto simp: abs_mult sgn_if)
  1606     also
  1607     have "\<dots> = \<lfloor>log 2 (1 + (r + sgn (sgn ai * b) / 2) / 2 powr k)\<rfloor>"
  1608       by (subst floor_eq) (auto simp: sgn_if)
  1609     also have "k + \<dots> = \<lfloor>log 2 (2 powr k * (1 + (r + sgn (sgn ai * b) / 2) / 2 powr k))\<rfloor>"
  1610       unfolding int_add_floor
  1611       using pos[OF less'] \<open>\<bar>b\<bar> \<le> _\<close>
  1612       by (simp add: field_simps add_log_eq_powr del: floor_add2)
  1613     also have "2 powr k * (1 + (r + sgn (sgn ai * b) / 2) / 2 powr k) =
  1614         2 powr k + r + sgn (sgn ai * b) / 2"
  1615       by (simp add: sgn_if field_simps)
  1616     finally show ?thesis .
  1617   qed
  1618   also have "2 powr k + r + sgn (sgn ai * b) / 2 = \<bar>ai + sgn b / 2\<bar>"
  1619     unfolding \<open>real_of_int \<bar>ai\<bar> = _\<close>[symmetric] using \<open>ai \<noteq> 0\<close>
  1620     by (auto simp: abs_if sgn_if algebra_simps)
  1621   finally show ?thesis .
  1622 qed
  1623 
  1624 context
  1625 begin
  1626 
  1627 qualified lemma compute_far_float_plus_down:
  1628   fixes m1 e1 m2 e2 :: int
  1629     and p :: nat
  1630   defines "k1 \<equiv> Suc p - nat (bitlen \<bar>m1\<bar>)"
  1631   assumes H: "bitlen \<bar>m2\<bar> \<le> e1 - e2 - k1 - 2" "m1 \<noteq> 0" "m2 \<noteq> 0" "e1 \<ge> e2"
  1632   shows "float_plus_down p (Float m1 e1) (Float m2 e2) =
  1633     float_round_down p (Float (m1 * 2 ^ (Suc (Suc k1)) + sgn m2) (e1 - int k1 - 2))"
  1634 proof -
  1635   let ?a = "real_of_float (Float m1 e1)"
  1636   let ?b = "real_of_float (Float m2 e2)"
  1637   let ?sum = "?a + ?b"
  1638   let ?shift = "real_of_int e2 - real_of_int e1 + real k1 + 1"
  1639   let ?m1 = "m1 * 2 ^ Suc k1"
  1640   let ?m2 = "m2 * 2 powr ?shift"
  1641   let ?m2' = "sgn m2 / 2"
  1642   let ?e = "e1 - int k1 - 1"
  1643 
  1644   have sum_eq: "?sum = (?m1 + ?m2) * 2 powr ?e"
  1645     by (auto simp: powr_add[symmetric] powr_mult[symmetric] algebra_simps
  1646       powr_realpow[symmetric] powr_mult_base)
  1647 
  1648   have "\<bar>?m2\<bar> * 2 < 2 powr (bitlen \<bar>m2\<bar> + ?shift + 1)"
  1649     by (auto simp: field_simps powr_add powr_mult_base powr_diff abs_mult)
  1650   also have "\<dots> \<le> 2 powr 0"
  1651     using H by (intro powr_mono) auto
  1652   finally have abs_m2_less_half: "\<bar>?m2\<bar> < 1 / 2"
  1653     by simp
  1654 
  1655   then have "\<bar>real_of_int m2\<bar> < 2 powr -(?shift + 1)"
  1656     unfolding powr_minus_divide by (auto simp: bitlen_alt_def field_simps powr_mult_base abs_mult)
  1657   also have "\<dots> \<le> 2 powr real_of_int (e1 - e2 - 2)"
  1658     by simp
  1659   finally have b_less_quarter: "\<bar>?b\<bar> < 1/4 * 2 powr real_of_int e1"
  1660     by (simp add: powr_add field_simps powr_diff abs_mult)
  1661   also have "1/4 < \<bar>real_of_int m1\<bar> / 2" using \<open>m1 \<noteq> 0\<close> by simp
  1662   finally have b_less_half_a: "\<bar>?b\<bar> < 1/2 * \<bar>?a\<bar>"
  1663     by (simp add: algebra_simps powr_mult_base abs_mult)
  1664   then have a_half_less_sum: "\<bar>?a\<bar> / 2 < \<bar>?sum\<bar>"
  1665     by (auto simp: field_simps abs_if split: if_split_asm)
  1666 
  1667   from b_less_half_a have "\<bar>?b\<bar> < \<bar>?a\<bar>" "\<bar>?b\<bar> \<le> \<bar>?a\<bar>"
  1668     by simp_all
  1669 
  1670   have "\<bar>real_of_float (Float m1 e1)\<bar> \<ge> 1/4 * 2 powr real_of_int e1"
  1671     using \<open>m1 \<noteq> 0\<close>
  1672     by (auto simp: powr_add powr_int bitlen_nonneg divide_right_mono abs_mult)
  1673   then have "?sum \<noteq> 0" using b_less_quarter
  1674     by (rule sum_neq_zeroI)
  1675   then have "?m1 + ?m2 \<noteq> 0"
  1676     unfolding sum_eq by (simp add: abs_mult zero_less_mult_iff)
  1677 
  1678   have "\<bar>real_of_int ?m1\<bar> \<ge> 2 ^ Suc k1" "\<bar>?m2'\<bar> < 2 ^ Suc k1"
  1679     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)
  1680   then have sum'_nz: "?m1 + ?m2' \<noteq> 0"
  1681     by (intro sum_neq_zeroI)
  1682 
  1683   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"
  1684     using \<open>?m1 + ?m2 \<noteq> 0\<close>
  1685     unfolding floor_add[symmetric] sum_eq
  1686     by (simp add: abs_mult log_mult) linarith
  1687   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>"
  1688     using abs_m2_less_half \<open>m1 \<noteq> 0\<close>
  1689     by (intro log2_abs_int_add_less_half_sgn_eq) (auto simp: abs_mult)
  1690   also have "sgn (real_of_int m2 * 2 powr ?shift) = sgn m2"
  1691     by (auto simp: sgn_if zero_less_mult_iff less_not_sym)
  1692   also
  1693   have "\<bar>?m1 + ?m2'\<bar> * 2 powr ?e = \<bar>?m1 * 2 + sgn m2\<bar> * 2 powr (?e - 1)"
  1694     by (auto simp: field_simps powr_minus[symmetric] powr_diff powr_mult_base)
  1695   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>"
  1696     using \<open>?m1 + ?m2' \<noteq> 0\<close>
  1697     unfolding floor_add_int
  1698     by (simp add: log_add_eq_powr abs_mult_pos del: floor_add2)
  1699   finally
  1700   have "\<lfloor>log 2 \<bar>?sum\<bar>\<rfloor> = \<lfloor>log 2 \<bar>real_of_float (Float (?m1*2 + sgn m2) (?e - 1))\<bar>\<rfloor>" .
  1701   then have "plus_down p (Float m1 e1) (Float m2 e2) =
  1702       truncate_down p (Float (?m1*2 + sgn m2) (?e - 1))"
  1703     unfolding plus_down_def
  1704   proof (rule truncate_down_log2_eqI)
  1705     let ?f = "(int p - \<lfloor>log 2 \<bar>real_of_float (Float m1 e1) + real_of_float (Float m2 e2)\<bar>\<rfloor>)"
  1706     let ?ai = "m1 * 2 ^ (Suc k1)"
  1707     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>"
  1708     proof (rule floor_sum_times_2_powr_sgn_eq)
  1709       show "?a * 2 powr real_of_int (-?e) = real_of_int ?ai"
  1710         by (simp add: powr_add powr_realpow[symmetric] powr_diff)
  1711       show "\<bar>?b * 2 powr real_of_int (-?e + 1)\<bar> \<le> 1"
  1712         using abs_m2_less_half
  1713         by (simp add: abs_mult powr_add[symmetric] algebra_simps powr_mult_base)
  1714     next
  1715       have "e1 + \<lfloor>log 2 \<bar>real_of_int m1\<bar>\<rfloor> - 1 = \<lfloor>log 2 \<bar>?a\<bar>\<rfloor> - 1"
  1716         using \<open>m1 \<noteq> 0\<close>
  1717         by (simp add: int_add_floor algebra_simps log_mult abs_mult del: floor_add2)
  1718       also have "\<dots> \<le> \<lfloor>log 2 \<bar>?a + ?b\<bar>\<rfloor>"
  1719         using a_half_less_sum \<open>m1 \<noteq> 0\<close> \<open>?sum \<noteq> 0\<close>
  1720         unfolding floor_diff_of_int[symmetric]
  1721         by (auto simp add: log_minus_eq_powr powr_minus_divide intro!: floor_mono)
  1722       finally
  1723       have "int p - \<lfloor>log 2 \<bar>?a + ?b\<bar>\<rfloor> \<le> p - (bitlen \<bar>m1\<bar>) - e1 + 2"
  1724         by (auto simp: algebra_simps bitlen_alt_def \<open>m1 \<noteq> 0\<close>)
  1725       also have "\<dots> \<le> - ?e"
  1726         using bitlen_nonneg[of "\<bar>m1\<bar>"] by (simp add: k1_def)
  1727       finally show "?f \<le> - ?e" by simp
  1728     qed
  1729     also have "sgn ?b = sgn m2"
  1730       using powr_gt_zero[of 2 e2]
  1731       by (auto simp add: sgn_if zero_less_mult_iff simp del: powr_gt_zero)
  1732     also have "\<lfloor>(real_of_int (2 * ?m1) + real_of_int (sgn m2)) * 2 powr real_of_int (?f - - ?e - 1)\<rfloor> =
  1733         \<lfloor>Float (?m1 * 2 + sgn m2) (?e - 1) * 2 powr ?f\<rfloor>"
  1734       by (simp add: powr_add[symmetric] algebra_simps powr_realpow[symmetric])
  1735     finally
  1736     show "\<lfloor>(?a + ?b) * 2 powr ?f\<rfloor> = \<lfloor>real_of_float (Float (?m1 * 2 + sgn m2) (?e - 1)) * 2 powr ?f\<rfloor>" .
  1737   qed
  1738   then show ?thesis
  1739     by transfer (simp add: plus_down_def ac_simps Let_def)
  1740 qed
  1741 
  1742 lemma compute_float_plus_down_naive[code]: "float_plus_down p x y = float_round_down p (x + y)"
  1743   by transfer (auto simp: plus_down_def)
  1744 
  1745 qualified lemma compute_float_plus_down[code]:
  1746   fixes p::nat and m1 e1 m2 e2::int
  1747   shows "float_plus_down p (Float m1 e1) (Float m2 e2) =
  1748     (if m1 = 0 then float_round_down p (Float m2 e2)
  1749     else if m2 = 0 then float_round_down p (Float m1 e1)
  1750     else
  1751       (if e1 \<ge> e2 then
  1752         (let k1 = Suc p - nat (bitlen \<bar>m1\<bar>) in
  1753           if bitlen \<bar>m2\<bar> > e1 - e2 - k1 - 2
  1754           then float_round_down p ((Float m1 e1) + (Float m2 e2))
  1755           else float_round_down p (Float (m1 * 2 ^ (Suc (Suc k1)) + sgn m2) (e1 - int k1 - 2)))
  1756     else float_plus_down p (Float m2 e2) (Float m1 e1)))"
  1757 proof -
  1758   {
  1759     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"
  1760     note compute_far_float_plus_down[OF this]
  1761   }
  1762   then show ?thesis
  1763     by transfer (simp add: Let_def plus_down_def ac_simps)
  1764 qed
  1765 
  1766 qualified lemma compute_float_plus_up[code]: "float_plus_up p x y = - float_plus_down p (-x) (-y)"
  1767   using truncate_down_uminus_eq[of p "x + y"]
  1768   by transfer (simp add: plus_down_def plus_up_def ac_simps)
  1769 
  1770 lemma mantissa_zero[simp]: "mantissa 0 = 0"
  1771   by (metis mantissa_0 zero_float.abs_eq)
  1772 
  1773 qualified lemma compute_float_less[code]: "a < b \<longleftrightarrow> is_float_pos (float_plus_down 0 b (- a))"
  1774   using truncate_down[of 0 "b - a"] truncate_down_pos[of "b - a" 0]
  1775   by transfer (auto simp: plus_down_def)
  1776 
  1777 qualified lemma compute_float_le[code]: "a \<le> b \<longleftrightarrow> is_float_nonneg (float_plus_down 0 b (- a))"
  1778   using truncate_down[of 0 "b - a"] truncate_down_nonneg[of "b - a" 0]
  1779   by transfer (auto simp: plus_down_def)
  1780 
  1781 end
  1782 
  1783 
  1784 subsection \<open>Lemmas needed by Approximate\<close>
  1785 
  1786 lemma Float_num[simp]:
  1787    "real_of_float (Float 1 0) = 1"
  1788    "real_of_float (Float 1 1) = 2"
  1789    "real_of_float (Float 1 2) = 4"
  1790    "real_of_float (Float 1 (- 1)) = 1/2"
  1791    "real_of_float (Float 1 (- 2)) = 1/4"
  1792    "real_of_float (Float 1 (- 3)) = 1/8"
  1793    "real_of_float (Float (- 1) 0) = -1"
  1794    "real_of_float (Float (numeral n) 0) = numeral n"
  1795    "real_of_float (Float (- numeral n) 0) = - numeral n"
  1796   using two_powr_int_float[of 2] two_powr_int_float[of "-1"] two_powr_int_float[of "-2"]
  1797     two_powr_int_float[of "-3"]
  1798   using powr_realpow[of 2 2] powr_realpow[of 2 3]
  1799   using powr_minus[of "2::real" 1] powr_minus[of "2::real" 2] powr_minus[of "2::real" 3]
  1800   by auto
  1801 
  1802 lemma real_of_Float_int[simp]: "real_of_float (Float n 0) = real n"
  1803   by simp
  1804 
  1805 lemma float_zero[simp]: "real_of_float (Float 0 e) = 0"
  1806   by simp
  1807 
  1808 lemma abs_div_2_less: "a \<noteq> 0 \<Longrightarrow> a \<noteq> -1 \<Longrightarrow> \<bar>(a::int) div 2\<bar> < \<bar>a\<bar>"
  1809   by arith
  1810 
  1811 lemma lapprox_rat: "real_of_float (lapprox_rat prec x y) \<le> real_of_int x / real_of_int y"
  1812   by (simp add: lapprox_rat.rep_eq truncate_down)
  1813 
  1814 lemma mult_div_le:
  1815   fixes a b :: int
  1816   assumes "b > 0"
  1817   shows "a \<ge> b * (a div b)"
  1818 proof -
  1819   from minus_div_mult_eq_mod [symmetric, of a b] have "a = b * (a div b) + a mod b"
  1820     by simp
  1821   also have "\<dots> \<ge> b * (a div b) + 0"
  1822     apply (rule add_left_mono)
  1823     apply (rule pos_mod_sign)
  1824     using assms
  1825     apply simp
  1826     done
  1827   finally show ?thesis
  1828     by simp
  1829 qed
  1830 
  1831 lemma lapprox_rat_nonneg:
  1832   assumes "0 \<le> x" and "0 \<le> y"
  1833   shows "0 \<le> real_of_float (lapprox_rat n x y)"
  1834   using assms
  1835   by transfer (simp add: truncate_down_nonneg)
  1836 
  1837 lemma rapprox_rat: "real_of_int x / real_of_int y \<le> real_of_float (rapprox_rat prec x y)"
  1838   by transfer (simp add: truncate_up)
  1839 
  1840 lemma rapprox_rat_le1:
  1841   assumes "0 \<le> x" "0 < y" "x \<le> y"
  1842   shows "real_of_float (rapprox_rat n x y) \<le> 1"
  1843   using assms
  1844   by transfer (simp add: truncate_up_le1)
  1845 
  1846 lemma rapprox_rat_nonneg_nonpos: "0 \<le> x \<Longrightarrow> y \<le> 0 \<Longrightarrow> real_of_float (rapprox_rat n x y) \<le> 0"
  1847   by transfer (simp add: truncate_up_nonpos divide_nonneg_nonpos)
  1848 
  1849 lemma rapprox_rat_nonpos_nonneg: "x \<le> 0 \<Longrightarrow> 0 \<le> y \<Longrightarrow> real_of_float (rapprox_rat n x y) \<le> 0"
  1850   by transfer (simp add: truncate_up_nonpos divide_nonpos_nonneg)
  1851 
  1852 lemma real_divl: "real_divl prec x y \<le> x / y"
  1853   by (simp add: real_divl_def truncate_down)
  1854 
  1855 lemma real_divr: "x / y \<le> real_divr prec x y"
  1856   by (simp add: real_divr_def truncate_up)
  1857 
  1858 lemma float_divl: "real_of_float (float_divl prec x y) \<le> x / y"
  1859   by transfer (rule real_divl)
  1860 
  1861 lemma real_divl_lower_bound: "0 \<le> x \<Longrightarrow> 0 \<le> y \<Longrightarrow> 0 \<le> real_divl prec x y"
  1862   by (simp add: real_divl_def truncate_down_nonneg)
  1863 
  1864 lemma float_divl_lower_bound: "0 \<le> x \<Longrightarrow> 0 \<le> y \<Longrightarrow> 0 \<le> real_of_float (float_divl prec x y)"
  1865   by transfer (rule real_divl_lower_bound)
  1866 
  1867 lemma exponent_1: "exponent 1 = 0"
  1868   using exponent_float[of 1 0] by (simp add: one_float_def)
  1869 
  1870 lemma mantissa_1: "mantissa 1 = 1"
  1871   using mantissa_float[of 1 0] by (simp add: one_float_def)
  1872 
  1873 lemma bitlen_1: "bitlen 1 = 1"
  1874   by (simp add: bitlen_alt_def)
  1875 
  1876 lemma mantissa_eq_zero_iff: "mantissa x = 0 \<longleftrightarrow> x = 0"
  1877   (is "?lhs \<longleftrightarrow> ?rhs")
  1878 proof
  1879   show ?rhs if ?lhs
  1880   proof -
  1881     from that have z: "0 = real_of_float x"
  1882       using mantissa_exponent by simp
  1883     show ?thesis
  1884       by (simp add: zero_float_def z)
  1885   qed
  1886   show ?lhs if ?rhs
  1887     using that by (simp add: zero_float_def)
  1888 qed
  1889 
  1890 lemma float_upper_bound: "x \<le> 2 powr (bitlen \<bar>mantissa x\<bar> + exponent x)"
  1891 proof (cases "x = 0")
  1892   case True
  1893   then show ?thesis by simp
  1894 next
  1895   case False
  1896   then have "mantissa x \<noteq> 0"
  1897     using mantissa_eq_zero_iff by auto
  1898   have "x = mantissa x * 2 powr (exponent x)"
  1899     by (rule mantissa_exponent)
  1900   also have "mantissa x \<le> \<bar>mantissa x\<bar>"
  1901     by simp
  1902   also have "\<dots> \<le> 2 powr (bitlen \<bar>mantissa x\<bar>)"
  1903     using bitlen_bounds[of "\<bar>mantissa x\<bar>"] bitlen_nonneg \<open>mantissa x \<noteq> 0\<close>
  1904     by (auto simp del: of_int_abs simp add: powr_int)
  1905   finally show ?thesis by (simp add: powr_add)
  1906 qed
  1907 
  1908 lemma real_divl_pos_less1_bound:
  1909   assumes "0 < x" "x \<le> 1"
  1910   shows "1 \<le> real_divl prec 1 x"
  1911   using assms
  1912   by (auto intro!: truncate_down_ge1 simp: real_divl_def)
  1913 
  1914 lemma float_divl_pos_less1_bound:
  1915   "0 < real_of_float x \<Longrightarrow> real_of_float x \<le> 1 \<Longrightarrow> prec \<ge> 1 \<Longrightarrow>
  1916     1 \<le> real_of_float (float_divl prec 1 x)"
  1917   by transfer (rule real_divl_pos_less1_bound)
  1918 
  1919 lemma float_divr: "real_of_float x / real_of_float y \<le> real_of_float (float_divr prec x y)"
  1920   by transfer (rule real_divr)
  1921 
  1922 lemma real_divr_pos_less1_lower_bound:
  1923   assumes "0 < x"
  1924     and "x \<le> 1"
  1925   shows "1 \<le> real_divr prec 1 x"
  1926 proof -
  1927   have "1 \<le> 1 / x"
  1928     using \<open>0 < x\<close> and \<open>x \<le> 1\<close> by auto
  1929   also have "\<dots> \<le> real_divr prec 1 x"
  1930     using real_divr[where x = 1 and y = x] by auto
  1931   finally show ?thesis by auto
  1932 qed
  1933 
  1934 lemma float_divr_pos_less1_lower_bound: "0 < x \<Longrightarrow> x \<le> 1 \<Longrightarrow> 1 \<le> float_divr prec 1 x"
  1935   by transfer (rule real_divr_pos_less1_lower_bound)
  1936 
  1937 lemma real_divr_nonpos_pos_upper_bound: "x \<le> 0 \<Longrightarrow> 0 \<le> y \<Longrightarrow> real_divr prec x y \<le> 0"
  1938   by (simp add: real_divr_def truncate_up_nonpos divide_le_0_iff)
  1939 
  1940 lemma float_divr_nonpos_pos_upper_bound:
  1941   "real_of_float x \<le> 0 \<Longrightarrow> 0 \<le> real_of_float y \<Longrightarrow> real_of_float (float_divr prec x y) \<le> 0"
  1942   by transfer (rule real_divr_nonpos_pos_upper_bound)
  1943 
  1944 lemma real_divr_nonneg_neg_upper_bound: "0 \<le> x \<Longrightarrow> y \<le> 0 \<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_nonneg_neg_upper_bound:
  1948   "0 \<le> real_of_float x \<Longrightarrow> real_of_float y \<le> 0 \<Longrightarrow> real_of_float (float_divr prec x y) \<le> 0"
  1949   by transfer (rule real_divr_nonneg_neg_upper_bound)
  1950 
  1951 lemma truncate_up_nonneg_mono:
  1952   assumes "0 \<le> x" "x \<le> y"
  1953   shows "truncate_up prec x \<le> truncate_up prec y"
  1954 proof -
  1955   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"
  1956     by arith
  1957   then show ?thesis
  1958   proof cases
  1959     case 1
  1960     then show ?thesis
  1961       using assms
  1962       by (auto simp: truncate_up_def round_up_def intro!: ceiling_mono)
  1963   next
  1964     case 2
  1965     from assms \<open>0 < x\<close> have "log 2 x \<le> log 2 y"
  1966       by auto
  1967     with \<open>\<lfloor>log 2 x\<rfloor> \<noteq> \<lfloor>log 2 y\<rfloor>\<close>
  1968     have logless: "log 2 x < log 2 y" and flogless: "\<lfloor>log 2 x\<rfloor> < \<lfloor>log 2 y\<rfloor>"
  1969       by (metis floor_less_cancel linorder_cases not_le)+
  1970     have "truncate_up prec x =
  1971       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>)"
  1972       using assms by (simp add: truncate_up_def round_up_def)
  1973     also have "\<lceil>x * 2 powr real_of_int (int prec - \<lfloor>log 2 x\<rfloor>)\<rceil> \<le> (2 ^ (Suc prec))"
  1974     proof (simp only: ceiling_le_iff)
  1975       have "x * 2 powr real_of_int (int prec - \<lfloor>log 2 x\<rfloor>) \<le>
  1976         x * (2 powr real (Suc prec) / (2 powr log 2 x))"
  1977         using real_of_int_floor_add_one_ge[of "log 2 x"] assms
  1978         by (auto simp add: algebra_simps powr_diff [symmetric]  intro!: mult_left_mono)
  1979       then show "x * 2 powr real_of_int (int prec - \<lfloor>log 2 x\<rfloor>) \<le> real_of_int ((2::int) ^ (Suc prec))"
  1980         using \<open>0 < x\<close> by (simp add: powr_realpow powr_add)
  1981     qed
  1982     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)"
  1983       by (auto simp: powr_realpow powr_add)
  1984         (metis power_Suc real_of_int_le_numeral_power_cancel_iff)
  1985     also
  1986     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)"
  1987       using logless flogless by (auto intro!: floor_mono)
  1988     also have "2 powr real_of_int (int (Suc prec)) \<le>
  1989         2 powr (log 2 y + real_of_int (int prec - \<lfloor>log 2 y\<rfloor> + 1))"
  1990       using assms \<open>0 < x\<close>
  1991       by (auto simp: algebra_simps)
  1992     finally have "truncate_up prec x \<le>
  1993         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)"
  1994       by simp
  1995     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>))"
  1996       by (subst powr_add[symmetric]) simp
  1997     also have "\<dots> = y"
  1998       using \<open>0 < x\<close> assms
  1999       by (simp add: powr_add)
  2000     also have "\<dots> \<le> truncate_up prec y"
  2001       by (rule truncate_up)
  2002     finally show ?thesis .
  2003   next
  2004     case 3
  2005     then show ?thesis
  2006       using assms
  2007       by (auto intro!: truncate_up_le)
  2008   qed
  2009 qed
  2010 
  2011 lemma truncate_up_switch_sign_mono:
  2012   assumes "x \<le> 0" "0 \<le> y"
  2013   shows "truncate_up prec x \<le> truncate_up prec y"
  2014 proof -
  2015   note truncate_up_nonpos[OF \<open>x \<le> 0\<close>]
  2016   also note truncate_up_le[OF \<open>0 \<le> y\<close>]
  2017   finally show ?thesis .
  2018 qed
  2019 
  2020 lemma truncate_down_switch_sign_mono:
  2021   assumes "x \<le> 0"
  2022     and "0 \<le> y"
  2023     and "x \<le> y"
  2024   shows "truncate_down prec x \<le> truncate_down prec y"
  2025 proof -
  2026   note truncate_down_le[OF \<open>x \<le> 0\<close>]
  2027   also note truncate_down_nonneg[OF \<open>0 \<le> y\<close>]
  2028   finally show ?thesis .
  2029 qed
  2030 
  2031 lemma truncate_down_nonneg_mono:
  2032   assumes "0 \<le> x" "x \<le> y"
  2033   shows "truncate_down prec x \<le> truncate_down prec y"
  2034 proof -
  2035   consider "x \<le> 0" | "\<lfloor>log 2 \<bar>x\<bar>\<rfloor> = \<lfloor>log 2 \<bar>y\<bar>\<rfloor>" |
  2036     "0 < x" "\<lfloor>log 2 \<bar>x\<bar>\<rfloor> \<noteq> \<lfloor>log 2 \<bar>y\<bar>\<rfloor>"
  2037     by arith
  2038   then show ?thesis
  2039   proof cases
  2040     case 1
  2041     with assms have "x = 0" "0 \<le> y" by simp_all
  2042     then show ?thesis
  2043       by (auto intro!: truncate_down_nonneg)
  2044   next
  2045     case 2
  2046     then show ?thesis
  2047       using assms
  2048       by (auto simp: truncate_down_def round_down_def intro!: floor_mono)
  2049   next
  2050     case 3
  2051     from \<open>0 < x\<close> have "log 2 x \<le> log 2 y" "0 < y" "0 \<le> y"
  2052       using assms by auto
  2053     with \<open>\<lfloor>log 2 \<bar>x\<bar>\<rfloor> \<noteq> \<lfloor>log 2 \<bar>y\<bar>\<rfloor>\<close>
  2054     have logless: "log 2 x < log 2 y" and flogless: "\<lfloor>log 2 x\<rfloor> < \<lfloor>log 2 y\<rfloor>"
  2055       unfolding atomize_conj abs_of_pos[OF \<open>0 < x\<close>] abs_of_pos[OF \<open>0 < y\<close>]
  2056       by (metis floor_less_cancel linorder_cases not_le)
  2057     have "2 powr prec \<le> y * 2 powr real prec / (2 powr log 2 y)"
  2058       using \<open>0 < y\<close> by simp
  2059     also have "\<dots> \<le> y * 2 powr real (Suc prec) / (2 powr (real_of_int \<lfloor>log 2 y\<rfloor> + 1))"
  2060       using \<open>0 \<le> y\<close> \<open>0 \<le> x\<close> assms(2)
  2061       by (auto intro!: powr_mono divide_left_mono
  2062           simp: of_nat_diff powr_add powr_diff)
  2063     also have "\<dots> = y * 2 powr real (Suc prec) / (2 powr real_of_int \<lfloor>log 2 y\<rfloor> * 2)"
  2064       by (auto simp: powr_add)
  2065     finally have "(2 ^ prec) \<le> \<lfloor>y * 2 powr real_of_int (int (Suc prec) - \<lfloor>log 2 \<bar>y\<bar>\<rfloor> - 1)\<rfloor>"
  2066       using \<open>0 \<le> y\<close>
  2067       by (auto simp: powr_diff le_floor_iff powr_realpow powr_add)
  2068     then have "(2 ^ (prec)) * 2 powr - real_of_int (int prec - \<lfloor>log 2 \<bar>y\<bar>\<rfloor>) \<le> truncate_down prec y"
  2069       by (auto simp: truncate_down_def round_down_def)
  2070     moreover have "x \<le> (2 ^ prec) * 2 powr - real_of_int (int prec - \<lfloor>log 2 \<bar>y\<bar>\<rfloor>)"
  2071     proof -
  2072       have "x = 2 powr (log 2 \<bar>x\<bar>)" using \<open>0 < x\<close> by simp
  2073       also have "\<dots> \<le> (2 ^ (Suc prec )) * 2 powr - real_of_int (int prec - \<lfloor>log 2 \<bar>x\<bar>\<rfloor>)"
  2074         using real_of_int_floor_add_one_ge[of "log 2 \<bar>x\<bar>"] \<open>0 < x\<close>
  2075         by (auto simp: powr_realpow[symmetric] powr_add[symmetric] algebra_simps
  2076           powr_mult_base le_powr_iff)
  2077       also
  2078       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)"
  2079         using logless flogless \<open>x > 0\<close> \<open>y > 0\<close>
  2080         by (auto intro!: floor_mono)
  2081       finally show ?thesis
  2082         by (auto simp: powr_realpow[symmetric] powr_diff assms of_nat_diff)
  2083     qed
  2084     ultimately show ?thesis
  2085       by (metis dual_order.trans truncate_down)
  2086   qed
  2087 qed
  2088 
  2089 lemma truncate_down_eq_truncate_up: "truncate_down p x = - truncate_up p (-x)"
  2090   and truncate_up_eq_truncate_down: "truncate_up p x = - truncate_down p (-x)"
  2091   by (auto simp: truncate_up_uminus_eq truncate_down_uminus_eq)
  2092 
  2093 lemma truncate_down_mono: "x \<le> y \<Longrightarrow> truncate_down p x \<le> truncate_down p y"
  2094   apply (cases "0 \<le> x")
  2095   apply (rule truncate_down_nonneg_mono, assumption+)
  2096   apply (simp add: truncate_down_eq_truncate_up)
  2097   apply (cases "0 \<le> y")
  2098   apply (auto intro: truncate_up_nonneg_mono truncate_up_switch_sign_mono)
  2099   done
  2100 
  2101 lemma truncate_up_mono: "x \<le> y \<Longrightarrow> truncate_up p x \<le> truncate_up p y"
  2102   by (simp add: truncate_up_eq_truncate_down truncate_down_mono)
  2103 
  2104 lemma Float_le_zero_iff: "Float a b \<le> 0 \<longleftrightarrow> a \<le> 0"
  2105  by (auto simp: zero_float_def mult_le_0_iff) (simp add: not_less [symmetric])
  2106 
  2107 lemma real_of_float_pprt[simp]:
  2108   fixes a :: float
  2109   shows "real_of_float (pprt a) = pprt (real_of_float a)"
  2110   unfolding pprt_def sup_float_def max_def sup_real_def by auto
  2111 
  2112 lemma real_of_float_nprt[simp]:
  2113   fixes a :: float
  2114   shows "real_of_float (nprt a) = nprt (real_of_float a)"
  2115   unfolding nprt_def inf_float_def min_def inf_real_def by auto
  2116 
  2117 context
  2118 begin
  2119 
  2120 lift_definition int_floor_fl :: "float \<Rightarrow> int" is floor .
  2121 
  2122 qualified lemma compute_int_floor_fl[code]:
  2123   "int_floor_fl (Float m e) = (if 0 \<le> e then m * 2 ^ nat e else m div (2 ^ (nat (-e))))"
  2124   apply transfer
  2125   apply (simp add: powr_int floor_divide_of_int_eq)
  2126   apply (metis (no_types, hide_lams)floor_divide_of_int_eq of_int_numeral of_int_power floor_of_int of_int_mult)
  2127   done
  2128 
  2129 lift_definition floor_fl :: "float \<Rightarrow> float" is "\<lambda>x. real_of_int \<lfloor>x\<rfloor>"
  2130   by simp
  2131 
  2132 qualified lemma compute_floor_fl[code]:
  2133   "floor_fl (Float m e) = (if 0 \<le> e then Float m e else Float (m div (2 ^ (nat (-e)))) 0)"
  2134   apply transfer
  2135   apply (simp add: powr_int floor_divide_of_int_eq)
  2136   apply (metis (no_types, hide_lams)floor_divide_of_int_eq of_int_numeral of_int_power of_int_mult)
  2137   done
  2138 
  2139 end
  2140 
  2141 lemma floor_fl: "real_of_float (floor_fl x) \<le> real_of_float x"
  2142   by transfer simp
  2143 
  2144 lemma int_floor_fl: "real_of_int (int_floor_fl x) \<le> real_of_float x"
  2145   by transfer simp
  2146 
  2147 lemma floor_pos_exp: "exponent (floor_fl x) \<ge> 0"
  2148 proof (cases "floor_fl x = float_of 0")
  2149   case True
  2150   then show ?thesis
  2151     by (simp add: floor_fl_def)
  2152 next
  2153   case False
  2154   have eq: "floor_fl x = Float \<lfloor>real_of_float x\<rfloor> 0"
  2155     by transfer simp
  2156   obtain i where "\<lfloor>real_of_float x\<rfloor> = mantissa (floor_fl x) * 2 ^ i" "0 = exponent (floor_fl x) - int i"
  2157     by (rule denormalize_shift[OF eq[THEN eq_reflection] False])
  2158   then show ?thesis
  2159     by simp
  2160 qed
  2161 
  2162 lemma compute_mantissa[code]:
  2163   "mantissa (Float m e) =
  2164     (if m = 0 then 0 else if 2 dvd m then mantissa (normfloat (Float m e)) else m)"
  2165   by (auto simp: mantissa_float Float.abs_eq)
  2166 
  2167 lemma compute_exponent[code]:
  2168   "exponent (Float m e) =
  2169     (if m = 0 then 0 else if 2 dvd m then exponent (normfloat (Float m e)) else e)"
  2170   by (auto simp: exponent_float Float.abs_eq)
  2171 
  2172 end