Theory Euler_MacLaurin.Euler_MacLaurin

section ‹The Euler--MacLaurin summation formula›
theory Euler_MacLaurin
imports 
  "HOL-Complex_Analysis.Complex_Analysis"
  Bernoulli.Periodic_Bernpoly
  Bernoulli.Bernoulli_FPS
begin
  
subsection ‹Auxiliary facts›

(* TODO Move? *)
lemma pbernpoly_of_int [simp]: "pbernpoly n (of_int a) = bernoulli n"
  by (simp add: pbernpoly_def)

lemma continuous_on_bernpoly' [continuous_intros]:
  assumes "continuous_on A f"
  shows   "continuous_on A (λx. bernpoly n (f x) :: 'a :: real_normed_algebra_1)"
  using continuous_on_compose2[OF continuous_on_bernpoly assms, of UNIV n] by auto

lemma sum_atLeastAtMost_int_last:
  assumes "a < (b :: int)"
  shows   "sum f {a..b} = sum f {a..<b} + f b"
proof -
  from assms have "{a..b} = insert b {a..<b}" by auto
  also have "sum f  = sum f {a..<b} + f b"
    by (subst sum.insert) (auto simp: add_ac)
  finally show ?thesis .
qed

lemma sum_atLeastAtMost_int_head:
  assumes "a < (b :: int)"
  shows   "sum f {a..b} = f a + sum f {a<..b}"
proof -
  from assms have "{a..b} = insert a {a<..b}" by auto
  also have "sum f  = f a + sum f {a<..b}"
    by (subst sum.insert) auto
  finally show ?thesis .
qed

lemma not_in_nonpos_Reals_imp_add_nonzero:
  assumes "z  0" "x  0"
  shows   "z + of_real x  0"
  using assms by (auto simp: add_eq_0_iff2)
(* END TODO *)

lemma negligible_atLeastAtMostI: "b  a  negligible {a..(b::real)}"
  by (cases "b < a") auto
    
lemma integrable_on_negligible: 
 "negligible A  (f :: 'n :: euclidean_space  'a :: banach) integrable_on A"
  by (subst integrable_spike_set_eq[of _ "{}"]) (simp_all add: integrable_on_empty)  
    
lemma Union_atLeastAtMost_real_of_int:
  assumes "a < b" 
  shows   "(n{a..<b}. {real_of_int n..real_of_int (n + 1)}) = {real_of_int a..real_of_int b}"
proof (intro equalityI subsetI)
  fix x assume x: "x  {real_of_int a..real_of_int b}"
  thus "x  (n{a..<b}. {real_of_int n..real_of_int (n + 1)})"
  proof (cases "x = real_of_int b")
    case True
    with assms show ?thesis by (auto intro!: bexI[of _ "b - 1"])
  next
    case False
    with x have x: "x  real_of_int a" "x < real_of_int b" by simp_all
    hence "x  of_int x" "x  of_int x + 1" by linarith+
    moreover from x have "x  a" "x < b" by linarith+ 
    ultimately have "n{a..<b}. x  {of_int n..of_int (n + 1)}"
      by (intro bexI[of _ "x"]) simp_all
    thus ?thesis by blast
  qed
qed auto


subsection ‹The remainder terms›

text ‹
  The following describes the remainder term in the classical version of the Euler--MacLaurin
  formula.
›
definition EM_remainder' :: "nat  (real  'a :: banach)  real  real  'a" where
  "EM_remainder' n f a b = ((-1) ^ Suc n / fact n) *R integral {a..b} (λt. pbernpoly n t *R f t)"


text ‹
  Next, we define the remainder term that occurs when one lets the right bound of summation 
  in the Euler--MacLaurin formula tend to infinity.
›
definition EM_remainder_converges :: "nat  (real  'a :: banach)  int  bool" where
  "EM_remainder_converges n f a  (L. ((λx. EM_remainder' n f a (of_int x))  L) at_top)"

