Theory Bernoulli.Bernoulli

(*  
  File:       Bernoulli.thy
  Author:     Lukas Bulwahn <lukas.bulwahn-at-gmail.com> 
  Author:     Manuel Eberl <manuel@pruvisto.org> 
*)
section ‹Bernoulli numbers›

theory Bernoulli
imports Complex_Main
begin

subsection ‹Preliminaries›
  
lemma power_numeral_reduce: "a ^ numeral n = a * a ^ pred_numeral n"
  by (simp only: numeral_eq_Suc power_Suc)

lemma fact_diff_Suc: "n < Suc m ⟹ fact (Suc m - n) = of_nat (Suc m - n) * fact (m - n)"
  by (subst fact_reduce) auto

lemma of_nat_binomial_Suc:
  assumes "k ≤ n"
  shows   "(of_nat (Suc n choose k) :: 'a :: field_char_0) = 
             of_nat (Suc n) / of_nat (Suc n - k) * of_nat (n choose k)"
  using assms by (simp add: binomial_fact divide_simps fact_diff_Suc of_nat_diff del: of_nat_Suc)

lemma integrals_eq:
  assumes "f 0 = g 0"
  assumes "⋀ x. ((λx. f x - g x) has_real_derivative 0) (at x)"
  shows "f x = g x"
proof -
  show "f x = g x"
  proof (cases "x ≠ 0")
    case True
    from assms DERIV_const_ratio_const[OF this, of "λx. f x - g x" 0]
    show ?thesis by auto
  qed (simp add: assms)
qed

lemma sum_diff: "((∑i≤n::nat. f (i + 1) - f i)::'a::field) = f (n + 1) - f 0"
  by (induct n) (auto simp add: field_simps)
    
lemma Rats_sum: "(⋀x. x ∈ A ⟹ f x ∈ ℚ) ⟹ sum f A ∈ ℚ"
  by (induction A rule: infinite_finite_induct) simp_all


subsection ‹Bernoulli Numbers and Bernoulli Polynomials›

declare sum.cong [fundef_cong]

fun bernoulli :: "nat ⇒ 'a :: field_char_0"
where
  "bernoulli 0 = 1"
| "bernoulli (Suc n) =  (-1 / (of_nat n + 2)) * (∑k ≤ n. (of_nat (n + 2 choose k) * bernoulli k))"
  
declare bernoulli.simps[simp del]
  
lemmas bernoulli_0 [simp] = bernoulli.simps(1)
lemmas bernoulli_Suc = bernoulli.simps(2)
lemma bernoulli_1 [simp]: "bernoulli 1 = -1/2" by (simp add: bernoulli_Suc)
lemma bernoulli_Suc_0 [simp]: "bernoulli (Suc 0) = -1/2" by (simp add: bernoulli_Suc)

lemma of_rat_bernoulli: "of_rat (bernoulli n) = bernoulli n"
  by (induction n rule: bernoulli.induct) 
     (auto simp: bernoulli_Suc of_rat_minus of_rat_sum of_rat_divide of_rat_mult of_rat_add)

(* TODO: Move *)
lemma of_real_of_rat: "of_real (of_rat x) = (of_rat x :: 'a :: real_field)"
  by (cases x) (auto simp: of_rat_rat)

lemma of_real_bernoulli: "of_real (bernoulli n) = (bernoulli n :: 'a :: real_field)"
  by (metis of_rat_bernoulli of_real_of_rat)

    
text ‹
  The ``normal'' Bernoulli numbers are the negative Bernoulli numbers $B_n^{-}$ we just defined
  (so called because $B_1^{-} = -\frac{1}{2}$). There is also another convention, the 
  positive Bernoulli numbers $B_n^{+}$, which differ from the negative ones only in that 
  $B_1^{+} = \frac{1}{2}$. Both conventions have their justification, since a number of theorems 
  are easier to state with one than the other.
›
definition bernoulli' :: "nat ⇒ 'a :: field_char_0" where
  "bernoulli' n = (if n = 1 then 1/2 else bernoulli n)"
  
lemma bernoulli'_0 [simp]: "bernoulli' 0 = 1" by (simp add: bernoulli'_def)
    
