(* Author:     Johannes Hoelzl <hoelzl@in.tum.de> 2008 / 2009 *)

header {* Prove Real Valued Inequalities by Computation *}

theory Approximation

imports Complex_Main Float Reflection Dense_Linear_Order Efficient_Nat

begin

section "Horner Scheme"

subsection {* Define auxiliary helper @{text horner} function *}

fun horner :: "(nat \<Rightarrow> nat) \<Rightarrow> (nat \<Rightarrow> nat \<Rightarrow> nat) \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> real \<Rightarrow> real" where

"horner F G 0 i k x       = 0" |

"horner F G (Suc n) i k x = 1 / real k - x * horner F G n (F i) (G i k) x"

lemma horner_schema': fixes x :: real  and a :: "nat \<Rightarrow> real"

  shows "a 0 - x * (\<Sum> i=0..<n. (-1)^i * a (Suc i) * x^i) = (\<Sum> i=0..<Suc n. (-1)^i * a i * x^i)"

proof -

  have shift_pow: "\<And>i. - (x * ((-1)^i * a (Suc i) * x ^ i)) = (-1)^(Suc i) * a (Suc i) * x ^ (Suc i)" by auto

  show ?thesis unfolding setsum_right_distrib shift_pow real_diff_def setsum_negf[symmetric] setsum_head_upt_Suc[OF zero_less_Suc]

    setsum_reindex[OF inj_Suc, unfolded comp_def, symmetric, of "\<lambda> n. (-1)^n  *a n * x^n"] by auto

qed

lemma horner_schema: fixes f :: "nat \<Rightarrow> nat" and G :: "nat \<Rightarrow> nat \<Rightarrow> nat" and F :: "nat \<Rightarrow> nat"

  assumes f_Suc: "\<And>n. f (Suc n) = G ((F o^ n) s) (f n)"

  shows "horner F G n ((F o^ j') s) (f j') x = (\<Sum> j = 0..< n. -1 ^ j * (1 / real (f (j' + j))) * x ^ j)"

proof (induct n arbitrary: i k j')

  case (Suc n)

  show ?case unfolding horner.simps Suc[where j'="Suc j'", unfolded funpow.simps comp_def f_Suc]

    using horner_schema'[of "\<lambda> j. 1 / real (f (j' + j))"] by auto

qed auto

