# HG changeset patch # User wenzelm # Date 1436222922 -7200 # Node ID 589ed01b94fe333a968e41e80a50304cb19fdc3f # Parent ade12ef2773c53605eec77a21b58bb34afd6c961 tuned proofs; diff -r ade12ef2773c -r 589ed01b94fe src/HOL/Decision_Procs/Approximation.thy --- a/src/HOL/Decision_Procs/Approximation.thy Mon Jul 06 22:57:34 2015 +0200 +++ b/src/HOL/Decision_Procs/Approximation.thy Tue Jul 07 00:48:42 2015 +0200 @@ -121,7 +121,7 @@ shows "(lb n ((F ^^ j') s) (f j') x) \ (\j=0..j=0.. (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 diff_mult_minus: "x - y * z = x + - y * z" for x y z :: float by simp have sum_eq: "(\j=0..j = 0..Selectors for next even or odd number\ text \ - The horner scheme computes alternating series. To get the upper and lower bounds we need to guarantee to access a even or odd member. To do this we use @{term get_odd} and @{term get_even}. - \ definition get_odd :: "nat \ nat" where @@ -149,16 +148,21 @@ definition get_even :: "nat \ nat" where "get_even n = (if even n then n else (Suc n))" -lemma get_odd[simp]: "odd (get_odd n)" unfolding get_odd_def by (cases "odd n", auto) -lemma get_even[simp]: "even (get_even n)" unfolding get_even_def by (cases "even n", auto) +lemma get_odd[simp]: "odd (get_odd n)" + unfolding get_odd_def by (cases "odd n") auto + +lemma get_even[simp]: "even (get_even n)" + unfolding get_even_def by (cases "even n") auto + lemma get_odd_ex: "\ k. Suc k = get_odd n \ odd (Suc k)" by (auto simp: get_odd_def odd_pos intro!: exI[of _ "n - 1"]) -lemma get_even_double: - "\i. get_even n = 2 * i" using get_even by (blast elim: evenE) - -lemma get_odd_double: - "\i. get_odd n = 2 * i + 1" using get_odd by (blast elim: oddE) +lemma get_even_double: "\i. get_even n = 2 * i" + using get_even by (blast elim: evenE) + +lemma get_odd_double: "\i. get_odd n = 2 * i + 1" + using get_odd by (blast elim: oddE) + section "Power function" @@ -183,17 +187,16 @@ intro: power_mono_odd power_mono power_mono_even zero_le_even_power) lemma bnds_power: - "\ (x::real) l u. (l1, u1) = float_power_bnds prec n l u \ x \ {l .. u} \ + "\(x::real) l u. (l1, u1) = float_power_bnds prec n l u \ x \ {l .. u} \ l1 \ x ^ n \ x ^ n \ u1" using float_power_bnds by auto + section "Square root" text \ - The square root computation is implemented as newton iteration. As first first step we use the nearest power of two greater than the square root. - \ fun sqrt_iteration :: "nat \ nat \ float \ float" where @@ -234,34 +237,40 @@ thus ?thesis by (simp add: field_simps) qed -lemma sqrt_iteration_bound: assumes "0 < real x" +lemma sqrt_iteration_bound: + assumes "0 < real x" shows "sqrt x < sqrt_iteration prec n x" proof (induct n) case 0 show ?case proof (cases x) case (Float m e) - hence "0 < m" using assms + hence "0 < m" + using assms apply (auto simp: sign_simps) by (meson not_less powr_ge_pzero) hence "0 < sqrt m" by auto - have int_nat_bl: "(nat (bitlen m)) = bitlen m" using bitlen_nonneg by auto + have int_nat_bl: "(nat (bitlen m)) = bitlen m" + using bitlen_nonneg by auto have "x = (m / 2^nat (bitlen m)) * 2 powr (e + (nat (bitlen m)))" unfolding Float by (auto simp: powr_realpow[symmetric] field_simps powr_add) also have "\ < 1 * 2 powr (e + nat (bitlen m))" proof (rule mult_strict_right_mono, auto) - show "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 + finally have "sqrt x < sqrt (2 powr (e + bitlen m))" + unfolding int_nat_bl by auto also have "\ \ 2 powr ((e + bitlen m) div 2 + 1)" proof - let ?E = "e + bitlen m" have E_mod_pow: "2 powr (?E mod 2) < 4" proof (cases "?E mod 2 = 1") - case True thus ?thesis by auto + case True + thus ?thesis by auto next case False have "0 \ ?E mod 2" by auto @@ -275,15 +284,16 @@ by (auto simp del: real_sqrt_four) hence E_mod_pow: "sqrt (2 powr (?E mod 2)) < 2" by auto - have E_eq: "2 powr ?E = 2 powr (?E div 2 + ?E div 2 + ?E mod 2)" by auto + have E_eq: "2 powr ?E = 2 powr (?E div 2 + ?E div 2 + ?E mod 2)" + by auto have "sqrt (2 powr ?E) = sqrt (2 powr (?E div 2) * 2 powr (?E div 2) * 2 powr (?E mod 2))" unfolding E_eq unfolding powr_add[symmetric] by (simp add: int_of_reals del: real_of_ints) also have "\ = 2 powr (?E div 2) * sqrt (2 powr (?E mod 2))" unfolding real_sqrt_mult[of _ "2 powr (?E mod 2)"] real_sqrt_abs2 by auto also have "\ < 2 powr (?E div 2) * 2 powr 1" - by (rule mult_strict_left_mono, auto intro: E_mod_pow) - also have "\ = 2 powr (?E div 2 + 1)" unfolding add.commute[of _ 1] powr_add[symmetric] - by simp + by (rule mult_strict_left_mono) (auto intro: E_mod_pow) + also have "\ = 2 powr (?E div 2 + 1)" + unfolding add.commute[of _ 1] powr_add[symmetric] by simp finally show ?thesis by auto qed finally show ?thesis using \0 < m\ @@ -293,18 +303,24 @@ next case (Suc n) let ?b = "sqrt_iteration prec n x" - have "0 < sqrt x" using \0 < real x\ by auto - also have "\ < 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 + have "0 < sqrt x" + using \0 < real x\ by auto + also have "\ < 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 "\ \ (?b + (float_divr prec x ?b))/2" by (rule divide_right_mono, auto simp add: float_divr) - also have "\ = (Float 1 (- 1)) * (?b + (float_divr prec x ?b))" by simp + also have "\ = (Float 1 (- 1)) * (?b + (float_divr prec x ?b))" + by simp also have "\ \ (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 . + finally show ?case + unfolding sqrt_iteration.simps Let_def distrib_left . qed -lemma sqrt_iteration_lower_bound: assumes "0 < real x" +lemma sqrt_iteration_lower_bound: + assumes "0 < real x" shows "0 < real (sqrt_iteration prec n x)" (is "0 < ?sqrt") proof - have "0 < sqrt x" using assms by auto @@ -312,25 +328,34 @@ finally show ?thesis . qed -lemma lb_sqrt_lower_bound: assumes "0 \ real x" +lemma lb_sqrt_lower_bound: + assumes "0 \ real x" shows "0 \ real (lb_sqrt prec x)" proof (cases "0 < x") - case True hence "0 < real x" and "0 \ x" using \0 \ real x\ by auto - hence "0 < sqrt_iteration prec prec x" using sqrt_iteration_lower_bound by auto - hence "0 \ real (float_divl prec x (sqrt_iteration prec prec x))" using float_divl_lower_bound[OF \0 \ x\] unfolding less_eq_float_def by auto - thus ?thesis unfolding lb_sqrt.simps using True by auto + case True + hence "0 < real x" and "0 \ x" + using \0 \ real x\ by auto + hence "0 < sqrt_iteration prec prec x" + using sqrt_iteration_lower_bound by auto + hence "0 \ real (float_divl prec x (sqrt_iteration prec prec x))" + using float_divl_lower_bound[OF \0 \ x\] unfolding less_eq_float_def by auto + thus ?thesis + unfolding lb_sqrt.simps using True by auto next - case False with \0 \ real x\ have "real x = 0" by auto - thus ?thesis unfolding lb_sqrt.simps by auto + case False + with \0 \ real x\ have "real x = 0" by auto + thus ?thesis + unfolding lb_sqrt.simps by auto qed lemma bnds_sqrt': "sqrt x \ {(lb_sqrt prec x) .. (ub_sqrt prec x)}" proof - - { fix x :: float assume "0 < x" - hence "0 < real x" and "0 \ real x" by auto + have lb: "lb_sqrt prec x \ sqrt x" if "0 < x" for x :: float + proof - + from that have "0 < real x" and "0 \ real x" by auto hence sqrt_gt0: "0 < sqrt x" by auto - hence sqrt_ub: "sqrt x < sqrt_iteration prec prec x" using sqrt_iteration_bound by auto - + hence sqrt_ub: "sqrt x < sqrt_iteration prec prec x" + using sqrt_iteration_bound by auto have "(float_divl prec x (sqrt_iteration prec prec x)) \ x / (sqrt_iteration prec prec x)" by (rule float_divl) also have "\ < x / sqrt x" @@ -339,27 +364,28 @@ also have "\ = sqrt x" unfolding inverse_eq_iff_eq[of _ "sqrt x", symmetric] sqrt_divide_self_eq[OF \0 \ real x\, symmetric] by auto - finally have "lb_sqrt prec x \ sqrt x" - unfolding lb_sqrt.simps if_P[OF \0 < x\] by auto } - note lb = this - - { fix x :: float assume "0 < x" - hence "0 < real x" by auto + finally show ?thesis + unfolding lb_sqrt.simps if_P[OF \0 < x\] by auto + qed + have ub: "sqrt x \ ub_sqrt prec x" if "0 < x" for x :: float + proof - + from that have "0 < real x" by auto hence "0 < sqrt x" by auto hence "sqrt x < sqrt_iteration prec prec x" using sqrt_iteration_bound by auto - hence "sqrt x \ ub_sqrt prec x" - unfolding ub_sqrt.simps if_P[OF \0 < x\] by auto } - note ub = this - + then show ?thesis + unfolding ub_sqrt.simps if_P[OF \0 < x\] by auto + qed show ?thesis using lb[of "-x"] ub[of "-x"] lb[of x] ub[of x] by (auto simp add: lb_sqrt.simps ub_sqrt.simps real_sqrt_minus) qed -lemma bnds_sqrt: "\ (x::real) lx ux. (l, u) = (lb_sqrt prec lx, ub_sqrt prec ux) \ x \ {lx .. ux} \ l \ sqrt x \ sqrt x \ u" +lemma bnds_sqrt: "\(x::real) lx ux. + (l, u) = (lb_sqrt prec lx, ub_sqrt prec ux) \ x \ {lx .. ux} \ l \ sqrt x \ sqrt x \ u" proof ((rule allI) +, rule impI, erule conjE, rule conjI) - fix x :: real fix lx ux + fix x :: real + fix lx ux assume "(l, u) = (lb_sqrt prec lx, ub_sqrt prec ux)" and x: "x \ {lx .. ux}" hence l: "l = lb_sqrt prec lx " and u: "u = ub_sqrt prec ux" by auto @@ -373,15 +399,14 @@ show "sqrt x \ u" unfolding u using bnds_sqrt'[of ux prec] by auto qed + section "Arcus tangens and \" subsection "Compute arcus tangens series" text \ - As first step we implement the computation of the arcus tangens series. This is only valid in the range @{term "{-1 :: real .. 1}"}. This is used to compute \ and then the entire arcus tangens. - \ fun ub_arctan_horner :: "nat \ nat \ nat \ float \ float" @@ -394,7 +419,8 @@ (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 \ real y" "real y \ 1" and "even n" + assumes "0 \ real y" "real y \ 1" + and "even n" shows "arctan (sqrt y) \ {(sqrt y * lb_arctan_horner prec n 1 y) .. (sqrt y * ub_arctan_horner prec (Suc n) 1 y)}" proof - @@ -407,6 +433,9 @@ have "arctan (sqrt y) \ { ?S n .. ?S (Suc n) }" proof (cases "sqrt y = 0") + case True + then show ?thesis by simp + next case False hence "0 < sqrt y" using \0 \ sqrt y\ by auto hence prem: "0 < 1 / (0 * 2 + (1::nat)) * sqrt y ^ (0 * 2 + 1)" by auto @@ -415,7 +444,7 @@ 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 \\ sqrt y \ \ 1\] Suc_eq_plus1 atLeast0LessThan . - qed auto + qed note arctan_bounds = this[unfolded atLeastAtMost_iff] have F: "\n. 2 * Suc n + 1 = 2 * n + 1 + 2" by auto @@ -425,25 +454,32 @@ and ub="\n i k x. ub_arctan_horner prec n k x", OF \0 \ real y\ F lb_arctan_horner.simps ub_arctan_horner.simps] - { have "(sqrt y * lb_arctan_horner prec n 1 y) \ ?S n" + have "(sqrt y * lb_arctan_horner prec n 1 y) \ arctan (sqrt y)" + proof - + have "(sqrt y * lb_arctan_horner prec n 1 y) \ ?S n" using bounds(1) \0 \ 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 by (auto intro!: mult_left_mono) also have "\ \ arctan (sqrt y)" using arctan_bounds .. - finally have "(sqrt y * lb_arctan_horner prec n 1 y) \ arctan (sqrt y)" . } + finally show ?thesis . + qed moreover - { have "arctan (sqrt y) \ ?S (Suc n)" using arctan_bounds .. + have "arctan (sqrt y) \ (sqrt y * ub_arctan_horner prec (Suc n) 1 y)" + proof - + have "arctan (sqrt y) \ ?S (Suc n)" using arctan_bounds .. also have "\ \ (sqrt y * ub_arctan_horner prec (Suc n) 1 y)" using bounds(2)[of "Suc n"] \0 \ 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 by (auto intro!: mult_left_mono) - finally have "arctan (sqrt y) \ (sqrt y * ub_arctan_horner prec (Suc n) 1 y)" . } + finally show ?thesis . + qed ultimately show ?thesis by auto qed -lemma arctan_0_1_bounds: assumes "0 \ real y" "real y \ 1" +lemma arctan_0_1_bounds: + assumes "0 \ real y" "real y \ 1" shows "arctan (sqrt y) \ {(sqrt y * lb_arctan_horner prec (get_even n) 1 y) .. (sqrt y * ub_arctan_horner prec (get_odd n) 1 y)}" @@ -451,8 +487,8 @@ 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 \ x" @@ -476,12 +512,15 @@ lemma arctan_mult_le: assumes "0 \ x" "x \ y" "y * z \ arctan y" shows "x * z \ arctan x" -proof cases - assume "x \ 0" +proof (cases "x = 0") + case True + then show ?thesis by simp +next + case False with assms have "z \ arctan y / y" by (simp add: field_simps) also have "\ \ arctan x / x" using assms \x \ 0\ by (auto intro!: arctan_divide_mono) finally show ?thesis using assms \x \ 0\ by (simp add: field_simps) -qed simp +qed lemma arctan_le_mult: assumes "0 < x" "x \ y" "arctan x \ x * z" @@ -514,7 +553,10 @@ ultimately show ?thesis by simp qed -lemma mult_nonneg_le_one: fixes a::real assumes "0 \ a" "0 \ b" "a \ 1" "b \ 1" shows "a * b \ 1" +lemma mult_nonneg_le_one: + fixes a :: real + assumes "0 \ a" "0 \ b" "a \ 1" "b \ 1" + shows "a * b \ 1" proof - have "a * b \ 1 * 1" by (intro mult_mono assms) simp_all @@ -564,7 +606,10 @@ 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 \ k" and "0 < k" and "1 \ k" by auto + { + fix prec n :: nat + fix k :: int + assume "1 < k" hence "0 \ k" and "0 < k" and "1 \ 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 @@ -581,7 +626,10 @@ finally have "arctan (1 / k) \ ?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 \ k" and "0 < k" by auto + { + fix prec n :: nat + fix k :: int + assume "1 < k" hence "0 \ 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 @@ -614,6 +662,7 @@ ultimately show ?thesis by auto qed + subsection "Compute arcus tangens in the entire domain" function lb_arctan :: "nat \ float \ float" and ub_arctan :: "nat \ float \ float" where @@ -661,7 +710,8 @@ declare ub_arctan_horner.simps[simp del] declare lb_arctan_horner.simps[simp del] -lemma lb_arctan_bound': assumes "0 \ real x" +lemma lb_arctan_bound': + assumes "0 \ real x" shows "lb_arctan prec x \ arctan x" proof - have "\ x < 0" and "0 \ x" @@ -674,13 +724,15 @@ show ?thesis proof (cases "x \ Float 1 (- 1)") - case True hence "real x \ 1" by simp + case True + hence "real x \ 1" by simp from arctan_0_1_bounds_round[OF \0 \ real x\ \real x \ 1\] show ?thesis unfolding lb_arctan.simps Let_def if_not_P[OF \\ x < 0\] if_P[OF True] using \0 \ x\ by (auto intro!: float_round_down_le) next - case False hence "0 < real x" by auto + case False + hence "0 < real x" by auto let ?R = "1 + sqrt (1 + real x * real 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)" @@ -707,44 +759,61 @@ show ?thesis proof (cases "x \ Float 1 1") case True - - have "x \ sqrt (1 + x * x)" using real_sqrt_sum_squares_ge2[where x=1, unfolded numeral_2_eq_2] by auto + have "x \ sqrt (1 + x * x)" + using real_sqrt_sum_squares_ge2[where x=1, unfolded numeral_2_eq_2] by auto also note \\ \ (ub_sqrt prec ?sxx)\ - finally have "real x \ ?fR" by (auto simp: float_plus_up.rep_eq plus_up_def intro!: truncate_up_le) - moreover have "?DIV \ real x / ?fR" by (rule float_divl) - ultimately have "real ?DIV \ 1" unfolding divide_le_eq_1_pos[OF \0 < real ?fR\, symmetric] by auto - - have "0 \ real ?DIV" using float_divl_lower_bound[OF \0 \ x\] \0 < ?fR\ unfolding less_eq_float_def by auto + finally have "real x \ ?fR" + by (auto simp: float_plus_up.rep_eq plus_up_def intro!: truncate_up_le) + moreover have "?DIV \ real x / ?fR" + by (rule float_divl) + ultimately have "real ?DIV \ 1" + unfolding divide_le_eq_1_pos[OF \0 < real ?fR\, symmetric] by auto + + have "0 \ real ?DIV" + using float_divl_lower_bound[OF \0 \ x\] \0 < ?fR\ + unfolding less_eq_float_def by auto from arctan_0_1_bounds_round[OF \0 \ real (?DIV)\ \real (?DIV) \ 1\] - have "Float 1 1 * ?lb_horner ?DIV \ 2 * arctan ?DIV" by simp + have "Float 1 1 * ?lb_horner ?DIV \ 2 * arctan ?DIV" + by simp also have "\ \ 2 * arctan (x / ?R)" 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 \\ x < 0\] if_not_P[OF \\ x \ 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]]) + 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 \\ x < 0\] + if_not_P[OF \\ x \ 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 hence "1 \ real x" by auto let "?invx" = "float_divr prec 1 x" - have "0 \ arctan x" using arctan_monotone'[OF \0 \ real x\] using arctan_tan[of 0, unfolded tan_zero] by auto + have "0 \ arctan x" using arctan_monotone'[OF \0 \ real x\] + using arctan_tan[of 0, unfolded tan_zero] by auto show ?thesis proof (cases "1 < ?invx") case True - show ?thesis unfolding lb_arctan.simps Let_def if_not_P[OF \\ x < 0\] if_not_P[OF \\ x \ Float 1 (- 1)\] if_not_P[OF False] if_P[OF True] + show ?thesis + unfolding lb_arctan.simps Let_def if_not_P[OF \\ x < 0\] + if_not_P[OF \\ x \ Float 1 (- 1)\] if_not_P[OF False] if_P[OF True] using \0 \ arctan x\ by auto next case False hence "real ?invx \ 1" by auto - have "0 \ real ?invx" by (rule order_trans[OF _ float_divr], auto simp add: \0 \ real x\) - - have "1 / x \ 0" and "0 < 1 / x" using \0 < real x\ by auto - - have "arctan (1 / x) \ arctan ?invx" unfolding one_float.rep_eq[symmetric] by (rule arctan_monotone', rule float_divr) - also have "\ \ ?ub_horner ?invx" using arctan_0_1_bounds_round[OF \0 \ real ?invx\ \real ?invx \ 1\] + have "0 \ real ?invx" + by (rule order_trans[OF _ float_divr]) (auto simp add: \0 \ real x\) + + have "1 / x \ 0" and "0 < 1 / x" + using \0 < real x\ by auto + + have "arctan (1 / x) \ arctan ?invx" + unfolding one_float.rep_eq[symmetric] by (rule arctan_monotone', rule float_divr) + also have "\ \ ?ub_horner ?invx" + using arctan_0_1_bounds_round[OF \0 \ real ?invx\ \real ?invx \ 1\] by (auto intro!: float_round_up_le) also note float_round_up finally have "pi / 2 - float_round_up prec (?ub_horner ?invx) \ arctan x" @@ -754,35 +823,45 @@ have "lb_pi prec * Float 1 (- 1) \ pi / 2" 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 \\ x < 0\] if_not_P[OF \\ x \ Float 1 (- 1)\] if_not_P[OF \\ x \ Float 1 1\] if_not_P[OF False] + show ?thesis + unfolding lb_arctan.simps Let_def if_not_P[OF \\ x < 0\] + if_not_P[OF \\ x \ Float 1 (- 1)\] if_not_P[OF \\ x \ Float 1 1\] if_not_P[OF False] by (auto intro!: float_plus_down_le) qed qed qed qed -lemma ub_arctan_bound': assumes "0 \ real x" +lemma ub_arctan_bound': + assumes "0 \ real x" shows "arctan x \ ub_arctan prec x" proof - - have "\ x < 0" and "0 \ x" using \0 \ real x\ by auto - - 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)))" + have "\ x < 0" and "0 \ x" + using \0 \ real x\ by auto + + 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)))" + let "?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 \ Float 1 (- 1)") - case True hence "real x \ 1" by auto - show ?thesis unfolding ub_arctan.simps Let_def if_not_P[OF \\ x < 0\] if_P[OF True] + case True + hence "real x \ 1" by auto + show ?thesis + unfolding ub_arctan.simps Let_def if_not_P[OF \\ x < 0\] if_P[OF True] using arctan_0_1_bounds_round[OF \0 \ real x\ \real x \ 1\] by (auto intro!: float_round_up_le) next - case False hence "0 < real x" by auto + case False + hence "0 < real x" by auto let ?R = "1 + sqrt (1 + real x * real 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 \ 1 + real x * real x" using sum_power2_ge_zero[of 1 "real x", unfolded numeral_2_eq_2] by auto + have sqr_ge0: "0 \ 1 + real x * real x" + using sum_power2_ge_zero[of 1 "real x", unfolded numeral_2_eq_2] by auto hence "0 \ real (1 + x*x)" by auto hence divisor_gt0: "0 < ?R" by (auto intro: add_pos_nonneg) @@ -792,7 +871,8 @@ also have "\ \ 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 \ sqrt (1 + x*x)" . - hence "?fR \ ?R" by (auto simp: float_plus_down.rep_eq plus_down_def truncate_down_le) + hence "?fR \ ?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] @@ -811,23 +891,31 @@ show ?thesis proof (cases "?DIV > 1") case True - have "pi / 2 \ ub_pi prec * Float 1 (- 1)" unfolding Float_num times_divide_eq_right mult_1_left using pi_boundaries by auto + have "pi / 2 \ ub_pi prec * Float 1 (- 1)" + unfolding Float_num times_divide_eq_right mult_1_left using pi_boundaries by auto from order_less_le_trans[OF arctan_ubound this, THEN less_imp_le] - show ?thesis unfolding ub_arctan.simps Let_def if_not_P[OF \\ x < 0\] if_not_P[OF \\ x \ Float 1 (- 1)\] if_P[OF \x \ Float 1 1\] if_P[OF True] . + show ?thesis + unfolding ub_arctan.simps Let_def if_not_P[OF \\ x < 0\] + if_not_P[OF \\ x \ Float 1 (- 1)\] if_P[OF \x \ Float 1 1\] if_P[OF True] . next case False hence "real ?DIV \ 1" by auto - have "0 \ x / ?R" using \0 \ real x\ \0 < ?R\ unfolding zero_le_divide_iff by auto - hence "0 \ real ?DIV" using monotone by (rule order_trans) - - have "arctan x = 2 * arctan (x / ?R)" using arctan_half unfolding numeral_2_eq_2 power_Suc2 power_0 mult_1_left . + have "0 \ x / ?R" + using \0 \ real x\ \0 < ?R\ unfolding zero_le_divide_iff by auto + hence "0 \ real ?DIV" + using monotone by (rule order_trans) + + have "arctan x = 2 * arctan (x / ?R)" + using arctan_half unfolding numeral_2_eq_2 power_Suc2 power_0 mult_1_left . also have "\ \ 2 * arctan (?DIV)" using arctan_monotone'[OF monotone] by (auto intro!: mult_left_mono) also have "\ \ (Float 1 1 * ?ub_horner ?DIV)" unfolding Float_num using arctan_0_1_bounds_round[OF \0 \ real ?DIV\ \real ?DIV \ 1\] by (auto intro!: float_round_up_le) - finally show ?thesis unfolding ub_arctan.simps Let_def if_not_P[OF \\ x < 0\] if_not_P[OF \\ x \ Float 1 (- 1)\] if_P[OF \x \ Float 1 1\] if_not_P[OF False] . + finally show ?thesis + unfolding ub_arctan.simps Let_def if_not_P[OF \\ x < 0\] + if_not_P[OF \\ x \ Float 1 (- 1)\] if_P[OF \x \ Float 1 1\] if_not_P[OF False] . qed next case False @@ -837,60 +925,89 @@ hence "0 < x" by auto let "?invx" = "float_divl prec 1 x" - have "0 \ arctan x" using arctan_monotone'[OF \0 \ real x\] using arctan_tan[of 0, unfolded tan_zero] by auto - - have "real ?invx \ 1" unfolding less_float_def by (rule order_trans[OF float_divl], auto simp add: \1 \ real x\ divide_le_eq_1_pos[OF \0 < real x\]) - have "0 \ real ?invx" using \0 < x\ by (intro float_divl_lower_bound) auto - - have "1 / x \ 0" and "0 < 1 / x" using \0 < real x\ by auto - - have "(?lb_horner ?invx) \ arctan (?invx)" using arctan_0_1_bounds_round[OF \0 \ real ?invx\ \real ?invx \ 1\] + have "0 \ arctan x" + using arctan_monotone'[OF \0 \ real x\] and arctan_tan[of 0, unfolded tan_zero] by auto + + have "real ?invx \ 1" + unfolding less_float_def + by (rule order_trans[OF float_divl]) + (auto simp add: \1 \ real x\ divide_le_eq_1_pos[OF \0 < real x\]) + have "0 \ real ?invx" + using \0 < x\ by (intro float_divl_lower_bound) auto + + have "1 / x \ 0" and "0 < 1 / x" + using \0 < real x\ by auto + + have "(?lb_horner ?invx) \ arctan (?invx)" + using arctan_0_1_bounds_round[OF \0 \ real ?invx\ \real ?invx \ 1\] by (auto intro!: float_round_down_le) - also have "\ \ arctan (1 / x)" unfolding one_float.rep_eq[symmetric] by (rule arctan_monotone', rule float_divl) + also have "\ \ arctan (1 / x)" + unfolding one_float.rep_eq[symmetric] by (rule arctan_monotone') (rule float_divl) finally have "arctan x \ pi / 2 - (?lb_horner ?invx)" using \0 \ arctan x\ arctan_inverse[OF \1 / x \ 0\] unfolding real_sgn_pos[OF \0 < 1 / x\] le_diff_eq by auto moreover - have "pi / 2 \ ub_pi prec * Float 1 (- 1)" unfolding Float_num times_divide_eq_right mult_1_right using pi_boundaries by auto + have "pi / 2 \ 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 \\ x < 0\]if_not_P[OF \\ x \ Float 1 (- 1)\] if_not_P[OF False] + show ?thesis + unfolding ub_arctan.simps Let_def if_not_P[OF \\ x < 0\] + if_not_P[OF \\ x \ Float 1 (- 1)\] if_not_P[OF False] by (auto intro!: float_round_up_le float_plus_up_le) qed qed qed -lemma arctan_boundaries: - "arctan x \ {(lb_arctan prec x) .. (ub_arctan prec x)}" +lemma arctan_boundaries: "arctan x \ {(lb_arctan prec x) .. (ub_arctan prec x)}" proof (cases "0 \ x") - case True hence "0 \ real x" by auto - show ?thesis using ub_arctan_bound'[OF \0 \ real x\] lb_arctan_bound'[OF \0 \ real x\] unfolding atLeastAtMost_iff by auto + case True + hence "0 \ real x" by auto + show ?thesis + using ub_arctan_bound'[OF \0 \ real x\] lb_arctan_bound'[OF \0 \ real x\] + unfolding atLeastAtMost_iff by auto next + case False let ?mx = "-x" - case False hence "x < 0" and "0 \ real ?mx" by auto + from False have "x < 0" and "0 \ real ?mx" + by auto hence bounds: "lb_arctan prec ?mx \ arctan ?mx \ arctan ?mx \ ub_arctan prec ?mx" using ub_arctan_bound'[OF \0 \ real ?mx\] lb_arctan_bound'[OF \0 \ real ?mx\] by auto - show ?thesis unfolding minus_float.rep_eq arctan_minus lb_arctan.simps[where x=x] ub_arctan.simps[where x=x] Let_def if_P[OF \x < 0\] + show ?thesis + unfolding minus_float.rep_eq arctan_minus lb_arctan.simps[where x=x] + ub_arctan.simps[where x=x] Let_def if_P[OF \x < 0\] unfolding atLeastAtMost_iff using bounds[unfolded minus_float.rep_eq arctan_minus] by (simp add: arctan_minus) qed lemma bnds_arctan: "\ (x::real) lx ux. (l, u) = (lb_arctan prec lx, ub_arctan prec ux) \ x \ {lx .. ux} \ l \ arctan x \ arctan x \ u" proof (rule allI, rule allI, rule allI, rule impI) - fix x :: real fix lx ux + fix x :: real + fix lx ux assume "(l, u) = (lb_arctan prec lx, ub_arctan prec ux) \ x \ {lx .. ux}" - hence l: "lb_arctan prec lx = l " and u: "ub_arctan prec ux = u" and x: "x \ {lx .. ux}" by auto - - { from arctan_boundaries[of lx prec, unfolded l] - have "l \ arctan lx" by (auto simp del: lb_arctan.simps) - also have "\ \ arctan x" using x by (auto intro: arctan_monotone') - finally have "l \ arctan x" . - } moreover - { have "arctan x \ arctan ux" using x by (auto intro: arctan_monotone') - also have "\ \ u" using arctan_boundaries[of ux prec, unfolded u] by (auto simp del: ub_arctan.simps) - finally have "arctan x \ u" . - } ultimately show "l \ arctan x \ arctan x \ u" .. + hence l: "lb_arctan prec lx = l " + and u: "ub_arctan prec ux = u" + and x: "x \ {lx .. ux}" + by auto + show "l \ arctan x \ arctan x \ u" + proof + show "l \ arctan x" + proof - + from arctan_boundaries[of lx prec, unfolded l] + have "l \ arctan lx" by (auto simp del: lb_arctan.simps) + also have "\ \ arctan x" using x by (auto intro: arctan_monotone') + finally show ?thesis . + qed + show "arctan x \ u" + proof - + have "arctan x \ arctan ux" using x by (auto intro: arctan_monotone') + also have "\ \ u" using arctan_boundaries[of ux prec, unfolded u] by (auto simp del: ub_arctan.simps) + finally show ?thesis . + qed + qed qed + section "Sinus and Cosinus" subsection "Compute the cosinus and sinus series" @@ -912,15 +1029,15 @@ proof - have "0 \ real (x * x)" by auto let "?f n" = "fact (2 * n) :: nat" - - { fix n - have F: "\m. ((\i. i + 2) ^^ n) m = m + 2 * n" by (induct n) auto - have "?f (Suc n) = ?f n * ((\i. i + 2) ^^ n) 1 * (((\i. i + 2) ^^ n) 1 + 1)" - unfolding F by auto } note f_eq = this - + have f_eq: "?f (Suc n) = ?f n * ((\i. i + 2) ^^ n) 1 * (((\i. i + 2) ^^ n) 1 + 1)" for n + proof - + have "\m. ((\i. i + 2) ^^ n) m = m + 2 * n" by (induct n) auto + then show ?thesis by auto + qed from horner_bounds[where lb="lb_sin_cos_aux prec" and ub="ub_sin_cos_aux prec" and j'=0, OF \0 \ real (x * x)\ f_eq lb_sin_cos_aux.simps ub_sin_cos_aux.simps] - show "?lb" and "?ub" by (auto simp add: power_mult power2_eq_square[of "real x"]) + 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 \ 1" @@ -930,15 +1047,20 @@ lemma one_le_ub_sin_cos_aux: "odd n \ 1 \ 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 \ real x" and "x \ pi / 2" +lemma cos_boundaries: + assumes "0 \ real x" and "x \ pi / 2" shows "cos x \ {(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") - case False hence "real x \ 0" by auto - hence "0 < x" and "0 < real x" using \0 \ real x\ by auto - have "0 < x * x" using \0 < x\ by simp - - { fix x n have "(\ i=0.. i = 0 ..< 2 * n. (if even(i) then ((- 1) ^ (i div 2))/((fact i)) else 0) * x ^ i)" (is "?sum = ?ifsum") + case False + hence "real x \ 0" by auto + hence "0 < x" and "0 < real x" + using \0 \ real x\ by auto + have "0 < x * x" + using \0 < x\ by simp + + have morph_to_if_power: "(\ i=0.. i = 0 ..< 2 * n. (if even(i) then ((- 1) ^ (i div 2))/((fact i)) else 0) * x ^ i)" + (is "?sum = ?ifsum") for x n proof - have "?sum = ?sum + (\ j = 0 ..< n. 0)" by auto also have "\ = @@ -947,9 +1069,8 @@ unfolding sum_split_even_odd atLeast0LessThan .. also have "\ = (\ i = 0 ..< 2 * n. (if even i then (- 1) ^ (i div 2) / ((fact i)) else 0) * x ^ i)" by (rule setsum.cong) auto - finally show ?thesis by assumption - qed } note morph_to_if_power = this - + finally show ?thesis . + qed { fix n :: nat assume "0 < n" hence "0 < 2 * n" by auto @@ -1000,16 +1121,20 @@ } note ub = this and lb } note ub = this(1) and lb = this(2) - have "cos x \ (ub_sin_cos_aux prec (get_odd n) 1 1 (x * x))" using ub[OF odd_pos[OF get_odd] get_odd] . + have "cos x \ (ub_sin_cos_aux prec (get_odd n) 1 1 (x * x))" + using ub[OF odd_pos[OF get_odd] get_odd] . moreover have "(lb_sin_cos_aux prec (get_even n) 1 1 (x * x)) \ cos x" proof (cases "0 < get_even n") - case True show ?thesis using lb[OF True get_even] . + case True + show ?thesis using lb[OF True get_even] . next case False hence "get_even n = 0" by auto - have "- (pi / 2) \ x" by (rule order_trans[OF _ \0 < real x\[THEN less_imp_le]], auto) - with \x \ pi / 2\ - show ?thesis unfolding \get_even n = 0\ lb_sin_cos_aux.simps minus_float.rep_eq zero_float.rep_eq using cos_ge_zero by auto + have "- (pi / 2) \ x" + by (rule order_trans[OF _ \0 < real x\[THEN less_imp_le]]) auto + with \x \ pi / 2\ show ?thesis + unfolding \get_even n = 0\ lb_sin_cos_aux.simps minus_float.rep_eq zero_float.rep_eq + using cos_ge_zero by auto qed ultimately show ?thesis by auto next @@ -1021,18 +1146,21 @@ by simp qed -lemma sin_aux: assumes "0 \ real x" - shows "(x * lb_sin_cos_aux prec n 2 1 (x * x)) \ (\ i=0.. i=0.. (x * ub_sin_cos_aux prec n 2 1 (x * x))" (is "?ub") +lemma sin_aux: + assumes "0 \ real x" + shows "(x * lb_sin_cos_aux prec n 2 1 (x * x)) \ + (\ i=0.. i=0.. + (x * ub_sin_cos_aux prec n 2 1 (x * x))" (is "?ub") proof - have "0 \ real (x * x)" by auto let "?f n" = "fact (2 * n + 1) :: nat" - - { fix n + have f_eq: "?f (Suc n) = ?f n * ((\i. i + 2) ^^ n) 2 * (((\i. i + 2) ^^ n) 2 + 1)" for n + proof - have F: "\m. ((\i. i + 2) ^^ n) m = m + 2 * n" by (induct n) auto - have "?f (Suc n) = ?f n * ((\i. i + 2) ^^ n) 2 * (((\i. i + 2) ^^ n) 2 + 1)" - unfolding F by auto } - note f_eq = this + show ?thesis + unfolding F by auto + qed from horner_bounds[where lb="lb_sin_cos_aux prec" and ub="ub_sin_cos_aux prec" and j'=0, OF \0 \ real (x * x)\ f_eq lb_sin_cos_aux.simps ub_sin_cos_aux.simps] show "?lb" and "?ub" using \0 \ real x\ @@ -1041,25 +1169,32 @@ by (auto intro!: mult_left_mono simp add: power_mult power2_eq_square[of "real x"]) qed -lemma sin_boundaries: assumes "0 \ real x" and "x \ pi / 2" +lemma sin_boundaries: + assumes "0 \ real x" + and "x \ pi / 2" shows "sin x \ {(x * lb_sin_cos_aux prec (get_even n) 2 1 (x * x)) .. (x * ub_sin_cos_aux prec (get_odd n) 2 1 (x * x))}" proof (cases "real x = 0") - case False hence "real x \ 0" by auto - hence "0 < x" and "0 < real x" using \0 \ real x\ by auto - have "0 < x * x" using \0 < x\ by simp - - { fix x::real and n - have "(\j = 0 ..< n. (- 1) ^ (((2 * j + 1) - Suc 0) div 2) / ((fact (2 * j + 1))) * x ^(2 * j + 1)) - = (\ i = 0 ..< 2 * n. (if even(i) then 0 else ((- 1) ^ ((i - Suc 0) div 2))/((fact i))) * x ^ i)" (is "?SUM = _") - proof - - have pow: "!!i. x ^ (2 * i + 1) = x * x ^ (2 * i)" by auto - have "?SUM = (\ j = 0 ..< n. 0) + ?SUM" by auto - also have "\ = (\ i = 0 ..< 2 * n. if even i then 0 else (- 1) ^ ((i - Suc 0) div 2) / ((fact i)) * x ^ i)" - unfolding sum_split_even_odd atLeast0LessThan .. - also have "\ = (\ i = 0 ..< 2 * n. (if even i then 0 else (- 1) ^ ((i - Suc 0) div 2) / ((fact i))) * x ^ i)" - by (rule setsum.cong) auto - finally show ?thesis by assumption - qed } note setsum_morph = this + case False + hence "real x \ 0" by auto + hence "0 < x" and "0 < real x" + using \0 \ real x\ by auto + have "0 < x * x" + using \0 < x\ by simp + + have setsum_morph: "(\j = 0 ..< n. (- 1) ^ (((2 * j + 1) - Suc 0) div 2) / ((fact (2 * j + 1))) * x ^(2 * j + 1)) = + (\ i = 0 ..< 2 * n. (if even(i) then 0 else ((- 1) ^ ((i - Suc 0) div 2))/((fact i))) * x ^ i)" + (is "?SUM = _") for x :: real and n + proof - + have pow: "!!i. x ^ (2 * i + 1) = x * x ^ (2 * i)" + by auto + have "?SUM = (\ j = 0 ..< n. 0) + ?SUM" + by auto + also have "\ = (\ i = 0 ..< 2 * n. if even i then 0 else (- 1) ^ ((i - Suc 0) div 2) / ((fact i)) * x ^ i)" + unfolding sum_split_even_odd atLeast0LessThan .. + also have "\ = (\ i = 0 ..< 2 * n. (if even i then 0 else (- 1) ^ ((i - Suc 0) div 2) / ((fact i))) * x ^ i)" + by (rule setsum.cong) auto + finally show ?thesis . + qed { fix n :: nat assume "0 < n" hence "0 < 2 * n + 1" by auto @@ -1070,14 +1205,20 @@ using Maclaurin_sin_expansion3[OF \0 < 2 * n + 1\ \0 < real x\] unfolding sin_coeff_def atLeast0LessThan by auto - have "?rest = cos t * (- 1) ^ n" unfolding sin_add cos_add real_of_nat_add distrib_right distrib_left by auto + have "?rest = cos t * (- 1) ^ n" + unfolding sin_add cos_add real_of_nat_add distrib_right distrib_left by auto moreover - have "t \ pi / 2" using \t < real x\ and \x \ pi / 2\ by auto - hence "0 \ cos t" using \0 < t\ and cos_ge_zero by auto - ultimately have even: "even n \ 0 \ ?rest" and odd: "odd n \ 0 \ - ?rest " by auto - - have "0 < ?fact" by (simp del: fact_Suc) - have "0 < ?pow" using \0 < real x\ by (rule zero_less_power) + have "t \ pi / 2" + using \t < real x\ and \x \ pi / 2\ by auto + hence "0 \ cos t" + using \0 < t\ and cos_ge_zero by auto + ultimately have even: "even n \ 0 \ ?rest" and odd: "odd n \ 0 \ - ?rest" + by auto + + have "0 < ?fact" + by (simp del: fact_Suc) + have "0 < ?pow" + using \0 < real x\ by (rule zero_less_power) { assume "even n" @@ -1111,15 +1252,20 @@ } note ub = this and lb } note ub = this(1) and lb = this(2) - have "sin x \ (x * ub_sin_cos_aux prec (get_odd n) 2 1 (x * x))" using ub[OF odd_pos[OF get_odd] get_odd] . + have "sin x \ (x * ub_sin_cos_aux prec (get_odd n) 2 1 (x * x))" + using ub[OF odd_pos[OF get_odd] get_odd] . moreover have "(x * lb_sin_cos_aux prec (get_even n) 2 1 (x * x)) \ sin x" proof (cases "0 < get_even n") - case True show ?thesis using lb[OF True get_even] . + case True + show ?thesis + using lb[OF True get_even] . next case False hence "get_even n = 0" by auto with \x \ pi / 2\ \0 \ real x\ - show ?thesis unfolding \get_even n = 0\ ub_sin_cos_aux.simps minus_float.rep_eq using sin_ge_zero by auto + show ?thesis + unfolding \get_even n = 0\ ub_sin_cos_aux.simps minus_float.rep_eq + using sin_ge_zero by auto qed ultimately show ?thesis by auto next @@ -1127,13 +1273,20 @@ 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 + 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 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) + 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 qed + subsection "Compute the cosinus in the entire domain" definition lb_cos :: "nat \ float \ float" where @@ -1152,16 +1305,20 @@ else if x < 1 then half (horner (x * Float 1 (- 1))) else half (half (horner (x * Float 1 (- 2)))))" -lemma lb_cos: assumes "0 \ real x" and "x \ pi" +lemma lb_cos: + assumes "0 \ real x" and "x \ pi" shows "cos x \ {(lb_cos prec x) .. (ub_cos prec x)}" (is "?cos x \ {(?lb x) .. (?ub x) }") proof - - { fix x :: real - have "cos x = cos (x / 2 + x / 2)" by auto + have x_half[symmetric]: "cos x = 2 * cos (x / 2) * cos (x / 2) - 1" for x :: real + proof - + have "cos x = cos (x / 2 + x / 2)" + by auto also have "\ = cos (x / 2) * cos (x / 2) + sin (x / 2) * sin (x / 2) - sin (x / 2) * sin (x / 2) + cos (x / 2) * cos (x / 2) - 1" unfolding cos_add by auto - also have "\ = 2 * cos (x / 2) * cos (x / 2) - 1" by algebra - finally have "cos x = 2 * cos (x / 2) * cos (x / 2) - 1" . - } note x_half = this[symmetric] + also have "\ = 2 * cos (x / 2) * cos (x / 2) - 1" + by algebra + finally show ?thesis . + qed have "\ x < 0" using \0 \ real x\ by auto let "?ub_horner x" = "ub_sin_cos_aux prec (get_odd (prec div 4 + 1)) 1 1 (x * x)" @@ -1171,44 +1328,60 @@ show ?thesis proof (cases "x < Float 1 (- 1)") - case True hence "x \ pi / 2" using pi_ge_two by auto - show ?thesis unfolding lb_cos_def[where x=x] ub_cos_def[where x=x] if_not_P[OF \\ x < 0\] if_P[OF \x < Float 1 (- 1)\] Let_def + case True + hence "x \ pi / 2" + using pi_ge_two by auto + show ?thesis + unfolding lb_cos_def[where x=x] ub_cos_def[where x=x] + if_not_P[OF \\ x < 0\] if_P[OF \x < Float 1 (- 1)\] Let_def using cos_boundaries[OF \0 \ real x\ \x \ pi / 2\] . next case False { fix y x :: float let ?x2 = "(x * Float 1 (- 1))" assume "y \ cos ?x2" and "-pi \ x" and "x \ pi" - hence "- (pi / 2) \ ?x2" and "?x2 \ pi / 2" using pi_ge_two unfolding Float_num by auto - hence "0 \ cos ?x2" by (rule cos_ge_zero) + hence "- (pi / 2) \ ?x2" and "?x2 \ pi / 2" + using pi_ge_two unfolding Float_num by auto + hence "0 \ cos ?x2" + by (rule cos_ge_zero) have "(?lb_half y) \ cos x" proof (cases "y < 0") - case True show ?thesis using cos_ge_minus_one unfolding if_P[OF True] by auto + case True + show ?thesis + using cos_ge_minus_one unfolding if_P[OF True] by auto next case False hence "0 \ real y" by auto from mult_mono[OF \y \ cos ?x2\ \y \ cos ?x2\ \0 \ cos ?x2\ this] have "real y * real y \ cos ?x2 * cos ?x2" . - hence "2 * real y * real y \ 2 * cos ?x2 * cos ?x2" by auto - hence "2 * real y * real y - 1 \ 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 + hence "2 * real y * real y \ 2 * cos ?x2 * cos ?x2" + by auto + hence "2 * real y * real y - 1 \ 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 intro!: float_plus_down_le) qed } note lb_half = this { fix y x :: float let ?x2 = "(x * Float 1 (- 1))" assume ub: "cos ?x2 \ y" and "- pi \ x" and "x \ pi" - hence "- (pi / 2) \ ?x2" and "?x2 \ pi / 2" using pi_ge_two unfolding Float_num by auto + hence "- (pi / 2) \ ?x2" and "?x2 \ pi / 2" + using pi_ge_two unfolding Float_num by auto hence "0 \ cos ?x2" by (rule cos_ge_zero) have "cos x \ (?ub_half y)" proof - - have "0 \ real y" using \0 \ cos ?x2\ ub by (rule order_trans) + have "0 \ real y" + using \0 \ cos ?x2\ ub by (rule order_trans) from mult_mono[OF ub ub this \0 \ cos ?x2\] have "cos ?x2 * cos ?x2 \ real y * real y" . - hence "2 * cos ?x2 * cos ?x2 \ 2 * real y * real y" by auto - hence "2 * cos (x / 2) * cos (x / 2) - 1 \ 2 * real y * real y - 1" unfolding Float_num by auto - thus ?thesis unfolding x_half Float_num + hence "2 * cos ?x2 * cos ?x2 \ 2 * real y * real y" + by auto + hence "2 * cos (x / 2) * cos (x / 2) - 1 \ 2 * real y * real y - 1" + unfolding Float_num by auto + thus ?thesis + unfolding x_half Float_num by (auto intro!: float_plus_up_le) qed } note ub_half = this @@ -1216,55 +1389,76 @@ let ?x2 = "x * Float 1 (- 1)" let ?x4 = "x * Float 1 (- 1) * Float 1 (- 1)" - have "-pi \ x" using pi_ge_zero[THEN le_imp_neg_le, unfolded minus_zero] \0 \ real x\ by (rule order_trans) + have "-pi \ x" + using pi_ge_zero[THEN le_imp_neg_le, unfolded minus_zero] \0 \ real x\ + by (rule order_trans) show ?thesis proof (cases "x < 1") - case True hence "real x \ 1" by auto - have "0 \ real ?x2" and "?x2 \ pi / 2" using pi_ge_two \0 \ real x\ using assms by auto + case True + hence "real x \ 1" by auto + have "0 \ real ?x2" and "?x2 \ pi / 2" + using pi_ge_two \0 \ real x\ using assms by auto from cos_boundaries[OF this] - have lb: "(?lb_horner ?x2) \ ?cos ?x2" and ub: "?cos ?x2 \ (?ub_horner ?x2)" by auto + have lb: "(?lb_horner ?x2) \ ?cos ?x2" and ub: "?cos ?x2 \ (?ub_horner ?x2)" + by auto have "(?lb x) \ ?cos x" proof - from lb_half[OF lb \-pi \ x\ \x \ pi\] - show ?thesis unfolding lb_cos_def[where x=x] Let_def using \\ x < 0\ \\ x < Float 1 (- 1)\ \x < 1\ by auto + show ?thesis + unfolding lb_cos_def[where x=x] Let_def + using \\ x < 0\ \\ x < Float 1 (- 1)\ \x < 1\ by auto qed moreover have "?cos x \ (?ub x)" proof - from ub_half[OF ub \-pi \ x\ \x \ pi\] - show ?thesis unfolding ub_cos_def[where x=x] Let_def using \\ x < 0\ \\ x < Float 1 (- 1)\ \x < 1\ by auto + show ?thesis + unfolding ub_cos_def[where x=x] Let_def + using \\ x < 0\ \\ x < Float 1 (- 1)\ \x < 1\ by auto qed ultimately show ?thesis by auto next case False - have "0 \ real ?x4" and "?x4 \ pi / 2" using pi_ge_two \0 \ real x\ \x \ pi\ unfolding Float_num by auto + have "0 \ real ?x4" and "?x4 \ pi / 2" + using pi_ge_two \0 \ real x\ \x \ pi\ unfolding Float_num by auto from cos_boundaries[OF this] - have lb: "(?lb_horner ?x4) \ ?cos ?x4" and ub: "?cos ?x4 \ (?ub_horner ?x4)" by auto - - have eq_4: "?x2 * Float 1 (- 1) = x * Float 1 (- 2)" by transfer simp + have lb: "(?lb_horner ?x4) \ ?cos ?x4" and ub: "?cos ?x4 \ (?ub_horner ?x4)" + by auto + + have eq_4: "?x2 * Float 1 (- 1) = x * Float 1 (- 2)" + by transfer simp have "(?lb x) \ ?cos x" proof - - have "-pi \ ?x2" and "?x2 \ pi" using pi_ge_two \0 \ real x\ \x \ pi\ by auto + have "-pi \ ?x2" and "?x2 \ pi" + using pi_ge_two \0 \ real x\ \x \ pi\ by auto from lb_half[OF lb_half[OF lb this] \-pi \ x\ \x \ pi\, unfolded eq_4] - show ?thesis unfolding lb_cos_def[where x=x] if_not_P[OF \\ x < 0\] if_not_P[OF \\ x < Float 1 (- 1)\] if_not_P[OF \\ x < 1\] Let_def . + show ?thesis + unfolding lb_cos_def[where x=x] if_not_P[OF \\ x < 0\] + if_not_P[OF \\ x < Float 1 (- 1)\] if_not_P[OF \\ x < 1\] Let_def . qed moreover have "?cos x \ (?ub x)" proof - - have "-pi \ ?x2" and "?x2 \ pi" using pi_ge_two \0 \ real x\ \ x \ pi\ by auto + have "-pi \ ?x2" and "?x2 \ pi" + using pi_ge_two \0 \ real x\ \ x \ pi\ by auto from ub_half[OF ub_half[OF ub this] \-pi \ x\ \x \ pi\, unfolded eq_4] - show ?thesis unfolding ub_cos_def[where x=x] if_not_P[OF \\ x < 0\] if_not_P[OF \\ x < Float 1 (- 1)\] if_not_P[OF \\ x < 1\] Let_def . + show ?thesis + unfolding ub_cos_def[where x=x] if_not_P[OF \\ x < 0\] + if_not_P[OF \\ x < Float 1 (- 1)\] if_not_P[OF \\ x < 1\] Let_def . qed ultimately show ?thesis by auto qed qed qed -lemma lb_cos_minus: assumes "-pi \ x" and "real x \ 0" +lemma lb_cos_minus: + assumes "-pi \ x" + and "real x \ 0" shows "cos (real(-x)) \ {(lb_cos prec (-x)) .. (ub_cos prec (-x))}" proof - - have "0 \ real (-x)" and "(-x) \ pi" using \-pi \ x\ \real x \ 0\ by auto + have "0 \ real (-x)" and "(-x) \ pi" + using \-pi \ x\ \real x \ 0\ by auto from lb_cos[OF this] show ?thesis . qed @@ -1282,33 +1476,46 @@ else if -2 * lpi \ lx \ ux \ 0 then (Float (- 1) 0, max (ub_cos prec (lx + 2 * lpi)) (ub_cos prec (-ux))) else (Float (- 1) 0, Float 1 0))" -lemma floor_int: - obtains k :: int where "real k = (floor_fl f)" +lemma floor_int: obtains k :: int where "real k = (floor_fl f)" by (simp add: floor_fl_def) -lemma cos_periodic_nat[simp]: fixes n :: nat shows "cos (x + n * (2 * pi)) = cos x" +lemma cos_periodic_nat[simp]: + fixes n :: nat + shows "cos (x + n * (2 * pi)) = cos x" proof (induct n arbitrary: x) + case 0 + then show ?case by simp +next case (Suc n) have split_pi_off: "x + (Suc n) * (2 * pi) = (x + n * (2 * pi)) + 2 * pi" unfolding Suc_eq_plus1 real_of_nat_add real_of_one distrib_right by auto - show ?case unfolding split_pi_off using Suc by auto -qed auto - -lemma cos_periodic_int[simp]: fixes i :: int shows "cos (x + i * (2 * pi)) = cos x" + show ?case + unfolding split_pi_off using Suc by auto +qed + +lemma cos_periodic_int[simp]: + fixes i :: int + shows "cos (x + i * (2 * pi)) = cos x" proof (cases "0 \ i") - case True hence i_nat: "real i = nat i" by auto - show ?thesis unfolding i_nat by auto + case True + hence i_nat: "real i = nat i" by auto + show ?thesis + unfolding i_nat by auto next - case False hence i_nat: "i = - real (nat (-i))" by auto - have "cos x = cos (x + i * (2 * pi) - i * (2 * pi))" by auto + case False + hence i_nat: "i = - real (nat (-i))" by auto + have "cos x = cos (x + i * (2 * pi) - i * (2 * pi))" + by auto also have "\ = cos (x + i * (2 * pi))" unfolding i_nat mult_minus_left diff_minus_eq_add by (rule cos_periodic_nat) finally show ?thesis by auto qed -lemma bnds_cos: "\ (x::real) lx ux. (l, u) = bnds_cos prec lx ux \ x \ {lx .. ux} \ l \ cos x \ cos x \ u" -proof ((rule allI | rule impI | erule conjE) +) - fix x :: real fix lx ux +lemma bnds_cos: "\(x::real) lx ux. (l, u) = + bnds_cos prec lx ux \ x \ {lx .. ux} \ l \ cos x \ cos x \ u" +proof (rule allI | rule impI | erule conjE)+ + fix x :: real + fix lx ux assume bnds: "(l, u) = bnds_cos prec lx ux" and x: "x \ {lx .. ux}" let ?lpi = "float_round_down prec (lb_pi prec)" @@ -1319,11 +1526,13 @@ 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 . + obtain k :: int where k: "k = real ?k" + by (rule floor_int) have upi: "pi \ ?upi" and lpi: "?lpi \ pi" using float_round_up[of "ub_pi prec" prec] pi_boundaries[of prec] - float_round_down[of prec "lb_pi prec"] by auto + float_round_down[of prec "lb_pi prec"] + by auto hence "lx + ?lx2 \ x - k * (2 * pi) \ x - k * (2 * pi) \ ux + ?ux2" using x by (cases "k = 0") @@ -1397,134 +1606,152 @@ show "l \ cos x \ cos x \ u" proof (cases "- ?lpi \ ?lx \ ?ux \ 0") - case True with bnds - have l: "l = lb_cos prec (-?lx)" - and u: "u = ub_cos prec (-?ux)" + case True + with bnds have l: "l = lb_cos prec (-?lx)" and u: "u = ub_cos prec (-?ux)" by (auto simp add: bnds_cos_def Let_def) - from True lpi[THEN le_imp_neg_le] lx ux - have "- pi \ x - k * (2 * pi)" - and "x - k * (2 * pi) \ 0" - by auto - with True negative_ux negative_lx - show ?thesis unfolding l u by simp - next case False note 1 = this show ?thesis - proof (cases "0 \ ?lx \ ?ux \ ?lpi") - case True with bnds 1 - have l: "l = lb_cos prec ?ux" - and u: "u = ub_cos prec ?lx" - by (auto simp add: bnds_cos_def Let_def) - - from True lpi lx ux - have "0 \ x - k * (2 * pi)" - and "x - k * (2 * pi) \ pi" + have "- pi \ x - k * (2 * pi)" and "x - k * (2 * pi) \ 0" by auto - with True positive_ux positive_lx - show ?thesis unfolding l u by simp - next case False note 2 = this show ?thesis - proof (cases "- ?lpi \ ?lx \ ?ux \ ?lpi") - case True note Cond = this with bnds 1 2 - have l: "l = min (lb_cos prec (-?lx)) (lb_cos prec ?ux)" - and u: "u = Float 1 0" - by (auto simp add: bnds_cos_def Let_def) - - show ?thesis unfolding u l using negative_lx positive_ux Cond - by (cases "x - k * (2 * pi) < 0") (auto simp add: real_of_float_min) - - next case False note 3 = this show ?thesis - proof (cases "0 \ ?lx \ ?ux \ 2 * ?lpi") - case True note Cond = this with bnds 1 2 3 - have l: "l = Float (- 1) 0" - and u: "u = max (ub_cos prec ?lx) (ub_cos prec (- (?ux - 2 * ?lpi)))" - by (auto simp add: bnds_cos_def Let_def) - - have "cos x \ real u" - proof (cases "x - k * (2 * pi) < pi") - case True hence "x - k * (2 * pi) \ pi" by simp - from positive_lx[OF Cond[THEN conjunct1] this] - show ?thesis unfolding u by (simp add: real_of_float_max) + with True negative_ux negative_lx show ?thesis + unfolding l u by simp + next + case 1: False + show ?thesis + proof (cases "0 \ ?lx \ ?ux \ ?lpi") + case True with bnds 1 + have l: "l = lb_cos prec ?ux" + and u: "u = ub_cos prec ?lx" + by (auto simp add: bnds_cos_def Let_def) + from True lpi lx ux + have "0 \ x - k * (2 * pi)" and "x - k * (2 * pi) \ pi" + by auto + with True positive_ux positive_lx show ?thesis + unfolding l u by simp next - case False hence "pi \ x - k * (2 * pi)" by simp - hence pi_x: "- pi \ x - k * (2 * pi) - 2 * pi" by simp - - have "?ux \ 2 * pi" using Cond lpi by auto - hence "x - k * (2 * pi) - 2 * pi \ 0" using ux by simp - - have ux_0: "real (?ux - 2 * ?lpi) \ 0" - using Cond by auto - - from 2 and Cond have "\ ?ux \ ?lpi" by auto - hence "- ?lpi \ ?ux - 2 * ?lpi" by auto - hence pi_ux: "- pi \ (?ux - 2 * ?lpi)" - using lpi[THEN le_imp_neg_le] by auto - - have x_le_ux: "x - k * (2 * pi) - 2 * pi \ (?ux - 2 * ?lpi)" - using ux lpi by auto - have "cos x = cos (x + (-k) * (2 * pi) + (-1::int) * (2 * pi))" - unfolding cos_periodic_int .. - also have "\ \ cos ((?ux - 2 * ?lpi))" - using cos_monotone_minus_pi_0'[OF pi_x x_le_ux ux_0] - by (simp only: minus_float.rep_eq real_of_int_minus real_of_one - mult_minus_left mult_1_left) simp - also have "\ = cos ((- (?ux - 2 * ?lpi)))" - unfolding uminus_float.rep_eq cos_minus .. - also have "\ \ (ub_cos prec (- (?ux - 2 * ?lpi)))" - using lb_cos_minus[OF pi_ux ux_0] by simp - finally show ?thesis unfolding u by (simp add: real_of_float_max) + case 2: False + show ?thesis + proof (cases "- ?lpi \ ?lx \ ?ux \ ?lpi") + case Cond: True + with bnds 1 2 have l: "l = min (lb_cos prec (-?lx)) (lb_cos prec ?ux)" + and u: "u = Float 1 0" + by (auto simp add: bnds_cos_def Let_def) + show ?thesis + unfolding u l using negative_lx positive_ux Cond + by (cases "x - k * (2 * pi) < 0") (auto simp add: real_of_float_min) + next + case 3: False + show ?thesis + proof (cases "0 \ ?lx \ ?ux \ 2 * ?lpi") + case Cond: True + with bnds 1 2 3 + have l: "l = Float (- 1) 0" + and u: "u = max (ub_cos prec ?lx) (ub_cos prec (- (?ux - 2 * ?lpi)))" + by (auto simp add: bnds_cos_def Let_def) + + have "cos x \ real u" + proof (cases "x - k * (2 * pi) < pi") + case True + hence "x - k * (2 * pi) \ pi" by simp + from positive_lx[OF Cond[THEN conjunct1] this] show ?thesis + unfolding u by (simp add: real_of_float_max) + next + case False + hence "pi \ x - k * (2 * pi)" by simp + hence pi_x: "- pi \ x - k * (2 * pi) - 2 * pi" by simp + + have "?ux \ 2 * pi" + using Cond lpi by auto + hence "x - k * (2 * pi) - 2 * pi \ 0" + using ux by simp + + have ux_0: "real (?ux - 2 * ?lpi) \ 0" + using Cond by auto + + from 2 and Cond have "\ ?ux \ ?lpi" by auto + hence "- ?lpi \ ?ux - 2 * ?lpi" by auto + hence pi_ux: "- pi \ (?ux - 2 * ?lpi)" + using lpi[THEN le_imp_neg_le] by auto + + have x_le_ux: "x - k * (2 * pi) - 2 * pi \ (?ux - 2 * ?lpi)" + using ux lpi by auto + have "cos x = cos (x + (-k) * (2 * pi) + (-1::int) * (2 * pi))" + unfolding cos_periodic_int .. + also have "\ \ cos ((?ux - 2 * ?lpi))" + using cos_monotone_minus_pi_0'[OF pi_x x_le_ux ux_0] + by (simp only: minus_float.rep_eq real_of_int_minus real_of_one + mult_minus_left mult_1_left) simp + also have "\ = cos ((- (?ux - 2 * ?lpi)))" + unfolding uminus_float.rep_eq cos_minus .. + also have "\ \ (ub_cos prec (- (?ux - 2 * ?lpi)))" + using lb_cos_minus[OF pi_ux ux_0] by simp + finally show ?thesis unfolding u by (simp add: real_of_float_max) + qed + thus ?thesis unfolding l by auto + next + case 4: False + show ?thesis + proof (cases "-2 * ?lpi \ ?lx \ ?ux \ 0") + case Cond: True + with bnds 1 2 3 4 have l: "l = Float (- 1) 0" + and u: "u = max (ub_cos prec (?lx + 2 * ?lpi)) (ub_cos prec (-?ux))" + by (auto simp add: bnds_cos_def Let_def) + + have "cos x \ u" + proof (cases "-pi < x - k * (2 * pi)") + case True + hence "-pi \ x - k * (2 * pi)" by simp + from negative_ux[OF this Cond[THEN conjunct2]] show ?thesis + unfolding u by (simp add: real_of_float_max) + next + case False + hence "x - k * (2 * pi) \ -pi" by simp + hence pi_x: "x - k * (2 * pi) + 2 * pi \ pi" by simp + + have "-2 * pi \ ?lx" using Cond lpi by auto + + hence "0 \ x - k * (2 * pi) + 2 * pi" using lx by simp + + have lx_0: "0 \ real (?lx + 2 * ?lpi)" + using Cond lpi by auto + + from 1 and Cond have "\ -?lpi \ ?lx" by auto + hence "?lx + 2 * ?lpi \ ?lpi" by auto + hence pi_lx: "(?lx + 2 * ?lpi) \ pi" + using lpi[THEN le_imp_neg_le] by auto + + have lx_le_x: "(?lx + 2 * ?lpi) \ x - k * (2 * pi) + 2 * pi" + using lx lpi by auto + + have "cos x = cos (x + (-k) * (2 * pi) + (1 :: int) * (2 * pi))" + unfolding cos_periodic_int .. + also have "\ \ cos ((?lx + 2 * ?lpi))" + using cos_monotone_0_pi_le[OF lx_0 lx_le_x pi_x] + by (simp only: minus_float.rep_eq real_of_int_minus real_of_one + mult_minus_left mult_1_left) simp + also have "\ \ (ub_cos prec (?lx + 2 * ?lpi))" + using lb_cos[OF lx_0 pi_lx] by simp + finally show ?thesis unfolding u by (simp add: real_of_float_max) + qed + thus ?thesis unfolding l by auto + next + case False + with bnds 1 2 3 4 show ?thesis + by (auto simp add: bnds_cos_def Let_def) + qed + qed + qed qed - thus ?thesis unfolding l by auto - next case False note 4 = this show ?thesis - proof (cases "-2 * ?lpi \ ?lx \ ?ux \ 0") - case True note Cond = this with bnds 1 2 3 4 - have l: "l = Float (- 1) 0" - and u: "u = max (ub_cos prec (?lx + 2 * ?lpi)) (ub_cos prec (-?ux))" - by (auto simp add: bnds_cos_def Let_def) - - have "cos x \ u" - proof (cases "-pi < x - k * (2 * pi)") - case True hence "-pi \ x - k * (2 * pi)" by simp - from negative_ux[OF this Cond[THEN conjunct2]] - show ?thesis unfolding u by (simp add: real_of_float_max) - next - case False hence "x - k * (2 * pi) \ -pi" by simp - hence pi_x: "x - k * (2 * pi) + 2 * pi \ pi" by simp - - have "-2 * pi \ ?lx" using Cond lpi by auto - - hence "0 \ x - k * (2 * pi) + 2 * pi" using lx by simp - - have lx_0: "0 \ real (?lx + 2 * ?lpi)" - using Cond lpi by auto - - from 1 and Cond have "\ -?lpi \ ?lx" by auto - hence "?lx + 2 * ?lpi \ ?lpi" by auto - hence pi_lx: "(?lx + 2 * ?lpi) \ pi" - using lpi[THEN le_imp_neg_le] by auto - - have lx_le_x: "(?lx + 2 * ?lpi) \ x - k * (2 * pi) + 2 * pi" - using lx lpi by auto - - have "cos x = cos (x + (-k) * (2 * pi) + (1 :: int) * (2 * pi))" - unfolding cos_periodic_int .. - also have "\ \ cos ((?lx + 2 * ?lpi))" - using cos_monotone_0_pi_le[OF lx_0 lx_le_x pi_x] - by (simp only: minus_float.rep_eq real_of_int_minus real_of_one - mult_minus_left mult_1_left) simp - also have "\ \ (ub_cos prec (?lx + 2 * ?lpi))" - using lb_cos[OF lx_0 pi_lx] by simp - finally show ?thesis unfolding u by (simp add: real_of_float_max) - qed - thus ?thesis unfolding l by auto - next - case False with bnds 1 2 3 4 show ?thesis by (auto simp add: bnds_cos_def Let_def) - qed qed qed qed qed + qed qed + section "Exponential function" subsection "Compute the series of the exponential function" -fun ub_exp_horner :: "nat \ nat \ nat \ nat \ float \ float" and lb_exp_horner :: "nat \ nat \ nat \ nat \ float \ float" where +fun ub_exp_horner :: "nat \ nat \ nat \ nat \ float \ float" + and lb_exp_horner :: "nat \ nat \ nat \ nat \ float \ float" +where "ub_exp_horner prec 0 i k x = 0" | "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))" | @@ -1532,18 +1759,24 @@ "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 \ 0" - shows "exp x \ { lb_exp_horner prec (get_even n) 1 1 x .. ub_exp_horner prec (get_odd n) 1 1 x }" +lemma bnds_exp_horner: + assumes "real x \ 0" + shows "exp x \ {lb_exp_horner prec (get_even n) 1 1 x .. ub_exp_horner prec (get_odd n) 1 1 x}" proof - - { fix n - have F: "\ m. ((\i. i + 1) ^^ n) m = n + m" by (induct n, auto) - have "fact (Suc n) = fact n * ((\i::nat. i + 1) ^^ n) 1" unfolding F by auto - } note f_eq = this + have f_eq: "fact (Suc n) = fact n * ((\i::nat. i + 1) ^^ n) 1" for n + proof - + have F: "\ m. ((\i. i + 1) ^^ n) m = n + m" + by (induct n) auto + show ?thesis + unfolding F by auto + qed note bounds = horner_bounds_nonpos[where f="fact" and lb="lb_exp_horner prec" and ub="ub_exp_horner prec" and j'=0 and s=1, OF assms f_eq lb_exp_horner.simps ub_exp_horner.simps] - { have "lb_exp_horner prec (get_even n) 1 1 x \ (\j = 0.. exp x" + proof - + have "lb_exp_horner prec (get_even n) 1 1 x \ (\j = 0.. \ exp x" proof - @@ -1553,9 +1786,11 @@ by (auto simp: zero_le_even_power) ultimately show ?thesis using get_odd exp_gt_zero by auto qed - finally have "lb_exp_horner prec (get_even n) 1 1 x \ exp x" . - } moreover - { + finally show ?thesis . + qed + moreover + have "exp x \ ub_exp_horner prec (get_odd n) 1 1 x" + proof - have x_less_zero: "real x ^ get_odd n \ 0" proof (cases "real x = 0") case True @@ -1565,8 +1800,8 @@ case False hence "real x < 0" using \real x \ 0\ by auto show ?thesis by (rule less_imp_le, auto simp add: power_less_zero_eq \real x < 0\) qed - - obtain t where "\t\ \ \real x\" and "exp x = (\m = 0..t\ \ \real x\" + and "exp x = (\m = 0.. 0" by (auto intro!: mult_nonneg_nonpos divide_nonpos_pos simp add: x_less_zero) @@ -1574,11 +1809,13 @@ using get_odd exp_gt_zero by auto also have "\ \ ub_exp_horner prec (get_odd n) 1 1 x" using bounds(2) by auto - finally have "exp x \ ub_exp_horner prec (get_odd n) 1 1 x" . - } ultimately show ?thesis by auto + finally show ?thesis . + qed + ultimately show ?thesis by auto qed -lemma ub_exp_horner_nonneg: "real x \ 0 \ 0 \ real (ub_exp_horner prec (get_odd n) (Suc 0) (Suc 0) x)" +lemma ub_exp_horner_nonneg: "real x \ 0 \ + 0 \ 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 @@ -1603,74 +1840,90 @@ (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 + by pat_completeness auto termination -by (relation "measure (\ v. let (prec, x) = case_sum id id v in (if 0 < x then 1 else 0))", auto) + by (relation "measure (\ 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) \ exp (- 1)" proof - have eq4: "4 = Suc (Suc (Suc (Suc 0)))" by auto - - have "1 / 4 = (Float 1 (- 2))" unfolding Float_num by auto + have "1 / 4 = (Float 1 (- 2))" + unfolding Float_num by auto also have "\ \ lb_exp_horner 3 (get_even 3) 1 1 (- 1)" by code_simp - also have "\ \ exp (- 1 :: float)" using bnds_exp_horner[where x="- 1"] by auto - finally show ?thesis by simp + also have "\ \ exp (- 1 :: float)" + using bnds_exp_horner[where x="- 1"] by auto + finally show ?thesis + by simp qed -lemma lb_exp_pos: assumes "\ 0 < x" shows "0 < lb_exp prec x" +lemma lb_exp_pos: + assumes "\ 0 < x" + shows "0 < lb_exp prec x" proof - let "?lb_horner x" = "lb_exp_horner prec (get_even (prec + 2)) 1 1 x" - let "?horner x" = "let y = ?lb_horner x in if y \ 0 then Float 1 (- 2) else y" - have pos_horner: "\ x. 0 < ?horner x" unfolding Let_def by (cases "?lb_horner x \ 0", auto) - moreover { fix x :: float fix num :: nat + let "?horner x" = "let y = ?lb_horner x in if y \ 0 then Float 1 (- 2) else y" + have pos_horner: "0 < ?horner x" for x + unfolding Let_def by (cases "?lb_horner x \ 0") auto + moreover have "0 < real ((?horner x) ^ num)" for x :: float and num :: nat + proof - have "0 < real (?horner x) ^ num" using \0 < ?horner x\ by simp also have "\ = (?horner x) ^ num" by auto - finally have "0 < real ((?horner x) ^ num)" . - } + finally show ?thesis . + qed ultimately show ?thesis unfolding lb_exp.simps if_not_P[OF \\ 0 < x\] Let_def - 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) + 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 \ 0" +lemma exp_boundaries': + assumes "x \ 0" shows "exp x \ { (lb_exp prec x) .. (ub_exp prec x)}" proof - let "?lb_exp_horner x" = "lb_exp_horner prec (get_even (prec + 2)) 1 1 x" let "?ub_exp_horner x" = "ub_exp_horner prec (get_odd (prec + 2)) 1 1 x" - have "real x \ 0" and "\ x > 0" using \x \ 0\ by auto + have "real x \ 0" and "\ x > 0" + using \x \ 0\ by auto show ?thesis proof (cases "x < - 1") - case False hence "- 1 \ real x" by auto + case False + hence "- 1 \ real x" by auto show ?thesis proof (cases "?lb_exp_horner x \ 0") - from \\ x < - 1\ have "- 1 \ real x" by auto - hence "exp (- 1) \ exp x" unfolding exp_le_cancel_iff . - from order_trans[OF exp_m1_ge_quarter this] - have "Float 1 (- 2) \ exp x" unfolding Float_num . - moreover case True - ultimately show ?thesis using bnds_exp_horner \real x \ 0\ \\ x > 0\ \\ x < - 1\ by auto + case True + from \\ x < - 1\ + have "- 1 \ real x" by auto + hence "exp (- 1) \ exp x" + unfolding exp_le_cancel_iff . + from order_trans[OF exp_m1_ge_quarter this] have "Float 1 (- 2) \ exp x" + unfolding Float_num . + with True show ?thesis + using bnds_exp_horner \real x \ 0\ \\ x > 0\ \\ x < - 1\ by auto next - case False thus ?thesis using bnds_exp_horner \real x \ 0\ \\ x > 0\ \\ x < - 1\ by (auto simp add: Let_def) + case False + thus ?thesis + using bnds_exp_horner \real x \ 0\ \\ x > 0\ \\ x < - 1\ by (auto simp add: Let_def) qed next case True - let ?num = "nat (- int_floor_fl x)" - have "real (int_floor_fl x) < - 1" using int_floor_fl[of x] \x < - 1\ - by simp + have "real (int_floor_fl x) < - 1" + using int_floor_fl[of x] \x < - 1\ by simp hence "real (int_floor_fl x) < 0" by simp hence "int_floor_fl x < 0" by auto hence "1 \ - int_floor_fl x" by auto hence "0 < nat (- int_floor_fl x)" by auto hence "0 < ?num" by auto hence "real ?num \ 0" by auto - have num_eq: "real ?num = - int_floor_fl x" using \0 < nat (- int_floor_fl x)\ by auto - have "0 < - int_floor_fl x" using \0 < ?num\[unfolded real_of_nat_less_iff[symmetric]] by simp - hence "real (int_floor_fl x) < 0" unfolding less_float_def by auto + have num_eq: "real ?num = - int_floor_fl x" + using \0 < nat (- int_floor_fl x)\ by auto + have "0 < - int_floor_fl x" + using \0 < ?num\[unfolded real_of_nat_less_iff[symmetric]] by simp + 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)" @@ -1683,17 +1936,20 @@ using float_divr_nonpos_pos_upper_bound[OF \real x \ 0\ \0 \ real (- floor_fl x)\] unfolding less_eq_float_def zero_float.rep_eq . - have "exp x = exp (?num * (x / ?num))" using \real ?num \ 0\ by auto - also have "\ = exp (x / ?num) ^ ?num" unfolding exp_real_of_nat_mult .. - also have "\ \ exp (float_divr prec x (- floor_fl x)) ^ ?num" unfolding num_eq fl_eq + have "exp x = exp (?num * (x / ?num))" + using \real ?num \ 0\ by auto + also have "\ = exp (x / ?num) ^ ?num" + unfolding exp_real_of_nat_mult .. + also have "\ \ exp (float_divr prec x (- floor_fl x)) ^ ?num" + unfolding num_eq fl_eq by (rule power_mono, rule exp_le_cancel_iff[THEN iffD2], rule float_divr) auto also have "\ \ (?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) also have "\ \ 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 \\ 0 < x\] if_P[OF \x < - 1\] floor_fl_def Let_def - . + finally show ?thesis + unfolding ub_exp.simps if_not_P[OF \\ 0 < x\] if_P[OF \x < - 1\] floor_fl_def Let_def . qed moreover have "lb_exp prec x \ exp x" @@ -1703,39 +1959,54 @@ show ?thesis proof (cases "?horner \ 0") - case False hence "0 \ real ?horner" by auto + case False + hence "0 \ real ?horner" by auto have div_less_zero: "real (float_divl prec x (- floor_fl x)) \ 0" - using \real (floor_fl x) < 0\ \real x \ 0\ by (auto intro!: order_trans[OF float_divl] divide_nonpos_neg) + using \real (floor_fl x) < 0\ \real x \ 0\ + by (auto intro!: order_trans[OF float_divl] divide_nonpos_neg) have "(?lb_exp_horner (float_divl prec x (- floor_fl x))) ^ ?num \ exp (float_divl prec x (- floor_fl x)) ^ ?num" - using \0 \ real ?horner\[unfolded floor_fl_def[symmetric]] bnds_exp_horner[OF div_less_zero, unfolded atLeastAtMost_iff, THEN conjunct1] by (auto intro!: power_mono) - also have "\ \ exp (x / ?num) ^ ?num" unfolding num_eq fl_eq + using \0 \ real ?horner\[unfolded floor_fl_def[symmetric]] + bnds_exp_horner[OF div_less_zero, unfolded atLeastAtMost_iff, THEN conjunct1] + by (auto intro!: power_mono) + also have "\ \ exp (x / ?num) ^ ?num" + unfolding num_eq fl_eq using float_divl by (auto intro!: power_mono simp del: uminus_float.rep_eq) - also have "\ = exp (?num * (x / ?num))" unfolding exp_real_of_nat_mult .. - also have "\ = exp x" using \real ?num \ 0\ by auto - finally show ?thesis using False - unfolding lb_exp.simps if_not_P[OF \\ 0 < x\] if_P[OF \x < - 1\] int_floor_fl_def Let_def if_not_P[OF False] + also have "\ = exp (?num * (x / ?num))" + unfolding exp_real_of_nat_mult .. + also have "\ = exp x" + using \real ?num \ 0\ by auto + finally show ?thesis + using False + unfolding lb_exp.simps if_not_P[OF \\ 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 \ (Float 1 (- 2)) ^ ?num" - by (metis Float_le_zero_iff less_imp_le linorder_not_less not_numeral_le_zero numeral_One power_down_fl) + by (metis Float_le_zero_iff less_imp_le linorder_not_less + not_numeral_le_zero numeral_One power_down_fl) then have "power_down_fl prec (Float 1 (- 2)) ?num \ real (Float 1 (- 2)) ^ ?num" by simp also - have "real (floor_fl x) \ 0" and "real (floor_fl x) \ 0" using \real (floor_fl x) < 0\ by auto + have "real (floor_fl x) \ 0" and "real (floor_fl x) \ 0" + using \real (floor_fl x) < 0\ by auto from divide_right_mono_neg[OF floor_fl[of x] \real (floor_fl x) \ 0\, unfolded divide_self[OF \real (floor_fl x) \ 0\]] - have "- 1 \ x / (- floor_fl x)" unfolding minus_float.rep_eq by auto + have "- 1 \ x / (- floor_fl x)" + unfolding minus_float.rep_eq by auto from order_trans[OF exp_m1_ge_quarter this[unfolded exp_le_cancel_iff[where x="- 1", symmetric]]] - have "Float 1 (- 2) \ exp (x / (- floor_fl x))" unfolding Float_num . + have "Float 1 (- 2) \ exp (x / (- floor_fl x))" + unfolding Float_num . hence "real (Float 1 (- 2)) ^ ?num \ exp (x / (- floor_fl x)) ^ ?num" by (metis Float_num(5) power_mono zero_le_divide_1_iff zero_le_numeral) - also have "\ = exp x" unfolding num_eq fl_eq exp_real_of_nat_mult[symmetric] using \real (floor_fl x) \ 0\ by auto + also have "\ = exp x" + unfolding num_eq fl_eq exp_real_of_nat_mult[symmetric] + using \real (floor_fl x) \ 0\ by auto finally show ?thesis - unfolding lb_exp.simps if_not_P[OF \\ 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 \\ 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 @@ -1746,21 +2017,28 @@ proof - show ?thesis proof (cases "0 < x") - case False hence "x \ 0" by auto + case False + hence "x \ 0" by auto from exp_boundaries'[OF this] show ?thesis . next - case True hence "-x \ 0" by auto + case True + hence "-x \ 0" by auto have "lb_exp prec x \ exp x" proof - from exp_boundaries'[OF \-x \ 0\] - have ub_exp: "exp (- real x) \ ub_exp prec (-x)" unfolding atLeastAtMost_iff minus_float.rep_eq by auto - - have "float_divl prec 1 (ub_exp prec (-x)) \ 1 / ub_exp prec (-x)" using float_divl[where x=1] by auto + have ub_exp: "exp (- real x) \ ub_exp prec (-x)" + unfolding atLeastAtMost_iff minus_float.rep_eq by auto + + have "float_divl prec 1 (ub_exp prec (-x)) \ 1 / ub_exp prec (-x)" + using float_divl[where x=1] by auto also have "\ \ exp x" - using ub_exp[unfolded inverse_le_iff_le[OF order_less_le_trans[OF exp_gt_zero ub_exp] exp_gt_zero, symmetric]] - unfolding exp_minus nonzero_inverse_inverse_eq[OF exp_not_eq_zero] inverse_eq_divide by auto - finally show ?thesis unfolding lb_exp.simps if_P[OF True] . + using ub_exp[unfolded inverse_le_iff_le[OF order_less_le_trans[OF exp_gt_zero ub_exp] + exp_gt_zero, symmetric]] + unfolding exp_minus nonzero_inverse_inverse_eq[OF exp_not_eq_zero] inverse_eq_divide + by auto + finally show ?thesis + unfolding lb_exp.simps if_P[OF True] . qed moreover have "exp x \ ub_exp prec x" @@ -1768,35 +2046,48 @@ have "\ 0 < -x" using \0 < x\ by auto from exp_boundaries'[OF \-x \ 0\] - have lb_exp: "lb_exp prec (-x) \ exp (- real x)" unfolding atLeastAtMost_iff minus_float.rep_eq by auto + have lb_exp: "lb_exp prec (-x) \ exp (- real x)" + unfolding atLeastAtMost_iff minus_float.rep_eq by auto have "exp x \ (1 :: float) / lb_exp prec (-x)" using lb_exp lb_exp_pos[OF \\ 0 < -x\, of prec] by (simp del: lb_exp.simps add: exp_minus inverse_eq_divide field_simps) - also have "\ \ float_divr prec 1 (lb_exp prec (-x))" using float_divr . - finally show ?thesis unfolding ub_exp.simps if_P[OF True] . + also have "\ \ float_divr prec 1 (lb_exp prec (-x))" + using float_divr . + finally show ?thesis + unfolding ub_exp.simps if_P[OF True] . qed - ultimately show ?thesis by auto + ultimately show ?thesis + by auto qed qed -lemma bnds_exp: "\ (x::real) lx ux. (l, u) = (lb_exp prec lx, ub_exp prec ux) \ x \ {lx .. ux} \ l \ exp x \ exp x \ u" +lemma bnds_exp: "\(x::real) lx ux. (l, u) = + (lb_exp prec lx, ub_exp prec ux) \ x \ {lx .. ux} \ l \ exp x \ exp x \ u" proof (rule allI, rule allI, rule allI, rule impI) - fix x::real and lx ux + fix x :: real and lx ux assume "(l, u) = (lb_exp prec lx, ub_exp prec ux) \ x \ {lx .. ux}" - hence l: "lb_exp prec lx = l " and u: "ub_exp prec ux = u" and x: "x \ {lx .. ux}" by auto - - { from exp_boundaries[of lx prec, unfolded l] - have "l \ exp lx" by (auto simp del: lb_exp.simps) - also have "\ \ exp x" using x by auto - finally have "l \ exp x" . - } moreover - { have "exp x \ exp ux" using x by auto - also have "\ \ u" using exp_boundaries[of ux prec, unfolded u] by (auto simp del: ub_exp.simps) - finally have "exp x \ u" . - } ultimately show "l \ exp x \ exp x \ u" .. + hence l: "lb_exp prec lx = l " and u: "ub_exp prec ux = u" and x: "x \ {lx .. ux}" + by auto + show "l \ exp x \ exp x \ u" + proof + show "l \ exp x" + proof - + from exp_boundaries[of lx prec, unfolded l] + have "l \ exp lx" by (auto simp del: lb_exp.simps) + also have "\ \ exp x" using x by auto + finally show ?thesis . + qed + show "exp x \ u" + proof - + have "exp x \ exp ux" using x by auto + also have "\ \ u" using exp_boundaries[of ux prec, unfolded u] by (auto simp del: ub_exp.simps) + finally show ?thesis . + qed + qed qed + section "Logarithm" subsection "Compute the logarithm series" @@ -1811,7 +2102,8 @@ (lapprox_rat prec 1 (int i)) (- float_round_up prec (x * ub_ln_horner prec n (Suc i) x))" lemma ln_bounds: - assumes "0 \ x" and "x < 1" + assumes "0 \ x" + and "x < 1" shows "(\i=0..<2*n. (- 1) ^ i * (1 / real (i + 1)) * x ^ (Suc i)) \ ln (x + 1)" (is "?lb") and "ln (x + 1) \ (\i=0..<2*n + 1. (- 1) ^ i * (1 / real (i + 1)) * x ^ (Suc i))" (is "?ub") proof - @@ -1823,54 +2115,72 @@ have "norm x < 1" using assms by auto have "?a ----> 0" unfolding Suc_eq_plus1[symmetric] inverse_eq_divide[symmetric] using tendsto_mult[OF LIMSEQ_inverse_real_of_nat LIMSEQ_Suc[OF LIMSEQ_power_zero[OF \norm x < 1\]]] by auto - { fix n have "0 \ ?a n" by (rule mult_nonneg_nonneg, auto simp: \0 \ x\) } - { fix n have "?a (Suc n) \ ?a n" unfolding inverse_eq_divide[symmetric] - proof (rule mult_mono) - show "0 \ x ^ Suc (Suc n)" by (auto simp add: \0 \ x\) - have "x ^ Suc (Suc n) \ x ^ Suc n * 1" unfolding power_Suc2 mult.assoc[symmetric] - by (rule mult_left_mono, fact less_imp_le[OF \x < 1\], auto simp: \0 \ x\) - thus "x ^ Suc (Suc n) \ x ^ Suc n" by auto - qed auto } + have "0 \ ?a n" for n + by (rule mult_nonneg_nonneg) (auto simp: \0 \ x\) + have "?a (Suc n) \ ?a n" for n + unfolding inverse_eq_divide[symmetric] + proof (rule mult_mono) + show "0 \ x ^ Suc (Suc n)" + by (auto simp add: \0 \ x\) + have "x ^ Suc (Suc n) \ x ^ Suc n * 1" + unfolding power_Suc2 mult.assoc[symmetric] + by (rule mult_left_mono, fact less_imp_le[OF \x < 1\]) (auto simp: \0 \ x\) + thus "x ^ Suc (Suc n) \ x ^ Suc n" by auto + qed auto from summable_Leibniz'(2,4)[OF \?a ----> 0\ \\n. 0 \ ?a n\, OF \\n. ?a (Suc n) \ ?a n\, unfolded ln_eq] - show "?lb" and "?ub" unfolding atLeast0LessThan by auto + show ?lb and ?ub + unfolding atLeast0LessThan by auto qed lemma ln_float_bounds: - assumes "0 \ real x" and "real x < 1" + assumes "0 \ real x" + and "real x < 1" shows "x * lb_ln_horner prec (get_even n) 1 x \ ln (x + 1)" (is "?lb \ ?ln") - and "ln (x + 1) \ x * ub_ln_horner prec (get_odd n) 1 x" (is "?ln \ ?ub") + and "ln (x + 1) \ x * ub_ln_horner prec (get_odd n) 1 x" (is "?ln \ ?ub") proof - obtain ev where ev: "get_even n = 2 * ev" using get_even_double .. obtain od where od: "get_odd n = 2 * od + 1" using get_odd_double .. let "?s n" = "(- 1) ^ n * (1 / real (1 + n)) * (real x)^(Suc n)" - have "?lb \ setsum ?s {0 ..< 2 * ev}" unfolding power_Suc2 mult.assoc[symmetric] times_float.rep_eq setsum_left_distrib[symmetric] unfolding mult.commute[of "real x"] ev + have "?lb \ setsum ?s {0 ..< 2 * ev}" + unfolding power_Suc2 mult.assoc[symmetric] times_float.rep_eq setsum_left_distrib[symmetric] + unfolding mult.commute[of "real x"] ev using horner_bounds(1)[where G="\ i k. Suc k" and F="\x. x" and f="\x. x" and lb="\n i k x. lb_ln_horner prec n k x" and ub="\n i k x. ub_ln_horner prec n k x" and j'=1 and n="2*ev", OF \0 \ real x\ refl lb_ln_horner.simps ub_ln_horner.simps] \0 \ real x\ by (rule mult_right_mono) - also have "\ \ ?ln" using ln_bounds(1)[OF \0 \ real x\ \real x < 1\] by auto + also have "\ \ ?ln" + using ln_bounds(1)[OF \0 \ real x\ \real x < 1\] by auto finally show "?lb \ ?ln" . - have "?ln \ setsum ?s {0 ..< 2 * od + 1}" using ln_bounds(2)[OF \0 \ real x\ \real x < 1\] by auto - also have "\ \ ?ub" unfolding power_Suc2 mult.assoc[symmetric] times_float.rep_eq setsum_left_distrib[symmetric] unfolding mult.commute[of "real x"] od + have "?ln \ setsum ?s {0 ..< 2 * od + 1}" + using ln_bounds(2)[OF \0 \ real x\ \real x < 1\] by auto + also have "\ \ ?ub" + unfolding power_Suc2 mult.assoc[symmetric] times_float.rep_eq setsum_left_distrib[symmetric] + unfolding mult.commute[of "real x"] od using horner_bounds(2)[where G="\ i k. Suc k" and F="\x. x" and f="\x. x" and lb="\n i k x. lb_ln_horner prec n k x" and ub="\n i k x. ub_ln_horner prec n k x" and j'=1 and n="2*od+1", OF \0 \ real x\ refl lb_ln_horner.simps ub_ln_horner.simps] \0 \ real x\ by (rule mult_right_mono) finally show "?ln \ ?ub" . qed -lemma ln_add: - fixes x::real assumes "0 < x" and "0 < y" shows "ln (x + y) = ln x + ln (1 + y / x)" +lemma ln_add: + fixes x :: real + assumes "0 < x" and "0 < y" + shows "ln (x + y) = ln x + ln (1 + y / x)" proof - have "x \ 0" using assms by auto - have "x + y = x * (1 + y / x)" unfolding distrib_left times_divide_eq_right nonzero_mult_divide_cancel_left[OF \x \ 0\] by auto + have "x + y = x * (1 + y / x)" + unfolding distrib_left times_divide_eq_right nonzero_mult_divide_cancel_left[OF \x \ 0\] + by auto moreover have "0 < y / x" using assms by auto hence "0 < 1 + y / x" by auto - ultimately show ?thesis using ln_mult assms by auto + ultimately show ?thesis + using ln_mult assms by auto qed + subsection "Compute the logarithm of 2" definition ub_ln2 where "ub_ln2 prec = (let third = rapprox_rat (max prec 1) 1 3 @@ -1896,7 +2206,8 @@ have ub3: "1 / 3 \ ?uthird" using rapprox_rat[of 1 3] by auto hence ub3_lb: "0 \ real ?uthird" by auto - have lb2: "0 \ real (Float 1 (- 1))" and ub2: "real (Float 1 (- 1)) < 1" unfolding Float_num by auto + have lb2: "0 \ real (Float 1 (- 1))" and ub2: "real (Float 1 (- 1)) < 1" + unfolding Float_num by auto have "0 \ (1::int)" and "0 < (3::int)" by auto have ub3_ub: "real ?uthird < 1" @@ -1906,24 +2217,31 @@ 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 ln2_sum Float_num(4)[symmetric] + 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) \ ln (real ?uthird + 1)" unfolding ln_le_cancel_iff[OF third_gt0 uthird_gt0] using ub3 by auto + have "ln (1 / 3 + 1) \ ln (real ?uthird + 1)" + unfolding ln_le_cancel_iff[OF third_gt0 uthird_gt0] using ub3 by auto also have "\ \ ?uthird * ub_ln_horner prec (get_odd prec) 1 ?uthird" using ln_float_bounds(2)[OF ub3_lb ub3_ub] . also note float_round_up finally show "ln (1 / 3 + 1) \ float_round_up prec (?uthird * ub_ln_horner prec (get_odd prec) 1 ?uthird)" . qed - show ?lb_ln2 unfolding lb_ln2_def Let_def ln2_sum Float_num(4)[symmetric] + 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 \ ln (real ?lthird + 1)" using ln_float_bounds(1)[OF lb3_lb lb3_ub] . note float_round_down_le[OF this] - also have "\ \ ln (1 / 3 + 1)" unfolding ln_le_cancel_iff[OF lthird_gt0 third_gt0] using lb3 by auto - finally show "float_round_down prec (?lthird * lb_ln_horner prec (get_even prec) 1 ?lthird) \ ln (1 / 3 + 1)" . + also have "\ \ ln (1 / 3 + 1)" + unfolding ln_le_cancel_iff[OF lthird_gt0 third_gt0] + using lb3 by auto + finally show "float_round_down prec (?lthird * lb_ln_horner prec (get_even prec) 1 ?lthird) \ + ln (1 / 3 + 1)" . qed qed + subsection "Compute the logarithm in the entire domain" function ub_ln :: "nat \ float \ float option" and lb_ln :: "nat \ float \ float option" where @@ -1942,26 +2260,32 @@ horner (max (x * lapprox_rat prec 2 3 - 1) 0))) else let l = bitlen (mantissa x) - 1 in 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 (\ v. let (prec, x) = case_sum id id v in (if x < 1 then 1 else 0))", auto) - fix prec and x :: float assume "\ real x \ 0" and "real x < 1" and "real (float_divl (max prec (Suc 0)) 1 x) < 1" - hence "0 < real x" "1 \ max prec (Suc 0)" "real x < 1" by auto + by pat_completeness auto + +termination +proof (relation "measure (\ v. let (prec, x) = case_sum id id v in (if x < 1 then 1 else 0))", auto) + fix prec and x :: float + assume "\ real x \ 0" and "real x < 1" and "real (float_divl (max prec (Suc 0)) 1 x) < 1" + hence "0 < real x" "1 \ max prec (Suc 0)" "real x < 1" + by auto from float_divl_pos_less1_bound[OF \0 < real x\ \real x < 1\[THEN less_imp_le] \1 \ max prec (Suc 0)\] - show False using \real (float_divl (max prec (Suc 0)) 1 x) < 1\ by auto + show False + using \real (float_divl (max prec (Suc 0)) 1 x) < 1\ by auto next - fix prec x assume "\ real x \ 0" and "real x < 1" and "real (float_divr prec 1 x) < 1" + fix prec x + assume "\ real x \ 0" and "real x < 1" and "real (float_divr prec 1 x) < 1" hence "0 < x" by auto - from float_divr_pos_less1_lower_bound[OF \0 < x\, of prec] \real x < 1\ - show False using \real (float_divr prec 1 x) < 1\ by auto + from float_divr_pos_less1_lower_bound[OF \0 < x\, of prec] \real x < 1\ show False + using \real (float_divr prec 1 x) < 1\ by auto qed -lemma float_pos_eq_mantissa_pos: "x > 0 \ mantissa x > 0" +lemma float_pos_eq_mantissa_pos: "x > 0 \ mantissa x > 0" apply (subst Float_mantissa_exponent[of x, symmetric]) apply (auto simp add: zero_less_mult_iff zero_float_def dest: less_zeroE) - by (metis not_le powr_ge_pzero) - -lemma Float_pos_eq_mantissa_pos: "Float m e > 0 \ m > 0" + apply (metis not_le powr_ge_pzero) + done + +lemma Float_pos_eq_mantissa_pos: "Float m e > 0 \ m > 0" using powr_gt_zero[of 2 "e"] by (auto simp add: zero_less_mult_iff zero_float_def simp del: powr_gt_zero dest: less_zeroE) @@ -1974,7 +2298,9 @@ proof - from assms have mantissa_pos: "m > 0" "mantissa x > 0" using Float_pos_eq_mantissa_pos[of m e] float_pos_eq_mantissa_pos[of x] by simp_all - thus ?th1 using bitlen_Float[of m e] assms by (auto simp add: zero_less_mult_iff intro!: arg_cong2[where f=Float]) + thus ?th1 + using bitlen_Float[of m e] assms + by (auto simp add: zero_less_mult_iff intro!: arg_cong2[where f=Float]) have "x \ float_of 0" unfolding zero_float_def[symmetric] using \0 < x\ by auto from denormalize_shift[OF assms(1) this] guess i . note i = this @@ -2010,16 +2336,22 @@ 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 \ m > 0" by simp - thus ?th1 ?th2 using Float_representation_aux[of m e] unfolding x_def[symmetric] + from assms Float_pos_eq_mantissa_pos have "x > 0 \ m > 0" + by simp + thus ?th1 ?th2 + using Float_representation_aux[of m e] + unfolding x_def[symmetric] by (auto dest: not_leE) qed -lemma ln_shifted_float: assumes "0 < m" shows "ln (Float m e) = ln 2 * (e + (bitlen m - 1)) + ln (Float m (- (bitlen m - 1)))" +lemma ln_shifted_float: + assumes "0 < m" + shows "ln (Float m e) = ln 2 * (e + (bitlen m - 1)) + ln (Float m (- (bitlen m - 1)))" proof - let ?B = "2^nat (bitlen m - 1)" def bl \ "bitlen m - 1" - have "0 < real m" and "\X. (0 :: real) < 2^X" and "0 < (2 :: real)" and "m \ 0" using assms by auto + have "0 < real m" and "\X. (0 :: real) < 2^X" and "0 < (2 :: real)" and "m \ 0" + using assms by auto hence "0 \ bl" by (simp add: bitlen_def bl_def) show ?thesis proof (cases "0 \ e") @@ -2032,18 +2364,25 @@ apply (cases "bl = 0", simp_all add: powr_minus ln_inverse ln_powr field_simps) done next - case False hence "0 < -e" by auto - have lne: "ln (2 powr real e) = ln (inverse (2 powr - e))" by (simp add: powr_minus) - hence pow_gt0: "(0::real) < 2^nat (-e)" by auto - hence inv_gt0: "(0::real) < inverse (2^nat (-e))" by auto - show ?thesis using False unfolding bl_def[symmetric] using \0 < real m\ \0 \ bl\ + case False + hence "0 < -e" by auto + have lne: "ln (2 powr real e) = ln (inverse (2 powr - e))" + by (simp add: powr_minus) + hence pow_gt0: "(0::real) < 2^nat (-e)" + by auto + hence inv_gt0: "(0::real) < inverse (2^nat (-e))" + by auto + show ?thesis + using False unfolding bl_def[symmetric] + using \0 < real m\ \0 \ bl\ by (auto simp add: lne ln_mult ln_powr ln_div field_simps) qed qed -lemma ub_ln_lb_ln_bounds': assumes "1 \ x" +lemma ub_ln_lb_ln_bounds': + assumes "1 \ x" shows "the (lb_ln prec x) \ ln x \ ln x \ the (ub_ln prec x)" - (is "?lb \ ?ln \ ?ln \ ?ub") + (is "?lb \ ?ln \ ?ln \ ?ub") proof (cases "x < Float 1 1") case True hence "real (x - 1) < 1" and "real x < 2" by auto @@ -2055,11 +2394,15 @@ show ?thesis proof (cases "x \ Float 3 (- 1)") case True - show ?thesis unfolding lb_ln.simps unfolding ub_ln.simps Let_def - using ln_float_bounds[OF \0 \ real (x - 1)\ \real (x - 1) < 1\, of prec] \\ x \ 0\ \\ x < 1\ True + show ?thesis + unfolding lb_ln.simps + unfolding ub_ln.simps Let_def + using ln_float_bounds[OF \0 \ real (x - 1)\ \real (x - 1) < 1\, of prec] + \\ x \ 0\ \\ x < 1\ True by (auto intro!: float_round_down_le float_round_up_le) next - case False hence *: "3 / 2 < x" by auto + case False + hence *: "3 / 2 < x" by auto with ln_add[of "3 / 2" "x - 3 / 2"] have add: "ln x = ln (3 / 2) + ln (real x * 2 / 3)" @@ -2137,15 +2480,20 @@ proof - def m \ "mantissa x" def e \ "exponent x" - from Float_mantissa_exponent[of x] have Float: "x = Float m e" by (simp add: m_def e_def) + from Float_mantissa_exponent[of x] have Float: "x = Float m e" + by (simp add: m_def e_def) let ?s = "Float (e + (bitlen m - 1)) 0" let ?x = "Float m (- (bitlen m - 1))" - have "0 < m" and "m \ 0" using \0 < x\ Float powr_gt_zero[of 2 e] + have "0 < m" and "m \ 0" using \0 < x\ Float powr_gt_zero[of 2 e] apply (auto simp add: zero_less_mult_iff) - using not_le powr_ge_pzero by blast - def bl \ "bitlen m - 1" hence "bl \ 0" using \m > 0\ by (simp add: bitlen_def) - have "1 \ Float m e" using \1 \ x\ Float unfolding less_eq_float_def by auto + using not_le powr_ge_pzero apply blast + done + def bl \ "bitlen m - 1" + hence "bl \ 0" + using \m > 0\ by (simp add: bitlen_def) + have "1 \ Float m e" + using \1 \ x\ Float unfolding less_eq_float_def by auto from bitlen_div[OF \0 < m\] float_gt1_scale[OF \1 \ Float m e\] \bl \ 0\ have x_bnds: "0 \ real (?x - 1)" "real (?x - 1) < 1" unfolding bl_def[symmetric] @@ -2153,7 +2501,8 @@ (auto simp : powr_minus field_simps inverse_eq_divide) { - have "float_round_down prec (lb_ln2 prec * ?s) \ ln 2 * (e + (bitlen m - 1))" (is "real ?lb2 \ _") + have "float_round_down prec (lb_ln2 prec * ?s) \ ln 2 * (e + (bitlen m - 1))" + (is "real ?lb2 \ _") apply (rule float_round_down_le) unfolding nat_0 power_0 mult_1_right times_float.rep_eq using lb_ln2[of prec] @@ -2171,17 +2520,18 @@ moreover { from ln_float_bounds(2)[OF x_bnds] - have "ln ?x \ float_round_up prec ((?x - 1) * ub_ln_horner prec (get_odd prec) 1 (?x - 1))" (is "_ \ real ?ub_horner") + have "ln ?x \ float_round_up prec ((?x - 1) * ub_ln_horner prec (get_odd prec) 1 (?x - 1))" + (is "_ \ real ?ub_horner") by (auto intro!: float_round_up_le) moreover - have "ln 2 * (e + (bitlen m - 1)) \ float_round_up prec (ub_ln2 prec * ?s)" (is "_ \ real ?ub2") + have "ln 2 * (e + (bitlen m - 1)) \ float_round_up prec (ub_ln2 prec * ?s)" + (is "_ \ 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) from float_gt1_scale[OF \1 \ Float m e\] show "0 \ real (e + (bitlen m - 1))" by auto - next have "0 \ ln (2 :: real)" by simp thus "0 \ real (ub_ln2 prec)" using ub_ln2[of prec] by arith qed auto @@ -2189,23 +2539,33 @@ 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 \\ x \ 0\] if_not_P[OF \\ x < 1\] if_not_P[OF False] if_not_P[OF \\ x \ Float 3 (- 1)\] Let_def - unfolding plus_float.rep_eq e_def[symmetric] m_def[symmetric] by simp + ultimately show ?thesis + unfolding lb_ln.simps + unfolding ub_ln.simps + unfolding if_not_P[OF \\ x \ 0\] if_not_P[OF \\ x < 1\] + if_not_P[OF False] if_not_P[OF \\ x \ Float 3 (- 1)\] Let_def + unfolding plus_float.rep_eq e_def[symmetric] m_def[symmetric] + by simp qed qed lemma ub_ln_lb_ln_bounds: assumes "0 < x" shows "the (lb_ln prec x) \ ln x \ ln x \ the (ub_ln prec x)" - (is "?lb \ ?ln \ ?ln \ ?ub") + (is "?lb \ ?ln \ ?ln \ ?ub") proof (cases "x < 1") - case False hence "1 \ x" unfolding less_float_def less_eq_float_def by auto - show ?thesis using ub_ln_lb_ln_bounds'[OF \1 \ x\] . + case False + hence "1 \ x" + unfolding less_float_def less_eq_float_def by auto + show ?thesis + using ub_ln_lb_ln_bounds'[OF \1 \ x\] . next - case True have "\ x \ 0" using \0 < x\ by auto - from True have "real x \ 1" "x \ 1" by simp_all - have "0 < real x" and "real x \ 0" using \0 < x\ by auto + case True + have "\ x \ 0" using \0 < x\ by auto + from True have "real x \ 1" "x \ 1" + by simp_all + have "0 < real x" and "real x \ 0" + using \0 < x\ by auto hence A: "0 < 1 / real x" by auto { @@ -2238,12 +2598,17 @@ proof - have "0 < x" proof (rule ccontr) - assume "\ 0 < x" hence "x \ 0" unfolding less_eq_float_def less_float_def by auto - thus False using assms by auto + assume "\ 0 < x" + hence "x \ 0" + unfolding less_eq_float_def less_float_def by auto + thus False + using assms by auto qed thus "0 < real x" by auto - have "the (lb_ln prec x) \ ln x" using ub_ln_lb_ln_bounds[OF \0 < x\] .. - thus "y \ ln x" unfolding assms[symmetric] by auto + have "the (lb_ln prec x) \ ln x" + using ub_ln_lb_ln_bounds[OF \0 < x\] .. + thus "y \ ln x" + unfolding assms[symmetric] by auto qed lemma ub_ln: @@ -2252,31 +2617,43 @@ proof - have "0 < x" proof (rule ccontr) - assume "\ 0 < x" hence "x \ 0" by auto - thus False using assms by auto + assume "\ 0 < x" + hence "x \ 0" by auto + thus False + using assms by auto qed thus "0 < real x" by auto - have "ln x \ the (ub_ln prec x)" using ub_ln_lb_ln_bounds[OF \0 < x\] .. - thus "ln x \ y" unfolding assms[symmetric] by auto + have "ln x \ the (ub_ln prec x)" + using ub_ln_lb_ln_bounds[OF \0 < x\] .. + thus "ln x \ y" + unfolding assms[symmetric] by auto qed -lemma bnds_ln: "\ (x::real) lx ux. (Some l, Some u) = (lb_ln prec lx, ub_ln prec ux) \ x \ {lx .. ux} \ l \ ln x \ ln x \ u" +lemma bnds_ln: "\(x::real) lx ux. (Some l, Some u) = + (lb_ln prec lx, ub_ln prec ux) \ x \ {lx .. ux} \ l \ ln x \ ln x \ u" proof (rule allI, rule allI, rule allI, rule impI) - fix x::real and lx ux + fix x :: real + fix lx ux assume "(Some l, Some u) = (lb_ln prec lx, ub_ln prec ux) \ x \ {lx .. ux}" - hence l: "Some l = lb_ln prec lx " and u: "Some u = ub_ln prec ux" and x: "x \ {lx .. ux}" by auto - - have "ln ux \ u" and "0 < real ux" using ub_ln u by auto - have "l \ ln lx" and "0 < real lx" and "0 < x" using lb_ln[OF l] x by auto + hence l: "Some l = lb_ln prec lx " and u: "Some u = ub_ln prec ux" and x: "x \ {lx .. ux}" + by auto + + have "ln ux \ u" and "0 < real ux" + using ub_ln u by auto + have "l \ ln lx" and "0 < real lx" and "0 < x" + using lb_ln[OF l] x by auto from ln_le_cancel_iff[OF \0 < real lx\ \0 < x\] \l \ ln lx\ - have "l \ ln x" using x unfolding atLeastAtMost_iff by auto + have "l \ ln x" + using x unfolding atLeastAtMost_iff by auto moreover from ln_le_cancel_iff[OF \0 < x\ \0 < real ux\] \ln ux \ real u\ - have "ln x \ u" using x unfolding atLeastAtMost_iff by auto + have "ln x \ u" + using x unfolding atLeastAtMost_iff by auto ultimately show "l \ ln x \ ln x \ u" .. qed + section "Implement floatarith" subsection "Define syntax and semantics" @@ -2317,35 +2694,42 @@ "interpret_floatarith (Num f) vs = f" | "interpret_floatarith (Var n) vs = vs ! n" -lemma interpret_floatarith_divide: "interpret_floatarith (Mult a (Inverse b)) vs = (interpret_floatarith a vs) / (interpret_floatarith b vs)" +lemma interpret_floatarith_divide: + "interpret_floatarith (Mult a (Inverse b)) vs = + (interpret_floatarith a vs) / (interpret_floatarith b vs)" unfolding divide_inverse interpret_floatarith.simps .. -lemma interpret_floatarith_diff: "interpret_floatarith (Add a (Minus b)) vs = (interpret_floatarith a vs) - (interpret_floatarith b vs)" +lemma interpret_floatarith_diff: + "interpret_floatarith (Add a (Minus b)) vs = + (interpret_floatarith a vs) - (interpret_floatarith b vs)" unfolding interpret_floatarith.simps by simp -lemma interpret_floatarith_sin: "interpret_floatarith (Cos (Add (Mult Pi (Num (Float 1 (- 1)))) (Minus a))) vs = - sin (interpret_floatarith a vs)" +lemma interpret_floatarith_sin: + "interpret_floatarith (Cos (Add (Mult Pi (Num (Float 1 (- 1)))) (Minus a))) vs = + sin (interpret_floatarith a vs)" unfolding sin_cos_eq interpret_floatarith.simps - interpret_floatarith_divide interpret_floatarith_diff + interpret_floatarith_divide interpret_floatarith_diff by auto lemma interpret_floatarith_tan: "interpret_floatarith (Mult (Cos (Add (Mult Pi (Num (Float 1 (- 1)))) (Minus a))) (Inverse (Cos a))) vs = - tan (interpret_floatarith a vs)" + tan (interpret_floatarith a vs)" unfolding interpret_floatarith.simps(3,4) interpret_floatarith_sin tan_def divide_inverse by auto -lemma interpret_floatarith_log: - "interpret_floatarith ((Mult (Ln x) (Inverse (Ln b)))) vs = - log (interpret_floatarith b vs) (interpret_floatarith x vs)" +lemma interpret_floatarith_log: + "interpret_floatarith ((Mult (Ln x) (Inverse (Ln b)))) vs = + log (interpret_floatarith b vs) (interpret_floatarith x vs)" unfolding log_def interpret_floatarith.simps divide_inverse .. lemma interpret_floatarith_num: shows "interpret_floatarith (Num (Float 0 0)) vs = 0" - and "interpret_floatarith (Num (Float 1 0)) vs = 1" - and "interpret_floatarith (Num (Float (- 1) 0)) vs = - 1" - and "interpret_floatarith (Num (Float (numeral a) 0)) vs = numeral a" - and "interpret_floatarith (Num (Float (- numeral a) 0)) vs = - numeral a" by auto + and "interpret_floatarith (Num (Float 1 0)) vs = 1" + and "interpret_floatarith (Num (Float (- 1) 0)) vs = - 1" + and "interpret_floatarith (Num (Float (numeral a) 0)) vs = numeral a" + and "interpret_floatarith (Num (Float (- numeral a) 0)) vs = - numeral a" + by auto + subsection "Implement approximation function" @@ -2362,8 +2746,7 @@ "lift_un' (Some (l1, u1)) f = Some (f l1 u1)" | "lift_un' b f = None" -definition -"bounded_by xs vs \ +definition "bounded_by xs vs \ (\ i < length vs. case vs ! i of None \ True | Some (l, u) \ xs ! i \ { real l .. real u })" @@ -2375,29 +2758,31 @@ lemma bounded_by_update: assumes "bounded_by xs vs" - and bnd: "xs ! i \ { real l .. real u }" + and bnd: "xs ! i \ { real l .. real u }" shows "bounded_by xs (vs[i := Some (l,u)])" proof - -{ fix j - let ?vs = "vs[i := Some (l,u)]" - assume "j < length ?vs" hence [simp]: "j < length vs" by simp - have "case ?vs ! j of None \ True | Some (l, u) \ xs ! j \ { real l .. real u }" - proof (cases "?vs ! j") - case (Some b) - thus ?thesis - proof (cases "i = j") - case True - thus ?thesis using \?vs ! j = Some b\ and bnd by auto - next - case False - thus ?thesis using \bounded_by xs vs\ unfolding bounded_by_def by auto - qed - qed auto } + { + fix j + let ?vs = "vs[i := Some (l,u)]" + assume "j < length ?vs" + hence [simp]: "j < length vs" by simp + have "case ?vs ! j of None \ True | Some (l, u) \ xs ! j \ { real l .. real u }" + proof (cases "?vs ! j") + case (Some b) + thus ?thesis + proof (cases "i = j") + case True + thus ?thesis using \?vs ! j = Some b\ and bnd by auto + next + case False + thus ?thesis using \bounded_by xs vs\ unfolding bounded_by_def by auto + qed + qed auto + } thus ?thesis unfolding bounded_by_def by auto qed -lemma bounded_by_None: - shows "bounded_by xs (replicate (length xs) None)" +lemma bounded_by_None: "bounded_by xs (replicate (length xs) None)" unfolding bounded_by_def by auto fun approx approx' :: "nat \ floatarith \ (float * float) option list \ (float * float) option" where @@ -2433,41 +2818,56 @@ assumes lift_bin'_Some: "Some (l, u) = lift_bin' a b f" shows "\ l1 u1 l2 u2. Some (l1, u1) = a \ Some (l2, u2) = b" proof (cases a) - case None hence "None = lift_bin' a b f" unfolding None lift_bin'.simps .. - thus ?thesis using lift_bin'_Some by auto + case None + hence "None = lift_bin' a b f" + unfolding None lift_bin'.simps .. + thus ?thesis + using lift_bin'_Some by auto next case (Some a') show ?thesis proof (cases b) - case None hence "None = lift_bin' a b f" unfolding None lift_bin'.simps .. + case None + hence "None = lift_bin' a b f" + unfolding None lift_bin'.simps .. thus ?thesis using lift_bin'_Some by auto next case (Some b') - obtain la ua where a': "a' = (la, ua)" by (cases a', auto) - obtain lb ub where b': "b' = (lb, ub)" by (cases b', auto) - thus ?thesis unfolding \a = Some a'\ \b = Some b'\ a' b' by auto + obtain la ua where a': "a' = (la, ua)" + by (cases a') auto + obtain lb ub where b': "b' = (lb, ub)" + by (cases b') auto + thus ?thesis + unfolding \a = Some a'\ \b = Some b'\ a' b' by auto qed qed lemma lift_bin'_f: assumes lift_bin'_Some: "Some (l, u) = lift_bin' (g a) (g b) f" - and Pa: "\l u. Some (l, u) = g a \ P l u a" and Pb: "\l u. Some (l, u) = g b \ P l u b" + and Pa: "\l u. Some (l, u) = g a \ P l u a" + and Pb: "\l u. Some (l, u) = g b \ P l u b" shows "\ l1 u1 l2 u2. P l1 u1 a \ P l2 u2 b \ l = fst (f l1 u1 l2 u2) \ u = snd (f l1 u1 l2 u2)" proof - obtain l1 u1 l2 u2 - where Sa: "Some (l1, u1) = g a" and Sb: "Some (l2, u2) = g b" using lift_bin'_ex[OF assms(1)] by auto - have lu: "(l, u) = f l1 u1 l2 u2" using lift_bin'_Some[unfolded Sa[symmetric] Sb[symmetric] lift_bin'.simps] by auto - have "l = fst (f l1 u1 l2 u2)" and "u = snd (f l1 u1 l2 u2)" unfolding lu[symmetric] by auto - thus ?thesis using Pa[OF Sa] Pb[OF Sb] by auto + where Sa: "Some (l1, u1) = g a" + and Sb: "Some (l2, u2) = g b" + using lift_bin'_ex[OF assms(1)] by auto + have lu: "(l, u) = f l1 u1 l2 u2" + using lift_bin'_Some[unfolded Sa[symmetric] Sb[symmetric] lift_bin'.simps] by auto + have "l = fst (f l1 u1 l2 u2)" and "u = snd (f l1 u1 l2 u2)" + unfolding lu[symmetric] by auto + thus ?thesis + using Pa[OF Sa] Pb[OF Sb] by auto qed lemma approx_approx': - assumes Pa: "\l u. Some (l, u) = approx prec a vs \ l \ interpret_floatarith a xs \ interpret_floatarith a xs \ u" - and approx': "Some (l, u) = approx' prec a vs" + assumes Pa: "\l u. Some (l, u) = approx prec a vs \ + l \ interpret_floatarith a xs \ interpret_floatarith a xs \ u" + and approx': "Some (l, u) = approx' prec a vs" shows "l \ interpret_floatarith a xs \ interpret_floatarith a xs \ u" proof - obtain l' u' where S: "Some (l', u') = approx prec a vs" - using approx' unfolding approx'.simps by (cases "approx prec a vs", auto) + using approx' unfolding approx'.simps by (cases "approx prec a vs") auto have l': "l = float_round_down prec l'" and u': "u = float_round_up prec u'" using approx' unfolding approx'.simps S[symmetric] by auto show ?thesis unfolding l' u' @@ -2477,11 +2877,13 @@ lemma lift_bin': assumes lift_bin'_Some: "Some (l, u) = lift_bin' (approx' prec a bs) (approx' prec b bs) f" - and Pa: "\l u. Some (l, u) = approx prec a bs \ l \ interpret_floatarith a xs \ interpret_floatarith a xs \ u" (is "\l u. _ = ?g a \ ?P l u a") - and Pb: "\l u. Some (l, u) = approx prec b bs \ l \ interpret_floatarith b xs \ interpret_floatarith b xs \ u" - shows "\ l1 u1 l2 u2. (l1 \ interpret_floatarith a xs \ interpret_floatarith a xs \ u1) \ - (l2 \ interpret_floatarith b xs \ interpret_floatarith b xs \ u2) \ - l = fst (f l1 u1 l2 u2) \ u = snd (f l1 u1 l2 u2)" + and Pa: "\l u. Some (l, u) = approx prec a bs \ + l \ interpret_floatarith a xs \ interpret_floatarith a xs \ u" (is "\l u. _ = ?g a \ ?P l u a") + and Pb: "\l u. Some (l, u) = approx prec b bs \ + l \ interpret_floatarith b xs \ interpret_floatarith b xs \ u" + shows "\l1 u1 l2 u2. (l1 \ interpret_floatarith a xs \ interpret_floatarith a xs \ u1) \ + (l2 \ interpret_floatarith b xs \ interpret_floatarith b xs \ u2) \ + l = fst (f l1 u1 l2 u2) \ u = snd (f l1 u1 l2 u2)" proof - { fix l u assume "Some (l, u) = approx' prec a bs" with approx_approx'[of prec a bs, OF _ this] Pa @@ -2498,68 +2900,93 @@ assumes lift_un'_Some: "Some (l, u) = lift_un' a f" shows "\ l u. Some (l, u) = a" proof (cases a) - case None hence "None = lift_un' a f" unfolding None lift_un'.simps .. - thus ?thesis using lift_un'_Some by auto + case None + hence "None = lift_un' a f" + unfolding None lift_un'.simps .. + thus ?thesis + using lift_un'_Some by auto next case (Some a') - obtain la ua where a': "a' = (la, ua)" by (cases a', auto) - thus ?thesis unfolding \a = Some a'\ a' by auto + obtain la ua where a': "a' = (la, ua)" + by (cases a') auto + thus ?thesis + unfolding \a = Some a'\ a' by auto qed lemma lift_un'_f: assumes lift_un'_Some: "Some (l, u) = lift_un' (g a) f" - and Pa: "\l u. Some (l, u) = g a \ P l u a" + and Pa: "\l u. Some (l, u) = g a \ P l u a" shows "\ l1 u1. P l1 u1 a \ l = fst (f l1 u1) \ u = snd (f l1 u1)" proof - - obtain l1 u1 where Sa: "Some (l1, u1) = g a" using lift_un'_ex[OF assms(1)] by auto - have lu: "(l, u) = f l1 u1" using lift_un'_Some[unfolded Sa[symmetric] lift_un'.simps] by auto - have "l = fst (f l1 u1)" and "u = snd (f l1 u1)" unfolding lu[symmetric] by auto - thus ?thesis using Pa[OF Sa] by auto + obtain l1 u1 where Sa: "Some (l1, u1) = g a" + using lift_un'_ex[OF assms(1)] by auto + have lu: "(l, u) = f l1 u1" + using lift_un'_Some[unfolded Sa[symmetric] lift_un'.simps] by auto + have "l = fst (f l1 u1)" and "u = snd (f l1 u1)" + unfolding lu[symmetric] by auto + thus ?thesis + using Pa[OF Sa] by auto qed lemma lift_un': assumes lift_un'_Some: "Some (l, u) = lift_un' (approx' prec a bs) f" - and Pa: "\l u. Some (l, u) = approx prec a bs \ l \ interpret_floatarith a xs \ interpret_floatarith a xs \ u" (is "\l u. _ = ?g a \ ?P l u a") - shows "\ l1 u1. (l1 \ interpret_floatarith a xs \ interpret_floatarith a xs \ u1) \ - l = fst (f l1 u1) \ u = snd (f l1 u1)" + and Pa: "\l u. Some (l, u) = approx prec a bs \ + l \ interpret_floatarith a xs \ interpret_floatarith a xs \ u" + (is "\l u. _ = ?g a \ ?P l u a") + shows "\l1 u1. (l1 \ interpret_floatarith a xs \ interpret_floatarith a xs \ u1) \ + l = fst (f l1 u1) \ u = snd (f l1 u1)" proof - - { fix l u assume "Some (l, u) = approx' prec a bs" - with approx_approx'[of prec a bs, OF _ this] Pa - have "l \ interpret_floatarith a xs \ interpret_floatarith a xs \ u" by auto } note Pa = this + have Pa: "l \ interpret_floatarith a xs \ interpret_floatarith a xs \ u" + if "Some (l, u) = approx' prec a bs" for l u + using approx_approx'[of prec a bs, OF _ that] Pa + by auto from lift_un'_f[where g="\a. approx' prec a bs" and P = ?P, OF lift_un'_Some, OF Pa] show ?thesis by auto qed lemma lift_un'_bnds: assumes bnds: "\ (x::real) lx ux. (l, u) = f lx ux \ x \ { lx .. ux } \ l \ f' x \ f' x \ u" - and lift_un'_Some: "Some (l, u) = lift_un' (approx' prec a bs) f" - and Pa: "\l u. Some (l, u) = approx prec a bs \ l \ interpret_floatarith a xs \ interpret_floatarith a xs \ u" + and lift_un'_Some: "Some (l, u) = lift_un' (approx' prec a bs) f" + and Pa: "\l u. Some (l, u) = approx prec a bs \ + l \ interpret_floatarith a xs \ interpret_floatarith a xs \ u" shows "real l \ f' (interpret_floatarith a xs) \ f' (interpret_floatarith a xs) \ real u" proof - from lift_un'[OF lift_un'_Some Pa] - obtain l1 u1 where "l1 \ interpret_floatarith a xs" and "interpret_floatarith a xs \ u1" and "l = fst (f l1 u1)" and "u = snd (f l1 u1)" by blast - hence "(l, u) = f l1 u1" and "interpret_floatarith a xs \ {l1 .. u1}" by auto - thus ?thesis using bnds by auto + obtain l1 u1 where "l1 \ interpret_floatarith a xs" + and "interpret_floatarith a xs \ u1" + and "l = fst (f l1 u1)" + and "u = snd (f l1 u1)" + by blast + hence "(l, u) = f l1 u1" and "interpret_floatarith a xs \ {l1 .. u1}" + by auto + thus ?thesis + using bnds by auto qed lemma lift_un_ex: assumes lift_un_Some: "Some (l, u) = lift_un a f" - shows "\ l u. Some (l, u) = a" + shows "\l u. Some (l, u) = a" proof (cases a) - case None hence "None = lift_un a f" unfolding None lift_un.simps .. - thus ?thesis using lift_un_Some by auto + case None + hence "None = lift_un a f" + unfolding None lift_un.simps .. + thus ?thesis + using lift_un_Some by auto next case (Some a') - obtain la ua where a': "a' = (la, ua)" by (cases a', auto) - thus ?thesis unfolding \a = Some a'\ a' by auto + obtain la ua where a': "a' = (la, ua)" + by (cases a') auto + thus ?thesis + unfolding \a = Some a'\ a' by auto qed lemma lift_un_f: assumes lift_un_Some: "Some (l, u) = lift_un (g a) f" - and Pa: "\l u. Some (l, u) = g a \ P l u a" + and Pa: "\l u. Some (l, u) = g a \ P l u a" shows "\ l1 u1. P l1 u1 a \ Some l = fst (f l1 u1) \ Some u = snd (f l1 u1)" proof - - obtain l1 u1 where Sa: "Some (l1, u1) = g a" using lift_un_ex[OF assms(1)] by auto + obtain l1 u1 where Sa: "Some (l1, u1) = g a" + using lift_un_ex[OF assms(1)] by auto have "fst (f l1 u1) \ None \ snd (f l1 u1) \ None" proof (rule ccontr) assume "\ (fst (f l1 u1) \ None \ snd (f l1 u1) \ None)" @@ -2567,64 +2994,88 @@ hence "lift_un (g a) f = None" proof (cases "fst (f l1 u1) = None") case True - then obtain b where b: "f l1 u1 = (None, b)" by (cases "f l1 u1", auto) - thus ?thesis unfolding Sa[symmetric] lift_un.simps b by auto + then obtain b where b: "f l1 u1 = (None, b)" + by (cases "f l1 u1") auto + thus ?thesis + unfolding Sa[symmetric] lift_un.simps b by auto next - case False hence "snd (f l1 u1) = None" using or by auto - with False obtain b where b: "f l1 u1 = (Some b, None)" by (cases "f l1 u1", auto) - thus ?thesis unfolding Sa[symmetric] lift_un.simps b by auto + case False + hence "snd (f l1 u1) = None" + using or by auto + with False obtain b where b: "f l1 u1 = (Some b, None)" + by (cases "f l1 u1") auto + thus ?thesis + unfolding Sa[symmetric] lift_un.simps b by auto qed - thus False using lift_un_Some by auto + thus False + using lift_un_Some by auto qed - then obtain a' b' where f: "f l1 u1 = (Some a', Some b')" by (cases "f l1 u1", auto) + then obtain a' b' where f: "f l1 u1 = (Some a', Some b')" + by (cases "f l1 u1") auto from lift_un_Some[unfolded Sa[symmetric] lift_un.simps f] - have "Some l = fst (f l1 u1)" and "Some u = snd (f l1 u1)" unfolding f by auto - thus ?thesis unfolding Sa[symmetric] lift_un.simps using Pa[OF Sa] by auto + have "Some l = fst (f l1 u1)" and "Some u = snd (f l1 u1)" + unfolding f by auto + thus ?thesis + unfolding Sa[symmetric] lift_un.simps using Pa[OF Sa] by auto qed lemma lift_un: assumes lift_un_Some: "Some (l, u) = lift_un (approx' prec a bs) f" - and Pa: "\l u. Some (l, u) = approx prec a bs \ l \ interpret_floatarith a xs \ interpret_floatarith a xs \ u" (is "\l u. _ = ?g a \ ?P l u a") - shows "\ l1 u1. (l1 \ interpret_floatarith a xs \ interpret_floatarith a xs \ u1) \ + and Pa: "\l u. Some (l, u) = approx prec a bs \ + l \ interpret_floatarith a xs \ interpret_floatarith a xs \ u" + (is "\l u. _ = ?g a \ ?P l u a") + shows "\l1 u1. (l1 \ interpret_floatarith a xs \ interpret_floatarith a xs \ u1) \ Some l = fst (f l1 u1) \ Some u = snd (f l1 u1)" proof - - { fix l u assume "Some (l, u) = approx' prec a bs" - with approx_approx'[of prec a bs, OF _ this] Pa - have "l \ interpret_floatarith a xs \ interpret_floatarith a xs \ u" by auto } note Pa = this + have Pa: "l \ interpret_floatarith a xs \ interpret_floatarith a xs \ u" + if "Some (l, u) = approx' prec a bs" for l u + using approx_approx'[of prec a bs, OF _ that] Pa by auto from lift_un_f[where g="\a. approx' prec a bs" and P = ?P, OF lift_un_Some, OF Pa] show ?thesis by auto qed lemma lift_un_bnds: - assumes bnds: "\ (x::real) lx ux. (Some l, Some u) = f lx ux \ x \ { lx .. ux } \ l \ f' x \ f' x \ u" - and lift_un_Some: "Some (l, u) = lift_un (approx' prec a bs) f" - and Pa: "\l u. Some (l, u) = approx prec a bs \ l \ interpret_floatarith a xs \ interpret_floatarith a xs \ u" + assumes bnds: "\(x::real) lx ux. (Some l, Some u) = f lx ux \ x \ { lx .. ux } \ l \ f' x \ f' x \ u" + and lift_un_Some: "Some (l, u) = lift_un (approx' prec a bs) f" + and Pa: "\l u. Some (l, u) = approx prec a bs \ + l \ interpret_floatarith a xs \ interpret_floatarith a xs \ u" shows "real l \ f' (interpret_floatarith a xs) \ f' (interpret_floatarith a xs) \ real u" proof - from lift_un[OF lift_un_Some Pa] - obtain l1 u1 where "l1 \ interpret_floatarith a xs" and "interpret_floatarith a xs \ u1" and "Some l = fst (f l1 u1)" and "Some u = snd (f l1 u1)" by blast - hence "(Some l, Some u) = f l1 u1" and "interpret_floatarith a xs \ {l1 .. u1}" by auto - thus ?thesis using bnds by auto + obtain l1 u1 where "l1 \ interpret_floatarith a xs" + and "interpret_floatarith a xs \ u1" + and "Some l = fst (f l1 u1)" + and "Some u = snd (f l1 u1)" + by blast + hence "(Some l, Some u) = f l1 u1" and "interpret_floatarith a xs \ {l1 .. u1}" + by auto + thus ?thesis + using bnds by auto qed lemma approx: assumes "bounded_by xs vs" - and "Some (l, u) = approx prec arith vs" (is "_ = ?g arith") + and "Some (l, u) = approx prec arith vs" (is "_ = ?g arith") shows "l \ interpret_floatarith arith xs \ interpret_floatarith arith xs \ u" (is "?P l u arith") using \Some (l, u) = approx prec arith vs\ 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 = float_plus_down prec l1 l2" and "u = float_plus_up prec u1 u2" - "l1 \ interpret_floatarith a xs" and "interpret_floatarith a xs \ u1" - "l2 \ interpret_floatarith b xs" and "interpret_floatarith b xs \ u2" unfolding fst_conv snd_conv by blast - thus ?case unfolding interpret_floatarith.simps by (auto intro!: float_plus_up_le float_plus_down_le) + obtain l1 u1 l2 u2 where "l = float_plus_down prec l1 l2" + and "u = float_plus_up prec u1 u2" "l1 \ interpret_floatarith a xs" + and "interpret_floatarith a xs \ u1" "l2 \ interpret_floatarith b xs" + and "interpret_floatarith b xs \ u2" + unfolding fst_conv snd_conv by blast + 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 - obtain l1 u1 where "l = -u1" and "u = -l1" - "l1 \ interpret_floatarith a xs" and "interpret_floatarith a xs \ u1" unfolding fst_conv snd_conv by blast - thus ?case unfolding interpret_floatarith.simps using minus_float.rep_eq by auto + obtain l1 u1 where "l = -u1" "u = -l1" + and "l1 \ interpret_floatarith a xs" "interpret_floatarith a xs \ u1" + unfolding fst_conv snd_conv by blast + thus ?case + unfolding interpret_floatarith.simps using minus_float.rep_eq by auto next case (Mult a b) from lift_bin'[OF Mult.prems[unfolded approx.simps]] Mult.hyps @@ -2634,34 +3085,49 @@ and "l1 \ interpret_floatarith a xs" and "interpret_floatarith a xs \ u1" and "l2 \ interpret_floatarith b xs" and "interpret_floatarith b xs \ u2" unfolding fst_conv snd_conv by blast hence bnds: - "nprt l1 * pprt u2 + nprt u1 * nprt u2 + pprt l1 * pprt l2 + pprt u1 * nprt l2 \ interpret_floatarith (Mult a b) xs" (is "?l \ _") - "interpret_floatarith (Mult a b) xs \ pprt u1 * pprt u2 + pprt l1 * nprt u2 + nprt u1 * pprt l2 + nprt l1 * nprt l2" (is "_ \ ?u") + "nprt l1 * pprt u2 + nprt u1 * nprt u2 + pprt l1 * pprt l2 + pprt u1 * nprt l2 \ + interpret_floatarith (Mult a b) xs" (is "?l \ _") + "interpret_floatarith (Mult a b) xs \ + pprt u1 * pprt u2 + pprt l1 * nprt u2 + nprt u1 * pprt l2 + nprt l1 * nprt l2" (is "_ \ ?u") unfolding interpret_floatarith.simps l u using mult_le_prts mult_ge_prts by auto from l u have "l \ ?l" "?u \ u" by (auto intro!: float_plus_up_le float_plus_down_le) - thus ?case using bnds by simp + 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 obtain l1 u1 where l': "Some l = (if 0 < l1 \ u1 < 0 then Some (float_divl prec 1 u1) else None)" and u': "Some u = (if 0 < l1 \ u1 < 0 then Some (float_divr prec 1 l1) else None)" - and l1: "l1 \ interpret_floatarith a xs" and u1: "interpret_floatarith a xs \ u1" by blast - have either: "0 < l1 \ u1 < 0" proof (rule ccontr) assume P: "\ (0 < l1 \ u1 < 0)" show False using l' unfolding if_not_P[OF P] by auto qed - moreover have l1_le_u1: "real l1 \ real u1" using l1 u1 by auto - ultimately have "real l1 \ 0" and "real u1 \ 0" by auto + and l1: "l1 \ interpret_floatarith a xs" + and u1: "interpret_floatarith a xs \ u1" + by blast + have either: "0 < l1 \ u1 < 0" + proof (rule ccontr) + assume P: "\ (0 < l1 \ u1 < 0)" + show False + using l' unfolding if_not_P[OF P] by auto + qed + moreover have l1_le_u1: "real l1 \ real u1" + using l1 u1 by auto + ultimately have "real l1 \ 0" and "real u1 \ 0" + by auto have inv: "inverse u1 \ inverse (interpret_floatarith a xs) \ inverse (interpret_floatarith a xs) \ inverse l1" proof (cases "0 < l1") - case True hence "0 < real u1" and "0 < real l1" "0 < interpret_floatarith a xs" + case True + hence "0 < real u1" and "0 < real l1" "0 < interpret_floatarith a xs" using l1_le_u1 l1 by auto show ?thesis unfolding inverse_le_iff_le[OF \0 < real u1\ \0 < interpret_floatarith a xs\] inverse_le_iff_le[OF \0 < interpret_floatarith a xs\ \0 < real l1\] using l1 u1 by auto next - case False hence "u1 < 0" using either by blast + case False + hence "u1 < 0" + using either by blast hence "real u1 < 0" and "real l1 < 0" "interpret_floatarith a xs < 0" using l1_le_u1 u1 by auto show ?thesis @@ -2670,47 +3136,81 @@ using l1 u1 by auto qed - from l' have "l = float_divl prec 1 u1" by (cases "0 < l1 \ u1 < 0", auto) - hence "l \ inverse u1" unfolding nonzero_inverse_eq_divide[OF \real u1 \ 0\] using float_divl[of prec 1 u1] by auto - also have "\ \ inverse (interpret_floatarith a xs)" using inv by auto + from l' have "l = float_divl prec 1 u1" + by (cases "0 < l1 \ u1 < 0") auto + hence "l \ inverse u1" + unfolding nonzero_inverse_eq_divide[OF \real u1 \ 0\] + using float_divl[of prec 1 u1] by auto + also have "\ \ inverse (interpret_floatarith a xs)" + using inv by auto finally have "l \ inverse (interpret_floatarith a xs)" . moreover - from u' have "u = float_divr prec 1 l1" by (cases "0 < l1 \ u1 < 0", auto) - hence "inverse l1 \ u" unfolding nonzero_inverse_eq_divide[OF \real l1 \ 0\] using float_divr[of 1 l1 prec] by auto - hence "inverse (interpret_floatarith a xs) \ u" by (rule order_trans[OF inv[THEN conjunct2]]) - ultimately show ?case unfolding interpret_floatarith.simps using l1 u1 by auto + from u' have "u = float_divr prec 1 l1" + by (cases "0 < l1 \ u1 < 0") auto + hence "inverse l1 \ u" + unfolding nonzero_inverse_eq_divide[OF \real l1 \ 0\] + using float_divr[of 1 l1 prec] by auto + hence "inverse (interpret_floatarith a xs) \ u" + by (rule order_trans[OF inv[THEN conjunct2]]) + ultimately show ?case + unfolding interpret_floatarith.simps using l1 u1 by auto next case (Abs x) from lift_un'[OF Abs.prems[unfolded approx.simps], unfolded fst_conv snd_conv] Abs.hyps - obtain l1 u1 where l': "l = (if l1 < 0 \ 0 < u1 then 0 else min \l1\ \u1\)" and u': "u = max \l1\ \u1\" - and l1: "l1 \ interpret_floatarith x xs" and u1: "interpret_floatarith x xs \ u1" by blast - thus ?case unfolding l' u' by (cases "l1 < 0 \ 0 < u1", auto simp add: real_of_float_min real_of_float_max) + obtain l1 u1 where l': "l = (if l1 < 0 \ 0 < u1 then 0 else min \l1\ \u1\)" + and u': "u = max \l1\ \u1\" + and l1: "l1 \ interpret_floatarith x xs" + and u1: "interpret_floatarith x xs \ u1" + by blast + thus ?case + unfolding l' u' + by (cases "l1 < 0 \ 0 < u1") (auto simp add: real_of_float_min real_of_float_max) next case (Min a b) from lift_bin'[OF Min.prems[unfolded approx.simps], unfolded fst_conv snd_conv] Min.hyps obtain l1 u1 l2 u2 where l': "l = min l1 l2" and u': "u = min u1 u2" and l1: "l1 \ interpret_floatarith a xs" and u1: "interpret_floatarith a xs \ u1" - and l1: "l2 \ interpret_floatarith b xs" and u1: "interpret_floatarith b xs \ u2" by blast - thus ?case unfolding l' u' by (auto simp add: real_of_float_min) + and l1: "l2 \ interpret_floatarith b xs" and u1: "interpret_floatarith b xs \ u2" + by blast + thus ?case + unfolding l' u' by (auto simp add: real_of_float_min) next case (Max a b) from lift_bin'[OF Max.prems[unfolded approx.simps], unfolded fst_conv snd_conv] Max.hyps obtain l1 u1 l2 u2 where l': "l = max l1 l2" and u': "u = max u1 u2" and l1: "l1 \ interpret_floatarith a xs" and u1: "interpret_floatarith a xs \ u1" - and l1: "l2 \ interpret_floatarith b xs" and u1: "interpret_floatarith b xs \ u2" by blast - thus ?case unfolding l' u' by (auto simp add: real_of_float_max) -next case (Cos a) with lift_un'_bnds[OF bnds_cos] show ?case by auto -next case (Arctan a) with lift_un'_bnds[OF bnds_arctan] show ?case by auto -next case Pi with pi_boundaries show ?case by auto -next case (Sqrt a) with lift_un'_bnds[OF bnds_sqrt] show ?case by auto -next case (Exp a) with lift_un'_bnds[OF bnds_exp] show ?case by auto -next case (Ln a) with lift_un_bnds[OF bnds_ln] show ?case by auto -next case (Power a n) with lift_un'_bnds[OF bnds_power] show ?case by auto -next case (Num f) thus ?case by auto + and l1: "l2 \ interpret_floatarith b xs" and u1: "interpret_floatarith b xs \ u2" + by blast + thus ?case + unfolding l' u' by (auto simp add: real_of_float_max) +next + case (Cos a) + with lift_un'_bnds[OF bnds_cos] show ?case by auto +next + case (Arctan a) + with lift_un'_bnds[OF bnds_arctan] show ?case by auto +next + case Pi + with pi_boundaries show ?case by auto +next + case (Sqrt a) + with lift_un'_bnds[OF bnds_sqrt] show ?case by auto +next + case (Exp a) + with lift_un'_bnds[OF bnds_exp] show ?case by auto +next + case (Ln a) + with lift_un_bnds[OF bnds_ln] show ?case by auto +next + case (Power a n) + with lift_un'_bnds[OF bnds_power] show ?case by auto +next + case (Num f) + thus ?case by auto next case (Var n) from this[symmetric] \bounded_by xs vs\[THEN bounded_byE, of n] - show ?case by (cases "n < length vs", auto) + show ?case by (cases "n < length vs") auto qed datatype form = Bound floatarith floatarith floatarith form @@ -2762,7 +3262,8 @@ lemma lazy_conj: "(if A then B else False) = (A \ B)" by simp lemma approx_form_approx_form': - assumes "approx_form' prec f s n l u bs ss" and "(x::real) \ { l .. u }" + assumes "approx_form' prec f s n l u bs ss" + and "(x::real) \ { l .. u }" obtains l' u' where "x \ { l' .. u' }" and "approx_form prec f (bs[n := Some (l', u')]) ss" using assms proof (induct s arbitrary: l u) @@ -2805,23 +3306,27 @@ and approx_form': "approx_form' prec f (ss ! n) n l u vs ss" by (cases "approx prec a vs", simp) (cases "approx prec b vs", auto) - { assume "xs ! n \ { interpret_floatarith a xs .. interpret_floatarith b xs }" - with approx[OF Bound.prems(2) l_eq] and approx[OF Bound.prems(2) u_eq] + have "interpret_form f xs" + if "xs ! n \ { interpret_floatarith a xs .. interpret_floatarith b xs }" + proof - + from approx[OF Bound.prems(2) l_eq] and approx[OF Bound.prems(2) u_eq] that have "xs ! n \ { l .. u}" by auto from approx_form_approx_form'[OF approx_form' this] obtain lx ux where bnds: "xs ! n \ { lx .. ux }" and approx_form: "approx_form prec f (vs[n := Some (lx, ux)]) ss" . - from \bounded_by xs vs\ bnds - have "bounded_by xs (vs[n := Some (lx, ux)])" by (rule bounded_by_update) - with Bound.hyps[OF approx_form] - have "interpret_form f xs" by blast } - thus ?case using interpret_form.simps x_eq and interpret_floatarith.simps by simp + from \bounded_by xs vs\ bnds have "bounded_by xs (vs[n := Some (lx, ux)])" + by (rule bounded_by_update) + with Bound.hyps[OF approx_form] show ?thesis + by blast + qed + thus ?case + using interpret_form.simps x_eq and interpret_floatarith.simps by simp next case (Assign x a f) - then obtain n - where x_eq: "x = Var n" by (cases x) auto + then obtain n where x_eq: "x = Var n" + by (cases x) auto with Assign.prems obtain l u where bnd_eq: "Some (l, u) = approx prec a vs" @@ -2829,26 +3334,29 @@ and approx_form': "approx_form' prec f (ss ! n) n l u vs ss" by (cases "approx prec a vs") auto - { assume bnds: "xs ! n = interpret_floatarith a xs" - with approx[OF Assign.prems(2) bnd_eq] + have "interpret_form f xs" + if bnds: "xs ! n = interpret_floatarith a xs" + proof - + from approx[OF Assign.prems(2) bnd_eq] bnds have "xs ! n \ { l .. u}" by auto from approx_form_approx_form'[OF approx_form' this] obtain lx ux where bnds: "xs ! n \ { lx .. ux }" and approx_form: "approx_form prec f (vs[n := Some (lx, ux)]) ss" . - from \bounded_by xs vs\ bnds - have "bounded_by xs (vs[n := Some (lx, ux)])" by (rule bounded_by_update) - with Assign.hyps[OF approx_form] - have "interpret_form f xs" by blast } - thus ?case using interpret_form.simps x_eq and interpret_floatarith.simps by simp + from \bounded_by xs vs\ bnds have "bounded_by xs (vs[n := Some (lx, ux)])" + by (rule bounded_by_update) + with Assign.hyps[OF approx_form] show ?thesis + by blast + qed + thus ?case + using interpret_form.simps x_eq and interpret_floatarith.simps by simp next case (Less 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: "real (float_plus_up prec u (-l')) < 0" - by (cases "approx prec a vs", auto, - cases "approx prec b vs", auto) + by (cases "approx prec a vs", auto, cases "approx prec b vs", auto) 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 @@ -2858,8 +3366,7 @@ where l_eq: "Some (l, u) = approx prec a vs" and u_eq: "Some (l', u') = approx prec b vs" and inequality: "real (float_plus_up prec u (-l')) \ 0" - by (cases "approx prec a vs", auto, - cases "approx prec b vs", auto) + by (cases "approx prec a vs", auto, cases "approx prec b vs", auto) 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 @@ -2880,10 +3387,11 @@ lemma approx_form: assumes "n = length xs" - assumes "approx_form prec f (replicate n None) ss" + and "approx_form prec f (replicate n None) ss" shows "interpret_form f xs" using approx_form_aux[OF _ bounded_by_None] assms by auto + subsection \Implementing Taylor series expansion\ fun isDERIV :: "nat \ floatarith \ real list \ bool" where @@ -3007,64 +3515,72 @@ have "0 < interpret_floatarith a xs" by auto thus ?case using Sqrt by auto next - case (Power a n) thus ?case by (cases n) auto + case (Power a n) + thus ?case by (cases n) auto qed auto lemma bounded_by_update_var: - assumes "bounded_by xs vs" and "vs ! i = Some (l, u)" + assumes "bounded_by xs vs" + and "vs ! i = Some (l, u)" and bnd: "x \ { real l .. real u }" shows "bounded_by (xs[i := x]) vs" proof (cases "i < length xs") case False - thus ?thesis using \bounded_by xs vs\ by auto + thus ?thesis + using \bounded_by xs vs\ by auto next + case True let ?xs = "xs[i := x]" - case True hence "i < length ?xs" by auto - { - fix j - assume "j < length vs" - have "case vs ! j of None \ True | Some (l, u) \ ?xs ! j \ { real l .. real u }" - proof (cases "vs ! j") - case (Some b) + from True have "i < length ?xs" by auto + have "case vs ! j of None \ True | Some (l, u) \ ?xs ! j \ {real l .. real u}" + if "j < length vs" for j + proof (cases "vs ! j") + case None + then show ?thesis by simp + next + case (Some b) + thus ?thesis + proof (cases "i = j") + case True + thus ?thesis using \vs ! i = Some (l, u)\ Some and bnd \i < length ?xs\ + by auto + next + case False thus ?thesis - proof (cases "i = j") - case True - thus ?thesis using \vs ! i = Some (l, u)\ Some and bnd \i < length ?xs\ - by auto - next - case False - thus ?thesis - using \bounded_by xs vs\[THEN bounded_byE, OF \j < length vs\] Some by auto - qed - qed auto - } - thus ?thesis unfolding bounded_by_def by auto + using \bounded_by xs vs\[THEN bounded_byE, OF \j < length vs\] Some by auto + qed + qed + thus ?thesis + unfolding bounded_by_def by auto qed lemma isDERIV_approx': assumes "bounded_by xs vs" - and vs_x: "vs ! x = Some (l, u)" and X_in: "X \ { real l .. real u }" + and vs_x: "vs ! x = Some (l, u)" + and X_in: "X \ {real l .. real u}" and approx: "isDERIV_approx prec x f vs" shows "isDERIV x f (xs[x := X])" proof - - note bounded_by_update_var[OF \bounded_by xs vs\ vs_x X_in] approx - thus ?thesis by (rule isDERIV_approx) + from bounded_by_update_var[OF \bounded_by xs vs\ vs_x X_in] approx + show ?thesis by (rule isDERIV_approx) qed lemma DERIV_approx: - assumes "n < length xs" and bnd: "bounded_by xs vs" + assumes "n < length xs" + and bnd: "bounded_by xs vs" and isD: "isDERIV_approx prec n f vs" and app: "Some (l, u) = approx prec (DERIV_floatarith n f) vs" (is "_ = approx _ ?D _") shows "\(x::real). l \ x \ x \ u \ DERIV (\ x. interpret_floatarith f (xs[n := x])) (xs!n) :> x" (is "\ x. _ \ _ \ DERIV (?i f) _ :> _") proof (rule exI[of _ "?i ?D (xs!n)"], rule conjI[OF _ conjI]) - let "?i f x" = "interpret_floatarith f (xs[n := x])" + let "?i f" = "\x. interpret_floatarith f (xs[n := x])" from approx[OF bnd app] show "l \ ?i ?D (xs!n)" and "?i ?D (xs!n) \ u" using \n < length xs\ by auto from DERIV_floatarith[OF \n < length xs\, of f "xs!n"] isDERIV_approx[OF bnd isD] - show "DERIV (?i f) (xs!n) :> (?i ?D (xs!n))" by simp + show "DERIV (?i f) (xs!n) :> (?i ?D (xs!n))" + by simp qed fun lift_bin :: "(float * float) option \ @@ -3099,23 +3615,25 @@ and x: "x \ { real l .. real u }" shows "bounded_by (x#xs) ((Some (l, u))#vs)" proof - - { - fix i assume *: "i < length ((Some (l, u))#vs)" - have "case ((Some (l,u))#vs) ! i of Some (l, u) \ (x#xs)!i \ { real l .. real u } | None \ True" - proof (cases i) - case 0 with x show ?thesis by auto - next - case (Suc i) with * have "i < length vs" by auto - from bnd[THEN bounded_byE, OF this] - show ?thesis unfolding Suc nth_Cons_Suc . - qed - } - thus ?thesis by (auto simp add: bounded_by_def) + have "case ((Some (l,u))#vs) ! i of Some (l, u) \ (x#xs)!i \ { real l .. real u } | None \ True" + if *: "i < length ((Some (l, u))#vs)" for i + proof (cases i) + case 0 + with x show ?thesis by auto + next + case (Suc i) + with * have "i < length vs" by auto + from bnd[THEN bounded_byE, OF this] + show ?thesis unfolding Suc nth_Cons_Suc . + qed + thus ?thesis + by (auto simp add: bounded_by_def) qed lemma approx_tse_generic: assumes "bounded_by xs vs" - and bnd_c: "bounded_by (xs[x := c]) vs" and "x < length vs" and "x < length xs" + and bnd_c: "bounded_by (xs[x := c]) vs" + and "x < length vs" and "x < length xs" and bnd_x: "vs ! x = Some (lx, ux)" and ate: "Some (l, u) = approx_tse prec x s c k f vs" shows "\ n. (\ m < n. \ (z::real) \ {lx .. ux}. @@ -3127,7 +3645,8 @@ inverse (real (\ j \ {k.. {l .. u})" (is "\ n. ?taylor f k l u n") -using ate proof (induct s arbitrary: k f l u) + using ate +proof (induct s arbitrary: k f l u) case 0 { fix t::real assume "t \ {lx .. ux}" @@ -3173,39 +3692,37 @@ (is "?f 0 (real c) \ _") by auto - { - fix f :: "'a \ 'a" fix n :: nat fix x :: 'a - have "(f ^^ Suc n) x = (f ^^ n) (f x)" - by (induct n) auto - } - note funpow_Suc = this[symmetric] - from Suc.hyps[OF ate, unfolded this] - obtain n - where DERIV_hyp: "\ m z. \ m < n ; (z::real) \ { lx .. ux } \ \ DERIV (?f (Suc m)) z :> ?f (Suc (Suc m)) z" - and hyp: "\ t \ {real lx .. real ux}. (\ i = 0.. j \ {Suc k.. j \ {Suc k.. {l2 .. u2}" + have funpow_Suc[symmetric]: "(f ^^ Suc n) x = (f ^^ n) (f x)" + for f :: "'a \ 'a" and n :: nat and x :: 'a + by (induct n) auto + from Suc.hyps[OF ate, unfolded this] obtain n + where DERIV_hyp: "\m z. \ m < n ; (z::real) \ { lx .. ux } \ \ + DERIV (?f (Suc m)) z :> ?f (Suc (Suc m)) z" + and hyp: "\t \ {real lx .. real ux}. + (\ i = 0.. j \ {Suc k.. j \ {Suc k.. {l2 .. u2}" (is "\ t \ _. ?X (Suc k) f n t \ _") by blast - { - fix m and z::real - assume "m < Suc n" and bnd_z: "z \ { lx .. ux }" - have "DERIV (?f m) z :> ?f (Suc m) z" - proof (cases m) - case 0 - with DERIV_floatarith[OF \x < length xs\ isDERIV_approx'[OF \bounded_by xs vs\ bnd_x bnd_z True]] - show ?thesis by simp - next - case (Suc m') - hence "m' < n" using \m < Suc n\ by auto - from DERIV_hyp[OF this bnd_z] - show ?thesis using Suc by simp - qed - } note DERIV = this - - have "\ k i. k < i \ {k ..< i} = insert k {Suc k ..< i}" by auto - hence setprod_head_Suc: "\ k i. \ {k ..< k + Suc i} = k * \ {Suc k ..< Suc k + i}" by auto - have setsum_move0: "\ k F. setsum F {0.. k. F (Suc k)) {0.. ?f (Suc m) z" + if "m < Suc n" and bnd_z: "z \ { lx .. ux }" for m and z::real + proof (cases m) + case 0 + with DERIV_floatarith[OF \x < length xs\ + isDERIV_approx'[OF \bounded_by xs vs\ bnd_x bnd_z True]] + show ?thesis by simp + next + case (Suc m') + hence "m' < n" + using \m < Suc n\ by auto + from DERIV_hyp[OF this bnd_z] show ?thesis + using Suc by simp + qed + + have "\k i. k < i \ {k ..< i} = insert k {Suc k ..< i}" by auto + hence setprod_head_Suc: "\k i. \{k ..< k + Suc i} = k * \{Suc k ..< Suc k + i}" + by auto + have setsum_move0: "\k F. setsum F {0.. k. F (Suc k)) {0.. "xs!x - c" @@ -3240,12 +3757,13 @@ lemma approx_tse: assumes "bounded_by xs vs" - and bnd_x: "vs ! x = Some (lx, ux)" and bnd_c: "real c \ {lx .. ux}" + and bnd_x: "vs ! x = Some (lx, ux)" + and bnd_c: "real c \ {lx .. ux}" and "x < length vs" and "x < length xs" and ate: "Some (l, u) = approx_tse prec x s c 1 f vs" - shows "interpret_floatarith f xs \ { l .. u }" + shows "interpret_floatarith f xs \ {l .. u}" proof - - def F \ "\ n z. interpret_floatarith ((DERIV_floatarith x ^^ n) f) (xs[x := z])" + def F \ "\n z. interpret_floatarith ((DERIV_floatarith x ^^ n) f) (xs[x := z])" hence F0: "F 0 = (\ z. interpret_floatarith f (xs[x := z]))" by auto hence "bounded_by (xs[x := c]) vs" and "x < length vs" "x < length xs" @@ -3267,17 +3785,19 @@ show ?thesis proof (cases n) - case 0 thus ?thesis using hyp[OF bnd_xs] unfolding F_def by auto + case 0 + thus ?thesis + using hyp[OF bnd_xs] unfolding F_def by auto next case (Suc n') show ?thesis proof (cases "xs ! x = c") case True from True[symmetric] hyp[OF bnd_xs] Suc show ?thesis - unfolding F_def Suc setsum_head_upt_Suc[OF zero_less_Suc] setsum_shift_bounds_Suc_ivl by auto + unfolding F_def Suc setsum_head_upt_Suc[OF zero_less_Suc] setsum_shift_bounds_Suc_ivl + by auto next case False - have "lx \ real c" "real c \ ux" "lx \ xs!x" "xs!x \ ux" using Suc bnd_c \bounded_by xs vs\[THEN bounded_byE, OF \x < length vs\] bnd_x by auto from Taylor.taylor[OF zero_less_Suc, of F, OF F0 DERIV[unfolded Suc] this False] @@ -3288,12 +3808,14 @@ unfolding atLeast0LessThan by blast from t_bnd bnd_xs bnd_c have *: "t \ {lx .. ux}" - by (cases "xs ! x < c", auto) + by (cases "xs ! x < c") auto have "interpret_floatarith f (xs[x := xs ! x]) = ?taylor t" unfolding fl_eq Suc by (auto simp add: algebra_simps divide_inverse) - also have "\ \ {l .. u}" using * by (rule hyp) - finally show ?thesis by simp + also have "\ \ {l .. u}" + using * by (rule hyp) + finally show ?thesis + by simp qed qed qed @@ -3309,10 +3831,12 @@ lemma approx_tse_form': fixes x :: real - assumes "approx_tse_form' prec t f s l u cmp" and "x \ {l .. u}" - shows "\ l' u' ly uy. x \ { l' .. u' } \ real l \ l' \ u' \ real u \ cmp ly uy \ - approx_tse prec 0 t ((l' + u') * Float 1 (- 1)) 1 f [Some (l', u')] = Some (ly, uy)" -using assms proof (induct s arbitrary: l u) + assumes "approx_tse_form' prec t f s l u cmp" + and "x \ {l .. u}" + shows "\l' u' ly uy. x \ {l' .. u'} \ real l \ l' \ u' \ real u \ cmp ly uy \ + approx_tse prec 0 t ((l' + u') * Float 1 (- 1)) 1 f [Some (l', u')] = Some (ly, uy)" + using assms +proof (induct s arbitrary: l u) case 0 then obtain ly uy where *: "approx_tse prec 0 t ((l + u) * Float 1 (- 1)) 1 f [Some (l, u)] = Some (ly, uy)" @@ -3328,69 +3852,73 @@ have m_l: "real l \ ?m" and m_u: "?m \ real u" unfolding less_eq_float_def using Suc.prems by auto - - with \x \ { l .. u }\ - have "x \ { l .. ?m} \ x \ { ?m .. u }" by auto + with \x \ { l .. u }\ consider "x \ { l .. ?m}" | "x \ {?m .. u}" + by atomize_elim auto thus ?case - proof (rule disjE) - assume "x \ { l .. ?m}" + proof cases + case 1 from Suc.hyps[OF l this] - obtain l' u' ly uy - where "x \ { l' .. u' } \ real l \ l' \ real u' \ ?m \ cmp ly uy \ - approx_tse prec 0 t ((l' + u') * Float 1 (- 1)) 1 f [Some (l', u')] = Some (ly, uy)" by blast - with m_u show ?thesis by (auto intro!: exI) + obtain l' u' ly uy where + "x \ {l' .. u'} \ real l \ l' \ real u' \ ?m \ cmp ly uy \ + approx_tse prec 0 t ((l' + u') * Float 1 (- 1)) 1 f [Some (l', u')] = Some (ly, uy)" + by blast + with m_u show ?thesis + by (auto intro!: exI) next - assume "x \ { ?m .. u }" + case 2 from Suc.hyps[OF u this] - obtain l' u' ly uy - where "x \ { l' .. u' } \ ?m \ real l' \ u' \ real u \ cmp ly uy \ - approx_tse prec 0 t ((l' + u') * Float 1 (- 1)) 1 f [Some (l', u')] = Some (ly, uy)" by blast - with m_u show ?thesis by (auto intro!: exI) + obtain l' u' ly uy where + "x \ { l' .. u' } \ ?m \ real l' \ u' \ real u \ cmp ly uy \ + approx_tse prec 0 t ((l' + u') * Float 1 (- 1)) 1 f [Some (l', u')] = Some (ly, uy)" + by blast + with m_u show ?thesis + by (auto intro!: exI) qed qed lemma approx_tse_form'_less: fixes x :: real assumes tse: "approx_tse_form' prec t (Add a (Minus b)) s l u (\ l u. 0 < l)" - and x: "x \ {l .. u}" + and x: "x \ {l .. u}" shows "interpret_floatarith b [x] < interpret_floatarith a [x]" proof - from approx_tse_form'[OF tse x] obtain l' u' ly uy - where x': "x \ { l' .. u' }" and "l \ real l'" + where x': "x \ {l' .. u'}" + and "l \ real l'" and "real u' \ u" and "0 < ly" and tse: "approx_tse prec 0 t ((l' + u') * Float 1 (- 1)) 1 (Add a (Minus b)) [Some (l', u')] = Some (ly, uy)" by blast - hence "bounded_by [x] [Some (l', u')]" by (auto simp add: bounded_by_def) - + hence "bounded_by [x] [Some (l', u')]" + by (auto simp add: bounded_by_def) from approx_tse[OF this _ _ _ _ tse[symmetric], of l' u'] x' have "ly \ interpret_floatarith a [x] - interpret_floatarith b [x]" by auto - from order_less_le_trans[OF _ this, of 0] \0 < ly\ - show ?thesis by auto + from order_less_le_trans[OF _ this, of 0] \0 < ly\ show ?thesis + by auto qed lemma approx_tse_form'_le: fixes x :: real assumes tse: "approx_tse_form' prec t (Add a (Minus b)) s l u (\ l u. 0 \ l)" - and x: "x \ {l .. u}" + and x: "x \ {l .. u}" shows "interpret_floatarith b [x] \ interpret_floatarith a [x]" proof - from approx_tse_form'[OF tse x] obtain l' u' ly uy - where x': "x \ { l' .. u' }" and "l \ real l'" + where x': "x \ {l' .. u'}" + and "l \ real l'" and "real u' \ u" and "0 \ ly" and tse: "approx_tse prec 0 t ((l' + u') * Float 1 (- 1)) 1 (Add a (Minus b)) [Some (l', u')] = Some (ly, uy)" by blast hence "bounded_by [x] [Some (l', u')]" by (auto simp add: bounded_by_def) - from approx_tse[OF this _ _ _ _ tse[symmetric], of l' u'] x' have "ly \ interpret_floatarith a [x] - interpret_floatarith b [x]" by auto - from order_trans[OF _ this, of 0] \0 \ ly\ - show ?thesis by auto + from order_trans[OF _ this, of 0] \0 \ ly\ show ?thesis + by auto qed fun approx_tse_concl where @@ -3408,34 +3936,37 @@ "approx_tse_concl _ _ _ _ _ _ _ _ \ False" definition -"approx_tse_form prec t s f = - (case f - of (Bound x a b f) \ x = Var 0 \ - (case (approx prec a [None], approx prec b [None]) - of (Some (l, u), Some (l', u')) \ approx_tse_concl prec t f s l u l' u' - | _ \ False) - | _ \ False)" + "approx_tse_form prec t s f = + (case f of + Bound x a b f \ + x = Var 0 \ + (case (approx prec a [None], approx prec b [None]) of + (Some (l, u), Some (l', u')) \ approx_tse_concl prec t f s l u l' u' + | _ \ False) + | _ \ False)" lemma approx_tse_form: assumes "approx_tse_form prec t s f" shows "interpret_form f [x]" proof (cases f) - case (Bound i a b f') note f_def = this + case f_def: (Bound i a b f') with assms obtain l u l' u' where a: "approx prec a [None] = Some (l, u)" and b: "approx prec b [None] = Some (l', u')" unfolding approx_tse_form_def by (auto elim!: case_optionE) - from Bound assms have "i = Var 0" unfolding approx_tse_form_def by auto + from f_def assms have "i = Var 0" + unfolding approx_tse_form_def by auto hence i: "interpret_floatarith i [x] = x" by auto - { let "?f z" = "interpret_floatarith z [x]" + { + let ?f = "\z. interpret_floatarith z [x]" assume "?f i \ { ?f a .. ?f b }" with approx[OF _ a[symmetric], of "[x]"] approx[OF _ b[symmetric], of "[x]"] have bnd: "x \ { l .. u'}" unfolding bounded_by_def i by auto have "interpret_form f' [x]" - using assms[unfolded Bound] + using assms[unfolded f_def] proof (induct f') case (Less lf rt) with a b @@ -3445,21 +3976,22 @@ show ?case using Less by auto next case (LessEqual lf rt) - with Bound a b assms + with f_def a b assms have "approx_tse_form' prec t (Add rt (Minus lf)) s l u' (\ l u. 0 \ l)" unfolding approx_tse_form_def by auto from approx_tse_form'_le[OF this bnd] show ?case using LessEqual by auto next case (AtLeastAtMost x lf rt) - with Bound a b assms + with f_def a b assms have "approx_tse_form' prec t (Add rt (Minus x)) s l u' (\ l u. 0 \ l)" and "approx_tse_form' prec t (Add x (Minus lf)) s l u' (\ l u. 0 \ l)" 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 ?case using AtLeastAtMost by auto qed (auto simp: f_def approx_tse_form_def elim!: case_optionE) - } thus ?thesis unfolding f_def by auto + } + 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.\ @@ -3479,6 +4011,7 @@ bs @ [approx prec x bs, approx prec a bs, approx prec b bs]" | "approx_form_eval _ _ bs = bs" + subsection \Implement proof method \texttt{approximation}\ lemmas interpret_form_equations = interpret_form.simps interpret_floatarith.simps interpret_floatarith_num @@ -3593,17 +4126,16 @@ ML_file "approximation.ML" method_setup approximation = \ - let val free = Args.context -- Args.term >> (fn (_, Free (n, _)) => n | (ctxt, t) => - error ("Bad free variable: " ^ Syntax.string_of_term ctxt t)); + let + val free = + Args.context -- Args.term >> (fn (_, Free (n, _)) => n | (ctxt, t) => + error ("Bad free variable: " ^ Syntax.string_of_term ctxt t)); in - Scan.lift Parse.nat - -- + Scan.lift Parse.nat -- Scan.optional (Scan.lift (Args.$$$ "splitting" |-- Args.colon) - |-- Parse.and_list' (free --| Scan.lift (Args.$$$ "=") -- Scan.lift Parse.nat)) [] - -- - Scan.option (Scan.lift (Args.$$$ "taylor" |-- Args.colon) - |-- (free |-- Scan.lift (Args.$$$ "=") |-- Scan.lift Parse.nat)) - >> + |-- Parse.and_list' (free --| Scan.lift (Args.$$$ "=") -- Scan.lift Parse.nat)) [] -- + Scan.option (Scan.lift (Args.$$$ "taylor" |-- Args.colon) |-- + (free |-- Scan.lift (Args.$$$ "=") |-- Scan.lift Parse.nat)) >> (fn ((prec, splitting), taylor) => fn ctxt => SIMPLE_METHOD' (Approximation.approximation_tac prec splitting taylor ctxt)) end @@ -3613,7 +4145,6 @@ section "Quickcheck Generator" ML_file "approximation_generator.ML" - setup "Approximation_Generator.setup" end