lemma bernoulli'_1 [simp]: "bernoulli' (Suc 0) = 1/2"
  by (simp add: bernoulli'_def)

lemma bernoulli_conv_bernoulli': "n ≠ 1 ⟹ bernoulli n = bernoulli' n"
  by (simp add: bernoulli'_def)
    
lemma bernoulli'_conv_bernoulli: "n ≠ 1 ⟹ bernoulli' n = bernoulli n"
  by (simp add: bernoulli'_def)
    
lemma bernoulli_conv_bernoulli'_if: 
    "n ≠ 1 ⟹ bernoulli n = (if n = 1 then -1/2 else bernoulli' n)"
  by (simp add: bernoulli'_def)

lemma of_rat_bernoulli': "of_rat (bernoulli' n) = bernoulli' n"
  by (auto simp: bernoulli'_def of_rat_bernoulli of_rat_divide)

lemma of_real_bernoulli': "of_real (bernoulli' n) = (bernoulli' n :: 'a :: real_field)"
  by (metis of_rat_bernoulli' of_real_of_rat)

lemma bernoulli_in_Rats: "bernoulli n ∈ ℚ"
proof (induction n rule: less_induct)
  case (less n)
  thus ?case
    by (cases n) (auto simp: bernoulli_Suc intro!: Rats_sum Rats_divide)
qed

lemma bernoulli'_in_Rats: "bernoulli' n ∈ ℚ"
  by (simp add: bernoulli'_def bernoulli_in_Rats)

definition bernpoly :: "nat ⇒ 'a ⇒ 'a :: real_algebra_1" where
  "bernpoly n = (λx. ∑k ≤ n. of_nat (n choose k) * of_real (bernoulli k) * x ^ (n - k))"
  
lemma bernpoly_altdef:
  "bernpoly n = (λx. ∑k≤n. of_nat (n choose k) * of_real (bernoulli (n - k)) * x ^ k)"
proof
  fix x :: 'a
  have "bernpoly n x = (∑k≤n. of_nat (n choose (n - k)) * of_real (bernoulli (n - k)) * x ^ (n - (n - k)))"
    unfolding bernpoly_def by (rule sum.reindex_bij_witness[of _ "λk. n - k" "λk. n - k"]) simp_all
  also have "… = (∑k≤n. of_nat (n choose k) * of_real (bernoulli (n - k)) * x ^ k)"
    by (intro sum.cong refl) (simp_all add: binomial_symmetric [symmetric])
  finally show "bernpoly n x = …" .
qed

lemma bernoulli_Suc': 
  "bernoulli (Suc n) = -1/(of_nat n + 2) * (∑k≤n. of_nat (n + 2 choose (k + 2)) * bernoulli (n - k))" 
proof -
  have "bernoulli (Suc n) = - 1 / (of_nat n + 2) * (∑k≤n. of_nat (n + 2 choose k) * bernoulli k)"
    unfolding bernoulli.simps ..
  also have "(∑k≤n. of_nat (n + 2 choose k) * bernoulli k) = 
               (∑k≤n. of_nat (n + 2 choose (n - k)) * bernoulli (n - k))"
    by (rule sum.reindex_bij_witness[of _ "λk. n - k" "λk. n - k"]) simp_all
  also have "… = (∑k≤n. of_nat (n + 2 choose (k + 2)) * bernoulli (n - k))"
    by (intro sum.cong refl, subst binomial_symmetric) simp_all
  finally show ?thesis .
qed
  

subsection ‹Basic Observations on Bernoulli Polynomials›

lemma bernpoly_0 [simp]: "bernpoly n 0 = of_real (bernoulli n)"
proof -
  have "bernpoly n 0 = (∑k≤n. of_nat (n choose k) * of_real (bernoulli (n - k)) * 0 ^ k :: 'a)"
    by (simp add: bernpoly_altdef)
  also have "… = (∑k∈{0}. of_nat (n choose k) * of_real (bernoulli (n - k)) * 0 ^ k)"
    by (rule sum.mono_neutral_right) (auto simp: power_0_left)
  also have "… = of_real (bernoulli n)"
    by simp
  finally show ?thesis .
qed

lemma continuous_on_bernpoly [continuous_intros]: 
  "continuous_on A (bernpoly n :: 'a ⇒ 'a :: real_normed_algebra_1)"
  unfolding bernpoly_def by (auto intro!: continuous_intros)

lemma isCont_bernpoly [continuous_intros]: 
  "isCont (bernpoly n :: 'a ⇒ 'a :: real_normed_algebra_1) x"
  unfolding bernpoly_def by (auto intro!: continuous_intros)  

lemma has_field_derivative_bernpoly:
  "(bernpoly (Suc n) has_field_derivative 
     (of_nat (n + 1) * bernpoly n x :: 'a :: real_normed_field)) (at x)"
proof -
  have "(bernpoly (Suc n) has_field_derivative 
          (∑k≤n. of_nat (Suc n - k) * x ^ (n - k) * (of_nat (Suc n choose k) * 
            of_real (bernoulli k)))) (at x)" (is "(_ has_field_derivative ?D) _")
    unfolding bernpoly_def by (rule DERIV_cong) (fast intro!: derivative_intros, simp)
  also have "?D = of_nat (n + 1) * bernpoly n x" 
    unfolding bernpoly_def sum_distrib_left
    by (force simp: of_nat_binomial_Suc nat_le_iff_add intro: sum.cong)
  ultimately show ?thesis by (auto simp del: of_nat_Suc One_nat_def)
qed
  
lemmas has_field_derivative_bernpoly' [derivative_intros] =
  DERIV_chain'[OF _ has_field_derivative_bernpoly]    

lemma sum_binomial_times_bernoulli:
  "(∑k≤n. of_nat ((Suc n) choose k) * bernoulli k) = (if n = 0 then 1 else 0)"
proof (cases n)
  case (Suc m)
  have "(∑k≤n. of_nat (Suc n choose k) * bernoulli k :: 'a) =
          (∑k≤m. of_nat (m + 2 choose k) * bernoulli k) * 
          (1 - of_nat (m + 2 choose Suc m) / of_nat (m + 2))"
    by (simp add: Suc bernoulli_Suc algebra_simps flip: add_2_eq_Suc')
  also have "of_nat (m + 2 choose Suc m) = (of_nat (m + 2) :: 'a)"
    by (simp add: eval_nat_numeral)
  also have "… / … = 1"
    by (metis of_nat_0 of_nat_eq_iff one_eq_divide_iff zero_eq_add_iff_both_eq_0 zero_neq_numeral)
  finally show ?thesis
    by (simp add: Suc)
qed auto
  
lemma sum_binomial_times_bernoulli':
  "(∑k<n. of_nat (n choose k) * bernoulli k :: 'a :: field_char_0) = (if n = 1 then 1 else 0)"
proof (cases n)
  case (Suc m)
  have "(∑k<n. of_nat (n choose k) * bernoulli k :: 'a) =
           (∑k≤m. of_nat (Suc m choose k) * bernoulli k)"
    unfolding Suc lessThan_Suc_atMost ..
  also have "… = (if n = 1 then 1 else 0)" 
    by (subst sum_binomial_times_bernoulli) (simp add: Suc)
  finally show ?thesis .
qed simp_all
  
lemma binomial_unroll:
  "n > 0 ⟹ (n choose k) = (if k = 0 then 1 else ((n - 1) choose (k - 1)) + ((n - 1) choose k))"
  by (auto simp add: gr0_conv_Suc)

lemma sum_unroll:
  "(∑k≤n::nat. f k) = (if n = 0 then f 0 else f n + (∑k≤n - 1. f k))"
  by (cases n) (simp_all add: add_ac)

lemma bernoulli_unroll:
  "n > 0 ⟹ bernoulli n = - 1 / (of_nat n + 1) * (∑k≤n - 1. of_nat (n + 1 choose k) * bernoulli k)"
  by (cases n) (simp add: bernoulli_Suc)+

lemmas bernoulli_unroll_all = binomial_unroll bernoulli_unroll sum_unroll bernpoly_def

lemma bernpoly_1_1: "bernpoly 1 1 = of_real (1/2)"
proof -
  have *: "(1 :: 'a) = of_real 1" by simp
  have "bernpoly 1 (1::'a) = 1 - of_real (1 / 2)"
    by (simp add: bernoulli_unroll_all)
  also have "… = of_real (1 - 1 / 2)"
    by (simp only: *  of_real_diff)
  also have "1 - 1 / 2 = (1 / 2 :: real)"
    by simp
  finally show ?thesis .
qed


subsection ‹Sum of Powers with Bernoulli Polynomials›

(* TODO: Generalisation not possible here because mean-value theorem 
   is only available for reals *)
lemma diff_bernpoly:
  fixes x :: real
  shows "bernpoly n (x + 1) - bernpoly n x = of_nat n * x ^ (n - 1)"
proof (induct n arbitrary: x)
  case 0
  show ?case unfolding bernpoly_def by auto
next
  case (Suc n)
  have "bernpoly (Suc n) (0 + 1) - bernpoly (Suc n) (0 :: real) = 
          (∑k≤n. of_real (real (Suc n choose k) * bernoulli k))"
    unfolding bernpoly_0 unfolding bernpoly_def by simp
  also have "… = of_nat (Suc n) * 0 ^ n"
    by (simp only: of_real_sum [symmetric] sum_binomial_times_bernoulli) simp
  finally have const: "bernpoly (Suc n) (0 + 1) - bernpoly (Suc n) 0 = …"
    by simp

  have hyps': "of_nat (Suc n) * bernpoly n (x + 1) - 
                  of_nat (Suc n) * bernpoly n x = 
                  of_nat n * of_nat (Suc n) * x ^ (n - Suc 0)" for x :: real
    unfolding right_diff_distrib[symmetric] 
    by (subst Suc) (simp_all add: algebra_simps)
  have "((λx. bernpoly (Suc n) (x + 1) - bernpoly (Suc n) x - of_nat (Suc n) * x ^ n) 
           has_field_derivative 0) (at x)" for x :: real
    by (rule derivative_eq_intros refl)+ (insert hyps'[of x], simp add: algebra_simps)
  from integrals_eq[OF const this] show ?case by simp
qed

lemma bernpoly_of_real: "bernpoly n (of_real x) = of_real (bernpoly n x)"
  by (simp add: bernpoly_def)
  
lemma bernpoly_1:
  assumes "n ≠ 1"
  shows   "bernpoly n 1 = of_real (bernoulli n)"
proof -
  have "bernpoly n 1 = (of_real (bernpoly n 1) :: 'a)"
    by (simp add: bernpoly_altdef)
  also have "bernpoly n 1 = (bernoulli n :: real)"
  proof (cases "n ≥ 2")
    case False
    with assms have "n = 0" by auto
    thus ?thesis by (simp add: bernpoly_def)
  next
    case True
    with diff_bernpoly[of n 0] show ?thesis
      by (simp add: power_0_left bernpoly_0)
  qed
  finally show ?thesis
    by simp
qed
  
lemma bernpoly_1': "bernpoly n 1 = of_real (bernoulli' n)"
  using bernpoly_1_1 [where ?'a = 'a]
  by (cases "n = 1") (simp_all add: bernpoly_1 bernoulli'_def)

theorem sum_of_powers: 
  "(∑k≤n::nat. (real k) ^ m) = (bernpoly (Suc m) (n + 1) - bernpoly (Suc m) 0) / (m + 1)"
proof -
  from diff_bernpoly[of "Suc m", simplified] have "(m + (1::real)) * (∑k≤n. (real k) ^ m) = (∑k≤n. bernpoly (Suc m) (real k + 1) - bernpoly (Suc m) (real k))"
    by (auto simp add: sum_distrib_left intro!: sum.cong)
  also have "... = (∑k≤n. bernpoly (Suc m) (real (k + 1)) - bernpoly (Suc m) (real k))"
    by (simp add: add_ac)
  also have "... = bernpoly (Suc m) (n + 1) - bernpoly (Suc m) 0"
    by (simp only: sum_diff[where f="λk. bernpoly (Suc m) (real k)"]) simp
  finally show ?thesis by (auto simp add: field_simps intro!: eq_divide_imp)
qed
  
lemma sum_of_powers_nat_aux: 
  assumes "real a = b / c" "real b' = b" "real c' = c"
  shows   "a = b' div c'"
proof (cases "c = 0")
  case False 
  with assms have "real (a * c') = real b'" by (simp add: field_simps)
  hence "b' = a * c'" by (subst (asm) of_nat_eq_iff) simp
  with False assms show ?thesis by simp
qed (insert assms, simp_all)


subsection ‹Instances for Square And Cubic Numbers›

theorem sum_of_squares: "real (∑k≤n::nat. k ^ 2) = real (2 * n ^ 3 + 3 * n ^ 2 + n) / 6"
  unfolding of_nat_sum of_nat_power sum_of_powers
  by (simp add: bernoulli_unroll_all field_simps power2_eq_square power_numeral_reduce)

corollary sum_of_squares_nat: "(∑k≤n::nat. k ^ 2) = (2 * n ^ 3 + 3 * n ^ 2 + n) div 6"
  by (rule sum_of_powers_nat_aux[OF sum_of_squares]) simp_all

theorem sum_of_cubes: "real (∑k≤n::nat. k ^ 3) = real (n ^ 2 + n) ^ 2 / 4"
  unfolding of_nat_sum of_nat_power sum_of_powers
  by (simp add: bernoulli_unroll_all field_simps power2_eq_square power_numeral_reduce)
                       
corollary sum_of_cubes_nat: "(∑k≤n::nat. k ^ 3) = (n ^ 2 + n) ^ 2 div 4"
  by (rule sum_of_powers_nat_aux[OF sum_of_cubes]) simp_all

end