# HG changeset patch
# User wenzelm
# Date 1415817056 -3600
# Node ID 87d4ce309e04b6ef1da6d519e8be33ff998c2b94
# Parent  b66788d19c0f0146b536cbab9237f42d1752419e# Parent  302104d8366b72e680aa8e859b16f9660ff80f06
merged

diff -r 302104d8366b -r 87d4ce309e04 NEWS
--- a/NEWS	Wed Nov 12 18:18:38 2014 +0100
+++ b/NEWS	Wed Nov 12 19:30:56 2014 +0100
@@ -168,6 +168,15 @@
 The "sos_cert" functionality is invoked as "sos" with additional
 argument. Minor INCOMPATIBILITY.
 
+* HOL-Decision_Procs:
+  - New counterexample generator quickcheck[approximation] for
+    inequalities of transcendental functions.
+    Uses hardware floating point arithmetic to randomly discover
+    potential counterexamples. Counterexamples are certified with the
+    'approximation' method.
+    See HOL/Decision_Procs/ex/Approximation_Quickcheck_Ex.thy for
+    examples.
+
 
 *** Document preparation ***
 
diff -r 302104d8366b -r 87d4ce309e04 src/HOL/Decision_Procs/Approximation.thy
--- a/src/HOL/Decision_Procs/Approximation.thy	Wed Nov 12 18:18:38 2014 +0100
+++ b/src/HOL/Decision_Procs/Approximation.thy	Wed Nov 12 19:30:56 2014 +0100
@@ -12,7 +12,6 @@
 keywords "approximate" :: diag
 begin
 
-declare powr_one [simp]
 declare powr_numeral [simp]
 declare powr_neg_one [simp]
 declare powr_neg_numeral [simp]
@@ -54,9 +53,13 @@
   fixes lb :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float" and ub :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float"
   assumes "0 \<le> real x" and f_Suc: "\<And>n. f (Suc n) = G ((F ^^ n) s) (f n)"
     and lb_0: "\<And> i k x. lb 0 i k x = 0"
-    and lb_Suc: "\<And> n i k x. lb (Suc n) i k x = lapprox_rat prec 1 k - x * (ub n (F i) (G i k) x)"
+    and lb_Suc: "\<And> n i k x. lb (Suc n) i k x = float_plus_down prec
+        (lapprox_rat prec 1 k)
+        (- float_round_up prec (x * (ub n (F i) (G i k) x)))"
     and ub_0: "\<And> i k x. ub 0 i k x = 0"
-    and ub_Suc: "\<And> n i k x. ub (Suc n) i k x = rapprox_rat prec 1 k - x * (lb n (F i) (G i k) x)"
+    and ub_Suc: "\<And> n i k x. ub (Suc n) i k x = float_plus_up prec
+        (rapprox_rat prec 1 k)
+        (- float_round_down prec (x * (lb n (F i) (G i k) x)))"
   shows "(lb n ((F ^^ j') s) (f j') x) \<le> horner F G n ((F ^^ j') s) (f j') x \<and>
          horner F G n ((F ^^ j') s) (f j') x \<le> (ub n ((F ^^ j') s) (f j') x)"
   (is "?lb n j' \<le> ?horner n j' \<and> ?horner n j' \<le> ?ub n j'")
@@ -67,7 +70,10 @@
   case (Suc n)
   thus ?case using lapprox_rat[of prec 1 "f j'"] using rapprox_rat[of 1 "f j'" prec]
     Suc[where j'="Suc j'"] `0 \<le> real x`
-    by (auto intro!: add_mono mult_left_mono simp add: lb_Suc ub_Suc field_simps f_Suc)
+    by (auto intro!: add_mono mult_left_mono float_round_down_le float_round_up_le
+      order_trans[OF add_mono[OF _ float_plus_down_le]]
+      order_trans[OF _ add_mono[OF _ float_plus_up_le]]
+      simp add: lb_Suc ub_Suc field_simps f_Suc)
 qed
 
 subsection "Theorems for floating point functions implementing the horner scheme"
@@ -83,11 +89,17 @@
   fixes F :: "nat \<Rightarrow> nat" and G :: "nat \<Rightarrow> nat \<Rightarrow> nat"
   assumes "0 \<le> real x" and f_Suc: "\<And>n. f (Suc n) = G ((F ^^ n) s) (f n)"
     and lb_0: "\<And> i k x. lb 0 i k x = 0"
-    and lb_Suc: "\<And> n i k x. lb (Suc n) i k x = lapprox_rat prec 1 k - x * (ub n (F i) (G i k) x)"
+    and lb_Suc: "\<And> n i k x. lb (Suc n) i k x = float_plus_down prec
+        (lapprox_rat prec 1 k)
+        (- float_round_up prec (x * (ub n (F i) (G i k) x)))"
     and ub_0: "\<And> i k x. ub 0 i k x = 0"
-    and ub_Suc: "\<And> n i k x. ub (Suc n) i k x = rapprox_rat prec 1 k - x * (lb n (F i) (G i k) x)"
-  shows "(lb n ((F ^^ j') s) (f j') x) \<le> (\<Sum>j=0..<n. (- 1) ^ j * (1 / (f (j' + j))) * (x ^ j))" (is "?lb") and
-    "(\<Sum>j=0..<n. (- 1) ^ j * (1 / (f (j' + j))) * (x ^ j)) \<le> (ub n ((F ^^ j') s) (f j') x)" (is "?ub")
+    and ub_Suc: "\<And> n i k x. ub (Suc n) i k x = float_plus_up prec
+        (rapprox_rat prec 1 k)
+        (- float_round_down prec (x * (lb n (F i) (G i k) x)))"
+  shows "(lb n ((F ^^ j') s) (f j') x) \<le> (\<Sum>j=0..<n. (- 1) ^ j * (1 / (f (j' + j))) * (x ^ j))"
+      (is "?lb")
+    and "(\<Sum>j=0..<n. (- 1) ^ j * (1 / (f (j' + j))) * (x ^ j)) \<le> (ub n ((F ^^ j') s) (f j') x)"
+      (is "?ub")
 proof -
   have "?lb  \<and> ?ub"
     using horner_bounds'[where lb=lb, OF `0 \<le> real x` f_Suc lb_0 lb_Suc ub_0 ub_Suc]
@@ -99,11 +111,15 @@
   fixes F :: "nat \<Rightarrow> nat" and G :: "nat \<Rightarrow> nat \<Rightarrow> nat"
   assumes "real x \<le> 0" and f_Suc: "\<And>n. f (Suc n) = G ((F ^^ n) s) (f n)"
     and lb_0: "\<And> i k x. lb 0 i k x = 0"
-    and lb_Suc: "\<And> n i k x. lb (Suc n) i k x = lapprox_rat prec 1 k + x * (ub n (F i) (G i k) x)"
+    and lb_Suc: "\<And> n i k x. lb (Suc n) i k x = float_plus_down prec
+        (lapprox_rat prec 1 k)
+        (float_round_down prec (x * (ub n (F i) (G i k) x)))"
     and ub_0: "\<And> i k x. ub 0 i k x = 0"
-    and ub_Suc: "\<And> n i k x. ub (Suc n) i k x = rapprox_rat prec 1 k + x * (lb n (F i) (G i k) x)"
-  shows "(lb n ((F ^^ j') s) (f j') x) \<le> (\<Sum>j=0..<n. (1 / (f (j' + j))) * real x ^ j)" (is "?lb") and
-    "(\<Sum>j=0..<n. (1 / (f (j' + j))) * real x ^ j) \<le> (ub n ((F ^^ j') s) (f j') x)" (is "?ub")
+    and ub_Suc: "\<And> n i k x. ub (Suc n) i k x = float_plus_up prec
+        (rapprox_rat prec 1 k)
+        (float_round_up prec (x * (lb n (F i) (G i k) x)))"
+  shows "(lb n ((F ^^ j') s) (f j') x) \<le> (\<Sum>j=0..<n. (1 / (f (j' + j))) * real x ^ j)" (is "?lb")
+    and "(\<Sum>j=0..<n. (1 / (f (j' + j))) * real x ^ j) \<le> (ub n ((F ^^ j') s) (f j') x)" (is "?ub")
 proof -
   { fix x y z :: float have "x - y * z = x + - y * z" by simp } note diff_mult_minus = this
   have sum_eq: "(\<Sum>j=0..<n. (1 / (f (j' + j))) * real x ^ j) =
@@ -111,10 +127,11 @@
     by (auto simp add: field_simps power_mult_distrib[symmetric])
   have "0 \<le> real (-x)" using assms by auto
   from horner_bounds[where G=G and F=F and f=f and s=s and prec=prec
-    and lb="\<lambda> n i k x. lb n i k (-x)" and ub="\<lambda> n i k x. ub n i k (-x)", unfolded lb_Suc ub_Suc diff_mult_minus,
-    OF this f_Suc lb_0 refl ub_0 refl]
+    and lb="\<lambda> n i k x. lb n i k (-x)" and ub="\<lambda> n i k x. ub n i k (-x)",
+    unfolded lb_Suc ub_Suc diff_mult_minus,
+    OF this f_Suc lb_0 _ ub_0 _]
   show "?lb" and "?ub" unfolding minus_minus sum_eq
-    by auto
+    by (auto simp: minus_float_round_up_eq minus_float_round_down_eq)
 qed
 
 subsection {* Selectors for next even or odd number *}
@@ -145,16 +162,29 @@
 
 section "Power function"
 
-definition float_power_bnds :: "nat \<Rightarrow> float \<Rightarrow> float \<Rightarrow> float * float" where
-"float_power_bnds n l u = (if odd n \<or> 0 < l then (l ^ n, u ^ n)
-                      else if u < 0         then (u ^ n, l ^ n)
-                                            else (0, (max (-l) u) ^ n))"
-
-lemma float_power_bnds: "(l1, u1) = float_power_bnds n l u \<Longrightarrow> x \<in> {l .. u} \<Longrightarrow> (x::real) ^ n \<in> {l1..u1}"
-  by (auto simp: float_power_bnds_def max_def split: split_if_asm
-           intro: power_mono_odd power_mono power_mono_even zero_le_even_power)
-
-lemma bnds_power: "\<forall> (x::real) l u. (l1, u1) = float_power_bnds n l u \<and> x \<in> {l .. u} \<longrightarrow> l1 \<le> x ^ n \<and> x ^ n \<le> u1"
+definition float_power_bnds :: "nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float \<Rightarrow> float * float" where
+"float_power_bnds prec n l u =
+  (if 0 < l then (power_down_fl prec l n, power_up_fl prec u n)
+  else if odd n then
+    (- power_up_fl prec (abs l) n,
+      if u < 0 then - power_down_fl prec (abs u) n else power_up_fl prec u n)
+  else if u < 0 then (power_down_fl prec (abs u) n, power_up_fl prec (abs l) n)
+  else (0, power_up_fl prec (max (abs l) (abs u)) n))"
+
+lemma le_minus_power_downI: "0 \<le> x \<Longrightarrow> x ^ n \<le> - a \<Longrightarrow> a \<le> - power_down prec x n"
+  by (subst le_minus_iff) (auto intro: power_down_le power_mono_odd)
+
+lemma float_power_bnds:
+  "(l1, u1) = float_power_bnds prec n l u \<Longrightarrow> x \<in> {l .. u} \<Longrightarrow> (x::real) ^ n \<in> {l1..u1}"
+  by (auto
+    simp: float_power_bnds_def max_def real_power_up_fl real_power_down_fl minus_le_iff
+    split: split_if_asm
+    intro!: power_up_le power_down_le le_minus_power_downI
+    intro: power_mono_odd power_mono power_mono_even zero_le_even_power)
+
+lemma bnds_power:
+  "\<forall> (x::real) l u. (l1, u1) = float_power_bnds prec n l u \<and> x \<in> {l .. u} \<longrightarrow>
+    l1 \<le> x ^ n \<and> x ^ n \<le> u1"
   using float_power_bnds by auto
 
 section "Square root"
@@ -169,13 +199,13 @@
 fun sqrt_iteration :: "nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float" where
 "sqrt_iteration prec 0 x = Float 1 ((bitlen \<bar>mantissa x\<bar> + exponent x) div 2 + 1)" |
 "sqrt_iteration prec (Suc m) x = (let y = sqrt_iteration prec m x
-                                  in Float 1 (- 1) * (y + float_divr prec x y))"
+                                  in Float 1 (- 1) * float_plus_up prec y (float_divr prec x y))"
 
 lemma compute_sqrt_iteration_base[code]:
   shows "sqrt_iteration prec n (Float m e) =
     (if n = 0 then Float 1 ((if m = 0 then 0 else bitlen \<bar>m\<bar> + e) div 2 + 1)
     else (let y = sqrt_iteration prec (n - 1) (Float m e) in
-      Float 1 (- 1) * (y + float_divr prec (Float m e) y)))"
+      Float 1 (- 1) * float_plus_up prec y (float_divr prec (Float m e) y)))"
   using bitlen_Float by (cases n) simp_all
 
 function ub_sqrt lb_sqrt :: "nat \<Rightarrow> float \<Rightarrow> float" where
@@ -220,7 +250,7 @@
       unfolding Float by (auto simp: powr_realpow[symmetric] field_simps powr_add)
     also have "\<dots> < 1 * 2 powr (e + nat (bitlen m))"
     proof (rule mult_strict_right_mono, auto)
-      show "real m < 2^nat (bitlen m)" using bitlen_bounds[OF `0 < m`, THEN conjunct2]
+      show "m < 2^nat (bitlen m)" using bitlen_bounds[OF `0 < m`, THEN conjunct2]
         unfolding real_of_int_less_iff[of m, symmetric] by auto
     qed
     finally have "sqrt x < sqrt (2 powr (e + bitlen m))" unfolding int_nat_bl by auto
@@ -264,8 +294,11 @@
   have "0 < sqrt x" using `0 < real x` by auto
   also have "\<dots> < real ?b" using Suc .
   finally have "sqrt x < (?b + x / ?b)/2" using sqrt_ub_pos_pos_1[OF Suc _ `0 < real x`] by auto
-  also have "\<dots> \<le> (?b + (float_divr prec x ?b))/2" by (rule divide_right_mono, auto simp add: float_divr)
+  also have "\<dots> \<le> (?b + (float_divr prec x ?b))/2"
+    by (rule divide_right_mono, auto simp add: float_divr)
   also have "\<dots> = (Float 1 (- 1)) * (?b + (float_divr prec x ?b))" by simp
+  also have "\<dots> \<le> (Float 1 (- 1)) * (float_plus_up prec ?b (float_divr prec x ?b))"
+    by (auto simp add: algebra_simps float_plus_up_le)
   finally show ?case unfolding sqrt_iteration.simps Let_def distrib_left .
 qed
 
@@ -352,31 +385,34 @@
 fun ub_arctan_horner :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float"
 and lb_arctan_horner :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float" where
   "ub_arctan_horner prec 0 k x = 0"
-| "ub_arctan_horner prec (Suc n) k x =
-    (rapprox_rat prec 1 k) - x * (lb_arctan_horner prec n (k + 2) x)"
+| "ub_arctan_horner prec (Suc n) k x = float_plus_up prec
+      (rapprox_rat prec 1 k) (- float_round_down prec (x * (lb_arctan_horner prec n (k + 2) x)))"
 | "lb_arctan_horner prec 0 k x = 0"
-| "lb_arctan_horner prec (Suc n) k x =
-    (lapprox_rat prec 1 k) - x * (ub_arctan_horner prec n (k + 2) x)"
+| "lb_arctan_horner prec (Suc n) k x = float_plus_down prec
+      (lapprox_rat prec 1 k) (- float_round_up prec (x * (ub_arctan_horner prec n (k + 2) x)))"
 
 lemma arctan_0_1_bounds':
-  assumes "0 \<le> real x" "real x \<le> 1" and "even n"
-  shows "arctan x \<in> {(x * lb_arctan_horner prec n 1 (x * x)) .. (x * ub_arctan_horner prec (Suc n) 1 (x * x))}"
+  assumes "0 \<le> real y" "real y \<le> 1" and "even n"
+  shows "arctan (sqrt y) \<in>
+      {(sqrt y * lb_arctan_horner prec n 1 y) .. (sqrt y * ub_arctan_horner prec (Suc n) 1 y)}"
 proof -
-  let ?c = "\<lambda>i. (- 1) ^ i * (1 / (i * 2 + (1::nat)) * real x ^ (i * 2 + 1))"
+  let ?c = "\<lambda>i. (- 1) ^ i * (1 / (i * 2 + (1::nat)) * sqrt y ^ (i * 2 + 1))"
   let ?S = "\<lambda>n. \<Sum> i=0..<n. ?c i"
 
-  have "0 \<le> real (x * x)" by auto
+  have "0 \<le> sqrt y" using assms by auto
+  have "sqrt y \<le> 1" using assms by auto
   from `even n` obtain m where "2 * m = n" by (blast elim: evenE)
 
-  have "arctan x \<in> { ?S n .. ?S (Suc n) }"
-  proof (cases "real x = 0")
+  have "arctan (sqrt y) \<in> { ?S n .. ?S (Suc n) }"
+  proof (cases "sqrt y = 0")
     case False
-    hence "0 < real x" using `0 \<le> real x` by auto
-    hence prem: "0 < 1 / (0 * 2 + (1::nat)) * real x ^ (0 * 2 + 1)" by auto
-
-    have "\<bar> real x \<bar> \<le> 1"  using `0 \<le> real x` `real x \<le> 1` by auto
-    from mp[OF summable_Leibniz(2)[OF zeroseq_arctan_series[OF this] monoseq_arctan_series[OF this]] prem, THEN spec, of m, unfolded `2 * m = n`]
-    show ?thesis unfolding arctan_series[OF `\<bar> real x \<bar> \<le> 1`] Suc_eq_plus1 atLeast0LessThan .
+    hence "0 < sqrt y" using `0 \<le> sqrt y` by auto
+    hence prem: "0 < 1 / (0 * 2 + (1::nat)) * sqrt y ^ (0 * 2 + 1)" by auto
+
+    have "\<bar> sqrt y \<bar> \<le> 1"  using `0 \<le> sqrt y` `sqrt y \<le> 1` by auto
+    from mp[OF summable_Leibniz(2)[OF zeroseq_arctan_series[OF this]
+      monoseq_arctan_series[OF this]] prem, THEN spec, of m, unfolded `2 * m = n`]
+    show ?thesis unfolding arctan_series[OF `\<bar> sqrt y \<bar> \<le> 1`] Suc_eq_plus1 atLeast0LessThan .
   qed auto
   note arctan_bounds = this[unfolded atLeastAtMost_iff]
 
