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