src/HOL/Library/Float.thy
author huffman
Fri Mar 30 12:32:35 2012 +0200 (2012-03-30)
changeset 47220 52426c62b5d0
parent 47165 9344891b504b
child 47230 6584098d5378
permissions -rw-r--r--
replace lemmas eval_nat_numeral with a simpler reformulation
     1 (*  Title:      HOL/Library/Float.thy
     2     Author:     Steven Obua 2008
     3     Author:     Johannes Hoelzl, TU Muenchen <hoelzl@in.tum.de> 2008 / 2009
     4 *)
     5 
     6 header {* Floating-Point Numbers *}
     7 
     8 theory Float
     9 imports Complex_Main Lattice_Algebras
    10 begin
    11 
    12 definition pow2 :: "int \<Rightarrow> real" where
    13   [simp]: "pow2 a = (if (0 <= a) then (2^(nat a)) else (inverse (2^(nat (-a)))))"
    14 
    15 datatype float = Float int int
    16 
    17 primrec of_float :: "float \<Rightarrow> real" where
    18   "of_float (Float a b) = real a * pow2 b"
    19 
    20 defs (overloaded)
    21   real_of_float_def [code_unfold]: "real == of_float"
    22 
    23 declare [[coercion "% x . Float x 0"]]
    24 declare [[coercion "real::float\<Rightarrow>real"]]
    25 
    26 primrec mantissa :: "float \<Rightarrow> int" where
    27   "mantissa (Float a b) = a"
    28 
    29 primrec scale :: "float \<Rightarrow> int" where
    30   "scale (Float a b) = b"
    31 
    32 instantiation float :: zero
    33 begin
    34 definition zero_float where "0 = Float 0 0"
    35 instance ..
    36 end
    37 
    38 instantiation float :: one
    39 begin
    40 definition one_float where "1 = Float 1 0"
    41 instance ..
    42 end
    43 
    44 lemma real_of_float_simp[simp]: "real (Float a b) = real a * pow2 b"
    45   unfolding real_of_float_def using of_float.simps .
    46 
    47 lemma real_of_float_neg_exp: "e < 0 \<Longrightarrow> real (Float m e) = real m * inverse (2^nat (-e))" by auto
    48 lemma real_of_float_nge0_exp: "\<not> 0 \<le> e \<Longrightarrow> real (Float m e) = real m * inverse (2^nat (-e))" by auto
    49 lemma real_of_float_ge0_exp: "0 \<le> e \<Longrightarrow> real (Float m e) = real m * (2^nat e)" by auto
    50 
    51 lemma Float_num[simp]: shows
    52    "real (Float 1 0) = 1" and "real (Float 1 1) = 2" and "real (Float 1 2) = 4" and
    53    "real (Float 1 -1) = 1/2" and "real (Float 1 -2) = 1/4" and "real (Float 1 -3) = 1/8" and
    54    "real (Float -1 0) = -1" and "real (Float (numeral n) 0) = numeral n"
    55   by auto
    56 
    57 lemma float_number_of_int[simp]: "real (Float n 0) = real n"
    58   by simp
    59 
    60 lemma pow2_0[simp]: "pow2 0 = 1" by simp
    61 lemma pow2_1[simp]: "pow2 1 = 2" by simp
    62 lemma pow2_neg: "pow2 x = inverse (pow2 (-x))" by simp
    63 
    64 lemma pow2_powr: "pow2 a = 2 powr a"
    65   by (simp add: powr_realpow[symmetric] powr_minus)
    66 
    67 declare pow2_def[simp del]
    68 
    69 lemma pow2_add: "pow2 (a+b) = (pow2 a) * (pow2 b)"
    70   by (simp add: pow2_powr powr_add)
    71 
    72 lemma pow2_diff: "pow2 (a - b) = pow2 a / pow2 b"
    73   by (simp add: pow2_powr powr_divide2)
    74   
    75 lemma pow2_add1: "pow2 (1 + a) = 2 * (pow2 a)"
    76   by (simp add: pow2_add)
    77 
    78 lemma float_components[simp]: "Float (mantissa f) (scale f) = f" by (cases f) auto
    79 
    80 lemma float_split: "\<exists> a b. x = Float a b" by (cases x) auto
    81 
    82 lemma float_split2: "(\<forall> a b. x \<noteq> Float a b) = False" by (auto simp add: float_split)
    83 
    84 lemma float_zero[simp]: "real (Float 0 e) = 0" by simp
    85 
    86 lemma abs_div_2_less: "a \<noteq> 0 \<Longrightarrow> a \<noteq> -1 \<Longrightarrow> abs((a::int) div 2) < abs a"
    87 by arith
    88 
    89 function normfloat :: "float \<Rightarrow> float" where
    90   "normfloat (Float a b) =
    91     (if a \<noteq> 0 \<and> even a then normfloat (Float (a div 2) (b+1))
    92      else if a=0 then Float 0 0 else Float a b)"
    93 by pat_completeness auto
    94 termination by (relation "measure (nat o abs o mantissa)") (auto intro: abs_div_2_less)
    95 declare normfloat.simps[simp del]
    96 
    97 theorem normfloat[symmetric, simp]: "real f = real (normfloat f)"
    98 proof (induct f rule: normfloat.induct)
    99   case (1 a b)
   100   have real2: "2 = real (2::int)"
   101     by auto
   102   show ?case
   103     apply (subst normfloat.simps)
   104     apply auto
   105     apply (subst 1[symmetric])
   106     apply (auto simp add: pow2_add even_def)
   107     done
   108 qed
   109 
   110 lemma pow2_neq_zero[simp]: "pow2 x \<noteq> 0"
   111   by (auto simp add: pow2_def)
   112 
   113 lemma pow2_int: "pow2 (int c) = 2^c"
   114   by (simp add: pow2_def)
   115 
   116 lemma zero_less_pow2[simp]: "0 < pow2 x"
   117   by (simp add: pow2_powr)
   118 
   119 lemma normfloat_imp_odd_or_zero:
   120   "normfloat f = Float a b \<Longrightarrow> odd a \<or> (a = 0 \<and> b = 0)"
   121 proof (induct f rule: normfloat.induct)
   122   case (1 u v)
   123   from 1 have ab: "normfloat (Float u v) = Float a b" by auto
   124   {
   125     assume eu: "even u"
   126     assume z: "u \<noteq> 0"
   127     have "normfloat (Float u v) = normfloat (Float (u div 2) (v + 1))"
   128       apply (subst normfloat.simps)
   129       by (simp add: eu z)
   130     with ab have "normfloat (Float (u div 2) (v + 1)) = Float a b" by simp
   131     with 1 eu z have ?case by auto
   132   }
   133   note case1 = this
   134   {
   135     assume "odd u \<or> u = 0"
   136     then have ou: "\<not> (u \<noteq> 0 \<and> even u)" by auto
   137     have "normfloat (Float u v) = (if u = 0 then Float 0 0 else Float u v)"
   138       apply (subst normfloat.simps)
   139       apply (simp add: ou)
   140       done
   141     with ab have "Float a b = (if u = 0 then Float 0 0 else Float u v)" by auto
   142     then have ?case
   143       apply (case_tac "u=0")
   144       apply (auto)
   145       by (insert ou, auto)
   146   }
   147   note case2 = this
   148   show ?case
   149     apply (case_tac "odd u \<or> u = 0")
   150     apply (rule case2)
   151     apply simp
   152     apply (rule case1)
   153     apply auto
   154     done
   155 qed
   156 
   157 lemma float_eq_odd_helper: 
   158   assumes odd: "odd a'"
   159     and floateq: "real (Float a b) = real (Float a' b')"
   160   shows "b \<le> b'"
   161 proof - 
   162   from odd have "a' \<noteq> 0" by auto
   163   with floateq have a': "real a' = real a * pow2 (b - b')"
   164     by (simp add: pow2_diff field_simps)
   165 
   166   {
   167     assume bcmp: "b > b'"
   168     then have "\<exists>c::nat. b - b' = int c + 1"
   169       by arith
   170     then guess c ..
   171     with a' have "real a' = real (a * 2^c * 2)"
   172       by (simp add: pow2_def nat_add_distrib)
   173     with odd have False
   174       unfolding real_of_int_inject by simp
   175   }
   176   then show ?thesis by arith
   177 qed
   178 
   179 lemma float_eq_odd: 
   180   assumes odd1: "odd a"
   181     and odd2: "odd a'"
   182     and floateq: "real (Float a b) = real (Float a' b')"
   183   shows "a = a' \<and> b = b'"
   184 proof -
   185   from 
   186      float_eq_odd_helper[OF odd2 floateq] 
   187      float_eq_odd_helper[OF odd1 floateq[symmetric]]
   188   have beq: "b = b'" by arith
   189   with floateq show ?thesis by auto
   190 qed
   191 
   192 theorem normfloat_unique:
   193   assumes real_of_float_eq: "real f = real g"
   194   shows "normfloat f = normfloat g"
   195 proof - 
   196   from float_split[of "normfloat f"] obtain a b where normf:"normfloat f = Float a b" by auto
   197   from float_split[of "normfloat g"] obtain a' b' where normg:"normfloat g = Float a' b'" by auto
   198   have "real (normfloat f) = real (normfloat g)"
   199     by (simp add: real_of_float_eq)
   200   then have float_eq: "real (Float a b) = real (Float a' b')"
   201     by (simp add: normf normg)
   202   have ab: "odd a \<or> (a = 0 \<and> b = 0)" by (rule normfloat_imp_odd_or_zero[OF normf])
   203   have ab': "odd a' \<or> (a' = 0 \<and> b' = 0)" by (rule normfloat_imp_odd_or_zero[OF normg])
   204   {
   205     assume odd: "odd a"
   206     then have "a \<noteq> 0" by (simp add: even_def) arith
   207     with float_eq have "a' \<noteq> 0" by auto
   208     with ab' have "odd a'" by simp
   209     from odd this float_eq have "a = a' \<and> b = b'" by (rule float_eq_odd)
   210   }
   211   note odd_case = this
   212   {
   213     assume even: "even a"
   214     with ab have a0: "a = 0" by simp
   215     with float_eq have a0': "a' = 0" by auto 
   216     from a0 a0' ab ab' have "a = a' \<and> b = b'" by auto
   217   }
   218   note even_case = this
   219   from odd_case even_case show ?thesis
   220     apply (simp add: normf normg)
   221     apply (case_tac "even a")
   222     apply auto
   223     done
   224 qed
   225 
   226 instantiation float :: plus
   227 begin
   228 fun plus_float where
   229 [simp del]: "(Float a_m a_e) + (Float b_m b_e) = 
   230      (if a_e \<le> b_e then Float (a_m + b_m * 2^(nat(b_e - a_e))) a_e 
   231                    else Float (a_m * 2^(nat (a_e - b_e)) + b_m) b_e)"
   232 instance ..
   233 end
   234 
   235 instantiation float :: uminus
   236 begin
   237 primrec uminus_float where [simp del]: "uminus_float (Float m e) = Float (-m) e"
   238 instance ..
   239 end
   240 
   241 instantiation float :: minus
   242 begin
   243 definition minus_float where "(z::float) - w = z + (- w)"
   244 instance ..
   245 end
   246 
   247 instantiation float :: times
   248 begin
   249 fun times_float where [simp del]: "(Float a_m a_e) * (Float b_m b_e) = Float (a_m * b_m) (a_e + b_e)"
   250 instance ..
   251 end
   252 
   253 primrec float_pprt :: "float \<Rightarrow> float" where
   254   "float_pprt (Float a e) = (if 0 <= a then (Float a e) else 0)"
   255 
   256 primrec float_nprt :: "float \<Rightarrow> float" where
   257   "float_nprt (Float a e) = (if 0 <= a then 0 else (Float a e))" 
   258 
   259 instantiation float :: ord
   260 begin
   261 definition le_float_def: "z \<le> (w :: float) \<equiv> real z \<le> real w"
   262 definition less_float_def: "z < (w :: float) \<equiv> real z < real w"
   263 instance ..
   264 end
   265 
   266 lemma real_of_float_add[simp]: "real (a + b) = real a + real (b :: float)"
   267   by (cases a, cases b) (simp add: algebra_simps plus_float.simps, 
   268       auto simp add: pow2_int[symmetric] pow2_add[symmetric])
   269 
   270 lemma real_of_float_minus[simp]: "real (- a) = - real (a :: float)"
   271   by (cases a) simp
   272 
   273 lemma real_of_float_sub[simp]: "real (a - b) = real a - real (b :: float)"
   274   by (cases a, cases b) (simp add: minus_float_def)
   275 
   276 lemma real_of_float_mult[simp]: "real (a*b) = real a * real (b :: float)"
   277   by (cases a, cases b) (simp add: times_float.simps pow2_add)
   278 
   279 lemma real_of_float_0[simp]: "real (0 :: float) = 0"
   280   by (auto simp add: zero_float_def)
   281 
   282 lemma real_of_float_1[simp]: "real (1 :: float) = 1"
   283   by (auto simp add: one_float_def)
   284 
   285 lemma zero_le_float:
   286   "(0 <= real (Float a b)) = (0 <= a)"
   287   apply auto
   288   apply (auto simp add: zero_le_mult_iff)
   289   apply (insert zero_less_pow2[of b])
   290   apply (simp_all)
   291   done
   292 
   293 lemma float_le_zero:
   294   "(real (Float a b) <= 0) = (a <= 0)"
   295   apply auto
   296   apply (auto simp add: mult_le_0_iff)
   297   apply (insert zero_less_pow2[of b])
   298   apply auto
   299   done
   300 
   301 lemma zero_less_float:
   302   "(0 < real (Float a b)) = (0 < a)"
   303   apply auto
   304   apply (auto simp add: zero_less_mult_iff)
   305   apply (insert zero_less_pow2[of b])
   306   apply (simp_all)
   307   done
   308 
   309 lemma float_less_zero:
   310   "(real (Float a b) < 0) = (a < 0)"
   311   apply auto
   312   apply (auto simp add: mult_less_0_iff)
   313   apply (insert zero_less_pow2[of b])
   314   apply (simp_all)
   315   done
   316 
   317 declare real_of_float_simp[simp del]
   318 
   319 lemma real_of_float_pprt[simp]: "real (float_pprt a) = pprt (real a)"
   320   by (cases a) (auto simp add: zero_le_float float_le_zero)
   321 
   322 lemma real_of_float_nprt[simp]: "real (float_nprt a) = nprt (real a)"
   323   by (cases a) (auto simp add: zero_le_float float_le_zero)
   324 
   325 instance float :: ab_semigroup_add
   326 proof (intro_classes)
   327   fix a b c :: float
   328   show "a + b + c = a + (b + c)"
   329     by (cases a, cases b, cases c)
   330       (auto simp add: algebra_simps power_add[symmetric] plus_float.simps)
   331 next
   332   fix a b :: float
   333   show "a + b = b + a"
   334     by (cases a, cases b) (simp add: plus_float.simps)
   335 qed
   336 
   337 instance float :: numeral ..
   338 
   339 lemma Float_add_same_scale: "Float x e + Float y e = Float (x + y) e"
   340   by (simp add: plus_float.simps)
   341 
   342 (* FIXME: define other constant for code_unfold_post *)
   343 lemma numeral_float_Float (*[code_unfold_post]*):
   344   "numeral k = Float (numeral k) 0"
   345   by (induct k, simp_all only: numeral.simps one_float_def
   346     Float_add_same_scale)
   347 
   348 lemma float_number_of[simp]: "real (numeral x :: float) = numeral x"
   349   by (simp only: numeral_float_Float Float_num)
   350 
   351 
   352 instance float :: comm_monoid_mult
   353 proof (intro_classes)
   354   fix a b c :: float
   355   show "a * b * c = a * (b * c)"
   356     by (cases a, cases b, cases c) (simp add: times_float.simps)
   357 next
   358   fix a b :: float
   359   show "a * b = b * a"
   360     by (cases a, cases b) (simp add: times_float.simps)
   361 next
   362   fix a :: float
   363   show "1 * a = a"
   364     by (cases a) (simp add: times_float.simps one_float_def)
   365 qed
   366 
   367 (* Floats do NOT form a cancel_semigroup_add: *)
   368 lemma "0 + Float 0 1 = 0 + Float 0 2"
   369   by (simp add: plus_float.simps zero_float_def)
   370 
   371 instance float :: comm_semiring
   372 proof (intro_classes)
   373   fix a b c :: float
   374   show "(a + b) * c = a * c + b * c"
   375     by (cases a, cases b, cases c) (simp add: plus_float.simps times_float.simps algebra_simps)
   376 qed
   377 
   378 (* Floats do NOT form an order, because "(x < y) = (x <= y & x <> y)" does NOT hold *)
   379 
   380 instance float :: zero_neq_one
   381 proof (intro_classes)
   382   show "(0::float) \<noteq> 1"
   383     by (simp add: zero_float_def one_float_def)
   384 qed
   385 
   386 lemma float_le_simp: "((x::float) \<le> y) = (0 \<le> y - x)"
   387   by (auto simp add: le_float_def)
   388 
   389 lemma float_less_simp: "((x::float) < y) = (0 < y - x)"
   390   by (auto simp add: less_float_def)
   391 
   392 lemma real_of_float_min: "real (min x y :: float) = min (real x) (real y)" unfolding min_def le_float_def by auto
   393 lemma real_of_float_max: "real (max a b :: float) = max (real a) (real b)" unfolding max_def le_float_def by auto
   394 
   395 lemma float_power: "real (x ^ n :: float) = real x ^ n"
   396   by (induct n) simp_all
   397 
   398 lemma zero_le_pow2[simp]: "0 \<le> pow2 s"
   399   apply (subgoal_tac "0 < pow2 s")
   400   apply (auto simp only:)
   401   apply auto
   402   done
   403 
   404 lemma pow2_less_0_eq_False[simp]: "(pow2 s < 0) = False"
   405   apply auto
   406   apply (subgoal_tac "0 \<le> pow2 s")
   407   apply simp
   408   apply simp
   409   done
   410 
   411 lemma pow2_le_0_eq_False[simp]: "(pow2 s \<le> 0) = False"
   412   apply auto
   413   apply (subgoal_tac "0 < pow2 s")
   414   apply simp
   415   apply simp
   416   done
   417 
   418 lemma float_pos_m_pos: "0 < Float m e \<Longrightarrow> 0 < m"
   419   unfolding less_float_def real_of_float_simp real_of_float_0 zero_less_mult_iff
   420   by auto
   421 
   422 lemma float_pos_less1_e_neg: assumes "0 < Float m e" and "Float m e < 1" shows "e < 0"
   423 proof -
   424   have "0 < m" using float_pos_m_pos `0 < Float m e` by auto
   425   hence "0 \<le> real m" and "1 \<le> real m" by auto
   426   
   427   show "e < 0"
   428   proof (rule ccontr)
   429     assume "\<not> e < 0" hence "0 \<le> e" by auto
   430     hence "1 \<le> pow2 e" unfolding pow2_def by auto
   431     from mult_mono[OF `1 \<le> real m` this `0 \<le> real m`]
   432     have "1 \<le> Float m e" by (simp add: le_float_def real_of_float_simp)
   433     thus False using `Float m e < 1` unfolding less_float_def le_float_def by auto
   434   qed
   435 qed
   436 
   437 lemma float_less1_mantissa_bound: assumes "0 < Float m e" "Float m e < 1" shows "m < 2^(nat (-e))"
   438 proof -
   439   have "e < 0" using float_pos_less1_e_neg assms by auto
   440   have "\<And>x. (0::real) < 2^x" by auto
   441   have "real m < 2^(nat (-e))" using `Float m e < 1`
   442     unfolding less_float_def real_of_float_neg_exp[OF `e < 0`] real_of_float_1
   443           real_mult_less_iff1[of _ _ 1, OF `0 < 2^(nat (-e))`, symmetric] 
   444           mult_assoc by auto
   445   thus ?thesis unfolding real_of_int_less_iff[symmetric] by auto
   446 qed
   447 
   448 function bitlen :: "int \<Rightarrow> int" where
   449 "bitlen 0 = 0" | 
   450 "bitlen -1 = 1" | 
   451 "0 < x \<Longrightarrow> bitlen x = 1 + (bitlen (x div 2))" | 
   452 "x < -1 \<Longrightarrow> bitlen x = 1 + (bitlen (x div 2))"
   453   apply (case_tac "x = 0 \<or> x = -1 \<or> x < -1 \<or> x > 0")
   454   apply auto
   455   done
   456 termination by (relation "measure (nat o abs)", auto)
   457 
   458 lemma bitlen_ge0: "0 \<le> bitlen x" by (induct x rule: bitlen.induct, auto)
   459 lemma bitlen_ge1: "x \<noteq> 0 \<Longrightarrow> 1 \<le> bitlen x" by (induct x rule: bitlen.induct, auto simp add: bitlen_ge0)
   460 
   461 lemma bitlen_bounds': assumes "0 < x" shows "2^nat (bitlen x - 1) \<le> x \<and> x + 1 \<le> 2^nat (bitlen x)" (is "?P x")
   462   using `0 < x`
   463 proof (induct x rule: bitlen.induct)
   464   fix x
   465   assume "0 < x" and hyp: "0 < x div 2 \<Longrightarrow> ?P (x div 2)" hence "0 \<le> x" and "x \<noteq> 0" by auto
   466   { fix x have "0 \<le> 1 + bitlen x" using bitlen_ge0[of x] by auto } note gt0_pls1 = this
   467 
   468   have "0 < (2::int)" by auto
   469 
   470   show "?P x"
   471   proof (cases "x = 1")
   472     case True show "?P x" unfolding True by auto
   473   next
   474     case False hence "2 \<le> x" using `0 < x` `x \<noteq> 1` by auto
   475     hence "2 div 2 \<le> x div 2" by (rule zdiv_mono1, auto)
   476     hence "0 < x div 2" and "x div 2 \<noteq> 0" by auto
   477     hence bitlen_s1_ge0: "0 \<le> bitlen (x div 2) - 1" using bitlen_ge1[OF `x div 2 \<noteq> 0`] by auto
   478 
   479     { from hyp[OF `0 < x div 2`]
   480       have "2 ^ nat (bitlen (x div 2) - 1) \<le> x div 2" by auto
   481       hence "2 ^ nat (bitlen (x div 2) - 1) * 2 \<le> x div 2 * 2" by (rule mult_right_mono, auto)
   482       also have "\<dots> \<le> x" using `0 < x` by auto
   483       finally have "2^nat (1 + bitlen (x div 2) - 1) \<le> x" unfolding power_Suc2[symmetric] Suc_nat_eq_nat_zadd1[OF bitlen_s1_ge0] by auto
   484     } moreover
   485     { have "x + 1 \<le> x - x mod 2 + 2"
   486       proof -
   487         have "x mod 2 < 2" using `0 < x` by auto
   488         hence "x < x - x mod 2 +  2" unfolding algebra_simps by auto
   489         thus ?thesis by auto
   490       qed
   491       also have "x - x mod 2 + 2 = (x div 2 + 1) * 2" unfolding algebra_simps using `0 < x` div_mod_equality[of x 2 0] by auto
   492       also have "\<dots> \<le> 2^nat (bitlen (x div 2)) * 2" using hyp[OF `0 < x div 2`, THEN conjunct2] by (rule mult_right_mono, auto)
   493       also have "\<dots> = 2^(1 + nat (bitlen (x div 2)))" unfolding power_Suc2[symmetric] by auto
   494       finally have "x + 1 \<le> 2^(1 + nat (bitlen (x div 2)))" .
   495     }
   496     ultimately show ?thesis
   497       unfolding bitlen.simps(3)[OF `0 < x`] nat_add_distrib[OF zero_le_one bitlen_ge0]
   498       unfolding add_commute nat_add_distrib[OF zero_le_one gt0_pls1]
   499       by auto
   500   qed
   501 next
   502   fix x :: int assume "x < -1" and "0 < x" hence False by auto
   503   thus "?P x" by auto
   504 qed auto
   505 
   506 lemma bitlen_bounds: assumes "0 < x" shows "2^nat (bitlen x - 1) \<le> x \<and> x < 2^nat (bitlen x)"
   507   using bitlen_bounds'[OF `0<x`] by auto
   508 
   509 lemma bitlen_div: assumes "0 < m" shows "1 \<le> real m / 2^nat (bitlen m - 1)" and "real m / 2^nat (bitlen m - 1) < 2"
   510 proof -
   511   let ?B = "2^nat(bitlen m - 1)"
   512 
   513   have "?B \<le> m" using bitlen_bounds[OF `0 <m`] ..
   514   hence "1 * ?B \<le> real m" unfolding real_of_int_le_iff[symmetric] by auto
   515   thus "1 \<le> real m / ?B" by auto
   516 
   517   have "m \<noteq> 0" using assms by auto
   518   have "0 \<le> bitlen m - 1" using bitlen_ge1[OF `m \<noteq> 0`] by auto
   519 
   520   have "m < 2^nat(bitlen m)" using bitlen_bounds[OF `0 <m`] ..
   521   also have "\<dots> = 2^nat(bitlen m - 1 + 1)" using bitlen_ge1[OF `m \<noteq> 0`] by auto
   522   also have "\<dots> = ?B * 2" unfolding nat_add_distrib[OF `0 \<le> bitlen m - 1` zero_le_one] by auto
   523   finally have "real m < 2 * ?B" unfolding real_of_int_less_iff[symmetric] by auto
   524   hence "real m / ?B < 2 * ?B / ?B" by (rule divide_strict_right_mono, auto)
   525   thus "real m / ?B < 2" by auto
   526 qed
   527 
   528 lemma float_gt1_scale: assumes "1 \<le> Float m e"
   529   shows "0 \<le> e + (bitlen m - 1)"
   530 proof (cases "0 \<le> e")
   531   have "0 < Float m e" using assms unfolding less_float_def le_float_def by auto
   532   hence "0 < m" using float_pos_m_pos by auto
   533   hence "m \<noteq> 0" by auto
   534   case True with bitlen_ge1[OF `m \<noteq> 0`] show ?thesis by auto
   535 next
   536   have "0 < Float m e" using assms unfolding less_float_def le_float_def by auto
   537   hence "0 < m" using float_pos_m_pos by auto
   538   hence "m \<noteq> 0" and "1 < (2::int)" by auto
   539   case False let ?S = "2^(nat (-e))"
   540   have "1 \<le> real m * inverse ?S" using assms unfolding le_float_def real_of_float_nge0_exp[OF False] by auto
   541   hence "1 * ?S \<le> real m * inverse ?S * ?S" by (rule mult_right_mono, auto)
   542   hence "?S \<le> real m" unfolding mult_assoc by auto
   543   hence "?S \<le> m" unfolding real_of_int_le_iff[symmetric] by auto
   544   from this bitlen_bounds[OF `0 < m`, THEN conjunct2]
   545   have "nat (-e) < (nat (bitlen m))" unfolding power_strict_increasing_iff[OF `1 < 2`, symmetric] by (rule order_le_less_trans)
   546   hence "-e < bitlen m" using False bitlen_ge0 by auto
   547   thus ?thesis by auto
   548 qed
   549 
   550 lemma normalized_float: assumes "m \<noteq> 0" shows "real (Float m (- (bitlen m - 1))) = real m / 2^nat (bitlen m - 1)"
   551 proof (cases "- (bitlen m - 1) = 0")
   552   case True show ?thesis unfolding real_of_float_simp pow2_def using True by auto
   553 next
   554   case False hence P: "\<not> 0 \<le> - (bitlen m - 1)" using bitlen_ge1[OF `m \<noteq> 0`] by auto
   555   show ?thesis unfolding real_of_float_nge0_exp[OF P] divide_inverse by auto
   556 qed
   557 
   558 (* BROKEN
   559 lemma bitlen_Pls: "bitlen (Int.Pls) = Int.Pls" by (subst Pls_def, subst Pls_def, simp)
   560 
   561 lemma bitlen_Min: "bitlen (Int.Min) = Int.Bit1 Int.Pls" by (subst Min_def, simp add: Bit1_def) 
   562 
   563 lemma bitlen_B0: "bitlen (Int.Bit0 b) = (if iszero b then Int.Pls else Int.succ (bitlen b))"
   564   apply (auto simp add: iszero_def succ_def)
   565   apply (simp add: Bit0_def Pls_def)
   566   apply (subst Bit0_def)
   567   apply simp
   568   apply (subgoal_tac "0 < 2 * b \<or> 2 * b < -1")
   569   apply auto
   570   done
   571 
   572 lemma bitlen_B1: "bitlen (Int.Bit1 b) = (if iszero (Int.succ b) then Int.Bit1 Int.Pls else Int.succ (bitlen b))"
   573 proof -
   574   have h: "! x. (2*x + 1) div 2 = (x::int)"
   575     by arith    
   576   show ?thesis
   577     apply (auto simp add: iszero_def succ_def)
   578     apply (subst Bit1_def)+
   579     apply simp
   580     apply (subgoal_tac "2 * b + 1 = -1")
   581     apply (simp only:)
   582     apply simp_all
   583     apply (subst Bit1_def)
   584     apply simp
   585     apply (subgoal_tac "0 < 2 * b + 1 \<or> 2 * b + 1 < -1")
   586     apply (auto simp add: h)
   587     done
   588 qed
   589 
   590 lemma bitlen_number_of: "bitlen (number_of w) = number_of (bitlen w)"
   591   by (simp add: number_of_is_id)
   592 BH *)
   593 
   594 lemma [code]: "bitlen x = 
   595      (if x = 0  then 0 
   596  else if x = -1 then 1 
   597                 else (1 + (bitlen (x div 2))))"
   598   by (cases "x = 0 \<or> x = -1 \<or> 0 < x") auto
   599 
   600 definition lapprox_posrat :: "nat \<Rightarrow> int \<Rightarrow> int \<Rightarrow> float"
   601 where
   602   "lapprox_posrat prec x y = 
   603    (let 
   604        l = nat (int prec + bitlen y - bitlen x) ;
   605        d = (x * 2^l) div y
   606     in normfloat (Float d (- (int l))))"
   607 
   608 lemma pow2_minus: "pow2 (-x) = inverse (pow2 x)"
   609   unfolding pow2_neg[of "-x"] by auto
   610 
   611 lemma lapprox_posrat: 
   612   assumes x: "0 \<le> x"
   613   and y: "0 < y"
   614   shows "real (lapprox_posrat prec x y) \<le> real x / real y"
   615 proof -
   616   let ?l = "nat (int prec + bitlen y - bitlen x)"
   617   
   618   have "real (x * 2^?l div y) * inverse (2^?l) \<le> (real (x * 2^?l) / real y) * inverse (2^?l)" 
   619     by (rule mult_right_mono, fact real_of_int_div4, simp)
   620   also have "\<dots> \<le> (real x / real y) * 2^?l * inverse (2^?l)" by auto
   621   finally have "real (x * 2^?l div y) * inverse (2^?l) \<le> real x / real y" unfolding mult_assoc by auto
   622   thus ?thesis unfolding lapprox_posrat_def Let_def normfloat real_of_float_simp
   623     unfolding pow2_minus pow2_int minus_minus .
   624 qed
   625 
   626 lemma real_of_int_div_mult: 
   627   fixes x y c :: int assumes "0 < y" and "0 < c"
   628   shows "real (x div y) \<le> real (x * c div y) * inverse (real c)"
   629 proof -
   630   have "c * (x div y) + 0 \<le> c * x div y" unfolding zdiv_zmult1_eq[of c x y]
   631     by (rule add_left_mono, 
   632         auto intro!: mult_nonneg_nonneg 
   633              simp add: pos_imp_zdiv_nonneg_iff[OF `0 < y`] `0 < c`[THEN less_imp_le] pos_mod_sign[OF `0 < y`])
   634   hence "real (x div y) * real c \<le> real (x * c div y)" 
   635     unfolding real_of_int_mult[symmetric] real_of_int_le_iff mult_commute by auto
   636   hence "real (x div y) * real c * inverse (real c) \<le> real (x * c div y) * inverse (real c)"
   637     using `0 < c` by auto
   638   thus ?thesis unfolding mult_assoc using `0 < c` by auto
   639 qed
   640 
   641 lemma lapprox_posrat_bottom: assumes "0 < y"
   642   shows "real (x div y) \<le> real (lapprox_posrat n x y)" 
   643 proof -
   644   have pow: "\<And>x. (0::int) < 2^x" by auto
   645   show ?thesis
   646     unfolding lapprox_posrat_def Let_def real_of_float_add normfloat real_of_float_simp pow2_minus pow2_int
   647     using real_of_int_div_mult[OF `0 < y` pow] by auto
   648 qed
   649 
   650 lemma lapprox_posrat_nonneg: assumes "0 \<le> x" and "0 < y"
   651   shows "0 \<le> real (lapprox_posrat n x y)" 
   652 proof -
   653   show ?thesis
   654     unfolding lapprox_posrat_def Let_def real_of_float_add normfloat real_of_float_simp pow2_minus pow2_int
   655     using pos_imp_zdiv_nonneg_iff[OF `0 < y`] assms by (auto intro!: mult_nonneg_nonneg)
   656 qed
   657 
   658 definition rapprox_posrat :: "nat \<Rightarrow> int \<Rightarrow> int \<Rightarrow> float"
   659 where
   660   "rapprox_posrat prec x y = (let
   661      l = nat (int prec + bitlen y - bitlen x) ;
   662      X = x * 2^l ;
   663      d = X div y ;
   664      m = X mod y
   665    in normfloat (Float (d + (if m = 0 then 0 else 1)) (- (int l))))"
   666 
   667 lemma rapprox_posrat:
   668   assumes x: "0 \<le> x"
   669   and y: "0 < y"
   670   shows "real x / real y \<le> real (rapprox_posrat prec x y)"
   671 proof -
   672   let ?l = "nat (int prec + bitlen y - bitlen x)" let ?X = "x * 2^?l"
   673   show ?thesis 
   674   proof (cases "?X mod y = 0")
   675     case True hence "y dvd ?X" using `0 < y` by auto
   676     from real_of_int_div[OF this]
   677     have "real (?X div y) * inverse (2 ^ ?l) = real ?X / real y * inverse (2 ^ ?l)" by auto
   678     also have "\<dots> = real x / real y * (2^?l * inverse (2^?l))" by auto
   679     finally have "real (?X div y) * inverse (2^?l) = real x / real y" by auto
   680     thus ?thesis unfolding rapprox_posrat_def Let_def normfloat if_P[OF True] 
   681       unfolding real_of_float_simp pow2_minus pow2_int minus_minus by auto
   682   next
   683     case False
   684     have "0 \<le> real y" and "real y \<noteq> 0" using `0 < y` by auto
   685     have "0 \<le> real y * 2^?l" by (rule mult_nonneg_nonneg, rule `0 \<le> real y`, auto)
   686 
   687     have "?X = y * (?X div y) + ?X mod y" by auto
   688     also have "\<dots> \<le> y * (?X div y) + y" by (rule add_mono, auto simp add: pos_mod_bound[OF `0 < y`, THEN less_imp_le])
   689     also have "\<dots> = y * (?X div y + 1)" unfolding right_distrib by auto
   690     finally have "real ?X \<le> real y * real (?X div y + 1)" unfolding real_of_int_le_iff real_of_int_mult[symmetric] .
   691     hence "real ?X / (real y * 2^?l) \<le> real y * real (?X div y + 1) / (real y * 2^?l)" 
   692       by (rule divide_right_mono, simp only: `0 \<le> real y * 2^?l`)
   693     also have "\<dots> = real y * real (?X div y + 1) / real y / 2^?l" by auto
   694     also have "\<dots> = real (?X div y + 1) * inverse (2^?l)" unfolding nonzero_mult_divide_cancel_left[OF `real y \<noteq> 0`] 
   695       unfolding divide_inverse ..
   696     finally show ?thesis unfolding rapprox_posrat_def Let_def normfloat real_of_float_simp if_not_P[OF False]
   697       unfolding pow2_minus pow2_int minus_minus by auto
   698   qed
   699 qed
   700 
   701 lemma rapprox_posrat_le1: assumes "0 \<le> x" and "0 < y" and "x \<le> y"
   702   shows "real (rapprox_posrat n x y) \<le> 1"
   703 proof -
   704   let ?l = "nat (int n + bitlen y - bitlen x)" let ?X = "x * 2^?l"
   705   show ?thesis
   706   proof (cases "?X mod y = 0")
   707     case True hence "y dvd ?X" using `0 < y` by auto
   708     from real_of_int_div[OF this]
   709     have "real (?X div y) * inverse (2 ^ ?l) = real ?X / real y * inverse (2 ^ ?l)" by auto
   710     also have "\<dots> = real x / real y * (2^?l * inverse (2^?l))" by auto
   711     finally have "real (?X div y) * inverse (2^?l) = real x / real y" by auto
   712     also have "real x / real y \<le> 1" using `0 \<le> x` and `0 < y` and `x \<le> y` by auto
   713     finally show ?thesis unfolding rapprox_posrat_def Let_def normfloat if_P[OF True]
   714       unfolding real_of_float_simp pow2_minus pow2_int minus_minus by auto
   715   next
   716     case False
   717     have "x \<noteq> y"
   718     proof (rule ccontr)
   719       assume "\<not> x \<noteq> y" hence "x = y" by auto
   720       have "?X mod y = 0" unfolding `x = y` using mod_mult_self1_is_0 by auto
   721       thus False using False by auto
   722     qed
   723     hence "x < y" using `x \<le> y` by auto
   724     hence "real x / real y < 1" using `0 < y` and `0 \<le> x` by auto
   725 
   726     from real_of_int_div4[of "?X" y]
   727     have "real (?X div y) \<le> (real x / real y) * 2^?l" unfolding real_of_int_mult times_divide_eq_left real_of_int_power real_numeral .
   728     also have "\<dots> < 1 * 2^?l" using `real x / real y < 1` by (rule mult_strict_right_mono, auto)
   729     finally have "?X div y < 2^?l" unfolding real_of_int_less_iff[of _ "2^?l", symmetric] by auto
   730     hence "?X div y + 1 \<le> 2^?l" by auto
   731     hence "real (?X div y + 1) * inverse (2^?l) \<le> 2^?l * inverse (2^?l)"
   732       unfolding real_of_int_le_iff[of _ "2^?l", symmetric] real_of_int_power real_numeral
   733       by (rule mult_right_mono, auto)
   734     hence "real (?X div y + 1) * inverse (2^?l) \<le> 1" by auto
   735     thus ?thesis unfolding rapprox_posrat_def Let_def normfloat real_of_float_simp if_not_P[OF False]
   736       unfolding pow2_minus pow2_int minus_minus by auto
   737   qed
   738 qed
   739 
   740 lemma zdiv_greater_zero: fixes a b :: int assumes "0 < a" and "a \<le> b"
   741   shows "0 < b div a"
   742 proof (rule ccontr)
   743   have "0 \<le> b" using assms by auto
   744   assume "\<not> 0 < b div a" hence "b div a = 0" using `0 \<le> b`[unfolded pos_imp_zdiv_nonneg_iff[OF `0<a`, of b, symmetric]] by auto
   745   have "b = a * (b div a) + b mod a" by auto
   746   hence "b = b mod a" unfolding `b div a = 0` by auto
   747   hence "b < a" using `0 < a`[THEN pos_mod_bound, of b] by auto
   748   thus False using `a \<le> b` by auto
   749 qed
   750 
   751 lemma rapprox_posrat_less1: assumes "0 \<le> x" and "0 < y" and "2 * x < y" and "0 < n"
   752   shows "real (rapprox_posrat n x y) < 1"
   753 proof (cases "x = 0")
   754   case True thus ?thesis unfolding rapprox_posrat_def True Let_def normfloat real_of_float_simp by auto
   755 next
   756   case False hence "0 < x" using `0 \<le> x` by auto
   757   hence "x < y" using assms by auto
   758   
   759   let ?l = "nat (int n + bitlen y - bitlen x)" let ?X = "x * 2^?l"
   760   show ?thesis
   761   proof (cases "?X mod y = 0")
   762     case True hence "y dvd ?X" using `0 < y` by auto
   763     from real_of_int_div[OF this]
   764     have "real (?X div y) * inverse (2 ^ ?l) = real ?X / real y * inverse (2 ^ ?l)" by auto
   765     also have "\<dots> = real x / real y * (2^?l * inverse (2^?l))" by auto
   766     finally have "real (?X div y) * inverse (2^?l) = real x / real y" by auto
   767     also have "real x / real y < 1" using `0 \<le> x` and `0 < y` and `x < y` by auto
   768     finally show ?thesis unfolding rapprox_posrat_def Let_def normfloat real_of_float_simp if_P[OF True]
   769       unfolding pow2_minus pow2_int minus_minus by auto
   770   next
   771     case False
   772     hence "(real x / real y) < 1 / 2" using `0 < y` and `0 \<le> x` `2 * x < y` by auto
   773 
   774     have "0 < ?X div y"
   775     proof -
   776       have "2^nat (bitlen x - 1) \<le> y" and "y < 2^nat (bitlen y)"
   777         using bitlen_bounds[OF `0 < x`, THEN conjunct1] bitlen_bounds[OF `0 < y`, THEN conjunct2] `x < y` by auto
   778       hence "(2::int)^nat (bitlen x - 1) < 2^nat (bitlen y)" by (rule order_le_less_trans)
   779       hence "bitlen x \<le> bitlen y" by auto
   780       hence len_less: "nat (bitlen x - 1) \<le> nat (int (n - 1) + bitlen y)" by auto
   781 
   782       have "x \<noteq> 0" and "y \<noteq> 0" using `0 < x` `0 < y` by auto
   783 
   784       have exp_eq: "nat (int (n - 1) + bitlen y) - nat (bitlen x - 1) = ?l"
   785         using `bitlen x \<le> bitlen y` bitlen_ge1[OF `x \<noteq> 0`] bitlen_ge1[OF `y \<noteq> 0`] `0 < n` by auto
   786 
   787       have "y * 2^nat (bitlen x - 1) \<le> y * x" 
   788         using bitlen_bounds[OF `0 < x`, THEN conjunct1] `0 < y`[THEN less_imp_le] by (rule mult_left_mono)
   789       also have "\<dots> \<le> 2^nat (bitlen y) * x" using bitlen_bounds[OF `0 < y`, THEN conjunct2, THEN less_imp_le] `0 \<le> x` by (rule mult_right_mono)
   790       also have "\<dots> \<le> x * 2^nat (int (n - 1) + bitlen y)" unfolding mult_commute[of x] by (rule mult_right_mono, auto simp add: `0 \<le> x`)
   791       finally have "real y * 2^nat (bitlen x - 1) * inverse (2^nat (bitlen x - 1)) \<le> real x * 2^nat (int (n - 1) + bitlen y) * inverse (2^nat (bitlen x - 1))"
   792         unfolding real_of_int_le_iff[symmetric] by auto
   793       hence "real y \<le> real x * (2^nat (int (n - 1) + bitlen y) / (2^nat (bitlen x - 1)))" 
   794         unfolding mult_assoc divide_inverse by auto
   795       also have "\<dots> = real x * (2^(nat (int (n - 1) + bitlen y) - nat (bitlen x - 1)))" using power_diff[of "2::real", OF _ len_less] by auto
   796       finally have "y \<le> x * 2^?l" unfolding exp_eq unfolding real_of_int_le_iff[symmetric] by auto
   797       thus ?thesis using zdiv_greater_zero[OF `0 < y`] by auto
   798     qed
   799 
   800     from real_of_int_div4[of "?X" y]
   801     have "real (?X div y) \<le> (real x / real y) * 2^?l" unfolding real_of_int_mult times_divide_eq_left real_of_int_power real_numeral .
   802     also have "\<dots> < 1/2 * 2^?l" using `real x / real y < 1/2` by (rule mult_strict_right_mono, auto)
   803     finally have "?X div y * 2 < 2^?l" unfolding real_of_int_less_iff[of _ "2^?l", symmetric] by auto
   804     hence "?X div y + 1 < 2^?l" using `0 < ?X div y` by auto
   805     hence "real (?X div y + 1) * inverse (2^?l) < 2^?l * inverse (2^?l)"
   806       unfolding real_of_int_less_iff[of _ "2^?l", symmetric] real_of_int_power real_numeral
   807       by (rule mult_strict_right_mono, auto)
   808     hence "real (?X div y + 1) * inverse (2^?l) < 1" by auto
   809     thus ?thesis unfolding rapprox_posrat_def Let_def normfloat real_of_float_simp if_not_P[OF False]
   810       unfolding pow2_minus pow2_int minus_minus by auto
   811   qed
   812 qed
   813 
   814 lemma approx_rat_pattern: fixes P and ps :: "nat * int * int"
   815   assumes Y: "\<And>y prec x. \<lbrakk>y = 0; ps = (prec, x, 0)\<rbrakk> \<Longrightarrow> P" 
   816   and A: "\<And>x y prec. \<lbrakk>0 \<le> x; 0 < y; ps = (prec, x, y)\<rbrakk> \<Longrightarrow> P"
   817   and B: "\<And>x y prec. \<lbrakk>x < 0; 0 < y; ps = (prec, x, y)\<rbrakk> \<Longrightarrow> P"
   818   and C: "\<And>x y prec. \<lbrakk>x < 0; y < 0; ps = (prec, x, y)\<rbrakk> \<Longrightarrow> P"
   819   and D: "\<And>x y prec. \<lbrakk>0 \<le> x; y < 0; ps = (prec, x, y)\<rbrakk> \<Longrightarrow> P"
   820   shows P
   821 proof -
   822   obtain prec x y where [simp]: "ps = (prec, x, y)" by (cases ps) auto
   823   from Y have "y = 0 \<Longrightarrow> P" by auto
   824   moreover {
   825     assume "0 < y"
   826     have P
   827     proof (cases "0 \<le> x")
   828       case True
   829       with A and `0 < y` show P by auto
   830     next
   831       case False
   832       with B and `0 < y` show P by auto
   833     qed
   834   } 
   835   moreover {
   836     assume "y < 0"
   837     have P
   838     proof (cases "0 \<le> x")
   839       case True
   840       with D and `y < 0` show P by auto
   841     next
   842       case False
   843       with C and `y < 0` show P by auto
   844     qed
   845   }
   846   ultimately show P by (cases "y = 0 \<or> 0 < y \<or> y < 0") auto
   847 qed
   848 
   849 function lapprox_rat :: "nat \<Rightarrow> int \<Rightarrow> int \<Rightarrow> float"
   850 where
   851   "y = 0 \<Longrightarrow> lapprox_rat prec x y = 0"
   852 | "0 \<le> x \<Longrightarrow> 0 < y \<Longrightarrow> lapprox_rat prec x y = lapprox_posrat prec x y"
   853 | "x < 0 \<Longrightarrow> 0 < y \<Longrightarrow> lapprox_rat prec x y = - (rapprox_posrat prec (-x) y)"
   854 | "x < 0 \<Longrightarrow> y < 0 \<Longrightarrow> lapprox_rat prec x y = lapprox_posrat prec (-x) (-y)"
   855 | "0 \<le> x \<Longrightarrow> y < 0 \<Longrightarrow> lapprox_rat prec x y = - (rapprox_posrat prec x (-y))"
   856 apply simp_all by (rule approx_rat_pattern)
   857 termination by lexicographic_order
   858 
   859 lemma compute_lapprox_rat[code]:
   860       "lapprox_rat prec x y = (if y = 0 then 0 else if 0 \<le> x then (if 0 < y then lapprox_posrat prec x y else - (rapprox_posrat prec x (-y))) 
   861                                                              else (if 0 < y then - (rapprox_posrat prec (-x) y) else lapprox_posrat prec (-x) (-y)))"
   862   by auto
   863             
   864 lemma lapprox_rat: "real (lapprox_rat prec x y) \<le> real x / real y"
   865 proof -      
   866   have h[rule_format]: "! a b b'. b' \<le> b \<longrightarrow> a \<le> b' \<longrightarrow> a \<le> (b::real)" by auto
   867   show ?thesis
   868     apply (case_tac "y = 0")
   869     apply simp
   870     apply (case_tac "0 \<le> x \<and> 0 < y")
   871     apply (simp add: lapprox_posrat)
   872     apply (case_tac "x < 0 \<and> 0 < y")
   873     apply simp
   874     apply (subst minus_le_iff)   
   875     apply (rule h[OF rapprox_posrat])
   876     apply (simp_all)
   877     apply (case_tac "x < 0 \<and> y < 0")
   878     apply simp
   879     apply (rule h[OF _ lapprox_posrat])
   880     apply (simp_all)
   881     apply (case_tac "0 \<le> x \<and> y < 0")
   882     apply (simp)
   883     apply (subst minus_le_iff)   
   884     apply (rule h[OF rapprox_posrat])
   885     apply simp_all
   886     apply arith
   887     done
   888 qed
   889 
   890 lemma lapprox_rat_bottom: assumes "0 \<le> x" and "0 < y"
   891   shows "real (x div y) \<le> real (lapprox_rat n x y)" 
   892   unfolding lapprox_rat.simps(2)[OF assms]  using lapprox_posrat_bottom[OF `0<y`] .
   893 
   894 function rapprox_rat :: "nat \<Rightarrow> int \<Rightarrow> int \<Rightarrow> float"
   895 where
   896   "y = 0 \<Longrightarrow> rapprox_rat prec x y = 0"
   897 | "0 \<le> x \<Longrightarrow> 0 < y \<Longrightarrow> rapprox_rat prec x y = rapprox_posrat prec x y"
   898 | "x < 0 \<Longrightarrow> 0 < y \<Longrightarrow> rapprox_rat prec x y = - (lapprox_posrat prec (-x) y)"
   899 | "x < 0 \<Longrightarrow> y < 0 \<Longrightarrow> rapprox_rat prec x y = rapprox_posrat prec (-x) (-y)"
   900 | "0 \<le> x \<Longrightarrow> y < 0 \<Longrightarrow> rapprox_rat prec x y = - (lapprox_posrat prec x (-y))"
   901 apply simp_all by (rule approx_rat_pattern)
   902 termination by lexicographic_order
   903 
   904 lemma compute_rapprox_rat[code]:
   905       "rapprox_rat prec x y = (if y = 0 then 0 else if 0 \<le> x then (if 0 < y then rapprox_posrat prec x y else - (lapprox_posrat prec x (-y))) else 
   906                                                                   (if 0 < y then - (lapprox_posrat prec (-x) y) else rapprox_posrat prec (-x) (-y)))"
   907   by auto
   908 
   909 lemma rapprox_rat: "real x / real y \<le> real (rapprox_rat prec x y)"
   910 proof -      
   911   have h[rule_format]: "! a b b'. b' \<le> b \<longrightarrow> a \<le> b' \<longrightarrow> a \<le> (b::real)" by auto
   912   show ?thesis
   913     apply (case_tac "y = 0")
   914     apply simp
   915     apply (case_tac "0 \<le> x \<and> 0 < y")
   916     apply (simp add: rapprox_posrat)
   917     apply (case_tac "x < 0 \<and> 0 < y")
   918     apply simp
   919     apply (subst le_minus_iff)   
   920     apply (rule h[OF _ lapprox_posrat])
   921     apply (simp_all)
   922     apply (case_tac "x < 0 \<and> y < 0")
   923     apply simp
   924     apply (rule h[OF rapprox_posrat])
   925     apply (simp_all)
   926     apply (case_tac "0 \<le> x \<and> y < 0")
   927     apply (simp)
   928     apply (subst le_minus_iff)   
   929     apply (rule h[OF _ lapprox_posrat])
   930     apply simp_all
   931     apply arith
   932     done
   933 qed
   934 
   935 lemma rapprox_rat_le1: assumes "0 \<le> x" and "0 < y" and "x \<le> y"
   936   shows "real (rapprox_rat n x y) \<le> 1"
   937   unfolding rapprox_rat.simps(2)[OF `0 \<le> x` `0 < y`] using rapprox_posrat_le1[OF assms] .
   938 
   939 lemma rapprox_rat_neg: assumes "x < 0" and "0 < y"
   940   shows "real (rapprox_rat n x y) \<le> 0"
   941   unfolding rapprox_rat.simps(3)[OF assms] using lapprox_posrat_nonneg[of "-x" y n] assms by auto
   942 
   943 lemma rapprox_rat_nonneg_neg: assumes "0 \<le> x" and "y < 0"
   944   shows "real (rapprox_rat n x y) \<le> 0"
   945   unfolding rapprox_rat.simps(5)[OF assms] using lapprox_posrat_nonneg[of x "-y" n] assms by auto
   946 
   947 lemma rapprox_rat_nonpos_pos: assumes "x \<le> 0" and "0 < y"
   948   shows "real (rapprox_rat n x y) \<le> 0"
   949 proof (cases "x = 0") 
   950   case True
   951   hence "0 \<le> x" by auto show ?thesis
   952     unfolding rapprox_rat.simps(2)[OF `0 \<le> x` `0 < y`]
   953     unfolding True rapprox_posrat_def Let_def
   954     by auto
   955 next
   956   case False
   957   hence "x < 0" using assms by auto
   958   show ?thesis using rapprox_rat_neg[OF `x < 0` `0 < y`] .
   959 qed
   960 
   961 fun float_divl :: "nat \<Rightarrow> float \<Rightarrow> float \<Rightarrow> float"
   962 where
   963   "float_divl prec (Float m1 s1) (Float m2 s2) = 
   964     (let
   965        l = lapprox_rat prec m1 m2;
   966        f = Float 1 (s1 - s2)
   967      in
   968        f * l)"     
   969 
   970 lemma float_divl: "real (float_divl prec x y) \<le> real x / real y"
   971   using lapprox_rat[of prec "mantissa x" "mantissa y"]
   972   by (cases x y rule: float.exhaust[case_product float.exhaust])
   973      (simp split: split_if_asm
   974            add: real_of_float_simp pow2_diff field_simps le_divide_eq mult_less_0_iff zero_less_mult_iff)
   975 
   976 lemma float_divl_lower_bound: assumes "0 \<le> x" and "0 < y" shows "0 \<le> float_divl prec x y"
   977 proof (cases x, cases y)
   978   fix xm xe ym ye :: int
   979   assume x_eq: "x = Float xm xe" and y_eq: "y = Float ym ye"
   980   have "0 \<le> xm"
   981     using `0 \<le> x`[unfolded x_eq le_float_def real_of_float_simp real_of_float_0 zero_le_mult_iff]
   982     by auto
   983   have "0 < ym"
   984     using `0 < y`[unfolded y_eq less_float_def real_of_float_simp real_of_float_0 zero_less_mult_iff]
   985     by auto
   986 
   987   have "\<And>n. 0 \<le> real (Float 1 n)"
   988     unfolding real_of_float_simp using zero_le_pow2 by auto
   989   moreover have "0 \<le> real (lapprox_rat prec xm ym)"
   990     apply (rule order_trans[OF _ lapprox_rat_bottom[OF `0 \<le> xm` `0 < ym`]])
   991     apply (auto simp add: `0 \<le> xm` pos_imp_zdiv_nonneg_iff[OF `0 < ym`])
   992     done
   993   ultimately show "0 \<le> float_divl prec x y"
   994     unfolding x_eq y_eq float_divl.simps Let_def le_float_def real_of_float_0
   995     by (auto intro!: mult_nonneg_nonneg)
   996 qed
   997 
   998 lemma float_divl_pos_less1_bound:
   999   assumes "0 < x" and "x < 1" and "0 < prec"
  1000   shows "1 \<le> float_divl prec 1 x"
  1001 proof (cases x)
  1002   case (Float m e)
  1003   from `0 < x` `x < 1` have "0 < m" "e < 0"
  1004     using float_pos_m_pos float_pos_less1_e_neg unfolding Float by auto
  1005   let ?b = "nat (bitlen m)" and ?e = "nat (-e)"
  1006   have "1 \<le> m" and "m \<noteq> 0" using `0 < m` by auto
  1007   with bitlen_bounds[OF `0 < m`] have "m < 2^?b" and "(2::int) \<le> 2^?b" by auto
  1008   hence "1 \<le> bitlen m" using power_le_imp_le_exp[of "2::int" 1 ?b] by auto
  1009   hence pow_split: "nat (int prec + bitlen m - 1) = (prec - 1) + ?b" using `0 < prec` by auto
  1010   
  1011   have pow_not0: "\<And>x. (2::real)^x \<noteq> 0" by auto
  1012 
  1013   from float_less1_mantissa_bound `0 < x` `x < 1` Float 
  1014   have "m < 2^?e" by auto
  1015   with bitlen_bounds[OF `0 < m`, THEN conjunct1] have "(2::int)^nat (bitlen m - 1) < 2^?e"
  1016     by (rule order_le_less_trans)
  1017   from power_less_imp_less_exp[OF _ this]
  1018   have "bitlen m \<le> - e" by auto
  1019   hence "(2::real)^?b \<le> 2^?e" by auto
  1020   hence "(2::real)^?b * inverse (2^?b) \<le> 2^?e * inverse (2^?b)"
  1021     by (rule mult_right_mono) auto
  1022   hence "(1::real) \<le> 2^?e * inverse (2^?b)" by auto
  1023   also
  1024   let ?d = "real (2 ^ nat (int prec + bitlen m - 1) div m) * inverse (2 ^ nat (int prec + bitlen m - 1))"
  1025   {
  1026     have "2^(prec - 1) * m \<le> 2^(prec - 1) * 2^?b"
  1027       using `m < 2^?b`[THEN less_imp_le] by (rule mult_left_mono) auto
  1028     also have "\<dots> = 2 ^ nat (int prec + bitlen m - 1)"
  1029       unfolding pow_split power_add by auto
  1030     finally have "2^(prec - 1) * m div m \<le> 2 ^ nat (int prec + bitlen m - 1) div m"
  1031       using `0 < m` by (rule zdiv_mono1)
  1032     hence "2^(prec - 1) \<le> 2 ^ nat (int prec + bitlen m - 1) div m"
  1033       unfolding div_mult_self2_is_id[OF `m \<noteq> 0`] .
  1034     hence "2^(prec - 1) * inverse (2 ^ nat (int prec + bitlen m - 1)) \<le> ?d"
  1035       unfolding real_of_int_le_iff[of "2^(prec - 1)", symmetric] by auto
  1036   }
  1037   from mult_left_mono[OF this [unfolded pow_split power_add inverse_mult_distrib mult_assoc[symmetric] right_inverse[OF pow_not0] mult_1_left], of "2^?e"]
  1038   have "2^?e * inverse (2^?b) \<le> 2^?e * ?d" unfolding pow_split power_add by auto
  1039   finally have "1 \<le> 2^?e * ?d" .
  1040   
  1041   have e_nat: "0 - e = int (nat (-e))" using `e < 0` by auto
  1042   have "bitlen 1 = 1" using bitlen.simps by auto
  1043   
  1044   show ?thesis 
  1045     unfolding one_float_def Float float_divl.simps Let_def
  1046       lapprox_rat.simps(2)[OF zero_le_one `0 < m`]
  1047       lapprox_posrat_def `bitlen 1 = 1`
  1048     unfolding le_float_def real_of_float_mult normfloat real_of_float_simp
  1049       pow2_minus pow2_int e_nat
  1050     using `1 \<le> 2^?e * ?d` by (auto simp add: pow2_def)
  1051 qed
  1052 
  1053 fun float_divr :: "nat \<Rightarrow> float \<Rightarrow> float \<Rightarrow> float"
  1054 where
  1055   "float_divr prec (Float m1 s1) (Float m2 s2) = 
  1056     (let
  1057        r = rapprox_rat prec m1 m2;
  1058        f = Float 1 (s1 - s2)
  1059      in
  1060        f * r)"  
  1061 
  1062 lemma float_divr: "real x / real y \<le> real (float_divr prec x y)"
  1063   using rapprox_rat[of "mantissa x" "mantissa y" prec]
  1064   by (cases x y rule: float.exhaust[case_product float.exhaust])
  1065      (simp split: split_if_asm
  1066            add: real_of_float_simp pow2_diff field_simps divide_le_eq mult_less_0_iff zero_less_mult_iff)
  1067 
  1068 lemma float_divr_pos_less1_lower_bound: assumes "0 < x" and "x < 1" shows "1 \<le> float_divr prec 1 x"
  1069 proof -
  1070   have "1 \<le> 1 / real x" using `0 < x` and `x < 1` unfolding less_float_def by auto
  1071   also have "\<dots> \<le> real (float_divr prec 1 x)" using float_divr[where x=1 and y=x] by auto
  1072   finally show ?thesis unfolding le_float_def by auto
  1073 qed
  1074 
  1075 lemma float_divr_nonpos_pos_upper_bound: assumes "x \<le> 0" and "0 < y" shows "float_divr prec x y \<le> 0"
  1076 proof (cases x, cases y)
  1077   fix xm xe ym ye :: int
  1078   assume x_eq: "x = Float xm xe" and y_eq: "y = Float ym ye"
  1079   have "xm \<le> 0" using `x \<le> 0`[unfolded x_eq le_float_def real_of_float_simp real_of_float_0 mult_le_0_iff] by auto
  1080   have "0 < ym" using `0 < y`[unfolded y_eq less_float_def real_of_float_simp real_of_float_0 zero_less_mult_iff] by auto
  1081 
  1082   have "\<And>n. 0 \<le> real (Float 1 n)" unfolding real_of_float_simp using zero_le_pow2 by auto
  1083   moreover have "real (rapprox_rat prec xm ym) \<le> 0" using rapprox_rat_nonpos_pos[OF `xm \<le> 0` `0 < ym`] .
  1084   ultimately show "float_divr prec x y \<le> 0"
  1085     unfolding x_eq y_eq float_divr.simps Let_def le_float_def real_of_float_0 real_of_float_mult by (auto intro!: mult_nonneg_nonpos)
  1086 qed
  1087 
  1088 lemma float_divr_nonneg_neg_upper_bound: assumes "0 \<le> x" and "y < 0" shows "float_divr prec x y \<le> 0"
  1089 proof (cases x, cases y)
  1090   fix xm xe ym ye :: int
  1091   assume x_eq: "x = Float xm xe" and y_eq: "y = Float ym ye"
  1092   have "0 \<le> xm" using `0 \<le> x`[unfolded x_eq le_float_def real_of_float_simp real_of_float_0 zero_le_mult_iff] by auto
  1093   have "ym < 0" using `y < 0`[unfolded y_eq less_float_def real_of_float_simp real_of_float_0 mult_less_0_iff] by auto
  1094   hence "0 < - ym" by auto
  1095 
  1096   have "\<And>n. 0 \<le> real (Float 1 n)" unfolding real_of_float_simp using zero_le_pow2 by auto
  1097   moreover have "real (rapprox_rat prec xm ym) \<le> 0" using rapprox_rat_nonneg_neg[OF `0 \<le> xm` `ym < 0`] .
  1098   ultimately show "float_divr prec x y \<le> 0"
  1099     unfolding x_eq y_eq float_divr.simps Let_def le_float_def real_of_float_0 real_of_float_mult by (auto intro!: mult_nonneg_nonpos)
  1100 qed
  1101 
  1102 primrec round_down :: "nat \<Rightarrow> float \<Rightarrow> float" where
  1103 "round_down prec (Float m e) = (let d = bitlen m - int prec in
  1104      if 0 < d then let P = 2^nat d ; n = m div P in Float n (e + d)
  1105               else Float m e)"
  1106 
  1107 primrec round_up :: "nat \<Rightarrow> float \<Rightarrow> float" where
  1108 "round_up prec (Float m e) = (let d = bitlen m - int prec in
  1109   if 0 < d then let P = 2^nat d ; n = m div P ; r = m mod P in Float (n + (if r = 0 then 0 else 1)) (e + d) 
  1110            else Float m e)"
  1111 
  1112 lemma round_up: "real x \<le> real (round_up prec x)"
  1113 proof (cases x)
  1114   case (Float m e)
  1115   let ?d = "bitlen m - int prec"
  1116   let ?p = "(2::int)^nat ?d"
  1117   have "0 < ?p" by auto
  1118   show "?thesis"
  1119   proof (cases "0 < ?d")
  1120     case True
  1121     hence pow_d: "pow2 ?d = real ?p" using pow2_int[symmetric] by simp
  1122     show ?thesis
  1123     proof (cases "m mod ?p = 0")
  1124       case True
  1125       have m: "m = m div ?p * ?p + 0" unfolding True[symmetric] using mod_div_equality [symmetric] .
  1126       have "real (Float m e) = real (Float (m div ?p) (e + ?d))" unfolding real_of_float_simp arg_cong[OF m, of real]
  1127         by (auto simp add: pow2_add `0 < ?d` pow_d)
  1128       thus ?thesis
  1129         unfolding Float round_up.simps Let_def if_P[OF `m mod ?p = 0`] if_P[OF `0 < ?d`]
  1130         by auto
  1131     next
  1132       case False
  1133       have "m = m div ?p * ?p + m mod ?p" unfolding mod_div_equality ..
  1134       also have "\<dots> \<le> (m div ?p + 1) * ?p" unfolding left_distrib mult_1 by (rule add_left_mono, rule pos_mod_bound[OF `0 < ?p`, THEN less_imp_le])
  1135       finally have "real (Float m e) \<le> real (Float (m div ?p + 1) (e + ?d))" unfolding real_of_float_simp add_commute[of e]
  1136         unfolding pow2_add mult_assoc[symmetric] real_of_int_le_iff[of m, symmetric]
  1137         by (auto intro!: mult_mono simp add: pow2_add `0 < ?d` pow_d)
  1138       thus ?thesis
  1139         unfolding Float round_up.simps Let_def if_not_P[OF `\<not> m mod ?p = 0`] if_P[OF `0 < ?d`] .
  1140     qed
  1141   next
  1142     case False
  1143     show ?thesis
  1144       unfolding Float round_up.simps Let_def if_not_P[OF False] .. 
  1145   qed
  1146 qed
  1147 
  1148 lemma round_down: "real (round_down prec x) \<le> real x"
  1149 proof (cases x)
  1150   case (Float m e)
  1151   let ?d = "bitlen m - int prec"
  1152   let ?p = "(2::int)^nat ?d"
  1153   have "0 < ?p" by auto
  1154   show "?thesis"
  1155   proof (cases "0 < ?d")
  1156     case True
  1157     hence pow_d: "pow2 ?d = real ?p" using pow2_int[symmetric] by simp
  1158     have "m div ?p * ?p \<le> m div ?p * ?p + m mod ?p" by (auto simp add: pos_mod_bound[OF `0 < ?p`, THEN less_imp_le])
  1159     also have "\<dots> \<le> m" unfolding mod_div_equality ..
  1160     finally have "real (Float (m div ?p) (e + ?d)) \<le> real (Float m e)" unfolding real_of_float_simp add_commute[of e]
  1161       unfolding pow2_add mult_assoc[symmetric] real_of_int_le_iff[of _ m, symmetric]
  1162       by (auto intro!: mult_mono simp add: pow2_add `0 < ?d` pow_d)
  1163     thus ?thesis
  1164       unfolding Float round_down.simps Let_def if_P[OF `0 < ?d`] .
  1165   next
  1166     case False
  1167     show ?thesis
  1168       unfolding Float round_down.simps Let_def if_not_P[OF False] .. 
  1169   qed
  1170 qed
  1171 
  1172 definition lb_mult :: "nat \<Rightarrow> float \<Rightarrow> float \<Rightarrow> float" where
  1173   "lb_mult prec x y =
  1174     (case normfloat (x * y) of Float m e \<Rightarrow>
  1175       let
  1176         l = bitlen m - int prec
  1177       in if l > 0 then Float (m div (2^nat l)) (e + l)
  1178                   else Float m e)"
  1179 
  1180 definition ub_mult :: "nat \<Rightarrow> float \<Rightarrow> float \<Rightarrow> float" where
  1181   "ub_mult prec x y =
  1182     (case normfloat (x * y) of Float m e \<Rightarrow>
  1183       let
  1184         l = bitlen m - int prec
  1185       in if l > 0 then Float (m div (2^nat l) + 1) (e + l)
  1186                   else Float m e)"
  1187 
  1188 lemma lb_mult: "real (lb_mult prec x y) \<le> real (x * y)"
  1189 proof (cases "normfloat (x * y)")
  1190   case (Float m e)
  1191   hence "odd m \<or> (m = 0 \<and> e = 0)" by (rule normfloat_imp_odd_or_zero)
  1192   let ?l = "bitlen m - int prec"
  1193   have "real (lb_mult prec x y) \<le> real (normfloat (x * y))"
  1194   proof (cases "?l > 0")
  1195     case False thus ?thesis unfolding lb_mult_def Float Let_def float.cases by auto
  1196   next
  1197     case True
  1198     have "real (m div 2^(nat ?l)) * pow2 ?l \<le> real m"
  1199     proof -
  1200       have "real (m div 2^(nat ?l)) * pow2 ?l = real (2^(nat ?l) * (m div 2^(nat ?l)))" unfolding real_of_int_mult real_of_int_power real_numeral unfolding pow2_int[symmetric] 
  1201         using `?l > 0` by auto
  1202       also have "\<dots> \<le> real (2^(nat ?l) * (m div 2^(nat ?l)) + m mod 2^(nat ?l))" unfolding real_of_int_add by auto
  1203       also have "\<dots> = real m" unfolding zmod_zdiv_equality[symmetric] ..
  1204       finally show ?thesis by auto
  1205     qed
  1206     thus ?thesis unfolding lb_mult_def Float Let_def float.cases if_P[OF True] real_of_float_simp pow2_add mult_commute mult_assoc by auto
  1207   qed
  1208   also have "\<dots> = real (x * y)" unfolding normfloat ..
  1209   finally show ?thesis .
  1210 qed
  1211 
  1212 lemma ub_mult: "real (x * y) \<le> real (ub_mult prec x y)"
  1213 proof (cases "normfloat (x * y)")
  1214   case (Float m e)
  1215   hence "odd m \<or> (m = 0 \<and> e = 0)" by (rule normfloat_imp_odd_or_zero)
  1216   let ?l = "bitlen m - int prec"
  1217   have "real (x * y) = real (normfloat (x * y))" unfolding normfloat ..
  1218   also have "\<dots> \<le> real (ub_mult prec x y)"
  1219   proof (cases "?l > 0")
  1220     case False thus ?thesis unfolding ub_mult_def Float Let_def float.cases by auto
  1221   next
  1222     case True
  1223     have "real m \<le> real (m div 2^(nat ?l) + 1) * pow2 ?l"
  1224     proof -
  1225       have "m mod 2^(nat ?l) < 2^(nat ?l)" by (rule pos_mod_bound) auto
  1226       hence mod_uneq: "real (m mod 2^(nat ?l)) \<le> 1 * 2^(nat ?l)" unfolding mult_1 real_of_int_less_iff[symmetric] by auto
  1227       
  1228       have "real m = real (2^(nat ?l) * (m div 2^(nat ?l)) + m mod 2^(nat ?l))" unfolding zmod_zdiv_equality[symmetric] ..
  1229       also have "\<dots> = real (m div 2^(nat ?l)) * 2^(nat ?l) + real (m mod 2^(nat ?l))" unfolding real_of_int_add by auto
  1230       also have "\<dots> \<le> (real (m div 2^(nat ?l)) + 1) * 2^(nat ?l)" unfolding left_distrib using mod_uneq by auto
  1231       finally show ?thesis unfolding pow2_int[symmetric] using True by auto
  1232     qed
  1233     thus ?thesis unfolding ub_mult_def Float Let_def float.cases if_P[OF True] real_of_float_simp pow2_add mult_commute mult_assoc by auto
  1234   qed
  1235   finally show ?thesis .
  1236 qed
  1237 
  1238 primrec float_abs :: "float \<Rightarrow> float" where
  1239   "float_abs (Float m e) = Float \<bar>m\<bar> e"
  1240 
  1241 instantiation float :: abs
  1242 begin
  1243 definition abs_float_def: "\<bar>x\<bar> = float_abs x"
  1244 instance ..
  1245 end
  1246 
  1247 lemma real_of_float_abs: "real \<bar>x :: float\<bar> = \<bar>real x\<bar>" 
  1248 proof (cases x)
  1249   case (Float m e)
  1250   have "\<bar>real m\<bar> * pow2 e = \<bar>real m * pow2 e\<bar>" unfolding abs_mult by auto
  1251   thus ?thesis unfolding Float abs_float_def float_abs.simps real_of_float_simp by auto
  1252 qed
  1253 
  1254 primrec floor_fl :: "float \<Rightarrow> float" where
  1255   "floor_fl (Float m e) = (if 0 \<le> e then Float m e
  1256                                   else Float (m div (2 ^ (nat (-e)))) 0)"
  1257 
  1258 lemma floor_fl: "real (floor_fl x) \<le> real x"
  1259 proof (cases x)
  1260   case (Float m e)
  1261   show ?thesis
  1262   proof (cases "0 \<le> e")
  1263     case False
  1264     hence me_eq: "pow2 (-e) = pow2 (int (nat (-e)))" by auto
  1265     have "real (Float (m div (2 ^ (nat (-e)))) 0) = real (m div 2 ^ (nat (-e)))" unfolding real_of_float_simp by auto
  1266     also have "\<dots> \<le> real m / real ((2::int) ^ (nat (-e)))" using real_of_int_div4 .
  1267     also have "\<dots> = real m * inverse (2 ^ (nat (-e)))" unfolding real_of_int_power real_numeral divide_inverse ..
  1268     also have "\<dots> = real (Float m e)" unfolding real_of_float_simp me_eq pow2_int pow2_neg[of e] ..
  1269     finally show ?thesis unfolding Float floor_fl.simps if_not_P[OF `\<not> 0 \<le> e`] .
  1270   next
  1271     case True thus ?thesis unfolding Float by auto
  1272   qed
  1273 qed
  1274 
  1275 lemma floor_pos_exp: assumes floor: "Float m e = floor_fl x" shows "0 \<le> e"
  1276 proof (cases x)
  1277   case (Float mx me)
  1278   from floor[unfolded Float floor_fl.simps] show ?thesis by (cases "0 \<le> me", auto)
  1279 qed
  1280 
  1281 declare floor_fl.simps[simp del]
  1282 
  1283 primrec ceiling_fl :: "float \<Rightarrow> float" where
  1284   "ceiling_fl (Float m e) = (if 0 \<le> e then Float m e
  1285                                     else Float (m div (2 ^ (nat (-e))) + 1) 0)"
  1286 
  1287 lemma ceiling_fl: "real x \<le> real (ceiling_fl x)"
  1288 proof (cases x)
  1289   case (Float m e)
  1290   show ?thesis
  1291   proof (cases "0 \<le> e")
  1292     case False
  1293     hence me_eq: "pow2 (-e) = pow2 (int (nat (-e)))" by auto
  1294     have "real (Float m e) = real m * inverse (2 ^ (nat (-e)))" unfolding real_of_float_simp me_eq pow2_int pow2_neg[of e] ..
  1295     also have "\<dots> = real m / real ((2::int) ^ (nat (-e)))" unfolding real_of_int_power real_numeral divide_inverse ..
  1296     also have "\<dots> \<le> 1 + real (m div 2 ^ (nat (-e)))" using real_of_int_div3[unfolded diff_le_eq] .
  1297     also have "\<dots> = real (Float (m div (2 ^ (nat (-e))) + 1) 0)" unfolding real_of_float_simp by auto
  1298     finally show ?thesis unfolding Float ceiling_fl.simps if_not_P[OF `\<not> 0 \<le> e`] .
  1299   next
  1300     case True thus ?thesis unfolding Float by auto
  1301   qed
  1302 qed
  1303 
  1304 declare ceiling_fl.simps[simp del]
  1305 
  1306 definition lb_mod :: "nat \<Rightarrow> float \<Rightarrow> float \<Rightarrow> float \<Rightarrow> float" where
  1307   "lb_mod prec x ub lb = x - ceiling_fl (float_divr prec x lb) * ub"
  1308 
  1309 definition ub_mod :: "nat \<Rightarrow> float \<Rightarrow> float \<Rightarrow> float \<Rightarrow> float" where
  1310   "ub_mod prec x ub lb = x - floor_fl (float_divl prec x ub) * lb"
  1311 
  1312 lemma lb_mod: fixes k :: int assumes "0 \<le> real x" and "real k * y \<le> real x" (is "?k * y \<le> ?x")
  1313   assumes "0 < real lb" "real lb \<le> y" (is "?lb \<le> y") "y \<le> real ub" (is "y \<le> ?ub")
  1314   shows "real (lb_mod prec x ub lb) \<le> ?x - ?k * y"
  1315 proof -
  1316   have "?lb \<le> ?ub" using assms by auto
  1317   have "0 \<le> ?lb" and "?lb \<noteq> 0" using assms by auto
  1318   have "?k * y \<le> ?x" using assms by auto
  1319   also have "\<dots> \<le> ?x / ?lb * ?ub" by (metis mult_left_mono[OF `?lb \<le> ?ub` `0 \<le> ?x`] divide_right_mono[OF _ `0 \<le> ?lb` ] times_divide_eq_left nonzero_mult_divide_cancel_right[OF `?lb \<noteq> 0`])
  1320   also have "\<dots> \<le> real (ceiling_fl (float_divr prec x lb)) * ?ub" by (metis mult_right_mono order_trans `0 \<le> ?lb` `?lb \<le> ?ub` float_divr ceiling_fl)
  1321   finally show ?thesis unfolding lb_mod_def real_of_float_sub real_of_float_mult by auto
  1322 qed
  1323 
  1324 lemma ub_mod:
  1325   fixes k :: int and x :: float
  1326   assumes "0 \<le> real x" and "real x \<le> real k * y" (is "?x \<le> ?k * y")
  1327   assumes "0 < real lb" "real lb \<le> y" (is "?lb \<le> y") "y \<le> real ub" (is "y \<le> ?ub")
  1328   shows "?x - ?k * y \<le> real (ub_mod prec x ub lb)"
  1329 proof -
  1330   have "?lb \<le> ?ub" using assms by auto
  1331   hence "0 \<le> ?lb" and "0 \<le> ?ub" and "?ub \<noteq> 0" using assms by auto
  1332   have "real (floor_fl (float_divl prec x ub)) * ?lb \<le> ?x / ?ub * ?lb" by (metis mult_right_mono order_trans `0 \<le> ?lb` `?lb \<le> ?ub` float_divl floor_fl)
  1333   also have "\<dots> \<le> ?x" by (metis mult_left_mono[OF `?lb \<le> ?ub` `0 \<le> ?x`] divide_right_mono[OF _ `0 \<le> ?ub` ] times_divide_eq_left nonzero_mult_divide_cancel_right[OF `?ub \<noteq> 0`])
  1334   also have "\<dots> \<le> ?k * y" using assms by auto
  1335   finally show ?thesis unfolding ub_mod_def real_of_float_sub real_of_float_mult by auto
  1336 qed
  1337 
  1338 lemma le_float_def'[code]: "f \<le> g = (case f - g of Float a b \<Rightarrow> a \<le> 0)"
  1339 proof -
  1340   have le_transfer: "(f \<le> g) = (real (f - g) \<le> 0)" by (auto simp add: le_float_def)
  1341   from float_split[of "f - g"] obtain a b where f_diff_g: "f - g = Float a b" by auto
  1342   with le_transfer have le_transfer': "f \<le> g = (real (Float a b) \<le> 0)" by simp
  1343   show ?thesis by (simp add: le_transfer' f_diff_g float_le_zero)
  1344 qed
  1345 
  1346 lemma less_float_def'[code]: "f < g = (case f - g of Float a b \<Rightarrow> a < 0)"
  1347 proof -
  1348   have less_transfer: "(f < g) = (real (f - g) < 0)" by (auto simp add: less_float_def)
  1349   from float_split[of "f - g"] obtain a b where f_diff_g: "f - g = Float a b" by auto
  1350   with less_transfer have less_transfer': "f < g = (real (Float a b) < 0)" by simp
  1351   show ?thesis by (simp add: less_transfer' f_diff_g float_less_zero)
  1352 qed
  1353 
  1354 end