@@ -385,111 +421,240 @@
   note bounds = horner_bounds[where s=1 and f="\<lambda>i. 2 * i + 1" and j'=0
     and lb="\<lambda>n i k x. lb_arctan_horner prec n k x"
     and ub="\<lambda>n i k x. ub_arctan_horner prec n k x",
-    OF `0 \<le> real (x*x)` F lb_arctan_horner.simps ub_arctan_horner.simps]
-
-  { have "(x * lb_arctan_horner prec n 1 (x*x)) \<le> ?S n"
-      using bounds(1) `0 \<le> real x`
+    OF `0 \<le> real y` F lb_arctan_horner.simps ub_arctan_horner.simps]
+
+  { have "(sqrt y * lb_arctan_horner prec n 1 y) \<le> ?S n"
+      using bounds(1) `0 \<le> sqrt y`
       unfolding power_add power_one_right mult.assoc[symmetric] setsum_left_distrib[symmetric]
-      unfolding mult.commute[where 'a=real] mult.commute[of _ "2::nat"] power_mult power2_eq_square[of "real x"]
+      unfolding mult.commute[where 'a=real] mult.commute[of _ "2::nat"] power_mult
       by (auto intro!: mult_left_mono)
-    also have "\<dots> \<le> arctan x" using arctan_bounds ..
-    finally have "(x * lb_arctan_horner prec n 1 (x*x)) \<le> arctan x" . }
+    also have "\<dots> \<le> arctan (sqrt y)" using arctan_bounds ..
+    finally have "(sqrt y * lb_arctan_horner prec n 1 y) \<le> arctan (sqrt y)" . }
   moreover
-  { have "arctan x \<le> ?S (Suc n)" using arctan_bounds ..
-    also have "\<dots> \<le> (x * ub_arctan_horner prec (Suc n) 1 (x*x))"
-      using bounds(2)[of "Suc n"] `0 \<le> real x`
+  { have "arctan (sqrt y) \<le> ?S (Suc n)" using arctan_bounds ..
+    also have "\<dots> \<le> (sqrt y * ub_arctan_horner prec (Suc n) 1 y)"
+      using bounds(2)[of "Suc n"] `0 \<le> sqrt y`
       unfolding power_add power_one_right mult.assoc[symmetric] setsum_left_distrib[symmetric]
-      unfolding mult.commute[where 'a=real] mult.commute[of _ "2::nat"] power_mult power2_eq_square[of "real x"]
+      unfolding mult.commute[where 'a=real] mult.commute[of _ "2::nat"] power_mult
       by (auto intro!: mult_left_mono)
-    finally have "arctan x \<le> (x * ub_arctan_horner prec (Suc n) 1 (x*x))" . }
+    finally have "arctan (sqrt y) \<le> (sqrt y * ub_arctan_horner prec (Suc n) 1 y)" . }
   ultimately show ?thesis by auto
 qed
 
-lemma arctan_0_1_bounds: assumes "0 \<le> real x" "real x \<le> 1"
-  shows "arctan x \<in> {(x * lb_arctan_horner prec (get_even n) 1 (x * x)) .. (x * ub_arctan_horner prec (get_odd n) 1 (x * x))}"
+lemma arctan_0_1_bounds: assumes "0 \<le> real y" "real y \<le> 1"
+  shows "arctan (sqrt y) \<in>
+    {(sqrt y * lb_arctan_horner prec (get_even n) 1 y) ..
+      (sqrt y * ub_arctan_horner prec (get_odd n) 1 y)}"
   using
     arctan_0_1_bounds'[OF assms, of n prec]
     arctan_0_1_bounds'[OF assms, of "n + 1" prec]
     arctan_0_1_bounds'[OF assms, of "n - 1" prec]
-  by (auto simp: get_even_def get_odd_def odd_pos simp del: ub_arctan_horner.simps lb_arctan_horner.simps)
+  by (auto simp: get_even_def get_odd_def odd_pos simp del: ub_arctan_horner.simps
+    lb_arctan_horner.simps)
+
+lemma arctan_lower_bound:
+  assumes "0 \<le> x"
+  shows "x / (1 + x\<^sup>2) \<le> arctan x" (is "?l x \<le> _")
+proof -
+  have "?l x - arctan x \<le> ?l 0 - arctan 0"
+    using assms
+    by (intro DERIV_nonpos_imp_nonincreasing[where f="\<lambda>x. ?l x - arctan x"])
+      (auto intro!: derivative_eq_intros simp: add_nonneg_eq_0_iff field_simps)
+  thus ?thesis by simp
+qed
+
+lemma arctan_divide_mono: "0 < x \<Longrightarrow> x \<le> y \<Longrightarrow> arctan y / y \<le> arctan x / x"
+  by (rule DERIV_nonpos_imp_nonincreasing[where f="\<lambda>x. arctan x / x"])
+    (auto intro!: derivative_eq_intros divide_nonpos_nonneg
+      simp: inverse_eq_divide arctan_lower_bound)
+
+lemma arctan_mult_mono: "0 \<le> x \<Longrightarrow> x \<le> y \<Longrightarrow> x * arctan y \<le> y * arctan x"
+  using arctan_divide_mono[of x y] by (cases "x = 0") (simp_all add: field_simps)
+
+lemma arctan_mult_le:
+  assumes "0 \<le> x" "x \<le> y" "y * z \<le> arctan y"
+  shows "x * z \<le> arctan x"
+proof cases
+  assume "x \<noteq> 0"
+  with assms have "z \<le> arctan y / y" by (simp add: field_simps)
+  also have "\<dots> \<le> arctan x / x" using assms `x \<noteq> 0` by (auto intro!: arctan_divide_mono)
+  finally show ?thesis using assms `x \<noteq> 0` by (simp add: field_simps)
+qed simp
+
+lemma arctan_le_mult:
+  assumes "0 < x" "x \<le> y" "arctan x \<le> x * z"
+  shows "arctan y \<le> y * z"
+proof -
+  from assms have "arctan y / y \<le> arctan x / x" by (auto intro!: arctan_divide_mono)
+  also have "\<dots> \<le> z" using assms by (auto simp: field_simps)
+  finally show ?thesis using assms by (simp add: field_simps)
+qed
+
+lemma arctan_0_1_bounds_le:
+  assumes "0 \<le> x" "x \<le> 1" "0 < real xl" "real xl \<le> x * x" "x * x \<le> real xu" "real xu \<le> 1"
+  shows "arctan x \<in>
+      {x * lb_arctan_horner p1 (get_even n) 1 xu .. x * ub_arctan_horner p2 (get_odd n) 1 xl}"
+proof -
+  from assms have "real xl \<le> 1" "sqrt (real xl) \<le> x" "x \<le> sqrt (real xu)" "0 \<le> real xu"
+    "0 \<le> real xl" "0 < sqrt (real xl)"
+    by (auto intro!: real_le_rsqrt real_le_lsqrt simp: power2_eq_square)
+  from arctan_0_1_bounds[OF `0 \<le> real xu`  `real xu \<le> 1`]
+  have "sqrt (real xu) * real (lb_arctan_horner p1 (get_even n) 1 xu) \<le> arctan (sqrt (real xu))"
+    by simp
+  from arctan_mult_le[OF `0 \<le> x` `x \<le> sqrt _`  this]
+  have "x * real (lb_arctan_horner p1 (get_even n) 1 xu) \<le> arctan x" .
+  moreover
+  from arctan_0_1_bounds[OF `0 \<le> real xl`  `real xl \<le> 1`]
+  have "arctan (sqrt (real xl)) \<le> sqrt (real xl) * real (ub_arctan_horner p2 (get_odd n) 1 xl)"
+    by simp
+  from arctan_le_mult[OF `0 < sqrt xl` `sqrt xl \<le> x` this]
+  have "arctan x \<le> x * real (ub_arctan_horner p2 (get_odd n) 1 xl)" .
+  ultimately show ?thesis by simp
+qed
+
+lemma mult_nonneg_le_one: fixes a::real assumes "0 \<le> a" "0 \<le> b" "a \<le> 1" "b \<le> 1" shows "a * b \<le> 1"
+proof -
+  have "a * b \<le> 1 * 1"
+    by (intro mult_mono assms) simp_all
+  thus ?thesis by simp
+qed
+
+lemma arctan_0_1_bounds_round:
+  assumes "0 \<le> real x" "real x \<le> 1"
+  shows "arctan x \<in>
+      {real x * lb_arctan_horner p1 (get_even n) 1 (float_round_up (Suc p2) (x * x)) ..
+        real x * ub_arctan_horner p3 (get_odd n) 1 (float_round_down (Suc p4) (x * x))}"
+  using assms
+  apply (cases "x > 0")
+   apply (intro arctan_0_1_bounds_le)
+   apply (auto simp: float_round_down.rep_eq float_round_up.rep_eq
+    intro!: truncate_up_le1 mult_nonneg_le_one truncate_down_le truncate_up_le truncate_down_pos
+      mult_pos_pos)
+  done
+
 
 subsection "Compute \<pi>"
 
 definition ub_pi :: "nat \<Rightarrow> float" where
-  "ub_pi prec = (let A = rapprox_rat prec 1 5 ;
-                     B = lapprox_rat prec 1 239
-                 in ((Float 1 2) * ((Float 1 2) * A * (ub_arctan_horner prec (get_odd (prec div 4 + 1)) 1 (A * A)) -
-                                                  B * (lb_arctan_horner prec (get_even (prec div 14 + 1)) 1 (B * B)))))"
+  "ub_pi prec =
+    (let
+      A = rapprox_rat prec 1 5 ;
+      B = lapprox_rat prec 1 239
+    in ((Float 1 2) * float_plus_up prec
+      ((Float 1 2) * float_round_up prec (A * (ub_arctan_horner prec (get_odd (prec div 4 + 1)) 1
+        (float_round_down (Suc prec) (A * A)))))
+      (- float_round_down prec (B * (lb_arctan_horner prec (get_even (prec div 14 + 1)) 1
+        (float_round_up (Suc prec) (B * B)))))))"
 
 definition lb_pi :: "nat \<Rightarrow> float" where
-  "lb_pi prec = (let A = lapprox_rat prec 1 5 ;
-                     B = rapprox_rat prec 1 239
-                 in ((Float 1 2) * ((Float 1 2) * A * (lb_arctan_horner prec (get_even (prec div 4 + 1)) 1 (A * A)) -
-                                                  B * (ub_arctan_horner prec (get_odd (prec div 14 + 1)) 1 (B * B)))))"
+  "lb_pi prec =
+    (let
+      A = lapprox_rat prec 1 5 ;
+      B = rapprox_rat prec 1 239
+    in ((Float 1 2) * float_plus_down prec
+      ((Float 1 2) * float_round_down prec (A * (lb_arctan_horner prec (get_even (prec div 4 + 1)) 1
+        (float_round_up (Suc prec) (A * A)))))
+      (- float_round_up prec (B * (ub_arctan_horner prec (get_odd (prec div 14 + 1)) 1
+        (float_round_down (Suc prec) (B * B)))))))"
 
 lemma pi_boundaries: "pi \<in> {(lb_pi n) .. (ub_pi n)}"
 proof -
-  have machin_pi: "pi = 4 * (4 * arctan (1 / 5) - arctan (1 / 239))" unfolding machin[symmetric] by auto
+  have machin_pi: "pi = 4 * (4 * arctan (1 / 5) - arctan (1 / 239))"
+    unfolding machin[symmetric] by auto
 
   { fix prec n :: nat fix k :: int assume "1 < k" hence "0 \<le> k" and "0 < k" and "1 \<le> k" by auto
     let ?k = "rapprox_rat prec 1 k"
+    let ?kl = "float_round_down (Suc prec) (?k * ?k)"
     have "1 div k = 0" using div_pos_pos_trivial[OF _ `1 < k`] by auto
 
-    have "0 \<le> real ?k" by (rule order_trans[OF _ rapprox_rat], auto simp add: `0 \<le> k`)
-    have "real ?k \<le> 1" 
-      by (rule rapprox_rat_le1, auto simp add: `0 < k` `1 \<le> k`)
-
+    have "0 \<le> real ?k" by (rule order_trans[OF _ rapprox_rat]) (auto simp add: `0 \<le> k`)
+    have "real ?k \<le> 1"
+      by (auto simp add: `0 < k` `1 \<le> k` less_imp_le
+        intro!: mult_nonneg_le_one order_trans[OF _ rapprox_rat] rapprox_rat_le1)
     have "1 / k \<le> ?k" using rapprox_rat[where x=1 and y=k] by auto
     hence "arctan (1 / k) \<le> arctan ?k" by (rule arctan_monotone')
-    also have "\<dots> \<le> (?k * ub_arctan_horner prec (get_odd n) 1 (?k * ?k))"
-      using arctan_0_1_bounds[OF `0 \<le> real ?k` `real ?k \<le> 1`] by auto
-    finally have "arctan (1 / k) \<le> ?k * ub_arctan_horner prec (get_odd n) 1 (?k * ?k)" .
+    also have "\<dots> \<le> (?k * ub_arctan_horner prec (get_odd n) 1 ?kl)"
+      using arctan_0_1_bounds_round[OF `0 \<le> real ?k` `real ?k \<le> 1`]
+      by auto
+    finally have "arctan (1 / k) \<le> ?k * ub_arctan_horner prec (get_odd n) 1 ?kl" .
   } note ub_arctan = this
 
   { fix prec n :: nat fix k :: int assume "1 < k" hence "0 \<le> k" and "0 < k" by auto
     let ?k = "lapprox_rat prec 1 k"
+    let ?ku = "float_round_up (Suc prec) (?k * ?k)"
     have "1 div k = 0" using div_pos_pos_trivial[OF _ `1 < k`] by auto
     have "1 / k \<le> 1" using `1 < k` by auto
-    have "\<And>n. 0 \<le> real ?k" using lapprox_rat_nonneg[where x=1 and y=k, OF zero_le_one `0 < k`] by (auto simp add: `1 div k = 0`)
-    have "\<And>n. real ?k \<le> 1" using lapprox_rat by (rule order_trans, auto simp add: `1 / k \<le> 1`)
+    have "0 \<le> real ?k" using lapprox_rat_nonneg[where x=1 and y=k, OF zero_le_one `0 \<le> k`]
+      by (auto simp add: `1 div k = 0`)
+    have "0 \<le> real (?k * ?k)" by simp
+    have "real ?k \<le> 1" using lapprox_rat by (rule order_trans, auto simp add: `1 / k \<le> 1`)
+    hence "real (?k * ?k) \<le> 1" using `0 \<le> real ?k` by (auto intro!: mult_nonneg_le_one)
 
     have "?k \<le> 1 / k" using lapprox_rat[where x=1 and y=k] by auto
 
-    have "?k * lb_arctan_horner prec (get_even n) 1 (?k * ?k) \<le> arctan ?k"
-      using arctan_0_1_bounds[OF `0 \<le> real ?k` `real ?k \<le> 1`] by auto
+    have "?k * lb_arctan_horner prec (get_even n) 1 ?ku \<le> arctan ?k"
+      using arctan_0_1_bounds_round[OF `0 \<le> real ?k` `real ?k \<le> 1`]
+      by auto
     also have "\<dots> \<le> arctan (1 / k)" using `?k \<le> 1 / k` by (rule arctan_monotone')
-    finally have "?k * lb_arctan_horner prec (get_even n) 1 (?k * ?k) \<le> arctan (1 / k)" .
+    finally have "?k * lb_arctan_horner prec (get_even n) 1 ?ku \<le> arctan (1 / k)" .
   } note lb_arctan = this
 
-  have "pi \<le> ub_pi n \<and> lb_pi n \<le> pi"
-    unfolding lb_pi_def ub_pi_def machin_pi Let_def unfolding Float_num
-    using lb_arctan[of 5] ub_arctan[of 239] lb_arctan[of 239] ub_arctan[of 5] powr_realpow[of 2 2]
-    by (auto intro!: mult_left_mono add_mono simp add: uminus_add_conv_diff [symmetric] simp del: uminus_add_conv_diff)
-  then show ?thesis by auto
+  have "pi \<le> ub_pi n "
+    unfolding ub_pi_def machin_pi Let_def times_float.rep_eq Float_num
+    using lb_arctan[of 239] ub_arctan[of 5] powr_realpow[of 2 2]
+    by (intro mult_left_mono float_plus_up_le float_plus_down_le)
+      (auto intro!: mult_left_mono float_round_down_le float_round_up_le diff_mono)
+  moreover have "lb_pi n \<le> pi"
+    unfolding lb_pi_def machin_pi Let_def times_float.rep_eq Float_num
+    using lb_arctan[of 5] ub_arctan[of 239]
+    by (intro mult_left_mono float_plus_up_le float_plus_down_le)
+      (auto intro!: mult_left_mono float_round_down_le float_round_up_le diff_mono)
+  ultimately show ?thesis by auto
 qed
 
 subsection "Compute arcus tangens in the entire domain"
 
 function lb_arctan :: "nat \<Rightarrow> float \<Rightarrow> float" and ub_arctan :: "nat \<Rightarrow> float \<Rightarrow> float" where
-  "lb_arctan prec x = (let ub_horner = \<lambda> x. x * ub_arctan_horner prec (get_odd (prec div 4 + 1)) 1 (x * x) ;
-                           lb_horner = \<lambda> x. x * lb_arctan_horner prec (get_even (prec div 4 + 1)) 1 (x * x)
-    in (if x < 0          then - ub_arctan prec (-x) else
-        if x \<le> Float 1 (- 1) then lb_horner x else
-        if x \<le> Float 1 1  then Float 1 1 * lb_horner (float_divl prec x (1 + ub_sqrt prec (1 + x * x)))
-                          else (let inv = float_divr prec 1 x
-                                in if inv > 1 then 0
-                                              else lb_pi prec * Float 1 (- 1) - ub_horner inv)))"
-
-| "ub_arctan prec x = (let lb_horner = \<lambda> x. x * lb_arctan_horner prec (get_even (prec div 4 + 1)) 1 (x * x) ;
-                           ub_horner = \<lambda> x. x * ub_arctan_horner prec (get_odd (prec div 4 + 1)) 1 (x * x)
-    in (if x < 0          then - lb_arctan prec (-x) else
-        if x \<le> Float 1 (- 1) then ub_horner x else
-        if x \<le> Float 1 1  then let y = float_divr prec x (1 + lb_sqrt prec (1 + x * x))
-                               in if y > 1 then ub_pi prec * Float 1 (- 1)
-                                           else Float 1 1 * ub_horner y
-                          else ub_pi prec * Float 1 (- 1) - lb_horner (float_divl prec 1 x)))"
+  "lb_arctan prec x =
+    (let
+      ub_horner = \<lambda> x. float_round_up prec
+        (x *
+          ub_arctan_horner prec (get_odd (prec div 4 + 1)) 1 (float_round_down (Suc prec) (x * x)));
+      lb_horner = \<lambda> x. float_round_down prec
+        (x *
+          lb_arctan_horner prec (get_even (prec div 4 + 1)) 1 (float_round_up (Suc prec) (x * x)))
+    in
+      if x < 0 then - ub_arctan prec (-x)
+      else if x \<le> Float 1 (- 1) then lb_horner x
+      else if x \<le> Float 1 1 then
+        Float 1 1 *
+        lb_horner
+          (float_divl prec x
+            (float_plus_up prec 1
+              (ub_sqrt prec (float_plus_up prec 1 (float_round_up prec (x * x))))))
+      else let inv = float_divr prec 1 x in
+        if inv > 1 then 0
+        else float_plus_down prec (lb_pi prec * Float 1 (- 1)) ( - ub_horner inv))"
+
+| "ub_arctan prec x =
+    (let
+      lb_horner = \<lambda> x. float_round_down prec
+        (x *
+          lb_arctan_horner prec (get_even (prec div 4 + 1)) 1 (float_round_up (Suc prec) (x * x))) ;
+      ub_horner = \<lambda> x. float_round_up prec
+        (x *
+          ub_arctan_horner prec (get_odd (prec div 4 + 1)) 1 (float_round_down (Suc prec) (x * x)))
+    in if x < 0 then - lb_arctan prec (-x)
+    else if x \<le> Float 1 (- 1) then ub_horner x
+    else if x \<le> Float 1 1 then
+      let y = float_divr prec x
+        (float_plus_down
+          (Suc prec) 1 (lb_sqrt prec (float_plus_down prec 1 (float_round_down prec (x * x)))))
+      in if y > 1 then ub_pi prec * Float 1 (- 1) else Float 1 1 * ub_horner y
+    else float_plus_up prec (ub_pi prec * Float 1 (- 1)) ( - lb_horner (float_divl prec 1 x)))"
 by pat_completeness auto
-termination by (relation "measure (\<lambda> v. let (prec, x) = case_sum id id v in (if x < 0 then 1 else 0))", auto)
+termination
+by (relation "measure (\<lambda> v. let (prec, x) = case_sum id id v in (if x < 0 then 1 else 0))", auto)
 
 declare ub_arctan_horner.simps[simp del]
 declare lb_arctan_horner.simps[simp del]
@@ -497,31 +662,40 @@
 lemma lb_arctan_bound': assumes "0 \<le> real x"
   shows "lb_arctan prec x \<le> arctan x"
 proof -
-  have "\<not> x < 0" and "0 \<le> x" using `0 \<le> real x` by auto
-  let "?ub_horner x" = "x * ub_arctan_horner prec (get_odd (prec div 4 + 1)) 1 (x * x)"
-    and "?lb_horner x" = "x * lb_arctan_horner prec (get_even (prec div 4 + 1)) 1 (x * x)"
+  have "\<not> x < 0" and "0 \<le> x"
+    using `0 \<le> real x` by (auto intro!: truncate_up_le )
+
+  let "?ub_horner x" =
+      "x * ub_arctan_horner prec (get_odd (prec div 4 + 1)) 1 (float_round_down (Suc prec) (x * x))"
+    and "?lb_horner x" =
+      "x * lb_arctan_horner prec (get_even (prec div 4 + 1)) 1 (float_round_up (Suc prec) (x * x))"
 
   show ?thesis
   proof (cases "x \<le> Float 1 (- 1)")
-    case True hence "real x \<le> 1" by auto
-    show ?thesis unfolding lb_arctan.simps Let_def if_not_P[OF `\<not> x < 0`] if_P[OF True]
-      using arctan_0_1_bounds[OF `0 \<le> real x` `real x \<le> 1`] by auto
+    case True hence "real x \<le> 1" by simp
+    from arctan_0_1_bounds_round[OF `0 \<le> real x` `real x \<le> 1`]
+    show ?thesis
+      unfolding lb_arctan.simps Let_def if_not_P[OF `\<not> x < 0`] if_P[OF True] using `0 \<le> x`
+      by (auto intro!: float_round_down_le)
   next
     case False hence "0 < real x" by auto
     let ?R = "1 + sqrt (1 + real x * real x)"
-    let ?fR = "1 + ub_sqrt prec (1 + x * x)"
+    let ?sxx = "float_plus_up prec 1 (float_round_up prec (x * x))"
+    let ?fR = "float_plus_up prec 1 (ub_sqrt prec ?sxx)"
     let ?DIV = "float_divl prec x ?fR"
 
-    have sqr_ge0: "0 \<le> 1 + real x * real x" using sum_power2_ge_zero[of 1 "real x", unfolded numeral_2_eq_2] by auto
-    hence divisor_gt0: "0 < ?R" by (auto intro: add_pos_nonneg)
-
-    have "sqrt (1 + x * x) \<le> ub_sqrt prec (1 + x * x)"
-      using bnds_sqrt'[of "1 + x * x"] by auto
-
-    hence "?R \<le> ?fR" by auto
+    have divisor_gt0: "0 < ?R" by (auto intro: add_pos_nonneg)
+
+    have "sqrt (1 + x*x) \<le> sqrt ?sxx"
+      by (auto simp: float_plus_up.rep_eq plus_up_def float_round_up.rep_eq intro!: truncate_up_le)
+    also have "\<dots> \<le> ub_sqrt prec ?sxx"
+      using bnds_sqrt'[of ?sxx prec] by auto
+    finally
+    have "sqrt (1 + x*x) \<le> ub_sqrt prec ?sxx" .
+    hence "?R \<le> ?fR" by (auto simp: float_plus_up.rep_eq plus_up_def intro!: truncate_up_le)
     hence "0 < ?fR" and "0 < real ?fR" using `0 < ?R` by auto
 
-    have monotone: "(float_divl prec x ?fR) \<le> x / ?R"
+    have monotone: "?DIV \<le> x / ?R"
     proof -
       have "?DIV \<le> real x / ?fR" by (rule float_divl)
       also have "\<dots> \<le> x / ?R" by (rule divide_left_mono[OF `?R \<le> ?fR` `0 \<le> real x` mult_pos_pos[OF order_less_le_trans[OF divisor_gt0 `?R \<le> real ?fR`] divisor_gt0]])
@@ -533,20 +707,20 @@
       case True
 
       have "x \<le> sqrt (1 + x * x)" using real_sqrt_sum_squares_ge2[where x=1, unfolded numeral_2_eq_2] by auto
-      also have "\<dots> \<le> (ub_sqrt prec (1 + x * x))"
-        using bnds_sqrt'[of "1 + x * x"] by auto
-      finally have "real x \<le> ?fR" by auto
+      also note `\<dots> \<le> (ub_sqrt prec ?sxx)`
+      finally have "real x \<le> ?fR" by (auto simp: float_plus_up.rep_eq plus_up_def intro!: truncate_up_le)
       moreover have "?DIV \<le> real x / ?fR" by (rule float_divl)
       ultimately have "real ?DIV \<le> 1" unfolding divide_le_eq_1_pos[OF `0 < real ?fR`, symmetric] by auto
 
       have "0 \<le> real ?DIV" using float_divl_lower_bound[OF `0 \<le> x`] `0 < ?fR` unfolding less_eq_float_def by auto
 
-      have "(Float 1 1 * ?lb_horner ?DIV) \<le> 2 * arctan (float_divl prec x ?fR)"
-        using arctan_0_1_bounds[OF `0 \<le> real ?DIV` `real ?DIV \<le> 1`] by auto
+      from arctan_0_1_bounds_round[OF `0 \<le> real (?DIV)` `real (?DIV) \<le> 1`]
+      have "Float 1 1 * ?lb_horner ?DIV \<le> 2 * arctan ?DIV" by simp
       also have "\<dots> \<le> 2 * arctan (x / ?R)"
-        using arctan_monotone'[OF monotone] by (auto intro!: mult_left_mono)
+        using arctan_monotone'[OF monotone] by (auto intro!: mult_left_mono arctan_monotone')
       also have "2 * arctan (x / ?R) = arctan x" using arctan_half[symmetric] unfolding numeral_2_eq_2 power_Suc2 power_0 mult_1_left .
-      finally show ?thesis unfolding lb_arctan.simps Let_def if_not_P[OF `\<not> x < 0`] if_not_P[OF `\<not> x \<le> Float 1 (- 1)`] if_P[OF True] .
+      finally show ?thesis unfolding lb_arctan.simps Let_def if_not_P[OF `\<not> x < 0`] if_not_P[OF `\<not> x \<le> Float 1 (- 1)`] if_P[OF True]
+        by (auto simp: float_round_down.rep_eq intro!: order_trans[OF mult_left_mono[OF truncate_down]])
     next
       case False
       hence "2 < real x" by auto
@@ -568,8 +742,10 @@
         have "1 / x \<noteq> 0" and "0 < 1 / x" using `0 < real x` by auto
 
         have "arctan (1 / x) \<le> arctan ?invx" unfolding one_float.rep_eq[symmetric] by (rule arctan_monotone', rule float_divr)
-        also have "\<dots> \<le> (?ub_horner ?invx)" using arctan_0_1_bounds[OF `0 \<le> real ?invx` `real ?invx \<le> 1`] by auto
-        finally have "pi / 2 - (?ub_horner ?invx) \<le> arctan x"
+        also have "\<dots> \<le> ?ub_horner ?invx" using arctan_0_1_bounds_round[OF `0 \<le> real ?invx` `real ?invx \<le> 1`]
+          by (auto intro!: float_round_up_le)
+        also note float_round_up
+        finally have "pi / 2 - float_round_up prec (?ub_horner ?invx) \<le> arctan x"
           using `0 \<le> arctan x` arctan_inverse[OF `1 / x \<noteq> 0`]
           unfolding real_sgn_pos[OF `0 < 1 / real x`] le_diff_eq by auto
         moreover
@@ -577,7 +753,7 @@
           unfolding Float_num times_divide_eq_right mult_1_left using pi_boundaries by simp
         ultimately
         show ?thesis unfolding lb_arctan.simps Let_def if_not_P[OF `\<not> x < 0`] if_not_P[OF `\<not> x \<le> Float 1 (- 1)`] if_not_P[OF `\<not> x \<le> Float 1 1`] if_not_P[OF False]
-          by auto
+          by (auto intro!: float_plus_down_le)
       qed
     qed
   qed
@@ -588,18 +764,20 @@
 proof -
   have "\<not> x < 0" and "0 \<le> x" using `0 \<le> real x` by auto
 
-  let "?ub_horner x" = "x * ub_arctan_horner prec (get_odd (prec div 4 + 1)) 1 (x * x)"
-    and "?lb_horner x" = "x * lb_arctan_horner prec (get_even (prec div 4 + 1)) 1 (x * x)"
+  let "?ub_horner x" = "float_round_up prec (x * ub_arctan_horner prec (get_odd (prec div 4 + 1)) 1 (float_round_down (Suc prec) (x * x)))"
+    and "?lb_horner x" = "float_round_down prec (x * lb_arctan_horner prec (get_even (prec div 4 + 1)) 1 (float_round_up (Suc prec) (x * x)))"
 
   show ?thesis
   proof (cases "x \<le> Float 1 (- 1)")
     case True hence "real x \<le> 1" by auto
     show ?thesis unfolding ub_arctan.simps Let_def if_not_P[OF `\<not> x < 0`] if_P[OF True]
-      using arctan_0_1_bounds[OF `0 \<le> real x` `real x \<le> 1`] by auto
+      using arctan_0_1_bounds_round[OF `0 \<le> real x` `real x \<le> 1`]
+      by (auto intro!: float_round_up_le)
   next
     case False hence "0 < real x" by auto
     let ?R = "1 + sqrt (1 + real x * real x)"
-    let ?fR = "1 + lb_sqrt prec (1 + x * x)"
+    let ?sxx = "float_plus_down prec 1 (float_round_down prec (x * x))"
+    let ?fR = "float_plus_down (Suc prec) 1 (lb_sqrt prec ?sxx)"
     let ?DIV = "float_divr prec x ?fR"
 
     have sqr_ge0: "0 \<le> 1 + real x * real x" using sum_power2_ge_zero[of 1 "real x", unfolded numeral_2_eq_2] by auto
@@ -607,11 +785,16 @@
 
     hence divisor_gt0: "0 < ?R" by (auto intro: add_pos_nonneg)
 
-    have "lb_sqrt prec (1 + x * x) \<le> sqrt (1 + x * x)"
-      using bnds_sqrt'[of "1 + x * x"] by auto
-    hence "?fR \<le> ?R" by auto
-    have "0 < real ?fR" by (rule order_less_le_trans[OF zero_less_one], auto simp add: lb_sqrt_lower_bound[OF `0 \<le> real (1 + x*x)`])
-
+    have "lb_sqrt prec ?sxx \<le> sqrt ?sxx"
+      using bnds_sqrt'[of ?sxx] by auto
+    also have "\<dots> \<le> sqrt (1 + x*x)"
+      by (auto simp: float_plus_down.rep_eq plus_down_def float_round_down.rep_eq truncate_down_le)
+    finally have "lb_sqrt prec ?sxx \<le> sqrt (1 + x*x)" .
+    hence "?fR \<le> ?R" by (auto simp: float_plus_down.rep_eq plus_down_def truncate_down_le)
+    have "0 < real ?fR"
+      by (auto simp: float_plus_down.rep_eq plus_down_def float_round_down.rep_eq
+        intro!: truncate_down_ge1 lb_sqrt_lower_bound order_less_le_trans[OF zero_less_one]
+        truncate_down_nonneg add_nonneg_nonneg)
     have monotone: "x / ?R \<le> (float_divr prec x ?fR)"
     proof -
       from divide_left_mono[OF `?fR \<le> ?R` `0 \<le> real x` mult_pos_pos[OF divisor_gt0 `0 < real ?fR`]]
@@ -640,7 +823,8 @@
         also have "\<dots> \<le> 2 * arctan (?DIV)"
           using arctan_monotone'[OF monotone] by (auto intro!: mult_left_mono)
         also have "\<dots> \<le> (Float 1 1 * ?ub_horner ?DIV)" unfolding Float_num
-          using arctan_0_1_bounds[OF `0 \<le> real ?DIV` `real ?DIV \<le> 1`] by auto
+          using arctan_0_1_bounds_round[OF `0 \<le> real ?DIV` `real ?DIV \<le> 1`]
+          by (auto intro!: float_round_up_le)
         finally show ?thesis unfolding ub_arctan.simps Let_def if_not_P[OF `\<not> x < 0`] if_not_P[OF `\<not> x \<le> Float 1 (- 1)`] if_P[OF `x \<le> Float 1 1`] if_not_P[OF False] .
       qed
     next
@@ -658,7 +842,8 @@
 
       have "1 / x \<noteq> 0" and "0 < 1 / x" using `0 < real x` by auto
 
-      have "(?lb_horner ?invx) \<le> arctan (?invx)" using arctan_0_1_bounds[OF `0 \<le> real ?invx` `real ?invx \<le> 1`] by auto
+      have "(?lb_horner ?invx) \<le> arctan (?invx)" using arctan_0_1_bounds_round[OF `0 \<le> real ?invx` `real ?invx \<le> 1`]
+        by (auto intro!: float_round_down_le)
       also have "\<dots> \<le> arctan (1 / x)" unfolding one_float.rep_eq[symmetric] by (rule arctan_monotone', rule float_divl)
       finally have "arctan x \<le> pi / 2 - (?lb_horner ?invx)"
         using `0 \<le> arctan x` arctan_inverse[OF `1 / x \<noteq> 0`]
@@ -667,7 +852,7 @@
       have "pi / 2 \<le> ub_pi prec * Float 1 (- 1)" unfolding Float_num times_divide_eq_right mult_1_right using pi_boundaries by auto
       ultimately
       show ?thesis unfolding ub_arctan.simps Let_def if_not_P[OF `\<not> x < 0`]if_not_P[OF `\<not> x \<le> Float 1 (- 1)`] if_not_P[OF False]
-        by auto
+        by (auto intro!: float_round_up_le float_plus_up_le)
     qed
   qed
 qed
@@ -711,11 +896,13 @@
 fun ub_sin_cos_aux :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float"
 and lb_sin_cos_aux :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float" where
   "ub_sin_cos_aux prec 0 i k x = 0"
-| "ub_sin_cos_aux prec (Suc n) i k x =
-    (rapprox_rat prec 1 k) - x * (lb_sin_cos_aux prec n (i + 2) (k * i * (i + 1)) x)"
+| "ub_sin_cos_aux prec (Suc n) i k x = float_plus_up prec
+    (rapprox_rat prec 1 k) (-
+      float_round_down prec (x * (lb_sin_cos_aux prec n (i + 2) (k * i * (i + 1)) x)))"
 | "lb_sin_cos_aux prec 0 i k x = 0"
-| "lb_sin_cos_aux prec (Suc n) i k x =
-    (lapprox_rat prec 1 k) - x * (ub_sin_cos_aux prec n (i + 2) (k * i * (i + 1)) x)"
+| "lb_sin_cos_aux prec (Suc n) i k x = float_plus_down prec
+    (lapprox_rat prec 1 k) (-
+      float_round_up prec (x * (ub_sin_cos_aux prec n (i + 2) (k * i * (i + 1)) x)))"
 
 lemma cos_aux:
   shows "(lb_sin_cos_aux prec n 1 1 (x * x)) \<le> (\<Sum> i=0..<n. (- 1) ^ i * (1/real (fact (2 * i))) * x ^(2 * i))" (is "?lb")
@@ -734,6 +921,13 @@
   show "?lb" and "?ub" by (auto simp add: power_mult power2_eq_square[of "real x"])
 qed
 
+lemma lb_sin_cos_aux_zero_le_one: "lb_sin_cos_aux prec n i j 0 \<le> 1"
+  by (cases j n rule: nat.exhaust[case_product nat.exhaust])
+    (auto intro!: float_plus_down_le order_trans[OF lapprox_rat])
+
+lemma one_le_ub_sin_cos_aux: "odd n \<Longrightarrow> 1 \<le> ub_sin_cos_aux prec n i (Suc 0) 0"
+  by (cases n) (auto intro!: float_plus_up_le order_trans[OF _ rapprox_rat])
+
 lemma cos_boundaries: assumes "0 \<le> real x" and "x \<le> pi / 2"
   shows "cos x \<in> {(lb_sin_cos_aux prec (get_even n) 1 1 (x * x)) .. (ub_sin_cos_aux prec (get_odd n) 1 1 (x * x))}"
 proof (cases "real x = 0")
@@ -818,16 +1012,11 @@
   ultimately show ?thesis by auto
 next
   case True
-  show ?thesis
-  proof (cases "n = 0")
-    case True
-    thus ?thesis unfolding `n = 0` get_even_def get_odd_def
-      using `real x = 0` lapprox_rat[where x="-1" and y=1]
-      by (auto simp: Float.compute_lapprox_rat Float.compute_rapprox_rat)
-  next
-    case False with not0_implies_Suc obtain m where "n = Suc m" by blast
-    thus ?thesis unfolding `n = Suc m` get_even_def get_odd_def using `real x = 0` rapprox_rat[where x=1 and y=1] lapprox_rat[where x=1 and y=1] by (cases "even (Suc m)", auto)
-  qed
+  hence "x = 0"
+    by transfer
+  thus ?thesis
+    using lb_sin_cos_aux_zero_le_one one_le_ub_sin_cos_aux
+    by simp
 qed
 
 lemma sin_aux: assumes "0 \<le> real x"
@@ -947,7 +1136,7 @@
 definition lb_cos :: "nat \<Rightarrow> float \<Rightarrow> float" where
 "lb_cos prec x = (let
     horner = \<lambda> x. lb_sin_cos_aux prec (get_even (prec div 4 + 1)) 1 1 (x * x) ;
-    half = \<lambda> x. if x < 0 then - 1 else Float 1 1 * x * x - 1
+    half = \<lambda> x. if x < 0 then - 1 else float_plus_down prec (Float 1 1 * x * x) (- 1)
   in if x < Float 1 (- 1) then horner x
 else if x < 1          then half (horner (x * Float 1 (- 1)))
                        else half (half (horner (x * Float 1 (- 2)))))"
@@ -955,7 +1144,7 @@
 definition ub_cos :: "nat \<Rightarrow> float \<Rightarrow> float" where
 "ub_cos prec x = (let
     horner = \<lambda> x. ub_sin_cos_aux prec (get_odd (prec div 4 + 1)) 1 1 (x * x) ;
-    half = \<lambda> x. Float 1 1 * x * x - 1
+    half = \<lambda> x. float_plus_up prec (Float 1 1 * x * x) (- 1)
   in if x < Float 1 (- 1) then horner x
 else if x < 1          then half (horner (x * Float 1 (- 1)))
                        else half (half (horner (x * Float 1 (- 2)))))"