lemma horner_bounds':

  assumes "0 \<le> Ifloat x" and f_Suc: "\<And>n. f (Suc n) = G ((F o^ n) s) (f n)"

  and lb_0: "\<And> i k x. lb 0 i k x = 0"

  and lb_Suc: "\<And> n i k x. lb (Suc n) i k x = lapprox_rat prec 1 (int k) - x * (ub n (F i) (G i k) x)"

  and ub_0: "\<And> i k x. ub 0 i k x = 0"

  and ub_Suc: "\<And> n i k x. ub (Suc n) i k x = rapprox_rat prec 1 (int k) - x * (lb n (F i) (G i k) x)"

  shows "Ifloat (lb n ((F o^ j') s) (f j') x) \<le> horner F G n ((F o^ j') s) (f j') (Ifloat x) \<and> 

         horner F G n ((F o^ j') s) (f j') (Ifloat x) \<le> Ifloat (ub n ((F o^ j') s) (f j') x)"

  (is "?lb n j' \<le> ?horner n j' \<and> ?horner n j' \<le> ?ub n j'")

proof (induct n arbitrary: j')

  case 0 thus ?case unfolding lb_0 ub_0 horner.simps by auto

next

  case (Suc n)

  have "?lb (Suc n) j' \<le> ?horner (Suc n) j'" unfolding lb_Suc ub_Suc horner.simps Ifloat_sub diff_def

  proof (rule add_mono)

    show "Ifloat (lapprox_rat prec 1 (int (f j'))) \<le> 1 / real (f j')" using lapprox_rat[of prec 1  "int (f j')"] by auto

    from Suc[where j'="Suc j'", unfolded funpow.simps comp_def f_Suc, THEN conjunct2] `0 \<le> Ifloat x`

    show "- Ifloat (x * ub n (F ((F o^ j') s)) (G ((F o^ j') s) (f j')) x) \<le> - (Ifloat x * horner F G n (F ((F o^ j') s)) (G ((F o^ j') s) (f j')) (Ifloat x))"

      unfolding Ifloat_mult neg_le_iff_le by (rule mult_left_mono)

  qed

  moreover have "?horner (Suc n) j' \<le> ?ub (Suc n) j'" unfolding ub_Suc ub_Suc horner.simps Ifloat_sub diff_def

  proof (rule add_mono)

    show "1 / real (f j') \<le> Ifloat (rapprox_rat prec 1 (int (f j')))" using rapprox_rat[of 1 "int (f j')" prec] by auto

    from Suc[where j'="Suc j'", unfolded funpow.simps comp_def f_Suc, THEN conjunct1] `0 \<le> Ifloat x`

    show "- (Ifloat x * horner F G n (F ((F o^ j') s)) (G ((F o^ j') s) (f j')) (Ifloat x)) \<le> 

          - Ifloat (x * lb n (F ((F o^ j') s)) (G ((F o^ j') s) (f j')) x)"

      unfolding Ifloat_mult neg_le_iff_le by (rule mult_left_mono)

  qed

  ultimately show ?case by blast

qed

subsection "Theorems for floating point functions implementing the horner scheme"

text {*

Here @{term_type "f :: nat \<Rightarrow> nat"} is the sequence defining the Taylor series, the coefficients are

all alternating and reciprocs. We use @{term G} and @{term F} to describe the computation of @{term f}.

*}

lemma horner_bounds: fixes F :: "nat \<Rightarrow> nat" and G :: "nat \<Rightarrow> nat \<Rightarrow> nat"

  assumes "0 \<le> Ifloat x" and f_Suc: "\<And>n. f (Suc n) = G ((F o^ n) s) (f n)"

  and lb_0: "\<And> i k x. lb 0 i k x = 0"

  and lb_Suc: "\<And> n i k x. lb (Suc n) i k x = lapprox_rat prec 1 (int k) - x * (ub n (F i) (G i k) x)"

  and ub_0: "\<And> i k x. ub 0 i k x = 0"

  and ub_Suc: "\<And> n i k x. ub (Suc n) i k x = rapprox_rat prec 1 (int k) - x * (lb n (F i) (G i k) x)"

  shows "Ifloat (lb n ((F o^ j') s) (f j') x) \<le> (\<Sum>j=0..<n. -1 ^ j * (1 / real (f (j' + j))) * Ifloat x ^ j)" (is "?lb") and 

    "(\<Sum>j=0..<n. -1 ^ j * (1 / real (f (j' + j))) * (Ifloat x ^ j)) \<le> Ifloat (ub n ((F o^ j') s) (f j') x)" (is "?ub")

proof -

  have "?lb  \<and> ?ub" 

    using horner_bounds'[where lb=lb, OF `0 \<le> Ifloat x` f_Suc lb_0 lb_Suc ub_0 ub_Suc]

    unfolding horner_schema[where f=f, OF f_Suc] .

  thus "?lb" and "?ub" by auto

qed

lemma horner_bounds_nonpos: fixes F :: "nat \<Rightarrow> nat" and G :: "nat \<Rightarrow> nat \<Rightarrow> nat"

  assumes "Ifloat x \<le> 0" and f_Suc: "\<And>n. f (Suc n) = G ((F o^ n) s) (f n)"

  and lb_0: "\<And> i k x. lb 0 i k x = 0"

  and lb_Suc: "\<And> n i k x. lb (Suc n) i k x = lapprox_rat prec 1 (int k) + x * (ub n (F i) (G i k) x)"

  and ub_0: "\<And> i k x. ub 0 i k x = 0"

  and ub_Suc: "\<And> n i k x. ub (Suc n) i k x = rapprox_rat prec 1 (int k) + x * (lb n (F i) (G i k) x)"

  shows "Ifloat (lb n ((F o^ j') s) (f j') x) \<le> (\<Sum>j=0..<n. (1 / real (f (j' + j))) * Ifloat x ^ j)" (is "?lb") and 

    "(\<Sum>j=0..<n. (1 / real (f (j' + j))) * (Ifloat x ^ j)) \<le> Ifloat (ub n ((F o^ j') s) (f j') x)" (is "?ub")

proof -

  { fix x y z :: float have "x - y * z = x + - y * z"

      by (cases x, cases y, cases z, simp add: plus_float.simps minus_float_def uminus_float.simps times_float.simps algebra_simps)

  } note diff_mult_minus = this

  { fix x :: float have "- (- x) = x" by (cases x, auto simp add: uminus_float.simps) } note minus_minus = this

  have move_minus: "Ifloat (-x) = -1 * Ifloat x" by auto

  have sum_eq: "(\<Sum>j=0..<n. (1 / real (f (j' + j))) * Ifloat x ^ j) = 

    (\<Sum>j = 0..<n. -1 ^ j * (1 / real (f (j' + j))) * Ifloat (- x) ^ j)"

  proof (rule setsum_cong, simp)

    fix j assume "j \<in> {0 ..< n}"

    show "1 / real (f (j' + j)) * Ifloat x ^ j = -1 ^ j * (1 / real (f (j' + j))) * Ifloat (- x) ^ j"

      unfolding move_minus power_mult_distrib real_mult_assoc[symmetric]

      unfolding real_mult_commute unfolding real_mult_assoc[of "-1 ^ j", symmetric] power_mult_distrib[symmetric]

      by auto

  qed

  have "0 \<le> Ifloat (-x)" using assms by auto

  from horner_bounds[where G=G and F=F and f=f and s=s and prec=prec

    and lb="\<lambda> n i k x. lb n i k (-x)" and ub="\<lambda> n i k x. ub n i k (-x)", unfolded lb_Suc ub_Suc diff_mult_minus,

    OF this f_Suc lb_0 refl ub_0 refl]

  show "?lb" and "?ub" unfolding minus_minus sum_eq

    by auto

qed

subsection {* 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 \<Rightarrow> nat" where

  "get_odd n = (if odd n then n else (Suc n))"

definition get_even :: "nat \<Rightarrow> 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_ex: "\<exists> k. Suc k = get_odd n \<and> odd (Suc k)"

proof (cases "odd n")

  case True hence "0 < n" by (rule odd_pos)

  from gr0_implies_Suc[OF this] obtain k where "Suc k = n" by auto 

  thus ?thesis unfolding get_odd_def if_P[OF True] using True[unfolded `Suc k = n`[symmetric]] by blast

next

  case False hence "odd (Suc n)" by auto

  thus ?thesis unfolding get_odd_def if_not_P[OF False] by blast

qed

lemma get_even_double: "\<exists>i. get_even n = 2 * i" using get_even[unfolded even_mult_two_ex] .

lemma get_odd_double: "\<exists>i. get_odd n = 2 * i + 1" using get_odd[unfolded odd_Suc_mult_two_ex] by auto

section "Power function"

definition float_power_bnds :: "nat \<Rightarrow> float \<Rightarrow> float \<Rightarrow> float * float" where

"float_power_bnds n l u = (if odd n \<or> 0 < l then (l ^ n, u ^ n)

                      else if u < 0         then (u ^ n, l ^ n)

                                            else (0, (max (-l) u) ^ n))"

lemma float_power_bnds: assumes "(l1, u1) = float_power_bnds n l u" and "x \<in> {Ifloat l .. Ifloat u}"

  shows "x ^ n \<in> {Ifloat l1..Ifloat u1}"

proof (cases "even n")

  case True 

  show ?thesis

  proof (cases "0 < l")

    case True hence "odd n \<or> 0 < l" and "0 \<le> Ifloat l" unfolding less_float_def by auto

    have u1: "u1 = u ^ n" and l1: "l1 = l ^ n" using assms unfolding float_power_bnds_def if_P[OF `odd n \<or> 0 < l`] by auto

    have "Ifloat l ^ n \<le> x ^ n" and "x ^ n \<le> Ifloat u ^ n " using `0 \<le> Ifloat l` and assms unfolding atLeastAtMost_iff using power_mono[of "Ifloat l" x] power_mono[of x "Ifloat u"] by auto

    thus ?thesis using assms `0 < l` unfolding atLeastAtMost_iff l1 u1 float_power less_float_def by auto

  next

    case False hence P: "\<not> (odd n \<or> 0 < l)" using `even n` by auto

    show ?thesis

    proof (cases "u < 0")

      case True hence "0 \<le> - Ifloat u" and "- Ifloat u \<le> - x" and "0 \<le> - x" and "-x \<le> - Ifloat l" using assms unfolding less_float_def by auto

      hence "Ifloat u ^ n \<le> x ^ n" and "x ^ n \<le> Ifloat l ^ n" using power_mono[of  "-x" "-Ifloat l" n] power_mono[of "-Ifloat u" "-x" n] 

	unfolding power_minus_even[OF `even n`] by auto

      moreover have u1: "u1 = l ^ n" and l1: "l1 = u ^ n" using assms unfolding float_power_bnds_def if_not_P[OF P] if_P[OF True] by auto

      ultimately show ?thesis using float_power by auto

    next

      case False 

      have "\<bar>x\<bar> \<le> Ifloat (max (-l) u)"

      proof (cases "-l \<le> u")

	case True thus ?thesis unfolding max_def if_P[OF True] using assms unfolding le_float_def by auto

      next

	case False thus ?thesis unfolding max_def if_not_P[OF False] using assms unfolding le_float_def by auto

      qed

      hence x_abs: "\<bar>x\<bar> \<le> \<bar>Ifloat (max (-l) u)\<bar>" by auto

      have u1: "u1 = (max (-l) u) ^ n" and l1: "l1 = 0" using assms unfolding float_power_bnds_def if_not_P[OF P] if_not_P[OF False] by auto

      show ?thesis unfolding atLeastAtMost_iff l1 u1 float_power using zero_le_even_power[OF `even n`] power_mono_even[OF `even n` x_abs] by auto

    qed

  qed

next

  case False hence "odd n \<or> 0 < l" by auto

  have u1: "u1 = u ^ n" and l1: "l1 = l ^ n" using assms unfolding float_power_bnds_def if_P[OF `odd n \<or> 0 < l`] by auto

  have "Ifloat l ^ n \<le> x ^ n" and "x ^ n \<le> Ifloat u ^ n " using assms unfolding atLeastAtMost_iff using power_mono_odd[OF False] by auto

  thus ?thesis unfolding atLeastAtMost_iff l1 u1 float_power less_float_def by auto

qed

lemma bnds_power: "\<forall> x l u. (l1, u1) = float_power_bnds n l u \<and> x \<in> {Ifloat l .. Ifloat u} \<longrightarrow> Ifloat l1 \<le> x ^ n \<and> x ^ n \<le> Ifloat 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 \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float" where

"sqrt_iteration prec 0 (Float m e) = Float 1 ((e + bitlen m) div 2 + 1)" |

"sqrt_iteration prec (Suc m) x = (let y = sqrt_iteration prec m x 

                                  in Float 1 -1 * (y + float_divr prec x y))"

definition ub_sqrt :: "nat \<Rightarrow> float \<Rightarrow> float option" where 

"ub_sqrt prec x = (if 0 < x then Some (sqrt_iteration prec prec x) else if x < 0 then None else Some 0)"

definition lb_sqrt :: "nat \<Rightarrow> float \<Rightarrow> float option" where

"lb_sqrt prec x = (if 0 < x then Some (float_divl prec x (sqrt_iteration prec prec x)) else if x < 0 then None else Some 0)"

lemma sqrt_ub_pos_pos_1:

  assumes "sqrt x < b" and "0 < b" and "0 < x"

  shows "sqrt x < (b + x / b)/2"

proof -

  from assms have "0 < (b - sqrt x) ^ 2 " by simp

  also have "\<dots> = b ^ 2 - 2 * b * sqrt x + (sqrt x) ^ 2" by algebra

  also have "\<dots> = b ^ 2 - 2 * b * sqrt x + x" using assms by (simp add: real_sqrt_pow2)

  finally have "0 < b ^ 2 - 2 * b * sqrt x + x" by assumption

  hence "0 < b / 2 - sqrt x + x / (2 * b)" using assms

    by (simp add: field_simps power2_eq_square)

  thus ?thesis by (simp add: field_simps)

qed

lemma sqrt_iteration_bound: assumes "0 < Ifloat x"

  shows "sqrt (Ifloat x) < Ifloat (sqrt_iteration prec n x)"

proof (induct n)

  case 0

  show ?case

  proof (cases x)

    case (Float m e)

    hence "0 < m" using float_pos_m_pos[unfolded less_float_def] assms by auto

    hence "0 < sqrt (real m)" by auto

    have int_nat_bl: "int (nat (bitlen m)) = bitlen m" using bitlen_ge0 by auto

    have "Ifloat x = (real m / 2^nat (bitlen m)) * pow2 (e + int (nat (bitlen m)))"

      unfolding pow2_add pow2_int Float Ifloat.simps by auto

    also have "\<dots> < 1 * pow2 (e + int (nat (bitlen m)))"

    proof (rule mult_strict_right_mono, auto)

      show "real m < 2^nat (bitlen m)" using bitlen_bounds[OF `0 < m`, THEN conjunct2] 

	unfolding real_of_int_less_iff[of m, symmetric] by auto

    qed

    finally have "sqrt (Ifloat x) < sqrt (pow2 (e + bitlen m))" unfolding int_nat_bl by auto

    also have "\<dots> \<le> pow2 ((e + bitlen m) div 2 + 1)"

    proof -

      let ?E = "e + bitlen m"

      have E_mod_pow: "pow2 (?E mod 2) < 4"

      proof (cases "?E mod 2 = 1")

	case True thus ?thesis by auto

      next

	case False 

	have "0 \<le> ?E mod 2" by auto 

	have "?E mod 2 < 2" by auto

	from this[THEN zless_imp_add1_zle]

	have "?E mod 2 \<le> 0" using False by auto

	from xt1(5)[OF `0 \<le> ?E mod 2` this]

	show ?thesis by auto

      qed

      hence "sqrt (pow2 (?E mod 2)) < sqrt (2 * 2)" by auto

      hence E_mod_pow: "sqrt (pow2 (?E mod 2)) < 2" unfolding real_sqrt_abs2 by auto

      have E_eq: "pow2 ?E = pow2 (?E div 2 + ?E div 2 + ?E mod 2)" by auto

      have "sqrt (pow2 ?E) = sqrt (pow2 (?E div 2) * pow2 (?E div 2) * pow2 (?E mod 2))"

	unfolding E_eq unfolding pow2_add ..

      also have "\<dots> = pow2 (?E div 2) * sqrt (pow2 (?E mod 2))"

	unfolding real_sqrt_mult[of _ "pow2 (?E mod 2)"] real_sqrt_abs2 by auto

      also have "\<dots> < pow2 (?E div 2) * 2" 

	by (rule mult_strict_left_mono, auto intro: E_mod_pow)

      also have "\<dots> = pow2 (?E div 2 + 1)" unfolding zadd_commute[of _ 1] pow2_add1 by auto

      finally show ?thesis by auto

    qed

    finally show ?thesis 

      unfolding Float sqrt_iteration.simps Ifloat.simps by auto

  qed

next

  case (Suc n)

  let ?b = "sqrt_iteration prec n x"

  have "0 < sqrt (Ifloat x)" using `0 < Ifloat x` by auto

  also have "\<dots> < Ifloat ?b" using Suc .

  finally have "sqrt (Ifloat x) < (Ifloat ?b + Ifloat x / Ifloat ?b)/2" using sqrt_ub_pos_pos_1[OF Suc _ `0 < Ifloat x`] by auto

  also have "\<dots> \<le> (Ifloat ?b + Ifloat (float_divr prec x ?b))/2" by (rule divide_right_mono, auto simp add: float_divr)

  also have "\<dots> = Ifloat (Float 1 -1) * (Ifloat ?b + Ifloat (float_divr prec x ?b))" by auto

  finally show ?case unfolding sqrt_iteration.simps Let_def Ifloat_mult Ifloat_add right_distrib .

qed

lemma sqrt_iteration_lower_bound: assumes "0 < Ifloat x"

  shows "0 < Ifloat (sqrt_iteration prec n x)" (is "0 < ?sqrt")

proof -

  have "0 < sqrt (Ifloat x)" using assms by auto

  also have "\<dots> < ?sqrt" using sqrt_iteration_bound[OF assms] .

  finally show ?thesis .

qed

lemma lb_sqrt_lower_bound: assumes "0 \<le> Ifloat x"

  shows "0 \<le> Ifloat (the (lb_sqrt prec x))"

proof (cases "0 < x")

  case True hence "0 < Ifloat x" and "0 \<le> x" using `0 \<le> Ifloat x` unfolding less_float_def le_float_def by auto

  hence "0 < sqrt_iteration prec prec x" unfolding less_float_def using sqrt_iteration_lower_bound by auto 

  hence "0 \<le> Ifloat (float_divl prec x (sqrt_iteration prec prec x))" using float_divl_lower_bound[OF `0 \<le> x`] unfolding le_float_def by auto

  thus ?thesis unfolding lb_sqrt_def using True by auto

next

  case False with `0 \<le> Ifloat x` have "Ifloat x = 0" unfolding less_float_def by auto

  thus ?thesis unfolding lb_sqrt_def less_float_def by auto

qed

lemma lb_sqrt_upper_bound: assumes "0 \<le> Ifloat x"

  shows "Ifloat (the (lb_sqrt prec x)) \<le> sqrt (Ifloat x)"

proof (cases "0 < x")

  case True hence "0 < Ifloat x" and "0 \<le> Ifloat x" unfolding less_float_def by auto

  hence sqrt_gt0: "0 < sqrt (Ifloat x)" by auto

  hence sqrt_ub: "sqrt (Ifloat x) < Ifloat (sqrt_iteration prec prec x)" using sqrt_iteration_bound by auto

  have "Ifloat (float_divl prec x (sqrt_iteration prec prec x)) \<le> Ifloat x / Ifloat (sqrt_iteration prec prec x)" by (rule float_divl)

  also have "\<dots> < Ifloat x / sqrt (Ifloat x)" 

    by (rule divide_strict_left_mono[OF sqrt_ub `0 < Ifloat x` mult_pos_pos[OF order_less_trans[OF sqrt_gt0 sqrt_ub] sqrt_gt0]])

  also have "\<dots> = sqrt (Ifloat x)" unfolding inverse_eq_iff_eq[of _ "sqrt (Ifloat x)", symmetric] sqrt_divide_self_eq[OF `0 \<le> Ifloat x`, symmetric] by auto

  finally show ?thesis unfolding lb_sqrt_def if_P[OF `0 < x`] by auto

next

  case False with `0 \<le> Ifloat x`

  have "\<not> x < 0" unfolding less_float_def le_float_def by auto

  show ?thesis unfolding lb_sqrt_def if_not_P[OF False] if_not_P[OF `\<not> x < 0`] using assms by auto

qed

lemma lb_sqrt: assumes "Some y = lb_sqrt prec x"

  shows "Ifloat y \<le> sqrt (Ifloat x)" and "0 \<le> Ifloat x"

proof -

  show "0 \<le> Ifloat x"

  proof (rule ccontr)

    assume "\<not> 0 \<le> Ifloat x"

    hence "lb_sqrt prec x = None" unfolding lb_sqrt_def less_float_def by auto

    thus False using assms by auto

  qed

  from lb_sqrt_upper_bound[OF this, of prec]

  show "Ifloat y \<le> sqrt (Ifloat x)" unfolding assms[symmetric] by auto

qed

lemma ub_sqrt_lower_bound: assumes "0 \<le> Ifloat x"

  shows "sqrt (Ifloat x) \<le> Ifloat (the (ub_sqrt prec x))"

proof (cases "0 < x")

  case True hence "0 < Ifloat x" unfolding less_float_def by auto

  hence "0 < sqrt (Ifloat x)" by auto

  hence "sqrt (Ifloat x) < Ifloat (sqrt_iteration prec prec x)" using sqrt_iteration_bound by auto

  thus ?thesis unfolding ub_sqrt_def if_P[OF `0 < x`] by auto

next

  case False with `0 \<le> Ifloat x`

  have "Ifloat x = 0" unfolding less_float_def le_float_def by auto

  thus ?thesis unfolding ub_sqrt_def less_float_def le_float_def by auto

qed

lemma ub_sqrt: assumes "Some y = ub_sqrt prec x"

  shows "sqrt (Ifloat x) \<le> Ifloat y" and "0 \<le> Ifloat x"

proof -

  show "0 \<le> Ifloat x"

  proof (rule ccontr)

    assume "\<not> 0 \<le> Ifloat x"

    hence "ub_sqrt prec x = None" unfolding ub_sqrt_def less_float_def by auto

    thus False using assms by auto

  qed

  from ub_sqrt_lower_bound[OF this, of prec]

  show "sqrt (Ifloat x) \<le> Ifloat y" unfolding assms[symmetric] by auto

qed

lemma bnds_sqrt: "\<forall> x lx ux. (Some l, Some u) = (lb_sqrt prec lx, ub_sqrt prec ux) \<and> x \<in> {Ifloat lx .. Ifloat ux} \<longrightarrow> Ifloat l \<le> sqrt x \<and> sqrt x \<le> Ifloat u"

proof (rule allI, rule allI, rule allI, rule impI)

  fix x lx ux

  assume "(Some l, Some u) = (lb_sqrt prec lx, ub_sqrt prec ux) \<and> x \<in> {Ifloat lx .. Ifloat ux}"

  hence l: "Some l = lb_sqrt prec lx " and u: "Some u = ub_sqrt prec ux" and x: "x \<in> {Ifloat lx .. Ifloat ux}" by auto

  have "Ifloat lx \<le> x" and "x \<le> Ifloat ux" using x by auto

  from lb_sqrt(1)[OF l] real_sqrt_le_mono[OF `Ifloat lx \<le> x`]

  have "Ifloat l \<le> sqrt x" by (rule order_trans)

  moreover

  from real_sqrt_le_mono[OF `x \<le> Ifloat ux`] ub_sqrt(1)[OF u]

  have "sqrt x \<le> Ifloat u" by (rule order_trans)

  ultimately show "Ifloat l \<le> sqrt x \<and> sqrt x \<le> Ifloat u" ..

qed

section "Arcus tangens and \<pi>"

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 \<pi> and then the entire arcus tangens.

*}

fun ub_arctan_horner :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float"

and lb_arctan_horner :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float" where

  "ub_arctan_horner prec 0 k x = 0"

| "ub_arctan_horner prec (Suc n) k x = 

    (rapprox_rat prec 1 (int k)) - x * (lb_arctan_horner prec n (k + 2) x)"

| "lb_arctan_horner prec 0 k x = 0"

| "lb_arctan_horner prec (Suc n) k x = 

    (lapprox_rat prec 1 (int k)) - x * (ub_arctan_horner prec n (k + 2) x)"

lemma arctan_0_1_bounds': assumes "0 \<le> Ifloat x" "Ifloat x \<le> 1" and "even n"

  shows "arctan (Ifloat x) \<in> {Ifloat (x * lb_arctan_horner prec n 1 (x * x)) .. Ifloat (x * ub_arctan_horner prec (Suc n) 1 (x * x))}"

proof -

  let "?c i" = "-1^i * (1 / real (i * 2 + 1) * Ifloat x ^ (i * 2 + 1))"

  let "?S n" = "\<Sum> i=0..<n. ?c i"

  have "0 \<le> Ifloat (x * x)" by auto

  from `even n` obtain m where "2 * m = n" unfolding even_mult_two_ex by auto

  have "arctan (Ifloat x) \<in> { ?S n .. ?S (Suc n) }"

  proof (cases "Ifloat x = 0")

    case False

    hence "0 < Ifloat x" using `0 \<le> Ifloat x` by auto

    hence prem: "0 < 1 / real (0 * 2 + (1::nat)) * Ifloat x ^ (0 * 2 + 1)" by auto 

    have "\<bar> Ifloat x \<bar> \<le> 1"  using `0 \<le> Ifloat x` `Ifloat x \<le> 1` by auto

    from mp[OF summable_Leibniz(2)[OF zeroseq_arctan_series[OF this] monoseq_arctan_series[OF this]] prem, THEN spec, of m, unfolded `2 * m = n`]

    show ?thesis unfolding arctan_series[OF `\<bar> Ifloat x \<bar> \<le> 1`] Suc_plus1  .

  qed auto

  note arctan_bounds = this[unfolded atLeastAtMost_iff]

  have F: "\<And>n. 2 * Suc n + 1 = 2 * n + 1 + 2" by auto

  note bounds = horner_bounds[where s=1 and f="\<lambda>i. 2 * i + 1" and j'=0 

    and lb="\<lambda>n i k x. lb_arctan_horner prec n k x"

    and ub="\<lambda>n i k x. ub_arctan_horner prec n k x", 

    OF `0 \<le> Ifloat (x*x)` F lb_arctan_horner.simps ub_arctan_horner.simps]

  { have "Ifloat (x * lb_arctan_horner prec n 1 (x*x)) \<le> ?S n"

      using bounds(1) `0 \<le> Ifloat x`

      unfolding Ifloat_mult power_add power_one_right real_mult_assoc[symmetric] setsum_left_distrib[symmetric]

      unfolding real_mult_commute mult_commute[of _ "2::nat"] power_mult power2_eq_square[of "Ifloat x"]

      by (auto intro!: mult_left_mono)

    also have "\<dots> \<le> arctan (Ifloat x)" using arctan_bounds ..

    finally have "Ifloat (x * lb_arctan_horner prec n 1 (x*x)) \<le> arctan (Ifloat x)" . }

  moreover

  { have "arctan (Ifloat x) \<le> ?S (Suc n)" using arctan_bounds ..

    also have "\<dots> \<le> Ifloat (x * ub_arctan_horner prec (Suc n) 1 (x*x))"

      using bounds(2)[of "Suc n"] `0 \<le> Ifloat x`

      unfolding Ifloat_mult power_add power_one_right real_mult_assoc[symmetric] setsum_left_distrib[symmetric]

      unfolding real_mult_commute mult_commute[of _ "2::nat"] power_mult power2_eq_square[of "Ifloat x"]

      by (auto intro!: mult_left_mono)

    finally have "arctan (Ifloat x) \<le> Ifloat (x * ub_arctan_horner prec (Suc n) 1 (x*x))" . }

  ultimately show ?thesis by auto

qed

lemma arctan_0_1_bounds: assumes "0 \<le> Ifloat x" "Ifloat x \<le> 1"

  shows "arctan (Ifloat x) \<in> {Ifloat (x * lb_arctan_horner prec (get_even n) 1 (x * x)) .. Ifloat (x * ub_arctan_horner prec (get_odd n) 1 (x * x))}"

proof (cases "even n")

  case True

  obtain n' where "Suc n' = get_odd n" and "odd (Suc n')" using get_odd_ex by auto

  hence "even n'" unfolding even_nat_Suc by auto

  have "arctan (Ifloat x) \<le> Ifloat (x * ub_arctan_horner prec (get_odd n) 1 (x * x))"

    unfolding `Suc n' = get_odd n`[symmetric] using arctan_0_1_bounds'[OF `0 \<le> Ifloat x` `Ifloat x \<le> 1` `even n'`] by auto

  moreover

  have "Ifloat (x * lb_arctan_horner prec (get_even n) 1 (x * x)) \<le> arctan (Ifloat x)"

    unfolding get_even_def if_P[OF True] using arctan_0_1_bounds'[OF `0 \<le> Ifloat x` `Ifloat x \<le> 1` `even n`] by auto

  ultimately show ?thesis by auto

next

  case False hence "0 < n" by (rule odd_pos)

  from gr0_implies_Suc[OF this] obtain n' where "n = Suc n'" ..

  from False[unfolded this even_nat_Suc]

  have "even n'" and "even (Suc (Suc n'))" by auto

  have "get_odd n = Suc n'" unfolding get_odd_def if_P[OF False] using `n = Suc n'` .

  have "arctan (Ifloat x) \<le> Ifloat (x * ub_arctan_horner prec (get_odd n) 1 (x * x))"

    unfolding `get_odd n = Suc n'` using arctan_0_1_bounds'[OF `0 \<le> Ifloat x` `Ifloat x \<le> 1` `even n'`] by auto

  moreover

  have "Ifloat (x * lb_arctan_horner prec (get_even n) 1 (x * x)) \<le> arctan (Ifloat x)"

    unfolding get_even_def if_not_P[OF False] unfolding `n = Suc n'` using arctan_0_1_bounds'[OF `0 \<le> Ifloat x` `Ifloat x \<le> 1` `even (Suc (Suc n'))`] by auto

  ultimately show ?thesis by auto

qed

subsection "Compute \<pi>"

definition ub_pi :: "nat \<Rightarrow> float" where

  "ub_pi prec = (let A = rapprox_rat prec 1 5 ; 

                     B = lapprox_rat prec 1 239

                 in ((Float 1 2) * ((Float 1 2) * A * (ub_arctan_horner prec (get_odd (prec div 4 + 1)) 1 (A * A)) - 

                                                  B * (lb_arctan_horner prec (get_even (prec div 14 + 1)) 1 (B * B)))))"

definition lb_pi :: "nat \<Rightarrow> float" where

  "lb_pi prec = (let A = lapprox_rat prec 1 5 ; 

                     B = rapprox_rat prec 1 239

                 in ((Float 1 2) * ((Float 1 2) * A * (lb_arctan_horner prec (get_even (prec div 4 + 1)) 1 (A * A)) - 

                                                  B * (ub_arctan_horner prec (get_odd (prec div 14 + 1)) 1 (B * B)))))"

lemma pi_boundaries: "pi \<in> {Ifloat (lb_pi n) .. Ifloat (ub_pi n)}"

proof -

  have machin_pi: "pi = 4 * (4 * arctan (1 / 5) - arctan (1 / 239))" unfolding machin[symmetric] by auto

  { fix prec n :: nat fix k :: int assume "1 < k" hence "0 \<le> k" and "0 < k" and "1 \<le> k" by auto

    let ?k = "rapprox_rat prec 1 k"

    have "1 div k = 0" using div_pos_pos_trivial[OF _ `1 < k`] by auto

    have "0 \<le> Ifloat ?k" by (rule order_trans[OF _ rapprox_rat], auto simp add: `0 \<le> k`)

    have "Ifloat ?k \<le> 1" unfolding rapprox_rat.simps(2)[OF zero_le_one `0 < k`]

      by (rule rapprox_posrat_le1, auto simp add: `0 < k` `1 \<le> k`)

    have "1 / real k \<le> Ifloat ?k" using rapprox_rat[where x=1 and y=k] by auto

    hence "arctan (1 / real k) \<le> arctan (Ifloat ?k)" by (rule arctan_monotone')

    also have "\<dots> \<le> Ifloat (?k * ub_arctan_horner prec (get_odd n) 1 (?k * ?k))"

      using arctan_0_1_bounds[OF `0 \<le> Ifloat ?k` `Ifloat ?k \<le> 1`] by auto

    finally have "arctan (1 / (real k)) \<le> Ifloat (?k * ub_arctan_horner prec (get_odd n) 1 (?k * ?k))" .

  } note ub_arctan = this

  { fix prec n :: nat fix k :: int assume "1 < k" hence "0 \<le> k" and "0 < k" by auto

    let ?k = "lapprox_rat prec 1 k"

    have "1 div k = 0" using div_pos_pos_trivial[OF _ `1 < k`] by auto

    have "1 / real k \<le> 1" using `1 < k` by auto

    have "\<And>n. 0 \<le> Ifloat ?k" using lapprox_rat_bottom[where x=1 and y=k, OF zero_le_one `0 < k`] by (auto simp add: `1 div k = 0`)

    have "\<And>n. Ifloat ?k \<le> 1" using lapprox_rat by (rule order_trans, auto simp add: `1 / real k \<le> 1`)

    have "Ifloat ?k \<le> 1 / real k" using lapprox_rat[where x=1 and y=k] by auto

    have "Ifloat (?k * lb_arctan_horner prec (get_even n) 1 (?k * ?k)) \<le> arctan (Ifloat ?k)"

      using arctan_0_1_bounds[OF `0 \<le> Ifloat ?k` `Ifloat ?k \<le> 1`] by auto

    also have "\<dots> \<le> arctan (1 / real k)" using `Ifloat ?k \<le> 1 / real k` by (rule arctan_monotone')

    finally have "Ifloat (?k * lb_arctan_horner prec (get_even n) 1 (?k * ?k)) \<le> arctan (1 / (real k))" .

  } note lb_arctan = this

  have "pi \<le> Ifloat (ub_pi n)"

    unfolding ub_pi_def machin_pi Let_def Ifloat_mult Ifloat_sub unfolding Float_num

    using lb_arctan[of 239] ub_arctan[of 5]

    by (auto intro!: mult_left_mono add_mono simp add: diff_minus simp del: lapprox_rat.simps rapprox_rat.simps)

  moreover

  have "Ifloat (lb_pi n) \<le> pi"

    unfolding lb_pi_def machin_pi Let_def Ifloat_mult Ifloat_sub Float_num

    using lb_arctan[of 5] ub_arctan[of 239]

    by (auto intro!: mult_left_mono add_mono simp add: diff_minus simp del: lapprox_rat.simps rapprox_rat.simps)

  ultimately show ?thesis by auto

qed

subsection "Compute arcus tangens in the entire domain"

function lb_arctan :: "nat \<Rightarrow> float \<Rightarrow> float" and ub_arctan :: "nat \<Rightarrow> float \<Rightarrow> float" where 

  "lb_arctan prec x = (let ub_horner = \<lambda> x. x * ub_arctan_horner prec (get_odd (prec div 4 + 1)) 1 (x * x) ;

                           lb_horner = \<lambda> x. x * lb_arctan_horner prec (get_even (prec div 4 + 1)) 1 (x * x)

    in (if x < 0          then - ub_arctan prec (-x) else

        if x \<le> Float 1 -1 then lb_horner x else

        if x \<le> Float 1 1  then Float 1 1 * lb_horner (float_divl prec x (1 + the (ub_sqrt prec (1 + x * x))))

                          else (let inv = float_divr prec 1 x 

                                in if inv > 1 then 0 

                                              else lb_pi prec * Float 1 -1 - ub_horner inv)))"

| "ub_arctan prec x = (let lb_horner = \<lambda> x. x * lb_arctan_horner prec (get_even (prec div 4 + 1)) 1 (x * x) ;

                           ub_horner = \<lambda> x. x * ub_arctan_horner prec (get_odd (prec div 4 + 1)) 1 (x * x)

    in (if x < 0          then - lb_arctan prec (-x) else

        if x \<le> Float 1 -1 then ub_horner x else

        if x \<le> Float 1 1  then let y = float_divr prec x (1 + the (lb_sqrt prec (1 + x * x)))

                               in if y > 1 then ub_pi prec * Float 1 -1 

                                           else Float 1 1 * ub_horner y 

                          else ub_pi prec * Float 1 -1 - lb_horner (float_divl prec 1 x)))"

by pat_completeness auto

termination by (relation "measure (\<lambda> v. let (prec, x) = sum_case id id v in (if x < 0 then 1 else 0))", auto simp add: less_float_def)

declare ub_arctan_horner.simps[simp del]

declare lb_arctan_horner.simps[simp del]

lemma lb_arctan_bound': assumes "0 \<le> Ifloat x"

  shows "Ifloat (lb_arctan prec x) \<le> arctan (Ifloat x)"

proof -

  have "\<not> x < 0" and "0 \<le> x" unfolding less_float_def le_float_def using `0 \<le> Ifloat x` by auto

  let "?ub_horner x" = "x * ub_arctan_horner prec (get_odd (prec div 4 + 1)) 1 (x * x)"

    and "?lb_horner x" = "x * lb_arctan_horner prec (get_even (prec div 4 + 1)) 1 (x * x)"

  show ?thesis

  proof (cases "x \<le> Float 1 -1")

    case True hence "Ifloat x \<le> 1" unfolding le_float_def Float_num by auto

    show ?thesis unfolding lb_arctan.simps Let_def if_not_P[OF `\<not> x < 0`] if_P[OF True]

      using arctan_0_1_bounds[OF `0 \<le> Ifloat x` `Ifloat x \<le> 1`] by auto

  next

    case False hence "0 < Ifloat x" unfolding le_float_def Float_num by auto

    let ?R = "1 + sqrt (1 + Ifloat x * Ifloat x)"

    let ?fR = "1 + the (ub_sqrt prec (1 + x * x))"

    let ?DIV = "float_divl prec x ?fR"

    have sqr_ge0: "0 \<le> 1 + Ifloat x * Ifloat x" using sum_power2_ge_zero[of 1 "Ifloat x", unfolded numeral_2_eq_2] by auto

    hence divisor_gt0: "0 < ?R" by (auto intro: add_pos_nonneg)

    have "sqrt (Ifloat (1 + x * x)) \<le> Ifloat (the (ub_sqrt prec (1 + x * x)))" by (rule ub_sqrt_lower_bound, auto simp add: sqr_ge0)

    hence "?R \<le> Ifloat ?fR" by auto

    hence "0 < ?fR" and "0 < Ifloat ?fR" unfolding less_float_def using `0 < ?R` by auto

    have monotone: "Ifloat (float_divl prec x ?fR) \<le> Ifloat x / ?R"

    proof -

      have "Ifloat ?DIV \<le> Ifloat x / Ifloat ?fR" by (rule float_divl)

      also have "\<dots> \<le> Ifloat x / ?R" by (rule divide_left_mono[OF `?R \<le> Ifloat ?fR` `0 \<le> Ifloat x` mult_pos_pos[OF order_less_le_trans[OF divisor_gt0 `?R \<le> Ifloat ?fR`] divisor_gt0]])

      finally show ?thesis .

    qed

    show ?thesis

    proof (cases "x \<le> Float 1 1")

      case True

      have "Ifloat x \<le> sqrt (Ifloat (1 + x * x))" using real_sqrt_sum_squares_ge2[where x=1, unfolded numeral_2_eq_2] by auto

      also have "\<dots> \<le> Ifloat (the (ub_sqrt prec (1 + x * x)))" by (rule ub_sqrt_lower_bound, auto simp add: sqr_ge0)

      finally have "Ifloat x \<le> Ifloat ?fR" by auto

      moreover have "Ifloat ?DIV \<le> Ifloat x / Ifloat ?fR" by (rule float_divl)

      ultimately have "Ifloat ?DIV \<le> 1" unfolding divide_le_eq_1_pos[OF `0 < Ifloat ?fR`, symmetric] by auto

      have "0 \<le> Ifloat ?DIV" using float_divl_lower_bound[OF `0 \<le> x` `0 < ?fR`] unfolding le_float_def by auto

      have "Ifloat (Float 1 1 * ?lb_horner ?DIV) \<le> 2 * arctan (Ifloat (float_divl prec x ?fR))" unfolding Ifloat_mult[of "Float 1 1"] Float_num

	using arctan_0_1_bounds[OF `0 \<le> Ifloat ?DIV` `Ifloat ?DIV \<le> 1`] by auto

      also have "\<dots> \<le> 2 * arctan (Ifloat x / ?R)"

	using arctan_monotone'[OF monotone] by (auto intro!: mult_left_mono)

      also have "2 * arctan (Ifloat x / ?R) = arctan (Ifloat x)" using arctan_half[symmetric] unfolding numeral_2_eq_2 power_Suc2 power_0 real_mult_1 . 

      finally show ?thesis unfolding lb_arctan.simps Let_def if_not_P[OF `\<not> x < 0`] if_not_P[OF `\<not> x \<le> Float 1 -1`] if_P[OF True] .

    next

      case False

      hence "2 < Ifloat x" unfolding le_float_def Float_num by auto

      hence "1 \<le> Ifloat x" by auto

      let "?invx" = "float_divr prec 1 x"

      have "0 \<le> arctan (Ifloat x)" using arctan_monotone'[OF `0 \<le> Ifloat 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 `\<not> x < 0`] if_not_P[OF `\<not> x \<le> Float 1 -1`] if_not_P[OF False] if_P[OF True] 

	  using `0 \<le> arctan (Ifloat x)` by auto

      next

	case False

	hence "Ifloat ?invx \<le> 1" unfolding less_float_def by auto

	have "0 \<le> Ifloat ?invx" by (rule order_trans[OF _ float_divr], auto simp add: `0 \<le> Ifloat x`)

	have "1 / Ifloat x \<noteq> 0" and "0 < 1 / Ifloat x" using `0 < Ifloat x` by auto

	have "arctan (1 / Ifloat x) \<le> arctan (Ifloat ?invx)" unfolding Ifloat_1[symmetric] by (rule arctan_monotone', rule float_divr)

	also have "\<dots> \<le> Ifloat (?ub_horner ?invx)" using arctan_0_1_bounds[OF `0 \<le> Ifloat ?invx` `Ifloat ?invx \<le> 1`] by auto

	finally have "pi / 2 - Ifloat (?ub_horner ?invx) \<le> arctan (Ifloat x)" 

	  using `0 \<le> arctan (Ifloat x)` arctan_inverse[OF `1 / Ifloat x \<noteq> 0`] 

	  unfolding real_sgn_pos[OF `0 < 1 / Ifloat x`] le_diff_eq by auto

	moreover

	have "Ifloat (lb_pi prec * Float 1 -1) \<le> pi / 2" unfolding Ifloat_mult Float_num times_divide_eq_right real_mult_1 using pi_boundaries by auto

	ultimately

	show ?thesis unfolding lb_arctan.simps Let_def if_not_P[OF `\<not> x < 0`] if_not_P[OF `\<not> x \<le> Float 1 -1`] if_not_P[OF `\<not> x \<le> Float 1 1`] if_not_P[OF False]

	  by auto

      qed

    qed

  qed

qed

lemma ub_arctan_bound': assumes "0 \<le> Ifloat x"

  shows "arctan (Ifloat x) \<le> Ifloat (ub_arctan prec x)"

proof -

  have "\<not> x < 0" and "0 \<le> x" unfolding less_float_def le_float_def using `0 \<le> Ifloat x` by auto

  let "?ub_horner x" = "x * ub_arctan_horner prec (get_odd (prec div 4 + 1)) 1 (x * x)"

    and "?lb_horner x" = "x * lb_arctan_horner prec (get_even (prec div 4 + 1)) 1 (x * x)"

  show ?thesis

  proof (cases "x \<le> Float 1 -1")

    case True hence "Ifloat x \<le> 1" unfolding le_float_def Float_num by auto

    show ?thesis unfolding ub_arctan.simps Let_def if_not_P[OF `\<not> x < 0`] if_P[OF True]

      using arctan_0_1_bounds[OF `0 \<le> Ifloat x` `Ifloat x \<le> 1`] by auto

  next

    case False hence "0 < Ifloat x" unfolding le_float_def Float_num by auto

    let ?R = "1 + sqrt (1 + Ifloat x * Ifloat x)"

    let ?fR = "1 + the (lb_sqrt prec (1 + x * x))"

    let ?DIV = "float_divr prec x ?fR"

    have sqr_ge0: "0 \<le> 1 + Ifloat x * Ifloat x" using sum_power2_ge_zero[of 1 "Ifloat x", unfolded numeral_2_eq_2] by auto

    hence "0 \<le> Ifloat (1 + x*x)" by auto

    hence divisor_gt0: "0 < ?R" by (auto intro: add_pos_nonneg)

    have "Ifloat (the (lb_sqrt prec (1 + x * x))) \<le> sqrt (Ifloat (1 + x * x))" by (rule lb_sqrt_upper_bound, auto simp add: sqr_ge0)

    hence "Ifloat ?fR \<le> ?R" by auto

    have "0 < Ifloat ?fR" unfolding Ifloat_add Ifloat_1 by (rule order_less_le_trans[OF zero_less_one], auto simp add: lb_sqrt_lower_bound[OF `0 \<le> Ifloat (1 + x*x)`])

    have monotone: "Ifloat x / ?R \<le> Ifloat (float_divr prec x ?fR)"

    proof -

      from divide_left_mono[OF `Ifloat ?fR \<le> ?R` `0 \<le> Ifloat x` mult_pos_pos[OF divisor_gt0 `0 < Ifloat ?fR`]]

      have "Ifloat x / ?R \<le> Ifloat x / Ifloat ?fR" .

      also have "\<dots> \<le> Ifloat ?DIV" by (rule float_divr)

      finally show ?thesis .

    qed

    show ?thesis

    proof (cases "x \<le> Float 1 1")

      case True

      show ?thesis

      proof (cases "?DIV > 1")

	case True

	have "pi / 2 \<le> Ifloat (ub_pi prec * Float 1 -1)" unfolding Ifloat_mult Float_num times_divide_eq_right real_mult_1 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 `\<not> x < 0`] if_not_P[OF `\<not> x \<le> Float 1 -1`] if_P[OF `x \<le> Float 1 1`] if_P[OF True] .

      next

	case False

	hence "Ifloat ?DIV \<le> 1" unfolding less_float_def by auto

	have "0 \<le> Ifloat x / ?R" using `0 \<le> Ifloat x` `0 < ?R` unfolding real_0_le_divide_iff by auto

	hence "0 \<le> Ifloat ?DIV" using monotone by (rule order_trans)

	have "arctan (Ifloat x) = 2 * arctan (Ifloat x / ?R)" using arctan_half unfolding numeral_2_eq_2 power_Suc2 power_0 real_mult_1 .

	also have "\<dots> \<le> 2 * arctan (Ifloat ?DIV)"

	  using arctan_monotone'[OF monotone] by (auto intro!: mult_left_mono)

	also have "\<dots> \<le> Ifloat (Float 1 1 * ?ub_horner ?DIV)" unfolding Ifloat_mult[of "Float 1 1"] Float_num

	  using arctan_0_1_bounds[OF `0 \<le> Ifloat ?DIV` `Ifloat ?DIV \<le> 1`] by auto

	finally show ?thesis unfolding ub_arctan.simps Let_def if_not_P[OF `\<not> x < 0`] if_not_P[OF `\<not> x \<le> Float 1 -1`] if_P[OF `x \<le> Float 1 1`] if_not_P[OF False] .

      qed

    next

      case False

      hence "2 < Ifloat x" unfolding le_float_def Float_num by auto

      hence "1 \<le> Ifloat x" by auto

      hence "0 < Ifloat x" by auto

      hence "0 < x" unfolding less_float_def by auto

      let "?invx" = "float_divl prec 1 x"

      have "0 \<le> arctan (Ifloat x)" using arctan_monotone'[OF `0 \<le> Ifloat x`] using arctan_tan[of 0, unfolded tan_zero] by auto

      have "Ifloat ?invx \<le> 1" unfolding less_float_def by (rule order_trans[OF float_divl], auto simp add: `1 \<le> Ifloat x` divide_le_eq_1_pos[OF `0 < Ifloat x`])

      have "0 \<le> Ifloat ?invx" unfolding Ifloat_0[symmetric] by (rule float_divl_lower_bound[unfolded le_float_def], auto simp add: `0 < x`)

      have "1 / Ifloat x \<noteq> 0" and "0 < 1 / Ifloat x" using `0 < Ifloat x` by auto

      have "Ifloat (?lb_horner ?invx) \<le> arctan (Ifloat ?invx)" using arctan_0_1_bounds[OF `0 \<le> Ifloat ?invx` `Ifloat ?invx \<le> 1`] by auto

      also have "\<dots> \<le> arctan (1 / Ifloat x)" unfolding Ifloat_1[symmetric] by (rule arctan_monotone', rule float_divl)

      finally have "arctan (Ifloat x) \<le> pi / 2 - Ifloat (?lb_horner ?invx)"

	using `0 \<le> arctan (Ifloat x)` arctan_inverse[OF `1 / Ifloat x \<noteq> 0`] 

	unfolding real_sgn_pos[OF `0 < 1 / Ifloat x`] le_diff_eq by auto

      moreover

      have "pi / 2 \<le> Ifloat (ub_pi prec * Float 1 -1)" unfolding Ifloat_mult Float_num times_divide_eq_right mult_1_right using pi_boundaries by auto

      ultimately

      show ?thesis unfolding ub_arctan.simps Let_def if_not_P[OF `\<not> x < 0`] if_not_P[OF `\<not> x \<le> Float 1 -1`] if_not_P[OF `\<not> x \<le> Float 1 1`] if_not_P[OF False]

	by auto

    qed

  qed

qed

lemma arctan_boundaries:

  "arctan (Ifloat x) \<in> {Ifloat (lb_arctan prec x) .. Ifloat (ub_arctan prec x)}"

proof (cases "0 \<le> x")

  case True hence "0 \<le> Ifloat x" unfolding le_float_def by auto

  show ?thesis using ub_arctan_bound'[OF `0 \<le> Ifloat x`] lb_arctan_bound'[OF `0 \<le> Ifloat x`] unfolding atLeastAtMost_iff by auto

next

  let ?mx = "-x"

  case False hence "x < 0" and "0 \<le> Ifloat ?mx" unfolding le_float_def less_float_def by auto

  hence bounds: "Ifloat (lb_arctan prec ?mx) \<le> arctan (Ifloat ?mx) \<and> arctan (Ifloat ?mx) \<le> Ifloat (ub_arctan prec ?mx)"

    using ub_arctan_bound'[OF `0 \<le> Ifloat ?mx`] lb_arctan_bound'[OF `0 \<le> Ifloat ?mx`] by auto

  show ?thesis unfolding Ifloat_minus 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 Ifloat_minus arctan_minus] by auto

qed

lemma bnds_arctan: "\<forall> x lx ux. (l, u) = (lb_arctan prec lx, ub_arctan prec ux) \<and> x \<in> {Ifloat lx .. Ifloat ux} \<longrightarrow> Ifloat l \<le> arctan x \<and> arctan x \<le> Ifloat u"

proof (rule allI, rule allI, rule allI, rule impI)

  fix x lx ux

  assume "(l, u) = (lb_arctan prec lx, ub_arctan prec ux) \<and> x \<in> {Ifloat lx .. Ifloat ux}"

  hence l: "lb_arctan prec lx = l " and u: "ub_arctan prec ux = u" and x: "x \<in> {Ifloat lx .. Ifloat ux}" by auto

  { from arctan_boundaries[of lx prec, unfolded l]

    have "Ifloat l \<le> arctan (Ifloat lx)" by (auto simp del: lb_arctan.simps)

    also have "\<dots> \<le> arctan x" using x by (auto intro: arctan_monotone')

    finally have "Ifloat l \<le> arctan x" .

  } moreover

  { have "arctan x \<le> arctan (Ifloat ux)" using x by (auto intro: arctan_monotone')

    also have "\<dots> \<le> Ifloat u" using arctan_boundaries[of ux prec, unfolded u] by (auto simp del: ub_arctan.simps)

    finally have "arctan x \<le> Ifloat u" .

  } ultimately show "Ifloat l \<le> arctan x \<and> arctan x \<le> Ifloat u" ..

qed

section "Sinus and Cosinus"

subsection "Compute the cosinus and sinus series"

fun ub_sin_cos_aux :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float"

and lb_sin_cos_aux :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float" where

  "ub_sin_cos_aux prec 0 i k x = 0"

| "ub_sin_cos_aux prec (Suc n) i k x = 

    (rapprox_rat prec 1 (int k)) - x * (lb_sin_cos_aux prec n (i + 2) (k * i * (i + 1)) x)"

| "lb_sin_cos_aux prec 0 i k x = 0"

| "lb_sin_cos_aux prec (Suc n) i k x = 

    (lapprox_rat prec 1 (int k)) - x * (ub_sin_cos_aux prec n (i + 2) (k * i * (i + 1)) x)"

lemma cos_aux:

  shows "Ifloat (lb_sin_cos_aux prec n 1 1 (x * x)) \<le> (\<Sum> i=0..<n. -1^i * (1/real (fact (2 * i))) * (Ifloat x)^(2 * i))" (is "?lb")

  and "(\<Sum> i=0..<n. -1^i * (1/real (fact (2 * i))) * (Ifloat x)^(2 * i)) \<le> Ifloat (ub_sin_cos_aux prec n 1 1 (x * x))" (is "?ub")

proof -

  have "0 \<le> Ifloat (x * x)" unfolding Ifloat_mult by auto

  let "?f n" = "fact (2 * n)"

  { fix n 

    have F: "\<And>m. ((\<lambda>i. i + 2) o^ n) m = m + 2 * n" by (induct n arbitrary: m, auto)

    have "?f (Suc n) = ?f n * ((\<lambda>i. i + 2) o^ n) 1 * (((\<lambda>i. i + 2) o^ n) 1 + 1)"

      unfolding F by auto } note f_eq = this

  from horner_bounds[where lb="lb_sin_cos_aux prec" and ub="ub_sin_cos_aux prec" and j'=0, 

    OF `0 \<le> Ifloat (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 "Ifloat x"])

qed

lemma cos_boundaries: assumes "0 \<le> Ifloat x" and "Ifloat x \<le> pi / 2"

  shows "cos (Ifloat x) \<in> {Ifloat (lb_sin_cos_aux prec (get_even n) 1 1 (x * x)) .. Ifloat (ub_sin_cos_aux prec (get_odd n) 1 1 (x * x))}"

proof (cases "Ifloat x = 0")

  case False hence "Ifloat x \<noteq> 0" by auto

  hence "0 < x" and "0 < Ifloat x" using `0 \<le> Ifloat x` unfolding less_float_def by auto

  have "0 < x * x" using `0 < x` unfolding less_float_def Ifloat_mult Ifloat_0

    using mult_pos_pos[where a="Ifloat x" and b="Ifloat x"] by auto

  { fix x n have "(\<Sum> i=0..<n. -1^i * (1/real (fact (2 * i))) * x ^ (2 * i))

    = (\<Sum> i = 0 ..< 2 * n. (if even(i) then (-1 ^ (i div 2))/(real (fact i)) else 0) * x ^ i)" (is "?sum = ?ifsum")

  proof -

    have "?sum = ?sum + (\<Sum> j = 0 ..< n. 0)" by auto

    also have "\<dots> = 

      (\<Sum> j = 0 ..< n. -1 ^ ((2 * j) div 2) / (real (fact (2 * j))) * x ^(2 * j)) + (\<Sum> j = 0 ..< n. 0)" by auto

    also have "\<dots> = (\<Sum> i = 0 ..< 2 * n. if even i then -1 ^ (i div 2) / (real (fact i)) * x ^ i else 0)"

      unfolding sum_split_even_odd ..

    also have "\<dots> = (\<Sum> i = 0 ..< 2 * n. (if even i then -1 ^ (i div 2) / (real (fact i)) else 0) * x ^ i)"

      by (rule setsum_cong2) auto

    finally show ?thesis by assumption

  qed } note morph_to_if_power = this

  { fix n :: nat assume "0 < n"

    hence "0 < 2 * n" by auto

    obtain t where "0 < t" and "t < Ifloat x" and

      cos_eq: "cos (Ifloat x) = (\<Sum> i = 0 ..< 2 * n. (if even(i) then (-1 ^ (i div 2))/(real (fact i)) else 0) * (Ifloat x) ^ i) 

      + (cos (t + 1/2 * real (2 * n) * pi) / real (fact (2*n))) * (Ifloat x)^(2*n)" 

      (is "_ = ?SUM + ?rest / ?fact * ?pow")

      using Maclaurin_cos_expansion2[OF `0 < Ifloat x` `0 < 2 * n`] by auto

    have "cos t * -1^n = cos t * cos (real n * pi) + sin t * sin (real n * pi)" by auto

    also have "\<dots> = cos (t + real n * pi)"  using cos_add by auto

    also have "\<dots> = ?rest" by auto

    finally have "cos t * -1^n = ?rest" .

    moreover

    have "t \<le> pi / 2" using `t < Ifloat x` and `Ifloat x \<le> pi / 2` by auto

    hence "0 \<le> cos t" using `0 < t` and cos_ge_zero by auto

    ultimately have even: "even n \<Longrightarrow> 0 \<le> ?rest" and odd: "odd n \<Longrightarrow> 0 \<le> - ?rest " by auto

    have "0 < ?fact" by auto

    have "0 < ?pow" using `0 < Ifloat x` by auto

    {

      assume "even n"

      have "Ifloat (lb_sin_cos_aux prec n 1 1 (x * x)) \<le> ?SUM"

	unfolding morph_to_if_power[symmetric] using cos_aux by auto 

      also have "\<dots> \<le> cos (Ifloat x)"

      proof -

	from even[OF `even n`] `0 < ?fact` `0 < ?pow`

	have "0 \<le> (?rest / ?fact) * ?pow" by (metis mult_nonneg_nonneg divide_nonneg_pos less_imp_le)

	thus ?thesis unfolding cos_eq by auto

      qed

      finally have "Ifloat (lb_sin_cos_aux prec n 1 1 (x * x)) \<le> cos (Ifloat x)" .

    } note lb = this

    {

      assume "odd n"

      have "cos (Ifloat x) \<le> ?SUM"

      proof -

	from `0 < ?fact` and `0 < ?pow` and odd[OF `odd n`]

	have "0 \<le> (- ?rest) / ?fact * ?pow"

	  by (metis mult_nonneg_nonneg divide_nonneg_pos less_imp_le)

	thus ?thesis unfolding cos_eq by auto

      qed

      also have "\<dots> \<le> Ifloat (ub_sin_cos_aux prec n 1 1 (x * x))"

	unfolding morph_to_if_power[symmetric] using cos_aux by auto

      finally have "cos (Ifloat x) \<le> Ifloat (ub_sin_cos_aux prec n 1 1 (x * x))" .

    } note ub = this and lb

  } note ub = this(1) and lb = this(2)

  have "cos (Ifloat x) \<le> Ifloat (ub_sin_cos_aux prec (get_odd n) 1 1 (x * x))" using ub[OF odd_pos[OF get_odd] get_odd] .

  moreover have "Ifloat (lb_sin_cos_aux prec (get_even n) 1 1 (x * x)) \<le> cos (Ifloat x)" 

  proof (cases "0 < get_even n")

    case True show ?thesis using lb[OF True get_even] .

  next

    case False

    hence "get_even n = 0" by auto

    have "- (pi / 2) \<le> Ifloat x" by (rule order_trans[OF _ `0 < Ifloat x`[THEN less_imp_le]], auto)

    with `Ifloat x \<le> pi / 2`

    show ?thesis unfolding `get_even n = 0` lb_sin_cos_aux.simps Ifloat_minus Ifloat_0 using cos_ge_zero by auto

  qed

  ultimately show ?thesis by auto

next

  case True

  show ?thesis

  proof (cases "n = 0")

    case True 

    thus ?thesis unfolding `n = 0` get_even_def get_odd_def using `Ifloat 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 `Ifloat 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

lemma sin_aux: assumes "0 \<le> Ifloat x"

  shows "Ifloat (x * lb_sin_cos_aux prec n 2 1 (x * x)) \<le> (\<Sum> i=0..<n. -1^i * (1/real (fact (2 * i + 1))) * (Ifloat x)^(2 * i + 1))" (is "?lb")

  and "(\<Sum> i=0..<n. -1^i * (1/real (fact (2 * i + 1))) * (Ifloat x)^(2 * i + 1)) \<le> Ifloat (x * ub_sin_cos_aux prec n 2 1 (x * x))" (is "?ub")

proof -

  have "0 \<le> Ifloat (x * x)" unfolding Ifloat_mult by auto

  let "?f n" = "fact (2 * n + 1)"

  { fix n 

    have F: "\<And>m. ((\<lambda>i. i + 2) o^ n) m = m + 2 * n" by (induct n arbitrary: m, auto)

    have "?f (Suc n) = ?f n * ((\<lambda>i. i + 2) o^ n) 2 * (((\<lambda>i. i + 2) o^ n) 2 + 1)"

      unfolding F by auto } note f_eq = this

  from horner_bounds[where lb="lb_sin_cos_aux prec" and ub="ub_sin_cos_aux prec" and j'=0,

    OF `0 \<le> Ifloat (x * x)` f_eq lb_sin_cos_aux.simps ub_sin_cos_aux.simps]

  show "?lb" and "?ub" using `0 \<le> Ifloat x` unfolding Ifloat_mult

    unfolding power_add power_one_right real_mult_assoc[symmetric] setsum_left_distrib[symmetric]

    unfolding real_mult_commute

    by (auto intro!: mult_left_mono simp add: power_mult power2_eq_square[of "Ifloat x"])

qed

lemma sin_boundaries: assumes "0 \<le> Ifloat x" and "Ifloat x \<le> pi / 2"

  shows "sin (Ifloat x) \<in> {Ifloat (x * lb_sin_cos_aux prec (get_even n) 2 1 (x * x)) .. Ifloat (x * ub_sin_cos_aux prec (get_odd n) 2 1 (x * x))}"

proof (cases "Ifloat x = 0")

  case False hence "Ifloat x \<noteq> 0" by auto

  hence "0 < x" and "0 < Ifloat x" using `0 \<le> Ifloat x` unfolding less_float_def by auto

  have "0 < x * x" using `0 < x` unfolding less_float_def Ifloat_mult Ifloat_0

    using mult_pos_pos[where a="Ifloat x" and b="Ifloat x"] by auto

  { fix x n have "(\<Sum> j = 0 ..< n. -1 ^ (((2 * j + 1) - Suc 0) div 2) / (real (fact (2 * j + 1))) * x ^(2 * j + 1))

    = (\<Sum> i = 0 ..< 2 * n. (if even(i) then 0 else (-1 ^ ((i - Suc 0) div 2))/(real (fact i))) * x ^ i)" (is "?SUM = _")

    proof -

      have pow: "!!i. x ^ (2 * i + 1) = x * x ^ (2 * i)" by auto

      have "?SUM = (\<Sum> j = 0 ..< n. 0) + ?SUM" by auto

      also have "\<dots> = (\<Sum> i = 0 ..< 2 * n. if even i then 0 else -1 ^ ((i - Suc 0) div 2) / (real (fact i)) * x ^ i)"

	unfolding sum_split_even_odd ..

      also have "\<dots> = (\<Sum> i = 0 ..< 2 * n. (if even i then 0 else -1 ^ ((i - Suc 0) div 2) / (real (fact i))) * x ^ i)"

	by (rule setsum_cong2) auto

      finally show ?thesis by assumption

    qed } note setsum_morph = this

  { fix n :: nat assume "0 < n"

    hence "0 < 2 * n + 1" by auto

    obtain t where "0 < t" and "t < Ifloat x" and

      sin_eq: "sin (Ifloat x) = (\<Sum> i = 0 ..< 2 * n + 1. (if even(i) then 0 else (-1 ^ ((i - Suc 0) div 2))/(real (fact i))) * (Ifloat x) ^ i) 

      + (sin (t + 1/2 * real (2 * n + 1) * pi) / real (fact (2*n + 1))) * (Ifloat x)^(2*n + 1)" 

      (is "_ = ?SUM + ?rest / ?fact * ?pow")

      using Maclaurin_sin_expansion3[OF `0 < 2 * n + 1` `0 < Ifloat x`] by auto

    have "?rest = cos t * -1^n" unfolding sin_add cos_add real_of_nat_add left_distrib right_distrib by auto

    moreover

    have "t \<le> pi / 2" using `t < Ifloat x` and `Ifloat x \<le> pi / 2` by auto

    hence "0 \<le> cos t" using `0 < t` and cos_ge_zero by auto

    ultimately have even: "even n \<Longrightarrow> 0 \<le> ?rest" and odd: "odd n \<Longrightarrow> 0 \<le> - ?rest " by auto

    have "0 < ?fact" by (rule real_of_nat_fact_gt_zero)

    have "0 < ?pow" using `0 < Ifloat x` by (rule zero_less_power)

    {

      assume "even n"

      have "Ifloat (x * lb_sin_cos_aux prec n 2 1 (x * x)) \<le> 

            (\<Sum> i = 0 ..< 2 * n. (if even(i) then 0 else (-1 ^ ((i - Suc 0) div 2))/(real (fact i))) * (Ifloat x) ^ i)"

	using sin_aux[OF `0 \<le> Ifloat x`] unfolding setsum_morph[symmetric] by auto

      also have "\<dots> \<le> ?SUM" by auto

      also have "\<dots> \<le> sin (Ifloat x)"

      proof -

	from even[OF `even n`] `0 < ?fact` `0 < ?pow`

	have "0 \<le> (?rest / ?fact) * ?pow" by (metis mult_nonneg_nonneg divide_nonneg_pos less_imp_le)

	thus ?thesis unfolding sin_eq by auto

      qed

      finally have "Ifloat (x * lb_sin_cos_aux prec n 2 1 (x * x)) \<le> sin (Ifloat x)" .

    } note lb = this

    {

      assume "odd n"

      have "sin (Ifloat x) \<le> ?SUM"

      proof -

	from `0 < ?fact` and `0 < ?pow` and odd[OF `odd n`]

	have "0 \<le> (- ?rest) / ?fact * ?pow"

	  by (metis mult_nonneg_nonneg divide_nonneg_pos less_imp_le)

	thus ?thesis unfolding sin_eq by auto

      qed

      also have "\<dots> \<le> (\<Sum> i = 0 ..< 2 * n. (if even(i) then 0 else (-1 ^ ((i - Suc 0) div 2))/(real (fact i))) * (Ifloat x) ^ i)"

	 by auto

      also have "\<dots> \<le> Ifloat (x * ub_sin_cos_aux prec n 2 1 (x * x))" 

	using sin_aux[OF `0 \<le> Ifloat x`] unfolding setsum_morph[symmetric] by auto

      finally have "sin (Ifloat x) \<le> Ifloat (x * ub_sin_cos_aux prec n 2 1 (x * x))" .

    } note ub = this and lb

  } note ub = this(1) and lb = this(2)

  have "sin (Ifloat x) \<le> Ifloat (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 "Ifloat (x * lb_sin_cos_aux prec (get_even n) 2 1 (x * x)) \<le> sin (Ifloat x)" 

  proof (cases "0 < get_even n")

    case True show ?thesis using lb[OF True get_even] .

  next