--- 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
--- 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