(* Author:     Johannes Hoelzl, TU Muenchen

   Coercions removed by Dmitriy Traytel *)

header {* Prove Real Valued Inequalities by Computation *}

theory Approximation

imports

  Complex_Main

  "~~/src/HOL/Library/Float"

  "~~/src/HOL/Library/Reflection"

  "~~/src/HOL/Decision_Procs/Dense_Linear_Order"

  "~~/src/HOL/Library/Efficient_Nat"

begin

section "Horner Scheme"

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

primrec 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 / 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 diff_minus 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 ^^ n) s) (f n)"

  shows "horner F G n ((F ^^ j') s) (f j') x = (\<Sum> j = 0..< n. -1 ^ j * (1 / (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 / (f (j' + j))"] by auto

qed auto

lemma horner_bounds':

  fixes lb :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float" and ub :: "nat \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> float \<Rightarrow> float"

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

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

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

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

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

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

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

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

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

  proof (rule add_mono)

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

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

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

          - (x * horner F G n (F ((F ^^ j') s)) (G ((F ^^ j') s) (f j')) x)"

      unfolding real_of_float_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 real_of_float_sub diff_minus

  proof (rule add_mono)

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

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

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

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

      unfolding real_of_float_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> real x" and f_Suc: "\<And>n. f (Suc n) = G ((F ^^ n) s) (f n)"

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

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

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

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

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

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

proof -

  have "?lb  \<and> ?ub"

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

    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 "real x \<le> 0" and f_Suc: "\<And>n. f (Suc n) = G ((F ^^ n) s) (f n)"

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

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

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

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

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

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

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: "(-x) = -1 * real x" by auto (* coercion "inside" is necessary *)

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

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

  proof (rule setsum_cong, simp)

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

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

      unfolding move_minus power_mult_distrib mult_assoc[symmetric]

      unfolding mult_commute unfolding mult_assoc[of "-1 ^ j", symmetric] power_mult_distrib[symmetric]

      by auto

  qed

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

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

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

    OF this f_Suc lb_0 refl ub_0 refl]

  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: fixes x :: real

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

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

proof (cases "even n")

  case True

  show ?thesis

  proof (cases "0 < l")

    case True hence "odd n \<or> 0 < l" and "0 \<le> real 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 "real l ^ n \<le> x ^ n" and "x ^ n \<le> real u ^ n " using `0 \<le> real l` and assms unfolding atLeastAtMost_iff using power_mono[of l x] power_mono[of x 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> - real u" and "- real u \<le> - x" and "0 \<le> - x" and "-x \<le> - real l" using assms unfolding less_float_def by auto

      hence "real u ^ n \<le> x ^ n" and "x ^ n \<le> real l ^ n" using power_mono[of  "-x" "-real l" n] power_mono[of "-real 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> real (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>real (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 "real l ^ n \<le> x ^ n" and "x ^ n \<le> real 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::real) l u. (l1, u1) = float_power_bnds n l u \<and> x \<in> {l .. u} \<longrightarrow> l1 \<le> x ^ n \<and> x ^ n \<le> u1"

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

function ub_sqrt lb_sqrt :: "nat \<Rightarrow> float \<Rightarrow> float" where

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

              else if x < 0 then - lb_sqrt prec (- x)

                            else 0)" |

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

              else if x < 0 then - ub_sqrt prec (- x)

                            else 0)"

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 lb_sqrt.simps[simp del]

declare ub_sqrt.simps[simp del]

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 < 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 float_pos_m_pos[unfolded less_float_def] assms by auto

    hence "0 < sqrt m" by auto

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

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

      unfolding pow2_add pow2_int Float real_of_float_simp by auto

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

    proof (rule mult_strict_right_mono, auto)

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

        unfolding real_of_int_less_iff[of m, symmetric] by auto

    qed

    finally have "sqrt 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 real_of_float_simp by auto

  qed

next

  case (Suc n)

  let ?b = "sqrt_iteration prec n x"

  have "0 < sqrt x" using `0 < real x` by auto

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

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

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

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

  finally show ?case unfolding sqrt_iteration.simps Let_def real_of_float_mult real_of_float_add right_distrib .

qed

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

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

  finally show ?thesis .

qed

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

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

