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