@@ -974,8 +1163,8 @@
   have "\<not> x < 0" using `0 \<le> real x` by auto
   let "?ub_horner x" = "ub_sin_cos_aux prec (get_odd (prec div 4 + 1)) 1 1 (x * x)"
   let "?lb_horner x" = "lb_sin_cos_aux prec (get_even (prec div 4 + 1)) 1 1 (x * x)"
-  let "?ub_half x" = "Float 1 1 * x * x - 1"
-  let "?lb_half x" = "if x < 0 then - 1 else Float 1 1 * x * x - 1"
+  let "?ub_half x" = "float_plus_up prec (Float 1 1 * x * x) (- 1)"
+  let "?lb_half x" = "if x < 0 then - 1 else float_plus_down prec (Float 1 1 * x * x) (- 1)"
 
   show ?thesis
   proof (cases "x < Float 1 (- 1)")
@@ -999,7 +1188,8 @@
         have "real y * real y \<le> cos ?x2 * cos ?x2" .
         hence "2 * real y * real y \<le> 2 * cos ?x2 * cos ?x2" by auto
         hence "2 * real y * real y - 1 \<le> 2 * cos (x / 2) * cos (x / 2) - 1" unfolding Float_num by auto
-        thus ?thesis unfolding if_not_P[OF False] x_half Float_num by auto
+        thus ?thesis unfolding if_not_P[OF False] x_half Float_num
+          by (auto intro!: float_plus_down_le)
       qed
     } note lb_half = this
 
@@ -1015,7 +1205,8 @@
         have "cos ?x2 * cos ?x2 \<le> real y * real y" .
         hence "2 * cos ?x2 * cos ?x2 \<le> 2 * real y * real y" by auto
         hence "2 * cos (x / 2) * cos (x / 2) - 1 \<le> 2 * real y * real y - 1" unfolding Float_num by auto
-        thus ?thesis unfolding x_half Float_num by auto
+        thus ?thesis unfolding x_half Float_num
+          by (auto intro!: float_plus_up_le)
       qed
     } note ub_half = this
 
@@ -1079,8 +1270,8 @@
     lpi = float_round_down prec (lb_pi prec) ;
     upi = float_round_up prec (ub_pi prec) ;
     k = floor_fl (float_divr prec (lx + lpi) (2 * lpi)) ;
-    lx = lx - k * 2 * (if k < 0 then lpi else upi) ;
-    ux = ux - k * 2 * (if k < 0 then upi else lpi)
+    lx = float_plus_down prec lx (- k * 2 * (if k < 0 then lpi else upi)) ;
+    ux = float_plus_up prec ux (- k * 2 * (if k < 0 then upi else lpi))
   in   if - lpi \<le> lx \<and> ux \<le> 0    then (lb_cos prec (-lx), ub_cos prec (-ux))
   else if 0 \<le> lx \<and> ux \<le> lpi      then (lb_cos prec ux, ub_cos prec lx)
   else if - lpi \<le> lx \<and> ux \<le> lpi  then (min (lb_cos prec (-lx)) (lb_cos prec ux), Float 1 0)
@@ -1120,20 +1311,24 @@
   let ?lpi = "float_round_down prec (lb_pi prec)"
   let ?upi = "float_round_up prec (ub_pi prec)"
   let ?k = "floor_fl (float_divr prec (lx + ?lpi) (2 * ?lpi))"
-  let ?lx = "lx - ?k * 2 * (if ?k < 0 then ?lpi else ?upi)"
-  let ?ux = "ux - ?k * 2 * (if ?k < 0 then ?upi else ?lpi)"
+  let ?lx2 = "(- ?k * 2 * (if ?k < 0 then ?lpi else ?upi))"
+  let ?ux2 = "(- ?k * 2 * (if ?k < 0 then ?upi else ?lpi))"
+  let ?lx = "float_plus_down prec lx ?lx2"
+  let ?ux = "float_plus_up prec ux ?ux2"
 
   obtain k :: int where k: "k = real ?k" using floor_int .
 
   have upi: "pi \<le> ?upi" and lpi: "?lpi \<le> pi"
     using float_round_up[of "ub_pi prec" prec] pi_boundaries[of prec]
           float_round_down[of prec "lb_pi prec"] by auto
-  hence "?lx \<le> x - k * (2 * pi) \<and> x - k * (2 * pi) \<le> ?ux"
-    using x unfolding k[symmetric]
+  hence "lx + ?lx2 \<le> x - k * (2 * pi) \<and> x - k * (2 * pi) \<le> ux + ?ux2"
+    using x
     by (cases "k = 0")
-       (auto intro!: add_mono
-                simp add: k [symmetric] uminus_add_conv_diff [symmetric]
-                simp del: float_of_numeral uminus_add_conv_diff)
+      (auto intro!: add_mono
+        simp add: k [symmetric] uminus_add_conv_diff [symmetric]
+        simp del: float_of_numeral uminus_add_conv_diff)
+  hence "?lx \<le> x - k * (2 * pi) \<and> x - k * (2 * pi) \<le> ?ux"
+    by (auto intro!: float_plus_down_le float_plus_up_le)
   note lx = this[THEN conjunct1] and ux = this[THEN conjunct2]
   hence lx_less_ux: "?lx \<le> real ?ux" by (rule order_trans)
 
@@ -1328,9 +1523,11 @@
 
 fun ub_exp_horner :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float" and lb_exp_horner :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float" where
 "ub_exp_horner prec 0 i k x       = 0" |
-"ub_exp_horner prec (Suc n) i k x = rapprox_rat prec 1 (int k) + x * lb_exp_horner prec n (i + 1) (k * i) x" |
+"ub_exp_horner prec (Suc n) i k x = float_plus_up prec
+    (rapprox_rat prec 1 (int k)) (float_round_up prec (x * lb_exp_horner prec n (i + 1) (k * i) x))" |
 "lb_exp_horner prec 0 i k x       = 0" |
-"lb_exp_horner prec (Suc n) i k x = lapprox_rat prec 1 (int k) + x * ub_exp_horner prec n (i + 1) (k * i) x"
+"lb_exp_horner prec (Suc n) i k x = float_plus_down prec
+    (lapprox_rat prec 1 (int k)) (float_round_down prec (x * ub_exp_horner prec n (i + 1) (k * i) x))"
 
 lemma bnds_exp_horner: assumes "real x \<le> 0"
   shows "exp x \<in> { lb_exp_horner prec (get_even n) 1 1 x .. ub_exp_horner prec (get_odd n) 1 1 x }"
@@ -1377,29 +1574,42 @@
   } ultimately show ?thesis by auto
 qed
 
+lemma ub_exp_horner_nonneg: "real x \<le> 0 \<Longrightarrow> 0 \<le> real (ub_exp_horner prec (get_odd n) (Suc 0) (Suc 0) x)"
+  using bnds_exp_horner[of x prec n]
+  by (intro order_trans[OF exp_ge_zero]) auto
+
+
 subsection "Compute the exponential function on the entire domain"
 
 function ub_exp :: "nat \<Rightarrow> float \<Rightarrow> float" and lb_exp :: "nat \<Rightarrow> float \<Rightarrow> float" where
-"lb_exp prec x = (if 0 < x then float_divl prec 1 (ub_exp prec (-x))
-             else let
-                horner = (\<lambda> x. let  y = lb_exp_horner prec (get_even (prec + 2)) 1 1 x  in if y \<le> 0 then Float 1 (- 2) else y)
-             in if x < - 1 then (horner (float_divl prec x (- floor_fl x))) ^ nat (- int_floor_fl x)
-                           else horner x)" |
-"ub_exp prec x = (if 0 < x    then float_divr prec 1 (lb_exp prec (-x))
-             else if x < - 1  then ub_exp_horner prec (get_odd (prec + 2)) 1 1 (float_divr prec x (- floor_fl x)) ^ (nat (- int_floor_fl x))
-                              else ub_exp_horner prec (get_odd (prec + 2)) 1 1 x)"
+"lb_exp prec x =
+  (if 0 < x then float_divl prec 1 (ub_exp prec (-x))
+  else
+    let
+      horner = (\<lambda> x. let  y = lb_exp_horner prec (get_even (prec + 2)) 1 1 x in
+        if y \<le> 0 then Float 1 (- 2) else y)
+    in
+      if x < - 1 then
+        power_down_fl prec (horner (float_divl prec x (- floor_fl x))) (nat (- int_floor_fl x))
+      else horner x)" |
+"ub_exp prec x =
+  (if 0 < x then float_divr prec 1 (lb_exp prec (-x))
+  else if x < - 1 then
+    power_up_fl prec
+      (ub_exp_horner prec (get_odd (prec + 2)) 1 1
+        (float_divr prec x (- floor_fl x))) (nat (- int_floor_fl x))
+  else ub_exp_horner prec (get_odd (prec + 2)) 1 1 x)"
 by pat_completeness auto
-termination by (relation "measure (\<lambda> v. let (prec, x) = case_sum id id v in (if 0 < x then 1 else 0))", auto)
+termination
+by (relation "measure (\<lambda> v. let (prec, x) = case_sum id id v in (if 0 < x then 1 else 0))", auto)
 
 lemma exp_m1_ge_quarter: "(1 / 4 :: real) \<le> exp (- 1)"
 proof -
   have eq4: "4 = Suc (Suc (Suc (Suc 0)))" by auto
 
   have "1 / 4 = (Float 1 (- 2))" unfolding Float_num by auto
-  also have "\<dots> \<le> lb_exp_horner 1 (get_even 4) 1 1 (- 1)"
-    unfolding get_even_def eq4
-    by (auto simp add: Float.compute_lapprox_rat Float.compute_rapprox_rat
-                  Float.compute_lapprox_posrat Float.compute_rapprox_posrat rat_precision_def Float.compute_bitlen)
+  also have "\<dots> \<le> lb_exp_horner 3 (get_even 3) 1 1 (- 1)"
+    by code_simp
   also have "\<dots> \<le> exp (- 1 :: float)" using bnds_exp_horner[where x="- 1"] by auto
   finally show ?thesis by simp
 qed
@@ -1416,7 +1626,8 @@
   }
   ultimately show ?thesis
     unfolding lb_exp.simps if_not_P[OF `\<not> 0 < x`] Let_def
-    by (cases "floor_fl x", cases "x < - 1", auto)
+    by (cases "floor_fl x", cases "x < - 1", auto simp: real_power_up_fl real_power_down_fl
+      intro!: power_up_less power_down_pos)
 qed
 
 lemma exp_boundaries': assumes "x \<le> 0"
@@ -1458,14 +1669,14 @@
     hence "real (int_floor_fl x) < 0" unfolding less_float_def by auto
     have fl_eq: "real (- int_floor_fl x) = real (- floor_fl x)"
       by (simp add: floor_fl_def int_floor_fl_def)
-    from `0 < - int_floor_fl x` have "0 < real (- floor_fl x)"
+    from `0 < - int_floor_fl x` have "0 \<le> real (- floor_fl x)"
       by (simp add: floor_fl_def int_floor_fl_def)
     from `real (int_floor_fl x) < 0` have "real (floor_fl x) < 0"
       by (simp add: floor_fl_def int_floor_fl_def)
     have "exp x \<le> ub_exp prec x"
     proof -
       have div_less_zero: "real (float_divr prec x (- floor_fl x)) \<le> 0"
-        using float_divr_nonpos_pos_upper_bound[OF `real x \<le> 0` `0 < real (- floor_fl x)`]
+        using float_divr_nonpos_pos_upper_bound[OF `real x \<le> 0` `0 \<le> real (- floor_fl x)`]
         unfolding less_eq_float_def zero_float.rep_eq .
 
       have "exp x = exp (?num * (x / ?num))" using `real ?num \<noteq> 0` by auto
@@ -1475,7 +1686,10 @@
       also have "\<dots> \<le> (?ub_exp_horner (float_divr prec x (- floor_fl x))) ^ ?num"
         unfolding real_of_float_power
         by (rule power_mono, rule bnds_exp_horner[OF div_less_zero, unfolded atLeastAtMost_iff, THEN conjunct2], auto)
-      finally show ?thesis unfolding ub_exp.simps if_not_P[OF `\<not> 0 < x`] if_P[OF `x < - 1`] floor_fl_def Let_def .
+      also have "\<dots> \<le> real (power_up_fl prec (?ub_exp_horner (float_divr prec x (- floor_fl x))) ?num)"
+        by (auto simp add: real_power_up_fl intro!: power_up ub_exp_horner_nonneg div_less_zero)
+      finally show ?thesis unfolding ub_exp.simps if_not_P[OF `\<not> 0 < x`] if_P[OF `x < - 1`] floor_fl_def Let_def
+        .
     qed
     moreover
     have "lb_exp prec x \<le> exp x"
@@ -1497,10 +1711,14 @@
           using float_divl by (auto intro!: power_mono simp del: uminus_float.rep_eq)
         also have "\<dots> = exp (?num * (x / ?num))" unfolding exp_real_of_nat_mult ..
         also have "\<dots> = exp x" using `real ?num \<noteq> 0` by auto
-        finally show ?thesis
-          unfolding lb_exp.simps if_not_P[OF `\<not> 0 < x`] if_P[OF `x < - 1`] int_floor_fl_def Let_def if_not_P[OF False] by auto
+        finally show ?thesis using False
+          unfolding lb_exp.simps if_not_P[OF `\<not> 0 < x`] if_P[OF `x < - 1`] int_floor_fl_def Let_def if_not_P[OF False]
+          by (auto simp: real_power_down_fl intro!: power_down_le)
       next
         case True
+        have "power_down_fl prec (Float 1 (- 2))  ?num \<le> real (Float 1 (- 2)) ^ ?num"
+          by (auto simp: real_power_down_fl power_down)
+        also
         have "real (floor_fl x) \<noteq> 0" and "real (floor_fl x) \<le> 0" using `real (floor_fl x) < 0` by auto
         from divide_right_mono_neg[OF floor_fl[of x] `real (floor_fl x) \<le> 0`, unfolded divide_self[OF `real (floor_fl x) \<noteq> 0`]]
         have "- 1 \<le> x / (- floor_fl x)" unfolding minus_float.rep_eq by auto
@@ -1510,7 +1728,8 @@
           by (auto intro!: power_mono)
         also have "\<dots> = exp x" unfolding num_eq fl_eq exp_real_of_nat_mult[symmetric] using `real (floor_fl x) \<noteq> 0` by auto
         finally show ?thesis
-          unfolding lb_exp.simps if_not_P[OF `\<not> 0 < x`] if_P[OF `x < - 1`] int_floor_fl_def Let_def if_P[OF True] real_of_float_power .
+          unfolding lb_exp.simps if_not_P[OF `\<not> 0 < x`] if_P[OF `x < - 1`] int_floor_fl_def Let_def if_P[OF True] real_of_float_power
+          .
       qed
     qed
     ultimately show ?thesis by auto
@@ -1579,9 +1798,11 @@
 fun ub_ln_horner :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float"
 and lb_ln_horner :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float" where
 "ub_ln_horner prec 0 i x       = 0" |
-"ub_ln_horner prec (Suc n) i x = rapprox_rat prec 1 (int i) - x * lb_ln_horner prec n (Suc i) x" |
+"ub_ln_horner prec (Suc n) i x = float_plus_up prec
+    (rapprox_rat prec 1 (int i)) (- float_round_down prec (x * lb_ln_horner prec n (Suc i) x))" |
 "lb_ln_horner prec 0 i x       = 0" |
-"lb_ln_horner prec (Suc n) i x = lapprox_rat prec 1 (int i) - x * ub_ln_horner prec n (Suc i) x"
+"lb_ln_horner prec (Suc n) i x = float_plus_down prec
+    (lapprox_rat prec 1 (int i)) (- float_round_up prec (x * ub_ln_horner prec n (Suc i) x))"
 
 lemma ln_bounds:
   assumes "0 \<le> x" and "x < 1"
@@ -1646,11 +1867,13 @@
 subsection "Compute the logarithm of 2"
 
 definition ub_ln2 where "ub_ln2 prec = (let third = rapprox_rat (max prec 1) 1 3
-                                        in (Float 1 (- 1) * ub_ln_horner prec (get_odd prec) 1 (Float 1 (- 1))) +
-                                           (third * ub_ln_horner prec (get_odd prec) 1 third))"
+                                        in float_plus_up prec
+                                          ((Float 1 (- 1) * ub_ln_horner prec (get_odd prec) 1 (Float 1 (- 1))))
+                                           (float_round_up prec (third * ub_ln_horner prec (get_odd prec) 1 third)))"
 definition lb_ln2 where "lb_ln2 prec = (let third = lapprox_rat prec 1 3
-                                        in (Float 1 (- 1) * lb_ln_horner prec (get_even prec) 1 (Float 1 (- 1))) +
-                                           (third * lb_ln_horner prec (get_even prec) 1 third))"
+                                        in float_plus_down prec
+                                          ((Float 1 (- 1) * lb_ln_horner prec (get_even prec) 1 (Float 1 (- 1))))
+                                           (float_round_down prec (third * lb_ln_horner prec (get_even prec) 1 third)))"
 
 lemma ub_ln2: "ln 2 \<le> ub_ln2 prec" (is "?ub_ln2")
   and lb_ln2: "lb_ln2 prec \<le> ln 2" (is "?lb_ln2")
@@ -1669,25 +1892,28 @@
   have lb2: "0 \<le> real (Float 1 (- 1))" and ub2: "real (Float 1 (- 1)) < 1" unfolding Float_num by auto
 
   have "0 \<le> (1::int)" and "0 < (3::int)" by auto
-  have ub3_ub: "real ?uthird < 1" by (simp add: Float.compute_rapprox_rat rapprox_posrat_less1)
+  have ub3_ub: "real ?uthird < 1"
+    by (simp add: Float.compute_rapprox_rat Float.compute_lapprox_rat rapprox_posrat_less1)
 
   have third_gt0: "(0 :: real) < 1 / 3 + 1" by auto
   have uthird_gt0: "0 < real ?uthird + 1" using ub3_lb by auto
   have lthird_gt0: "0 < real ?lthird + 1" using lb3_lb by auto
 
-  show ?ub_ln2 unfolding ub_ln2_def Let_def plus_float.rep_eq ln2_sum Float_num(4)[symmetric]
-  proof (rule add_mono, fact ln_float_bounds(2)[OF lb2 ub2])
+  show ?ub_ln2 unfolding ub_ln2_def Let_def ln2_sum Float_num(4)[symmetric]
+  proof (rule float_plus_up_le, rule add_mono, fact ln_float_bounds(2)[OF lb2 ub2])
     have "ln (1 / 3 + 1) \<le> ln (real ?uthird + 1)" unfolding ln_le_cancel_iff[OF third_gt0 uthird_gt0] using ub3 by auto
     also have "\<dots> \<le> ?uthird * ub_ln_horner prec (get_odd prec) 1 ?uthird"
       using ln_float_bounds(2)[OF ub3_lb ub3_ub] .
-    finally show "ln (1 / 3 + 1) \<le> ?uthird * ub_ln_horner prec (get_odd prec) 1 ?uthird" .
+    also note float_round_up
+    finally show "ln (1 / 3 + 1) \<le> float_round_up prec (?uthird * ub_ln_horner prec (get_odd prec) 1 ?uthird)" .
   qed
-  show ?lb_ln2 unfolding lb_ln2_def Let_def plus_float.rep_eq ln2_sum Float_num(4)[symmetric]
-  proof (rule add_mono, fact ln_float_bounds(1)[OF lb2 ub2])
+  show ?lb_ln2 unfolding lb_ln2_def Let_def ln2_sum Float_num(4)[symmetric]
+  proof (rule float_plus_down_le, rule add_mono, fact ln_float_bounds(1)[OF lb2 ub2])
     have "?lthird * lb_ln_horner prec (get_even prec) 1 ?lthird \<le> ln (real ?lthird + 1)"
       using ln_float_bounds(1)[OF lb3_lb lb3_ub] .
+    note float_round_down_le[OF this]
     also have "\<dots> \<le> ln (1 / 3 + 1)" unfolding ln_le_cancel_iff[OF lthird_gt0 third_gt0] using lb3 by auto
-    finally show "?lthird * lb_ln_horner prec (get_even prec) 1 ?lthird \<le> ln (1 / 3 + 1)" .
+    finally show "float_round_down prec (?lthird * lb_ln_horner prec (get_even prec) 1 ?lthird) \<le> ln (1 / 3 + 1)" .
   qed
 qed
 
@@ -1696,25 +1922,25 @@
 function ub_ln :: "nat \<Rightarrow> float \<Rightarrow> float option" and lb_ln :: "nat \<Rightarrow> float \<Rightarrow> float option" where
 "ub_ln prec x = (if x \<le> 0          then None
             else if x < 1          then Some (- the (lb_ln prec (float_divl (max prec 1) 1 x)))
-            else let horner = \<lambda>x. x * ub_ln_horner prec (get_odd prec) 1 x in
+            else let horner = \<lambda>x. float_round_up prec (x * ub_ln_horner prec (get_odd prec) 1 x) in
                  if x \<le> Float 3 (- 1) then Some (horner (x - 1))
-            else if x < Float 1 1  then Some (horner (Float 1 (- 1)) + horner (x * rapprox_rat prec 2 3 - 1))
+            else if x < Float 1 1  then Some (float_round_up prec (horner (Float 1 (- 1)) + horner (x * rapprox_rat prec 2 3 - 1)))
                                    else let l = bitlen (mantissa x) - 1 in
-                                        Some (ub_ln2 prec * (Float (exponent x + l) 0) + horner (Float (mantissa x) (- l) - 1)))" |
+                                        Some (float_plus_up prec (float_round_up prec (ub_ln2 prec * (Float (exponent x + l) 0))) (horner (Float (mantissa x) (- l) - 1))))" |
 "lb_ln prec x = (if x \<le> 0          then None
             else if x < 1          then Some (- the (ub_ln prec (float_divr prec 1 x)))
-            else let horner = \<lambda>x. x * lb_ln_horner prec (get_even prec) 1 x in
+            else let horner = \<lambda>x. float_round_down prec (x * lb_ln_horner prec (get_even prec) 1 x) in
                  if x \<le> Float 3 (- 1) then Some (horner (x - 1))
-            else if x < Float 1 1  then Some (horner (Float 1 (- 1)) +
-                                              horner (max (x * lapprox_rat prec 2 3 - 1) 0))
+            else if x < Float 1 1  then Some (float_round_down prec (horner (Float 1 (- 1)) +
+                                              horner (max (x * lapprox_rat prec 2 3 - 1) 0)))
                                    else let l = bitlen (mantissa x) - 1 in
-                                        Some (lb_ln2 prec * (Float (exponent x + l) 0) + horner (Float (mantissa x) (- l) - 1)))"
+                                        Some (float_plus_down prec (float_round_down prec (lb_ln2 prec * (Float (exponent x + l) 0))) (horner (Float (mantissa x) (- l) - 1))))"
 by pat_completeness auto
 
 termination proof (relation "measure (\<lambda> v. let (prec, x) = case_sum id id v in (if x < 1 then 1 else 0))", auto)
   fix prec and x :: float assume "\<not> real x \<le> 0" and "real x < 1" and "real (float_divl (max prec (Suc 0)) 1 x) < 1"
   hence "0 < real x" "1 \<le> max prec (Suc 0)" "real x < 1" by auto
