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