definition EM_remainder :: "nat  (real  'a :: banach)  int  'a" where
  "EM_remainder n f a = 
     (if EM_remainder_converges n f a then
        Lim at_top (λx. EM_remainder' n f a (of_int x)) else 0)"

text ‹
  The following lemmas are fairly obvious -- but tedious to prove -- 
  properties of the remainder terms.
›

lemma EM_remainder_eqI: 
  fixes L
  assumes "((λx. EM_remainder' n f b (of_int x))  L) at_top"
  shows   "EM_remainder n f b = L"
  using assms by (auto simp: EM_remainder_def EM_remainder_converges_def intro!: tendsto_Lim)

lemma integrable_EM_remainder'_int:
  fixes a b :: int and f :: "real  'a :: banach"
  assumes "continuous_on {of_int a..of_int b} f"
  shows   "(λt. pbernpoly n t *R f t) integrable_on {a..b}"
proof -
  have [continuous_intros]: "continuous_on A f" if "A  {of_int a..of_int b}" for A
    using continuous_on_subset[OF assms that] .
  consider "a > b" | "a = b" | "a < b" "n = 1" | "a < b" "n  1"
    by (cases a b rule: linorder_cases) auto
  thus ?thesis
  proof cases
    assume "a < b" and "n  1"
    thus ?thesis by (intro integrable_continuous_real continuous_intros) auto
  next
    assume ab: "a < b" and [simp]: "n = 1"
    let ?A = "(λn. {real_of_int n..real_of_int (n+1)}) ` {a..<b}"
    show ?thesis
    proof (rule integrable_combine_division; (intro ballI)?)
      show "?A division_of {of_int a..of_int b}"
        using Union_atLeastAtMost_real_of_int[OF ab] by (simp add: division_of_def)
    next
      fix I assume "I  ?A"
      then obtain i where i: "i  {a..<b}" "I = {of_int i..of_int (i + 1)}" by auto
      show "(λt. pbernpoly n t *R f t) integrable_on I"
      proof (rule integrable_spike)
        show "(λt. (t - of_int i - 1/2) *R f t) integrable_on I"
          using i by (auto intro!: integrable_continuous_real continuous_intros)
      next
        fix x assume "x  I - {of_int (i + 1)}"
        with i have "of_int i  x" "x < of_int i + 1" by simp_all
        hence "floor x = i" by linarith
        thus "pbernpoly n x *R f x = (x - of_int i - 1 / 2) *R f x"
          by (simp add: pbernpoly_def bernpoly_def frac_def)
      qed simp_all
    qed
  qed (simp_all add: integrable_on_negligible)
qed

lemma integrable_EM_remainder':
  fixes a b :: real and f :: "real  'a :: banach"
  assumes "continuous_on {a..b} f"
  shows   "(λt. pbernpoly n t *R f t) integrable_on {a..b}"
proof (cases "a  b")
  case True
  define a' b' where "a' = a" and "b' = b"
  from True have *: "a'  b'" "a'  a" "b'  b" by (auto simp: a'_def b'_def)
  from * have A: "(λt. pbernpoly n t *R f t) integrable_on ({a'..b'})"
    by (intro integrable_EM_remainder'_int continuous_on_subset[OF assms]) auto
  have B: "(λt. pbernpoly n t *R f t) integrable_on ({a..a'})"
  proof (rule integrable_spike)
    show "pbernpoly n x *R f x = bernpoly n (x - of_int (floor a)) *R f x" 
         if "x  {a..real_of_int a'} - {real_of_int a'}"for x
    proof -
      have "x  a" "x <real_of_int a'" using that by auto
      with True have "floor x = floor a" unfolding a'_def
        using ceiling_diff_floor_le_1[of a] by (intro floor_unique; linarith) 
      thus ?thesis by (simp add: pbernpoly_def frac_def)
    qed
  qed (insert *, auto intro!: continuous_intros integrable_continuous_real 
         continuous_on_subset[OF assms])
  have C: "(λt. pbernpoly n t *R f t) integrable_on ({b'..b})"
  proof (rule integrable_spike)
    show "pbernpoly n x *R f x = bernpoly n (x - of_int b') *R f x"
         if "x  {real_of_int b'..b} - {real_of_int b'}" for x
    proof -
      have "x  b" "x > real_of_int b'" using that by auto
      with True have "floor x = b'" unfolding b'_def by (intro floor_unique; linarith) 
      thus ?thesis by (simp add: pbernpoly_def frac_def)
    qed
  qed (insert *, auto intro!: continuous_intros integrable_continuous_real 
         continuous_on_subset[OF assms])
  have "(λt. pbernpoly n t *R f t) integrable_on ({a..a'}  {a'..b'}  {b'..b})" using * A B C
    by (intro integrable_Un; (subst ivl_disj_un)?)
       (auto simp: ivl_disj_un max_def min_def)
  also have "{a..a'}  {a'..b'}  {b'..b} = {a..b}" using * by auto
  finally show ?thesis .
next
  assume *: "¬ceiling a  floor b"
  show ?thesis
  proof (rule integrable_spike)
    show "(λt. bernpoly n (t - floor a) *R f t) integrable_on {a..b}" using *
      by (auto intro!: integrable_continuous_real continuous_intros assms)
  next
    show "pbernpoly n x *R f x = bernpoly n (x - floor a) *R f x"
         if "x  {a..b} - {}" for x
    proof -
      from * have **: "b < floor a + 1" 
        unfolding ceiling_altdef by (auto split: if_splits simp: le_floor_iff)
      from that have x: "x  a" "x  b" by simp_all
      with * ** have "floor x = floor a" by linarith
      thus ?thesis by (simp add: pbernpoly_def frac_def)
    qed
  qed simp_all
qed

lemma EM_remainder'_bounded_linear:
  assumes "bounded_linear h"
  assumes "continuous_on {a..b} f"
  shows   "EM_remainder' n (λx. h (f x)) a b = h (EM_remainder' n f a b)"
proof -
  have "integral {a..b} (λt. pbernpoly n t *R h (f t)) = 
           integral {a..b} (λt. h (pbernpoly n t *R f t))" using assms
    by (simp add: linear_simps)
  also have " = h (integral {a..b} (λt. pbernpoly n t *R f t))"
    by (subst integral_linear [OF _ assms(1), symmetric])
       (auto intro!: integrable_EM_remainder' assms(2) simp: o_def)
  finally show ?thesis using assms(1)
    by (simp add: EM_remainder'_def linear_simps)
qed

lemma EM_remainder_converges_of_real:
  assumes "EM_remainder_converges n f a" "continuous_on {of_int a..} f"
  shows   "EM_remainder_converges n (λx. of_real (f x)) a"
proof -
  from assms obtain L 
    where L: "((λb. EM_remainder' n f (real_of_int a) (real_of_int b))  L) at_top"
    by (auto simp: EM_remainder_converges_def)
  have "((λb. EM_remainder' n (λx. of_real (f x)) (of_int a) (of_int b))  of_real L) at_top"
  proof (rule Lim_transform_eventually)
    show "eventually (λb. of_real (EM_remainder' n f (of_int a) (of_int b)) = 
            EM_remainder' n (λx. of_real (f x)) (of_int a) (of_int b)) at_top"
      using eventually_ge_at_top[of a]
      by eventually_elim
         (intro EM_remainder'_bounded_linear [OF bounded_linear_of_real, symmetric]
                continuous_on_subset[OF assms(2)], auto)
  qed (intro tendsto_intros L)
  thus ?thesis unfolding EM_remainder_converges_def ..
qed

lemma EM_remainder_converges_of_real_iff:
  fixes f :: "real  real"
  assumes "continuous_on {of_int a..} f"
  shows   "EM_remainder_converges n (λx. of_real (f x) :: 
            'a :: {banach, real_normed_algebra_1, real_inner}) a  EM_remainder_converges n f a"
proof
  assume "EM_remainder_converges n (λx. of_real (f x) :: 'a) a"
  then obtain L :: 'a
    where L: "((λb. EM_remainder' n (λx. of_real (f x)) (of_int a) (of_int b))  L) at_top"
    by (auto simp: EM_remainder_converges_def)
  have "((λb. EM_remainder' n f (of_int a) (of_int b))  L  1) at_top"
  proof (rule Lim_transform_eventually)
    show "eventually (λb. EM_remainder' n (λx. of_real (f x) :: 'a) (of_int a) (of_int b)  1 =
            EM_remainder' n f (of_int a) (of_int b)) at_top" using eventually_ge_at_top[of a]
      by eventually_elim
         (subst EM_remainder'_bounded_linear [OF bounded_linear_of_real],
          auto intro!: continuous_on_subset[OF assms])
  qed (intro tendsto_intros L)
  thus "EM_remainder_converges n f a" unfolding EM_remainder_converges_def ..
qed (intro EM_remainder_converges_of_real assms)

lemma EM_remainder_of_real:
  assumes "continuous_on {a..} f"
  shows   "EM_remainder n (λx. of_real (f x) :: 
             'a :: {banach, real_normed_algebra_1, real_inner}) a = 
             of_real (EM_remainder n f a)"
proof -
  have eq: "EM_remainder' n (λx. of_real (f x) :: 'a) (real_of_int a) = 
              (λx::int. of_real (EM_remainder' n f a x))"
      by (intro ext EM_remainder'_bounded_linear[OF bounded_linear_of_real]
            continuous_on_subset[OF assms]) auto
  show ?thesis
  proof (cases "EM_remainder_converges n f a")
    case False
    with EM_remainder_converges_of_real_iff[OF assms, of n] show ?thesis
      by (auto simp: EM_remainder_def)
  next
    case True
    then obtain L where L: "((λx. EM_remainder' n f a (real_of_int x))  L) at_top" 
      by (auto simp: EM_remainder_converges_def)
    have L': "((λx. EM_remainder' n (λx. of_real (f x) :: 'a) a 
                (real_of_int x))  of_real L) at_top" unfolding eq by (intro tendsto_of_real L)
    from L L' tendsto_Lim[OF _ L] tendsto_Lim[OF _ L'] show ?thesis
      by (auto simp: EM_remainder_def EM_remainder_converges_def)
  qed
qed

lemma EM_remainder'_cong:
  assumes "x. x  {a..b}  f x = g x" "n = n'" "a = a'" "b = b'"
  shows   "EM_remainder' n f a b = EM_remainder' n' g a' b'"
proof -
  have "integral {a..b} (λt. pbernpoly n t *R f t) = integral {a'..b'} (λt. pbernpoly n' t *R g t)"
    unfolding assms using assms by (intro integral_cong) auto
  with assms show ?thesis by (simp add: EM_remainder'_def)
qed

lemma EM_remainder_converges_cong:
  assumes "x. x  of_int a  f x = g x" "n = n'" "a = a'"
  shows   "EM_remainder_converges n f a = EM_remainder_converges n' g a'"
  unfolding EM_remainder_converges_def
  by (subst EM_remainder'_cong[OF _ refl refl refl, of _ _ f g]) (use assms in auto)

lemma EM_remainder_cong:
  assumes "x. x  of_int a  f x = g x" "n = n'" "a = a'"
  shows   "EM_remainder n f a = EM_remainder n' g a'"
proof -
  have *: "EM_remainder_converges n f a = EM_remainder_converges n' g a'"
    using assms by (intro EM_remainder_converges_cong) auto
  show ?thesis unfolding EM_remainder_def
    by (subst EM_remainder'_cong[OF _ refl refl refl, of _ _ f g]) (use assms * in auto)
qed

lemma EM_remainder_converges_cnj: 
  assumes "continuous_on {a..} f" and "EM_remainder_converges n f a"
  shows   "EM_remainder_converges n (λx. cnj (f x)) a"
proof -
  interpret bounded_linear cnj by (rule bounded_linear_cnj)
  obtain L where L: "((λx. EM_remainder' n f (real_of_int a) (real_of_int x))  L) at_top"
    using assms unfolding EM_remainder_converges_def by blast
  note tendsto_cnj [OF this]
  also have "(λx. cnj (EM_remainder' n f (real_of_int a) (real_of_int x))) =
               (λx. EM_remainder' n (λx. cnj (f x)) (real_of_int a) (real_of_int x))"
    by (subst EM_remainder'_bounded_linear [OF bounded_linear_cnj])
       (rule continuous_on_subset [OF assms(1)], auto)
  finally have L': "(  cnj L) at_top" .
  thus "EM_remainder_converges n (λx. cnj (f x)) a"
    by (auto simp: EM_remainder_converges_def)
qed

lemma EM_remainder_converges_cnj_iff:
  assumes "continuous_on {of_int a..} f"
  shows   "EM_remainder_converges n (λx. cnj (f x)) a  EM_remainder_converges n f a"
proof
  assume "EM_remainder_converges n (λx. cnj (f x)) a"
  hence "EM_remainder_converges n (λx. cnj (cnj (f x))) a"
    by (rule EM_remainder_converges_cnj [rotated]) (auto intro: continuous_intros assms)
  thus "EM_remainder_converges n f a" by simp
qed (intro EM_remainder_converges_cnj assms)

lemma EM_remainder_cnj: 
  assumes "continuous_on {a..} f"
  shows   "EM_remainder n (λx. cnj (f x)) a = cnj (EM_remainder n f a)"
proof (cases "EM_remainder_converges n f a")
  case False
  hence "¬EM_remainder_converges n (λx. cnj (f x)) a"
    by (subst EM_remainder_converges_cnj_iff [OF assms])
  with False show ?thesis by (simp add: EM_remainder_def)
next
  case True
  then obtain L where L: "((λx. EM_remainder' n f (real_of_int a) (real_of_int x))  L) at_top"
    unfolding EM_remainder_converges_def by blast
  note tendsto_cnj [OF this]
  also have "(λx. cnj (EM_remainder' n f (real_of_int a) (real_of_int x))) =
               (λx. EM_remainder' n (λx. cnj (f x)) (real_of_int a) (real_of_int x))"
    by (subst EM_remainder'_bounded_linear [OF bounded_linear_cnj])
       (rule continuous_on_subset [OF assms(1)], auto)
  finally have L': "(  cnj L) at_top" .
  moreover from assms and L have "EM_remainder n f a = L"
    by (intro EM_remainder_eqI)
  ultimately show "EM_remainder n (λx. cnj (f x)) a = cnj (EM_remainder n f a)"
    using L' by (intro EM_remainder_eqI) simp_all
qed

lemma EM_remainder'_combine:
  fixes f :: "real  'a :: banach"
  assumes [continuous_intros]: "continuous_on {a..c} f"
  assumes "a  b" "b  c"
  shows   "EM_remainder' n f a b + EM_remainder' n f b c = EM_remainder' n f a c"
proof -
  have "integral {a..b} (λt. pbernpoly n t *R f t) + integral {b..c} (λt. pbernpoly n t *R f t) =
          integral {a..c} (λt. pbernpoly n t *R f t)"
    by (intro Henstock_Kurzweil_Integration.integral_combine assms integrable_EM_remainder')
  from this [symmetric] show ?thesis by (simp add: EM_remainder'_def algebra_simps) 
qed

lemma uniformly_convergent_EM_remainder':
  fixes f :: "'a  real  'b :: {banach,real_normed_algebra}"
  assumes deriv: "y. a  y  (G has_real_derivative g y) (at y within {a..})"
  assumes integrable: "a' b y. y  A   a  a'  a'  b  
                         (λt. pbernpoly n t *R f y t) integrable_on {a'..b}"
  assumes conv: "convergent (λy. G (real y))"
  assumes bound: "eventually (λx. yA. norm (f y x)  g x) at_top"
  shows   "uniformly_convergent_on A (λb s. EM_remainder' n (f s) a b)"
proof -
  interpret bounded_linear "λx::'b. ((- 1) ^ Suc n / fact n) *R x"
    by (rule bounded_linear_scaleR_right)

  from bounded_pbernpoly obtain C where C: "x. norm (pbernpoly n x)  C" by auto
  from C[of 0] have [simp]: "C  0" by simp

  show ?thesis unfolding EM_remainder'_def
  proof (intro uniformly_convergent_on uniformly_convergent_improper_integral')
    fix x assume "x  a"
    thus "((λx. C * G x) has_real_derivative C * g x) (at x within {a..})"
      by (intro DERIV_cmult deriv)
  next
    fix y a' b assume "y  A" "a  a'" "a'  b"
    thus "(λt. pbernpoly n t *R f y t) integrable_on {a'..b}"
      by (rule integrable)
  next
    from conv obtain L where "(λy. G (real y))  L"
      by (auto simp: convergent_def)
    from tendsto_mult[OF tendsto_const[of C] this]
      show "convergent (λy. C * G (real y))"
        by (auto simp: convergent_def)
  next
    show "F x in at_top. yA. norm (pbernpoly n x *R f y x)  C * g x"
      using C unfolding norm_scaleR 
      by (intro eventually_mono[OF bound] ballI mult_mono) auto
  qed
qed

lemma uniform_limit_EM_remainder:
  fixes f :: "'a  real  'b :: {banach,real_normed_algebra}"
  assumes deriv: "y. a  y  (G has_real_derivative g y) (at y within {a..})"
  assumes integrable: "a' b y. y  A   a  a'  a'  b  
                         (λt. pbernpoly n t *R f y t) integrable_on {a'..b}"
  assumes conv: "convergent (λy. G (real y))"
  assumes bound: "eventually (λx. yA. norm (f y x)  g x) at_top"
  shows   "uniform_limit A (λb s. EM_remainder' n (f s) a b) 
             (λs. EM_remainder n (f s) a) sequentially"
proof -
  have *: "uniformly_convergent_on A (λb s. EM_remainder' n (f s) a b)"
    by (rule uniformly_convergent_EM_remainder'[OF assms])
  also have "?this  ?thesis"
    unfolding uniformly_convergent_uniform_limit_iff
  proof (intro uniform_limit_cong refl always_eventually allI ballI)
    fix s assume "s  A"
    with * have **: "convergent (λb. EM_remainder' n (f s) a b)"
      by (rule uniformly_convergent_imp_convergent)
    show "lim (λb. EM_remainder' n (f s) a b) = EM_remainder n (f s) a"
    proof (rule sym, rule EM_remainder_eqI)
      have "((λx. EM_remainder' n (f s) (real_of_int a) (real x)) 
               lim (λx. EM_remainder' n (f s) (real_of_int a) (real x))) at_top"
        (is "(_  ?L) _") using ** unfolding convergent_LIMSEQ_iff by blast
      hence "((λx. EM_remainder' n (f s) (real_of_int a) (real (nat x)))  ?L) at_top"
        by (rule filterlim_compose) (fact filterlim_nat_sequentially)
      thus "((λx. EM_remainder' n (f s) (real_of_int a) (real_of_int x))  ?L) at_top"
        by (rule Lim_transform_eventually)
           (auto intro: eventually_mono[OF eventually_ge_at_top[of 0]])
    qed
  qed
  finally show  .
qed

lemma tendsto_EM_remainder:
  fixes f :: "real  'b :: {banach,real_normed_algebra}"
  assumes deriv: "y. a  y  (G has_real_derivative g y) (at y within {a..})"
  assumes integrable: "a' b . a  a'  a'  b  
                         (λt. pbernpoly n t *R f t) integrable_on {a'..b}"
  assumes conv: "convergent (λy. G (real y))"
  assumes bound: "eventually (λx. norm (f x)  g x) at_top"
  shows   "filterlim (λb. EM_remainder' n f a b) 
             (nhds (EM_remainder n f a)) sequentially"
proof -
  have "uniform_limit {()} (λb s. EM_remainder' n f a b) 
             (λs. EM_remainder n f a) sequentially"
    using assms by (intro uniform_limit_EM_remainder[where G = G and g = g]) auto
  moreover have "()  {()}" by simp
  ultimately show ?thesis by (rule tendsto_uniform_limitI)
qed

lemma EM_remainder_0 [simp]: "EM_remainder n (λx. 0) a = 0"
  by (rule EM_remainder_eqI) (simp add: EM_remainder'_def)

lemma holomorphic_EM_remainder':
  assumes deriv: 
    "z t. z  U  t  {a..x}  
       ((λz. f z t) has_field_derivative f' z t) (at z within U)"
  assumes int: "b c z e. a  b  c  x  z  U 
                  (λt. of_real (bernpoly n (t - e)) * f z t) integrable_on {b..c}"
  assumes cont: "continuous_on (U × {a..x}) (λ(z, t). f' z t)"
  assumes "convex U"
  shows "(λs. EM_remainder' n (f s) a x) holomorphic_on U"
  unfolding EM_remainder'_def scaleR_conv_of_real
proof (intro holomorphic_intros)
  have holo: "(λz. integral (cbox b c) (λt. of_real (bernpoly n (t - e)) * f z t)) holomorphic_on U"
    if "b  a" "c  x" for b c e :: real
  proof (rule leibniz_rule_holomorphic)
    fix z t assume "z  U" "t  cbox b c"
    thus "((λz. complex_of_real (bernpoly n (t - e)) * f z t) has_field_derivative
             complex_of_real (bernpoly n (t - e)) * f' z t) (at z within U)"
      using that by (intro DERIV_cmult deriv) auto
  next
    fix z assume "z  U"
    thus "(λt. complex_of_real (bernpoly n (t - e)) * f z t) integrable_on cbox b c"
      using that int[of b c z] by auto
  next
    have "continuous_on (U × {b..c}) (λ(z, t). f' z t)"
      using cont by (rule continuous_on_subset) (insert that, auto)
    thus "continuous_on (U × cbox b c) (λ(z, t).
            complex_of_real (bernpoly n (t - e)) * f' z t)"
      by (auto simp: case_prod_unfold intro!: continuous_intros)
  qed fact+

  consider "a > x" | "a  x" "floor x  a" | "a  x" "floor x > a" by force
  hence "(λz. integral (cbox a x) (λt. of_real (pbernpoly n t) * f z t)) holomorphic_on U"
    (is "?f a x holomorphic_on _")
  proof cases
    case 2
    have "(λz. integral (cbox a x) (λt. of_real (bernpoly n (t - of_int x)) * f z t))
                   holomorphic_on U" 
      by (intro holo) auto
    also have "(λz. integral (cbox a x) (λt. of_real (bernpoly n (t - of_int x)) * f z t)) = ?f a x"
    proof (intro ext integral_cong, goal_cases)
      case (1 z t)
      hence "t  a" "t  x" by auto
      hence "floor t = floor x" using 2 by linarith
      thus ?case by (simp add: pbernpoly_def frac_def)
    qed
    finally show ?thesis .
  next
    case 3
          define N :: "int set" where "N = {a..<x}"
    define A where "A = insert {a..of_int a} (insert {of_int x..x} 
                          ((λn. {of_int n..of_int n + 1}) ` N))"
    {
      fix X assume "X  A"
      then consider "X = {a..of_int a}" | "X = {of_int x..x}" |
             n where "X = {of_int n..of_int n + 1}" "n  N" by (auto simp: A_def)
    } note A_cases = this

    have division: "A division_of {a..x}"
    proof (rule division_ofI)
      show "finite A" by (auto simp: A_def N_def)
      fix K assume K: "K  A"
      from 3 have "of_int a  x" 
        using ceiling_le[of a "floor x"] by linarith
      moreover from 3 have "of_int x  a" by linarith
      ultimately show "K  {a..x}" using K 3 by (auto simp: A_def N_def) linarith+
      from K show "K  {}" and "a b. K = cbox a b" by (auto simp: A_def)
    next
      fix K1 K2 assume K: "K1  A" "K2  A" "K1  K2"
      have F1: "interior {a..a}  interior {x..x} = {}" using 3 ceiling_le[of a "floor x"]
        by (auto simp: min_def max_def) 
      hence F2: "interior {x..x}  interior {a..a} = {}" by simp
      have F3: "interior {a..a}  interior {of_int n..of_int n+1} = {}"
               "interior {x..x}  interior {of_int n..of_int n+1} = {}"
               "interior {of_int n..of_int n+1}  interior {a..a} = {}"
               "interior {of_int n..of_int n+1}  interior {x..x} = {}"if "n  N" for n
        using 3 ceiling_le[of a "floor x"] that by (auto simp: min_def max_def N_def)
      have F4: "interior {real_of_int n..of_int n+1}  interior {of_int m..of_int m+1} = {}"
        if "{real_of_int n..of_int n+1}  {of_int m..of_int m+1}" for m n
      proof -
        from that have "n  m" by auto
        thus ?thesis by simp
      qed
      from F1 F2 F3 F4 K show "interior K1  interior K2 = {}"
        by (elim A_cases) (simp_all only: not_False_eq_True)
    next
      show "A = {a..x}"
      proof (cases "a = x")
        case True
        thus ?thesis using 3 by (auto simp: A_def N_def intro: order.trans) linarith+
      next
        case False
        with 3 have *: "a < x" by linarith
        have "A = {a..of_int a}  (nN. {of_int n..of_int (n + 1)})  {of_int x..x}"
          by (simp add: A_def Un_ac)
        also have "(nN. {of_int n..of_int (n + 1)}) = {of_int a..real_of_int x}"
          using * unfolding N_def by (intro Union_atLeastAtMost_real_of_int)
        also have "{a..of_int a}   = {a..real_of_int x}"
          using 3 * by (intro ivl_disj_un) auto
        also have "  {of_int x..x} = {a..x}"
          using 3 * by (intro ivl_disj_un) auto
        finally show ?thesis .
      qed
    qed


    have "(λz. XA. integral X (λt. of_real (bernpoly n (t - Inf X)) * f z t)) 
            holomorphic_on U"
    proof (intro holomorphic_on_sum holo, goal_cases)
      case (1 X)
      from 1 and division have subset: "X  {a..x}" by (auto simp: division_of_def)
      from 1 obtain b c where [simp]: "X = cbox b c" "b  c" by (auto simp: A_def)
      from subset have "b  a" "c  x" by auto
      hence "(λx. integral (cbox b c) (λt. of_real (bernpoly n (t - Inf {b..c})) * f x t)) 
               holomorphic_on U" by (intro holo) auto
      thus ?case by simp
    qed
    also have "?this  (λz. integral {a..x} (λt. of_real (pbernpoly n t) * f z t)) 
                           holomorphic_on U"
    proof (intro holomorphic_cong refl, goal_cases)
      case (1 z)
      have "((λt. of_real (pbernpoly n t) * f z t) has_integral
              (XA. integral X (λt. of_real (bernpoly n (t - Inf X)) * f z t))) {a..x}"
        using division
      proof (rule has_integral_combine_division)
        fix X assume X: "X  A"
        then obtain b c where X': "X = {b..c}" "b  c" by (elim A_cases) auto
        from X and division have "X  {a..x}" by (auto simp: division_of_def)
        with X' have bc: "b  a" "c  x" by auto
        have "((λt. of_real (bernpoly n (t - of_int Inf X)) * f z t) has_integral
                 integral X (λt. of_real (bernpoly n (t - of_int Inf X)) * f z t)) X"
          unfolding X' using z  U bc by (intro integrable_integral int)
        also have "?this  ((λt. of_real (pbernpoly n t) * f z t) has_integral
                 integral X (λt. of_real (bernpoly n (t - of_int Inf X)) * f z t)) X"
        proof (rule has_integral_spike_eq[of "{Sup X}"], goal_cases)
          case (2 t)
          note t = this
          from X  A have "t = Inf X"
          proof (cases rule: A_cases [consumes 1])
            case 1
            with t show ?thesis
              by (intro floor_unique) (auto simp: ceiling_altdef split: if_splits, (linarith+)?)
          next
            case 2
            with t show ?thesis
              by (intro floor_unique) (auto simp: ceiling_altdef split: if_splits, (linarith+)?)
          next
            case 3
            with t show ?thesis
              by (intro floor_unique) (auto simp: ceiling_altdef N_def split: if_splits)
          qed
          thus ?case by (simp add: pbernpoly_def frac_def)
        qed auto
        finally show  .
      qed
      thus ?case by (simp add: has_integral_iff)
    qed
    finally show ?thesis by simp
  qed auto
  thus "(λz. integral {a..x} (λt. of_real (pbernpoly n t) * f z t)) holomorphic_on U"
    by simp
qed

lemma
  assumes deriv: "y. a  y  (G has_real_derivative g y) (at y within {a..})"
  assumes deriv': 
    "z t x. z  U  x  a  t  {a..x}  
       ((λz. f z t) has_field_derivative f' z t) (at z within U)"
  assumes cont: "continuous_on (U × {of_int a..}) (λ(z, t). f' z t)"
  assumes int: "b c z e. a  b  z  U 
                  (λt. of_real (bernpoly n (t - e)) * f z t) integrable_on {b..c}"
  assumes int': "a' b y. y  U   a  a'  a'  b  
                         (λt. pbernpoly n t *R f y t) integrable_on {a'..b}"
  assumes conv: "convergent (λy. G (real y))"
  assumes bound: "eventually (λx. yU. norm (f y x)  g x) at_top"
  assumes "open U" 
  shows analytic_EM_remainder: "(λs::complex. EM_remainder n (f s) a) analytic_on U"
    and holomorphic_EM_remainder: "(λs::complex. EM_remainder n (f s) a) holomorphic_on U"
proof -
  show "(λs::complex. EM_remainder n (f s) a) analytic_on U"
  unfolding analytic_on_def
  proof
    fix z assume "z  U"
    from z  U and open U obtain ε where ε: "ε > 0" "ball z ε  U"
      by (auto simp: open_contains_ball)
    have "(λs. EM_remainder n (f s) a) holomorphic_on ball z ε"
    proof (rule holomorphic_uniform_sequence)
      fix x :: nat
      show "(λs. EM_remainder' n (f s) a x) holomorphic_on ball z ε"
      proof (rule holomorphic_EM_remainder', goal_cases)
        fix s t assume "s  ball z ε" "t  {real_of_int a..real x}"
        thus "((λz. f z t) has_field_derivative f' s t) (at s within ball z ε)"
          using ε by (intro DERIV_subset[OF deriv'[of _ x]]) auto
      next
        case (2 b c s e)
        with ε have "s  U" by blast
        with 2 show ?case using ε int[of b s e c] by (cases "a  x") auto
      next
        from cont show "continuous_on (ball z ε × {real_of_int a..real x}) (λ(z, t). f' z t)"
          by (rule continuous_on_subset) (insert ε, auto)
      qed (auto)
    next
      fix s assume s: "s  ball z ε"
      have "open (ball z ε)" by simp
      with s obtain δ where δ: "δ > 0" "cball s δ  ball z ε" 
        unfolding open_contains_cball by blast
      moreover have bound': "eventually (λx. ycball s δ. norm (f y x)  g x) at_top"
        by (intro eventually_mono [OF bound]) (insert δ ε, auto)
      have "uniform_limit (cball s δ) (λx s. EM_remainder' n (f s) (real_of_int a) (real x))
                        (λs. EM_remainder n (f s) a) sequentially"
        by (rule uniform_limit_EM_remainder[OF deriv int' conv bound']) (insert δ ε s, auto)
      ultimately show "δ>0. cball s δ  ball z ε  uniform_limit (cball s δ)
                          (λx s. EM_remainder' n (f s) (real_of_int a) (real x))
                          (λs. EM_remainder n (f s) a) sequentially" by blast
    qed auto
    with ε show "ε>0. (λs. EM_remainder n (f s) a) holomorphic_on ball z ε"
      by blast
  qed
  thus "(λs::complex. EM_remainder n (f s) a) holomorphic_on U"
    by (rule analytic_imp_holomorphic)
qed



text ‹
  The following lemma is the first step in the proof of the Euler--MacLaurin formula: 
  We show the relationship between the first-order remainder term and the difference of 
  the integral and the sum.
›
  
context
  fixes f f' :: "real  'a :: banach"
  fixes a b :: int and I S :: 'a
  fixes Y :: "real set"
  assumes "a  b"
  assumes fin: "finite Y"
  assumes cont: "continuous_on {real_of_int a..real_of_int b} f"
  assumes deriv [derivative_intros]: 
            "x::real. x  {a..b} - Y  (f has_vector_derivative f' x) (at x)"
  defines S_def: "S  (i{a<..b}. f i)" and I_def: "I  integral {a..b} f"
begin
  
lemma
  diff_sum_integral_has_integral_int:
    "((λt. (frac t - 1/2) *R f' t) has_integral (S - I - (f b - f a) /R 2)) {a..b}"
proof (cases "a = b")
  case True
  thus ?thesis by (simp_all add: S_def I_def has_integral_refl)
next
  case False
  with a  b have ab: "a < b" by simp
  let ?A = "(λn. {real_of_int n..real_of_int (n+1)}) ` {a..<b}"
  have division: "?A division_of {of_int a..of_int b}"
    using Union_atLeastAtMost_real_of_int[OF ab] by (simp add: division_of_def)
  have cont' [continuous_intros]: "continuous_on A f" if "A  {of_int a..of_int b}" for A
    using continuous_on_subset[OF cont that] .
  
  define d where "d = (λx. (f x + f (x + 1)) /R 2 - integral {x..x+1} f)"    
  have "((λt. (frac t - 1/2) *R f' t) has_integral d i) {of_int i..of_int (i+1)}"
    if i: "i  {a..<b}" for i
  proof (rule has_integral_spike)
    show "(frac x - 1 / 2) *R f' x = (x - of_int i - 1 / 2) *R f' x"
         if "x  {of_int i..of_int (i + 1)} - {of_int (i + 1)}" for x
    proof -
      have "x  of_int i" "x < of_int (i + 1)" using that by auto
      hence "floor x = of_int i" by (subst floor_unique) auto
      thus ?thesis by (simp add: frac_def)
    qed
  next
    define h where "h = (λx::real. (x - of_int i - 1 / 2) *R f' x)"
    define g where "g = (λx::real. (x - of_int i - 1/2) *R f x - integral {of_int i..x} f)"
    have *: "((λx. integral {real_of_int i..x} f) has_vector_derivative f x) (at x within {i..i+1})"
      if "x  {of_int i<..<of_int i + 1}" for x using that i
      by (intro integral_has_vector_derivative cont') auto
    have "((λx. integral {real_of_int i..x} f) has_vector_derivative f x) (at x)" 
      if "x  {of_int i<..<of_int i + 1}" for x 
      using that i at_within_interior[of x "{of_int i..of_int (i + 1)}"] *[of x] by simp
    hence "(h has_integral g (of_int (i + 1)) - g (of_int i)) {of_int i..of_int (i+1)}"
      unfolding g_def h_def using that 
      by (intro fundamental_theorem_of_calculus_interior_strong[OF fin])
         (auto intro!: derivative_eq_intros continuous_intros indefinite_integral_continuous_1
             integrable_continuous_real)
    also have "g (of_int (i + 1)) - g (of_int i) = d i"
      by (simp add: g_def scaleR_add_right [symmetric] d_def)
    finally show "(h has_integral d i) {of_int i..of_int (i + 1)}" .
  qed simp_all
  hence *: "I. I?A  ((λx. (frac x - 1 / 2) *R f' x) has_integral d (Inf I)) I"
    by (auto simp: add_ac)
  have "((λx::real. (frac x - 1 / 2) *R f' x) has_integral (I?A. d (Inf I))) (?A)"
    by (intro has_integral_Union * finite_imageI) (force intro!: negligible_atLeastAtMostI pairwiseI)+
  also have "?A = {of_int a..of_int b}"
    by (intro Union_atLeastAtMost_real_of_int ab)
  also have "(I?A. d (Inf I)) = (i=a..<b. d i)"
    by (subst sum.reindex) (auto simp: inj_on_def)
  also have " = (1 / 2) *R ((i = a..<b. f (real_of_int i)) +
                    (i = a..<b. f (real_of_int (i + 1)))) -
                  (i = a..<b. integral {real_of_int i..1 + real_of_int i} f)"
    (is "_ = _ *R (?S1 + ?S2) - ?S3") 
    by (simp add: d_def algebra_simps sum.distrib sum_subtractf scaleR_sum_right)
  also have "?S1 = (i = a..b. f (real_of_int i)) - f b"
    unfolding S_def using ab by (subst sum_atLeastAtMost_int_last) auto
  also have "(i = a..b. f (real_of_int i)) = S + f a"
    unfolding S_def using ab by (subst sum_atLeastAtMost_int_head) auto
  also have "?S2 = S" unfolding S_def
    by (intro sum.reindex_bij_witness[of _ "λi. i-1" "λi. i+1"]) auto
  also have "(1 / 2) *R (S + f a - f b + S) = 
                (1/2) *R S + (1/2) *R S - (f b - f a) /R 2"
    by (simp add: algebra_simps)
  also have "(1/2) *R S + (1/2) *R S = S" by (simp add: scaleR_add_right [symmetric])
      
  also have "?S3 = (I?A. integral I f)"
    by (subst sum.reindex) (auto simp: inj_on_def add_ac)
  also have " = I" unfolding I_def
    by (intro integral_combine_division_topdown [symmetric] division integrable_continuous_real 
              continuous_intros) simp_all
  finally show ?thesis by (simp add: algebra_simps)
qed

lemma diff_sum_integral_has_integral_int':
  "((λt. pbernpoly 1 t *R f' t) has_integral (S - I - (f b - f a) /R 2 )) {a..b}"
  using diff_sum_integral_has_integral_int by (simp add: pbernpoly_def bernpoly_def)
  
lemma EM_remainder'_Suc_0: "EM_remainder' (Suc 0) f' a b = S - I - (f b - f a) /R 2"
  using diff_sum_integral_has_integral_int' by (simp add: has_integral_iff EM_remainder'_def)

end


text ‹
  Next, we show that the $n$-th-order remainder can be expressed in terms of the $n+1$-th-order 
  remainder term. Iterating this essentially yields the Euler--MacLaurin formula.
›

context
  fixes f f' :: "real  'a :: banach" and a b :: int and n :: nat and A :: "real set"
  assumes ab: "a  b" and n: "n > 0"
  assumes fin:   "finite A"
  assumes cont:  "continuous_on {of_int a..of_int b} f"
  assumes cont': "continuous_on {of_int a..of_int b} f'"
  assumes deriv: "x. x  {of_int a<..<of_int b} - A  (f has_vector_derivative f' x) (at x)"
begin
  
lemma EM_remainder'_integral_conv_Suc:
  shows   "integral {a..b} (λt. pbernpoly n t *R f t) =
              (bernoulli (Suc n) / real (Suc n)) *R (f b - f a) -
              integral {a..b} (λt. pbernpoly (Suc n) t *R f' t) /R real (Suc n)"
  unfolding EM_remainder'_def
proof -
  let ?h = "λi. (pbernpoly (Suc n) (real_of_int i) / real (Suc n)) *R f (real_of_int i)"
  define T where "T = integral {a..b} (λt. (pbernpoly (Suc n) t / real (Suc n)) *R f' t)"
  note [derivative_intros] = has_field_derivative_pbernpoly_Suc'
  let ?A = "real_of_int ` {a..b}  A"
    
  have "((λt. pbernpoly n t *R f t) has_integral (-T + (?h b - ?h a))) {a..b}"
  proof (rule integration_by_parts_interior_strong[OF bounded_bilinear_scaleR])
    from fin show "finite ?A" by simp
    from n > 0 show "continuous_on {of_int a..of_int b} (λt. pbernpoly (Suc n) t / real (Suc n))"
      by (intro continuous_intros) auto
    show "continuous_on {of_int a..of_int b} f" by fact
    show "(f has_vector_derivative f' t) (at t)" if "t  {of_int a<..<of_int b} - ?A" for t
      using deriv[of t] that by auto
    have "(λt. pbernpoly (Suc n) t *R f' t) integrable_on {a..b}"
      by (intro integrable_EM_remainder' cont')
    hence "(λt. (1 / real (Suc n)) *R pbernpoly (Suc n) t *R f' t) integrable_on {a..b}" 
      by (rule integrable_cmul)
    also have "(λt. (1 / real (Suc n)) *R pbernpoly (Suc n) t *R f' t) =
                 (λt. (pbernpoly (Suc n) t / real (Suc n)) *R f' t)"
      by (rule ext) (simp add: algebra_simps)
    finally show "((λt. (pbernpoly (Suc n) t / real (Suc n)) *R f' t) 
                    has_integral ?h b - ?h a - (- T + (?h b - ?h a))) {a..b}"
      using integrable_EM_remainder'[of a b f' "Suc n"]
        by (simp add: has_integral_iff T_def)
    qed (insert ab n, auto intro!: derivative_eq_intros
           simp: has_real_derivative_iff_has_vector_derivative [symmetric] not_le elim!: Ints_cases)
    also have "?h b - ?h a = (bernoulli (Suc n) / real (Suc n)) *R (f b - f a)"
      using n by (simp add: algebra_simps bernoulli'_def)
    finally have "integral {a..b} (λt. pbernpoly n t *R f t) =  - T"
      by (simp add: has_integral_iff)
    also have "T = integral {a..b} (λt. (1 / real (Suc n)) *R (pbernpoly (Suc n) t) *R f' t)"
      by (simp add: T_def)
    also have " = integral {a..b} (λt. pbernpoly (Suc n) t *R f' t) /R real (Suc n)"
      by (subst integral_cmul) (simp_all add: divide_simps)
    finally show ?thesis .
qed 

lemma EM_remainder'_conv_Suc: 
  "EM_remainder' n f a b = 
     ((-1) ^ Suc n * bernoulli (Suc n) / fact (Suc n)) *R (f b - f a) + 
     EM_remainder' (Suc n) f' a b"
  by (simp add: EM_remainder'_def EM_remainder'_integral_conv_Suc scaleR_diff_right 
                scaleR_add_right field_simps del: of_nat_Suc)

end

context
  fixes f f' :: "real  'a :: banach" and a :: int and n :: nat and A :: "real set" and C
  assumes n: "n > 0"
  assumes fin:   "finite A"
  assumes cont:  "continuous_on {of_int a..} f"
  assumes cont': "continuous_on {of_int a..} f'"
  assumes lim:   "(f  C) at_top"
  assumes deriv: "x. x  {of_int a<..} - A  (f has_vector_derivative f' x) (at x)"
begin  
  
lemma
  shows EM_remainder_converges_iff_Suc_converges:
          "EM_remainder_converges n f a  EM_remainder_converges (Suc n) f' a"
  and   EM_remainder_conv_Suc: 
           "EM_remainder_converges n f a  
              EM_remainder n f a = 
                  ((-1) ^ Suc n * bernoulli (Suc n) / fact (Suc n)) *R (C - f a) + 
                  EM_remainder (Suc n) f' a"
proof (rule iffI)
  define g where "g = (λx. ((-1) ^ Suc n * bernoulli (Suc n) / fact (Suc n)) *R (f x - f a))"
  define G where "G = ((-1) ^ Suc n * bernoulli (Suc n) / fact (Suc n)) *R (C - f a)"
  have limit_g: "(g  G) at_top" unfolding g_def G_def by (intro tendsto_intros lim)
  have *: "eventually (λx. EM_remainder' n f (real_of_int a) (real_of_int x) = 
              g x + EM_remainder' (Suc n) f' (real_of_int a) (real_of_int x)) at_top"
    using eventually_ge_at_top[of a]
  proof eventually_elim
    case (elim b)
    thus ?case
    using EM_remainder'_conv_Suc[OF elim n fin continuous_on_subset[OF cont] 
            continuous_on_subset[OF cont'] deriv] by (auto simp: g_def)
  qed
  
  {
    assume "EM_remainder_converges n f a"
    then obtain L
      where L: "((λb. EM_remainder' n f (real_of_int a) (real_of_int b))  L) at_top"
      by (auto simp: EM_remainder_converges_def)
    have *: "((λb. EM_remainder' (Suc n) f' (real_of_int a) (real_of_int b))  L - G) at_top"
    proof (rule Lim_transform_eventually)
      show "F x in at_top. EM_remainder' n f (real_of_int a) (real_of_int x) - g x =
               EM_remainder' (Suc n) f' (real_of_int a) (real_of_int x)"
        using * by (simp add: algebra_simps)
      show "((λx. EM_remainder' n f (real_of_int a) (real_of_int x) - g x)  L - G) at_top"
        by (intro tendsto_intros filterlim_compose[OF limit_g] L)
    qed
    from * show "EM_remainder_converges (Suc n) f' a" unfolding EM_remainder_converges_def ..
    from * have "EM_remainder (Suc n) f' a = L - G" by (rule EM_remainder_eqI)
    moreover from L have "EM_remainder n f a = L" by (rule EM_remainder_eqI)
    ultimately show "EM_remainder n f a = G + EM_remainder (Suc n) f' a" by (simp add: G_def)
  }
  {
    assume "EM_remainder_converges (Suc n) f' a" 
    then obtain L
      where L: "((λb. EM_remainder' (Suc n) f' (real_of_int a) (real_of_int b))  L) at_top"
      by (auto simp: EM_remainder_converges_def)
    have *: "((λb. EM_remainder' n f (real_of_int a) (real_of_int b))  G + L) at_top"
    proof (rule Lim_transform_eventually)
      show "F x in at_top. g x + EM_remainder' (Suc n) f' (real_of_int a) (real_of_int x) =
                EM_remainder' n f (real_of_int a) (real_of_int x)"
        using * by (subst eq_commute)
      show "((λx. g x + EM_remainder' (Suc n) f' (real_of_int a) (real_of_int x))  G + L) at_top"
        by (intro tendsto_intros filterlim_compose[OF limit_g] L)
    qed
    thus "EM_remainder_converges n f a" unfolding EM_remainder_converges_def ..
  }
qed

end
  


subsection ‹The conventional version of the Euler--MacLaurin formula›

text ‹
  The following theorems are the classic Euler--MacLaurin formula that can be found,
  with slight variations, in many sources (e.\,g.\ cite"AS_HMF" and "apostol99" and "GKP_CM").
›

context
  fixes f :: "real  'a :: banach"
  fixes fs :: "nat  real  'a"
  fixes a b :: int assumes ab: "a  b"
  fixes N :: nat assumes N: "N > 0"
  fixes Y :: "real set" assumes fin: "finite Y"
  assumes fs_0 [simp]: "fs 0 = f"
  assumes fs_cont [continuous_intros]:  
    "k. k  N  continuous_on {real_of_int a..real_of_int b} (fs k)"
  assumes fs_deriv [derivative_intros]: 
    "k x. k < N  x  {a..b} - Y  (fs k has_vector_derivative fs (Suc k) x) (at x)"
begin

theorem euler_maclaurin_raw_strong_int:
  defines "S  (i{a<..b}. f (of_int i))"
  defines "I  integral {of_int a..of_int b} f"
  defines "c'  λk. (bernoulli' (Suc k) / fact (Suc k)) *R (fs k b - fs k a)"
  shows   "S - I = (k<N. c' k) + EM_remainder' N (fs N) a b"
proof -
  define c :: "nat  'a" 
    where "c = (λk. ((-1) ^ (Suc k) * bernoulli (Suc k) / fact (Suc k)) *R (fs k b - fs k a))"
  have "S - I = (k<m. c k) + EM_remainder' m (fs m) a b" if "m  1" "m  N" for m
  using that
  proof (induction m rule: dec_induct)
    case base
    with ab fin fs_cont[of 0] show ?case using fs_deriv[of 0] N unfolding One_nat_def
      by (subst EM_remainder'_Suc_0[of _ _ Y f]) (simp_all add: algebra_simps S_def I_def c_def)
  next
    case (step n)
    from step.prems have "S - I = (k<n. c k) + EM_remainder' n (fs n) a b"
      by (intro step.IH) simp_all
    also have "(k<n. c k) = (k<Suc n. c k) +
                  (((-1) ^ n * bernoulli (Suc n) / fact (Suc n)) *R (fs n b - fs n a))"
      (is "_ = _ + ?c") by (simp add: EM_remainder'_Suc_0 c_def)
    also have " + EM_remainder' n (fs n) a b = (k<Suc n. c k) + (?c + EM_remainder' n (fs n) a b)"
      by (simp add: add.assoc)
    also from step.prems step.hyps ab fin
      have "?c + EM_remainder' n (fs n) a b = EM_remainder' (Suc n) (fs (Suc n)) a b"
      by (subst EM_remainder'_conv_Suc [where A = Y]) 
         (auto intro!: fs_deriv fs_cont)
    finally show ?case .
  qed
  from this[of N] and N 
    have "S - I = sum c {..<N} + EM_remainder' N (fs N) (real_of_int a) (real_of_int b)" by simp
  also have "sum c {..<N} = sum c' {..<N}"
  proof (intro sum.cong refl)
    fix k :: nat
    show "c k = c' k"
      by (cases "even k")
         (auto simp: c_def c'_def bernoulli'_def algebra_simps bernoulli_odd_eq_0)
  qed
  finally show ?thesis .
qed

end

theorem euler_maclaurin_strong_raw_nat:
  assumes "a  b" "0 < N" "finite Y" "fs 0 = f"
    "(k. k  N  continuous_on {real a..real b} (fs k))"
    "(k x. k < N  x  {real a..real b} - Y 
       (fs k has_vector_derivative fs (Suc k) x) (at x))"
  shows "(i{a<..b}. f (real i)) - integral {real a..real b} f =
           (k<N. (bernoulli' (Suc k) / fact (Suc k)) *R (fs k (real b) - fs k (real a))) +
           EM_remainder' N (fs N) (real a) (real b)"
proof -
  have "(i{int a<..int b}. f (real_of_int i)) - 
           integral {real_of_int (int a)..real_of_int (int b)} f =
           (k<N. (bernoulli' (Suc k) / fact (Suc k)) *R 
             (fs k (real_of_int (int b)) - fs k (real_of_int (int a)))) +
           EM_remainder' N (fs N) (real_of_int (int a)) (real_of_int (int b))"
    using assms by (intro euler_maclaurin_raw_strong_int[where Y = Y] assms) simp_all
  also have "(i{int a<..int b}. f (real_of_int i)) = (i{a<..b}. f (real i))"
    by (intro sum.reindex_bij_witness[of _ int nat]) auto
  finally show ?thesis by simp
qed

  
subsection ‹The ``Concrete Mathematics'' version of the Euler--MacLaurin formula›

text ‹
  As explained in \textit{Concrete Mathematics}~cite"GKP_CM", the above form of the 
  formula has some drawbacks: When applying it to determine the asymptotics of some concrete 
  function, one is usually left with several different unwieldy constant terms that are difficult 
  to get rid of.

  There is no general way to determine what these constant terms are, but in concrete applications, 
  they can often be determined or estimated by other means. We can therefore simply group all the 
  constant terms into a single constant and have the user provide a proof of what it is.
›

locale euler_maclaurin_int =
  fixes F f :: "real  'a :: banach"
  fixes fs :: "nat  real  'a"
  fixes a :: int
  fixes N :: nat assumes N: "N > 0"
  fixes C :: 'a
  fixes Y :: "real set" assumes fin: "finite Y"
  assumes fs_0 [simp]: "fs 0 = f"
  assumes fs_cont [continuous_intros]:  
    "k. k  N  continuous_on {real_of_int a..} (fs k)"
  assumes fs_deriv [derivative_intros]: 
    "k x. k < N  x  {of_int a..} - Y  (fs k has_vector_derivative fs (Suc k) x) (at x)"
  assumes F_cont [continuous_intros]: "continuous_on {of_int a..} F"
  assumes F_deriv [derivative_intros]: 
    "x. x  {of_int a..} - Y  (F has_vector_derivative f x) (at x)" 
  assumes limit: 
    "((λb. (k=a..b. f k) - F (of_int b) - 
         (i<N. (bernoulli' (Suc i) / fact (Suc i)) *R fs i (of_int b)))  C) at_top"
begin

context
  fixes C' T
  defines "C'  -f a + F a + C + (k<N. (bernoulli' (Suc k) / fact (Suc k)) *R (fs k (of_int a)))"
      and "T  (λx. i<N. (bernoulli' (Suc i) / fact (Suc i)) *R fs i x)"
begin
  
lemma euler_maclaurin_strong_int_aux:
  assumes ab: "a  b"
  defines "S  (k=a..b. f (of_int k))"
  shows "S - F (of_int b) - T (of_int b) = EM_remainder' N (fs N) (of_int a) (of_int b) + (C - C')"
proof (cases "a = b")
  case True
  thus ?thesis unfolding C'_def by (simp add: S_def EM_remainder'_def T_def)
next
  case False
  with assms have ab: "a < b" by simp
  define T' where "T' = (k<N. (bernoulli' (Suc k) / fact (Suc k)) *R (fs k (of_int a)))"
  have "(i{a<..b}. f (of_int i)) - integral {of_int a..of_int b} f =
          (k<N. (bernoulli' (Suc k) / fact (Suc k)) *R (fs k (of_int b) - fs k (of_int a))) +
          EM_remainder' N (fs N) (of_int a) (of_int b)" using ab
    by (intro euler_maclaurin_raw_strong_int [where Y = Y] N fin fs_0
              continuous_on_subset[OF fs_cont] fs_deriv) auto
  also have "(f has_integral (F b - F a)) {of_int a..of_int b}" using ab
    by (intro fundamental_theorem_of_calculus_strong[OF fin])
       (auto intro!: continuous_on_subset[OF F_cont] derivative_intros)
  hence "integral {of_int a..of_int b} f = F (of_int b) - F (of_int a)"
    by (simp add: has_integral_iff)
  also have "(k<N. (bernoulli' (Suc k) / fact (Suc k)) *R (fs k (of_int b) - fs k (of_int a))) =
               T (of_int b) - T'"
    by (simp add: T_def T'_def algebra_simps sum_subtractf)
  also have "(i{a<..b}. f (of_int i)) = S - f (of_int a)"
    unfolding S_def using ab by (subst sum_atLeastAtMost_int_head) auto
  finally show ?thesis by (simp add: algebra_simps C'_def T'_def)  
qed

lemma EM_remainder_limit:
  assumes ab: "a  b"
  defines "D  EM_remainder' N (fs N) (of_int a) (of_int b)"
  shows "EM_remainder N (fs N) b = C' - D"
    and EM_remainder_converges: "EM_remainder_converges N (fs N) b"
proof -
  note limit
  also have "((λb. (k = a..b. f (of_int k)) - F (of_int b) -
               (i<N. (bernoulli' (Suc i) / fact (Suc i)) *R fs i (of_int b)))  C) at_top =
             ((λb. (k = a..b. f (of_int k)) - F (of_int b) - T (of_int b))  C) at_top"
    unfolding T_def ..
  also have "eventually (λx. (k=a..x. f k) - F (of_int x) - T (of_int x) = 
               EM_remainder' N (fs N) (of_int a) (of_int x) + (C - C')) at_top"
    (is "eventually (λx. ?f x = ?g x) _")
    using eventually_gt_at_top[of b]
    by eventually_elim (rule euler_maclaurin_strong_int_aux, insert ab, simp_all)
  hence "(?f  C) at_top  (?g  C) at_top" by (intro filterlim_cong refl)
  finally have "((λx. ?g x - (C - C'))  (C - (C - C'))) at_top"
    by (rule tendsto_diff[OF _ tendsto_const])
  hence *: "((λx. EM_remainder' N (fs N) (of_int a) (of_int x))  C') at_top"
    by simp