-  from float_divl_pos_less1_bound[OF `0 < real x` `real x < 1` `1 \<le> max prec (Suc 0)`]
+  from float_divl_pos_less1_bound[OF `0 < real x` `real x < 1`[THEN less_imp_le] `1 \<le> max prec (Suc 0)`]
   show False using `real (float_divl (max prec (Suc 0)) 1 x) < 1` by auto
 next
   fix prec x assume "\<not> real x \<le> 0" and "real x < 1" and "real (float_divr prec 1 x) < 1"
@@ -1763,20 +1989,20 @@
   defines "x \<equiv> Float m e"
   shows "ub_ln prec x = (if x \<le> 0          then None
               else if x < 1          then Some (- the (lb_ln prec (float_divl (max prec 1) 1 x)))
-            else let horner = \<lambda>x. x * ub_ln_horner prec (get_odd prec) 1 x in
+            else let horner = \<lambda>x. float_round_up prec (x * ub_ln_horner prec (get_odd prec) 1 x) in
                  if x \<le> Float 3 (- 1) then Some (horner (x - 1))
-            else if x < Float 1 1  then Some (horner (Float 1 (- 1)) + horner (x * rapprox_rat prec 2 3 - 1))
+            else if x < Float 1 1  then Some (float_round_up prec (horner (Float 1 (- 1)) + horner (x * rapprox_rat prec 2 3 - 1)))
                                    else let l = bitlen m - 1 in
-                                        Some (ub_ln2 prec * (Float (e + l) 0) + horner (Float m (- l) - 1)))"
+                                        Some (float_plus_up prec (float_round_up prec (ub_ln2 prec * (Float (e + l) 0))) (horner (Float m (- l) - 1))))"
     (is ?th1)
   and "lb_ln prec x = (if x \<le> 0          then None
             else if x < 1          then Some (- the (ub_ln prec (float_divr prec 1 x)))
-            else let horner = \<lambda>x. x * lb_ln_horner prec (get_even prec) 1 x in
+            else let horner = \<lambda>x. float_round_down prec (x * lb_ln_horner prec (get_even prec) 1 x) in
                  if x \<le> Float 3 (- 1) then Some (horner (x - 1))
-            else if x < Float 1 1  then Some (horner (Float 1 (- 1)) +
-                                              horner (max (x * lapprox_rat prec 2 3 - 1) 0))
+            else if x < Float 1 1  then Some (float_round_down prec (horner (Float 1 (- 1)) +
+                                              horner (max (x * lapprox_rat prec 2 3 - 1) 0)))
                                    else let l = bitlen m - 1 in
-                                        Some (lb_ln2 prec * (Float (e + l) 0) + horner (Float m (- l) - 1)))"
+                                        Some (float_plus_down prec (float_round_down prec (lb_ln2 prec * (Float (e + l) 0))) (horner (Float m (- l) - 1))))"
     (is ?th2)
 proof -
   from assms Float_pos_eq_mantissa_pos have "x > 0 \<Longrightarrow> m > 0" by simp
@@ -1826,7 +2052,7 @@
     case True
     show ?thesis unfolding lb_ln.simps unfolding ub_ln.simps Let_def
       using ln_float_bounds[OF `0 \<le> real (x - 1)` `real (x - 1) < 1`, of prec] `\<not> x \<le> 0` `\<not> x < 1` True
-      by auto
+      by (auto intro!: float_round_down_le float_round_up_le)
   next
     case False hence *: "3 / 2 < x" by auto
 
@@ -1834,8 +2060,8 @@
     have add: "ln x = ln (3 / 2) + ln (real x * 2 / 3)"
       by (auto simp add: algebra_simps diff_divide_distrib)
 
-    let "?ub_horner x" = "x * ub_ln_horner prec (get_odd prec) 1 x"
-    let "?lb_horner x" = "x * lb_ln_horner prec (get_even prec) 1 x"
+    let "?ub_horner x" = "float_round_up prec (x * ub_ln_horner prec (get_odd prec) 1 x)"
+    let "?lb_horner x" = "float_round_down prec (x * lb_ln_horner prec (get_even prec) 1 x)"
 
     { have up: "real (rapprox_rat prec 2 3) \<le> 1"
         by (rule rapprox_rat_le1) simp_all
@@ -1853,15 +2079,17 @@
         show "0 < real (x * rapprox_rat prec 2 3 - 1) + 1" using pos by auto
       qed
       also have "\<dots> \<le> ?ub_horner (x * rapprox_rat prec 2 3 - 1)"
-      proof (rule ln_float_bounds(2))
+      proof (rule float_round_up_le, rule ln_float_bounds(2))
         from mult_less_le_imp_less[OF `real x < 2` up] low *
         show "real (x * rapprox_rat prec 2 3 - 1) < 1" by auto
         show "0 \<le> real (x * rapprox_rat prec 2 3 - 1)" using pos by auto
       qed
-      finally have "ln x
-        \<le> ?ub_horner (Float 1 (- 1))
-          + ?ub_horner (x * rapprox_rat prec 2 3 - 1)"
-        using ln_float_bounds(2)[of "Float 1 (- 1)" prec prec] add by auto }
+     finally have "ln x \<le> ?ub_horner (Float 1 (-1))
+          + ?ub_horner ((x * rapprox_rat prec 2 3 - 1))"
+        using ln_float_bounds(2)[of "Float 1 (- 1)" prec prec] add
+        by (auto intro!: add_mono float_round_up_le)
+      note float_round_up_le[OF this, of prec]
+    }
     moreover
     { let ?max = "max (x * lapprox_rat prec 2 3 - 1) 0"
 
@@ -1873,7 +2101,7 @@
 
       have "?lb_horner ?max
         \<le> ln (real ?max + 1)"
-      proof (rule ln_float_bounds(1))
+      proof (rule float_round_down_le, rule ln_float_bounds(1))
         from mult_less_le_imp_less[OF `real x < 2` up] * low
         show "real ?max < 1" by (cases "real (lapprox_rat prec 2 3) = 0",
           auto simp add: real_of_float_max)
@@ -1887,9 +2115,11 @@
           by (cases "0 < real x * real (lapprox_posrat prec 2 3) - 1",
               auto simp add: max_def)
       qed
-      finally have "?lb_horner (Float 1 (- 1)) + ?lb_horner ?max
-        \<le> ln x"
-        using ln_float_bounds(1)[of "Float 1 (- 1)" prec prec] add by auto }
+      finally have "?lb_horner (Float 1 (- 1)) + ?lb_horner ?max \<le> ln x"
+        using ln_float_bounds(1)[of "Float 1 (- 1)" prec prec] add
+        by (auto intro!: add_mono float_round_down_le)
+      note float_round_down_le[OF this, of prec]
+    }
     ultimately
     show ?thesis unfolding lb_ln.simps unfolding ub_ln.simps Let_def
       using `\<not> x \<le> 0` `\<not> x < 1` True False by auto
@@ -1917,7 +2147,8 @@
          (auto simp : powr_minus field_simps inverse_eq_divide)
 
     {
-      have "lb_ln2 prec * ?s \<le> ln 2 * (e + (bitlen m - 1))" (is "?lb2 \<le> _")
+      have "float_round_down prec (lb_ln2 prec * ?s) \<le> ln 2 * (e + (bitlen m - 1))" (is "real ?lb2 \<le> _")
+        apply (rule float_round_down_le)
         unfolding nat_0 power_0 mult_1_right times_float.rep_eq
         using lb_ln2[of prec]
       proof (rule mult_mono)
@@ -1926,16 +2157,19 @@
       qed auto
       moreover
       from ln_float_bounds(1)[OF x_bnds]
-      have "(?x - 1) * lb_ln_horner prec (get_even prec) 1 (?x - 1) \<le> ln ?x" (is "?lb_horner \<le> _") by auto
-      ultimately have "?lb2 + ?lb_horner \<le> ln x"
-        unfolding Float ln_shifted_float[OF `0 < m`, of e] by auto
+      have "float_round_down prec ((?x - 1) * lb_ln_horner prec (get_even prec) 1 (?x - 1)) \<le> ln ?x" (is "real ?lb_horner \<le> _")
+        by (auto intro!: float_round_down_le)
+      ultimately have "float_plus_down prec ?lb2 ?lb_horner \<le> ln x"
+        unfolding Float ln_shifted_float[OF `0 < m`, of e] by (auto intro!: float_plus_down_le)
     }
     moreover
     {
       from ln_float_bounds(2)[OF x_bnds]
-      have "ln ?x \<le> (?x - 1) * ub_ln_horner prec (get_odd prec) 1 (?x - 1)" (is "_ \<le> ?ub_horner") by auto
+      have "ln ?x \<le> float_round_up prec ((?x - 1) * ub_ln_horner prec (get_odd prec) 1 (?x - 1))" (is "_ \<le> real ?ub_horner")
+        by (auto intro!: float_round_up_le)
       moreover
-      have "ln 2 * (e + (bitlen m - 1)) \<le> ub_ln2 prec * ?s" (is "_ \<le> ?ub2")
+      have "ln 2 * (e + (bitlen m - 1)) \<le> float_round_up prec (ub_ln2 prec * ?s)" (is "_ \<le> real ?ub2")
+        apply (rule float_round_up_le)
         unfolding nat_0 power_0 mult_1_right times_float.rep_eq
         using ub_ln2[of prec]
       proof (rule mult_mono)
@@ -1945,8 +2179,9 @@
         have "0 \<le> ln 2" by simp
         thus "0 \<le> real (ub_ln2 prec)" using ub_ln2[of prec] by arith
       qed auto
-      ultimately have "ln x \<le> ?ub2 + ?ub_horner"
-        unfolding Float ln_shifted_float[OF `0 < m`, of e] by auto
+      ultimately have "ln x \<le> float_plus_up prec ?ub2 ?ub_horner"
+        unfolding Float ln_shifted_float[OF `0 < m`, of e]
+        by (auto intro!: float_plus_up_le)
     }
     ultimately show ?thesis unfolding lb_ln.simps unfolding ub_ln.simps
       unfolding if_not_P[OF `\<not> x \<le> 0`] if_not_P[OF `\<not> x < 1`] if_not_P[OF False] if_not_P[OF `\<not> x \<le> Float 3 (- 1)`] Let_def
@@ -1963,13 +2198,13 @@
   show ?thesis using ub_ln_lb_ln_bounds'[OF `1 \<le> x`] .
 next
   case True have "\<not> x \<le> 0" using `0 < x` by auto