proof (cases "0 < x")

  case True hence "0 < real x" and "0 \<le> x" using `0 \<le> real 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> real (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.simps using True by auto

next

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

  thus ?thesis unfolding lb_sqrt.simps less_float_def by auto

qed

lemma bnds_sqrt':

  shows "sqrt x \<in> {(lb_sqrt prec x) .. (ub_sqrt prec x) }"

proof -

  { fix x :: float assume "0 < x"

    hence "0 < real x" and "0 \<le> real x" unfolding less_float_def 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

    have "(float_divl prec x (sqrt_iteration prec prec x)) \<le>

          x / (sqrt_iteration prec prec x)" by (rule float_divl)

    also have "\<dots> < x / sqrt x"

      by (rule divide_strict_left_mono[OF sqrt_ub `0 < real x`

               mult_pos_pos[OF order_less_trans[OF sqrt_gt0 sqrt_ub] sqrt_gt0]])

    also have "\<dots> = sqrt x"

      unfolding inverse_eq_iff_eq[of _ "sqrt x", symmetric]

                sqrt_divide_self_eq[OF `0 \<le> real x`, symmetric] by auto

    finally have "lb_sqrt prec x \<le> 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" unfolding less_float_def 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 \<le> ub_sqrt prec x"

      unfolding ub_sqrt.simps if_P[OF `0 < x`] by auto }

  note ub = this

  show ?thesis

  proof (cases "0 < x")

    case True with lb ub show ?thesis by auto

  next case False show ?thesis

  proof (cases "real x = 0")

    case True thus ?thesis

      by (auto simp add: less_float_def lb_sqrt.simps ub_sqrt.simps)

  next

    case False with `\<not> 0 < x` have "x < 0" and "0 < -x"

      by (auto simp add: less_float_def)

    with `\<not> 0 < x`

    show ?thesis using lb[OF `0 < -x`] ub[OF `0 < -x`]

      by (auto simp add: real_sqrt_minus lb_sqrt.simps ub_sqrt.simps)

  qed qed

qed

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

proof ((rule allI) +, rule impI, erule conjE, rule conjI)

  fix x :: real fix lx ux

  assume "(l, u) = (lb_sqrt prec lx, ub_sqrt prec ux)"

    and x: "x \<in> {lx .. ux}"

  hence l: "l = lb_sqrt prec lx " and u: "u = ub_sqrt prec ux" by auto

  have "sqrt lx \<le> sqrt x" using x by auto

  from order_trans[OF _ this]

  show "l \<le> sqrt x" unfolding l using bnds_sqrt'[of lx prec] by auto

  have "sqrt x \<le> sqrt ux" using x by auto

  from order_trans[OF this]

  show "sqrt x \<le> u" unfolding u using bnds_sqrt'[of ux prec] by auto

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 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 k) - x * (ub_arctan_horner prec n (k + 2) x)"

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

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

proof -

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

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

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

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

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

  proof (cases "real x = 0")

    case False

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

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

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

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

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

  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> real (x*x)` F lb_arctan_horner.simps ub_arctan_horner.simps]

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

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

      unfolding real_of_float_mult power_add power_one_right mult_assoc[symmetric] setsum_left_distrib[symmetric]

      unfolding mult_commute[where 'a=real] mult_commute[of _ "2::nat"] power_mult power2_eq_square[of "real x"]

      by (auto intro!: mult_left_mono)

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

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

  moreover

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

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

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

      unfolding real_of_float_mult power_add power_one_right mult_assoc[symmetric] setsum_left_distrib[symmetric]

      unfolding mult_commute[where 'a=real] mult_commute[of _ "2::nat"] power_mult power2_eq_square[of "real x"]

      by (auto intro!: mult_left_mono)

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

  ultimately show ?thesis by auto

qed

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

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

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_Suc by auto

  have "arctan x \<le> 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> real x` `real x \<le> 1` `even n'`] by auto

  moreover

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

    unfolding get_even_def if_P[OF True] using arctan_0_1_bounds'[OF `0 \<le> real x` `real 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_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 x \<le> 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> real x` `real x \<le> 1` `even n'`] by auto

  moreover

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

    unfolding get_even_def if_not_P[OF False] unfolding `n = Suc n'` using arctan_0_1_bounds'[OF `0 \<le> real x` `real 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> {(lb_pi n) .. (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> real ?k" by (rule order_trans[OF _ rapprox_rat], auto simp add: `0 \<le> k`)

    have "real ?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 / k \<le> ?k" using rapprox_rat[where x=1 and y=k] by auto

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

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

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

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

  } 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 / k \<le> 1" using `1 < k` by auto

    have "\<And>n. 0 \<le> real ?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. real ?k \<le> 1" using lapprox_rat by (rule order_trans, auto simp add: `1 / k \<le> 1`)

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

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

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

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

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

  } note lb_arctan = this

  have "pi \<le> ub_pi n"

    unfolding ub_pi_def machin_pi Let_def real_of_float_mult real_of_float_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 "lb_pi n \<le> pi"

    unfolding lb_pi_def machin_pi Let_def real_of_float_mult real_of_float_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 + ub_sqrt prec (1 + x * x)))

                          else (let inv = float_divr prec 1 x

                                in if inv > 1 then 0

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

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

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

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

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

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

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

                                           else Float 1 1 * ub_horner y

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

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> real x"

  shows "lb_arctan prec x \<le> arctan x"

proof -

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

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

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

  show ?thesis

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

    case True hence "real 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> real x` `real x \<le> 1`] by auto

  next

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

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

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

    let ?DIV = "float_divl prec x ?fR"

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

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

    have "sqrt (1 + x * x) \<le> ub_sqrt prec (1 + x * x)"

      using bnds_sqrt'[of "1 + x * x"] by auto

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

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

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

    proof -

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

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

      finally show ?thesis .

    qed

    show ?thesis

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

      case True

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

      also have "\<dots> \<le> (ub_sqrt prec (1 + x * x))"

        using bnds_sqrt'[of "1 + x * x"] by auto

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

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

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

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

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

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

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

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

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

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

    next

      case False

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

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

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

      have "0 \<le> arctan x" using arctan_monotone'[OF `0 \<le> 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 `\<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 x` by auto

      next

        case False

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

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

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

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

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

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

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

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

        moreover

        have "lb_pi prec * Float 1 -1 \<le> pi / 2" unfolding real_of_float_mult Float_num times_divide_eq_right mult_1_left 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> real x"

  shows "arctan x \<le> ub_arctan prec x"

proof -

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

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

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

  show ?thesis

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

    case True hence "real 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> real x` `real x \<le> 1`] by auto

  next

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

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

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

    let ?DIV = "float_divr prec x ?fR"

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

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

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

    have "lb_sqrt prec (1 + x * x) \<le> sqrt (1 + x * x)"

      using bnds_sqrt'[of "1 + x * x"] by auto

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

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

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

    proof -

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

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

      also have "\<dots> \<le> ?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> ub_pi prec * Float 1 -1" unfolding real_of_float_mult 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 `\<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 "real ?DIV \<le> 1" unfolding less_float_def by auto

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

        hence "0 \<le> 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 "\<dots> \<le> 2 * arctan (?DIV)"

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

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

          using arctan_0_1_bounds[OF `0 \<le> real ?DIV` `real ?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 < real x" unfolding le_float_def Float_num by auto

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

      hence "0 < real x" by auto

      hence "0 < x" unfolding less_float_def by auto

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

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

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

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

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

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

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

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

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

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

      moreover

      have "pi / 2 \<le> ub_pi prec * Float 1 -1" unfolding real_of_float_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 x \<in> {(lb_arctan prec x) .. (ub_arctan prec x)}"

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

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

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

next

  let ?mx = "-x"

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

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

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

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

qed

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

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

  fix x :: real fix lx ux

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

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

  { from arctan_boundaries[of lx prec, unfolded l]

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

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

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

  } moreover

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

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

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

  } ultimately show "l \<le> arctan x \<and> arctan x \<le> 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 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 k) - x * (ub_sin_cos_aux prec n (i + 2) (k * i * (i + 1)) x)"

lemma cos_aux:

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

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

proof -

  have "0 \<le> real (x * x)" unfolding real_of_float_mult by auto

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

  { fix n

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

    have "?f (Suc n) = ?f n * ((\<lambda>i. i + 2) ^^ n) 1 * (((\<lambda>i. i + 2) ^^ 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> 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"])

