Quellcode-Bibliothek

^© Kompilation durch diese Firma

[Weder Korrektheit noch Funktionsfähigkeit der Software werden zugesichert.]

Datei: Approximation_Bounds.thy Sprache: Isabelle

Original von: Isabelle^©

(*
  Author:     Johannes Hoelzl, TU Muenchen
  Coercions removed by Dmitriy Traytel

  This file contains only general material about computing lower/upper bounds
  on real functions. Approximation.thy contains the actual approximation algorithm
  and the approximation oracle. This is in order to make a clear separation between
  "morally immaculate" material about upper/lower bounds and the trusted oracle/reflection.
*)

theory Approximation_Bounds
imports
  Complex_Main
  "HOL-Library.Interval_Float"
  Dense_Linear_Order
begin

declare powr_neg_one [simp]
declare powr_neg_numeral [simp]

context includes interval.lifting begin

section "Horner Scheme"

subsection \<open>Define auxiliary helper \<open>horner\<close> function\<close>

primrec horner :: "(nat \ nat) \ (nat \ nat \ nat) \ nat \ nat \ nat \ real \ 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 \ real"
  shows "a 0 - x * (\ i=0.. i=0..
proof -
  have shift_pow: "\i. - (x * ((-1)^i * a (Suc i) * x ^ i)) = (-1)^(Suc i) * a (Suc i) * x ^ (Suc i)"
    by auto
  show ?thesis
    unfolding sum_distrib_left shift_pow uminus_add_conv_diff [symmetric] sum_negf[symmetric]
    sum.atLeast_Suc_lessThan[OF zero_less_Suc]
    sum.reindex[OF inj_Suc, unfolded comp_def, symmetric, of "\ n. (-1)^n *a n * x^n"] by auto
qed

lemma horner_schema:
  fixes f :: "nat \ nat" and G :: "nat \ nat \ nat" and F :: "nat \ nat"
  assumes f_Suc: "\n. f (Suc n) = G ((F ^^ n) s) (f n)"
  shows "horner F G n ((F ^^ j') s) (f j') x = (\ j = 0..< n. (- 1) ^ j * (1 / (f (j' + j))) * x ^ j)"
proof (induct n arbitrary: j')
  case 0
  then show ?case by auto
next
  case (Suc n)
  show ?case unfolding horner.simps Suc[where j'="Suc j'", unfolded funpow.simps comp_def f_Suc]
    using horner_schema'[of "\ j. 1 / (f (j' + j))"] by auto
qed

