feat(number_theory/bernoulli): Faulhaber's theorem (#6409)

Co-authored-by Fabian Kruse




Co-authored-by: Moritz Firsching <firsching@google.com>
Co-authored-by: Benjamin Davidson <68528197+benjamindavidson@users.noreply.github.com>

docs/100.yaml

@@ -271,6 +271,8 @@
   title  : Fourier Series
 77:
   title  : Sum of kth powers
+  decl   : sum_range_pow
+  author : mathlib (Moritz Firsching, Fabian Kruse, Ashvni Narayanan)
 78:
   title  : The Cauchy-Schwarz Inequality
   decls  :

docs/references.bib

@@ -508,6 +508,17 @@ @Book{            MM92
   publisher     = {Springer Science \& Business Media}
 }
 
+@Article{         MR0236876,
+  title         = {A new proof that metric spaces are paracompact},
+  author        = {Mary Ellen Rudin},
+  year          = 1969,
+  journal       = {Proc. Amer. Math. Soc.},
+  volume        = {20},
+  pages         = {603},
+  mrnumber      = {0236876},
+  doi           = {10.1090/S0002-9939-1969-0236876-3}
+}
+
 @Article{         MR1167694,
   author        = {Blass, Andreas},
   title         = {A game semantics for linear logic},
@@ -620,17 +631,6 @@ @InCollection{    MR541021
   mrreviewer    = {Daniel Lazard}
 }
 
-@Article{         MR0236876,
-  title         = {A new proof that metric spaces are paracompact},
-  author        = {Mary Ellen Rudin},
-  year          = 1969,
-  journal       = {Proc. Amer. Math. Soc.},
-  volume        = {20},
-  pages         = {603},
-  mrnumber      = {0236876},
-  doi           = {10.1090/S0002-9939-1969-0236876-3}
-}
-
 @Article{         MR577178,
   author        = {Cirel\cprime son, B. S.},
   title         = {Quantum generalizations of {B}ell's inequality},
@@ -649,6 +649,21 @@ @Article{         MR577178
   url           = {https://doi.org/10.1007/BF00417500}
 }
 
+@Article{         orosi2018faulhaber,
+  author        = {Greg {Orosi}},
+  title         = {{A simple derivation of Faulhaber's formula}},
+  fjournal      = {{Applied Mathematics E-Notes}},
+  journal       = {{Appl. Math. E-Notes}},
+  issn          = {1607-2510/e},
+  volume        = {18},
+  pages         = {124--126},
+  year          = {2018},
+  publisher     = {Tsing Hua University, Department of Mathematics, Hsinchu},
+  language      = {English},
+  msc2010       = {41A58 30K05},
+  zbl           = {1411.41023}
+}
+
 @Misc{            ponton2020chebyshev,
   title         = {Roots of {C}hebyshev polynomials: a purely algebraic
                   approach},

src/number_theory/bernoulli.lean

@@ -52,9 +52,9 @@ then defined as `bernoulli = (-1)^n * bernoulli'`.
 -/
 
 open_locale big_operators
+open_locale nat
 open nat
 open finset
-open_locale nat
 
 /-!
 
@@ -66,27 +66,27 @@ open_locale nat
 the $n$-th Bernoulli number $B_n$ is defined recursively via
 $$B_n = 1 - \sum_{k < n} \binom{n}{k}\frac{B_k}{n+1-k}$$ -/
 def bernoulli' : ℕ → ℚ :=
-well_founded.fix nat.lt_wf
+well_founded.fix lt_wf
   (λ n bernoulli', 1 - ∑ k : fin n, n.choose k / (n - k + 1) * bernoulli' k k.2)
 
 lemma bernoulli'_def' (n : ℕ) :
   bernoulli' n = 1 - ∑ k : fin n, (n.choose k) / (n - k + 1) * bernoulli' k :=
 well_founded.fix_eq _ _ _
 
 lemma bernoulli'_def (n : ℕ) :
-  bernoulli' n = 1 - ∑ k in finset.range n, (n.choose k) / (n - k + 1) * bernoulli' k :=
+  bernoulli' n = 1 - ∑ k in range n, (n.choose k) / (n - k + 1) * bernoulli' k :=
 by { rw [bernoulli'_def', ← fin.sum_univ_eq_sum_range], refl }
 
 lemma bernoulli'_spec (n : ℕ) :
-  ∑ k in finset.range n.succ, (n.choose (n - k) : ℚ) / (n - k + 1) * bernoulli' k = 1 :=
+  ∑ k in range n.succ, (n.choose (n - k) : ℚ) / (n - k + 1) * bernoulli' k = 1 :=
 begin
-  rw [finset.sum_range_succ, bernoulli'_def n, nat.sub_self],
-  conv in (nat.choose _ (_ - _)) { rw choose_symm (le_of_lt (finset.mem_range.1 H)) },
+  rw [sum_range_succ, bernoulli'_def n, nat.sub_self],
+  conv in (nat.choose _ (_ - _)) { rw choose_symm (le_of_lt (mem_range.1 H)) },
   simp only [one_mul, cast_one, sub_self, sub_add_cancel, choose_zero_right, zero_add, div_one],
 end
 
 lemma bernoulli'_spec' (n : ℕ) :
-  ∑ k in finset.nat.antidiagonal n,
+  ∑ k in nat.antidiagonal n,
   ((k.1 + k.2).choose k.2 : ℚ) / (k.2 + 1) * bernoulli' k.1 = 1 :=
 begin
   refine ((nat.sum_antidiagonal_eq_sum_range_succ_mk _ n).trans _).trans (bernoulli'_spec n),
@@ -105,24 +105,24 @@ section examples
 
 open finset
 
-@[simp] lemma bernoulli'_zero  : bernoulli' 0 = 1   := rfl
+@[simp] lemma bernoulli'_zero : bernoulli' 0 = 1 := rfl
 
-@[simp] lemma bernoulli'_one   : bernoulli' 1 = 1/2 :=
+@[simp] lemma bernoulli'_one : bernoulli' 1 = 1/2 :=
 begin
     rw [bernoulli'_def, sum_range_one], norm_num
 end
 
-@[simp] lemma bernoulli'_two   : bernoulli' 2 = 1/6 :=
+@[simp] lemma bernoulli'_two : bernoulli' 2 = 1/6 :=
 begin
   rw [bernoulli'_def, sum_range_succ, sum_range_one], norm_num
 end
 
-@[simp] lemma bernoulli'_three : bernoulli' 3 = 0   :=
+@[simp] lemma bernoulli'_three : bernoulli' 3 = 0 :=
 begin
   rw [bernoulli'_def, sum_range_succ, sum_range_succ, sum_range_one], norm_num
 end
 
-@[simp] lemma bernoulli'_four  : bernoulli' 4 = -1/30 :=
+@[simp] lemma bernoulli'_four : bernoulli' 4 = -1/30 :=
 begin
   rw [bernoulli'_def, sum_range_succ, sum_range_succ, sum_range_succ, sum_range_one],
   rw (show nat.choose 4 2 = 6, from dec_trivial), -- shrug
@@ -134,7 +134,7 @@ end examples
 open nat finset
 
 @[simp] lemma sum_bernoulli' (n : ℕ) :
-  ∑ k in finset.range n, (n.choose k : ℚ) * bernoulli' k = n :=
+  ∑ k in range n, (n.choose k : ℚ) * bernoulli' k = n :=
 begin
   cases n with n, { simp },
   rw [sum_range_succ, bernoulli'_def],
@@ -204,15 +204,16 @@ open ring_hom
 /-- Odd Bernoulli numbers (greater than 1) are zero. -/
 theorem bernoulli'_odd_eq_zero {n : ℕ} (h_odd : odd n) (hlt : 1 < n) : bernoulli' n = 0 :=
 begin
+  --have f := bernoulli'_power_series,
   have f : power_series.mk (λ n, (bernoulli' n / n!)) * (exp ℚ - 1) = X * exp ℚ,
   { simpa [bernoulli'_power_series] using bernoulli'_power_series_mul_exp_sub_one ℚ },
   have g : eval_neg_hom (mk (λ (n : ℕ), bernoulli' n / ↑(n!)) * (exp ℚ - 1)) * (exp ℚ) =
     (eval_neg_hom (X * exp ℚ)) * (exp ℚ) := by congr',
   rw [map_mul, map_sub, map_one, map_mul, mul_assoc, sub_mul, mul_assoc (eval_neg_hom X) _ _,
     mul_comm (eval_neg_hom (exp ℚ)) (exp ℚ), exp_mul_exp_neg_eq_one, eval_neg_hom_X, mul_one,
     one_mul] at g,
-  suffices h : (mk (λ (n : ℕ), bernoulli' n / ↑(n!)) - eval_neg_hom (mk (λ (n : ℕ),
-    bernoulli' n / ↑(n!))) ) * (exp ℚ - 1) = X * (exp ℚ - 1),
+  suffices h : (mk (λ (n : ℕ), bernoulli' n / ↑n!) - eval_neg_hom (mk (λ (n : ℕ),
+    bernoulli' n / ↑n!)) ) * (exp ℚ - 1) = X * (exp ℚ - 1),
   { rw [mul_eq_mul_right_iff] at h,
     cases h,
     { simp only [eval_neg_hom, rescale, coeff_mk, coe_mk, power_series.ext_iff,
@@ -239,9 +240,9 @@ def bernoulli (n : ℕ) : ℚ := (-1)^n * (bernoulli' n)
 lemma bernoulli'_eq_bernoulli (n : ℕ) : bernoulli' n = (-1)^n * bernoulli n :=
 by simp [bernoulli, ← mul_assoc, ← pow_two, ← pow_mul, mul_comm n 2, pow_mul]
 
-@[simp] lemma bernoulli_zero  : bernoulli 0 = 1 := rfl
+@[simp] lemma bernoulli_zero : bernoulli 0 = 1 := rfl
 
-@[simp] lemma bernoulli_one   : bernoulli 1 = -1/2 :=
+@[simp] lemma bernoulli_one : bernoulli 1 = -1/2 :=
 by norm_num [bernoulli, bernoulli'_one]
 
 theorem bernoulli_eq_bernoulli'_of_ne_one {n : ℕ} (hn : n ≠ 1) : bernoulli n = bernoulli' n :=
@@ -328,3 +329,65 @@ begin
   rw [cast_mul, cast_mul, mul_div_mul_right _ _ hj, add_choose, cast_dvd_char_zero],
   exact factorial_mul_factorial_dvd_factorial_add i j,
 end
+
+
+section faulhaber
+
+/-- Faulhaber's theorem relating the sum of of p-th powers to the Bernoulli numbers.
+See https://proofwiki.org/wiki/Faulhaber%27s_Formula and [orosi2018faulhaber] for
+the proof provided here. -/
+theorem sum_range_pow (n p : ℕ) :
+  ∑ k in range n, (k : ℚ) ^ p =
+    ∑ i in range (p + 1), bernoulli i * (p + 1).choose i * n ^ (p + 1 - i) / (p + 1) :=
+begin
+  -- trivial fact about cast factorials
+  have hne : ∀ m : ℕ, (m! : ℚ) ≠ 0 := λ m, by exact_mod_cast factorial_ne_zero m,
+  -- the Cauchy product of two power series
+  have h_cauchy : mk (λ p, bernoulli p / ↑p!) * mk (λ q, coeff ℚ (q + 1) (exp ℚ ^ n))
+                = mk (λ p, ∑ i in range (p + 1),
+                      bernoulli i * ↑((p + 1).choose i) * ↑n ^ (p + 1 - i) / ↑(p + 1)!),
+  { ext q,
+    let f := λ a, λ b, bernoulli a / ↑a! * coeff ℚ (b + 1) (exp ℚ ^ n),
+    -- key step: using `power_series.coeff_mul` and then rewriting sums
+    simp only [coeff_mul, coeff_mk, cast_mul, nat.sum_antidiagonal_eq_sum_range_succ f],
+    refine sum_congr rfl _,
+    simp_intros m h only [finset.mem_range],
+    simp only [f, exp_pow_eq_rescale_exp, rescale, one_div, coeff_mk, ring_hom.coe_mk, coeff_exp,
+              ring_hom.id_apply, cast_mul, rat.algebra_map_rat_rat],
+    -- manipulate factorials and binomial coefficients
+    rw [choose_eq_factorial_div_factorial h.le, eq_comm, div_eq_iff (hne q.succ), succ_eq_add_one,
+        mul_assoc _ _ ↑q.succ!, mul_comm _ ↑q.succ!, ←mul_assoc, div_mul_eq_mul_div,
+        mul_comm (n ^ (q - m + 1) : ℚ), ←mul_assoc _ _ (n ^ (q - m + 1) : ℚ), ←one_div, mul_one_div,
+        div_div_eq_div_mul, ←nat.sub_add_comm (le_of_lt_succ h), cast_dvd, cast_mul],
+    { ring },
+    { exact factorial_mul_factorial_dvd_factorial h.le },
+    { simp [hne] } },
+  -- same as our goal except we pull out `p!` for convenience
+  have hps : ∑ k in range n, ↑k ^ p
+          = (∑ i in range (p + 1), bernoulli i * ↑((p + 1).choose i) * ↑n ^ (p + 1 - i) / ↑(p + 1)!)
+            * ↑p!,
+  { suffices : power_series.mk (λ p, ∑ k in range n, ↑k ^ p * algebra_map ℚ ℚ (↑p!)⁻¹)
+             = power_series.mk (λ p, ∑ i in range (p + 1),
+                                bernoulli i * ↑((p + 1).choose i) * ↑n ^ (p + 1 - i) / ↑(p + 1)!),
+    { rw [← div_eq_iff (hne p), div_eq_mul_inv, sum_mul],
+      rw power_series.ext_iff at this,
+      simpa using this p },
+    -- the power series `exp ℚ - 1` is non-zero, a fact we need in order to use `mul_right_inj'`
+    have hexp : exp ℚ - 1 ≠ 0,
+    { simp only [exp, power_series.ext_iff, ne, not_forall],
+      use 1,
+      simp },
+    have h_r : exp ℚ ^ n - 1 = X * mk (λ p, coeff ℚ (p + 1) (exp ℚ ^ n)),
+    { have h_const : C ℚ (constant_coeff ℚ (exp ℚ ^ n)) = 1 := by simp,
+      rw [←h_const, sub_const_eq_X_mul_shift] },
+    -- key step: a chain of equalities of power series
+    rw [←mul_right_inj' hexp, mul_comm, ←exp_pow_sum, ←geom_series_def, geom_sum_mul, h_r,
+        ←bernoulli_power_series_mul_exp_sub_one, bernoulli_power_series, mul_right_comm],
+    simp [h_cauchy, mul_comm] },
+  -- the rest is showing that `hps` can be massaged into our goal
+  rw [hps, sum_mul],
+  refine sum_congr rfl (λ x hx, _),
+  field_simp [mul_right_comm _ ↑p!, ←mul_assoc _ _ ↑p!, cast_add_one_ne_zero, hne],
+end
+
+end faulhaber