qed

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

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

proof (cases "real x = 0")

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

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

  have "0 < x * x" using `0 < x` unfolding less_float_def real_of_float_mult real_of_float_0

    using mult_pos_pos[where a="real x" and b="real 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 < real x" and

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

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

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

      using Maclaurin_cos_expansion2[OF `0 < real x` `0 < 2 * n`]

      unfolding cos_coeff_def by auto

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

    also have "\<dots> = cos (t + 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 < real x` and `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 < real x` by auto

    {

      assume "even n"

      have "(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 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 "(lb_sin_cos_aux prec n 1 1 (x * x)) \<le> cos x" .

    } note lb = this

    {

      assume "odd n"

      have "cos 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> (ub_sin_cos_aux prec n 1 1 (x * x))"

        unfolding morph_to_if_power[symmetric] using cos_aux by auto

      finally have "cos x \<le> (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 x \<le> (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)) \<le> cos 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> x" by (rule order_trans[OF _ `0 < real x`[THEN less_imp_le]], auto)

    with `x \<le> pi / 2`

    show ?thesis unfolding `get_even n = 0` lb_sin_cos_aux.simps real_of_float_minus real_of_float_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 `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)

  qed

qed

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

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

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

proof -

  have "0 \<le> real (x * x)" unfolding real_of_float_mult by auto

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

  { fix n

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

    have "?f (Suc n) = ?f n * ((\<lambda>i. i + 2) ^^ n) 2 * (((\<lambda>i. i + 2) ^^ 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> real (x * x)` f_eq lb_sin_cos_aux.simps ub_sin_cos_aux.simps]

  show "?lb" and "?ub" using `0 \<le> real x` unfolding real_of_float_mult

    unfolding power_add power_one_right mult_assoc[symmetric] setsum_left_distrib[symmetric]

    unfolding mult_commute[where 'a=real]

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

qed

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

  shows "sin x \<in> {(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 \<noteq> 0" by auto

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

  have "0 < x * x" using `0 < x` unfolding less_float_def real_of_float_mult real_of_float_0

    using mult_pos_pos[where a="real x" and b="real 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 < real x" and

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

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

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

      using Maclaurin_sin_expansion3[OF `0 < 2 * n + 1` `0 < real x`]

      unfolding sin_coeff_def 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 < real x` and `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 (simp del: fact_Suc)

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

    {

      assume "even n"

      have "(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))) * (real x) ^ i)"

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

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

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

    } note lb = this

    {

      assume "odd n"

      have "sin 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))) * (real x) ^ i)"

         by auto

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