-  from True have "real x < 1" by simp
+  from True have "real x \<le> 1" "x \<le> 1" by simp_all
   have "0 < real x" and "real x \<noteq> 0" using `0 < x` by auto
   hence A: "0 < 1 / real x" by auto
 
   {
     let ?divl = "float_divl (max prec 1) 1 x"
-    have A': "1 \<le> ?divl" using float_divl_pos_less1_bound[OF `0 < real x` `real x < 1`] by auto
+    have A': "1 \<le> ?divl" using float_divl_pos_less1_bound[OF `0 < real x` `real x \<le> 1`] by auto
     hence B: "0 < real ?divl" by auto
 
     have "ln ?divl \<le> ln (1 / x)" unfolding ln_le_cancel_iff[OF B A] using float_divl[of _ 1 x] by auto
@@ -1979,7 +2214,7 @@
   } moreover
   {
     let ?divr = "float_divr prec 1 x"
-    have A': "1 \<le> ?divr" using float_divr_pos_less1_lower_bound[OF `0 < x` `x < 1`] unfolding less_eq_float_def less_float_def by auto
+    have A': "1 \<le> ?divr" using float_divr_pos_less1_lower_bound[OF `0 < x` `x \<le> 1`] unfolding less_eq_float_def less_float_def by auto
     hence B: "0 < real ?divr" by auto
 
     have "ln (1 / x) \<le> ln ?divr" unfolding ln_le_cancel_iff[OF A B] using float_divr[of 1 x] by auto
@@ -2162,11 +2397,19 @@
 
 fun approx approx' :: "nat \<Rightarrow> floatarith \<Rightarrow> (float * float) option list \<Rightarrow> (float * float) option" where
 "approx' prec a bs          = (case (approx prec a bs) of Some (l, u) \<Rightarrow> Some (float_round_down prec l, float_round_up prec u) | None \<Rightarrow> None)" |
-"approx prec (Add a b) bs   = lift_bin' (approx' prec a bs) (approx' prec b bs) (\<lambda> l1 u1 l2 u2. (l1 + l2, u1 + u2))" |
+"approx prec (Add a b) bs   =
+  lift_bin' (approx' prec a bs) (approx' prec b bs)
+    (\<lambda> l1 u1 l2 u2. (float_plus_down prec l1 l2, float_plus_up prec u1 u2))" |
 "approx prec (Minus a) bs   = lift_un' (approx' prec a bs) (\<lambda> l u. (-u, -l))" |
-"approx prec (Mult a b) bs  = lift_bin' (approx' prec a bs) (approx' prec b bs)
-                                    (\<lambda> a1 a2 b1 b2. (nprt a1 * pprt b2 + nprt a2 * nprt b2 + pprt a1 * pprt b1 + pprt a2 * nprt b1,
-                                                     pprt a2 * pprt b2 + pprt a1 * nprt b2 + nprt a2 * pprt b1 + nprt a1 * nprt b1))" |
+"approx prec (Mult a b) bs  =
+  lift_bin' (approx' prec a bs) (approx' prec b bs)
+    (\<lambda> a1 a2 b1 b2.
+      (float_plus_down prec (nprt a1 * pprt b2)
+          (float_plus_down prec (nprt a2 * nprt b2)
+            (float_plus_down prec (pprt a1 * pprt b1) (pprt a2 * nprt b1))),
+        float_plus_up prec (pprt a2 * pprt b2)
+            (float_plus_up prec (pprt a1 * nprt b2)
+              (float_plus_up prec (nprt a2 * pprt b1) (nprt a1 * nprt b1)))))" |
 "approx prec (Inverse a) bs = lift_un (approx' prec a bs) (\<lambda> l u. if (0 < l \<or> u < 0) then (Some (float_divl prec 1 u), Some (float_divr prec 1 l)) else (None, None))" |
 "approx prec (Cos a) bs     = lift_un' (approx' prec a bs) (bnds_cos prec)" |
 "approx prec Pi bs          = Some (lb_pi prec, ub_pi prec)" |
@@ -2177,7 +2420,7 @@
 "approx prec (Sqrt a) bs    = lift_un' (approx' prec a bs) (\<lambda> l u. (lb_sqrt prec l, ub_sqrt prec u))" |
 "approx prec (Exp a) bs     = lift_un' (approx' prec a bs) (\<lambda> l u. (lb_exp prec l, ub_exp prec u))" |
 "approx prec (Ln a) bs      = lift_un (approx' prec a bs) (\<lambda> l u. (lb_ln prec l, ub_ln prec u))" |
-"approx prec (Power a n) bs = lift_un' (approx' prec a bs) (float_power_bnds n)" |
+"approx prec (Power a n) bs = lift_un' (approx' prec a bs) (float_power_bnds prec n)" |
 "approx prec (Num f) bs     = Some (f, f)" |
 "approx prec (Var i) bs    = (if i < length bs then bs ! i else None)"
 
@@ -2367,10 +2610,10 @@
 proof (induct arith arbitrary: l u)
   case (Add a b)
   from lift_bin'[OF Add.prems[unfolded approx.simps]] Add.hyps
-  obtain l1 u1 l2 u2 where "l = l1 + l2" and "u = u1 + u2"
+  obtain l1 u1 l2 u2 where "l = float_plus_down prec l1 l2" and "u = float_plus_up prec u1 u2"
     "l1 \<le> interpret_floatarith a xs" and "interpret_floatarith a xs \<le> u1"
     "l2 \<le> interpret_floatarith b xs" and "interpret_floatarith b xs \<le> u2" unfolding fst_conv snd_conv by blast
-  thus ?case unfolding interpret_floatarith.simps by auto
+  thus ?case unfolding interpret_floatarith.simps by (auto intro!: float_plus_up_le float_plus_down_le)
 next
   case (Minus a)
   from lift_un'[OF Minus.prems[unfolded approx.simps]] Minus.hyps
@@ -2381,12 +2624,18 @@
   case (Mult a b)
   from lift_bin'[OF Mult.prems[unfolded approx.simps]] Mult.hyps
   obtain l1 u1 l2 u2
-    where l: "l = nprt l1 * pprt u2 + nprt u1 * nprt u2 + pprt l1 * pprt l2 + pprt u1 * nprt l2"
-    and u: "u = pprt u1 * pprt u2 + pprt l1 * nprt u2 + nprt u1 * pprt l2 + nprt l1 * nprt l2"
+    where l: "l = float_plus_down prec (nprt l1 * pprt u2) (float_plus_down prec (nprt u1 * nprt u2) (float_plus_down prec (pprt l1 * pprt l2) (pprt u1 * nprt l2)))"
+    and u: "u = float_plus_up prec (pprt u1 * pprt u2) (float_plus_up prec (pprt l1 * nprt u2) (float_plus_up prec (nprt u1 * pprt l2) (nprt l1 * nprt l2)))"
     and "l1 \<le> interpret_floatarith a xs" and "interpret_floatarith a xs \<le> u1"
     and "l2 \<le> interpret_floatarith b xs" and "interpret_floatarith b xs \<le> u2" unfolding fst_conv snd_conv by blast
-  thus ?case unfolding interpret_floatarith.simps l u
+  hence bnds:
+    "nprt l1 * pprt u2 + nprt u1 * nprt u2 + pprt l1 * pprt l2 + pprt u1 * nprt l2 \<le> interpret_floatarith (Mult a b) xs" (is "?l \<le> _")
+    "interpret_floatarith (Mult a b) xs \<le> pprt u1 * pprt u2 + pprt l1 * nprt u2 + nprt u1 * pprt l2 + nprt l1 * nprt l2" (is "_ \<le> ?u")
+    unfolding interpret_floatarith.simps l u
     using mult_le_prts mult_ge_prts by auto
+  from l u have "l \<le> ?l" "?u \<le> u"
+    by (auto intro!: float_plus_up_le float_plus_down_le)
+  thus ?case using bnds by simp
 next
   case (Inverse a)
   from lift_un[OF Inverse.prems[unfolded approx.simps], unfolded if_distrib[of fst] if_distrib[of snd] fst_conv snd_conv] Inverse.hyps
@@ -2464,13 +2713,17 @@
               | Less floatarith floatarith
               | LessEqual floatarith floatarith
               | AtLeastAtMost floatarith floatarith floatarith
+              | Conj form form
+              | Disj form form
 
 fun interpret_form :: "form \<Rightarrow> real list \<Rightarrow> bool" where
 "interpret_form (Bound x a b f) vs = (interpret_floatarith x vs \<in> { interpret_floatarith a vs .. interpret_floatarith b vs } \<longrightarrow> interpret_form f vs)" |
 "interpret_form (Assign x a f) vs  = (interpret_floatarith x vs = interpret_floatarith a vs \<longrightarrow> interpret_form f vs)" |
 "interpret_form (Less a b) vs      = (interpret_floatarith a vs < interpret_floatarith b vs)" |
 "interpret_form (LessEqual a b) vs = (interpret_floatarith a vs \<le> interpret_floatarith b vs)" |
-"interpret_form (AtLeastAtMost x a b) vs = (interpret_floatarith x vs \<in> { interpret_floatarith a vs .. interpret_floatarith b vs })"
+"interpret_form (AtLeastAtMost x a b) vs = (interpret_floatarith x vs \<in> { interpret_floatarith a vs .. interpret_floatarith b vs })" |
+"interpret_form (Conj f g) vs \<longleftrightarrow> interpret_form f vs \<and> interpret_form g vs" |
+"interpret_form (Disj f g) vs \<longleftrightarrow> interpret_form f vs \<or> interpret_form g vs"
 
 fun approx_form' and approx_form :: "nat \<Rightarrow> form \<Rightarrow> (float * float) option list \<Rightarrow> nat list \<Rightarrow> bool" where
 "approx_form' prec f 0 n l u bs ss = approx_form prec f (bs[n := Some (l, u)]) ss" |
@@ -2487,16 +2740,18 @@
     | _ \<Rightarrow> False)" |
 "approx_form prec (Less a b) bs ss =
    (case (approx prec a bs, approx prec b bs)
-   of (Some (l, u), Some (l', u')) \<Rightarrow> u < l'
+   of (Some (l, u), Some (l', u')) \<Rightarrow> float_plus_up prec u (-l') < 0
     | _ \<Rightarrow> False)" |
 "approx_form prec (LessEqual a b) bs ss =
    (case (approx prec a bs, approx prec b bs)
-   of (Some (l, u), Some (l', u')) \<Rightarrow> u \<le> l'
+   of (Some (l, u), Some (l', u')) \<Rightarrow> float_plus_up prec u (-l') \<le> 0
     | _ \<Rightarrow> False)" |
 "approx_form prec (AtLeastAtMost x a b) bs ss =
    (case (approx prec x bs, approx prec a bs, approx prec b bs)
-   of (Some (lx, ux), Some (l, u), Some (l', u')) \<Rightarrow> u \<le> lx \<and> ux \<le> l'
+   of (Some (lx, ux), Some (l, u), Some (l', u')) \<Rightarrow> float_plus_up prec u (-lx) \<le> 0 \<and> float_plus_up prec ux (-l') \<le> 0
     | _ \<Rightarrow> False)" |
+"approx_form prec (Conj a b) bs ss \<longleftrightarrow> approx_form prec a bs ss \<and> approx_form prec b bs ss" |
+"approx_form prec (Disj a b) bs ss \<longleftrightarrow> approx_form prec a bs ss \<or> approx_form prec b bs ss" |
 "approx_form _ _ _ _ = False"
 
 lemma lazy_conj: "(if A then B else False) = (A \<and> B)" by simp
@@ -2586,20 +2841,22 @@
   then obtain l u l' u'
     where l_eq: "Some (l, u) = approx prec a vs"
       and u_eq: "Some (l', u') = approx prec b vs"
-      and inequality: "u < l'"
+      and inequality: "real (float_plus_up prec u (-l')) < 0"
     by (cases "approx prec a vs", auto,
       cases "approx prec b vs", auto)
-  from inequality approx[OF Less.prems(2) l_eq] approx[OF Less.prems(2) u_eq]
+  from le_less_trans[OF float_plus_up inequality]
+    approx[OF Less.prems(2) l_eq] approx[OF Less.prems(2) u_eq]
   show ?case by auto
 next
   case (LessEqual a b)
   then obtain l u l' u'
     where l_eq: "Some (l, u) = approx prec a vs"
       and u_eq: "Some (l', u') = approx prec b vs"
-      and inequality: "u \<le> l'"
+      and inequality: "real (float_plus_up prec u (-l')) \<le> 0"
     by (cases "approx prec a vs", auto,
       cases "approx prec b vs", auto)
-  from inequality approx[OF LessEqual.prems(2) l_eq] approx[OF LessEqual.prems(2) u_eq]
+  from order_trans[OF float_plus_up inequality]
+    approx[OF LessEqual.prems(2) l_eq] approx[OF LessEqual.prems(2) u_eq]
   show ?case by auto
 next
   case (AtLeastAtMost x a b)
@@ -2607,13 +2864,14 @@
     where x_eq: "Some (lx, ux) = approx prec x vs"
     and l_eq: "Some (l, u) = approx prec a vs"
     and u_eq: "Some (l', u') = approx prec b vs"
-    and inequality: "u \<le> lx \<and> ux \<le> l'"
+    and inequality: "real (float_plus_up prec u (-lx)) \<le> 0" "real (float_plus_up prec ux (-l')) \<le> 0"
     by (cases "approx prec x vs", auto,
       cases "approx prec a vs", auto,
       cases "approx prec b vs", auto)
-  from inequality approx[OF AtLeastAtMost.prems(2) l_eq] approx[OF AtLeastAtMost.prems(2) u_eq] approx[OF AtLeastAtMost.prems(2) x_eq]
+  from order_trans[OF float_plus_up inequality(1)] order_trans[OF float_plus_up inequality(2)]
+    approx[OF AtLeastAtMost.prems(2) l_eq] approx[OF AtLeastAtMost.prems(2) u_eq] approx[OF AtLeastAtMost.prems(2) x_eq]
   show ?case by auto
-qed
+qed auto
 
 lemma approx_form:
   assumes "n = length xs"
@@ -3136,19 +3394,26 @@
   show ?thesis by auto
 qed
 
+fun approx_tse_concl where
+"approx_tse_concl prec t (Less lf rt) s l u l' u' \<longleftrightarrow>
+    approx_tse_form' prec t (Add rt (Minus lf)) s l u' (\<lambda> l u. 0 < l)" |
+"approx_tse_concl prec t (LessEqual lf rt) s l u l' u' \<longleftrightarrow>
+    approx_tse_form' prec t (Add rt (Minus lf)) s l u' (\<lambda> l u. 0 \<le> l)" |
+"approx_tse_concl prec t (AtLeastAtMost x lf rt) s l u l' u' \<longleftrightarrow>
+    (if approx_tse_form' prec t (Add x (Minus lf)) s l u' (\<lambda> l u. 0 \<le> l) then
+      approx_tse_form' prec t (Add rt (Minus x)) s l u' (\<lambda> l u. 0 \<le> l) else False)" |
+"approx_tse_concl prec t (Conj f g) s l u l' u' \<longleftrightarrow>
+    approx_tse_concl prec t f s l u l' u' \<and> approx_tse_concl prec t g s l u l' u'" |
+"approx_tse_concl prec t (Disj f g) s l u l' u' \<longleftrightarrow>
+    approx_tse_concl prec t f s l u l' u' \<or> approx_tse_concl prec t g s l u l' u'" |
+"approx_tse_concl _ _ _ _ _ _ _ _ \<longleftrightarrow> False"
+
 definition
 "approx_tse_form prec t s f =
   (case f
    of (Bound x a b f) \<Rightarrow> x = Var 0 \<and>
      (case (approx prec a [None], approx prec b [None])
-      of (Some (l, u), Some (l', u')) \<Rightarrow>
-        (case f
-         of Less lf rt \<Rightarrow> approx_tse_form' prec t (Add rt (Minus lf)) s l u' (\<lambda> l u. 0 < l)
-          | LessEqual lf rt \<Rightarrow> approx_tse_form' prec t (Add rt (Minus lf)) s l u' (\<lambda> l u. 0 \<le> l)
-          | AtLeastAtMost x lf rt \<Rightarrow>
-            (if approx_tse_form' prec t (Add x (Minus lf)) s l u' (\<lambda> l u. 0 \<le> l) then
-            approx_tse_form' prec t (Add rt (Minus x)) s l u' (\<lambda> l u. 0 \<le> l) else False)
-          | _ \<Rightarrow> False)
+      of (Some (l, u), Some (l', u')) \<Rightarrow> approx_tse_concl prec t f s l u l' u'
        | _ \<Rightarrow> False)
    | _ \<Rightarrow> False)"
 
@@ -3171,48 +3436,32 @@
     have bnd: "x \<in> { l .. u'}" unfolding bounded_by_def i by auto
 
     have "interpret_form f' [x]"
-    proof (cases f')
+      using assms[unfolded Bound]
+    proof (induct f')
       case (Less lf rt)
-      with Bound a b assms
+      with a b
       have "approx_tse_form' prec t (Add rt (Minus lf)) s l u' (\<lambda> l u. 0 < l)"
         unfolding approx_tse_form_def by auto
       from approx_tse_form'_less[OF this bnd]
-      show ?thesis using Less by auto
+      show ?case using Less by auto
     next
       case (LessEqual lf rt)
       with Bound a b assms
       have "approx_tse_form' prec t (Add rt (Minus lf)) s l u' (\<lambda> l u. 0 \<le> l)"
         unfolding approx_tse_form_def by auto
       from approx_tse_form'_le[OF this bnd]
-      show ?thesis using LessEqual by auto
+      show ?case using LessEqual by auto
     next
       case (AtLeastAtMost x lf rt)
       with Bound a b assms
       have "approx_tse_form' prec t (Add rt (Minus x)) s l u' (\<lambda> l u. 0 \<le> l)"
         and "approx_tse_form' prec t (Add x (Minus lf)) s l u' (\<lambda> l u. 0 \<le> l)"
-        unfolding approx_tse_form_def lazy_conj by auto
+        unfolding approx_tse_form_def lazy_conj by (auto split: split_if_asm)
       from approx_tse_form'_le[OF this(1) bnd] approx_tse_form'_le[OF this(2) bnd]
-      show ?thesis using AtLeastAtMost by auto
-    next
-      case (Bound x a b f') with assms
-      show ?thesis by (auto elim!: case_optionE simp add: f_def approx_tse_form_def)
-    next
-      case (Assign x a f') with assms
-      show ?thesis by (auto elim!: case_optionE simp add: f_def approx_tse_form_def)
-    qed } thus ?thesis unfolding f_def by auto
-next
-  case Assign
-  with assms show ?thesis by (auto simp add: approx_tse_form_def)
-next
-  case LessEqual
-  with assms show ?thesis by (auto simp add: approx_tse_form_def)
-next
-  case Less
-  with assms show ?thesis by (auto simp add: approx_tse_form_def)
-next
-  case AtLeastAtMost
-  with assms show ?thesis by (auto simp add: approx_tse_form_def)
-qed
+      show ?case using AtLeastAtMost by auto
+    qed (auto simp: f_def approx_tse_form_def elim!: case_optionE)
+  } thus ?thesis unfolding f_def by auto
+qed (insert assms, auto simp add: approx_tse_form_def)
 
 text {* @{term approx_form_eval} is only used for the {\tt value}-command. *}
 
@@ -3245,7 +3494,8 @@
     | term_of_bool false = @{term False};
 
   val mk_int = HOLogic.mk_number @{typ int} o @{code integer_of_int};
-  val dest_int = @{code int_of_integer} o snd o HOLogic.dest_number;
+  fun dest_int (@{term int_of_integer} $ j) = @{code int_of_integer} (snd (HOLogic.dest_number j))
+    | dest_int i = @{code int_of_integer} (snd (HOLogic.dest_number i));
 
   fun term_of_float (@{code Float} (k, l)) =
     @{term Float} $ mk_int k $ mk_int l;
@@ -3291,6 +3541,10 @@
         (floatarith_of_term a, floatarith_of_term b)
     | form_of_term (@{term LessEqual} $ a $ b) = @{code LessEqual}
         (floatarith_of_term a, floatarith_of_term b)
+    | form_of_term (@{term Conj} $ a $ b) = @{code Conj}
+        (form_of_term a, form_of_term b)
+    | form_of_term (@{term Disj} $ a $ b) = @{code Disj}
+        (form_of_term a, form_of_term b)
     | form_of_term (@{term AtLeastAtMost} $ a $ b $ c) = @{code AtLeastAtMost}
         (floatarith_of_term a, floatarith_of_term b, floatarith_of_term c)
     | form_of_term t = bad t;
@@ -3453,4 +3707,11 @@
 
 ML_file "approximation.ML"
 
+
+section "Quickcheck Generator"
+
+ML_file "approximation_generator.ML"
+
+setup "Approximation_Generator.setup"
+
 end
diff -r 302104d8366b -r 87d4ce309e04 src/HOL/Decision_Procs/Decision_Procs.thy
--- a/src/HOL/Decision_Procs/Decision_Procs.thy	Wed Nov 12 18:18:38 2014 +0100
+++ b/src/HOL/Decision_Procs/Decision_Procs.thy	Wed Nov 12 19:30:56 2014 +0100
@@ -10,6 +10,7 @@
   Commutative_Ring_Complete
   "ex/Commutative_Ring_Ex"
   "ex/Approximation_Ex"
+  "ex/Approximation_Quickcheck_Ex"
   "ex/Dense_Linear_Order_Ex"
 begin
 
diff -r 302104d8366b -r 87d4ce309e04 src/HOL/Decision_Procs/approximation_generator.ML
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/HOL/Decision_Procs/approximation_generator.ML	Wed Nov 12 19:30:56 2014 +0100
@@ -0,0 +1,211 @@
+(*  Title:      HOL/Decision_Procs/approximation_generator.ML
+    Author:     Fabian Immler, TU Muenchen
+*)
+
+signature APPROXIMATION_GENERATOR =
+sig
+  val custom_seed: int Config.T
+  val precision: int Config.T
+  val epsilon: real Config.T
+  val approximation_generator:
+    Proof.context ->
+    (term * term list) list ->
+    bool -> int list -> (bool * term list) option * Quickcheck.report option
+  val setup: theory -> theory
+end;
+
+structure Approximation_Generator : APPROXIMATION_GENERATOR =
+struct
+
+val custom_seed = Attrib.setup_config_int @{binding quickcheck_approximation_custom_seed} (K ~1)
+
+val precision = Attrib.setup_config_int @{binding quickcheck_approximation_precision} (K 30)
+
+val epsilon = Attrib.setup_config_real @{binding quickcheck_approximation_epsilon} (K 0.0)
+
+val random_float = @{code "random_class.random::_ \<Rightarrow> _ \<Rightarrow> (float \<times> (unit \<Rightarrow> term)) \<times> _"}
+
+fun nat_of_term t =
+  (HOLogic.dest_nat t handle TERM _ => snd (HOLogic.dest_number t)
+    handle TERM _ => raise TERM ("nat_of_term", [t]));
+
+fun int_of_term t = snd (HOLogic.dest_number t) handle TERM _ => raise TERM ("int_of_term", [t]);
+
+fun real_of_man_exp m e = Real.fromManExp {man = Real.fromInt m, exp = e}
+
+fun mapprox_float (@{term Float} $ m $ e) = real_of_man_exp (int_of_term m) (int_of_term e)
+  | mapprox_float t = Real.fromInt (snd (HOLogic.dest_number t))
+      handle TERM _ => raise TERM ("mapprox_float", [t]);
+
+(* TODO: define using compiled terms? *)
+fun mapprox_floatarith (@{term Add} $ a $ b) xs = mapprox_floatarith a xs + mapprox_floatarith b xs
+  | mapprox_floatarith (@{term Minus} $ a) xs = ~ (mapprox_floatarith a xs)
+  | mapprox_floatarith (@{term Mult} $ a $ b) xs = mapprox_floatarith a xs * mapprox_floatarith b xs
+  | mapprox_floatarith (@{term Inverse} $ a) xs = 1.0 / mapprox_floatarith a xs
+  | mapprox_floatarith (@{term Cos} $ a) xs = Math.cos (mapprox_floatarith a xs)
+  | mapprox_floatarith (@{term Arctan} $ a) xs = Math.atan (mapprox_floatarith a xs)
+  | mapprox_floatarith (@{term Abs} $ a) xs = abs (mapprox_floatarith a xs)
+  | mapprox_floatarith (@{term Max} $ a $ b) xs =
+      Real.max (mapprox_floatarith a xs, mapprox_floatarith b xs)
+  | mapprox_floatarith (@{term Min} $ a $ b) xs =
+      Real.min (mapprox_floatarith a xs, mapprox_floatarith b xs)
+  | mapprox_floatarith @{term Pi} _ = Math.pi
+  | mapprox_floatarith (@{term Sqrt} $ a) xs = Math.sqrt (mapprox_floatarith a xs)
+  | mapprox_floatarith (@{term Exp} $ a) xs = Math.exp (mapprox_floatarith a xs)
+  | mapprox_floatarith (@{term Ln} $ a) xs = Math.ln (mapprox_floatarith a xs)
+  | mapprox_floatarith (@{term Power} $ a $ n) xs =
+      Math.pow (mapprox_floatarith a xs, Real.fromInt (nat_of_term n))
+  | mapprox_floatarith (@{term Var} $ n) xs = nth xs (nat_of_term n)
+  | mapprox_floatarith (@{term Num} $ m) _ = mapprox_float m
+  | mapprox_floatarith t _ = raise TERM ("mapprox_floatarith", [t])
+
+fun mapprox_atLeastAtMost eps x a b xs =
+    let
+      val x' = mapprox_floatarith x xs
+    in
+      mapprox_floatarith a xs + eps <= x' andalso x' + eps <= mapprox_floatarith b xs
+    end
+
+fun mapprox_form eps (@{term Bound} $ x $ a $ b $ f) xs =
+    (not (mapprox_atLeastAtMost eps x a b xs)) orelse mapprox_form eps f xs
+| mapprox_form eps (@{term Assign} $ x $ a $ f) xs =
+    (Real.!= (mapprox_floatarith x xs, mapprox_floatarith a xs)) orelse mapprox_form eps f xs
+| mapprox_form eps (@{term Less} $ a $ b) xs = mapprox_floatarith a xs + eps < mapprox_floatarith b xs
+| mapprox_form eps (@{term LessEqual} $ a $ b) xs = mapprox_floatarith a xs + eps <= mapprox_floatarith b xs
+| mapprox_form eps (@{term AtLeastAtMost} $ x $ a $ b) xs = mapprox_atLeastAtMost eps x a b xs
+| mapprox_form eps (@{term Conj} $ f $ g) xs = mapprox_form eps f xs andalso mapprox_form eps g xs
+| mapprox_form eps (@{term Disj} $ f $ g) xs = mapprox_form eps f xs orelse mapprox_form eps g xs
+| mapprox_form _ t _ = raise TERM ("mapprox_form", [t])
+
+fun dest_interpret_form (@{const "interpret_form"} $ b $ xs) = (b, xs)
+  | dest_interpret_form t = raise TERM ("dest_interpret_form", [t])
+
+fun optionT t = Type (@{type_name "option"}, [t])
+fun mk_Some t = Const (@{const_name "Some"}, t --> optionT t)
+
+fun random_float_list size xs seed =
+  fold (K (apsnd (random_float size) #-> (fn c => apfst (fn b => b::c)))) xs ([],seed)
+
+fun real_of_Float (@{code Float} (m, e)) =
+    real_of_man_exp (@{code integer_of_int} m) (@{code integer_of_int} e)
+
+fun is_True @{term True} = true
+  | is_True _ = false
+
+val postproc_form_eqs =
+  @{lemma
+    "real (Float 0 a) = 0"
+    "real (Float (numeral m) 0) = numeral m"
+    "real (Float 1 0) = 1"
+    "real (Float (- 1) 0) = - 1"
+    "real (Float 1 (numeral e)) = 2 ^ numeral e"
+    "real (Float 1 (- numeral e)) = 1 / 2 ^ numeral e"
+    "real (Float a 1) = a * 2"
+    "real (Float a (-1)) = a / 2"
+    "real (Float (- a) b) = - real (Float a b)"
+    "real (Float (numeral m) (numeral e)) = numeral m * 2 ^ (numeral e)"
+    "real (Float (numeral m) (- numeral e)) = numeral m / 2 ^ (numeral e)"
+    "- (c * d::real) = -c * d"
+    "- (c / d::real) = -c / d"
+    "- (0::real) = 0"
+    "int_of_integer (numeral k) = numeral k"
+    "int_of_integer (- numeral k) = - numeral k"
+    "int_of_integer 0 = 0"
+    "int_of_integer 1 = 1"
+    "int_of_integer (- 1) = - 1"
+    by auto
+  }
+
+fun rewrite_with ctxt thms = Simplifier.rewrite (put_simpset HOL_basic_ss ctxt addsimps thms)
+fun conv_term thy conv r = cterm_of thy r |> conv |> Thm.prop_of |> Logic.dest_equals |> snd
+
+fun approx_random ctxt prec eps frees e xs genuine_only size seed =
+  let
+    val (rs, seed') = random_float_list size xs seed
+    fun mk_approx_form e ts =
+      @{const "approx_form"} $
+        HOLogic.mk_number @{typ nat} prec $
+        e $
+        (HOLogic.mk_list @{typ "(float * float) option"}
+          (map (fn t => mk_Some @{typ "float * float"} $ HOLogic.mk_prod (t, t)) ts)) $
+        @{term "[] :: nat list"}
+  in
+    (if mapprox_form eps e (map (real_of_Float o fst) rs)
+    then
+      let
+        val ts = map (fn x => snd x ()) rs
+        val ts' = map
+          (AList.lookup op = (map dest_Free xs ~~ ts)
+            #> the_default Term.dummy
+            #> curry op $ @{term "real::float\<Rightarrow>_"}
+            #> conv_term (Proof_Context.theory_of ctxt) (rewrite_with ctxt postproc_form_eqs))
+          frees
+      in
+        if approximate ctxt (mk_approx_form e ts) |> is_True
+        then SOME (true, ts')
+        else (if genuine_only then NONE else SOME (false, ts'))
+      end
+    else NONE, seed')
+  end
+
+val preproc_form_eqs =
+  @{lemma
+    "(a::real) \<in> {b .. c} \<longleftrightarrow> b \<le> a \<and> a \<le> c"
+    "a = b \<longleftrightarrow> a \<le> b \<and> b \<le> a"
+    "(p \<longrightarrow> q) \<longleftrightarrow> \<not>p \<or> q"
+    "(p \<longleftrightarrow> q) \<longleftrightarrow> (p \<longrightarrow> q) \<and> (q \<longrightarrow> p)"
+    "\<not> (a < b) \<longleftrightarrow> b \<le> a"
+    "\<not> (a \<le> b) \<longleftrightarrow> b < a"
+    "\<not> (p \<and> q) \<longleftrightarrow> \<not> p \<or> \<not> q"
+    "\<not> (p \<or> q) \<longleftrightarrow> \<not> p \<and> \<not> q"
+    "\<not> \<not> q \<longleftrightarrow> q"
+    by auto
+  }
+
+fun reify_goal ctxt t =
+  HOLogic.mk_not t
+    |> conv_term (Proof_Context.theory_of ctxt) (rewrite_with ctxt preproc_form_eqs)
+    |> conv_term (Proof_Context.theory_of ctxt) (Reification.conv ctxt form_equations)
+    |> dest_interpret_form
+    ||> HOLogic.dest_list
+
+fun approximation_generator_raw ctxt t =
+  let
+    val iterations = Config.get ctxt Quickcheck.iterations
+    val prec = Config.get ctxt precision
+    val eps = Config.get ctxt epsilon
+    val cs = Config.get ctxt custom_seed
+    val seed = (Code_Numeral.natural_of_integer (cs + 1), Code_Numeral.natural_of_integer 1)
+    val run = if cs < 0
+      then (fn f => fn seed => (Random_Engine.run f, seed))
+      else (fn f => fn seed => f seed)
+    val frees = Term.add_frees t []
+    val (e, xs) = reify_goal ctxt t
+    fun single_tester b s =
+      approx_random ctxt prec eps frees e xs b s |> run
+    fun iterate _ _ 0 _ = NONE
+      | iterate genuine_only size j seed =
+        case single_tester genuine_only size seed of
+          (NONE, seed') => iterate genuine_only size (j - 1) seed'
+        | (SOME q, _) => SOME q
+  in
+    fn genuine_only => fn size => (iterate genuine_only size iterations seed, NONE)
+  end
+
+fun approximation_generator ctxt [(t, _)] =
+  (fn genuine_only =>
+    fn [_, size] =>
+      approximation_generator_raw ctxt t genuine_only
+        (Code_Numeral.natural_of_integer size))
+  | approximation_generator _ _ =
+      error "Quickcheck-approximation does not support type variables (or finite instantiations)"
+
+val test_goals =
+  Quickcheck_Common.generator_test_goal_terms
+    ("approximation", (fn _ => fn _ => false, approximation_generator))
+
+val active = Attrib.setup_config_bool @{binding quickcheck_approximation_active} (K false)
+
+val setup = Context.theory_map (Quickcheck.add_tester ("approximation", (active, test_goals)))
+
+end
diff -r 302104d8366b -r 87d4ce309e04 src/HOL/Decision_Procs/ex/Approximation_Ex.thy
--- a/src/HOL/Decision_Procs/ex/Approximation_Ex.thy	Wed Nov 12 18:18:38 2014 +0100
+++ b/src/HOL/Decision_Procs/ex/Approximation_Ex.thy	Wed Nov 12 19:30:56 2014 +0100
@@ -65,7 +65,7 @@
   by (approximation 10)
 
 lemma "3.2 \<le> x \<and> x \<le> 6.2 \<Longrightarrow> sin x \<le> 0"
-  by (approximation 9)
+  by (approximation 10)
 
 lemma
   defines "g \<equiv> 9.80665" and "v \<equiv> 128.61" and "d \<equiv> pi / 180"
diff -r 302104d8366b -r 87d4ce309e04 src/HOL/Decision_Procs/ex/Approximation_Quickcheck_Ex.thy
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/HOL/Decision_Procs/ex/Approximation_Quickcheck_Ex.thy	Wed Nov 12 19:30:56 2014 +0100
@@ -0,0 +1,39 @@
+theory Approximation_Quickcheck_Ex
+imports "../Approximation"
+begin
+
+lemma
+  fixes x::real and y::real
+  shows "sin x \<le> tan x"
+  using [[quickcheck_approximation_custom_seed = 1]]
+  quickcheck[approximation, expect=counterexample]
+  oops
+
+lemma "x \<le> y \<Longrightarrow> arctan y / y \<le> arctan x / x"
+  using [[quickcheck_approximation_custom_seed = 1]]
+  quickcheck[approximation, expect=counterexample]
+  oops
+
+lemma "0 < x \<Longrightarrow> x \<le> y \<Longrightarrow> arctan y / y \<le> arctan x / x"
+  using [[quickcheck_approximation_custom_seed = 1]]
+  quickcheck[approximation, expect=no_counterexample]
+  by (rule arctan_divide_mono)
+
+lemma
+  fixes x::real
+  shows "exp (exp x + exp y + sin x * sin y) - 0.4 > 0 \<or> 0.98 - sin x / (sin x * sin y + 2)^2 > 0"
+  using [[quickcheck_approximation_custom_seed = 1]]
+  quickcheck[approximation, expect=counterexample, size=10, iterations=1000, verbose]
+  oops
+
+lemma
+  fixes x::real
+  shows "x > 1 \<Longrightarrow> x \<le> 2 powr 20 * log 2 x + 1 \<and> (sin x)\<^sup>2 + (cos x)\<^sup>2 = 1"
+  using [[quickcheck_approximation_custom_seed = 1]]
+  using [[quickcheck_approximation_epsilon = 0.00000001]]
+    --\<open>avoids spurious counterexamples in approximate computation of @{term "(sin x)\<^sup>2 + (cos x)\<^sup>2"}
+      and therefore avoids expensive failing attempts for certification\<close>
+  quickcheck[approximation, expect=counterexample, size=20]
+  oops
+
+end
diff -r 302104d8366b -r 87d4ce309e04 src/HOL/Library/Float.thy
--- a/src/HOL/Library/Float.thy	Wed Nov 12 18:18:38 2014 +0100
+++ b/src/HOL/Library/Float.thy	Wed Nov 12 19:30:56 2014 +0100
@@ -139,6 +139,15 @@
   finally show ?thesis .
 qed
 
+lemma power_float[simp]: assumes "a \<in> float" shows "a ^ b \<in> float"
+proof -
+  from assms obtain m e::int where "a = m * 2 powr e"
+    by (auto simp: float_def)
+  thus ?thesis
+    by (auto intro!: floatI[where m="m^b" and e = "e*b"]
+      simp: power_mult_distrib powr_realpow[symmetric] powr_powr)
+qed
+
 lift_definition Float :: "int \<Rightarrow> int \<Rightarrow> float" is "\<lambda>(m::int) (e::int). m * 2 powr e" by simp
 declare Float.rep_eq[simp]
 
@@ -182,6 +191,9 @@
   proof qed (transfer, fastforce simp add: field_simps intro: mult_left_mono mult_right_mono)+
 end
 
+lemma Float_0_eq_0[simp]: "Float 0 e = 0"
+  by transfer simp
+
 lemma real_of_float_power[simp]: fixes f::float shows "real (f^n) = real f^n"
   by (induct n) simp_all
 
@@ -248,6 +260,46 @@
     and float_of_neg_numeral[simp]: "- numeral k = float_of (- numeral k)"
   unfolding real_of_float_eq by simp_all
 
+subsection {* Quickcheck *}
+
+instantiation float :: exhaustive
+begin
+
+definition exhaustive_float where
+  "exhaustive_float f d =
+    Quickcheck_Exhaustive.exhaustive (%x. Quickcheck_Exhaustive.exhaustive (%y. f (Float x y)) d) d"
+
+instance ..
+
+end
+
+definition (in term_syntax) [code_unfold]:
+  "valtermify_float x y = Code_Evaluation.valtermify Float {\<cdot>} x {\<cdot>} y"
+
+instantiation float :: full_exhaustive
+begin
+
+definition full_exhaustive_float where
+  "full_exhaustive_float f d =
+    Quickcheck_Exhaustive.full_exhaustive
+      (\<lambda>x. Quickcheck_Exhaustive.full_exhaustive (\<lambda>y. f (valtermify_float x y)) d) d"
+
+instance ..
+
+end
+
+instantiation float :: random
+begin
+
+definition "Quickcheck_Random.random i =
+  scomp (Quickcheck_Random.random (2 ^ nat_of_natural i))
+    (\<lambda>man. scomp (Quickcheck_Random.random i) (\<lambda>exp. Pair (valtermify_float man exp)))"
+
+instance ..
+
+end
+
+
 subsection {* Represent floats as unique mantissa and exponent *}
 
 lemma int_induct_abs[case_names less]:
@@ -435,13 +487,12 @@
   by transfer simp
 hide_fact (open) compute_float_one
 
-definition normfloat :: "float \<Rightarrow> float" where
-  [simp]: "normfloat x = x"
+lift_definition normfloat :: "float \<Rightarrow> float" is "\<lambda>x. x" .
+lemma normloat_id[simp]: "normfloat x = x" by transfer rule
 
 lemma compute_normfloat[code]: "normfloat (Float m e) =
   (if m mod 2 = 0 \<and> m \<noteq> 0 then normfloat (Float (m div 2) (e + 1))
                            else if m = 0 then 0 else Float m e)"
-  unfolding normfloat_def
   by transfer (auto simp add: powr_add zmod_eq_0_iff)
 hide_fact (open) compute_normfloat
 
@@ -510,7 +561,32 @@
   by transfer simp
 hide_fact (open) compute_float_eq
 
-subsection {* Rounding Real numbers *}
+
+subsection {* Lemmas for types @{typ real}, @{typ nat}, @{typ int}*}
+
+lemmas real_of_ints =
+  real_of_int_zero
+  real_of_one
+  real_of_int_add
+  real_of_int_minus
+  real_of_int_diff
+  real_of_int_mult
+  real_of_int_power
+  real_numeral
+lemmas real_of_nats =
+  real_of_nat_zero
+  real_of_nat_one
+  real_of_nat_1
+  real_of_nat_add
+  real_of_nat_mult
+  real_of_nat_power
+  real_of_nat_numeral
+
+lemmas int_of_reals = real_of_ints[symmetric]
+lemmas nat_of_reals = real_of_nats[symmetric]
+
+
+subsection {* Rounding Real Numbers *}
 
 definition round_down :: "int \<Rightarrow> real \<Rightarrow> real" where
   "round_down prec x = floor (x * 2 powr prec) * 2 powr -prec"
@@ -561,8 +637,85 @@
   by (simp add: powr_add powr_mult field_simps powr_divide2[symmetric])
     (simp add: powr_add[symmetric])
 
+lemma round_up_uminus_eq: "round_up p (-x) = - round_down p x"
+  and round_down_uminus_eq: "round_down p (-x) = - round_up p x"
+  by (auto simp: round_up_def round_down_def ceiling_def)
+
+lemma round_up_mono: "x \<le> y \<Longrightarrow> round_up p x \<le> round_up p y"
+  by (auto intro!: ceiling_mono simp: round_up_def)
+
+lemma round_up_le1:
+  assumes "x \<le> 1" "prec \<ge> 0"
+  shows "round_up prec x \<le> 1"
+proof -
+  have "real \<lceil>x * 2 powr prec\<rceil> \<le> real \<lceil>2 powr real prec\<rceil>"
+    using assms by (auto intro!: ceiling_mono)
+  also have "\<dots> = 2 powr prec" using assms by (auto simp: powr_int intro!: exI[where x="2^nat prec"])
+  finally show ?thesis
+    by (simp add: round_up_def) (simp add: powr_minus inverse_eq_divide)
+qed
+
+lemma round_up_less1:
+  assumes "x < 1 / 2" "p > 0"
+  shows "round_up p x < 1"
+proof -
+  have "x * 2 powr p < 1 / 2 * 2 powr p"
+    using assms by simp
+  also have "\<dots> \<le> 2 powr p - 1" using `p > 0`
+    by (auto simp: powr_divide2[symmetric] powr_int field_simps self_le_power)
+  finally show ?thesis using `p > 0`
+    by (simp add: round_up_def field_simps powr_minus powr_int ceiling_less_eq)
+qed
+
+lemma round_down_ge1:
+  assumes x: "x \<ge> 1"
+  assumes prec: "p \<ge> - log 2 x"
+  shows "1 \<le> round_down p x"
+proof cases
+  assume nonneg: "0 \<le> p"
+  have "2 powr p = real \<lfloor>2 powr real p\<rfloor>"
+    using nonneg by (auto simp: powr_int)
+  also have "\<dots> \<le> real \<lfloor>x * 2 powr p\<rfloor>"
+    using assms by (auto intro!: floor_mono)
+  finally show ?thesis
+    by (simp add: round_down_def) (simp add: powr_minus inverse_eq_divide)
+next
+  assume neg: "\<not> 0 \<le> p"
+  have "x = 2 powr (log 2 x)"
+    using x by simp
+  also have "2 powr (log 2 x) \<ge> 2 powr - p"
+    using prec by auto
+  finally have x_le: "x \<ge> 2 powr -p" .
+
+  from neg have "2 powr real p \<le> 2 powr 0"
+    by (intro powr_mono) auto
+  also have "\<dots> \<le> \<lfloor>2 powr 0\<rfloor>" by simp
+  also have "\<dots> \<le> \<lfloor>x * 2 powr real p\<rfloor>" unfolding real_of_int_le_iff
+    using x x_le by (intro floor_mono) (simp add: powr_minus_divide field_simps)
+  finally show ?thesis
+    using prec x
+    by (simp add: round_down_def powr_minus_divide pos_le_divide_eq)
+qed
+
+lemma round_up_le0: "x \<le> 0 \<Longrightarrow> round_up p x \<le> 0"
+  unfolding round_up_def
+  by (auto simp: field_simps mult_le_0_iff zero_le_mult_iff)
+
+
 subsection {* Rounding Floats *}
 
+definition div_twopow::"int \<Rightarrow> nat \<Rightarrow> int" where [simp]: "div_twopow x n = x div (2 ^ n)"
+
+definition mod_twopow::"int \<Rightarrow> nat \<Rightarrow> int" where [simp]: "mod_twopow x n = x mod (2 ^ n)"
+
+lemma compute_div_twopow[code]:
+  "div_twopow x n = (if x = 0 \<or> x = -1 \<or> n = 0 then x else div_twopow (x div 2) (n - 1))"
+  by (cases n) (auto simp: zdiv_zmult2_eq div_eq_minus1)
+
+lemma compute_mod_twopow[code]:
+  "mod_twopow x n = (if n = 0 then 0 else x mod 2 + 2 * mod_twopow (x div 2) (n - 1))"
+  by (cases n) (auto simp: zmod_zmult2_eq)
+
 lift_definition float_up :: "int \<Rightarrow> float \<Rightarrow> float" is round_up by simp
 declare float_up.rep_eq[simp]
 
@@ -599,7 +752,7 @@
 
 lemma compute_float_down[code]:
   "float_down p (Float m e) =
-    (if p + e < 0 then Float (m div 2^nat (-(p + e))) (-p) else Float m e)"
+    (if p + e < 0 then Float (div_twopow m (nat (-(p + e)))) (-p) else Float m e)"
 proof cases
   assume "p + e < 0"
   hence "real ((2::int) ^ nat (-(p + e))) = 2 powr (-(p + e))"
@@ -649,72 +802,10 @@
   floor_divide_eq_div dvd_neg_div del: divide_minus_left real_of_int_minus)
 
 lemma compute_float_up[code]:
-  "float_up p (Float m e) =
-    (let P = 2^nat (-(p + e)); r = m mod P in
-      if p + e < 0 then Float (m div P + (if r = 0 then 0 else 1)) (-p) else Float m e)"
-proof cases
-  assume "p + e < 0"
-  hence "real ((2::int) ^ nat (-(p + e))) = 2 powr (-(p + e))"
-    using powr_realpow[of 2 "nat (-(p + e))"] by simp
-  also have "... = 1 / 2 powr p / 2 powr e"
-  unfolding powr_minus_divide real_of_int_minus by (simp add: powr_add)
-  finally have twopow_rewrite:
-    "real ((2::int) ^ nat (- (p + e))) = 1 / 2 powr real p / 2 powr real e" .
-  with `p + e < 0` have powr_rewrite:
-    "2 powr real e * 2 powr real p = 1 / real ((2::int) ^ nat (- (p + e)))"
-    unfolding powr_divide2 by simp
-  show ?thesis
-  proof cases
-    assume "2^nat (-(p + e)) dvd m"
-    with `p + e < 0` twopow_rewrite show ?thesis
-      by transfer (auto simp: ac_simps round_up_def floor_divide_eq_div dvd_eq_mod_eq_0)
-  next
-    assume ndvd: "\<not> 2 ^ nat (- (p + e)) dvd m"
-    have one_div: "real m * (1 / real ((2::int) ^ nat (- (p + e)))) =
-      real m / real ((2::int) ^ nat (- (p + e)))"
-      by (simp add: field_simps)
-    have "real \<lceil>real m * (2 powr real e * 2 powr real p)\<rceil> =
-      real \<lfloor>real m * (2 powr real e * 2 powr real p)\<rfloor> + 1"
-      using ndvd unfolding powr_rewrite one_div
-      by (subst ceil_divide_floor_conv) (auto simp: field_simps)
-    thus ?thesis using `p + e < 0` twopow_rewrite
-      by transfer (auto simp: ac_simps round_up_def floor_divide_eq_div[symmetric])
-  qed
-next
-  assume "\<not> p + e < 0"
-  then have r1: "real e + real p = real (nat (e + p))" by simp
-  have r: "\<lceil>(m * 2 powr e) * 2 powr real p\<rceil> = (m * 2 powr e) * 2 powr real p"
-    by (auto simp add: ac_simps powr_add[symmetric] r1 powr_realpow
-      intro: exI[where x="m*2^nat (e+p)"])
-  then show ?thesis using `\<not> p + e < 0`
-    by transfer (simp add: round_up_def floor_divide_eq_div field_simps powr_add powr_minus)
-qed
+  "float_up p x = - float_down p (-x)"
+  by transfer (simp add: round_down_uminus_eq)
 hide_fact (open) compute_float_up
 
-lemmas real_of_ints =
-  real_of_int_zero
-  real_of_one
-  real_of_int_add
-  real_of_int_minus
-  real_of_int_diff
-  real_of_int_mult
-  real_of_int_power
-  real_numeral
-lemmas real_of_nats =
-  real_of_nat_zero
-  real_of_nat_one
-  real_of_nat_1
-  real_of_nat_add
-  real_of_nat_mult
-  real_of_nat_power
-
-lemmas int_of_reals = real_of_ints[symmetric]
-lemmas nat_of_reals = real_of_nats[symmetric]
-
-lemma two_real_int: "(2::real) = real (2::int)" by simp
-lemma two_real_nat: "(2::real) = real (2::nat)" by simp
-
-lemma mult_cong: "a = c ==> b = d ==> a*b = c*d" by simp
 
 subsection {* Compute bitlen of integers *}
 
@@ -752,7 +843,7 @@
     apply (simp add: powr_realpow[symmetric])
     using `x > 0` by simp
   finally show "x < 2 ^ nat (bitlen x)" using `x > 0`
-    by (simp add: bitlen_def ac_simps int_of_reals del: real_of_ints)
+    by (simp add: bitlen_def ac_simps)
 qed
 
 lemma bitlen_pow2[simp]:
@@ -851,13 +942,16 @@
     hence "?S \<le> real m" unfolding mult.assoc by auto
     hence "?S \<le> m" unfolding real_of_int_le_iff[symmetric] by auto
     from this bitlen_bounds[OF `0 < m`, THEN conjunct2]
-    have "nat (-e) < (nat (bitlen m))" unfolding power_strict_increasing_iff[OF `1 < 2`, symmetric] by (rule order_le_less_trans)
+    have "nat (-e) < (nat (bitlen m))" unfolding power_strict_increasing_iff[OF `1 < 2`, symmetric]
+      by (rule order_le_less_trans)
     hence "-e < bitlen m" using False by auto
     thus ?thesis by auto
   qed
 qed
 
-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"
+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"
 proof -
   let ?B = "2^nat(bitlen m - 1)"
 
@@ -876,10 +970,126 @@
   thus "real m / ?B < 2" by auto
 qed
 
-subsection {* Approximation of positive rationals *}
+subsection {* Truncating Real Numbers*}
+
+definition truncate_down::"nat \<Rightarrow> real \<Rightarrow> real" where
+  "truncate_down prec x = round_down (prec - \<lfloor>log 2 \<bar>x\<bar>\<rfloor> - 1) x"
+
+lemma truncate_down: "truncate_down prec x \<le> x"
+  using round_down by (simp add: truncate_down_def)
+
+lemma truncate_down_le: "x \<le> y \<Longrightarrow> truncate_down prec x \<le> y"
+  by (rule order_trans[OF truncate_down])
+
+lemma truncate_down_zero[simp]: "truncate_down prec 0 = 0"
+  by (simp add: truncate_down_def)
+
+lemma truncate_down_float[simp]: "truncate_down p x \<in> float"
+  by (auto simp: truncate_down_def)
+
+definition truncate_up::"nat \<Rightarrow> real \<Rightarrow> real" where
+  "truncate_up prec x = round_up (prec - \<lfloor>log 2 \<bar>x\<bar>\<rfloor> - 1) x"
+
+lemma truncate_up: "x \<le> truncate_up prec x"
+  using round_up by (simp add: truncate_up_def)
+
+lemma truncate_up_le: "x \<le> y \<Longrightarrow> x \<le> truncate_up prec y"
+  by (rule order_trans[OF _ truncate_up])
+
+lemma truncate_up_zero[simp]: "truncate_up prec 0 = 0"
+  by (simp add: truncate_up_def)
+
+lemma truncate_up_uminus_eq: "truncate_up prec (-x) = - truncate_down prec x"
+  and truncate_down_uminus_eq: "truncate_down prec (-x) = - truncate_up prec x"
+  by (auto simp: truncate_up_def round_up_def truncate_down_def round_down_def ceiling_def)
+
+lemma truncate_up_float[simp]: "truncate_up p x \<in> float"
+  by (auto simp: truncate_up_def)
+
+lemma mult_powr_eq: "0 < b \<Longrightarrow> b \<noteq> 1 \<Longrightarrow> 0 < x \<Longrightarrow> x * b powr y = b powr (y + log b x)"
+  by (simp_all add: powr_add)
+
+lemma truncate_down_pos:
+  assumes "x > 0" "p > 0"
+  shows "truncate_down p x > 0"
+proof -
+  have "0 \<le> log 2 x - real \<lfloor>log 2 x\<rfloor>"
+    by (simp add: algebra_simps)
+  from this assms
+  show ?thesis
+    by (auto simp: truncate_down_def round_down_def mult_powr_eq
+      intro!: ge_one_powr_ge_zero mult_pos_pos)
+qed
+
+lemma truncate_down_nonneg: "0 \<le> y \<Longrightarrow> 0 \<le> truncate_down prec y"
+  by (auto simp: truncate_down_def round_down_def)
+
+lemma truncate_down_ge1: "1 \<le> x \<Longrightarrow> 1 \<le> p \<Longrightarrow> 1 \<le> truncate_down p x"
+  by (auto simp: truncate_down_def algebra_simps intro!: round_down_ge1 add_mono)
+
+lemma truncate_up_nonpos: "x \<le> 0 \<Longrightarrow> truncate_up prec x \<le> 0"
+  by (auto simp: truncate_up_def round_up_def intro!: mult_nonpos_nonneg)
 
-lemma zdiv_zmult_twopow_eq: fixes a b::int shows "a div b div (2 ^ n) = a div (b * 2 ^ n)"
-by (simp add: zdiv_zmult2_eq)
+lemma truncate_up_le1:
+  assumes "x \<le> 1" "1 \<le> p" shows "truncate_up p x \<le> 1"
+proof -
+  {
+    assume "x \<le> 0"
+    with truncate_up_nonpos[OF this, of p] have ?thesis by simp
+  } moreover {
+    assume "x > 0"
+    hence le: "\<lfloor>log 2 \<bar>x\<bar>\<rfloor> \<le> 0"
+      using assms by (auto simp: log_less_iff)
+    from assms have "1 \<le> int p" by simp
+    from add_mono[OF this le]
+    have ?thesis using assms
+      by (simp add: truncate_up_def round_up_le1 add_mono)
+  } ultimately show ?thesis by arith
+qed
+
+subsection {* Truncating Floats*}
+
+lift_definition float_round_up :: "nat \<Rightarrow> float \<Rightarrow> float" is truncate_up
+  by (simp add: truncate_up_def)
+
+lemma float_round_up: "real x \<le> real (float_round_up prec x)"
+  using truncate_up by transfer simp
+
+lemma float_round_up_zero[simp]: "float_round_up prec 0 = 0"
+  by transfer simp
+
+lift_definition float_round_down :: "nat \<Rightarrow> float \<Rightarrow> float" is truncate_down
+  by (simp add: truncate_down_def)
+
+lemma float_round_down: "real (float_round_down prec x) \<le> real x"
+  using truncate_down by transfer simp
+
+lemma float_round_down_zero[simp]: "float_round_down prec 0 = 0"
+  by transfer simp
+
+lemmas float_round_up_le = order_trans[OF _ float_round_up]
+  and float_round_down_le = order_trans[OF float_round_down]
+
+lemma minus_float_round_up_eq: "- float_round_up prec x = float_round_down prec (- x)"
+  and minus_float_round_down_eq: "- float_round_down prec x = float_round_up prec (- x)"
+  by (transfer, simp add: truncate_down_uminus_eq truncate_up_uminus_eq)+
+
+lemma compute_float_round_down[code]:
+  "float_round_down prec (Float m e) = (let d = bitlen (abs m) - int prec in
+    if 0 < d then Float (div_twopow m (nat d)) (e + d)
+             else Float m e)"
+  using Float.compute_float_down[of "prec - bitlen \<bar>m\<bar> - e" m e, symmetric]
+  by transfer (simp add: field_simps abs_mult log_mult bitlen_def truncate_down_def
+    cong del: if_weak_cong)
+hide_fact (open) compute_float_round_down
+
+lemma compute_float_round_up[code]:
+  "float_round_up prec x = - float_round_down prec (-x)"
+  by transfer (simp add: truncate_down_uminus_eq)
+hide_fact (open) compute_float_round_up
+
+
+subsection {* Approximation of positive rationals *}
 
 lemma div_mult_twopow_eq: fixes a b::nat shows "a div ((2::nat) ^ n) div b = a div (b * 2 ^ n)"
   by (cases "b=0") (simp_all add: div_mult2_eq[symmetric] ac_simps)
@@ -901,7 +1111,7 @@
        l = rat_precision prec x y;
        d = if 0 \<le> l then x * 2^nat l div y else x div 2^nat (- l) div y
     in normfloat (Float d (- l)))"
-    unfolding div_mult_twopow_eq normfloat_def
+    unfolding div_mult_twopow_eq
     by transfer
        (simp add: round_down_def powr_int real_div_nat_eq_floor_of_divide field_simps Let_def
              del: two_powr_minus_int_float)
@@ -910,18 +1120,17 @@
 lift_definition rapprox_posrat :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float"
   is "\<lambda>prec (x::nat) (y::nat). round_up (rat_precision prec x y) (x / y)" by simp
 
-(* TODO: optimize using zmod_zmult2_eq, pdivmod ? *)
 lemma compute_rapprox_posrat[code]:
   fixes prec x y
+  notes divmod_int_mod_div[simp]
   defines "l \<equiv> rat_precision prec x y"
   shows "rapprox_posrat prec x y = (let
      l = l ;
      X = if 0 \<le> l then (x * 2^nat l, y) else (x, y * 2^nat(-l)) ;
-     d = fst X div snd X ;
-     m = fst X mod snd X
+     (d, m) = divmod_int (fst X) (snd X)
    in normfloat (Float (d + (if m = 0 \<or> y = 0 then 0 else 1)) (- l)))"
 proof (cases "y = 0")
-  assume "y = 0" thus ?thesis unfolding normfloat_def by transfer simp
+  assume "y = 0" thus ?thesis by transfer simp
 next
   assume "y \<noteq> 0"
   show ?thesis
@@ -932,7 +1141,6 @@
     moreover have "real x * 2 powr real l = real x'"
       by (simp add: powr_realpow[symmetric] `0 \<le> l` x'_def)
     ultimately show ?thesis
-      unfolding normfloat_def
       using ceil_divide_floor_conv[of y x'] powr_realpow[of 2 "nat l"] `0 \<le> l` `y \<noteq> 0`
         l_def[symmetric, THEN meta_eq_to_obj_eq]
       by transfer (auto simp add: floor_divide_eq_div [symmetric] round_up_def)
@@ -945,7 +1153,6 @@
       using `\<not> 0 \<le> l`
       by (simp add: powr_realpow[symmetric] powr_minus y'_def field_simps)
     ultimately show ?thesis
-      unfolding normfloat_def
       using ceil_divide_floor_conv[of y' x] `\<not> 0 \<le> l` `y' \<noteq> 0` `y \<noteq> 0`
         l_def[symmetric, THEN meta_eq_to_obj_eq]
       by transfer
@@ -966,41 +1173,9 @@
     using assms by (auto intro!: pos_add_strict simp add: field_simps rat_precision_def)
 qed
 
-lemma power_aux:
-  assumes "x > 0"
-  shows "(2::int) ^ nat (x - 1) \<le> 2 ^ nat x - 1"
-proof -
-  def y \<equiv> "nat (x - 1)"
-  moreover
-  have "(2::int) ^ y \<le> (2 ^ (y + 1)) - 1" by simp
-  ultimately show ?thesis using assms by simp
-qed
-
 lemma rapprox_posrat_less1:
-  assumes "0 \<le> x" and "0 < y" and "2 * x < y" and "0 < n"
-  shows "real (rapprox_posrat n x y) < 1"
-proof -
-  have powr1: "2 powr real (rat_precision n (int x) (int y)) =
-    2 ^ nat (rat_precision n (int x) (int y))" using rat_precision_pos[of x y n] assms
-    by (simp add: powr_realpow[symmetric])
-  have "x * 2 powr real (rat_precision n (int x) (int y)) / y = (x / y) *
-     2 powr real (rat_precision n (int x) (int y))" by simp
-  also have "... < (1 / 2) * 2 powr real (rat_precision n (int x) (int y))"
-    apply (rule mult_strict_right_mono) by (insert assms) auto
-  also have "\<dots> = 2 powr real (rat_precision n (int x) (int y) - 1)"
-    using powr_add [of 2 _ "- 1", simplified add_uminus_conv_diff] by (simp add: powr_minus)
-  also have "\<dots> = 2 ^ nat (rat_precision n (int x) (int y) - 1)"
-    using rat_precision_pos[of x y n] assms by (simp add: powr_realpow[symmetric])
-  also have "\<dots> \<le> 2 ^ nat (rat_precision n (int x) (int y)) - 1"
-    unfolding int_of_reals real_of_int_le_iff
-    using rat_precision_pos[OF assms] by (rule power_aux)
-  finally show ?thesis
-    apply (transfer fixing: n x y)
-    apply (simp add: round_up_def field_simps powr_minus powr1)
-    unfolding int_of_reals real_of_int_less_iff
-    apply (simp add: ceiling_less_eq)
-    done
-qed
+  shows "0 \<le> x \<Longrightarrow> 0 < y \<Longrightarrow> 2 * x < y \<Longrightarrow> 0 < n \<Longrightarrow> real (rapprox_posrat n x y) < 1"
+  by transfer (simp add: rat_precision_pos round_up_less1)
 
 lift_definition lapprox_rat :: "nat \<Rightarrow> int \<Rightarrow> int \<Rightarrow> float" is
   "\<lambda>prec (x::int) (y::int). round_down (rat_precision prec \<bar>x\<bar> \<bar>y\<bar>) (x / y)" by simp
@@ -1020,16 +1195,15 @@
 lift_definition rapprox_rat :: "nat \<Rightarrow> int \<Rightarrow> int \<Rightarrow> float" is
   "\<lambda>prec (x::int) (y::int). round_up (rat_precision prec \<bar>x\<bar> \<bar>y\<bar>) (x / y)" by simp
 
+lemma "rapprox_rat = rapprox_posrat"
+  by transfer auto
+
+lemma "lapprox_rat = lapprox_posrat"
+  by transfer auto
+
 lemma compute_rapprox_rat[code]:
-  "rapprox_rat prec x y =
-    (if y = 0 then 0
-    else if 0 \<le> x then
-      (if 0 < y then rapprox_posrat prec (nat x) (nat y)
-      else - (lapprox_posrat prec (nat x) (nat (-y))))
-      else (if 0 < y
-        then - (lapprox_posrat prec (nat (-x)) (nat y))
-        else rapprox_posrat prec (nat (-x)) (nat (-y))))"
-  by transfer (auto simp: round_up_def round_down_def ceiling_def ac_simps)
+  "rapprox_rat prec x y = - lapprox_rat prec (-x) y"
+  by transfer (simp add: round_down_uminus_eq)
 hide_fact (open) compute_rapprox_rat
 
 subsection {* Division *}
@@ -1063,22 +1237,467 @@
   by (simp add: real_divr_def)
 
 lemma compute_float_divr[code]:
-  "float_divr prec (Float m1 s1) (Float m2 s2) = rapprox_rat prec m1 m2 * Float 1 (s1 - s2)"
+  "float_divr prec x y = - float_divl prec (-x) y"
+  by transfer (simp add: real_divr_def real_divl_def round_down_uminus_eq)
+hide_fact (open) compute_float_divr
+
+
+subsection {* Approximate Power *}
+
+lemma div2_less_self[termination_simp]: fixes n::nat shows "odd n \<Longrightarrow> n div 2 < n"
+  by (simp add: odd_pos)
+
+fun power_down :: "nat \<Rightarrow> real \<Rightarrow> nat \<Rightarrow> real" where
+  "power_down p x 0 = 1"
+| "power_down p x (Suc n) =
+    (if odd n then truncate_down (Suc p) ((power_down p x (Suc n div 2))\<^sup>2) else truncate_down (Suc p) (x * power_down p x n))"
+
+fun power_up :: "nat \<Rightarrow> real \<Rightarrow> nat \<Rightarrow> real" where
+  "power_up p x 0 = 1"
+| "power_up p x (Suc n) =
+    (if odd n then truncate_up p ((power_up p x (Suc n div 2))\<^sup>2) else truncate_up p (x * power_up p x n))"
+
+lift_definition power_up_fl :: "nat \<Rightarrow> float \<Rightarrow> nat \<Rightarrow> float" is power_up
+  by (induct_tac rule: power_up.induct) simp_all
+
+lift_definition power_down_fl :: "nat \<Rightarrow> float \<Rightarrow> nat \<Rightarrow> float" is power_down
+  by (induct_tac rule: power_down.induct) simp_all
+
+lemma power_float_transfer[transfer_rule]:
+  "(rel_fun pcr_float (rel_fun op = pcr_float)) op ^ op ^"
+  unfolding power_def
+  by transfer_prover
+
+lemma compute_power_up_fl[code]:
+  "power_up_fl p x 0 = 1"
+  "power_up_fl p x (Suc n) =
+    (if odd n then float_round_up p ((power_up_fl p x (Suc n div 2))\<^sup>2) else float_round_up p (x * power_up_fl p x n))"
+  and compute_power_down_fl[code]:
+  "power_down_fl p x 0 = 1"
+  "power_down_fl p x (Suc n) =
+    (if odd n then float_round_down (Suc p) ((power_down_fl p x (Suc n div 2))\<^sup>2) else float_round_down (Suc p) (x * power_down_fl p x n))"
+  unfolding atomize_conj
+  by transfer simp
+
+lemma power_down_pos: "0 < x \<Longrightarrow> 0 < power_down p x n"
+  by (induct p x n rule: power_down.induct)
+    (auto simp del: odd_Suc_div_two intro!: truncate_down_pos)
+
+lemma power_down_nonneg: "0 \<le> x \<Longrightarrow> 0 \<le> power_down p x n"
+  by (induct p x n rule: power_down.induct)
+    (auto simp del: odd_Suc_div_two intro!: truncate_down_nonneg mult_nonneg_nonneg)
+
+lemma power_down: "0 \<le> x \<Longrightarrow> power_down p x n \<le> x ^ n"
+proof (induct p x n rule: power_down.induct)
+  case (2 p x n)
+  {
+    assume "odd n"
+    hence "(power_down p x (Suc n div 2)) ^ 2 \<le> (x ^ (Suc n div 2)) ^ 2"
+      using 2
+      by (auto intro: power_mono power_down_nonneg simp del: odd_Suc_div_two)
+    also have "\<dots> = x ^ (Suc n div 2 * 2)"
+      by (simp add: power_mult[symmetric])
+    also have "Suc n div 2 * 2 = Suc n"
+      using `odd n` by presburger
+    finally have ?case
+      using `odd n`
+      by (auto intro!: truncate_down_le simp del: odd_Suc_div_two)
+  } thus ?case
+    by (auto intro!: truncate_down_le mult_left_mono 2 mult_nonneg_nonneg power_down_nonneg)
+qed simp
+
+lemma power_up: "0 \<le> x \<Longrightarrow> x ^ n \<le> power_up p x n"
+proof (induct p x n rule: power_up.induct)
+  case (2 p x n)
+  {
+    assume "odd n"
+    hence "Suc n = Suc n div 2 * 2"
+      using `odd n` even_Suc by presburger
+    hence "x ^ Suc n \<le> (x ^ (Suc n div 2))\<^sup>2"
+      by (simp add: power_mult[symmetric])
+    also have "\<dots> \<le> (power_up p x (Suc n div 2))\<^sup>2"
+      using 2 `odd n`
+      by (auto intro: power_mono simp del: odd_Suc_div_two )
+    finally have ?case
+      using `odd n`
+      by (auto intro!: truncate_up_le simp del: odd_Suc_div_two )
+  } thus ?case
+    by (auto intro!: truncate_up_le mult_left_mono 2)
+qed simp
+
+lemmas power_up_le = order_trans[OF _ power_up]
+  and power_up_less = less_le_trans[OF _ power_up]
+  and power_down_le = order_trans[OF power_down]
+
+lemma power_down_fl: "0 \<le> x \<Longrightarrow> power_down_fl p x n \<le> x ^ n"
+  by transfer (rule power_down)
+
+lemma power_up_fl: "0 \<le> x \<Longrightarrow> x ^ n \<le> power_up_fl p x n"
+  by transfer (rule power_up)
+
+lemma real_power_up_fl: "real (power_up_fl p x n) = power_up p x n"
+  by transfer simp
+
+lemma real_power_down_fl: "real (power_down_fl p x n) = power_down p x n"
+  by transfer simp
+
+
+subsection {* Approximate Addition *}
+
+definition "plus_down prec x y = truncate_down prec (x + y)"
+
+definition "plus_up prec x y = truncate_up prec (x + y)"
+
+lemma float_plus_down_float[intro, simp]: "x \<in> float \<Longrightarrow> y \<in> float \<Longrightarrow> plus_down p x y \<in> float"
+  by (simp add: plus_down_def)
+
+lemma float_plus_up_float[intro, simp]: "x \<in> float \<Longrightarrow> y \<in> float \<Longrightarrow> plus_up p x y \<in> float"
+  by (simp add: plus_up_def)
+
+lift_definition float_plus_down::"nat \<Rightarrow> float \<Rightarrow> float \<Rightarrow> float" is plus_down ..
+
+lift_definition float_plus_up::"nat \<Rightarrow> float \<Rightarrow> float \<Rightarrow> float" is plus_up ..
+
+lemma plus_down: "plus_down prec x y \<le> x + y"
+  and plus_up: "x + y \<le> plus_up prec x y"
+  by (auto simp: plus_down_def truncate_down plus_up_def truncate_up)
+
+lemma float_plus_down: "real (float_plus_down prec x y) \<le> x + y"
+  and float_plus_up: "x + y \<le> real (float_plus_up prec x y)"
+  by (transfer, rule plus_down plus_up)+
+
+lemmas plus_down_le = order_trans[OF plus_down]
+  and plus_up_le = order_trans[OF _ plus_up]
+  and float_plus_down_le = order_trans[OF float_plus_down]
+  and float_plus_up_le = order_trans[OF _ float_plus_up]
+
+lemma compute_plus_up[code]: "plus_up p x y = - plus_down p (-x) (-y)"
+  using truncate_down_uminus_eq[of p "x + y"]
+  by (auto simp: plus_down_def plus_up_def)
+
+lemma
+  truncate_down_log2_eqI:
+  assumes "\<lfloor>log 2 \<bar>x\<bar>\<rfloor> = \<lfloor>log 2 \<bar>y\<bar>\<rfloor>"
+  assumes "\<lfloor>x * 2 powr (p - \<lfloor>log 2 \<bar>x\<bar>\<rfloor> - 1)\<rfloor> = \<lfloor>y * 2 powr (p - \<lfloor>log 2 \<bar>x\<bar>\<rfloor> - 1)\<rfloor>"
+  shows "truncate_down p x = truncate_down p y"
+  using assms by (auto simp: truncate_down_def round_down_def)
+
+lemma bitlen_eq_zero_iff: "bitlen x = 0 \<longleftrightarrow> x \<le> 0"
+  by (clarsimp simp add: bitlen_def)
+    (metis Float.compute_bitlen add.commute bitlen_def bitlen_nonneg less_add_same_cancel2 not_less
+      zero_less_one)
+
+lemma
+  sum_neq_zeroI:
+  fixes a k::real
+  shows "abs a \<ge> k \<Longrightarrow> abs b < k \<Longrightarrow> a + b \<noteq> 0"
+    and "abs a > k \<Longrightarrow> abs b \<le> k \<Longrightarrow> a + b \<noteq> 0"
+  by auto
+
+lemma
+  abs_real_le_2_powr_bitlen[simp]:
+  "\<bar>real m2\<bar> < 2 powr real (bitlen \<bar>m2\<bar>)"
 proof cases
-  let ?f1 = "real m1 * 2 powr real s1" and ?f2 = "real m2 * 2 powr real s2"
-  let ?m = "real m1 / real m2" and ?s = "2 powr real (s1 - s2)"
-  assume not_0: "m1 \<noteq> 0 \<and> m2 \<noteq> 0"
-  then have eq2: "(int prec + \<lfloor>log 2 \<bar>?f2\<bar>\<rfloor> - \<lfloor>log 2 \<bar>?f1\<bar>\<rfloor>) = rat_precision prec \<bar>m1\<bar> \<bar>m2\<bar> + (s2 - s1)"
-    by (simp add: abs_mult log_mult rat_precision_def bitlen_def)
-  have eq1: "real m1 * 2 powr real s1 / (real m2 * 2 powr real s2) = ?m * ?s"
-    by (simp add: field_simps powr_divide2[symmetric])
+  assume "m2 \<noteq> 0"
+  hence "\<bar>m2\<bar> < 2 ^ nat (bitlen \<bar>m2\<bar>)"
+    using bitlen_bounds[of "\<bar>m2\<bar>"]
+    by (auto simp: powr_add bitlen_nonneg)
+  thus ?thesis
+    by (simp add: powr_int bitlen_nonneg real_of_int_less_iff[symmetric])
+qed simp
+
+lemma floor_sum_times_2_powr_sgn_eq:
+  fixes ai p q::int
+  and a b::real
+  assumes "a * 2 powr p = ai"
+  assumes b_le_1: "abs (b * 2 powr (p + 1)) \<le> 1"
+  assumes leqp: "q \<le> p"
+  shows "\<lfloor>(a + b) * 2 powr q\<rfloor> = \<lfloor>(2 * ai + sgn b) * 2 powr (q - p - 1)\<rfloor>"
+proof -
+  {
+    assume "b = 0"
+    hence ?thesis
+      by (simp add: assms(1)[symmetric] powr_add[symmetric] algebra_simps powr_mult_base)
+  } moreover {
+    assume "b > 0"
+    hence "b * 2 powr p < abs (b * 2 powr (p + 1))" by simp
+    also note b_le_1
+    finally have b_less_1: "b * 2 powr real p < 1" .
+
+    from b_less_1 `b > 0` have floor_eq: "\<lfloor>b * 2 powr real p\<rfloor> = 0" "\<lfloor>sgn b / 2\<rfloor> = 0"
+      by (simp_all add: floor_eq_iff)
+
+    have "\<lfloor>(a + b) * 2 powr q\<rfloor> = \<lfloor>(a + b) * 2 powr p * 2 powr (q - p)\<rfloor>"
+      by (simp add: algebra_simps powr_realpow[symmetric] powr_add[symmetric])
+    also have "\<dots> = \<lfloor>(ai + b * 2 powr p) * 2 powr (q - p)\<rfloor>"
+      by (simp add: assms algebra_simps)
+    also have "\<dots> = \<lfloor>(ai + b * 2 powr p) / real ((2::int) ^ nat (p - q))\<rfloor>"
+      using assms
+      by (simp add: algebra_simps powr_realpow[symmetric] divide_powr_uminus powr_add[symmetric])
+    also have "\<dots> = \<lfloor>ai / real ((2::int) ^ nat (p - q))\<rfloor>"
+      by (simp del: real_of_int_power add: floor_divide_real_eq_div floor_eq)
+    finally have "\<lfloor>(a + b) * 2 powr real q\<rfloor> = \<lfloor>real ai / real ((2::int) ^ nat (p - q))\<rfloor>" .
+    moreover
+    {
+      have "\<lfloor>(2 * ai + sgn b) * 2 powr (real (q - p) - 1)\<rfloor> = \<lfloor>(ai + sgn b / 2) * 2 powr (q - p)\<rfloor>"
+        by (subst powr_divide2[symmetric]) (simp add: field_simps)
+      also have "\<dots> = \<lfloor>(ai + sgn b / 2) / real ((2::int) ^ nat (p - q))\<rfloor>"
+        using leqp by (simp add: powr_realpow[symmetric] powr_divide2[symmetric])
+      also have "\<dots> = \<lfloor>ai / real ((2::int) ^ nat (p - q))\<rfloor>"
+        by (simp del: real_of_int_power add: floor_divide_real_eq_div floor_eq)
+      finally
+      have "\<lfloor>(2 * ai + (sgn b)) * 2 powr (real (q - p) - 1)\<rfloor> =
+          \<lfloor>real ai / real ((2::int) ^ nat (p - q))\<rfloor>"
+        .
+    } ultimately have ?thesis by simp
+  } moreover {
+    assume "\<not> 0 \<le> b"
+    hence "0 > b" by simp
+    hence floor_eq: "\<lfloor>b * 2 powr (real p + 1)\<rfloor> = -1"
+      using b_le_1
+      by (auto simp: floor_eq_iff algebra_simps pos_divide_le_eq[symmetric] abs_if divide_powr_uminus
+        intro!: mult_neg_pos split: split_if_asm)
+    have "\<lfloor>(a + b) * 2 powr q\<rfloor> = \<lfloor>(2*a + 2*b) * 2 powr p * 2 powr (q - p - 1)\<rfloor>"
+      by (simp add: algebra_simps powr_realpow[symmetric] powr_add[symmetric] powr_mult_base)
+    also have "\<dots> = \<lfloor>(2 * (a * 2 powr p) + 2 * b * 2 powr p) * 2 powr (q - p - 1)\<rfloor>"
+      by (simp add: algebra_simps)
+    also have "\<dots> = \<lfloor>(2 * ai + b * 2 powr (p + 1)) / 2 powr (1 - q + p)\<rfloor>"
+      using assms by (simp add: algebra_simps powr_mult_base divide_powr_uminus)
+    also have "\<dots> = \<lfloor>(2 * ai + b * 2 powr (p + 1)) / real ((2::int) ^ nat (p - q + 1))\<rfloor>"
+      using assms by (simp add: algebra_simps powr_realpow[symmetric])
+    also have "\<dots> = \<lfloor>(2 * ai - 1) / real ((2::int) ^ nat (p - q + 1))\<rfloor>"
+      using `b < 0` assms
+      by (simp add: floor_divide_eq_div floor_eq floor_divide_real_eq_div
+        del: real_of_int_mult real_of_int_power real_of_int_diff)
+    also have "\<dots> = \<lfloor>(2 * ai - 1) * 2 powr (q - p - 1)\<rfloor>"
+      using assms by (simp add: algebra_simps divide_powr_uminus powr_realpow[symmetric])
+    finally have ?thesis using `b < 0` by simp
+  } ultimately show ?thesis by arith
+qed
+
+lemma
+  log2_abs_int_add_less_half_sgn_eq:
+  fixes ai::int and b::real
+  assumes "abs b \<le> 1/2" "ai \<noteq> 0"
+  shows "\<lfloor>log 2 \<bar>real ai + b\<bar>\<rfloor> = \<lfloor>log 2 \<bar>ai + sgn b / 2\<bar>\<rfloor>"
+proof cases
+  assume "b = 0" thus ?thesis by simp
+next
+  assume "b \<noteq> 0"
+  def k \<equiv> "\<lfloor>log 2 \<bar>ai\<bar>\<rfloor>"
+  hence "\<lfloor>log 2 \<bar>ai\<bar>\<rfloor> = k" by simp
+  hence k: "2 powr k \<le> \<bar>ai\<bar>" "\<bar>ai\<bar> < 2 powr (k + 1)"
+    by (simp_all add: floor_log_eq_powr_iff `ai \<noteq> 0`)
+  have "k \<ge> 0"
+    using assms by (auto simp: k_def)
+  def r \<equiv> "\<bar>ai\<bar> - 2 ^ nat k"
+  have r: "0 \<le> r" "r < 2 powr k"
+    using `k \<ge> 0` k
+    by (auto simp: r_def k_def algebra_simps powr_add abs_if powr_int)
+  hence "r \<le> (2::int) ^ nat k - 1"
+    using `k \<ge> 0` by (auto simp: powr_int)
+  from this[simplified real_of_int_le_iff[symmetric]] `0 \<le> k`
+  have r_le: "r \<le> 2 powr k - 1"
+    by (auto simp: algebra_simps powr_int simp del: real_of_int_le_iff)
+
+  have "\<bar>ai\<bar> = 2 powr k + r"
+    using `k \<ge> 0` by (auto simp: k_def r_def powr_realpow[symmetric])
+
+  have pos: "\<And>b::real. abs b < 1 \<Longrightarrow> 0 < 2 powr k + (r + b)"
+    using `0 \<le> k` `ai \<noteq> 0`
+    by (auto simp add: r_def powr_realpow[symmetric] abs_if sgn_if algebra_simps
+      split: split_if_asm)
+  have less: "\<bar>sgn ai * b\<bar> < 1"
+    and less': "\<bar>sgn (sgn ai * b) / 2\<bar> < 1"
+    using `abs b \<le> _` by (auto simp: abs_if sgn_if split: split_if_asm)
+
+  have floor_eq: "\<And>b::real. abs b \<le> 1 / 2 \<Longrightarrow>
+      \<lfloor>log 2 (1 + (r + b) / 2 powr k)\<rfloor> = (if r = 0 \<and> b < 0 then -1 else 0)"
+    using `k \<ge> 0` r r_le
+    by (auto simp: floor_log_eq_powr_iff powr_minus_divide field_simps sgn_if)
 
-  show ?thesis
-    using not_0
-    by (transfer fixing: m1 s1 m2 s2 prec) (unfold eq1 eq2 round_up_shift real_divr_def,
-      simp add: field_simps)
-qed (transfer, auto simp: real_divr_def)
-hide_fact (open) compute_float_divr
+  from `real \<bar>ai\<bar> = _` have "\<bar>ai + b\<bar> = 2 powr k + (r + sgn ai * b)"
+    using `abs b <= _` `0 \<le> k` r
+    by (auto simp add: sgn_if abs_if)
+  also have "\<lfloor>log 2 \<dots>\<rfloor> = \<lfloor>log 2 (2 powr k + r + sgn (sgn ai * b) / 2)\<rfloor>"
+  proof -
+    have "2 powr k + (r + (sgn ai) * b) = 2 powr k * (1 + (r + sgn ai * b) / 2 powr k)"
+      by (simp add: field_simps)
+    also have "\<lfloor>log 2 \<dots>\<rfloor> = k + \<lfloor>log 2 (1 + (r + sgn ai * b) / 2 powr k)\<rfloor>"
+      using pos[OF less]
+      by (subst log_mult) (simp_all add: log_mult powr_mult field_simps)
+    also
+    let ?if = "if r = 0 \<and> sgn ai * b < 0 then -1 else 0"
+    have "\<lfloor>log 2 (1 + (r + sgn ai * b) / 2 powr k)\<rfloor> = ?if"
+      using `abs b <= _`
+      by (intro floor_eq) (auto simp: abs_mult sgn_if)
+    also
+    have "\<dots> = \<lfloor>log 2 (1 + (r + sgn (sgn ai * b) / 2) / 2 powr k)\<rfloor>"
+      by (subst floor_eq) (auto simp: sgn_if)
+    also have "k + \<dots> = \<lfloor>log 2 (2 powr k * (1 + (r + sgn (sgn ai * b) / 2) / 2 powr k))\<rfloor>"
+      unfolding floor_add2[symmetric]
+      using pos[OF less'] `abs b \<le> _`
+      by (simp add: field_simps add_log_eq_powr)
+    also have "2 powr k * (1 + (r + sgn (sgn ai * b) / 2) / 2 powr k) =
+        2 powr k + r + sgn (sgn ai * b) / 2"
+      by (simp add: sgn_if field_simps)
+    finally show ?thesis .
+  qed
+  also have "2 powr k + r + sgn (sgn ai * b) / 2 = \<bar>ai + sgn b / 2\<bar>"
+    unfolding `real \<bar>ai\<bar> = _`[symmetric] using `ai \<noteq> 0`
+    by (auto simp: abs_if sgn_if algebra_simps)
+  finally show ?thesis .
+qed
+
+lemma compute_far_float_plus_down:
+  fixes m1 e1 m2 e2::int and p::nat
+  defines "k1 \<equiv> p - nat (bitlen \<bar>m1\<bar>)"
+  assumes H: "bitlen \<bar>m2\<bar> \<le> e1 - e2 - k1 - 2" "m1 \<noteq> 0" "m2 \<noteq> 0" "e1 \<ge> e2"
+  shows "float_plus_down p (Float m1 e1) (Float m2 e2) =
+    float_round_down p (Float (m1 * 2 ^ (Suc (Suc k1)) + sgn m2) (e1 - int k1 - 2))"
+proof -
+  let ?a = "real (Float m1 e1)"
+  let ?b = "real (Float m2 e2)"
+  let ?sum = "?a + ?b"
+  let ?shift = "real e2 - real e1 + real k1 + 1"
+  let ?m1 = "m1 * 2 ^ Suc k1"
+  let ?m2 = "m2 * 2 powr ?shift"
+  let ?m2' = "sgn m2 / 2"
+  let ?e = "e1 - int k1 - 1"
+
+  have sum_eq: "?sum = (?m1 + ?m2) * 2 powr ?e"
+    by (auto simp: powr_add[symmetric] powr_mult[symmetric] algebra_simps
+      powr_realpow[symmetric] powr_mult_base)
+
+  have "\<bar>?m2\<bar> * 2 < 2 powr (bitlen \<bar>m2\<bar> + ?shift + 1)"
+    by (auto simp: field_simps powr_add powr_mult_base powr_numeral powr_divide2[symmetric] abs_mult)
+  also have "\<dots> \<le> 2 powr 0"
+    using H by (intro powr_mono) auto
+  finally have abs_m2_less_half: "\<bar>?m2\<bar> < 1 / 2"
+    by simp
+
+  hence "\<bar>real m2\<bar> < 2 powr -(?shift + 1)"
+    unfolding powr_minus_divide by (auto simp: bitlen_def field_simps powr_mult_base abs_mult)
+  also have "\<dots> \<le> 2 powr real (e1 - e2 - 2)"
+    by simp
+  finally have b_less_quarter: "\<bar>?b\<bar> < 1/4 * 2 powr real e1"
+    by (simp add: powr_add field_simps powr_divide2[symmetric] powr_numeral abs_mult)
+  also have "1/4 < \<bar>real m1\<bar> / 2" using `m1 \<noteq> 0` by simp
+  finally have b_less_half_a: "\<bar>?b\<bar> < 1/2 * \<bar>?a\<bar>"
+    by (simp add: algebra_simps powr_mult_base abs_mult)
+  hence a_half_less_sum: "\<bar>?a\<bar> / 2 < \<bar>?sum\<bar>"
+    by (auto simp: field_simps abs_if split: split_if_asm)
+
+  from b_less_half_a have "\<bar>?b\<bar> < \<bar>?a\<bar>" "\<bar>?b\<bar> \<le> \<bar>?a\<bar>"
+    by simp_all
+
+  have "\<bar>real (Float m1 e1)\<bar> \<ge> 1/4 * 2 powr real e1"
+    using `m1 \<noteq> 0`
+    by (auto simp: powr_add powr_int bitlen_nonneg divide_right_mono abs_mult)
+  hence "?sum \<noteq> 0" using b_less_quarter
+    by (rule sum_neq_zeroI)
+  hence "?m1 + ?m2 \<noteq> 0"
+    unfolding sum_eq by (simp add: abs_mult zero_less_mult_iff)
+
+  have "\<bar>real ?m1\<bar> \<ge> 2 ^ Suc k1" "\<bar>?m2'\<bar> < 2 ^ Suc k1"
+    using `m1 \<noteq> 0` `m2 \<noteq> 0` by (auto simp: sgn_if less_1_mult abs_mult simp del: power.simps)
+  hence sum'_nz: "?m1 + ?m2' \<noteq> 0"
+    by (intro sum_neq_zeroI)
+
+  have "\<lfloor>log 2 \<bar>real (Float m1 e1) + real (Float m2 e2)\<bar>\<rfloor> = \<lfloor>log 2 \<bar>?m1 + ?m2\<bar>\<rfloor> + ?e"
+    using `?m1 + ?m2 \<noteq> 0`
+    unfolding floor_add[symmetric] sum_eq
+    by (simp add: abs_mult log_mult)
+  also have "\<lfloor>log 2 \<bar>?m1 + ?m2\<bar>\<rfloor> = \<lfloor>log 2 \<bar>?m1 + sgn (real m2 * 2 powr ?shift) / 2\<bar>\<rfloor>"
+    using abs_m2_less_half `m1 \<noteq> 0`
+    by (intro log2_abs_int_add_less_half_sgn_eq) (auto simp: abs_mult)
+  also have "sgn (real m2 * 2 powr ?shift) = sgn m2"
+    by (auto simp: sgn_if zero_less_mult_iff less_not_sym)
+  also
+  have "\<bar>?m1 + ?m2'\<bar> * 2 powr ?e = \<bar>?m1 * 2 + sgn m2\<bar> * 2 powr (?e - 1)"
+    by (auto simp: field_simps powr_minus[symmetric] powr_divide2[symmetric] powr_mult_base)
+  hence "\<lfloor>log 2 \<bar>?m1 + ?m2'\<bar>\<rfloor> + ?e = \<lfloor>log 2 \<bar>real (Float (?m1 * 2 + sgn m2) (?e - 1))\<bar>\<rfloor>"
+    using `?m1 + ?m2' \<noteq> 0`
+    unfolding floor_add[symmetric]
+    by (simp add: log_add_eq_powr abs_mult_pos)
+  finally
+  have "\<lfloor>log 2 \<bar>?sum\<bar>\<rfloor> = \<lfloor>log 2 \<bar>real (Float (?m1*2 + sgn m2) (?e - 1))\<bar>\<rfloor>" .
+  hence "plus_down p (Float m1 e1) (Float m2 e2) =
+      truncate_down p (Float (?m1*2 + sgn m2) (?e - 1))"
+    unfolding plus_down_def
+  proof (rule truncate_down_log2_eqI)
+    let ?f = "(int p - \<lfloor>log 2 \<bar>real (Float m1 e1) + real (Float m2 e2)\<bar>\<rfloor> - 1)"
+    let ?ai = "m1 * 2 ^ (Suc k1)"
+    have "\<lfloor>(?a + ?b) * 2 powr real ?f\<rfloor> = \<lfloor>(real (2 * ?ai) + sgn ?b) * 2 powr real (?f - - ?e - 1)\<rfloor>"
+    proof (rule floor_sum_times_2_powr_sgn_eq)
+      show "?a * 2 powr real (-?e) = real ?ai"
+        by (simp add: powr_add powr_realpow[symmetric] powr_divide2[symmetric])
+      show "\<bar>?b * 2 powr real (-?e + 1)\<bar> \<le> 1"
+        using abs_m2_less_half
+        by (simp add: abs_mult powr_add[symmetric] algebra_simps powr_mult_base)
+    next
+      have "e1 + \<lfloor>log 2 \<bar>real m1\<bar>\<rfloor> - 1 = \<lfloor>log 2 \<bar>?a\<bar>\<rfloor> - 1"
+        using `m1 \<noteq> 0`
+        by (simp add: floor_add2[symmetric] algebra_simps log_mult abs_mult del: floor_add2)
+      also have "\<dots> \<le> \<lfloor>log 2 \<bar>?a + ?b\<bar>\<rfloor>"
+        using a_half_less_sum `m1 \<noteq> 0` `?sum \<noteq> 0`
+        unfolding floor_subtract[symmetric]
+        by (auto simp add: log_minus_eq_powr powr_minus_divide
+          intro!: floor_mono)
+      finally
+      have "int p - \<lfloor>log 2 \<bar>?a + ?b\<bar>\<rfloor> \<le> p - (bitlen \<bar>m1\<bar>) - e1 + 2"
+        by (auto simp: algebra_simps bitlen_def `m1 \<noteq> 0`)
+      also have "\<dots> \<le> 1 - ?e"
+        using bitlen_nonneg[of "\<bar>m1\<bar>"] by (simp add: k1_def)
+      finally show "?f \<le> - ?e" by simp
+    qed
+    also have "sgn ?b = sgn m2"
+      using powr_gt_zero[of 2 e2]
+      by (auto simp add: sgn_if zero_less_mult_iff simp del: powr_gt_zero)
+    also have "\<lfloor>(real (2 * ?m1) + real (sgn m2)) * 2 powr real (?f - - ?e - 1)\<rfloor> =
+        \<lfloor>Float (?m1 * 2 + sgn m2) (?e - 1) * 2 powr ?f\<rfloor>"
+      by (simp add: powr_add[symmetric] algebra_simps powr_realpow[symmetric])
+    finally
+    show "\<lfloor>(?a + ?b) * 2 powr ?f\<rfloor> = \<lfloor>real (Float (?m1 * 2 + sgn m2) (?e - 1)) * 2 powr ?f\<rfloor>" .
+  qed
+  thus ?thesis
+    by transfer (simp add: plus_down_def ac_simps Let_def)
+qed
+
+lemma compute_float_plus_down_naive[code]: "float_plus_down p x y = float_round_down p (x + y)"
+  by transfer (auto simp: plus_down_def)
+
+lemma compute_float_plus_down[code]:
+  fixes p::nat and m1 e1 m2 e2::int
+  shows "float_plus_down p (Float m1 e1) (Float m2 e2) =
+    (if m1 = 0 then float_round_down p (Float m2 e2)
+    else if m2 = 0 then float_round_down p (Float m1 e1)
+    else (if e1 \<ge> e2 then
+      (let
+        k1 = p - nat (bitlen \<bar>m1\<bar>)
+      in
+        if bitlen \<bar>m2\<bar> > e1 - e2 - k1 - 2 then float_round_down p ((Float m1 e1) + (Float m2 e2))
+        else float_round_down p (Float (m1 * 2 ^ (Suc (Suc k1)) + sgn m2) (e1 - int k1 - 2)))
+    else float_plus_down p (Float m2 e2) (Float m1 e1)))"
+proof -
+  {
+    assume H: "bitlen \<bar>m2\<bar> \<le> e1 - e2 - (p - nat (bitlen \<bar>m1\<bar>)) - 2" "m1 \<noteq> 0" "m2 \<noteq> 0" "e1 \<ge> e2"
+    note compute_far_float_plus_down[OF H]
+  }
+  thus ?thesis
+    by transfer (simp add: Let_def plus_down_def ac_simps)
+qed
+hide_fact (open) compute_far_float_plus_down
+hide_fact (open) compute_float_plus_down
+
+lemma compute_float_plus_up[code]: "float_plus_up p x y = - float_plus_down p (-x) (-y)"
+  using truncate_down_uminus_eq[of p "x + y"]
+  by transfer (simp add: plus_down_def plus_up_def ac_simps)
+hide_fact (open) compute_float_plus_up
+
+lemma mantissa_zero[simp]: "mantissa 0 = 0"
+by (metis mantissa_0 zero_float.abs_eq)
+
 
 subsection {* Lemmas needed by Approximate *}
 
@@ -1113,12 +1732,9 @@
 
 lemma lapprox_rat_nonneg:
   fixes n x y
-  defines "p \<equiv> int n - ((bitlen \<bar>x\<bar>) - (bitlen \<bar>y\<bar>))"
-  assumes "0 \<le> x" and "0 < y"
+  assumes "0 \<le> x" and "0 \<le> y"
   shows "0 \<le> real (lapprox_rat n x y)"
-using assms unfolding lapprox_rat_def p_def[symmetric] round_down_def real_of_int_minus[symmetric]
-   powr_int[of 2, simplified]
-  by auto
+  using assms by (auto simp: lapprox_rat_def simp: round_down_nonneg)
 
 lemma rapprox_rat: "real x / real y \<le> real (rapprox_rat prec x y)"
   using round_up by (simp add: rapprox_rat_def)
@@ -1130,32 +1746,17 @@
 proof -
   have "bitlen \<bar>x\<bar> \<le> bitlen \<bar>y\<bar>"
     using xy unfolding bitlen_def by (auto intro!: floor_mono)
-  then have "0 \<le> rat_precision n \<bar>x\<bar> \<bar>y\<bar>" by (simp add: rat_precision_def)
-  have "real \<lceil>real x / real y * 2 powr real (rat_precision n \<bar>x\<bar> \<bar>y\<bar>)\<rceil>
-      \<le> real \<lceil>2 powr real (rat_precision n \<bar>x\<bar> \<bar>y\<bar>)\<rceil>"
-    using xy by (auto intro!: ceiling_mono simp: field_simps)
-  also have "\<dots> = 2 powr real (rat_precision n \<bar>x\<bar> \<bar>y\<bar>)"
-    using `0 \<le> rat_precision n \<bar>x\<bar> \<bar>y\<bar>`
-    by (auto intro!: exI[of _ "2^nat (rat_precision n \<bar>x\<bar> \<bar>y\<bar>)"] simp: powr_int)
-  finally show ?thesis
-    by (simp add: rapprox_rat_def round_up_def)
-       (simp add: powr_minus inverse_eq_divide)
+  from this assms show ?thesis
+    by transfer (auto intro!: round_up_le1 simp: rat_precision_def)
 qed
 
-lemma rapprox_rat_nonneg_neg:
-  "0 \<le> x \<Longrightarrow> y < 0 \<Longrightarrow> real (rapprox_rat n x y) \<le> 0"
-  unfolding rapprox_rat_def round_up_def
-  by (auto simp: field_simps mult_le_0_iff zero_le_mult_iff)
+lemma rapprox_rat_nonneg_nonpos:
+  "0 \<le> x \<Longrightarrow> y \<le> 0 \<Longrightarrow> real (rapprox_rat n x y) \<le> 0"
+  by transfer (simp add: round_up_le0 divide_nonneg_nonpos)
 
-lemma rapprox_rat_neg:
-  "x < 0 \<Longrightarrow> 0 < y \<Longrightarrow> real (rapprox_rat n x y) \<le> 0"
-  unfolding rapprox_rat_def round_up_def
-  by (auto simp: field_simps mult_le_0_iff)
-
-lemma rapprox_rat_nonpos_pos:
-  "x \<le> 0 \<Longrightarrow> 0 < y \<Longrightarrow> real (rapprox_rat n x y) \<le> 0"
-  unfolding rapprox_rat_def round_up_def
-  by (auto simp: field_simps mult_le_0_iff)
+lemma rapprox_rat_nonpos_nonneg:
+  "x \<le> 0 \<Longrightarrow> 0 \<le> y \<Longrightarrow> real (rapprox_rat n x y) \<le> 0"
+  by transfer (simp add: round_up_le0 divide_nonpos_nonneg)
 
 lemma real_divl: "real_divl prec x y \<le> x / y"
   by (simp add: real_divl_def round_down)
@@ -1168,7 +1769,7 @@
 
 lemma real_divl_lower_bound:
   "0 \<le> x \<Longrightarrow> 0 \<le> y \<Longrightarrow> 0 \<le> real_divl prec x y"
-  by (simp add: real_divl_def round_down_def zero_le_mult_iff zero_le_divide_iff)
+  by (simp add: real_divl_def round_down_nonneg)
 
 lemma float_divl_lower_bound:
   "0 \<le> x \<Longrightarrow> 0 \<le> y \<Longrightarrow> 0 \<le> real (float_divl prec x y)"
@@ -1196,149 +1797,52 @@
   also have "mantissa x \<le> \<bar>mantissa x\<bar>" by simp
   also have "... \<le> 2 powr (bitlen \<bar>mantissa x\<bar>)"
     using bitlen_bounds[of "\<bar>mantissa x\<bar>"] bitlen_nonneg `mantissa x \<noteq> 0`
-    by (simp add: powr_int) (simp only: two_real_int int_of_reals real_of_int_abs[symmetric]
-      real_of_int_le_iff less_imp_le)
+    by (auto simp del: real_of_int_abs simp add: powr_int)
   finally show ?thesis by (simp add: powr_add)
 qed
 
 lemma real_divl_pos_less1_bound:
-  "0 < x \<Longrightarrow> x < 1 \<Longrightarrow> prec \<ge> 1 \<Longrightarrow> 1 \<le> real_divl prec 1 x"
-proof (unfold real_divl_def)
-  fix prec :: nat and x :: real assume x: "0 < x" "x < 1" and prec: "1 \<le> prec"
-  def p \<equiv> "int prec + \<lfloor>log 2 \<bar>x\<bar>\<rfloor>"
-  show "1 \<le> round_down (int prec + \<lfloor>log 2 \<bar>x\<bar>\<rfloor> - \<lfloor>log 2 \<bar>1\<bar>\<rfloor>) (1 / x) "
-  proof cases
-    assume nonneg: "0 \<le> p"
-    hence "2 powr real (p) = floor (real ((2::int) ^ nat p)) * floor (1::real)"
-      by (simp add: powr_int del: real_of_int_power) simp
-    also have "floor (1::real) \<le> floor (1 / x)" using x prec by simp
-    also have "floor (real ((2::int) ^ nat p)) * floor (1 / x) \<le>
-      floor (real ((2::int) ^ nat p) * (1 / x))"
-      by (rule le_mult_floor) (auto simp: x prec less_imp_le)
-    finally have "2 powr real p \<le> floor (2 powr nat p / x)" by (simp add: powr_realpow)
-    thus ?thesis unfolding p_def[symmetric]
-      using x prec nonneg by (simp add: powr_minus inverse_eq_divide round_down_def)
-  next
-    assume neg: "\<not> 0 \<le> p"
-
-    have "x = 2 powr (log 2 x)"
-      using x by simp
-    also have "2 powr (log 2 x) \<le> 2 powr p"
-    proof (rule powr_mono)
-      have "log 2 x \<le> \<lceil>log 2 x\<rceil>"
-        by simp
-      also have "\<dots> \<le> \<lfloor>log 2 x\<rfloor> + 1"
-        using ceiling_diff_floor_le_1[of "log 2 x"] by simp
-      also have "\<dots> \<le> \<lfloor>log 2 x\<rfloor> + prec"
-        using prec by simp
-      finally show "log 2 x \<le> real p"
-        using x by (simp add: p_def)
-    qed simp
-    finally have x_le: "x \<le> 2 powr p" .
-
-    from neg have "2 powr real p \<le> 2 powr 0"
-      by (intro powr_mono) auto
-    also have "\<dots> \<le> \<lfloor>2 powr 0\<rfloor>" by simp
-    also have "\<dots> \<le> \<lfloor>2 powr real p / x\<rfloor>" unfolding real_of_int_le_iff
-      using x x_le by (intro floor_mono) (simp add:  pos_le_divide_eq)
-    finally show ?thesis
-      using prec x unfolding p_def[symmetric]
-      by (simp add: round_down_def powr_minus_divide pos_le_divide_eq)
-  qed
+  assumes "0 < x" "x \<le> 1" "prec \<ge> 1"
+  shows "1 \<le> real_divl prec 1 x"
+proof -
+  have "log 2 x \<le> real prec + real \<lfloor>log 2 x\<rfloor>" using `prec \<ge> 1` by arith
+  from this assms show ?thesis
+    by (simp add: real_divl_def log_divide round_down_ge1)
 qed
 
 lemma float_divl_pos_less1_bound:
-  "0 < real x \<Longrightarrow> real x < 1 \<Longrightarrow> prec \<ge> 1 \<Longrightarrow> 1 \<le> real (float_divl prec 1 x)"
+  "0 < real x \<Longrightarrow> real x \<le> 1 \<Longrightarrow> prec \<ge> 1 \<Longrightarrow> 1 \<le> real (float_divl prec 1 x)"
   by (transfer, rule real_divl_pos_less1_bound)
 
 lemma float_divr: "real x / real y \<le> real (float_divr prec x y)"
   by transfer (rule real_divr)
 
-lemma real_divr_pos_less1_lower_bound: assumes "0 < x" and "x < 1" shows "1 \<le> real_divr prec 1 x"
+lemma real_divr_pos_less1_lower_bound: assumes "0 < x" and "x \<le> 1" shows "1 \<le> real_divr prec 1 x"
 proof -
-  have "1 \<le> 1 / x" using `0 < x` and `x < 1` by auto
+  have "1 \<le> 1 / x" using `0 < x` and `x <= 1` by auto
   also have "\<dots> \<le> real_divr prec 1 x" using real_divr[where x=1 and y=x] by auto
   finally show ?thesis by auto
 qed
 
-lemma float_divr_pos_less1_lower_bound: "0 < x \<Longrightarrow> x < 1 \<Longrightarrow> 1 \<le> float_divr prec 1 x"
+lemma float_divr_pos_less1_lower_bound: "0 < x \<Longrightarrow> x \<le> 1 \<Longrightarrow> 1 \<le> float_divr prec 1 x"
   by transfer (rule real_divr_pos_less1_lower_bound)
 
 lemma real_divr_nonpos_pos_upper_bound:
-  "x \<le> 0 \<Longrightarrow> 0 < y \<Longrightarrow> real_divr prec x y \<le> 0"
-  by (auto simp: field_simps mult_le_0_iff divide_le_0_iff round_up_def real_divr_def)
+  "x \<le> 0 \<Longrightarrow> 0 \<le> y \<Longrightarrow> real_divr prec x y \<le> 0"
+  by (simp add: real_divr_def round_up_le0 divide_le_0_iff)
 
 lemma float_divr_nonpos_pos_upper_bound:
-  "real x \<le> 0 \<Longrightarrow> 0 < real y \<Longrightarrow> real (float_divr prec x y) \<le> 0"
+  "real x \<le> 0 \<Longrightarrow> 0 \<le> real y \<Longrightarrow> real (float_divr prec x y) \<le> 0"
   by transfer (rule real_divr_nonpos_pos_upper_bound)
 
 lemma real_divr_nonneg_neg_upper_bound:
-  "0 \<le> x \<Longrightarrow> y < 0 \<Longrightarrow> real_divr prec x y \<le> 0"
-  by (auto simp: field_simps mult_le_0_iff zero_le_mult_iff divide_le_0_iff round_up_def real_divr_def)
+  "0 \<le> x \<Longrightarrow> y \<le> 0 \<Longrightarrow> real_divr prec x y \<le> 0"
+  by (simp add: real_divr_def round_up_le0 divide_le_0_iff)
 
 lemma float_divr_nonneg_neg_upper_bound:
-  "0 \<le> real x \<Longrightarrow> real y < 0 \<Longrightarrow> real (float_divr prec x y) \<le> 0"
+  "0 \<le> real x \<Longrightarrow> real y \<le> 0 \<Longrightarrow> real (float_divr prec x y) \<le> 0"
   by transfer (rule real_divr_nonneg_neg_upper_bound)
 
-definition truncate_down::"nat \<Rightarrow> real \<Rightarrow> real" where
-  "truncate_down prec x = round_down (prec - \<lfloor>log 2 \<bar>x\<bar>\<rfloor> - 1) x"
-
-lemma truncate_down: "truncate_down prec x \<le> x"
-  using round_down by (simp add: truncate_down_def)
-
-lemma truncate_down_le: "x \<le> y \<Longrightarrow> truncate_down prec x \<le> y"
-  by (rule order_trans[OF truncate_down])
-
-definition truncate_up::"nat \<Rightarrow> real \<Rightarrow> real" where
-  "truncate_up prec x = round_up (prec - \<lfloor>log 2 \<bar>x\<bar>\<rfloor> - 1) x"
-
-lemma truncate_up: "x \<le> truncate_up prec x"
-  using round_up by (simp add: truncate_up_def)
-
-lemma truncate_up_le: "x \<le> y \<Longrightarrow> x \<le> truncate_up prec y"
-  by (rule order_trans[OF _ truncate_up])
-
-lemma truncate_up_zero[simp]: "truncate_up prec 0 = 0"
-  by (simp add: truncate_up_def)
-
-lift_definition float_round_up :: "nat \<Rightarrow> float \<Rightarrow> float" is truncate_up
-  by (simp add: truncate_up_def)
-
-lemma float_round_up: "real x \<le> real (float_round_up prec x)"
-  using truncate_up by transfer simp
-
-lift_definition float_round_down :: "nat \<Rightarrow> float \<Rightarrow> float" is truncate_down
-  by (simp add: truncate_down_def)
-
-lemma float_round_down: "real (float_round_down prec x) \<le> real x"
-  using truncate_down by transfer simp
-
-lemma floor_add2[simp]: "\<lfloor> real i + x \<rfloor> = i + \<lfloor> x \<rfloor>"
-  using floor_add[of x i] by (simp del: floor_add add: ac_simps)
-
-lemma compute_float_round_down[code]:
-  "float_round_down prec (Float m e) = (let d = bitlen (abs m) - int prec in
-    if 0 < d then let P = 2^nat d ; n = m div P in Float n (e + d)
-             else Float m e)"
-  using Float.compute_float_down[of "prec - bitlen \<bar>m\<bar> - e" m e, symmetric]
-  by transfer (simp add: field_simps abs_mult log_mult bitlen_def truncate_down_def
-    cong del: if_weak_cong)
-hide_fact (open) compute_float_round_down
-
-lemma compute_float_round_up[code]:
-  "float_round_up prec (Float m e) = (let d = (bitlen (abs m) - int prec) in
-     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)
-              else Float m e)"
-  using Float.compute_float_up[of "prec - bitlen \<bar>m\<bar> - e" m e, symmetric]
-  unfolding Let_def
-  by transfer (simp add: field_simps abs_mult log_mult bitlen_def truncate_up_def
-    cong del: if_weak_cong)
-hide_fact (open) compute_float_round_up
-
-lemma round_up_mono: "x \<le> y \<Longrightarrow> round_up p x \<le> round_up p y"
-  by (auto intro!: ceiling_mono simp: round_up_def)
-
 lemma truncate_up_nonneg_mono:
   assumes "0 \<le> x" "x \<le> y"
   shows "truncate_up prec x \<le> truncate_up prec y"
@@ -1394,13 +1898,6 @@
     by blast
 qed
 
-lemma truncate_up_nonpos: "x \<le> 0 \<Longrightarrow> truncate_up prec x \<le> 0"
-  by (auto simp: truncate_up_def round_up_def intro!: mult_nonpos_nonneg)
-
-lemma truncate_down_nonpos: "x \<le> 0 \<Longrightarrow> truncate_down prec x \<le> 0"
-  by (auto simp: truncate_down_def round_down_def intro!: mult_nonpos_nonneg
-    order_le_less_trans[of _ 0, OF mult_nonpos_nonneg])
-
 lemma truncate_up_switch_sign_mono:
   assumes "x \<le> 0" "0 \<le> y"
   shows "truncate_up prec x \<le> truncate_up prec y"
@@ -1437,32 +1934,16 @@
     by (auto simp: truncate_down_def round_down_def)
 qed
 
-lemma truncate_down_nonneg: "0 \<le> y \<Longrightarrow> 0 \<le> truncate_down prec y"
-  by (auto simp: truncate_down_def round_down_def)
-
-lemma truncate_down_zero: "truncate_down prec 0 = 0"
-  by (auto simp: truncate_down_def round_down_def)
-
 lemma truncate_down_switch_sign_mono:
   assumes "x \<le> 0" "0 \<le> y"
   assumes "x \<le> y"
   shows "truncate_down prec x \<le> truncate_down prec y"
 proof -
-  note truncate_down_nonpos[OF `x \<le> 0`]
+  note truncate_down_le[OF `x \<le> 0`]
   also note truncate_down_nonneg[OF `0 \<le> y`]
   finally show ?thesis .
 qed
 
-lemma truncate_up_uminus_truncate_down:
-  "truncate_up prec x = - truncate_down prec (- x)"
-  "truncate_up prec (-x) = - truncate_down prec x"
-  by (auto simp: truncate_up_def round_up_def truncate_down_def round_down_def ceiling_def)
-
-lemma truncate_down_uminus_truncate_up:
-  "truncate_down prec x = - truncate_up prec (- x)"
-  "truncate_down prec (-x) = - truncate_up prec x"
-  by (auto simp: truncate_up_def round_up_def truncate_down_def round_down_def ceiling_def)
-
 lemma truncate_down_nonneg_mono:
   assumes "0 \<le> x" "x \<le> y"
   shows "truncate_down prec x \<le> truncate_down prec y"
@@ -1475,7 +1956,7 @@
     assume "~ 0 < x"
     with assms have "x = 0" "0 \<le> y" by simp_all
     hence ?thesis
-      by (auto simp add: truncate_down_zero intro!: truncate_down_nonneg)
+      by (auto intro!: truncate_down_nonneg)
   } moreover {
     assume "\<lfloor>log 2 \<bar>x\<bar>\<rfloor> = \<lfloor>log 2 \<bar>y\<bar>\<rfloor>"
     hence ?thesis
@@ -1522,16 +2003,20 @@
   } ultimately show ?thesis by blast
 qed
 
+lemma truncate_down_eq_truncate_up: "truncate_down p x = - truncate_up p (-x)"
+  and truncate_up_eq_truncate_down: "truncate_up p x = - truncate_down p (-x)"
+  by (auto simp: truncate_up_uminus_eq truncate_down_uminus_eq)
+
 lemma truncate_down_mono: "x \<le> y \<Longrightarrow> truncate_down p x \<le> truncate_down p y"
   apply (cases "0 \<le> x")
   apply (rule truncate_down_nonneg_mono, assumption+)
-  apply (simp add: truncate_down_uminus_truncate_up)
+  apply (simp add: truncate_down_eq_truncate_up)
   apply (cases "0 \<le> y")
   apply (auto intro: truncate_up_nonneg_mono truncate_up_switch_sign_mono)
   done
 
 lemma truncate_up_mono: "x \<le> y \<Longrightarrow> truncate_up p x \<le> truncate_up p y"
-  by (simp add: truncate_up_uminus_truncate_down truncate_down_mono)
+  by (simp add: truncate_up_eq_truncate_down truncate_down_mono)
 
 lemma Float_le_zero_iff: "Float a b \<le> 0 \<longleftrightarrow> a \<le> 0"
  apply (auto simp: zero_float_def mult_le_0_iff)
@@ -1573,5 +2058,13 @@
   then show ?thesis by simp
 qed
 
+lemma compute_mantissa[code]:
+  "mantissa (Float m e) = (if m = 0 then 0 else if 2 dvd m then mantissa (normfloat (Float m e)) else m)"
+  by (auto simp: mantissa_float Float.abs_eq)
+
+lemma compute_exponent[code]:
+  "exponent (Float m e) = (if m = 0 then 0 else if 2 dvd m then exponent (normfloat (Float m e)) else e)"
+  by (auto simp: exponent_float Float.abs_eq)
+
 end
 
diff -r 302104d8366b -r 87d4ce309e04 src/HOL/Real.thy
--- a/src/HOL/Real.thy	Wed Nov 12 18:18:38 2014 +0100
+++ b/src/HOL/Real.thy	Wed Nov 12 19:30:56 2014 +0100
@@ -1550,6 +1550,25 @@
 by (auto simp add: power2_eq_square)
 
 
+lemma numeral_power_eq_real_of_int_cancel_iff[simp]:
+  "numeral x ^ n = real (y::int) \<longleftrightarrow> numeral x ^ n = y"
+  by (metis real_numeral(1) real_of_int_inject real_of_int_power)
+
+lemma real_of_int_eq_numeral_power_cancel_iff[simp]:
+  "real (y::int) = numeral x ^ n \<longleftrightarrow> y = numeral x ^ n"
+  using numeral_power_eq_real_of_int_cancel_iff[of x n y]
+  by metis
+
+lemma numeral_power_eq_real_of_nat_cancel_iff[simp]:
+  "numeral x ^ n = real (y::nat) \<longleftrightarrow> numeral x ^ n = y"
+  by (metis of_nat_eq_iff of_nat_numeral real_of_int_eq_numeral_power_cancel_iff
+    real_of_int_of_nat_eq zpower_int)
+
+lemma real_of_nat_eq_numeral_power_cancel_iff[simp]:
+  "real (y::nat) = numeral x ^ n \<longleftrightarrow> y = numeral x ^ n"
+  using numeral_power_eq_real_of_nat_cancel_iff[of x n y]
+  by metis
+
 lemma numeral_power_le_real_of_nat_cancel_iff[simp]:
   "(numeral x::real) ^ n \<le> real a \<longleftrightarrow> (numeral x::nat) ^ n \<le> a"
   unfolding real_of_nat_le_iff[symmetric] by simp
@@ -1566,6 +1585,22 @@
   "real a \<le> (numeral x::real) ^ n \<longleftrightarrow> a \<le> (numeral x::int) ^ n"
   unfolding real_of_int_le_iff[symmetric] by simp
 
+lemma numeral_power_less_real_of_nat_cancel_iff[simp]:
+  "(numeral x::real) ^ n < real a \<longleftrightarrow> (numeral x::nat) ^ n < a"
+  unfolding real_of_nat_less_iff[symmetric] by simp
+
+lemma real_of_nat_less_numeral_power_cancel_iff[simp]:
+  "real a < (numeral x::real) ^ n \<longleftrightarrow> a < (numeral x::nat) ^ n"
+  unfolding real_of_nat_less_iff[symmetric] by simp
+
+lemma numeral_power_less_real_of_int_cancel_iff[simp]:
+  "(numeral x::real) ^ n < real a \<longleftrightarrow> (numeral x::int) ^ n < a"
+  unfolding real_of_int_less_iff[symmetric] by simp
+
+lemma real_of_int_less_numeral_power_cancel_iff[simp]:
+  "real a < (numeral x::real) ^ n \<longleftrightarrow> a < (numeral x::int) ^ n"
+  unfolding real_of_int_less_iff[symmetric] by simp
+
 lemma neg_numeral_power_le_real_of_int_cancel_iff[simp]:
   "(- numeral x::real) ^ n \<le> real a \<longleftrightarrow> (- numeral x::int) ^ n \<le> a"
   unfolding real_of_int_le_iff[symmetric] by simp
@@ -1719,9 +1754,15 @@
 lemma floor_le_eq: "(floor x <= a) = (x < real a + 1)"
   by linarith
 
+lemma floor_eq_iff: "floor x = b \<longleftrightarrow> real b \<le> x \<and> x < real (b + 1)"
+  by linarith
+
 lemma floor_add [simp]: "floor (x + real a) = floor x + a"
   by linarith
 
+lemma floor_add2[simp]: "floor (real a + x) = a + floor x"
+  by linarith
+
 lemma floor_subtract [simp]: "floor (x - real a) = floor x - a"
   by linarith
 
@@ -2015,6 +2056,14 @@
     by simp
 qed
 
+lemma floor_numeral_power[simp]:
+  "\<lfloor>numeral x ^ n\<rfloor> = numeral x ^ n"
+  by (metis floor_of_int of_int_numeral of_int_power)
+
+lemma ceiling_numeral_power[simp]:
+  "\<lceil>numeral x ^ n\<rceil> = numeral x ^ n"
+  by (metis ceiling_of_int of_int_numeral of_int_power)
+
 
 subsection {* Implementation of rational real numbers *}
 
diff -r 302104d8366b -r 87d4ce309e04 src/HOL/Transcendental.thy
--- a/src/HOL/Transcendental.thy	Wed Nov 12 18:18:38 2014 +0100
+++ b/src/HOL/Transcendental.thy	Wed Nov 12 19:30:56 2014 +0100
@@ -1861,6 +1861,9 @@
 lemma powr_minus_divide: "x powr (-a) = 1/(x powr a)"
   by (simp add: divide_inverse powr_minus)
 
+lemma divide_powr_uminus: "a / b powr c = a * b powr (- c)"
+  by (simp add: powr_minus_divide)
+
 lemma powr_less_mono: "a < b \<Longrightarrow> 1 < x \<Longrightarrow> x powr a < x powr b"
   by (simp add: powr_def)
 
@@ -1926,6 +1929,12 @@
 lemma log_divide: "0 < a \<Longrightarrow> a \<noteq> 1 \<Longrightarrow> 0 < x \<Longrightarrow> 0 < y \<Longrightarrow> log a (x/y) = log a x - log a y"
   by (simp add: log_mult divide_inverse log_inverse)
 
+lemma log_add_eq_powr: "0 < b \<Longrightarrow> b \<noteq> 1 \<Longrightarrow> 0 < x \<Longrightarrow> log b x + y = log b (x * b powr y)"
+  and add_log_eq_powr: "0 < b \<Longrightarrow> b \<noteq> 1 \<Longrightarrow> 0 < x \<Longrightarrow> y + log b x = log b (b powr y * x)"
+  and log_minus_eq_powr: "0 < b \<Longrightarrow> b \<noteq> 1 \<Longrightarrow> 0 < x \<Longrightarrow> log b x - y = log b (x * b powr -y)"
+  and minus_log_eq_powr: "0 < b \<Longrightarrow> b \<noteq> 1 \<Longrightarrow> 0 < x \<Longrightarrow> y - log b x = log b (b powr y / x)"
+  by (simp_all add: log_mult log_divide)
+
 lemma log_less_cancel_iff [simp]:
   "1 < a \<Longrightarrow> 0 < x \<Longrightarrow> 0 < y \<Longrightarrow> log a x < log a y \<longleftrightarrow> x < y"
   apply safe
@@ -1982,6 +1991,34 @@
 lemma log_le_one_cancel_iff[simp]: "1 < a \<Longrightarrow> 0 < x \<Longrightarrow> log a x \<le> 1 \<longleftrightarrow> x \<le> a"
   using log_le_cancel_iff[of a x a] by simp
 
+lemma le_log_iff:
+  assumes "1 < b" "x > 0"
+  shows "y \<le> log b x \<longleftrightarrow> b powr y \<le> x"
+  by (metis assms(1) assms(2) dual_order.strict_trans powr_le_cancel_iff powr_log_cancel
+    powr_one_eq_one powr_one_gt_zero_iff)
+
+lemma less_log_iff:
+  assumes "1 < b" "x > 0"
+  shows "y < log b x \<longleftrightarrow> b powr y < x"
+  by (metis assms(1) assms(2) dual_order.strict_trans less_irrefl powr_less_cancel_iff
+    powr_log_cancel zero_less_one)
+
+lemma
+  assumes "1 < b" "x > 0"
+  shows log_less_iff: "log b x < y \<longleftrightarrow> x < b powr y"
+    and log_le_iff: "log b x \<le> y \<longleftrightarrow> x \<le> b powr y"
+  using le_log_iff[OF assms, of y] less_log_iff[OF assms, of y]
+  by auto
+
+lemmas powr_le_iff = le_log_iff[symmetric]
+  and powr_less_iff = le_log_iff[symmetric]
+  and less_powr_iff = log_less_iff[symmetric]
+  and le_powr_iff = log_le_iff[symmetric]
+
+lemma
+  floor_log_eq_powr_iff: "x > 0 \<Longrightarrow> b > 1 \<Longrightarrow> \<lfloor>log b x\<rfloor> = k \<longleftrightarrow> b powr k \<le> x \<and> x < b powr (k + 1)"
+  by (auto simp add: floor_eq_iff powr_le_iff less_powr_iff)
+
 lemma powr_realpow: "0 < x ==> x powr (real n) = x^n"
   apply (induct n)
   apply simp
@@ -2013,6 +2050,14 @@
   then show ?thesis by (simp add: assms powr_realpow[symmetric])
 qed
 
+lemma compute_powr[code]:
+  fixes i::real
+  shows "b powr i =
+    (if b \<le> 0 then Code.abort (STR ''op powr with nonpositive base'') (\<lambda>_. b powr i)
+    else if floor i = i then (if 0 \<le> i then b ^ natfloor i else 1 / b ^ natfloor (- i))
+    else Code.abort (STR ''op powr with non-integer exponent'') (\<lambda>_. b powr i))"
+  by (auto simp: natfloor_def powr_int)
+
 lemma powr_one: "0 < x \<Longrightarrow> x powr 1 = x"
   using powr_realpow [of x 1] by simp