src/HOL/Library/Float.thy
 author wenzelm Tue May 15 13:57:39 2018 +0200 (16 months ago) changeset 68189 6163c90694ef parent 67573 ed0a7090167d child 68406 6beb45f6cf67 permissions -rw-r--r--
1 (*  Title:      HOL/Library/Float.thy
2     Author:     Johannes Hölzl, Fabian Immler
4 *)
6 section \<open>Floating-Point Numbers\<close>
8 theory Float
9 imports Log_Nat Lattice_Algebras
10 begin
12 definition "float = {m * 2 powr e | (m :: int) (e :: int). True}"
14 typedef float = float
15   morphisms real_of_float float_of
16   unfolding float_def by auto
18 setup_lifting type_definition_float
20 declare real_of_float [code_unfold]
22 lemmas float_of_inject[simp]
24 declare [[coercion "real_of_float :: float \<Rightarrow> real"]]
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 ..
29 declare real_of_float_inverse[simp] float_of_inverse [simp]
30 declare real_of_float [simp]
33 subsection \<open>Real operations preserving the representation as floating point number\<close>
35 lemma floatI: "m * 2 powr e = x \<Longrightarrow> x \<in> float" for m e :: int
36   by (auto simp: float_def)
38 lemma zero_float[simp]: "0 \<in> float"
39   by (auto simp: float_def)
41 lemma one_float[simp]: "1 \<in> float"
42   by (intro floatI[of 1 0]) simp
44 lemma numeral_float[simp]: "numeral i \<in> float"
45   by (intro floatI[of "numeral i" 0]) simp
47 lemma neg_numeral_float[simp]: "- numeral i \<in> float"
48   by (intro floatI[of "- numeral i" 0]) simp
50 lemma real_of_int_float[simp]: "real_of_int x \<in> float" for x :: int
51   by (intro floatI[of x 0]) simp
53 lemma real_of_nat_float[simp]: "real x \<in> float" for x :: nat
54   by (intro floatI[of x 0]) simp
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
59 lemma two_powr_nat_float[simp]: "2 powr (real i) \<in> float" for i :: nat
60   by (intro floatI[of 1 i]) simp
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
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
68 lemma two_powr_numeral_float[simp]: "2 powr numeral i \<in> float"
69   by (intro floatI[of 1 "numeral i"]) simp
71 lemma two_powr_neg_numeral_float[simp]: "2 powr - numeral i \<in> float"
72   by (intro floatI[of 1 "- numeral i"]) simp
74 lemma two_pow_float[simp]: "2 ^ n \<in> float"
75   by (intro floatI[of 1 n]) (simp add: powr_realpow)
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
96   next
97     case 2
98     then show ?thesis by (rule *)
99   qed
100 qed
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)
109   done
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)
118   done
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
123 lemma abs_float[simp]: "x \<in> float \<Longrightarrow> \<bar>x\<bar> \<in> float"
124   by (cases x rule: linorder_cases[of 0]) auto
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)
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)
136   done
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)
145   done
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
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
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
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]
184 code_datatype Float
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)))"
191 subsection \<open>Arithmetic operations on floating point numbers\<close>
193 instantiation float :: "{ring_1,linorder,linordered_ring,linordered_idom,numeral,equal}"
194 begin
196 lift_definition zero_float :: float is 0 by simp
197 declare zero_float.rep_eq[simp]
199 lift_definition one_float :: float is 1 by simp
200 declare one_float.rep_eq[simp]
202 lift_definition plus_float :: "float \<Rightarrow> float \<Rightarrow> float" is "(+)" by simp
203 declare plus_float.rep_eq[simp]
205 lift_definition times_float :: "float \<Rightarrow> float \<Rightarrow> float" is "( * )" by simp
206 declare times_float.rep_eq[simp]
208 lift_definition minus_float :: "float \<Rightarrow> float \<Rightarrow> float" is "(-)" by simp
209 declare minus_float.rep_eq[simp]
211 lift_definition uminus_float :: "float \<Rightarrow> float" is "uminus" by simp
212 declare uminus_float.rep_eq[simp]
214 lift_definition abs_float :: "float \<Rightarrow> float" is abs by simp
215 declare abs_float.rep_eq[simp]
217 lift_definition sgn_float :: "float \<Rightarrow> float" is sgn by simp
218 declare sgn_float.rep_eq[simp]
220 lift_definition equal_float :: "float \<Rightarrow> float \<Rightarrow> bool" is "(=) :: real \<Rightarrow> real \<Rightarrow> bool" .
222 lift_definition less_eq_float :: "float \<Rightarrow> float \<Rightarrow> bool" is "(\<le>)" .
223 declare less_eq_float.rep_eq[simp]
225 lift_definition less_float :: "float \<Rightarrow> float \<Rightarrow> bool" is "(<)" .
226 declare less_float.rep_eq[simp]
228 instance
229   by standard (transfer; fastforce simp add: field_simps intro: mult_left_mono mult_right_mono)+
231 end
233 lemma real_of_float [simp]: "real_of_float (of_nat n) = of_nat n"
234   by (induct n) simp_all
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)
239 lemma Float_0_eq_0[simp]: "Float 0 e = 0"
240   by transfer simp
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
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)
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
268     done
269 qed
272 begin
274 definition inf_float :: "float \<Rightarrow> float \<Rightarrow> float"
275   where "inf_float a b = min a b"
277 definition sup_float :: "float \<Rightarrow> float \<Rightarrow> float"
278   where "sup_float a b = max a b"
280 instance
281   by standard (transfer; simp add: inf_float_def sup_float_def real_of_float_min real_of_float_max)+
283 end
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
292 lemma transfer_numeral [transfer_rule]:
293   "rel_fun (=) pcr_float (numeral :: _ \<Rightarrow> real) (numeral :: _ \<Rightarrow> float)"
294   by (simp add: rel_fun_def float.pcr_cr_eq cr_float_def)
296 lemma float_neg_numeral[simp]: "real_of_float (- numeral x :: float) = - numeral x"
297   by simp
299 lemma transfer_neg_numeral [transfer_rule]:
300   "rel_fun (=) pcr_float (- numeral :: _ \<Rightarrow> real) (- numeral :: _ \<Rightarrow> float)"
301   by (simp add: rel_fun_def float.pcr_cr_eq cr_float_def)
303 lemma float_of_numeral: "numeral k = float_of (numeral k)"
304   and float_of_neg_numeral: "- numeral k = float_of (- numeral k)"
305   unfolding real_of_float_eq by simp_all
308 subsection \<open>Quickcheck\<close>
310 instantiation float :: exhaustive
311 begin
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"
317 instance ..
319 end
321 definition (in term_syntax) [code_unfold]:
322   "valtermify_float x y = Code_Evaluation.valtermify Float {\<cdot>} x {\<cdot>} y"
324 instantiation float :: full_exhaustive
325 begin
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"
332 instance ..
334 end
336 instantiation float :: random
337 begin
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)))"
343 instance ..
345 end
348 subsection \<open>Represent floats as unique mantissa and exponent\<close>
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
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
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)"
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
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
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"])
423       done
424   qed
425   with that show thesis by blast
426 qed
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
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))"
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))"
448 lemma exponent_0[simp]: "exponent 0 = 0" (is ?E)
449   and mantissa_0[simp]: "mantissa 0 = 0" (is ?M)
450 proof -
451   have "\<And>p::int \<times> int. fst p = 0 \<and> snd p = 0 \<longleftrightarrow> p = (0, 0)"
452     by auto
453   then show ?E ?M
454     by (auto simp add: mantissa_def exponent_def zero_float_def)
455 qed
457 lemma mantissa_exponent: "real_of_float f = mantissa f * 2 powr exponent f" (is ?E)
458   and mantissa_not_dvd: "f \<noteq> 0 \<Longrightarrow> \<not> 2 dvd mantissa f" (is "_ \<Longrightarrow> ?D")
459 proof cases
460   assume [simp]: "f \<noteq> 0"
461   have "f = mantissa f * 2 powr exponent f \<and> \<not> 2 dvd mantissa f"
462   proof (cases f rule: float_normed_cases)
463     case zero
464     then show ?thesis by simp
465   next
466     case (powr m e)
467     then have "\<exists>p::int \<times> int. (f = 0 \<and> fst p = 0 \<and> snd p = 0) \<or>
468       (f \<noteq> 0 \<and> real_of_float f = real_of_int (fst p) * 2 powr real_of_int (snd p) \<and> \<not> 2 dvd fst p)"
469       by auto
470     then show ?thesis
471       unfolding exponent_def mantissa_def
472       by (rule someI2_ex) simp
473   qed
474   then show ?E ?D by auto
475 qed simp
477 lemma mantissa_noteq_0: "f \<noteq> 0 \<Longrightarrow> mantissa f \<noteq> 0"
478   using mantissa_not_dvd[of f] by auto
480 lemma mantissa_eq_zero_iff: "mantissa x = 0 \<longleftrightarrow> x = 0"
481   (is "?lhs \<longleftrightarrow> ?rhs")
482 proof
483   show ?rhs if ?lhs
484   proof -
485     from that have z: "0 = real_of_float x"
486       using mantissa_exponent by simp
487     show ?thesis
488       by (simp add: zero_float_def z)
489   qed
490   show ?lhs if ?rhs
491     using that by simp
492 qed
494 lemma mantissa_pos_iff: "0 < mantissa x \<longleftrightarrow> 0 < x"
495   by (auto simp: mantissa_exponent sign_simps)
497 lemma mantissa_nonneg_iff: "0 \<le> mantissa x \<longleftrightarrow> 0 \<le> x"
498   by (auto simp: mantissa_exponent sign_simps zero_le_mult_iff)
500 lemma mantissa_neg_iff: "0 > mantissa x \<longleftrightarrow> 0 > x"
501   by (auto simp: mantissa_exponent sign_simps zero_le_mult_iff)
503 lemma
504   fixes m e :: int
505   defines "f \<equiv> float_of (m * 2 powr e)"
506   assumes dvd: "\<not> 2 dvd m"
507   shows mantissa_float: "mantissa f = m" (is "?M")
508     and exponent_float: "m \<noteq> 0 \<Longrightarrow> exponent f = e" (is "_ \<Longrightarrow> ?E")
509 proof cases
510   assume "m = 0"
511   with dvd show "mantissa f = m" by auto
512 next
513   assume "m \<noteq> 0"
514   then have f_not_0: "f \<noteq> 0" by (simp add: f_def zero_float_def)
515   from mantissa_exponent[of f] have "m * 2 powr e = mantissa f * 2 powr exponent f"
516     by (auto simp add: f_def)
517   then show ?M ?E
518     using mantissa_not_dvd[OF f_not_0] dvd
519     by (auto simp: mult_powr_eq_mult_powr_iff)
520 qed
523 subsection \<open>Compute arithmetic operations\<close>
525 lemma Float_mantissa_exponent: "Float (mantissa f) (exponent f) = f"
526   unfolding real_of_float_eq mantissa_exponent[of f] by simp
528 lemma Float_cases [cases type: float]:
529   fixes f :: float
530   obtains (Float) m e :: int where "f = Float m e"
531   using Float_mantissa_exponent[symmetric]
532   by (atomize_elim) auto
534 lemma denormalize_shift:
535   assumes f_def: "f = Float m e"
536     and not_0: "f \<noteq> 0"
537   obtains i where "m = mantissa f * 2 ^ i" "e = exponent f - i"
538 proof
539   from mantissa_exponent[of f] f_def
540   have "m * 2 powr e = mantissa f * 2 powr exponent f"
541     by simp
542   then have eq: "m = mantissa f * 2 powr (exponent f - e)"
543     by (simp add: powr_diff field_simps)
544   moreover
545   have "e \<le> exponent f"
546   proof (rule ccontr)
547     assume "\<not> e \<le> exponent f"
548     then have pos: "exponent f < e" by simp
549     then have "2 powr (exponent f - e) = 2 powr - real_of_int (e - exponent f)"
550       by simp
551     also have "\<dots> = 1 / 2^nat (e - exponent f)"
552       using pos by (simp add: powr_realpow[symmetric] powr_diff)
553     finally have "m * 2^nat (e - exponent f) = real_of_int (mantissa f)"
554       using eq by simp
555     then have "mantissa f = m * 2^nat (e - exponent f)"
556       by linarith
557     with \<open>exponent f < e\<close> have "2 dvd mantissa f"
558       apply (intro dvdI[where k="m * 2^(nat (e-exponent f)) div 2"])
559       apply (cases "nat (e - exponent f)")
560       apply auto
561       done
562     then show False using mantissa_not_dvd[OF not_0] by simp
563   qed
564   ultimately have "real_of_int m = mantissa f * 2^nat (exponent f - e)"
566   with \<open>e \<le> exponent f\<close>
567   show "m = mantissa f * 2 ^ nat (exponent f - e)"
568     by linarith
569   show "e = exponent f - nat (exponent f - e)"
570     using \<open>e \<le> exponent f\<close> by auto
571 qed
573 context
574 begin
576 qualified lemma compute_float_zero[code_unfold, code]: "0 = Float 0 0"
577   by transfer simp
579 qualified lemma compute_float_one[code_unfold, code]: "1 = Float 1 0"
580   by transfer simp
582 lift_definition normfloat :: "float \<Rightarrow> float" is "\<lambda>x. x" .
583 lemma normloat_id[simp]: "normfloat x = x" by transfer rule
585 qualified lemma compute_normfloat[code]:
586   "normfloat (Float m e) =
587     (if m mod 2 = 0 \<and> m \<noteq> 0 then normfloat (Float (m div 2) (e + 1))
588      else if m = 0 then 0 else Float m e)"
591 qualified lemma compute_float_numeral[code_abbrev]: "Float (numeral k) 0 = numeral k"
592   by transfer simp
594 qualified lemma compute_float_neg_numeral[code_abbrev]: "Float (- numeral k) 0 = - numeral k"
595   by transfer simp
597 qualified lemma compute_float_uminus[code]: "- Float m1 e1 = Float (- m1) e1"
598   by transfer simp
600 qualified lemma compute_float_times[code]: "Float m1 e1 * Float m2 e2 = Float (m1 * m2) (e1 + e2)"
603 qualified lemma compute_float_plus[code]:
604   "Float m1 e1 + Float m2 e2 =
605     (if m1 = 0 then Float m2 e2
606      else if m2 = 0 then Float m1 e1
607      else if e1 \<le> e2 then Float (m1 + m2 * 2^nat (e2 - e1)) e1
608      else Float (m2 + m1 * 2^nat (e1 - e2)) e2)"
609   by transfer (simp add: field_simps powr_realpow[symmetric] powr_diff)
611 qualified lemma compute_float_minus[code]: "f - g = f + (-g)" for f g :: float
612   by simp
614 qualified lemma compute_float_sgn[code]:
615   "sgn (Float m1 e1) = (if 0 < m1 then 1 else if m1 < 0 then -1 else 0)"
616   by transfer (simp add: sgn_mult)
618 lift_definition is_float_pos :: "float \<Rightarrow> bool" is "(<) 0 :: real \<Rightarrow> bool" .
620 qualified lemma compute_is_float_pos[code]: "is_float_pos (Float m e) \<longleftrightarrow> 0 < m"
621   by transfer (auto simp add: zero_less_mult_iff not_le[symmetric, of _ 0])
623 lift_definition is_float_nonneg :: "float \<Rightarrow> bool" is "(\<le>) 0 :: real \<Rightarrow> bool" .
625 qualified lemma compute_is_float_nonneg[code]: "is_float_nonneg (Float m e) \<longleftrightarrow> 0 \<le> m"
626   by transfer (auto simp add: zero_le_mult_iff not_less[symmetric, of _ 0])
628 lift_definition is_float_zero :: "float \<Rightarrow> bool"  is "(=) 0 :: real \<Rightarrow> bool" .
630 qualified lemma compute_is_float_zero[code]: "is_float_zero (Float m e) \<longleftrightarrow> 0 = m"
631   by transfer (auto simp add: is_float_zero_def)
633 qualified lemma compute_float_abs[code]: "\<bar>Float m e\<bar> = Float \<bar>m\<bar> e"
634   by transfer (simp add: abs_mult)
636 qualified lemma compute_float_eq[code]: "equal_class.equal f g = is_float_zero (f - g)"
637   by transfer simp
639 end
642 subsection \<open>Lemmas for types @{typ real}, @{typ nat}, @{typ int}\<close>
644 lemmas real_of_ints =
646   of_int_minus
647   of_int_diff
648   of_int_mult
649   of_int_power
650   of_int_numeral of_int_neg_numeral
652 lemmas int_of_reals = real_of_ints[symmetric]
655 subsection \<open>Rounding Real Numbers\<close>
657 definition round_down :: "int \<Rightarrow> real \<Rightarrow> real"
658   where "round_down prec x = \<lfloor>x * 2 powr prec\<rfloor> * 2 powr -prec"
660 definition round_up :: "int \<Rightarrow> real \<Rightarrow> real"
661   where "round_up prec x = \<lceil>x * 2 powr prec\<rceil> * 2 powr -prec"
663 lemma round_down_float[simp]: "round_down prec x \<in> float"
664   unfolding round_down_def
665   by (auto intro!: times_float simp: of_int_minus[symmetric] simp del: of_int_minus)
667 lemma round_up_float[simp]: "round_up prec x \<in> float"
668   unfolding round_up_def
669   by (auto intro!: times_float simp: of_int_minus[symmetric] simp del: of_int_minus)
671 lemma round_up: "x \<le> round_up prec x"
672   by (simp add: powr_minus_divide le_divide_eq round_up_def ceiling_correct)
674 lemma round_down: "round_down prec x \<le> x"
675   by (simp add: powr_minus_divide divide_le_eq round_down_def)
677 lemma round_up_0[simp]: "round_up p 0 = 0"
678   unfolding round_up_def by simp
680 lemma round_down_0[simp]: "round_down p 0 = 0"
681   unfolding round_down_def by simp
683 lemma round_up_diff_round_down: "round_up prec x - round_down prec x \<le> 2 powr -prec"
684 proof -
685   have "round_up prec x - round_down prec x = (\<lceil>x * 2 powr prec\<rceil> - \<lfloor>x * 2 powr prec\<rfloor>) * 2 powr -prec"
686     by (simp add: round_up_def round_down_def field_simps)
687   also have "\<dots> \<le> 1 * 2 powr -prec"
688     by (rule mult_mono)
689       (auto simp del: of_int_diff simp: of_int_diff[symmetric] ceiling_diff_floor_le_1)
690   finally show ?thesis by simp
691 qed
693 lemma round_down_shift: "round_down p (x * 2 powr k) = 2 powr k * round_down (p + k) x"
694   unfolding round_down_def
698 lemma round_up_shift: "round_up p (x * 2 powr k) = 2 powr k * round_up (p + k) x"
699   unfolding round_up_def
703 lemma round_up_uminus_eq: "round_up p (-x) = - round_down p x"
704   and round_down_uminus_eq: "round_down p (-x) = - round_up p x"
705   by (auto simp: round_up_def round_down_def ceiling_def)
707 lemma round_up_mono: "x \<le> y \<Longrightarrow> round_up p x \<le> round_up p y"
708   by (auto intro!: ceiling_mono simp: round_up_def)
710 lemma round_up_le1:
711   assumes "x \<le> 1" "prec \<ge> 0"
712   shows "round_up prec x \<le> 1"
713 proof -
714   have "real_of_int \<lceil>x * 2 powr prec\<rceil> \<le> real_of_int \<lceil>2 powr real_of_int prec\<rceil>"
715     using assms by (auto intro!: ceiling_mono)
716   also have "\<dots> = 2 powr prec" using assms by (auto simp: powr_int intro!: exI[where x="2^nat prec"])
717   finally show ?thesis
719 qed
721 lemma round_up_less1:
722   assumes "x < 1 / 2" "p > 0"
723   shows "round_up p x < 1"
724 proof -
725   have "x * 2 powr p < 1 / 2 * 2 powr p"
726     using assms by simp
727   also have "\<dots> \<le> 2 powr p - 1" using \<open>p > 0\<close>
728     by (auto simp: powr_diff powr_int field_simps self_le_power)
729   finally show ?thesis using \<open>p > 0\<close>
730     by (simp add: round_up_def field_simps powr_minus powr_int ceiling_less_iff)
731 qed
733 lemma round_down_ge1:
734   assumes x: "x \<ge> 1"
735   assumes prec: "p \<ge> - log 2 x"
736   shows "1 \<le> round_down p x"
737 proof cases
738   assume nonneg: "0 \<le> p"
739   have "2 powr p = real_of_int \<lfloor>2 powr real_of_int p\<rfloor>"
740     using nonneg by (auto simp: powr_int)
741   also have "\<dots> \<le> real_of_int \<lfloor>x * 2 powr p\<rfloor>"
742     using assms by (auto intro!: floor_mono)
743   finally show ?thesis
745 next
746   assume neg: "\<not> 0 \<le> p"
747   have "x = 2 powr (log 2 x)"
748     using x by simp
749   also have "2 powr (log 2 x) \<ge> 2 powr - p"
750     using prec by auto
751   finally have x_le: "x \<ge> 2 powr -p" .
753   from neg have "2 powr real_of_int p \<le> 2 powr 0"
754     by (intro powr_mono) auto
755   also have "\<dots> \<le> \<lfloor>2 powr 0::real\<rfloor>" by simp
756   also have "\<dots> \<le> \<lfloor>x * 2 powr (real_of_int p)\<rfloor>"
757     unfolding of_int_le_iff
758     using x x_le by (intro floor_mono) (simp add: powr_minus_divide field_simps)
759   finally show ?thesis
760     using prec x
761     by (simp add: round_down_def powr_minus_divide pos_le_divide_eq)
762 qed
764 lemma round_up_le0: "x \<le> 0 \<Longrightarrow> round_up p x \<le> 0"
765   unfolding round_up_def
766   by (auto simp: field_simps mult_le_0_iff zero_le_mult_iff)
769 subsection \<open>Rounding Floats\<close>
771 definition div_twopow :: "int \<Rightarrow> nat \<Rightarrow> int"
772   where [simp]: "div_twopow x n = x div (2 ^ n)"
774 definition mod_twopow :: "int \<Rightarrow> nat \<Rightarrow> int"
775   where [simp]: "mod_twopow x n = x mod (2 ^ n)"
777 lemma compute_div_twopow[code]:
778   "div_twopow x n = (if x = 0 \<or> x = -1 \<or> n = 0 then x else div_twopow (x div 2) (n - 1))"
779   by (cases n) (auto simp: zdiv_zmult2_eq div_eq_minus1)
781 lemma compute_mod_twopow[code]:
782   "mod_twopow x n = (if n = 0 then 0 else x mod 2 + 2 * mod_twopow (x div 2) (n - 1))"
783   by (cases n) (auto simp: zmod_zmult2_eq)
785 lift_definition float_up :: "int \<Rightarrow> float \<Rightarrow> float" is round_up by simp
786 declare float_up.rep_eq[simp]
788 lemma round_up_correct: "round_up e f - f \<in> {0..2 powr -e}"
789   unfolding atLeastAtMost_iff
790 proof
791   have "round_up e f - f \<le> round_up e f - round_down e f"
792     using round_down by simp
793   also have "\<dots> \<le> 2 powr -e"
794     using round_up_diff_round_down by simp
795   finally show "round_up e f - f \<le> 2 powr - (real_of_int e)"
796     by simp
797 qed (simp add: algebra_simps round_up)
799 lemma float_up_correct: "real_of_float (float_up e f) - real_of_float f \<in> {0..2 powr -e}"
800   by transfer (rule round_up_correct)
802 lift_definition float_down :: "int \<Rightarrow> float \<Rightarrow> float" is round_down by simp
803 declare float_down.rep_eq[simp]
805 lemma round_down_correct: "f - (round_down e f) \<in> {0..2 powr -e}"
806   unfolding atLeastAtMost_iff
807 proof
808   have "f - round_down e f \<le> round_up e f - round_down e f"
809     using round_up by simp
810   also have "\<dots> \<le> 2 powr -e"
811     using round_up_diff_round_down by simp
812   finally show "f - round_down e f \<le> 2 powr - (real_of_int e)"
813     by simp
814 qed (simp add: algebra_simps round_down)
816 lemma float_down_correct: "real_of_float f - real_of_float (float_down e f) \<in> {0..2 powr -e}"
817   by transfer (rule round_down_correct)
819 context
820 begin
822 qualified lemma compute_float_down[code]:
823   "float_down p (Float m e) =
824     (if p + e < 0 then Float (div_twopow m (nat (-(p + e)))) (-p) else Float m e)"
825 proof (cases "p + e < 0")
826   case True
827   then have "real_of_int ((2::int) ^ nat (-(p + e))) = 2 powr (-(p + e))"
828     using powr_realpow[of 2 "nat (-(p + e))"] by simp
829   also have "\<dots> = 1 / 2 powr p / 2 powr e"
831   finally show ?thesis
832     using \<open>p + e < 0\<close>
833     apply transfer
834     apply (simp add: ac_simps round_down_def floor_divide_of_int_eq[symmetric])
835     proof - (*FIXME*)
836       fix pa :: int and ea :: int and ma :: int
837       assume a1: "2 ^ nat (- pa - ea) = 1 / (2 powr real_of_int pa * 2 powr real_of_int ea)"
838       assume "pa + ea < 0"
839       have "\<lfloor>real_of_int ma / real_of_int (int 2 ^ nat (- (pa + ea)))\<rfloor> =
840           \<lfloor>real_of_float (Float ma (pa + ea))\<rfloor>"
842       then show "\<lfloor>real_of_int ma * (2 powr real_of_int pa * 2 powr real_of_int ea)\<rfloor> =
843           ma div 2 ^ nat (- pa - ea)"
844         by (metis Float.rep_eq add_uminus_conv_diff floor_divide_of_int_eq
846     qed
847 next
848   case False
849   then have r: "real_of_int e + real_of_int p = real (nat (e + p))"
850     by simp
851   have r: "\<lfloor>(m * 2 powr e) * 2 powr real_of_int p\<rfloor> = (m * 2 powr e) * 2 powr real_of_int p"
852     by (auto intro: exI[where x="m*2^nat (e+p)"]
854   with \<open>\<not> p + e < 0\<close> show ?thesis
856 qed
858 lemma abs_round_down_le: "\<bar>f - (round_down e f)\<bar> \<le> 2 powr -e"
859   using round_down_correct[of f e] by simp
861 lemma abs_round_up_le: "\<bar>f - (round_up e f)\<bar> \<le> 2 powr -e"
862   using round_up_correct[of e f] by simp
864 lemma round_down_nonneg: "0 \<le> s \<Longrightarrow> 0 \<le> round_down p s"
865   by (auto simp: round_down_def)
867 lemma ceil_divide_floor_conv:
868   assumes "b \<noteq> 0"
869   shows "\<lceil>real_of_int a / real_of_int b\<rceil> =
870     (if b dvd a then a div b else \<lfloor>real_of_int a / real_of_int b\<rfloor> + 1)"
871 proof (cases "b dvd a")
872   case True
873   then show ?thesis
874     by (simp add: ceiling_def of_int_minus[symmetric] divide_minus_left[symmetric]
875       floor_divide_of_int_eq dvd_neg_div del: divide_minus_left of_int_minus)
876 next
877   case False
878   then have "a mod b \<noteq> 0"
879     by auto
880   then have ne: "real_of_int (a mod b) / real_of_int b \<noteq> 0"
881     using \<open>b \<noteq> 0\<close> by auto
882   have "\<lceil>real_of_int a / real_of_int b\<rceil> = \<lfloor>real_of_int a / real_of_int b\<rfloor> + 1"
883     apply (rule ceiling_eq)
884     apply (auto simp: floor_divide_of_int_eq[symmetric])
885   proof -
886     have "real_of_int \<lfloor>real_of_int a / real_of_int b\<rfloor> \<le> real_of_int a / real_of_int b"
887       by simp
888     moreover have "real_of_int \<lfloor>real_of_int a / real_of_int b\<rfloor> \<noteq> real_of_int a / real_of_int b"
889       apply (subst (2) real_of_int_div_aux)
890       unfolding floor_divide_of_int_eq
891       using ne \<open>b \<noteq> 0\<close> apply auto
892       done
893     ultimately show "real_of_int \<lfloor>real_of_int a / real_of_int b\<rfloor> < real_of_int a / real_of_int b" by arith
894   qed
895   then show ?thesis
896     using \<open>\<not> b dvd a\<close> by simp
897 qed
899 qualified lemma compute_float_up[code]: "float_up p x = - float_down p (-x)"
900   by transfer (simp add: round_down_uminus_eq)
902 end
905 lemma bitlen_Float:
906   fixes m e
907   defines [THEN meta_eq_to_obj_eq]: "f \<equiv> Float m e"
908   shows "bitlen \<bar>mantissa f\<bar> + exponent f = (if m = 0 then 0 else bitlen \<bar>m\<bar> + e)"
909 proof (cases "m = 0")
910   case True
911   then show ?thesis by (simp add: f_def bitlen_alt_def)
912 next
913   case False
914   then have "f \<noteq> 0"
915     unfolding real_of_float_eq by (simp add: f_def)
916   then have "mantissa f \<noteq> 0"
918   moreover
919   obtain i where "m = mantissa f * 2 ^ i" "e = exponent f - int i"
920     by (rule f_def[THEN denormalize_shift, OF \<open>f \<noteq> 0\<close>])
921   ultimately show ?thesis by (simp add: abs_mult)
922 qed
924 lemma float_gt1_scale:
925   assumes "1 \<le> Float m e"
926   shows "0 \<le> e + (bitlen m - 1)"
927 proof -
928   have "0 < Float m e" using assms by auto
929   then have "0 < m" using powr_gt_zero[of 2 e]
930     by (auto simp: zero_less_mult_iff)
931   then have "m \<noteq> 0" by auto
932   show ?thesis
933   proof (cases "0 \<le> e")
934     case True
935     then show ?thesis
936       using \<open>0 < m\<close> by (simp add: bitlen_alt_def)
937   next
938     case False
939     have "(1::int) < 2" by simp
940     let ?S = "2^(nat (-e))"
941     have "inverse (2 ^ nat (- e)) = 2 powr e"
942       using assms False powr_realpow[of 2 "nat (-e)"]
943       by (auto simp: powr_minus field_simps)
944     then have "1 \<le> real_of_int m * inverse ?S"
945       using assms False powr_realpow[of 2 "nat (-e)"]
946       by (auto simp: powr_minus)
947     then have "1 * ?S \<le> real_of_int m * inverse ?S * ?S"
948       by (rule mult_right_mono) auto
949     then have "?S \<le> real_of_int m"
950       unfolding mult.assoc by auto
951     then have "?S \<le> m"
952       unfolding of_int_le_iff[symmetric] by auto
953     from this bitlen_bounds[OF \<open>0 < m\<close>, THEN conjunct2]
954     have "nat (-e) < (nat (bitlen m))"
955       unfolding power_strict_increasing_iff[OF \<open>1 < 2\<close>, symmetric]
956       by (rule order_le_less_trans)
957     then have "-e < bitlen m"
958       using False by auto
959     then show ?thesis
960       by auto
961   qed
962 qed
965 subsection \<open>Truncating Real Numbers\<close>
967 definition truncate_down::"nat \<Rightarrow> real \<Rightarrow> real"
968   where "truncate_down prec x = round_down (prec - \<lfloor>log 2 \<bar>x\<bar>\<rfloor>) x"
970 lemma truncate_down: "truncate_down prec x \<le> x"
971   using round_down by (simp add: truncate_down_def)
973 lemma truncate_down_le: "x \<le> y \<Longrightarrow> truncate_down prec x \<le> y"
974   by (rule order_trans[OF truncate_down])
976 lemma truncate_down_zero[simp]: "truncate_down prec 0 = 0"
979 lemma truncate_down_float[simp]: "truncate_down p x \<in> float"
980   by (auto simp: truncate_down_def)
982 definition truncate_up::"nat \<Rightarrow> real \<Rightarrow> real"
983   where "truncate_up prec x = round_up (prec - \<lfloor>log 2 \<bar>x\<bar>\<rfloor>) x"
985 lemma truncate_up: "x \<le> truncate_up prec x"
986   using round_up by (simp add: truncate_up_def)
988 lemma truncate_up_le: "x \<le> y \<Longrightarrow> x \<le> truncate_up prec y"
989   by (rule order_trans[OF _ truncate_up])
991 lemma truncate_up_zero[simp]: "truncate_up prec 0 = 0"
994 lemma truncate_up_uminus_eq: "truncate_up prec (-x) = - truncate_down prec x"
995   and truncate_down_uminus_eq: "truncate_down prec (-x) = - truncate_up prec x"
996   by (auto simp: truncate_up_def round_up_def truncate_down_def round_down_def ceiling_def)
998 lemma truncate_up_float[simp]: "truncate_up p x \<in> float"
999   by (auto simp: truncate_up_def)
1001 lemma mult_powr_eq: "0 < b \<Longrightarrow> b \<noteq> 1 \<Longrightarrow> 0 < x \<Longrightarrow> x * b powr y = b powr (y + log b x)"
1004 lemma truncate_down_pos:
1005   assumes "x > 0"
1006   shows "truncate_down p x > 0"
1007 proof -
1008   have "0 \<le> log 2 x - real_of_int \<lfloor>log 2 x\<rfloor>"
1010   with assms
1011   show ?thesis
1012     apply (auto simp: truncate_down_def round_down_def mult_powr_eq
1013       intro!: ge_one_powr_ge_zero mult_pos_pos)
1014     by linarith
1015 qed
1017 lemma truncate_down_nonneg: "0 \<le> y \<Longrightarrow> 0 \<le> truncate_down prec y"
1018   by (auto simp: truncate_down_def round_down_def)
1020 lemma truncate_down_ge1: "1 \<le> x \<Longrightarrow> 1 \<le> truncate_down p x"
1021   apply (auto simp: truncate_down_def algebra_simps intro!: round_down_ge1)
1022   apply linarith
1023   done
1025 lemma truncate_up_nonpos: "x \<le> 0 \<Longrightarrow> truncate_up prec x \<le> 0"
1026   by (auto simp: truncate_up_def round_up_def intro!: mult_nonpos_nonneg)
1028 lemma truncate_up_le1:
1029   assumes "x \<le> 1"
1030   shows "truncate_up p x \<le> 1"
1031 proof -
1032   consider "x \<le> 0" | "x > 0"
1033     by arith
1034   then show ?thesis
1035   proof cases
1036     case 1
1037     with truncate_up_nonpos[OF this, of p] show ?thesis
1038       by simp
1039   next
1040     case 2
1041     then have le: "\<lfloor>log 2 \<bar>x\<bar>\<rfloor> \<le> 0"
1042       using assms by (auto simp: log_less_iff)
1043     from assms have "0 \<le> int p" by simp
1045     show ?thesis
1047   qed
1048 qed
1050 lemma truncate_down_shift_int:
1051   "truncate_down p (x * 2 powr real_of_int k) = truncate_down p x * 2 powr k"
1052   by (cases "x = 0")
1053     (simp_all add: algebra_simps abs_mult log_mult truncate_down_def
1054       round_down_shift[of _ _ k, simplified])
1056 lemma truncate_down_shift_nat: "truncate_down p (x * 2 powr real k) = truncate_down p x * 2 powr k"
1057   by (metis of_int_of_nat_eq truncate_down_shift_int)
1059 lemma truncate_up_shift_int: "truncate_up p (x * 2 powr real_of_int k) = truncate_up p x * 2 powr k"
1060   by (cases "x = 0")
1061     (simp_all add: algebra_simps abs_mult log_mult truncate_up_def
1062       round_up_shift[of _ _ k, simplified])
1064 lemma truncate_up_shift_nat: "truncate_up p (x * 2 powr real k) = truncate_up p x * 2 powr k"
1065   by (metis of_int_of_nat_eq truncate_up_shift_int)
1068 subsection \<open>Truncating Floats\<close>
1070 lift_definition float_round_up :: "nat \<Rightarrow> float \<Rightarrow> float" is truncate_up
1073 lemma float_round_up: "real_of_float x \<le> real_of_float (float_round_up prec x)"
1074   using truncate_up by transfer simp
1076 lemma float_round_up_zero[simp]: "float_round_up prec 0 = 0"
1077   by transfer simp
1079 lift_definition float_round_down :: "nat \<Rightarrow> float \<Rightarrow> float" is truncate_down
1082 lemma float_round_down: "real_of_float (float_round_down prec x) \<le> real_of_float x"
1083   using truncate_down by transfer simp
1085 lemma float_round_down_zero[simp]: "float_round_down prec 0 = 0"
1086   by transfer simp
1088 lemmas float_round_up_le = order_trans[OF _ float_round_up]
1089   and float_round_down_le = order_trans[OF float_round_down]
1091 lemma minus_float_round_up_eq: "- float_round_up prec x = float_round_down prec (- x)"
1092   and minus_float_round_down_eq: "- float_round_down prec x = float_round_up prec (- x)"
1093   by (transfer; simp add: truncate_down_uminus_eq truncate_up_uminus_eq)+
1095 context
1096 begin
1098 qualified lemma compute_float_round_down[code]:
1099   "float_round_down prec (Float m e) =
1100     (let d = bitlen \<bar>m\<bar> - int prec - 1 in
1101       if 0 < d then Float (div_twopow m (nat d)) (e + d)
1102       else Float m e)"
1103   using Float.compute_float_down[of "Suc prec - bitlen \<bar>m\<bar> - e" m e, symmetric]
1104   by transfer
1105     (simp add: field_simps abs_mult log_mult bitlen_alt_def truncate_down_def
1106       cong del: if_weak_cong)
1108 qualified lemma compute_float_round_up[code]:
1109   "float_round_up prec x = - float_round_down prec (-x)"
1110   by transfer (simp add: truncate_down_uminus_eq)
1112 end
1115 subsection \<open>Approximation of positive rationals\<close>
1117 lemma div_mult_twopow_eq: "a div ((2::nat) ^ n) div b = a div (b * 2 ^ n)" for a b :: nat
1118   by (cases "b = 0") (simp_all add: div_mult2_eq[symmetric] ac_simps)
1120 lemma real_div_nat_eq_floor_of_divide: "a div b = real_of_int \<lfloor>a / b\<rfloor>" for a b :: nat
1121   by (simp add: floor_divide_of_nat_eq [of a b])
1123 definition "rat_precision prec x y =
1124   (let d = bitlen x - bitlen y
1125    in int prec - d + (if Float (abs x) 0 < Float (abs y) d then 1 else 0))"
1127 lemma floor_log_divide_eq:
1128   assumes "i > 0" "j > 0" "p > 1"
1129   shows "\<lfloor>log p (i / j)\<rfloor> = floor (log p i) - floor (log p j) -
1130     (if i \<ge> j * p powr (floor (log p i) - floor (log p j)) then 0 else 1)"
1131 proof -
1132   let ?l = "log p"
1133   let ?fl = "\<lambda>x. floor (?l x)"
1134   have "\<lfloor>?l (i / j)\<rfloor> = \<lfloor>?l i - ?l j\<rfloor>" using assms
1135     by (auto simp: log_divide)
1136   also have "\<dots> = floor (real_of_int (?fl i - ?fl j) + (?l i - ?fl i - (?l j - ?fl j)))"
1137     (is "_ = floor (_ + ?r)")
1140   also note \<open>p > 1\<close>
1141   note powr = powr_le_cancel_iff[symmetric, OF \<open>1 < p\<close>, THEN iffD2]
1142   note powr_strict = powr_less_cancel_iff[symmetric, OF \<open>1 < p\<close>, THEN iffD2]
1143   have "floor ?r = (if i \<ge> j * p powr (?fl i - ?fl j) then 0 else -1)" (is "_ = ?if")
1144     using assms
1145     by (linarith |
1146       auto
1147         intro!: floor_eq2
1148         intro: powr_strict powr
1149         simp: powr_diff powr_add divide_simps algebra_simps)+
1150   finally
1151   show ?thesis by simp
1152 qed
1154 lemma truncate_down_rat_precision:
1155     "truncate_down prec (real x / real y) = round_down (rat_precision prec x y) (real x / real y)"
1156   and truncate_up_rat_precision:
1157     "truncate_up prec (real x / real y) = round_up (rat_precision prec x y) (real x / real y)"
1158   unfolding truncate_down_def truncate_up_def rat_precision_def
1159   by (cases x; cases y) (auto simp: floor_log_divide_eq algebra_simps bitlen_alt_def)
1161 lift_definition lapprox_posrat :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float"
1162   is "\<lambda>prec (x::nat) (y::nat). truncate_down prec (x / y)"
1163   by simp
1165 context
1166 begin
1168 qualified lemma compute_lapprox_posrat[code]:
1169   "lapprox_posrat prec x y =
1170    (let
1171       l = rat_precision prec x y;
1172       d = if 0 \<le> l then x * 2^nat l div y else x div 2^nat (- l) div y
1173     in normfloat (Float d (- l)))"
1174     unfolding div_mult_twopow_eq
1175     by transfer
1176       (simp add: round_down_def powr_int real_div_nat_eq_floor_of_divide field_simps Let_def
1177         truncate_down_rat_precision del: two_powr_minus_int_float)
1179 end
1181 lift_definition rapprox_posrat :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float"
1182   is "\<lambda>prec (x::nat) (y::nat). truncate_up prec (x / y)"
1183   by simp
1185 context
1186 begin
1188 qualified lemma compute_rapprox_posrat[code]:
1189   fixes prec x y
1190   defines "l \<equiv> rat_precision prec x y"
1191   shows "rapprox_posrat prec x y =
1192    (let
1193       l = l;
1194       (r, s) = if 0 \<le> l then (x * 2^nat l, y) else (x, y * 2^nat(-l));
1195       d = r div s;
1196       m = r mod s
1197     in normfloat (Float (d + (if m = 0 \<or> y = 0 then 0 else 1)) (- l)))"
1198 proof (cases "y = 0")
1199   assume "y = 0"
1200   then show ?thesis by transfer simp
1201 next
1202   assume "y \<noteq> 0"
1203   show ?thesis
1204   proof (cases "0 \<le> l")
1205     case True
1206     define x' where "x' = x * 2 ^ nat l"
1207     have "int x * 2 ^ nat l = x'"
1209     moreover have "real x * 2 powr l = real x'"
1210       by (simp add: powr_realpow[symmetric] \<open>0 \<le> l\<close> x'_def)
1211     ultimately show ?thesis
1212       using ceil_divide_floor_conv[of y x'] powr_realpow[of 2 "nat l"] \<open>0 \<le> l\<close> \<open>y \<noteq> 0\<close>
1213         l_def[symmetric, THEN meta_eq_to_obj_eq]
1214       apply transfer
1215       apply (auto simp add: round_up_def truncate_up_rat_precision)
1216       apply (metis dvd_triv_left of_nat_dvd_iff)
1217       apply (metis floor_divide_of_int_eq of_int_of_nat_eq)
1218       done
1219    next
1220     case False
1221     define y' where "y' = y * 2 ^ nat (- l)"
1222     from \<open>y \<noteq> 0\<close> have "y' \<noteq> 0" by (simp add: y'_def)
1223     have "int y * 2 ^ nat (- l) = y'"
1225     moreover have "real x * real_of_int (2::int) powr real_of_int l / real y = x / real y'"
1226       using \<open>\<not> 0 \<le> l\<close> by (simp add: powr_realpow[symmetric] powr_minus y'_def field_simps)
1227     ultimately show ?thesis
1228       using ceil_divide_floor_conv[of y' x] \<open>\<not> 0 \<le> l\<close> \<open>y' \<noteq> 0\<close> \<open>y \<noteq> 0\<close>
1229         l_def[symmetric, THEN meta_eq_to_obj_eq]
1230       apply transfer
1231       apply (auto simp add: round_up_def ceil_divide_floor_conv truncate_up_rat_precision)
1232       apply (metis dvd_triv_left of_nat_dvd_iff)
1233       apply (metis floor_divide_of_int_eq of_int_of_nat_eq)
1234       done
1235   qed
1236 qed
1238 end
1240 lemma rat_precision_pos:
1241   assumes "0 \<le> x"
1242     and "0 < y"
1243     and "2 * x < y"
1244   shows "rat_precision n (int x) (int y) > 0"
1245 proof -
1246   have "0 < x \<Longrightarrow> log 2 x + 1 = log 2 (2 * x)"
1248   then have "bitlen (int x) < bitlen (int y)"
1249     using assms
1252   then show ?thesis
1253     using assms
1255 qed
1257 lemma rapprox_posrat_less1:
1258   "0 \<le> x \<Longrightarrow> 0 < y \<Longrightarrow> 2 * x < y \<Longrightarrow> real_of_float (rapprox_posrat n x y) < 1"
1259   by transfer (simp add: rat_precision_pos round_up_less1 truncate_up_rat_precision)
1261 lift_definition lapprox_rat :: "nat \<Rightarrow> int \<Rightarrow> int \<Rightarrow> float" is
1262   "\<lambda>prec (x::int) (y::int). truncate_down prec (x / y)"
1263   by simp
1265 context
1266 begin
1268 qualified lemma compute_lapprox_rat[code]:
1269   "lapprox_rat prec x y =
1270    (if y = 0 then 0
1271     else if 0 \<le> x then
1272      (if 0 < y then lapprox_posrat prec (nat x) (nat y)
1273       else - (rapprox_posrat prec (nat x) (nat (-y))))
1274       else
1275        (if 0 < y
1276         then - (rapprox_posrat prec (nat (-x)) (nat y))
1277         else lapprox_posrat prec (nat (-x)) (nat (-y))))"
1278   by transfer (simp add: truncate_up_uminus_eq)
1280 lift_definition rapprox_rat :: "nat \<Rightarrow> int \<Rightarrow> int \<Rightarrow> float" is
1281   "\<lambda>prec (x::int) (y::int). truncate_up prec (x / y)"
1282   by simp
1284 lemma "rapprox_rat = rapprox_posrat"
1285   by transfer auto
1287 lemma "lapprox_rat = lapprox_posrat"
1288   by transfer auto
1290 qualified lemma compute_rapprox_rat[code]:
1291   "rapprox_rat prec x y = - lapprox_rat prec (-x) y"
1292   by transfer (simp add: truncate_down_uminus_eq)
1294 qualified lemma compute_truncate_down[code]:
1295   "truncate_down p (Ratreal r) = (let (a, b) = quotient_of r in lapprox_rat p a b)"
1296   by transfer (auto split: prod.split simp: of_rat_divide dest!: quotient_of_div)
1298 qualified lemma compute_truncate_up[code]:
1299   "truncate_up p (Ratreal r) = (let (a, b) = quotient_of r in rapprox_rat p a b)"
1300   by transfer (auto split: prod.split simp:  of_rat_divide dest!: quotient_of_div)
1302 end
1305 subsection \<open>Division\<close>
1307 definition "real_divl prec a b = truncate_down prec (a / b)"
1309 definition "real_divr prec a b = truncate_up prec (a / b)"
1311 lift_definition float_divl :: "nat \<Rightarrow> float \<Rightarrow> float \<Rightarrow> float" is real_divl
1314 context
1315 begin
1317 qualified lemma compute_float_divl[code]:
1318   "float_divl prec (Float m1 s1) (Float m2 s2) = lapprox_rat prec m1 m2 * Float 1 (s1 - s2)"
1319   apply transfer
1320   unfolding real_divl_def of_int_1 mult_1 truncate_down_shift_int[symmetric]
1321   apply (simp add: powr_diff powr_minus)
1322   done
1324 lift_definition float_divr :: "nat \<Rightarrow> float \<Rightarrow> float \<Rightarrow> float" is real_divr
1327 qualified lemma compute_float_divr[code]:
1328   "float_divr prec x y = - float_divl prec (-x) y"
1329   by transfer (simp add: real_divr_def real_divl_def truncate_down_uminus_eq)
1331 end
1334 subsection \<open>Approximate Power\<close>
1336 lemma div2_less_self[termination_simp]: "odd n \<Longrightarrow> n div 2 < n" for n :: nat
1339 fun power_down :: "nat \<Rightarrow> real \<Rightarrow> nat \<Rightarrow> real"
1340 where
1341   "power_down p x 0 = 1"
1342 | "power_down p x (Suc n) =
1343     (if odd n then truncate_down (Suc p) ((power_down p x (Suc n div 2))\<^sup>2)
1344      else truncate_down (Suc p) (x * power_down p x n))"
1346 fun power_up :: "nat \<Rightarrow> real \<Rightarrow> nat \<Rightarrow> real"
1347 where
1348   "power_up p x 0 = 1"
1349 | "power_up p x (Suc n) =
1350     (if odd n then truncate_up p ((power_up p x (Suc n div 2))\<^sup>2)
1351      else truncate_up p (x * power_up p x n))"
1353 lift_definition power_up_fl :: "nat \<Rightarrow> float \<Rightarrow> nat \<Rightarrow> float" is power_up
1354   by (induct_tac rule: power_up.induct) simp_all
1356 lift_definition power_down_fl :: "nat \<Rightarrow> float \<Rightarrow> nat \<Rightarrow> float" is power_down
1357   by (induct_tac rule: power_down.induct) simp_all
1359 lemma power_float_transfer[transfer_rule]:
1360   "(rel_fun pcr_float (rel_fun (=) pcr_float)) (^) (^)"
1361   unfolding power_def
1362   by transfer_prover
1364 lemma compute_power_up_fl[code]:
1365     "power_up_fl p x 0 = 1"
1366     "power_up_fl p x (Suc n) =
1367       (if odd n then float_round_up p ((power_up_fl p x (Suc n div 2))\<^sup>2)
1368        else float_round_up p (x * power_up_fl p x n))"
1369   and compute_power_down_fl[code]:
1370     "power_down_fl p x 0 = 1"
1371     "power_down_fl p x (Suc n) =
1372       (if odd n then float_round_down (Suc p) ((power_down_fl p x (Suc n div 2))\<^sup>2)
1373        else float_round_down (Suc p) (x * power_down_fl p x n))"
1374   unfolding atomize_conj by transfer simp
1376 lemma power_down_pos: "0 < x \<Longrightarrow> 0 < power_down p x n"
1377   by (induct p x n rule: power_down.induct)
1378     (auto simp del: odd_Suc_div_two intro!: truncate_down_pos)
1380 lemma power_down_nonneg: "0 \<le> x \<Longrightarrow> 0 \<le> power_down p x n"
1381   by (induct p x n rule: power_down.induct)
1382     (auto simp del: odd_Suc_div_two intro!: truncate_down_nonneg mult_nonneg_nonneg)
1384 lemma power_down: "0 \<le> x \<Longrightarrow> power_down p x n \<le> x ^ n"
1385 proof (induct p x n rule: power_down.induct)
1386   case (2 p x n)
1387   have ?case if "odd n"
1388   proof -
1389     from that 2 have "(power_down p x (Suc n div 2)) ^ 2 \<le> (x ^ (Suc n div 2)) ^ 2"
1390       by (auto intro: power_mono power_down_nonneg simp del: odd_Suc_div_two)
1391     also have "\<dots> = x ^ (Suc n div 2 * 2)"
1393     also have "Suc n div 2 * 2 = Suc n"
1394       using \<open>odd n\<close> by presburger
1395     finally show ?thesis
1396       using that by (auto intro!: truncate_down_le simp del: odd_Suc_div_two)
1397   qed
1398   then show ?case
1399     by (auto intro!: truncate_down_le mult_left_mono 2 mult_nonneg_nonneg power_down_nonneg)
1400 qed simp
1402 lemma power_up: "0 \<le> x \<Longrightarrow> x ^ n \<le> power_up p x n"
1403 proof (induct p x n rule: power_up.induct)
1404   case (2 p x n)
1405   have ?case if "odd n"
1406   proof -
1407     from that even_Suc have "Suc n = Suc n div 2 * 2"
1408       by presburger
1409     then have "x ^ Suc n \<le> (x ^ (Suc n div 2))\<^sup>2"
1411     also from that 2 have "\<dots> \<le> (power_up p x (Suc n div 2))\<^sup>2"
1412       by (auto intro: power_mono simp del: odd_Suc_div_two)
1413     finally show ?thesis
1414       using that by (auto intro!: truncate_up_le simp del: odd_Suc_div_two)
1415   qed
1416   then show ?case
1417     by (auto intro!: truncate_up_le mult_left_mono 2)
1418 qed simp
1420 lemmas power_up_le = order_trans[OF _ power_up]
1421   and power_up_less = less_le_trans[OF _ power_up]
1422   and power_down_le = order_trans[OF power_down]
1424 lemma power_down_fl: "0 \<le> x \<Longrightarrow> power_down_fl p x n \<le> x ^ n"
1425   by transfer (rule power_down)
1427 lemma power_up_fl: "0 \<le> x \<Longrightarrow> x ^ n \<le> power_up_fl p x n"
1428   by transfer (rule power_up)
1430 lemma real_power_up_fl: "real_of_float (power_up_fl p x n) = power_up p x n"
1431   by transfer simp
1433 lemma real_power_down_fl: "real_of_float (power_down_fl p x n) = power_down p x n"
1434   by transfer simp
1439 definition "plus_down prec x y = truncate_down prec (x + y)"
1441 definition "plus_up prec x y = truncate_up prec (x + y)"
1443 lemma float_plus_down_float[intro, simp]: "x \<in> float \<Longrightarrow> y \<in> float \<Longrightarrow> plus_down p x y \<in> float"
1446 lemma float_plus_up_float[intro, simp]: "x \<in> float \<Longrightarrow> y \<in> float \<Longrightarrow> plus_up p x y \<in> float"
1449 lift_definition float_plus_down :: "nat \<Rightarrow> float \<Rightarrow> float \<Rightarrow> float" is plus_down ..
1451 lift_definition float_plus_up :: "nat \<Rightarrow> float \<Rightarrow> float \<Rightarrow> float" is plus_up ..
1453 lemma plus_down: "plus_down prec x y \<le> x + y"
1454   and plus_up: "x + y \<le> plus_up prec x y"
1455   by (auto simp: plus_down_def truncate_down plus_up_def truncate_up)
1457 lemma float_plus_down: "real_of_float (float_plus_down prec x y) \<le> x + y"
1458   and float_plus_up: "x + y \<le> real_of_float (float_plus_up prec x y)"
1459   by (transfer; rule plus_down plus_up)+
1461 lemmas plus_down_le = order_trans[OF plus_down]
1462   and plus_up_le = order_trans[OF _ plus_up]
1463   and float_plus_down_le = order_trans[OF float_plus_down]
1464   and float_plus_up_le = order_trans[OF _ float_plus_up]
1466 lemma compute_plus_up[code]: "plus_up p x y = - plus_down p (-x) (-y)"
1467   using truncate_down_uminus_eq[of p "x + y"]
1468   by (auto simp: plus_down_def plus_up_def)
1470 lemma truncate_down_log2_eqI:
1471   assumes "\<lfloor>log 2 \<bar>x\<bar>\<rfloor> = \<lfloor>log 2 \<bar>y\<bar>\<rfloor>"
1472   assumes "\<lfloor>x * 2 powr (p - \<lfloor>log 2 \<bar>x\<bar>\<rfloor>)\<rfloor> = \<lfloor>y * 2 powr (p - \<lfloor>log 2 \<bar>x\<bar>\<rfloor>)\<rfloor>"
1473   shows "truncate_down p x = truncate_down p y"
1474   using assms by (auto simp: truncate_down_def round_down_def)
1476 lemma sum_neq_zeroI:
1477   "\<bar>a\<bar> \<ge> k \<Longrightarrow> \<bar>b\<bar> < k \<Longrightarrow> a + b \<noteq> 0"
1478   "\<bar>a\<bar> > k \<Longrightarrow> \<bar>b\<bar> \<le> k \<Longrightarrow> a + b \<noteq> 0"
1479   for a k :: real
1480   by auto
1482 lemma abs_real_le_2_powr_bitlen[simp]: "\<bar>real_of_int m2\<bar> < 2 powr real_of_int (bitlen \<bar>m2\<bar>)"
1483 proof (cases "m2 = 0")
1484   case True
1485   then show ?thesis by simp
1486 next
1487   case False
1488   then have "\<bar>m2\<bar> < 2 ^ nat (bitlen \<bar>m2\<bar>)"
1489     using bitlen_bounds[of "\<bar>m2\<bar>"]
1490     by (auto simp: powr_add bitlen_nonneg)
1491   then show ?thesis
1492     by (metis bitlen_nonneg powr_int of_int_abs of_int_less_numeral_power_cancel_iff
1493         zero_less_numeral)
1494 qed
1496 lemma floor_sum_times_2_powr_sgn_eq:
1497   fixes ai p q :: int
1498     and a b :: real
1499   assumes "a * 2 powr p = ai"
1500     and b_le_1: "\<bar>b * 2 powr (p + 1)\<bar> \<le> 1"
1501     and leqp: "q \<le> p"
1502   shows "\<lfloor>(a + b) * 2 powr q\<rfloor> = \<lfloor>(2 * ai + sgn b) * 2 powr (q - p - 1)\<rfloor>"
1503 proof -
1504   consider "b = 0" | "b > 0" | "b < 0" by arith
1505   then show ?thesis
1506   proof cases
1507     case 1
1508     then show ?thesis
1510   next
1511     case 2
1512     then have "b * 2 powr p < \<bar>b * 2 powr (p + 1)\<bar>"
1513       by simp
1514     also note b_le_1
1515     finally have b_less_1: "b * 2 powr real_of_int p < 1" .
1517     from b_less_1 \<open>b > 0\<close> have floor_eq: "\<lfloor>b * 2 powr real_of_int p\<rfloor> = 0" "\<lfloor>sgn b / 2\<rfloor> = 0"
1520     have "\<lfloor>(a + b) * 2 powr q\<rfloor> = \<lfloor>(a + b) * 2 powr p * 2 powr (q - p)\<rfloor>"
1522     also have "\<dots> = \<lfloor>(ai + b * 2 powr p) * 2 powr (q - p)\<rfloor>"
1523       by (simp add: assms algebra_simps)
1524     also have "\<dots> = \<lfloor>(ai + b * 2 powr p) / real_of_int ((2::int) ^ nat (p - q))\<rfloor>"
1525       using assms
1527     also have "\<dots> = \<lfloor>ai / real_of_int ((2::int) ^ nat (p - q))\<rfloor>"
1528       by (simp del: of_int_power add: floor_divide_real_eq_div floor_eq)
1529     finally have "\<lfloor>(a + b) * 2 powr real_of_int q\<rfloor> = \<lfloor>real_of_int ai / real_of_int ((2::int) ^ nat (p - q))\<rfloor>" .
1530     moreover
1531     have "\<lfloor>(2 * ai + (sgn b)) * 2 powr (real_of_int (q - p) - 1)\<rfloor> =
1532         \<lfloor>real_of_int ai / real_of_int ((2::int) ^ nat (p - q))\<rfloor>"
1533     proof -
1534       have "\<lfloor>(2 * ai + sgn b) * 2 powr (real_of_int (q - p) - 1)\<rfloor> = \<lfloor>(ai + sgn b / 2) * 2 powr (q - p)\<rfloor>"
1535         by (subst powr_diff) (simp add: field_simps)
1536       also have "\<dots> = \<lfloor>(ai + sgn b / 2) / real_of_int ((2::int) ^ nat (p - q))\<rfloor>"
1537         using leqp by (simp add: powr_realpow[symmetric] powr_diff)
1538       also have "\<dots> = \<lfloor>ai / real_of_int ((2::int) ^ nat (p - q))\<rfloor>"
1539         by (simp del: of_int_power add: floor_divide_real_eq_div floor_eq)
1540       finally show ?thesis .
1541     qed
1542     ultimately show ?thesis by simp
1543   next
1544     case 3
1545     then have floor_eq: "\<lfloor>b * 2 powr (real_of_int p + 1)\<rfloor> = -1"
1546       using b_le_1
1547       by (auto simp: floor_eq_iff algebra_simps pos_divide_le_eq[symmetric] abs_if divide_powr_uminus
1548         intro!: mult_neg_pos split: if_split_asm)
1549     have "\<lfloor>(a + b) * 2 powr q\<rfloor> = \<lfloor>(2*a + 2*b) * 2 powr p * 2 powr (q - p - 1)\<rfloor>"
1551     also have "\<dots> = \<lfloor>(2 * (a * 2 powr p) + 2 * b * 2 powr p) * 2 powr (q - p - 1)\<rfloor>"
1553     also have "\<dots> = \<lfloor>(2 * ai + b * 2 powr (p + 1)) / 2 powr (1 - q + p)\<rfloor>"
1554       using assms by (simp add: algebra_simps powr_mult_base divide_powr_uminus)
1555     also have "\<dots> = \<lfloor>(2 * ai + b * 2 powr (p + 1)) / real_of_int ((2::int) ^ nat (p - q + 1))\<rfloor>"
1556       using assms by (simp add: algebra_simps powr_realpow[symmetric])
1557     also have "\<dots> = \<lfloor>(2 * ai - 1) / real_of_int ((2::int) ^ nat (p - q + 1))\<rfloor>"
1558       using \<open>b < 0\<close> assms
1559       by (simp add: floor_divide_of_int_eq floor_eq floor_divide_real_eq_div
1560         del: of_int_mult of_int_power of_int_diff)
1561     also have "\<dots> = \<lfloor>(2 * ai - 1) * 2 powr (q - p - 1)\<rfloor>"
1562       using assms by (simp add: algebra_simps divide_powr_uminus powr_realpow[symmetric])
1563     finally show ?thesis
1564       using \<open>b < 0\<close> by simp
1565   qed
1566 qed
1569   fixes ai :: int
1570     and b :: real
1571   assumes "\<bar>b\<bar> \<le> 1/2"
1572     and "ai \<noteq> 0"
1573   shows "\<lfloor>log 2 \<bar>real_of_int ai + b\<bar>\<rfloor> = \<lfloor>log 2 \<bar>ai + sgn b / 2\<bar>\<rfloor>"
1574 proof (cases "b = 0")
1575   case True
1576   then show ?thesis by simp
1577 next
1578   case False
1579   define k where "k = \<lfloor>log 2 \<bar>ai\<bar>\<rfloor>"
1580   then have "\<lfloor>log 2 \<bar>ai\<bar>\<rfloor> = k"
1581     by simp
1582   then have k: "2 powr k \<le> \<bar>ai\<bar>" "\<bar>ai\<bar> < 2 powr (k + 1)"
1583     by (simp_all add: floor_log_eq_powr_iff \<open>ai \<noteq> 0\<close>)
1584   have "k \<ge> 0"
1585     using assms by (auto simp: k_def)
1586   define r where "r = \<bar>ai\<bar> - 2 ^ nat k"
1587   have r: "0 \<le> r" "r < 2 powr k"
1588     using \<open>k \<ge> 0\<close> k
1589     by (auto simp: r_def k_def algebra_simps powr_add abs_if powr_int)
1590   then have "r \<le> (2::int) ^ nat k - 1"
1591     using \<open>k \<ge> 0\<close> by (auto simp: powr_int)
1592   from this[simplified of_int_le_iff[symmetric]] \<open>0 \<le> k\<close>
1593   have r_le: "r \<le> 2 powr k - 1"
1594     by (auto simp: algebra_simps powr_int)
1597   have "\<bar>ai\<bar> = 2 powr k + r"
1598     using \<open>k \<ge> 0\<close> by (auto simp: k_def r_def powr_realpow[symmetric])
1600   have pos: "\<bar>b\<bar> < 1 \<Longrightarrow> 0 < 2 powr k + (r + b)" for b :: real
1601     using \<open>0 \<le> k\<close> \<open>ai \<noteq> 0\<close>
1602     by (auto simp add: r_def powr_realpow[symmetric] abs_if sgn_if algebra_simps
1603       split: if_split_asm)
1604   have less: "\<bar>sgn ai * b\<bar> < 1"
1605     and less': "\<bar>sgn (sgn ai * b) / 2\<bar> < 1"
1606     using \<open>\<bar>b\<bar> \<le> _\<close> by (auto simp: abs_if sgn_if split: if_split_asm)
1608   have floor_eq: "\<And>b::real. \<bar>b\<bar> \<le> 1 / 2 \<Longrightarrow>
1609       \<lfloor>log 2 (1 + (r + b) / 2 powr k)\<rfloor> = (if r = 0 \<and> b < 0 then -1 else 0)"
1610     using \<open>k \<ge> 0\<close> r r_le
1611     by (auto simp: floor_log_eq_powr_iff powr_minus_divide field_simps sgn_if)
1613   from \<open>real_of_int \<bar>ai\<bar> = _\<close> have "\<bar>ai + b\<bar> = 2 powr k + (r + sgn ai * b)"
1614     using \<open>\<bar>b\<bar> \<le> _\<close> \<open>0 \<le> k\<close> r
1615     by (auto simp add: sgn_if abs_if)
1616   also have "\<lfloor>log 2 \<dots>\<rfloor> = \<lfloor>log 2 (2 powr k + r + sgn (sgn ai * b) / 2)\<rfloor>"
1617   proof -
1618     have "2 powr k + (r + (sgn ai) * b) = 2 powr k * (1 + (r + sgn ai * b) / 2 powr k)"
1620     also have "\<lfloor>log 2 \<dots>\<rfloor> = k + \<lfloor>log 2 (1 + (r + sgn ai * b) / 2 powr k)\<rfloor>"
1621       using pos[OF less]
1622       by (subst log_mult) (simp_all add: log_mult powr_mult field_simps)
1623     also
1624     let ?if = "if r = 0 \<and> sgn ai * b < 0 then -1 else 0"
1625     have "\<lfloor>log 2 (1 + (r + sgn ai * b) / 2 powr k)\<rfloor> = ?if"
1626       using \<open>\<bar>b\<bar> \<le> _\<close>
1627       by (intro floor_eq) (auto simp: abs_mult sgn_if)
1628     also
1629     have "\<dots> = \<lfloor>log 2 (1 + (r + sgn (sgn ai * b) / 2) / 2 powr k)\<rfloor>"
1630       by (subst floor_eq) (auto simp: sgn_if)
1631     also have "k + \<dots> = \<lfloor>log 2 (2 powr k * (1 + (r + sgn (sgn ai * b) / 2) / 2 powr k))\<rfloor>"
1633       using pos[OF less'] \<open>\<bar>b\<bar> \<le> _\<close>
1635     also have "2 powr k * (1 + (r + sgn (sgn ai * b) / 2) / 2 powr k) =
1636         2 powr k + r + sgn (sgn ai * b) / 2"
1637       by (simp add: sgn_if field_simps)
1638     finally show ?thesis .
1639   qed
1640   also have "2 powr k + r + sgn (sgn ai * b) / 2 = \<bar>ai + sgn b / 2\<bar>"
1641     unfolding \<open>real_of_int \<bar>ai\<bar> = _\<close>[symmetric] using \<open>ai \<noteq> 0\<close>
1642     by (auto simp: abs_if sgn_if algebra_simps)
1643   finally show ?thesis .
1644 qed
1646 context
1647 begin
1649 qualified lemma compute_far_float_plus_down:
1650   fixes m1 e1 m2 e2 :: int
1651     and p :: nat
1652   defines "k1 \<equiv> Suc p - nat (bitlen \<bar>m1\<bar>)"
1653   assumes H: "bitlen \<bar>m2\<bar> \<le> e1 - e2 - k1 - 2" "m1 \<noteq> 0" "m2 \<noteq> 0" "e1 \<ge> e2"
1654   shows "float_plus_down p (Float m1 e1) (Float m2 e2) =
1655     float_round_down p (Float (m1 * 2 ^ (Suc (Suc k1)) + sgn m2) (e1 - int k1 - 2))"
1656 proof -
1657   let ?a = "real_of_float (Float m1 e1)"
1658   let ?b = "real_of_float (Float m2 e2)"
1659   let ?sum = "?a + ?b"
1660   let ?shift = "real_of_int e2 - real_of_int e1 + real k1 + 1"
1661   let ?m1 = "m1 * 2 ^ Suc k1"
1662   let ?m2 = "m2 * 2 powr ?shift"
1663   let ?m2' = "sgn m2 / 2"
1664   let ?e = "e1 - int k1 - 1"
1666   have sum_eq: "?sum = (?m1 + ?m2) * 2 powr ?e"
1667     by (auto simp: powr_add[symmetric] powr_mult[symmetric] algebra_simps
1668       powr_realpow[symmetric] powr_mult_base)
1670   have "\<bar>?m2\<bar> * 2 < 2 powr (bitlen \<bar>m2\<bar> + ?shift + 1)"
1671     by (auto simp: field_simps powr_add powr_mult_base powr_diff abs_mult)
1672   also have "\<dots> \<le> 2 powr 0"
1673     using H by (intro powr_mono) auto
1674   finally have abs_m2_less_half: "\<bar>?m2\<bar> < 1 / 2"
1675     by simp
1677   then have "\<bar>real_of_int m2\<bar> < 2 powr -(?shift + 1)"
1678     unfolding powr_minus_divide by (auto simp: bitlen_alt_def field_simps powr_mult_base abs_mult)
1679   also have "\<dots> \<le> 2 powr real_of_int (e1 - e2 - 2)"
1680     by simp
1681   finally have b_less_quarter: "\<bar>?b\<bar> < 1/4 * 2 powr real_of_int e1"
1683   also have "1/4 < \<bar>real_of_int m1\<bar> / 2" using \<open>m1 \<noteq> 0\<close> by simp
1684   finally have b_less_half_a: "\<bar>?b\<bar> < 1/2 * \<bar>?a\<bar>"
1685     by (simp add: algebra_simps powr_mult_base abs_mult)
1686   then have a_half_less_sum: "\<bar>?a\<bar> / 2 < \<bar>?sum\<bar>"
1687     by (auto simp: field_simps abs_if split: if_split_asm)
1689   from b_less_half_a have "\<bar>?b\<bar> < \<bar>?a\<bar>" "\<bar>?b\<bar> \<le> \<bar>?a\<bar>"
1690     by simp_all
1692   have "\<bar>real_of_float (Float m1 e1)\<bar> \<ge> 1/4 * 2 powr real_of_int e1"
1693     using \<open>m1 \<noteq> 0\<close>
1694     by (auto simp: powr_add powr_int bitlen_nonneg divide_right_mono abs_mult)
1695   then have "?sum \<noteq> 0" using b_less_quarter
1696     by (rule sum_neq_zeroI)
1697   then have "?m1 + ?m2 \<noteq> 0"
1698     unfolding sum_eq by (simp add: abs_mult zero_less_mult_iff)
1700   have "\<bar>real_of_int ?m1\<bar> \<ge> 2 ^ Suc k1" "\<bar>?m2'\<bar> < 2 ^ Suc k1"
1701     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)
1702   then have sum'_nz: "?m1 + ?m2' \<noteq> 0"
1703     by (intro sum_neq_zeroI)
1705   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"
1706     using \<open>?m1 + ?m2 \<noteq> 0\<close>
1708     by (simp add: abs_mult log_mult) linarith
1709   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>"
1710     using abs_m2_less_half \<open>m1 \<noteq> 0\<close>
1711     by (intro log2_abs_int_add_less_half_sgn_eq) (auto simp: abs_mult)
1712   also have "sgn (real_of_int m2 * 2 powr ?shift) = sgn m2"
1713     by (auto simp: sgn_if zero_less_mult_iff less_not_sym)
1714   also
1715   have "\<bar>?m1 + ?m2'\<bar> * 2 powr ?e = \<bar>?m1 * 2 + sgn m2\<bar> * 2 powr (?e - 1)"
1716     by (auto simp: field_simps powr_minus[symmetric] powr_diff powr_mult_base)
1717   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>"
1718     using \<open>?m1 + ?m2' \<noteq> 0\<close>
1721   finally
1722   have "\<lfloor>log 2 \<bar>?sum\<bar>\<rfloor> = \<lfloor>log 2 \<bar>real_of_float (Float (?m1*2 + sgn m2) (?e - 1))\<bar>\<rfloor>" .
1723   then have "plus_down p (Float m1 e1) (Float m2 e2) =
1724       truncate_down p (Float (?m1*2 + sgn m2) (?e - 1))"
1725     unfolding plus_down_def
1726   proof (rule truncate_down_log2_eqI)
1727     let ?f = "(int p - \<lfloor>log 2 \<bar>real_of_float (Float m1 e1) + real_of_float (Float m2 e2)\<bar>\<rfloor>)"
1728     let ?ai = "m1 * 2 ^ (Suc k1)"
1729     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>"
1730     proof (rule floor_sum_times_2_powr_sgn_eq)
1731       show "?a * 2 powr real_of_int (-?e) = real_of_int ?ai"
1733       show "\<bar>?b * 2 powr real_of_int (-?e + 1)\<bar> \<le> 1"
1734         using abs_m2_less_half
1736     next
1737       have "e1 + \<lfloor>log 2 \<bar>real_of_int m1\<bar>\<rfloor> - 1 = \<lfloor>log 2 \<bar>?a\<bar>\<rfloor> - 1"
1738         using \<open>m1 \<noteq> 0\<close>
1740       also have "\<dots> \<le> \<lfloor>log 2 \<bar>?a + ?b\<bar>\<rfloor>"
1741         using a_half_less_sum \<open>m1 \<noteq> 0\<close> \<open>?sum \<noteq> 0\<close>
1742         unfolding floor_diff_of_int[symmetric]
1743         by (auto simp add: log_minus_eq_powr powr_minus_divide intro!: floor_mono)
1744       finally
1745       have "int p - \<lfloor>log 2 \<bar>?a + ?b\<bar>\<rfloor> \<le> p - (bitlen \<bar>m1\<bar>) - e1 + 2"
1746         by (auto simp: algebra_simps bitlen_alt_def \<open>m1 \<noteq> 0\<close>)
1747       also have "\<dots> \<le> - ?e"
1748         using bitlen_nonneg[of "\<bar>m1\<bar>"] by (simp add: k1_def)
1749       finally show "?f \<le> - ?e" by simp
1750     qed
1751     also have "sgn ?b = sgn m2"
1752       using powr_gt_zero[of 2 e2]
1753       by (auto simp add: sgn_if zero_less_mult_iff simp del: powr_gt_zero)
1754     also have "\<lfloor>(real_of_int (2 * ?m1) + real_of_int (sgn m2)) * 2 powr real_of_int (?f - - ?e - 1)\<rfloor> =
1755         \<lfloor>Float (?m1 * 2 + sgn m2) (?e - 1) * 2 powr ?f\<rfloor>"
1757     finally
1758     show "\<lfloor>(?a + ?b) * 2 powr ?f\<rfloor> = \<lfloor>real_of_float (Float (?m1 * 2 + sgn m2) (?e - 1)) * 2 powr ?f\<rfloor>" .
1759   qed
1760   then show ?thesis
1761     by transfer (simp add: plus_down_def ac_simps Let_def)
1762 qed
1764 lemma compute_float_plus_down_naive[code]: "float_plus_down p x y = float_round_down p (x + y)"
1765   by transfer (auto simp: plus_down_def)
1767 qualified lemma compute_float_plus_down[code]:
1768   fixes p::nat and m1 e1 m2 e2::int
1769   shows "float_plus_down p (Float m1 e1) (Float m2 e2) =
1770     (if m1 = 0 then float_round_down p (Float m2 e2)
1771     else if m2 = 0 then float_round_down p (Float m1 e1)
1772     else
1773       (if e1 \<ge> e2 then
1774         (let k1 = Suc p - nat (bitlen \<bar>m1\<bar>) in
1775           if bitlen \<bar>m2\<bar> > e1 - e2 - k1 - 2
1776           then float_round_down p ((Float m1 e1) + (Float m2 e2))
1777           else float_round_down p (Float (m1 * 2 ^ (Suc (Suc k1)) + sgn m2) (e1 - int k1 - 2)))
1778     else float_plus_down p (Float m2 e2) (Float m1 e1)))"
1779 proof -
1780   {
1781     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"
1782     note compute_far_float_plus_down[OF this]
1783   }
1784   then show ?thesis
1785     by transfer (simp add: Let_def plus_down_def ac_simps)
1786 qed
1788 qualified lemma compute_float_plus_up[code]: "float_plus_up p x y = - float_plus_down p (-x) (-y)"
1789   using truncate_down_uminus_eq[of p "x + y"]
1790   by transfer (simp add: plus_down_def plus_up_def ac_simps)
1792 lemma mantissa_zero[simp]: "mantissa 0 = 0"
1793   by (metis mantissa_0 zero_float.abs_eq)
1795 qualified lemma compute_float_less[code]: "a < b \<longleftrightarrow> is_float_pos (float_plus_down 0 b (- a))"
1796   using truncate_down[of 0 "b - a"] truncate_down_pos[of "b - a" 0]
1797   by transfer (auto simp: plus_down_def)
1799 qualified lemma compute_float_le[code]: "a \<le> b \<longleftrightarrow> is_float_nonneg (float_plus_down 0 b (- a))"
1800   using truncate_down[of 0 "b - a"] truncate_down_nonneg[of "b - a" 0]
1801   by transfer (auto simp: plus_down_def)
1803 end
1806 subsection \<open>Lemmas needed by Approximate\<close>
1808 lemma Float_num[simp]:
1809    "real_of_float (Float 1 0) = 1"
1810    "real_of_float (Float 1 1) = 2"
1811    "real_of_float (Float 1 2) = 4"
1812    "real_of_float (Float 1 (- 1)) = 1/2"
1813    "real_of_float (Float 1 (- 2)) = 1/4"
1814    "real_of_float (Float 1 (- 3)) = 1/8"
1815    "real_of_float (Float (- 1) 0) = -1"
1816    "real_of_float (Float (numeral n) 0) = numeral n"
1817    "real_of_float (Float (- numeral n) 0) = - numeral n"
1818   using two_powr_int_float[of 2] two_powr_int_float[of "-1"] two_powr_int_float[of "-2"]
1819     two_powr_int_float[of "-3"]
1820   using powr_realpow[of 2 2] powr_realpow[of 2 3]
1821   using powr_minus[of "2::real" 1] powr_minus[of "2::real" 2] powr_minus[of "2::real" 3]
1822   by auto
1824 lemma real_of_Float_int[simp]: "real_of_float (Float n 0) = real n"
1825   by simp
1827 lemma float_zero[simp]: "real_of_float (Float 0 e) = 0"
1828   by simp
1830 lemma abs_div_2_less: "a \<noteq> 0 \<Longrightarrow> a \<noteq> -1 \<Longrightarrow> \<bar>(a::int) div 2\<bar> < \<bar>a\<bar>"
1831   by arith
1833 lemma lapprox_rat: "real_of_float (lapprox_rat prec x y) \<le> real_of_int x / real_of_int y"
1834   by (simp add: lapprox_rat.rep_eq truncate_down)
1836 lemma mult_div_le:
1837   fixes a b :: int
1838   assumes "b > 0"
1839   shows "a \<ge> b * (a div b)"
1840 proof -
1841   from minus_div_mult_eq_mod [symmetric, of a b] have "a = b * (a div b) + a mod b"
1842     by simp
1843   also have "\<dots> \<ge> b * (a div b) + 0"
1845     apply (rule pos_mod_sign)
1846     using assms
1847     apply simp
1848     done
1849   finally show ?thesis
1850     by simp
1851 qed
1853 lemma lapprox_rat_nonneg:
1854   assumes "0 \<le> x" and "0 \<le> y"
1855   shows "0 \<le> real_of_float (lapprox_rat n x y)"
1856   using assms
1857   by transfer (simp add: truncate_down_nonneg)
1859 lemma rapprox_rat: "real_of_int x / real_of_int y \<le> real_of_float (rapprox_rat prec x y)"
1860   by transfer (simp add: truncate_up)
1862 lemma rapprox_rat_le1:
1863   assumes "0 \<le> x" "0 < y" "x \<le> y"
1864   shows "real_of_float (rapprox_rat n x y) \<le> 1"
1865   using assms
1866   by transfer (simp add: truncate_up_le1)
1868 lemma rapprox_rat_nonneg_nonpos: "0 \<le> x \<Longrightarrow> y \<le> 0 \<Longrightarrow> real_of_float (rapprox_rat n x y) \<le> 0"
1869   by transfer (simp add: truncate_up_nonpos divide_nonneg_nonpos)
1871 lemma rapprox_rat_nonpos_nonneg: "x \<le> 0 \<Longrightarrow> 0 \<le> y \<Longrightarrow> real_of_float (rapprox_rat n x y) \<le> 0"
1872   by transfer (simp add: truncate_up_nonpos divide_nonpos_nonneg)
1874 lemma real_divl: "real_divl prec x y \<le> x / y"
1875   by (simp add: real_divl_def truncate_down)
1877 lemma real_divr: "x / y \<le> real_divr prec x y"
1878   by (simp add: real_divr_def truncate_up)
1880 lemma float_divl: "real_of_float (float_divl prec x y) \<le> x / y"
1881   by transfer (rule real_divl)
1883 lemma real_divl_lower_bound: "0 \<le> x \<Longrightarrow> 0 \<le> y \<Longrightarrow> 0 \<le> real_divl prec x y"
1884   by (simp add: real_divl_def truncate_down_nonneg)
1886 lemma float_divl_lower_bound: "0 \<le> x \<Longrightarrow> 0 \<le> y \<Longrightarrow> 0 \<le> real_of_float (float_divl prec x y)"
1887   by transfer (rule real_divl_lower_bound)
1889 lemma exponent_1: "exponent 1 = 0"
1890   using exponent_float[of 1 0] by (simp add: one_float_def)
1892 lemma mantissa_1: "mantissa 1 = 1"
1893   using mantissa_float[of 1 0] by (simp add: one_float_def)
1895 lemma bitlen_1: "bitlen 1 = 1"
1898 lemma float_upper_bound: "x \<le> 2 powr (bitlen \<bar>mantissa x\<bar> + exponent x)"
1899 proof (cases "x = 0")
1900   case True
1901   then show ?thesis by simp
1902 next
1903   case False
1904   then have "mantissa x \<noteq> 0"
1905     using mantissa_eq_zero_iff by auto
1906   have "x = mantissa x * 2 powr (exponent x)"
1907     by (rule mantissa_exponent)
1908   also have "mantissa x \<le> \<bar>mantissa x\<bar>"
1909     by simp
1910   also have "\<dots> \<le> 2 powr (bitlen \<bar>mantissa x\<bar>)"
1911     using bitlen_bounds[of "\<bar>mantissa x\<bar>"] bitlen_nonneg \<open>mantissa x \<noteq> 0\<close>
1912     by (auto simp del: of_int_abs simp add: powr_int)
1914 qed
1916 lemma real_divl_pos_less1_bound:
1917   assumes "0 < x" "x \<le> 1"
1918   shows "1 \<le> real_divl prec 1 x"
1919   using assms
1920   by (auto intro!: truncate_down_ge1 simp: real_divl_def)
1922 lemma float_divl_pos_less1_bound:
1923   "0 < real_of_float x \<Longrightarrow> real_of_float x \<le> 1 \<Longrightarrow> prec \<ge> 1 \<Longrightarrow>
1924     1 \<le> real_of_float (float_divl prec 1 x)"
1925   by transfer (rule real_divl_pos_less1_bound)
1927 lemma float_divr: "real_of_float x / real_of_float y \<le> real_of_float (float_divr prec x y)"
1928   by transfer (rule real_divr)
1930 lemma real_divr_pos_less1_lower_bound:
1931   assumes "0 < x"
1932     and "x \<le> 1"
1933   shows "1 \<le> real_divr prec 1 x"
1934 proof -
1935   have "1 \<le> 1 / x"
1936     using \<open>0 < x\<close> and \<open>x \<le> 1\<close> by auto
1937   also have "\<dots> \<le> real_divr prec 1 x"
1938     using real_divr[where x = 1 and y = x] by auto
1939   finally show ?thesis by auto
1940 qed
1942 lemma float_divr_pos_less1_lower_bound: "0 < x \<Longrightarrow> x \<le> 1 \<Longrightarrow> 1 \<le> float_divr prec 1 x"
1943   by transfer (rule real_divr_pos_less1_lower_bound)
1945 lemma real_divr_nonpos_pos_upper_bound: "x \<le> 0 \<Longrightarrow> 0 \<le> y \<Longrightarrow> real_divr prec x y \<le> 0"
1946   by (simp add: real_divr_def truncate_up_nonpos divide_le_0_iff)
1948 lemma float_divr_nonpos_pos_upper_bound:
1949   "real_of_float x \<le> 0 \<Longrightarrow> 0 \<le> real_of_float y \<Longrightarrow> real_of_float (float_divr prec x y) \<le> 0"
1950   by transfer (rule real_divr_nonpos_pos_upper_bound)
1952 lemma real_divr_nonneg_neg_upper_bound: "0 \<le> x \<Longrightarrow> y \<le> 0 \<Longrightarrow> real_divr prec x y \<le> 0"
1953   by (simp add: real_divr_def truncate_up_nonpos divide_le_0_iff)
1955 lemma float_divr_nonneg_neg_upper_bound:
1956   "0 \<le> real_of_float x \<Longrightarrow> real_of_float y \<le> 0 \<Longrightarrow> real_of_float (float_divr prec x y) \<le> 0"
1957   by transfer (rule real_divr_nonneg_neg_upper_bound)
1959 lemma truncate_up_nonneg_mono:
1960   assumes "0 \<le> x" "x \<le> y"
1961   shows "truncate_up prec x \<le> truncate_up prec y"
1962 proof -
1963   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"
1964     by arith
1965   then show ?thesis
1966   proof cases
1967     case 1
1968     then show ?thesis
1969       using assms
1970       by (auto simp: truncate_up_def round_up_def intro!: ceiling_mono)
1971   next
1972     case 2
1973     from assms \<open>0 < x\<close> have "log 2 x \<le> log 2 y"
1974       by auto
1975     with \<open>\<lfloor>log 2 x\<rfloor> \<noteq> \<lfloor>log 2 y\<rfloor>\<close>
1976     have logless: "log 2 x < log 2 y" and flogless: "\<lfloor>log 2 x\<rfloor> < \<lfloor>log 2 y\<rfloor>"
1977       by (metis floor_less_cancel linorder_cases not_le)+
1978     have "truncate_up prec x =
1979       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>)"
1980       using assms by (simp add: truncate_up_def round_up_def)
1981     also have "\<lceil>x * 2 powr real_of_int (int prec - \<lfloor>log 2 x\<rfloor>)\<rceil> \<le> (2 ^ (Suc prec))"
1982     proof (simp only: ceiling_le_iff)
1983       have "x * 2 powr real_of_int (int prec - \<lfloor>log 2 x\<rfloor>) \<le>
1984         x * (2 powr real (Suc prec) / (2 powr log 2 x))"
1985         using real_of_int_floor_add_one_ge[of "log 2 x"] assms
1986         by (auto simp add: algebra_simps powr_diff [symmetric]  intro!: mult_left_mono)
1987       then show "x * 2 powr real_of_int (int prec - \<lfloor>log 2 x\<rfloor>) \<le> real_of_int ((2::int) ^ (Suc prec))"
1989     qed
1990     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)"
1991       by (auto simp: powr_realpow powr_add)
1992         (metis power_Suc of_int_le_numeral_power_cancel_iff)
1993     also
1994     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)"
1995       using logless flogless by (auto intro!: floor_mono)
1996     also have "2 powr real_of_int (int (Suc prec)) \<le>
1997         2 powr (log 2 y + real_of_int (int prec - \<lfloor>log 2 y\<rfloor> + 1))"
1998       using assms \<open>0 < x\<close>
1999       by (auto simp: algebra_simps)
2000     finally have "truncate_up prec x \<le>
2001         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)"
2002       by simp
2003     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>))"
2005     also have "\<dots> = y"
2006       using \<open>0 < x\<close> assms
2008     also have "\<dots> \<le> truncate_up prec y"
2009       by (rule truncate_up)
2010     finally show ?thesis .
2011   next
2012     case 3
2013     then show ?thesis
2014       using assms
2015       by (auto intro!: truncate_up_le)
2016   qed
2017 qed
2019 lemma truncate_up_switch_sign_mono:
2020   assumes "x \<le> 0" "0 \<le> y"
2021   shows "truncate_up prec x \<le> truncate_up prec y"
2022 proof -
2023   note truncate_up_nonpos[OF \<open>x \<le> 0\<close>]
2024   also note truncate_up_le[OF \<open>0 \<le> y\<close>]
2025   finally show ?thesis .
2026 qed
2028 lemma truncate_down_switch_sign_mono:
2029   assumes "x \<le> 0"
2030     and "0 \<le> y"
2031     and "x \<le> y"
2032   shows "truncate_down prec x \<le> truncate_down prec y"
2033 proof -
2034   note truncate_down_le[OF \<open>x \<le> 0\<close>]
2035   also note truncate_down_nonneg[OF \<open>0 \<le> y\<close>]
2036   finally show ?thesis .
2037 qed
2039 lemma truncate_down_nonneg_mono:
2040   assumes "0 \<le> x" "x \<le> y"
2041   shows "truncate_down prec x \<le> truncate_down prec y"
2042 proof -
2043   consider "x \<le> 0" | "\<lfloor>log 2 \<bar>x\<bar>\<rfloor> = \<lfloor>log 2 \<bar>y\<bar>\<rfloor>" |
2044     "0 < x" "\<lfloor>log 2 \<bar>x\<bar>\<rfloor> \<noteq> \<lfloor>log 2 \<bar>y\<bar>\<rfloor>"
2045     by arith
2046   then show ?thesis
2047   proof cases
2048     case 1
2049     with assms have "x = 0" "0 \<le> y" by simp_all
2050     then show ?thesis
2051       by (auto intro!: truncate_down_nonneg)
2052   next
2053     case 2
2054     then show ?thesis
2055       using assms
2056       by (auto simp: truncate_down_def round_down_def intro!: floor_mono)
2057   next
2058     case 3
2059     from \<open>0 < x\<close> have "log 2 x \<le> log 2 y" "0 < y" "0 \<le> y"
2060       using assms by auto
2061     with \<open>\<lfloor>log 2 \<bar>x\<bar>\<rfloor> \<noteq> \<lfloor>log 2 \<bar>y\<bar>\<rfloor>\<close>
2062     have logless: "log 2 x < log 2 y" and flogless: "\<lfloor>log 2 x\<rfloor> < \<lfloor>log 2 y\<rfloor>"
2063       unfolding atomize_conj abs_of_pos[OF \<open>0 < x\<close>] abs_of_pos[OF \<open>0 < y\<close>]
2064       by (metis floor_less_cancel linorder_cases not_le)
2065     have "2 powr prec \<le> y * 2 powr real prec / (2 powr log 2 y)"
2066       using \<open>0 < y\<close> by simp
2067     also have "\<dots> \<le> y * 2 powr real (Suc prec) / (2 powr (real_of_int \<lfloor>log 2 y\<rfloor> + 1))"
2068       using \<open>0 \<le> y\<close> \<open>0 \<le> x\<close> assms(2)
2069       by (auto intro!: powr_mono divide_left_mono
2071     also have "\<dots> = y * 2 powr real (Suc prec) / (2 powr real_of_int \<lfloor>log 2 y\<rfloor> * 2)"
2073     finally have "(2 ^ prec) \<le> \<lfloor>y * 2 powr real_of_int (int (Suc prec) - \<lfloor>log 2 \<bar>y\<bar>\<rfloor> - 1)\<rfloor>"
2074       using \<open>0 \<le> y\<close>
2075       by (auto simp: powr_diff le_floor_iff powr_realpow powr_add)
2076     then have "(2 ^ (prec)) * 2 powr - real_of_int (int prec - \<lfloor>log 2 \<bar>y\<bar>\<rfloor>) \<le> truncate_down prec y"
2077       by (auto simp: truncate_down_def round_down_def)
2078     moreover have "x \<le> (2 ^ prec) * 2 powr - real_of_int (int prec - \<lfloor>log 2 \<bar>y\<bar>\<rfloor>)"
2079     proof -
2080       have "x = 2 powr (log 2 \<bar>x\<bar>)" using \<open>0 < x\<close> by simp
2081       also have "\<dots> \<le> (2 ^ (Suc prec )) * 2 powr - real_of_int (int prec - \<lfloor>log 2 \<bar>x\<bar>\<rfloor>)"
2082         using real_of_int_floor_add_one_ge[of "log 2 \<bar>x\<bar>"] \<open>0 < x\<close>
2083         by (auto simp: powr_realpow[symmetric] powr_add[symmetric] algebra_simps
2084           powr_mult_base le_powr_iff)
2085       also
2086       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)"
2087         using logless flogless \<open>x > 0\<close> \<open>y > 0\<close>
2088         by (auto intro!: floor_mono)
2089       finally show ?thesis
2090         by (auto simp: powr_realpow[symmetric] powr_diff assms of_nat_diff)
2091     qed
2092     ultimately show ?thesis
2093       by (metis dual_order.trans truncate_down)
2094   qed
2095 qed
2097 lemma truncate_down_eq_truncate_up: "truncate_down p x = - truncate_up p (-x)"
2098   and truncate_up_eq_truncate_down: "truncate_up p x = - truncate_down p (-x)"
2099   by (auto simp: truncate_up_uminus_eq truncate_down_uminus_eq)
2101 lemma truncate_down_mono: "x \<le> y \<Longrightarrow> truncate_down p x \<le> truncate_down p y"
2102   apply (cases "0 \<le> x")
2103   apply (rule truncate_down_nonneg_mono, assumption+)
2105   apply (cases "0 \<le> y")
2106   apply (auto intro: truncate_up_nonneg_mono truncate_up_switch_sign_mono)
2107   done
2109 lemma truncate_up_mono: "x \<le> y \<Longrightarrow> truncate_up p x \<le> truncate_up p y"
2110   by (simp add: truncate_up_eq_truncate_down truncate_down_mono)
2112 lemma Float_le_zero_iff: "Float a b \<le> 0 \<longleftrightarrow> a \<le> 0"
2113   by (auto simp: zero_float_def mult_le_0_iff)
2115 lemma real_of_float_pprt[simp]:
2116   fixes a :: float
2117   shows "real_of_float (pprt a) = pprt (real_of_float a)"
2118   unfolding pprt_def sup_float_def max_def sup_real_def by auto
2120 lemma real_of_float_nprt[simp]:
2121   fixes a :: float
2122   shows "real_of_float (nprt a) = nprt (real_of_float a)"
2123   unfolding nprt_def inf_float_def min_def inf_real_def by auto
2125 context
2126 begin
2128 lift_definition int_floor_fl :: "float \<Rightarrow> int" is floor .
2130 qualified lemma compute_int_floor_fl[code]:
2131   "int_floor_fl (Float m e) = (if 0 \<le> e then m * 2 ^ nat e else m div (2 ^ (nat (-e))))"
2132   apply transfer
2133   apply (simp add: powr_int floor_divide_of_int_eq)
2134   apply (metis (no_types, hide_lams)floor_divide_of_int_eq of_int_numeral of_int_power floor_of_int of_int_mult)
2135   done
2137 lift_definition floor_fl :: "float \<Rightarrow> float" is "\<lambda>x. real_of_int \<lfloor>x\<rfloor>"
2138   by simp
2140 qualified lemma compute_floor_fl[code]:
2141   "floor_fl (Float m e) = (if 0 \<le> e then Float m e else Float (m div (2 ^ (nat (-e)))) 0)"
2142   apply transfer
2143   apply (simp add: powr_int floor_divide_of_int_eq)
2144   apply (metis (no_types, hide_lams)floor_divide_of_int_eq of_int_numeral of_int_power of_int_mult)
2145   done
2147 end
2149 lemma floor_fl: "real_of_float (floor_fl x) \<le> real_of_float x"
2150   by transfer simp
2152 lemma int_floor_fl: "real_of_int (int_floor_fl x) \<le> real_of_float x"
2153   by transfer simp
2155 lemma floor_pos_exp: "exponent (floor_fl x) \<ge> 0"
2156 proof (cases "floor_fl x = 0")
2157   case True
2158   then show ?thesis
2160 next
2161   case False
2162   have eq: "floor_fl x = Float \<lfloor>real_of_float x\<rfloor> 0"
2163     by transfer simp
2164   obtain i where "\<lfloor>real_of_float x\<rfloor> = mantissa (floor_fl x) * 2 ^ i" "0 = exponent (floor_fl x) - int i"
2165     by (rule denormalize_shift[OF eq False])
2166   then show ?thesis
2167     by simp
2168 qed
2170 lemma compute_mantissa[code]:
2171   "mantissa (Float m e) =
2172     (if m = 0 then 0 else if 2 dvd m then mantissa (normfloat (Float m e)) else m)"
2173   by (auto simp: mantissa_float Float.abs_eq zero_float_def[symmetric])
2175 lemma compute_exponent[code]:
2176   "exponent (Float m e) =
2177     (if m = 0 then 0 else if 2 dvd m then exponent (normfloat (Float m e)) else e)"
2178   by (auto simp: exponent_float Float.abs_eq zero_float_def[symmetric])
2180 lifting_update Float.float.lifting
2181 lifting_forget Float.float.lifting
2183 end