feat(Algebra/Group /ForwardDiff.lean): add five theorems for forward difference

Mathlib/Algebra/Group/ForwardDiff.lean

@@ -1,7 +1,7 @@
 /-
 Copyright (c) 2024 David Loeffler. All rights reserved.
 Released under Apache 2.0 license as described in the file LICENSE.
-Authors: Giulio Caflisch, David Loeffler
+Authors: Giulio Caflisch, David Loeffler, Yu Shao, Beibei Xiong, Weijie Jiang
 -/
 import Mathlib.Algebra.BigOperators.Pi
 import Mathlib.Algebra.Group.AddChar
@@ -193,3 +193,204 @@ lemma fwdDiff_addChar_eq {M R : Type*} [AddCommMonoid M] [Ring R]
   | succ n IH =>
     simp only [pow_succ, Function.iterate_succ_apply', fwdDiff, IH, ← mul_sub, mul_assoc]
     rw [sub_mul, ← AddChar.map_add_eq_mul, add_comm h x, one_mul]
+
+
+/-!
+## Forward differences of Polynomials
+
+This section develops the theory of forward differences for polynomial functions `P : R → R`,
+where the step size `h` is `1`. We prove several key results:
+
+* `fwdDiff_iter_pow_eq_zero_of_lt`: The `n`-th difference of `x ↦ x^j` is zero if `j < n`.
+* `fwdDiff_iter_eq_factorial`: The `n`-th difference of `x ↦ x^n` is the constant `n!`.
+* `fwdDiff_iter_succ_sum_eq_zero`: The `(d+1)`-th difference of a polynomial of degree `d` is zero.
+* `sum_range_choose_mul_fwdDiff_iter_at_zero`: **Newton's series**
+  for a polynomial, expressing `P(x)` as a sum
+  of its forward differences at `0` weighted by binomial coefficients.
+* `sum_sum_range_eq_sum_range_choose_mul_fwdDiff_iter_at_zero`: A formula for the sum of a
+  polynomial sequence `∑_{i=0..p} P(i)`, which generalizes **Faulhaber's formula**.
+-/
+
+open fwdDiff
+variable {R : Type*} [CommRing R]
+
+/--
+The `n`-th forward difference of the function `x ↦ x^j` is zero if `j < n`.
+This is a building block for showing that the `(p+1)`-th difference of a polynomial of
+degree `p` is zero.
+-/
+theorem fwdDiff_iter_pow_eq_zero_of_lt {j n : ℕ} (h : j < n) :
+    ((fwdDiffₗ R R 1 ^ n) fun x ↦ x ^ j) = 0 := by
+  induction n generalizing j with
+  | zero => contradiction
+  | succ n ih =>
+    rw [pow_succ, Module.End.mul_apply]
+    have : ((fwdDiffₗ R R 1) fun x ↦ x ^ j) =
+        ∑ i ∈ Finset.range j, j.choose i • fun x : R ↦ x ^ i := by
+      ext x
+      simp only [fwdDiffₗ_apply, nsmul_eq_mul, sum_apply, Pi.mul_apply, Pi.natCast_apply, fwdDiff,
+        add_pow, one_pow, Finset.sum_range_succ, mul_one, choose_self, cast_one,
+        add_sub_cancel_right, mul_comm]
+    rw [this, map_sum]
+    exact Finset.sum_eq_zero fun i hi ↦ have _ := Finset.mem_range.1 hi; by
+      rw [map_nsmul, ih (by omega)]; rw [nsmul_zero]
+
+
+/-- The `n`-th forward difference of `x ↦ x^n` is the constant function `n!`. -/
+theorem fwdDiff_iter_eq_factorial {n : ℕ} :
+    ((fwdDiffₗ R R 1 ^ n) fun x ↦ x ^ n) = (fun _ ↦ (n.factorial : R))  := by
+  induction n with
+  | zero => simp only [pow_zero, Module.End.one_apply, factorial_zero, cast_one]
+  | succ n ih =>
+    simp at ih
+    rw [pow_succ, Module.End.mul_apply]
+    have : ((fwdDiffₗ R R 1) fun x ↦ x ^ (n + 1)) =
+        ∑ k ∈ Finset.range (n + 1), (n + 1).choose k • fun x : R ↦ x ^ k := by
+      ext x
+      simp only [fwdDiffₗ_apply, nsmul_eq_mul, sum_apply, Pi.mul_apply, Pi.natCast_apply,
+        fwdDiff, add_pow, one_pow, Finset.sum_range_succ, Nat.choose_self, cast_one, mul_one,
+        add_sub_assoc, sub_self, add_zero]
+      simp only [choose_succ_self_right, cast_add, cast_one, mul_comm]
+    rw [this, map_sum, Nat.factorial_succ, Nat.cast_mul]; ext x
+    rw [funext_iff] at ih
+    simp only [← ih x, Finset.sum_range_succ, choose_succ_self_right, cast_add, cast_one]
+    have : (fwdDiffₗ R R 1 ^ n) ((n + 1) • fun x ↦ x ^ n) =
+      ((n : R → R) + 1) * (fwdDiffₗ R R 1 ^ n) fun x ↦ x ^ n := by
+      rw [map_nsmul]
+      simp only [nsmul_eq_mul, cast_add, cast_one]
+    simp only [this, Pi.add_apply, sum_apply, Pi.mul_apply,
+      Pi.natCast_apply, Pi.one_apply, add_eq_right]
+    exact Finset.sum_eq_zero fun i hi ↦ have _ := Finset.mem_range.1 hi; by
+      rw [map_nsmul, fwdDiff_iter_pow_eq_zero_of_lt
+        (by linarith), Pi.smul_apply, Pi.zero_apply, smul_zero]
+
+
+/--
+The `(n+1)`-th forward difference of a polynomial of degree at most `n` is zero.
+A polynomial `P(x) = ∑_{k=0..n} aₖ xᵏ` has `Δ^[n+1] P = 0`.
+-/
+theorem fwdDiff_iter_succ_sum_eq_zero {n : ℕ} (a : ℕ → R):
+    ((fwdDiffₗ R R 1 ^ (n + 1)) fun x ↦ ∑ k ∈ Finset.range (n + 1), a k * (x ^ k)) = 0 := by
+  induction n with
+  | zero =>
+    unfold fwdDiffₗ
+    simp; ext x
+    simp only [Pi.zero_apply]
+  | succ n ih =>
+    rw [pow_succ, Module.End.mul_apply]
+    have : ((fwdDiffₗ R R 1 ^ (n + 1)) (fun x => ∑ k ∈ range (n + 1), a k * x ^ k)) =
+      ∑ k ∈ range (n + 1), (fwdDiffₗ R R 1 ^ (n + 1))
+      ((fun x : R ↦ a k) * (fun x : R ↦  x ^ k)) := by
+      ext x; simp
+      have : (fun x => ∑ k ∈ range (n + 1), a k * x ^ k) =
+        (∑ k ∈ range (n + 1), fun x => a k * x ^ k) := by
+        ext x; simp only [Finset.sum_apply]
+      simp only [this, map_sum, sum_apply]
+      congr 1
+    rw [this] at ih
+    have : ((fwdDiffₗ R R 1) fun x ↦ ∑ k ∈ range (n + 1 + 1), a k * x ^ k) =
+      ∑ k ∈ range (n + 1 + 1), a k • (fun x : R ↦  (x + 1) ^ k) -
+      ∑ k ∈ range (n + 1 + 1),a k • fun x : R ↦ x ^ k := by
+      ext x; simp
+      unfold fwdDiff
+      rfl
+    rw [this]
+    simp only [map_sub, map_sum, coe_fwdDiffₗ_pow]
+    ext x
+    rw [Finset.sum_range_succ]
+    nth_rw 2 [Finset.sum_range_succ]
+    simp only [fwdDiff_iter_const_smul, Pi.sub_apply, Pi.add_apply, sum_apply, Pi.smul_apply,
+      smul_eq_mul]
+    rw [← add_sub, ← coe_fwdDiffₗ_pow]
+    have : ((fwdDiffₗ R R 1 ^ (n + 1)) fun x => x ^ (n + 1)) (x + 1) =
+        ((fwdDiffₗ R R 1 ^ (n + 1)) fun x => (x + 1) ^ (n + 1)) x := by
+      simp only [coe_fwdDiffₗ_pow, fwdDiff_iter_eq_sum_shift, fwdDiff_iter_eq_sum_shift,
+        Int.reduceNeg, nsmul_eq_mul, mul_one, zsmul_eq_mul, Int.cast_mul, Int.cast_pow,
+        Int.cast_neg, Int.cast_one, Int.cast_natCast,
+        add_assoc (b := (1 : R)), add_comm (a := (1 : R)), ← add_assoc]
+    simp only [← this, fwdDiff_iter_eq_factorial (n := n + 1) (R := R), sub_add_cancel_right, ←
+      sub_eq_add_neg, ← Finset.sum_sub_distrib]
+    exact Finset.sum_eq_zero fun k hk ↦ have _ := Finset.mem_range.1 hk; by
+      rw [← mul_sub, fwdDiff_iter_pow_eq_zero_of_lt (by linarith) ]
+      simp only [Pi.zero_apply, sub_zero]
+      have : (fwdDiffₗ R R 1 ^ (n + 1)) (fun x => (x + 1) ^ k) x =
+      (fwdDiffₗ R R 1 ^ (n + 1)) (fun x => x ^ k) (x + 1) := by
+        rw [coe_fwdDiffₗ_pow, fwdDiff_iter_eq_sum_shift, fwdDiff_iter_eq_sum_shift]
+        congr 1; ext k
+        simp only [Int.reduceNeg, nsmul_eq_mul, mul_one,
+          zsmul_eq_mul, Int.cast_mul, Int.cast_pow, Int.cast_neg, Int.cast_one, Int.cast_natCast]
+        ring
+      rw [this, fwdDiff_iter_pow_eq_zero_of_lt (by linarith)]
+      simp only [Pi.zero_apply, mul_zero]
+
+/--
+**Newton's series** for a polynomial function.
+Any function `f` defined by a polynomial can be expressed as a sum of its forward
+differences at `0`, weighted by binomial coefficients.
+`f(x) = ∑_{k=0..p} (p choose k) * (Δ^k f)(0)`.
+-/
+theorem sum_range_choose_mul_fwdDiff_iter_at_zero {n p : ℕ} (a : ℕ → R):
+    ∑ k ∈ Finset.range (n + 1), a k * (p ^ k) =
+    ∑ k ∈ Finset.range (p + 1), p.choose k *
+      ((fwdDiffₗ R R 1 ^ k) fun x ↦ ∑ k ∈ Finset.range (n + 1), a k * (x ^ k)) 0 := by
+  obtain h := shift_eq_sum_fwdDiff_iter (n := p) (y := 0)
+    (f := (fun x => ∑ i ∈ Finset.range (n + 1), a i * (x ^ i))) (h := 1)
+  simp only [mul_one, zero_add, nsmul_eq_mul] at h
+  rw [h]
+  exact Finset.sum_congr rfl fun k hk ↦ by
+    have _ := Finset.mem_range.1 hk
+    congr 1
+    simp only [coe_fwdDiffₗ_pow, fwdDiff_iter_eq_sum_shift, Int.reduceNeg, nsmul_eq_mul,
+      mul_one, zsmul_eq_mul, Int.cast_mul, Int.cast_pow, Int.cast_neg, Int.cast_one,
+      Int.cast_natCast]
+
+/--
+A formula for the sum of a polynomial sequence `∑_{k=0..p} P(k)`, which
+generalizes **Faulhaber's formula**.
+-/
+theorem sum_sum_range_eq_sum_range_choose_mul_fwdDiff_iter_at_zero {p n : ℕ} (a : ℕ → R) :
+    ∑ k ∈ Finset.range (p + 1), (∑ m ∈ Finset.range (n + 1), a m * k ^ m) =
+    ∑ k ∈ Finset.range (p + 1), ((p + 1).choose (k + 1)) *
+      (((fwdDiffₗ R R 1 ^ k) fun y ↦ ∑ i ∈ Finset.range (n + 1), a i * y ^ i) 0) := by
+  conv => enter [1, 2, x]; rw [sum_range_choose_mul_fwdDiff_iter_at_zero]; simp
+  have sum_extend_inner_range : ∑ x ∈ Finset.range (p + 1), ∑ k ∈ Finset.range (x + 1),
+    ↑(x.choose k) * ((fwdDiffₗ R R 1 ^ k) fun x ↦ ∑ m ∈ Finset.range (n + 1), a m * ↑x ^ m) 0 =
+    ∑ x ∈ Finset.range (p + 1), ∑ k ∈ Finset.range (p + 1),
+    ↑(x.choose k) * ((fwdDiffₗ R R 1 ^ k) fun x ↦ ∑ m ∈ Finset.range (n + 1), a m * ↑x ^ m) 0 := by
+    apply Finset.sum_congr rfl
+    intro x hx
+    have sum_sum_eq_zero : ∑ k ∈ Finset.Ico (x + 1) (p + 1), ↑(x.choose k) *
+      ((fwdDiffₗ R R 1 ^ k) fun x ↦ ∑ m ∈ Finset.range (n + 1), a m * x ^ m) 0 = 0 := by
+      rw [Finset.sum_Ico_eq_sum_range]
+      simp
+      simp at hx
+      have : ∑ k ∈ Finset.range (p - x), 0 = (0 : R) := by simp only [Finset.sum_const_zero]
+      rw [← this]
+      apply Finset.sum_congr rfl
+      intro y hy; simp only [mem_range] at hy
+      have : x + 1 + y > x := by omega
+      rw [Nat.choose_eq_zero_of_lt this]
+      simp
+    nth_rw 1 3 [Finset.range_eq_Ico]
+    have hx' : 0 ≤ (x + 1) := by omega
+    have hxp' : x + 1 ≤ p + 1 := by
+      simp only [mem_range] at hx
+      omega
+    rw [← Finset.sum_Ico_consecutive _ hx' hxp', sum_sum_eq_zero, add_zero]
+  rw [sum_extend_inner_range, Finset.sum_comm]
+  simp_rw [← Finset.sum_mul]
+  apply Finset.sum_congr rfl
+  intro k hk; simp only [mem_range] at hk
+  congr 1
+  norm_cast
+  have hk1 : 0 ≤ k := by omega
+  have hk2 : k ≤ p + 1 := by omega
+  simp_rw [← Ico_zero_eq_range, ← Finset.sum_Ico_consecutive _ hk1 hk2]
+  have l1 : ∑ x ∈ Ico 0 k, x.choose k = 0 := by
+    simp only [Ico_zero_eq_range, sum_eq_zero_iff, mem_range]
+    intro _ _
+    exact choose_eq_zero_iff.mpr (by omega)
+  have l2 : ∑ x ∈ Ico k (p + 1), x.choose k = (p + 1).choose (k + 1) := by
+    rw [Finset.sum_Ico_succ_top (by omega), Finset.sum_Ico_add_eq_sum_Icc (by omega),
+      Nat.sum_Icc_choose]
+  simp_rw [l1, l2, Nat.zero_add]