lemma horner_bounds':
  fixes lb :: "nat \ nat \ nat \ float \ float" and ub :: "nat \ nat \ nat \ float \ float"
  assumes "0 \ real_of_float x" and f_Suc: "\n. f (Suc n) = G ((F ^^ n) s) (f n)"
    and lb_0: "\ i k x. lb 0 i k x = 0"
    and lb_Suc: "\ n i k x. lb (Suc n) i k x = float_plus_down prec
        (lapprox_rat prec 1 k)
        (- float_round_up prec (x * (ub n (F i) (G i k) x)))"
    and ub_0: "\ i k x. ub 0 i k x = 0"
    and ub_Suc: "\ n i k x. ub (Suc n) i k x = float_plus_up prec
        (rapprox_rat prec 1 k)
        (- float_round_down prec (x * (lb n (F i) (G i k) x)))"
  shows "(lb n ((F ^^ j') s) (f j') x) \ horner F G n ((F ^^ j') s) (f j') x \
         horner F G n ((F ^^ j') s) (f j') x \<le> (ub n ((F ^^ j') s) (f j') x)"
  (is "?lb n j' \ ?horner n j' \ ?horner n j' \ ?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)
  thus ?case using lapprox_rat[of prec 1 "f j'"] using rapprox_rat[of 1 "f j'" prec]
    Suc[where j'="Suc j'"] \0 \ real_of_float x\
    by (auto intro!: add_mono mult_left_mono float_round_down_le float_round_up_le
      order_trans[OF add_mono[OF _ float_plus_down_le]]
      order_trans[OF _ add_mono[OF _ float_plus_up_le]]
      simp add: lb_Suc ub_Suc field_simps f_Suc)
qed

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

text \<open>

Here \<^term_type>\<open>f :: nat \<Rightarrow> nat\<close> is the sequence defining the Taylor series, the coefficients are
all alternating and reciprocs. We use \<^term>\<open>G\<close> and \<^term>\<open>F\<close> to describe the computation of \<^term>\<open>f\<close>.

\<close>

lemma horner_bounds:
  fixes F :: "nat \ nat" and G :: "nat \ nat \ nat"
  assumes "0 \ real_of_float x" and f_Suc: "\n. f (Suc n) = G ((F ^^ n) s) (f n)"
    and lb_0: "\ i k x. lb 0 i k x = 0"
    and lb_Suc: "\ n i k x. lb (Suc n) i k x = float_plus_down prec
        (lapprox_rat prec 1 k)
        (- float_round_up prec (x * (ub n (F i) (G i k) x)))"
    and ub_0: "\ i k x. ub 0 i k x = 0"
    and ub_Suc: "\ n i k x. ub (Suc n) i k x = float_plus_up prec
        (rapprox_rat prec 1 k)
        (- float_round_down prec (x * (lb n (F i) (G i k) x)))"
  shows "(lb n ((F ^^ j') s) (f j') x) \ (\j=0..
      (is "?lb")
    and "(\j=0.. (ub n ((F ^^ j') s) (f j') x)"
      (is "?ub")
proof -
  have "?lb \ ?ub"
    using horner_bounds'[where lb=lb, OF \0 \ real_of_float x\ f_Suc lb_0 lb_Suc ub_0 ub_Suc]
    unfolding horner_schema[where f=f, OF f_Suc] by simp
  thus "?lb" and "?ub" by auto
qed

lemma horner_bounds_nonpos:
  fixes F :: "nat \ nat" and G :: "nat \ nat \ nat"
  assumes "real_of_float x \ 0" and f_Suc: "\n. f (Suc n) = G ((F ^^ n) s) (f n)"
    and lb_0: "\ i k x. lb 0 i k x = 0"
    and lb_Suc: "\ n i k x. lb (Suc n) i k x = float_plus_down prec
        (lapprox_rat prec 1 k)
        (float_round_down prec (x * (ub n (F i) (G i k) x)))"
    and ub_0: "\ i k x. ub 0 i k x = 0"
    and ub_Suc: "\ n i k x. ub (Suc n) i k x = float_plus_up prec
        (rapprox_rat prec 1 k)
        (float_round_up prec (x * (lb n (F i) (G i k) x)))"
  shows "(lb n ((F ^^ j') s) (f j') x) \ (\j=0..
    and "(\j=0.. (ub n ((F ^^ j') s) (f j') x)" (is "?ub")
proof -
  have diff_mult_minus: "x - y * z = x + - y * z" for x y z :: float by simp
  have sum_eq: "(\j=0..
    (\<Sum>j = 0..<n. (- 1) ^ j * (1 / (f (j' + j))) * real_of_float (- x) ^ j)"
    by (auto simp add: field_simps power_mult_distrib[symmetric])
  have "0 \ real_of_float (-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="\ n i k x. lb n i k (-x)" and ub="\ n i k x. ub n i k (-x)",
    unfolded lb_Suc ub_Suc diff_mult_minus,
    OF this f_Suc lb_0 _ ub_0 _]
  show "?lb" and "?ub" unfolding minus_minus sum_eq
    by (auto simp: minus_float_round_up_eq minus_float_round_down_eq)
qed

subsection \<open>Selectors for next even or odd number\<close>

text \<open>
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>\<open>get_odd\<close> and \<^term>\<open>get_even\<close>.
\<close>

definition get_odd :: "nat \ nat" where
  "get_odd n = (if odd n then n else (Suc n))"

definition get_even :: "nat \ nat" where
  "get_even n = (if even n then n else (Suc n))"

lemma get_odd[simp]: "odd (get_odd n)"
  unfolding get_odd_def by (cases "odd n") auto

lemma get_even[simp]: "even (get_even n)"
  unfolding get_even_def by (cases "even n") auto

lemma get_odd_ex: "\ k. Suc k = get_odd n \ odd (Suc k)"
  by (auto simp: get_odd_def odd_pos intro!: exI[of _ "n - 1"])

lemma get_even_double: "\i. get_even n = 2 * i"
  using get_even by (blast elim: evenE)

lemma get_odd_double: "\i. get_odd n = 2 * i + 1"
  using get_odd by (blast elim: oddE)

section "Power function"

definition float_power_bnds :: "nat \ nat \ float \ float \ float * float" where
"float_power_bnds prec n l u =
  (if 0 < l then (power_down_fl prec l n, power_up_fl prec u n)
  else if odd n then
    (- power_up_fl prec \<bar>l\<bar> n,
      if u < 0 then - power_down_fl prec \<bar>u\<bar> n else power_up_fl prec u n)
  else if u < 0 then (power_down_fl prec \<bar>u\<bar> n, power_up_fl prec \<bar>l\<bar> n)
  else (0, power_up_fl prec (max \<bar>l\<bar> \<bar>u\<bar>) n))"

lemma le_minus_power_downI: "0 \ x \ x ^ n \ - a \ a \ - power_down prec x n"
  by (subst le_minus_iff) (auto intro: power_down_le power_mono_odd)

lemma float_power_bnds:
  "(l1, u1) = float_power_bnds prec n l u \ x \ {l .. u} \ (x::real) ^ n \ {l1..u1}"
  by (auto
    simp: float_power_bnds_def max_def real_power_up_fl real_power_down_fl minus_le_iff
    split: if_split_asm
    intro!: power_up_le power_down_le le_minus_power_downI
    intro: power_mono_odd power_mono power_mono_even zero_le_even_power)

lemma bnds_power:
  "\(x::real) l u. (l1, u1) = float_power_bnds prec n l u \ x \ {l .. u} \
    l1 \<le> x ^ n \<and> x ^ n \<le> u1"
  using float_power_bnds by auto

lift_definition power_float_interval :: "nat \ nat \ float interval \ float interval"
  is "\p n (l, u). float_power_bnds p n l u"
  using float_power_bnds
  by (auto simp: bnds_power dest!: float_power_bnds[OF sym])

lemma lower_power_float_interval:
  "lower (power_float_interval p n x) = fst (float_power_bnds p n (lower x) (upper x))"
  by transfer auto
lemma upper_power_float_interval:
  "upper (power_float_interval p n x) = snd (float_power_bnds p n (lower x) (upper x))"
  by transfer auto

lemma power_float_intervalI: "x \\<^sub>r X \ x ^ n \\<^sub>r power_float_interval p n X"
  using float_power_bnds[OF prod.collapse]
  by (auto simp: set_of_eq lower_power_float_interval upper_power_float_interval)

lemma power_float_interval_mono:
  "set_of (power_float_interval prec n A)
    \<subseteq> set_of (power_float_interval prec n B)"
  if "set_of A \ set_of B"
proof -
  define la where "la = real_of_float (lower A)"
  define ua where "ua = real_of_float (upper A)"
  define lb where "lb = real_of_float (lower B)"
  define ub where "ub = real_of_float (upper B)"
  have ineqs: "lb \ la" "la \ ua" "ua \ ub" "lb \ ub"
    using that lower_le_upper[of A] lower_le_upper[of B]
    by (auto simp: la_def ua_def lb_def ub_def set_of_eq)
  show ?thesis
    using ineqs
    by (simp add: set_of_subset_iff float_power_bnds_def max_def
        power_down_fl.rep_eq power_up_fl.rep_eq
        lower_power_float_interval upper_power_float_interval
        la_def[symmetric] ua_def[symmetric] lb_def[symmetric] ub_def[symmetric])
      (auto intro!: power_down_mono power_up_mono intro: order_trans[where y=0])
qed

section \<open>Approximation utility functions\<close>

definition bnds_mult :: "nat \ float \ float \ float \ float \ float \ float" where
  "bnds_mult prec a1 a2 b1 b2 =
      (float_plus_down prec (nprt a1 * pprt b2)
          (float_plus_down prec (nprt a2 * nprt b2)
            (float_plus_down prec (pprt a1 * pprt b1) (pprt a2 * nprt b1))),
        float_plus_up prec (pprt a2 * pprt b2)
            (float_plus_up prec (pprt a1 * nprt b2)
              (float_plus_up prec (nprt a2 * pprt b1) (nprt a1 * nprt b1))))"

lemma bnds_mult:
  fixes prec :: nat and a1 aa2 b1 b2 :: float
  assumes "(l, u) = bnds_mult prec a1 a2 b1 b2"
  assumes "a \ {real_of_float a1..real_of_float a2}"
  assumes "b \ {real_of_float b1..real_of_float b2}"
  shows   "a * b \ {real_of_float l..real_of_float u}"
proof -
  from assms have "real_of_float l \ a * b"
    by (intro order.trans[OF _ mult_ge_prts[of a1 a a2 b1 b b2]])
       (auto simp: bnds_mult_def intro!: float_plus_down_le)
  moreover from assms have "real_of_float u \ a * b"
    by (intro order.trans[OF mult_le_prts[of a1 a a2 b1 b b2]])
       (auto simp: bnds_mult_def intro!: float_plus_up_le)
  ultimately show ?thesis by simp
qed

lift_definition mult_float_interval::"nat \ float interval \ float interval \ float interval"
  is "\prec. \(a1, a2). \(b1, b2). bnds_mult prec a1 a2 b1 b2"
  by (auto dest!: bnds_mult[OF sym])

lemma lower_mult_float_interval:
  "lower (mult_float_interval p x y) = fst (bnds_mult p (lower x) (upper x) (lower y) (upper y))"
  by transfer auto
lemma upper_mult_float_interval:
  "upper (mult_float_interval p x y) = snd (bnds_mult p (lower x) (upper x) (lower y) (upper y))"
  by transfer auto

lemma mult_float_interval:
  "set_of (real_interval A) * set_of (real_interval B) \
    set_of (real_interval (mult_float_interval prec A B))"
proof -
  let ?bm = "bnds_mult prec (lower A) (upper A) (lower B) (upper B)"
  show ?thesis
    using bnds_mult[of "fst ?bm" "snd ?bm", simplified, OF refl]
    by (auto simp: set_of_eq set_times_def upper_mult_float_interval lower_mult_float_interval)
qed

lemma mult_float_intervalI:
  "x * y \\<^sub>r mult_float_interval prec A B"
  if "x \\<^sub>i real_interval A" "y \\<^sub>i real_interval B"
  using mult_float_interval[of A B] that
  by (auto simp: )

lemma mult_float_interval_mono:
  "set_of (mult_float_interval prec A B) \ set_of (mult_float_interval prec X Y)"
  if "set_of A \ set_of X" "set_of B \ set_of Y"
  using that
  apply transfer
  unfolding bnds_mult_def atLeastatMost_subset_iff float_plus_down.rep_eq float_plus_up.rep_eq
  by (auto simp: float_plus_down.rep_eq float_plus_up.rep_eq mult_float_mono1 mult_float_mono2)

definition map_bnds :: "(nat \ float \ float) \ (nat \ float \ float) \
                           nat \<Rightarrow> (float \<times> float) \<Rightarrow> (float \<times> float)" where
  "map_bnds lb ub prec = (\(l,u). (lb prec l, ub prec u))"

lemma map_bnds:
  assumes "(lf, uf) = map_bnds lb ub prec (l, u)"
  assumes "mono f"
  assumes "x \ {real_of_float l..real_of_float u}"
  assumes "real_of_float (lb prec l) \ f (real_of_float l)"
  assumes "real_of_float (ub prec u) \ f (real_of_float u)"
  shows   "f x \ {real_of_float lf..real_of_float uf}"
proof -
  from assms have "real_of_float lf = real_of_float (lb prec l)"
    by (simp add: map_bnds_def)
  also have "real_of_float (lb prec l) \ f (real_of_float l)" by fact
  also from assms have "\ \ f x"
    by (intro monoD[OF \<open>mono f\<close>]) auto
  finally have lf: "real_of_float lf \ f x" .

  from assms have "f x \ f (real_of_float u)"
    by (intro monoD[OF \<open>mono f\<close>]) auto
  also have "\ \ real_of_float (ub prec u)" by fact
  also from assms have "\ = real_of_float uf"
    by (simp add: map_bnds_def)
  finally have uf: "f x \ real_of_float uf" .

  from lf uf show ?thesis by simp
qed

section "Square root"

text \<open>
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.
\<close>

fun sqrt_iteration :: "nat \ nat \ float \ float" where
"sqrt_iteration prec 0 x = Float 1 ((bitlen \mantissa x\ + exponent x) div 2 + 1)" |
"sqrt_iteration prec (Suc m) x = (let y = sqrt_iteration prec m x
                                  in Float 1 (- 1) * float_plus_up prec y (float_divr prec x y))"

lemma compute_sqrt_iteration_base[code]:
  shows "sqrt_iteration prec n (Float m e) =
    (if n = 0 then Float 1 ((if m = 0 then 0 else bitlen \<bar>m\<bar> + e) div 2 + 1)
    else (let y = sqrt_iteration prec (n - 1) (Float m e) in
      Float 1 (- 1) * float_plus_up prec y (float_divr prec (Float m e) y)))"
  using bitlen_Float by (cases n) simp_all

function ub_sqrt lb_sqrt :: "nat \ float \ 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 (\ v. let (prec, x) = case_sum id id v in (if x < 0 then 1 else 0))", auto)

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)\<^sup>2 " by simp
  also have "\ = b\<^sup>2 - 2 * b * sqrt x + (sqrt x)\<^sup>2" by algebra
  also have "\ = b\<^sup>2 - 2 * b * sqrt x + x" using assms by simp
  finally have "0 < b\<^sup>2 - 2 * b * sqrt x + x" .
  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_of_float x"
  shows "sqrt x < sqrt_iteration prec n x"
proof (induct n)
  case 0
  show ?case
  proof (cases x)
    case (Float m e)
    hence "0 < m"
      using assms
      by (auto simp: algebra_split_simps)
    hence "0 < sqrt m" by auto

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

    have "x = (m / 2^nat (bitlen m)) * 2 powr (e + (nat (bitlen m)))"
      unfolding Float by (auto simp: powr_realpow[symmetric] field_simps powr_add)
    also have "\ < 1 * 2 powr (e + nat (bitlen m))"
    proof (rule mult_strict_right_mono, auto)
      show "m < 2^nat (bitlen m)"
        using bitlen_bounds[OF \<open>0 < m\<close>, THEN conjunct2]
        unfolding of_int_less_iff[of m, symmetric] by auto
    qed
    finally have "sqrt x < sqrt (2 powr (e + bitlen m))"
      unfolding int_nat_bl by auto
    also have "\ \ 2 powr ((e + bitlen m) div 2 + 1)"
    proof -
      let ?E = "e + bitlen m"
      have E_mod_pow: "2 powr (?E mod 2) < 4"
      proof (cases "?E mod 2 = 1")
        case True
        thus ?thesis by auto
      next
        case False
        have "0 \ ?E mod 2" by auto
        have "?E mod 2 < 2" by auto
        from this[THEN zless_imp_add1_zle]
        have "?E mod 2 \ 0" using False by auto
        from xt1(5)[OF \<open>0 \<le> ?E mod 2\<close> this]
        show ?thesis by auto
      qed
      hence "sqrt (2 powr (?E mod 2)) < sqrt (2 * 2)"
        by (intro real_sqrt_less_mono) auto
      hence E_mod_pow: "sqrt (2 powr (?E mod 2)) < 2" by auto

      have E_eq: "2 powr ?E = 2 powr (?E div 2 + ?E div 2 + ?E mod 2)"
        by auto
      have "sqrt (2 powr ?E) = sqrt (2 powr (?E div 2) * 2 powr (?E div 2) * 2 powr (?E mod 2))"
        unfolding E_eq unfolding powr_add[symmetric] by (metis of_int_add)
      also have "\ = 2 powr (?E div 2) * sqrt (2 powr (?E mod 2))"
        unfolding real_sqrt_mult[of _ "2 powr (?E mod 2)"] real_sqrt_abs2 by auto
      also have "\ < 2 powr (?E div 2) * 2 powr 1"
        by (rule mult_strict_left_mono) (auto intro: E_mod_pow)
      also have "\ = 2 powr (?E div 2 + 1)"
        unfolding add.commute[of _ 1] powr_add[symmetric] by simp
      finally show ?thesis by auto
    qed
    finally show ?thesis using \<open>0 < m\<close>
      unfolding Float
      by (subst compute_sqrt_iteration_base) (simp add: ac_simps)
  qed
next
  case (Suc n)
  let ?b = "sqrt_iteration prec n x"
  have "0 < sqrt x"
    using \<open>0 < real_of_float x\<close> by auto
  also have "\ < real_of_float ?b"
    using Suc .
  finally have "sqrt x < (?b + x / ?b)/2"
    using sqrt_ub_pos_pos_1[OF Suc _ \<open>0 < real_of_float x\<close>] by auto
  also have "\ \ (?b + (float_divr prec x ?b))/2"
    by (rule divide_right_mono, auto simp add: float_divr)
  also have "\ = (Float 1 (- 1)) * (?b + (float_divr prec x ?b))"
    by simp
  also have "\ \ (Float 1 (- 1)) * (float_plus_up prec ?b (float_divr prec x ?b))"
    by (auto simp add: algebra_simps float_plus_up_le)
  finally show ?case
    unfolding sqrt_iteration.simps Let_def distrib_left .
qed

lemma sqrt_iteration_lower_bound:
  assumes "0 < real_of_float x"
  shows "0 < real_of_float (sqrt_iteration prec n x)" (is "0 < ?sqrt")
proof -
  have "0 < sqrt x" using assms by auto
  also have "\ < ?sqrt" using sqrt_iteration_bound[OF assms] .
  finally show ?thesis .
qed

lemma lb_sqrt_lower_bound:
  assumes "0 \ real_of_float x"
  shows "0 \ real_of_float (lb_sqrt prec x)"
proof (cases "0 < x")
  case True
  hence "0 < real_of_float x" and "0 \ x"
    using \<open>0 \<le> real_of_float x\<close> by auto
  hence "0 < sqrt_iteration prec prec x"
    using sqrt_iteration_lower_bound by auto
  hence "0 \ real_of_float (float_divl prec x (sqrt_iteration prec prec x))"
    using float_divl_lower_bound[OF \<open>0 \<le> x\<close>] unfolding less_eq_float_def by auto
  thus ?thesis
    unfolding lb_sqrt.simps using True by auto
next
  case False
  with \<open>0 \<le> real_of_float x\<close> have "real_of_float x = 0" by auto
  thus ?thesis
    unfolding lb_sqrt.simps by auto
qed

lemma bnds_sqrt': "sqrt x \ {(lb_sqrt prec x) .. (ub_sqrt prec x)}"
proof -
  have lb: "lb_sqrt prec x \ sqrt x" if "0 < x" for x :: float
  proof -
    from that have "0 < real_of_float x" and "0 \ real_of_float x" by auto
    hence sqrt_gt0: "0 < sqrt x" by auto
    hence sqrt_ub: "sqrt x < sqrt_iteration prec prec x"
      using sqrt_iteration_bound by auto
    have "(float_divl prec x (sqrt_iteration prec prec x)) \
          x / (sqrt_iteration prec prec x)" by (rule float_divl)
    also have "\ < x / sqrt x"
      by (rule divide_strict_left_mono[OF sqrt_ub \<open>0 < real_of_float x\<close>
               mult_pos_pos[OF order_less_trans[OF sqrt_gt0 sqrt_ub] sqrt_gt0]])
    also have "\ = sqrt x"
      unfolding inverse_eq_iff_eq[of _ "sqrt x", symmetric]
                sqrt_divide_self_eq[OF \<open>0 \<le> real_of_float x\<close>, symmetric] by auto
    finally show ?thesis
      unfolding lb_sqrt.simps if_P[OF \<open>0 < x\<close>] by auto
  qed
  have ub: "sqrt x \ ub_sqrt prec x" if "0 < x" for x :: float
  proof -
    from that have "0 < real_of_float x" by auto
    hence "0 < sqrt x" by auto
    hence "sqrt x < sqrt_iteration prec prec x"
      using sqrt_iteration_bound by auto
    then show ?thesis
      unfolding ub_sqrt.simps if_P[OF \<open>0 < x\<close>] by auto
  qed
  show ?thesis
    using lb[of "-x"] ub[of "-x"] lb[of x] ub[of x]
    by (auto simp add: lb_sqrt.simps ub_sqrt.simps real_sqrt_minus)
qed

lemma bnds_sqrt: "\(x::real) lx ux.
  (l, u) = (lb_sqrt prec lx, ub_sqrt prec ux) \<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 \ {lx .. ux}"
  hence l: "l = lb_sqrt prec lx " and u: "u = ub_sqrt prec ux" by auto

  have "sqrt lx \ sqrt x" using x by auto
  from order_trans[OF _ this]
  show "l \ sqrt x" unfolding l using bnds_sqrt'[of lx prec] by auto

  have "sqrt x \ sqrt ux" using x by auto
  from order_trans[OF this]
  show "sqrt x \ u" unfolding u using bnds_sqrt'[of ux prec] by auto
qed

lift_definition sqrt_float_interval::"nat \ float interval \ float interval"
  is "\prec. \(lx, ux). (lb_sqrt prec lx, ub_sqrt prec ux)"
  using bnds_sqrt'
  by auto (meson order_trans real_sqrt_le_iff)

lemma lower_float_interval: "lower (sqrt_float_interval prec X) = lb_sqrt prec (lower X)"
  by transfer auto

lemma upper_float_interval: "upper (sqrt_float_interval prec X) = ub_sqrt prec (upper X)"
  by transfer auto

lemma sqrt_float_interval:
  "sqrt ` set_of (real_interval X) \ set_of (real_interval (sqrt_float_interval prec X))"
  using bnds_sqrt
  by (auto simp: set_of_eq lower_float_interval upper_float_interval)

lemma sqrt_float_intervalI: "sqrt x \\<^sub>r sqrt_float_interval p X" if "x \\<^sub>r X"
  using sqrt_float_interval[of X p] that
  by auto

section "Arcus tangens and \"

subsection "Compute arcus tangens series"

text \<open>
As first step we implement the computation of the arcus tangens series. This is only valid in the range
\<^term>\<open>{-1 :: real .. 1}\<close>. This is used to compute \<pi> and then the entire arcus tangens.
\<close>

fun ub_arctan_horner :: "nat \ nat \ nat \ float \ float"
and lb_arctan_horner :: "nat \ nat \ nat \ float \ float" where
  "ub_arctan_horner prec 0 k x = 0"
| "ub_arctan_horner prec (Suc n) k x = float_plus_up prec
      (rapprox_rat prec 1 k) (- float_round_down prec (x * (lb_arctan_horner prec n (k + 2) x)))"
| "lb_arctan_horner prec 0 k x = 0"
| "lb_arctan_horner prec (Suc n) k x = float_plus_down prec
      (lapprox_rat prec 1 k) (- float_round_up prec (x * (ub_arctan_horner prec n (k + 2) x)))"

lemma arctan_0_1_bounds':
  assumes "0 \ real_of_float y" "real_of_float y \ 1"
    and "even n"
  shows "arctan (sqrt y) \
      {(sqrt y * lb_arctan_horner prec n 1 y) .. (sqrt y * ub_arctan_horner prec (Suc n) 1 y)}"
proof -
  let ?c = "\i. (- 1) ^ i * (1 / (i * 2 + (1::nat)) * sqrt y ^ (i * 2 + 1))"
  let ?S = "\n. \ i=0..

  have "0 \ sqrt y" using assms by auto
  have "sqrt y \ 1" using assms by auto
  from \<open>even n\<close> obtain m where "2 * m = n" by (blast elim: evenE)

  have "arctan (sqrt y) \ { ?S n .. ?S (Suc n) }"
  proof (cases "sqrt y = 0")
    case True
    then show ?thesis by simp
  next
    case False
    hence "0 < sqrt y" using \<open>0 \<le> sqrt y\<close> by auto
    hence prem: "0 < 1 / (0 * 2 + (1::nat)) * sqrt y ^ (0 * 2 + 1)" by auto

    have "\ sqrt y \ \ 1" using \0 \ sqrt y\ \sqrt y \ 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 \<open>2 * m = n\<close>]
    show ?thesis unfolding arctan_series[OF \<open>\<bar> sqrt y \<bar> \<le> 1\<close>] Suc_eq_plus1 atLeast0LessThan .
  qed
  note arctan_bounds = this[unfolded atLeastAtMost_iff]

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

  note bounds = horner_bounds[where s=1 and f="\i. 2 * i + 1" and j'=0
    and lb="\n i k x. lb_arctan_horner prec n k x"
    and ub="\n i k x. ub_arctan_horner prec n k x",
    OF \<open>0 \<le> real_of_float y\<close> F lb_arctan_horner.simps ub_arctan_horner.simps]

  have "(sqrt y * lb_arctan_horner prec n 1 y) \ arctan (sqrt y)"
  proof -
    have "(sqrt y * lb_arctan_horner prec n 1 y) \ ?S n"
      using bounds(1) \<open>0 \<le> sqrt y\<close>
      apply (simp only: power_add power_one_right mult.assoc[symmetric] sum_distrib_right[symmetric])
      apply (simp only: mult.commute[where 'a=real] mult.commute[of _ "2::nat"] power_mult)
      apply (auto intro!: mult_left_mono)
      done
    also have "\ \ arctan (sqrt y)" using arctan_bounds ..
    finally show ?thesis .
  qed
  moreover
  have "arctan (sqrt y) \ (sqrt y * ub_arctan_horner prec (Suc n) 1 y)"
  proof -
    have "arctan (sqrt y) \ ?S (Suc n)" using arctan_bounds ..
    also have "\ \ (sqrt y * ub_arctan_horner prec (Suc n) 1 y)"
      using bounds(2)[of "Suc n"] \<open>0 \<le> sqrt y\<close>
      apply (simp only: power_add power_one_right mult.assoc[symmetric] sum_distrib_right[symmetric])
      apply (simp only: mult.commute[where 'a=real] mult.commute[of _ "2::nat"] power_mult)
      apply (auto intro!: mult_left_mono)
      done
    finally show ?thesis .
  qed
  ultimately show ?thesis by auto
qed

lemma arctan_0_1_bounds:
  assumes "0 \ real_of_float y" "real_of_float y \ 1"
  shows "arctan (sqrt y) \
    {(sqrt y * lb_arctan_horner prec (get_even n) 1 y) ..
      (sqrt y * ub_arctan_horner prec (get_odd n) 1 y)}"
  using
    arctan_0_1_bounds'[OF assms, of n prec]
    arctan_0_1_bounds'[OF assms, of "n + 1" prec]
    arctan_0_1_bounds'[OF assms, of "n - 1" prec]
  by (auto simp: get_even_def get_odd_def odd_pos
    simp del: ub_arctan_horner.simps lb_arctan_horner.simps)

lemma arctan_lower_bound:
  assumes "0 \ x"
  shows "x / (1 + x\<^sup>2) \ arctan x" (is "?l x \ _")
proof -
  have "?l x - arctan x \ ?l 0 - arctan 0"
    using assms
    by (intro DERIV_nonpos_imp_nonincreasing[where f="\x. ?l x - arctan x"])
      (auto intro!: derivative_eq_intros simp: add_nonneg_eq_0_iff field_simps)
  thus ?thesis by simp
qed

lemma arctan_divide_mono: "0 < x \ x \ y \ arctan y / y \ arctan x / x"
  by (rule DERIV_nonpos_imp_nonincreasing[where f="\x. arctan x / x"])
    (auto intro!: derivative_eq_intros divide_nonpos_nonneg
      simp: inverse_eq_divide arctan_lower_bound)

lemma arctan_mult_mono: "0 \ x \ x \ y \ x * arctan y \ y * arctan x"
  using arctan_divide_mono[of x y] by (cases "x = 0") (simp_all add: field_simps)

lemma arctan_mult_le:
  assumes "0 \ x" "x \ y" "y * z \ arctan y"
  shows "x * z \ arctan x"
proof (cases "x = 0")
  case True
  then show ?thesis by simp
next
  case False
  with assms have "z \ arctan y / y" by (simp add: field_simps)
  also have "\ \ arctan x / x" using assms \x \ 0\ by (auto intro!: arctan_divide_mono)
  finally show ?thesis using assms \<open>x \<noteq> 0\<close> by (simp add: field_simps)
qed

lemma arctan_le_mult:
  assumes "0 < x" "x \ y" "arctan x \ x * z"
  shows "arctan y \ y * z"
proof -
  from assms have "arctan y / y \ arctan x / x" by (auto intro!: arctan_divide_mono)
  also have "\ \ z" using assms by (auto simp: field_simps)
  finally show ?thesis using assms by (simp add: field_simps)
qed

lemma arctan_0_1_bounds_le:
  assumes "0 \ x" "x \ 1" "0 < real_of_float xl" "real_of_float xl \ x * x" "x * x \ real_of_float xu" "real_of_float xu \ 1"
  shows "arctan x \
      {x * lb_arctan_horner p1 (get_even n) 1 xu .. x * ub_arctan_horner p2 (get_odd n) 1 xl}"
proof -
  from assms have "real_of_float xl \ 1" "sqrt (real_of_float xl) \ x" "x \ sqrt (real_of_float xu)" "0 \ real_of_float xu"
    "0 \ real_of_float xl" "0 < sqrt (real_of_float xl)"
    by (auto intro!: real_le_rsqrt real_le_lsqrt simp: power2_eq_square)
  from arctan_0_1_bounds[OF \<open>0 \<le> real_of_float xu\<close>  \<open>real_of_float xu \<le> 1\<close>]
  have "sqrt (real_of_float xu) * real_of_float (lb_arctan_horner p1 (get_even n) 1 xu) \ arctan (sqrt (real_of_float xu))"
    by simp
  from arctan_mult_le[OF \<open>0 \<le> x\<close> \<open>x \<le> sqrt _\<close>  this]
  have "x * real_of_float (lb_arctan_horner p1 (get_even n) 1 xu) \ arctan x" .
  moreover
  from arctan_0_1_bounds[OF \<open>0 \<le> real_of_float xl\<close>  \<open>real_of_float xl \<le> 1\<close>]
  have "arctan (sqrt (real_of_float xl)) \ sqrt (real_of_float xl) * real_of_float (ub_arctan_horner p2 (get_odd n) 1 xl)"
    by simp
  from arctan_le_mult[OF \<open>0 < sqrt xl\<close> \<open>sqrt xl \<le> x\<close> this]
  have "arctan x \ x * real_of_float (ub_arctan_horner p2 (get_odd n) 1 xl)" .
  ultimately show ?thesis by simp
qed

lemma arctan_0_1_bounds_round:
  assumes "0 \ real_of_float x" "real_of_float x \ 1"
  shows "arctan x \
      {real_of_float x * lb_arctan_horner p1 (get_even n) 1 (float_round_up (Suc p2) (x * x)) ..
        real_of_float x * ub_arctan_horner p3 (get_odd n) 1 (float_round_down (Suc p4) (x * x))}"
  using assms
  apply (cases "x > 0")
   apply (intro arctan_0_1_bounds_le)
   apply (auto simp: float_round_down.rep_eq float_round_up.rep_eq
    intro!: truncate_up_le1 mult_le_one truncate_down_le truncate_up_le truncate_down_pos
      mult_pos_pos)
  done

subsection "Compute \"

definition ub_pi :: "nat \ float" where
  "ub_pi prec =
    (let
      A = rapprox_rat prec 1 5 ;
      B = lapprox_rat prec 1 239
    in ((Float 1 2) * float_plus_up prec
      ((Float 1 2) * float_round_up prec (A * (ub_arctan_horner prec (get_odd (prec div 4 + 1)) 1
        (float_round_down (Suc prec) (A * A)))))
      (- float_round_down prec (B * (lb_arctan_horner prec (get_even (prec div 14 + 1)) 1
        (float_round_up (Suc prec) (B * B)))))))"

definition lb_pi :: "nat \ float" where
  "lb_pi prec =
    (let
      A = lapprox_rat prec 1 5 ;
      B = rapprox_rat prec 1 239
    in ((Float 1 2) * float_plus_down prec
      ((Float 1 2) * float_round_down prec (A * (lb_arctan_horner prec (get_even (prec div 4 + 1)) 1
        (float_round_up (Suc prec) (A * A)))))
      (- float_round_up prec (B * (ub_arctan_horner prec (get_odd (prec div 14 + 1)) 1
        (float_round_down (Suc prec) (B * B)))))))"

lemma pi_boundaries: "pi \ {(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 \ k" and "0 < k" and "1 \ k" by auto
    let ?k = "rapprox_rat prec 1 k"
    let ?kl = "float_round_down (Suc prec) (?k * ?k)"
    have "1 div k = 0" using div_pos_pos_trivial[OF _ \<open>1 < k\<close>] by auto

    have "0 \ real_of_float ?k" by (rule order_trans[OF _ rapprox_rat]) (auto simp add: \0 \ k\)
    have "real_of_float ?k \ 1"
      by (auto simp add: \<open>0 < k\<close> \<open>1 \<le> k\<close> less_imp_le
        intro!: mult_le_one order_trans[OF _ rapprox_rat] rapprox_rat_le1)
    have "1 / k \ ?k" using rapprox_rat[where x=1 and y=k] by auto
    hence "arctan (1 / k) \ arctan ?k" by (rule arctan_monotone')
    also have "\ \ (?k * ub_arctan_horner prec (get_odd n) 1 ?kl)"
      using arctan_0_1_bounds_round[OF \<open>0 \<le> real_of_float ?k\<close> \<open>real_of_float ?k \<le> 1\<close>]
      by auto
    finally have "arctan (1 / k) \ ?k * ub_arctan_horner prec (get_odd n) 1 ?kl" .
  } note ub_arctan = this

  {
    fix prec n :: nat
    fix k :: int
    assume "1 < k" hence "0 \ k" and "0 < k" by auto
    let ?k = "lapprox_rat prec 1 k"
    let ?ku = "float_round_up (Suc prec) (?k * ?k)"
    have "1 div k = 0" using div_pos_pos_trivial[OF _ \<open>1 < k\<close>] by auto
    have "1 / k \ 1" using \1 < k\ by auto
    have "0 \ real_of_float ?k" using lapprox_rat_nonneg[where x=1 and y=k, OF zero_le_one \0 \ k\]
      by (auto simp add: \<open>1 div k = 0\<close>)
    have "0 \ real_of_float (?k * ?k)" by simp
    have "real_of_float ?k \ 1" using lapprox_rat by (rule order_trans, auto simp add: \1 / k \ 1\)
    hence "real_of_float (?k * ?k) \ 1" using \0 \ real_of_float ?k\ by (auto intro!: mult_le_one)

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

    have "?k * lb_arctan_horner prec (get_even n) 1 ?ku \ arctan ?k"
      using arctan_0_1_bounds_round[OF \<open>0 \<le> real_of_float ?k\<close> \<open>real_of_float ?k \<le> 1\<close>]
      by auto
    also have "\ \ arctan (1 / k)" using \?k \ 1 / k\ by (rule arctan_monotone')
    finally have "?k * lb_arctan_horner prec (get_even n) 1 ?ku \ arctan (1 / k)" .
  } note lb_arctan = this

  have "pi \ ub_pi n "
    unfolding ub_pi_def machin_pi Let_def times_float.rep_eq Float_num
    using lb_arctan[of 239] ub_arctan[of 5] powr_realpow[of 2 2]
    by (intro mult_left_mono float_plus_up_le float_plus_down_le)
      (auto intro!: mult_left_mono float_round_down_le float_round_up_le diff_mono)
  moreover have "lb_pi n \ pi"
    unfolding lb_pi_def machin_pi Let_def times_float.rep_eq Float_num
    using lb_arctan[of 5] ub_arctan[of 239]
    by (intro mult_left_mono float_plus_up_le float_plus_down_le)
      (auto intro!: mult_left_mono float_round_down_le float_round_up_le diff_mono)
  ultimately show ?thesis by auto
qed

lift_definition pi_float_interval::"nat \ float interval" is "\prec. (lb_pi prec, ub_pi prec)"
  using pi_boundaries
  by (auto intro: order_trans)

lemma lower_pi_float_interval: "lower (pi_float_interval prec) = lb_pi prec"
  by transfer auto
lemma upper_pi_float_interval: "upper (pi_float_interval prec) = ub_pi prec"
  by transfer auto
lemma pi_float_interval: "pi \ set_of (real_interval (pi_float_interval prec))"
  using pi_boundaries
  by (auto simp: set_of_eq lower_pi_float_interval upper_pi_float_interval)

subsection "Compute arcus tangens in the entire domain"

function lb_arctan :: "nat \ float \ float" and ub_arctan :: "nat \ float \ float" where
  "lb_arctan prec x =
    (let
      ub_horner = \<lambda> x. float_round_up prec
        (x *
          ub_arctan_horner prec (get_odd (prec div 4 + 1)) 1 (float_round_down (Suc prec) (x * x)));
      lb_horner = \<lambda> x. float_round_down prec
        (x *
          lb_arctan_horner prec (get_even (prec div 4 + 1)) 1 (float_round_up (Suc prec) (x * x)))
    in
      if x < 0 then - ub_arctan prec (-x)
      else if x \<le> Float 1 (- 1) then lb_horner x
      else if x \<le> Float 1 1 then
        Float 1 1 *
        lb_horner
          (float_divl prec x
            (float_plus_up prec 1
              (ub_sqrt prec (float_plus_up prec 1 (float_round_up prec (x * x))))))
      else let inv = float_divr prec 1 x in
        if inv > 1 then 0
        else float_plus_down prec (lb_pi prec * Float 1 (- 1)) ( - ub_horner inv))"

| "ub_arctan prec x =
    (let
      lb_horner = \<lambda> x. float_round_down prec
        (x *
          lb_arctan_horner prec (get_even (prec div 4 + 1)) 1 (float_round_up (Suc prec) (x * x))) ;
      ub_horner = \<lambda> x. float_round_up prec
        (x *
          ub_arctan_horner prec (get_odd (prec div 4 + 1)) 1 (float_round_down (Suc prec) (x * x)))
    in if x < 0 then - lb_arctan prec (-x)
    else if x \<le> Float 1 (- 1) then ub_horner x
    else if x \<le> Float 1 1 then
      let y = float_divr prec x
        (float_plus_down
          (Suc prec) 1 (lb_sqrt prec (float_plus_down prec 1 (float_round_down prec (x * x)))))
      in if y > 1 then ub_pi prec * Float 1 (- 1) else Float 1 1 * ub_horner y
    else float_plus_up prec (ub_pi prec * Float 1 (- 1)) ( - lb_horner (float_divl prec 1 x)))"
by pat_completeness auto
termination
by (relation "measure (\ v. let (prec, x) = case_sum id id v in (if x < 0 then 1 else 0))", auto)

declare ub_arctan_horner.simps[simp del]
declare lb_arctan_horner.simps[simp del]

lemma lb_arctan_bound':
  assumes "0 \ real_of_float x"
  shows "lb_arctan prec x \ arctan x"
proof -
  have "\ x < 0" and "0 \ x"
    using \<open>0 \<le> real_of_float x\<close> by (auto intro!: truncate_up_le )

  let "?ub_horner x" =
      "x * ub_arctan_horner prec (get_odd (prec div 4 + 1)) 1 (float_round_down (Suc prec) (x * x))"
    and "?lb_horner x" =
      "x * lb_arctan_horner prec (get_even (prec div 4 + 1)) 1 (float_round_up (Suc prec) (x * x))"

  show ?thesis
  proof (cases "x \ Float 1 (- 1)")
    case True
    hence "real_of_float x \ 1" by simp
    from arctan_0_1_bounds_round[OF \<open>0 \<le> real_of_float x\<close> \<open>real_of_float x \<le> 1\<close>]
    show ?thesis
      unfolding lb_arctan.simps Let_def if_not_P[OF \<open>\<not> x < 0\<close>] if_P[OF True] using \<open>0 \<le> x\<close>
      by (auto intro!: float_round_down_le)
  next
    case False
    hence "0 < real_of_float x" by auto
    let ?R = "1 + sqrt (1 + real_of_float x * real_of_float x)"
    let ?sxx = "float_plus_up prec 1 (float_round_up prec (x * x))"
    let ?fR = "float_plus_up prec 1 (ub_sqrt prec ?sxx)"
    let ?DIV = "float_divl prec x ?fR"

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

    have "sqrt (1 + x*x) \ sqrt ?sxx"
      by (auto simp: float_plus_up.rep_eq plus_up_def float_round_up.rep_eq intro!: truncate_up_le)
    also have "\ \ ub_sqrt prec ?sxx"
      using bnds_sqrt'[of ?sxx prec] by auto
    finally
    have "sqrt (1 + x*x) \ ub_sqrt prec ?sxx" .
    hence "?R \ ?fR" by (auto simp: float_plus_up.rep_eq plus_up_def intro!: truncate_up_le)
    hence "0 < ?fR" and "0 < real_of_float ?fR" using \<open>0 < ?R\<close> by auto

    have monotone: "?DIV \ x / ?R"
    proof -
      have "?DIV \ real_of_float x / ?fR" by (rule float_divl)
      also have "\ \ x / ?R" by (rule divide_left_mono[OF \?R \ ?fR\ \0 \ real_of_float x\ mult_pos_pos[OF order_less_le_trans[OF divisor_gt0 \?R \ real_of_float ?fR\] divisor_gt0]])
      finally show ?thesis .
    qed

    show ?thesis
    proof (cases "x \ Float 1 1")
      case True
      have "x \ sqrt (1 + x * x)"
        using real_sqrt_sum_squares_ge2[where x=1, unfolded numeral_2_eq_2] by auto
      also note \<open>\<dots> \<le> (ub_sqrt prec ?sxx)\<close>
      finally have "real_of_float x \ ?fR"
        by (auto simp: float_plus_up.rep_eq plus_up_def intro!: truncate_up_le)
      moreover have "?DIV \ real_of_float x / ?fR"
        by (rule float_divl)
      ultimately have "real_of_float ?DIV \ 1"
        unfolding divide_le_eq_1_pos[OF \<open>0 < real_of_float ?fR\<close>, symmetric] by auto

      have "0 \ real_of_float ?DIV"
        using float_divl_lower_bound[OF \<open>0 \<le> x\<close>] \<open>0 < ?fR\<close>
        unfolding less_eq_float_def by auto

      from arctan_0_1_bounds_round[OF \<open>0 \<le> real_of_float (?DIV)\<close> \<open>real_of_float (?DIV) \<le> 1\<close>]
      have "Float 1 1 * ?lb_horner ?DIV \ 2 * arctan ?DIV"
        by simp
      also have "\ \ 2 * arctan (x / ?R)"
        using arctan_monotone'[OF monotone] by (auto intro!: mult_left_mono arctan_monotone')
      also have "2 * arctan (x / ?R) = arctan x"
        using arctan_half[symmetric] unfolding numeral_2_eq_2 power_Suc2 power_0 mult_1_left .
      finally show ?thesis
        unfolding lb_arctan.simps Let_def if_not_P[OF \<open>\<not> x < 0\<close>]
          if_not_P[OF \<open>\<not> x \<le> Float 1 (- 1)\<close>] if_P[OF True]
        by (auto simp: float_round_down.rep_eq
          intro!: order_trans[OF mult_left_mono[OF truncate_down]])
    next
      case False
      hence "2 < real_of_float x" by auto
      hence "1 \ real_of_float x" by auto

      let "?invx" = "float_divr prec 1 x"
      have "0 \ arctan x" using arctan_monotone'[OF \0 \ real_of_float 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 \<open>\<not> x < 0\<close>]
            if_not_P[OF \<open>\<not> x \<le> Float 1 (- 1)\<close>] if_not_P[OF False] if_P[OF True]
          using \<open>0 \<le> arctan x\<close> by auto
      next
        case False
        hence "real_of_float ?invx \ 1" by auto
        have "0 \ real_of_float ?invx"
          by (rule order_trans[OF _ float_divr]) (auto simp add: \<open>0 \<le> real_of_float x\<close>)

        have "1 / x \ 0" and "0 < 1 / x"
          using \<open>0 < real_of_float x\<close> by auto

        have "arctan (1 / x) \ arctan ?invx"
          unfolding one_float.rep_eq[symmetric] by (rule arctan_monotone', rule float_divr)
        also have "\ \ ?ub_horner ?invx"
          using arctan_0_1_bounds_round[OF \<open>0 \<le> real_of_float ?invx\<close> \<open>real_of_float ?invx \<le> 1\<close>]
          by (auto intro!: float_round_up_le)
        also note float_round_up
        finally have "pi / 2 - float_round_up prec (?ub_horner ?invx) \ arctan x"
          using \<open>0 \<le> arctan x\<close> arctan_inverse[OF \<open>1 / x \<noteq> 0\<close>]
          unfolding sgn_pos[OF \<open>0 < 1 / real_of_float x\<close>] le_diff_eq by auto
        moreover
        have "lb_pi prec * Float 1 (- 1) \ pi / 2"
          unfolding Float_num times_divide_eq_right mult_1_left using pi_boundaries by simp
        ultimately
        show ?thesis
          unfolding lb_arctan.simps Let_def if_not_P[OF \<open>\<not> x < 0\<close>]
            if_not_P[OF \<open>\<not> x \<le> Float 1 (- 1)\<close>] if_not_P[OF \<open>\<not> x \<le> Float 1 1\<close>] if_not_P[OF False]
          by (auto intro!: float_plus_down_le)
      qed
    qed
  qed
qed

lemma ub_arctan_bound':
  assumes "0 \ real_of_float x"
  shows "arctan x \ ub_arctan prec x"
proof -
  have "\ x < 0" and "0 \ x"
    using \<open>0 \<le> real_of_float x\<close> by auto

  let "?ub_horner x" =
    "float_round_up prec (x * ub_arctan_horner prec (get_odd (prec div 4 + 1)) 1 (float_round_down (Suc prec) (x * x)))"
  let "?lb_horner x" =
    "float_round_down prec (x * lb_arctan_horner prec (get_even (prec div 4 + 1)) 1 (float_round_up (Suc prec) (x * x)))"

  show ?thesis
  proof (cases "x \ Float 1 (- 1)")
    case True
    hence "real_of_float x \ 1" by auto
    show ?thesis
      unfolding ub_arctan.simps Let_def if_not_P[OF \<open>\<not> x < 0\<close>] if_P[OF True]
      using arctan_0_1_bounds_round[OF \<open>0 \<le> real_of_float x\<close> \<open>real_of_float x \<le> 1\<close>]
      by (auto intro!: float_round_up_le)
  next
    case False
    hence "0 < real_of_float x" by auto
    let ?R = "1 + sqrt (1 + real_of_float x * real_of_float x)"
    let ?sxx = "float_plus_down prec 1 (float_round_down prec (x * x))"
    let ?fR = "float_plus_down (Suc prec) 1 (lb_sqrt prec ?sxx)"
    let ?DIV = "float_divr prec x ?fR"

    have sqr_ge0: "0 \ 1 + real_of_float x * real_of_float x"
      using sum_power2_ge_zero[of 1 "real_of_float x", unfolded numeral_2_eq_2] by auto
    hence "0 \ real_of_float (1 + x*x)" by auto

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

    have "lb_sqrt prec ?sxx \ sqrt ?sxx"
      using bnds_sqrt'[of ?sxx] by auto
    also have "\ \ sqrt (1 + x*x)"
      by (auto simp: float_plus_down.rep_eq plus_down_def float_round_down.rep_eq truncate_down_le)
    finally have "lb_sqrt prec ?sxx \ sqrt (1 + x*x)" .
    hence "?fR \ ?R"
      by (auto simp: float_plus_down.rep_eq plus_down_def truncate_down_le)
    have "0 < real_of_float ?fR"
      by (auto simp: float_plus_down.rep_eq plus_down_def float_round_down.rep_eq
        intro!: truncate_down_ge1 lb_sqrt_lower_bound order_less_le_trans[OF zero_less_one]
        truncate_down_nonneg add_nonneg_nonneg)
    have monotone: "x / ?R \ (float_divr prec x ?fR)"
    proof -
      from divide_left_mono[OF \<open>?fR \<le> ?R\<close> \<open>0 \<le> real_of_float x\<close> mult_pos_pos[OF divisor_gt0 \<open>0 < real_of_float ?fR\<close>]]
      have "x / ?R \ x / ?fR" .
      also have "\ \ ?DIV" by (rule float_divr)
      finally show ?thesis .
    qed

    show ?thesis
    proof (cases "x \ Float 1 1")
      case True
      show ?thesis
      proof (cases "?DIV > 1")
        case True
        have "pi / 2 \ ub_pi prec * Float 1 (- 1)"
          unfolding Float_num times_divide_eq_right mult_1_left using pi_boundaries by auto
        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 \<open>\<not> x < 0\<close>]
            if_not_P[OF \<open>\<not> x \<le> Float 1 (- 1)\<close>] if_P[OF \<open>x \<le> Float 1 1\<close>] if_P[OF True] .
      next
        case False
        hence "real_of_float ?DIV \ 1" by auto

        have "0 \ x / ?R"
          using \<open>0 \<le> real_of_float x\<close> \<open>0 < ?R\<close> unfolding zero_le_divide_iff by auto
        hence "0 \ real_of_float ?DIV"
          using monotone by (rule order_trans)

        have "arctan x = 2 * arctan (x / ?R)"
          using arctan_half unfolding numeral_2_eq_2 power_Suc2 power_0 mult_1_left .
        also have "\ \ 2 * arctan (?DIV)"
          using arctan_monotone'[OF monotone] by (auto intro!: mult_left_mono)
        also have "\ \ (Float 1 1 * ?ub_horner ?DIV)" unfolding Float_num
          using arctan_0_1_bounds_round[OF \<open>0 \<le> real_of_float ?DIV\<close> \<open>real_of_float ?DIV \<le> 1\<close>]
          by (auto intro!: float_round_up_le)
        finally show ?thesis
          unfolding ub_arctan.simps Let_def if_not_P[OF \<open>\<not> x < 0\<close>]
            if_not_P[OF \<open>\<not> x \<le> Float 1 (- 1)\<close>] if_P[OF \<open>x \<le> Float 1 1\<close>] if_not_P[OF False] .
      qed
    next
      case False
      hence "2 < real_of_float x" by auto
      hence "1 \ real_of_float x" by auto
      hence "0 < real_of_float x" by auto
      hence "0 < x" by auto

      let "?invx" = "float_divl prec 1 x"
      have "0 \ arctan x"
        using arctan_monotone'[OF \0 \ real_of_float x\] and arctan_tan[of 0, unfolded tan_zero] by auto

      have "real_of_float ?invx \ 1"
        unfolding less_float_def
        by (rule order_trans[OF float_divl])
          (auto simp add: \<open>1 \<le> real_of_float x\<close> divide_le_eq_1_pos[OF \<open>0 < real_of_float x\<close>])
      have "0 \ real_of_float ?invx"
        using \<open>0 < x\<close> by (intro float_divl_lower_bound) auto

      have "1 / x \ 0" and "0 < 1 / x"
        using \<open>0 < real_of_float x\<close> by auto

      have "(?lb_horner ?invx) \ arctan (?invx)"
        using arctan_0_1_bounds_round[OF \<open>0 \<le> real_of_float ?invx\<close> \<open>real_of_float ?invx \<le> 1\<close>]
        by (auto intro!: float_round_down_le)
      also have "\ \ arctan (1 / x)"
        unfolding one_float.rep_eq[symmetric] by (rule arctan_monotone') (rule float_divl)
      finally have "arctan x \ pi / 2 - (?lb_horner ?invx)"
        using \<open>0 \<le> arctan x\<close> arctan_inverse[OF \<open>1 / x \<noteq> 0\<close>]
        unfolding sgn_pos[OF \<open>0 < 1 / x\<close>] le_diff_eq by auto
      moreover
      have "pi / 2 \ ub_pi prec * Float 1 (- 1)"
        unfolding Float_num times_divide_eq_right mult_1_right
        using pi_boundaries by auto
      ultimately
      show ?thesis
        unfolding ub_arctan.simps Let_def if_not_P[OF \<open>\<not> x < 0\<close>]
          if_not_P[OF \<open>\<not> x \<le> Float 1 (- 1)\<close>] if_not_P[OF False]
        by (auto intro!: float_round_up_le float_plus_up_le)
    qed
  qed
qed

lemma arctan_boundaries: "arctan x \ {(lb_arctan prec x) .. (ub_arctan prec x)}"
proof (cases "0 \ x")
  case True
  hence "0 \ real_of_float x" by auto
  show ?thesis
    using ub_arctan_bound'[OF \0 \ real_of_float x\] lb_arctan_bound'[OF \0 \ real_of_float x\]
    unfolding atLeastAtMost_iff by auto
next
  case False
  let ?mx = "-x"
  from False have "x < 0" and "0 \ real_of_float ?mx"
    by auto
  hence bounds: "lb_arctan prec ?mx \ arctan ?mx \ arctan ?mx \ ub_arctan prec ?mx"
    using ub_arctan_bound'[OF \0 \ real_of_float ?mx\] lb_arctan_bound'[OF \0 \ real_of_float ?mx\] by auto
  show ?thesis
    unfolding minus_float.rep_eq arctan_minus lb_arctan.simps[where x=x]
      ub_arctan.simps[where x=x] Let_def if_P[OF \<open>x < 0\<close>]
    unfolding atLeastAtMost_iff using bounds[unfolded minus_float.rep_eq arctan_minus]
    by (simp add: arctan_minus)
qed

lemma bnds_arctan: "\ (x::real) lx ux. (l, u) = (lb_arctan prec lx, ub_arctan prec ux) \ x \ {lx .. ux} \ l \ arctan x \ arctan x \ u"
proof (rule allI, rule allI, rule allI, rule impI)
  fix x :: real
  fix lx ux
  assume "(l, u) = (lb_arctan prec lx, ub_arctan prec ux) \ x \ {lx .. ux}"
  hence l: "lb_arctan prec lx = l "
    and u: "ub_arctan prec ux = u"
    and x: "x \ {lx .. ux}"
    by auto
  show "l \ arctan x \ arctan x \ u"
  proof
    show "l \ arctan x"
    proof -
      from arctan_boundaries[of lx prec, unfolded l]
      have "l \ arctan lx" by (auto simp del: lb_arctan.simps)
      also have "\ \ arctan x" using x by (auto intro: arctan_monotone')
      finally show ?thesis .
    qed
    show "arctan x \ u"
    proof -
      have "arctan x \ arctan ux" using x by (auto intro: arctan_monotone')
      also have "\ \ u" using arctan_boundaries[of ux prec, unfolded u] by (auto simp del: ub_arctan.simps)
      finally show ?thesis .
    qed
  qed
qed

lemmas [simp del] = lb_arctan.simps ub_arctan.simps

lemma lb_arctan: "arctan (real_of_float x) \ y \ real_of_float (lb_arctan prec x) \ y"
  and ub_arctan: "y \ arctan x \ y \ ub_arctan prec x"
  for x::float and y::real
  using arctan_boundaries[of x prec] by auto

lift_definition arctan_float_interval :: "nat \ float interval \ float interval"
  is "\prec. \(lx, ux). (lb_arctan prec lx, ub_arctan prec ux)"
  by (auto intro!: lb_arctan ub_arctan arctan_monotone')

lemma lower_arctan_float_interval: "lower (arctan_float_interval p x) = lb_arctan p (lower x)"
  by transfer auto
lemma upper_arctan_float_interval: "upper (arctan_float_interval p x) = ub_arctan p (upper x)"
  by transfer auto

lemma arctan_float_interval:
  "arctan ` set_of (real_interval x) \ set_of (real_interval (arctan_float_interval p x))"
  by (auto simp: set_of_eq lower_arctan_float_interval upper_arctan_float_interval
      intro!: lb_arctan ub_arctan arctan_monotone')

lemma arctan_float_intervalI:
  "arctan x \\<^sub>r arctan_float_interval p X" if "x \\<^sub>r X"
  using arctan_float_interval[of X p] that
  by auto

section "Sinus and Cosinus"

subsection "Compute the cosinus and sinus series"

fun ub_sin_cos_aux :: "nat \ nat \ nat \ nat \ float \ float"
and lb_sin_cos_aux :: "nat \ nat \ nat \ nat \ float \ float" where
  "ub_sin_cos_aux prec 0 i k x = 0"
| "ub_sin_cos_aux prec (Suc n) i k x = float_plus_up prec
    (rapprox_rat prec 1 k) (-
      float_round_down prec (x * (lb_sin_cos_aux prec n (i + 2) (k * i * (i + 1)) x)))"
| "lb_sin_cos_aux prec 0 i k x = 0"
| "lb_sin_cos_aux prec (Suc n) i k x = float_plus_down prec
    (lapprox_rat prec 1 k) (-
      float_round_up prec (x * (ub_sin_cos_aux prec n (i + 2) (k * i * (i + 1)) x)))"

lemma cos_aux:
  shows "(lb_sin_cos_aux prec n 1 1 (x * x)) \ (\ i=0..
  and "(\ i=0.. (ub_sin_cos_aux prec n 1 1 (x * x))" (is "?ub")
proof -
  have "0 \ real_of_float (x * x)" by auto
  let "?f n" = "fact (2 * n) :: nat"
  have f_eq: "?f (Suc n) = ?f n * ((\i. i + 2) ^^ n) 1 * (((\i. i + 2) ^^ n) 1 + 1)" for n
  proof -
    have "\m. ((\i. i + 2) ^^ n) m = m + 2 * n" by (induct n) auto
    then show ?thesis by auto
  qed
  from horner_bounds[where lb="lb_sin_cos_aux prec" and ub="ub_sin_cos_aux prec" and j'=0,
    OF \<open>0 \<le> real_of_float (x * x)\<close> 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_of_float x"])
qed

lemma lb_sin_cos_aux_zero_le_one: "lb_sin_cos_aux prec n i j 0 \ 1"
  by (cases j n rule: nat.exhaust[case_product nat.exhaust])
    (auto intro!: float_plus_down_le order_trans[OF lapprox_rat])

lemma one_le_ub_sin_cos_aux: "odd n \ 1 \ ub_sin_cos_aux prec n i (Suc 0) 0"
  by (cases n) (auto intro!: float_plus_up_le order_trans[OF _ rapprox_rat])

lemma cos_boundaries:
  assumes "0 \ real_of_float x" and "x \ pi / 2"
  shows "cos x \ {(lb_sin_cos_aux prec (get_even n) 1 1 (x * x)) .. (ub_sin_cos_aux prec (get_odd n) 1 1 (x * x))}"
proof (cases "real_of_float x = 0")
  case False
  hence "real_of_float x \ 0" by auto
  hence "0 < x" and "0 < real_of_float x"
    using \<open>0 \<le> real_of_float x\<close> by auto
  have "0 < x * x"
    using \<open>0 < x\<close> by simp

  have morph_to_if_power: "(\ i=0..
    (\<Sum> i = 0 ..< 2 * n. (if even(i) then ((- 1) ^ (i div 2))/((fact i)) else 0) * x ^ i)"
    (is "?sum = ?ifsum") for x n
  proof -
    have "?sum = ?sum + (\ j = 0 ..< n. 0)" by auto
    also have "\ =
      (\<Sum> j = 0 ..< n. (- 1) ^ ((2 * j) div 2) / ((fact (2 * j))) * x ^(2 * j)) + (\<Sum> j = 0 ..< n. 0)" by auto
    also have "\ = (\ i = 0 ..< 2 * n. if even i then (- 1) ^ (i div 2) / ((fact i)) * x ^ i else 0)"
      unfolding sum_split_even_odd atLeast0LessThan ..
    also have "\ = (\ i = 0 ..< 2 * n. (if even i then (- 1) ^ (i div 2) / ((fact i)) else 0) * x ^ i)"
      by (rule sum.cong) auto
    finally show ?thesis .
  qed

  { fix n :: nat assume "0 < n"
    hence "0 < 2 * n" by auto
    obtain t where "0 < t" and "t < real_of_float x" and
      cos_eq: "cos x = (\ i = 0 ..< 2 * n. (if even(i) then ((- 1) ^ (i div 2))/((fact i)) else 0) * (real_of_float x) ^ i)
      + (cos (t + 1/2 * (2 * n) * pi) / (fact (2*n))) * (real_of_float x)^(2*n)"
      (is "_ = ?SUM + ?rest / ?fact * ?pow")
      using Maclaurin_cos_expansion2[OF \<open>0 < real_of_float x\<close> \<open>0 < 2 * n\<close>]
      unfolding cos_coeff_def atLeast0LessThan by auto

    have "cos t * (- 1) ^ n = cos t * cos (n * pi) + sin t * sin (n * pi)" by auto
    also have "\ = cos (t + n * pi)" by (simp add: cos_add)
    also have "\ = ?rest" by auto
    finally have "cos t * (- 1) ^ n = ?rest" .
    moreover
    have "t \ pi / 2" using \t < real_of_float x\ and \x \ pi / 2\ by auto
    hence "0 \ cos t" using \0 < t\ and cos_ge_zero by auto
    ultimately have even: "even n \ 0 \ ?rest" and odd: "odd n \ 0 \ - ?rest " by auto

    have "0 < ?fact" by auto
    have "0 < ?pow" using \<open>0 < real_of_float x\<close> by auto

    {
      assume "even n"
      have "(lb_sin_cos_aux prec n 1 1 (x * x)) \ ?SUM"
        unfolding morph_to_if_power[symmetric] using cos_aux by auto
      also have "\ \ cos x"
      proof -
        from even[OF \<open>even n\<close>] \<open>0 < ?fact\<close> \<open>0 < ?pow\<close>
        have "0 \ (?rest / ?fact) * ?pow" by simp
        thus ?thesis unfolding cos_eq by auto
      qed
      finally have "(lb_sin_cos_aux prec n 1 1 (x * x)) \ cos x" .
    } note lb = this

    {
      assume "odd n"
      have "cos x \ ?SUM"
      proof -
        from \<open>0 < ?fact\<close> and \<open>0 < ?pow\<close> and odd[OF \<open>odd n\<close>]
        have "0 \ (- ?rest) / ?fact * ?pow"
          by (metis mult_nonneg_nonneg divide_nonneg_pos less_imp_le)
        thus ?thesis unfolding cos_eq by auto
      qed
      also have "\ \ (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 \ (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 \ (ub_sin_cos_aux prec (get_odd n) 1 1 (x * x))"
    using ub[OF odd_pos[OF get_odd] get_odd] .
  moreover have "(lb_sin_cos_aux prec (get_even n) 1 1 (x * x)) \ cos x"
  proof (cases "0 < get_even n")
    case True
    show ?thesis using lb[OF True get_even] .
  next
    case False
    hence "get_even n = 0" by auto
    have "- (pi / 2) \ x"
      by (rule order_trans[OF _ \<open>0 < real_of_float x\<close>[THEN less_imp_le]]) auto
    with \<open>x \<le> pi / 2\<close> show ?thesis
      unfolding \<open>get_even n = 0\<close> lb_sin_cos_aux.simps minus_float.rep_eq zero_float.rep_eq
      using cos_ge_zero by auto
  qed
  ultimately show ?thesis by auto
next
  case True
  hence "x = 0"
    by (simp add: real_of_float_eq)
  thus ?thesis
    using lb_sin_cos_aux_zero_le_one one_le_ub_sin_cos_aux
    by simp
qed

lemma sin_aux:
  assumes "0 \ real_of_float x"
  shows "(x * lb_sin_cos_aux prec n 2 1 (x * x)) \
      (\<Sum> i=0..<n. (- 1) ^ i * (1/(fact (2 * i + 1))) * x^(2 * i + 1))" (is "?lb")
    and "(\ i=0..
      (x * ub_sin_cos_aux prec n 2 1 (x * x))" (is "?ub")
proof -
  have "0 \ real_of_float (x * x)" by auto
  let "?f n" = "fact (2 * n + 1) :: nat"
  have f_eq: "?f (Suc n) = ?f n * ((\i. i + 2) ^^ n) 2 * (((\i. i + 2) ^^ n) 2 + 1)" for n
  proof -
    have F: "\m. ((\i. i + 2) ^^ n) m = m + 2 * n" by (induct n) auto
    show ?thesis
      unfolding F by auto
  qed
  from horner_bounds[where lb="lb_sin_cos_aux prec" and ub="ub_sin_cos_aux prec" and j'=0,
    OF \<open>0 \<le> real_of_float (x * x)\<close> f_eq lb_sin_cos_aux.simps ub_sin_cos_aux.simps]
  show "?lb" and "?ub" using \<open>0 \<le> real_of_float x\<close>
    apply (simp_all only: power_add power_one_right mult.assoc[symmetric] sum_distrib_right[symmetric])
    apply (simp_all only: mult.commute[where 'a=real] of_nat_fact)
    apply (auto intro!: mult_left_mono simp add: power_mult power2_eq_square[of "real_of_float x"])
    done
qed

lemma sin_boundaries:
  assumes "0 \ real_of_float x"
    and "x \ pi / 2"
  shows "sin x \ {(x * lb_sin_cos_aux prec (get_even n) 2 1 (x * x)) .. (x * ub_sin_cos_aux prec (get_odd n) 2 1 (x * x))}"
proof (cases "real_of_float x = 0")
  case False
  hence "real_of_float x \ 0" by auto
  hence "0 < x" and "0 < real_of_float x"
    using \<open>0 \<le> real_of_float x\<close> by auto
  have "0 < x * x"
    using \<open>0 < x\<close> by simp

  have sum_morph: "(\j = 0 ..< n. (- 1) ^ (((2 * j + 1) - Suc 0) div 2) / ((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))/((fact i))) * x ^ i)"
    (is "?SUM = _") for x :: real and n
  proof -
    have pow: "!!i. x ^ (2 * i + 1) = x * x ^ (2 * i)"
      by auto
    have "?SUM = (\ j = 0 ..< n. 0) + ?SUM"
      by auto
    also have "\ = (\ i = 0 ..< 2 * n. if even i then 0 else (- 1) ^ ((i - Suc 0) div 2) / ((fact i)) * x ^ i)"
      unfolding sum_split_even_odd atLeast0LessThan ..
    also have "\ = (\ i = 0 ..< 2 * n. (if even i then 0 else (- 1) ^ ((i - Suc 0) div 2) / ((fact i))) * x ^ i)"
      by (rule sum.cong) auto
    finally show ?thesis .
  qed

  { fix n :: nat assume "0 < n"
    hence "0 < 2 * n + 1" by auto
    obtain t where "0 < t" and "t < real_of_float x" and
      sin_eq: "sin x = (\ i = 0 ..< 2 * n + 1. (if even(i) then 0 else ((- 1) ^ ((i - Suc 0) div 2))/((fact i))) * (real_of_float x) ^ i)
      + (sin (t + 1/2 * (2 * n + 1) * pi) / (fact (2*n + 1))) * (real_of_float x)^(2*n + 1)"
      (is "_ = ?SUM + ?rest / ?fact * ?pow")
      using Maclaurin_sin_expansion3[OF \<open>0 < 2 * n + 1\<close> \<open>0 < real_of_float x\<close>]
      unfolding sin_coeff_def atLeast0LessThan by auto

    have "?rest = cos t * (- 1) ^ n"
      unfolding sin_add cos_add of_nat_add distrib_right distrib_left by auto
    moreover
    have "t \ pi / 2"
      using \<open>t < real_of_float x\<close> and \<open>x \<le> pi / 2\<close> by auto
    hence "0 \ cos t"
      using \<open>0 < t\<close> and cos_ge_zero by auto
    ultimately have even: "even n \ 0 \ ?rest" and odd: "odd n \ 0 \ - ?rest"
      by auto

    have "0 < ?fact"
      by (simp del: fact_Suc)
    have "0 < ?pow"
      using \<open>0 < real_of_float x\<close> by (rule zero_less_power)

    {
      assume "even n"
      have "(x * lb_sin_cos_aux prec n 2 1 (x * x)) \
            (\<Sum> i = 0 ..< 2 * n. (if even(i) then 0 else ((- 1) ^ ((i - Suc 0) div 2))/((fact i))) * (real_of_float x) ^ i)"
        using sin_aux[OF \<open>0 \<le> real_of_float x\<close>] unfolding sum_morph[symmetric] by auto
      also have "\ \ ?SUM" by auto
      also have "\ \ sin x"
      proof -
        from even[OF \<open>even n\<close>] \<open>0 < ?fact\<close> \<open>0 < ?pow\<close>
        have "0 \ (?rest / ?fact) * ?pow" by simp
        thus ?thesis unfolding sin_eq by auto
      qed
      finally have "(x * lb_sin_cos_aux prec n 2 1 (x * x)) \ sin x" .
    } note lb = this

    {
      assume "odd n"
      have "sin x \ ?SUM"
      proof -
        from \<open>0 < ?fact\<close> and \<open>0 < ?pow\<close> and odd[OF \<open>odd n\<close>]
        have "0 \ (- ?rest) / ?fact * ?pow"
          by (metis mult_nonneg_nonneg divide_nonneg_pos less_imp_le)
        thus ?thesis unfolding sin_eq by auto
      qed
      also have "\ \ (\ i = 0 ..< 2 * n. (if even(i) then 0 else ((- 1) ^ ((i - Suc 0) div 2))/((fact i))) * (real_of_float x) ^ i)"
         by auto
--> --------------------

--> maximum size reached

--> --------------------

¤ Dauer der Verarbeitung: 0.28 Sekunden (vorverarbeitet) ¤

Download des Quellennavigators
Download des sprechenden Kalenders
in der Quellcodebibliothek suchen

Haftungshinweis

Die Informationen auf dieser Webseite wurden nach bestem Wissen sorgfältig zusammengestellt. Es wird jedoch weder Vollständigkeit, noch Richtigkeit, noch Qualität der bereit gestellten Informationen zugesichert.

Bemerkung:

Die farbliche Syntaxdarstellung ist noch experimentell.