NewtonDeterminant.lean

1 /-
2 Copyright (c) 2026 Komorebi P. All rights reserved.
3 Authors: Komorebi P.
4 -/
5 import KomorebiPaperIV.ReciprocalCompletion
6 import Mathlib.LinearAlgebra.Matrix.Block
7 import Mathlib.Data.Nat.Factorial.BigOperators
8 
9 namespace Komorebi.PaperIV
10 noncomputable section
11 open Matrix Finset
12 variable {K : Type*} [Field K] [CharZero K]
13 
14 /-- Newton's lower Hessenberg matrix, including its negative superdiagonal. -/
15 def newtonHessenberg (p : ℕ → K) (n : ℕ) : Matrix (Fin n) (Fin n) K :=
16   fun i j => (if j.val ≤ i.val then p (i.val+1-j.val) else 0) -
17     (if j.val = i.val+1 then (i.val+1 : K) else 0)
18 
19 omit [CharZero K] in
20 private theorem newton_column_sum (p b : ℕ → K) (n : ℕ)
21     (hb : ∀ i, ∑ j ∈ range (i+1), b j * p (i+1-j) = (i+1 : K)*b (i+1))
22     (i : Fin (n+1)) :
23     (∑ j : Fin (n+1), b j.val * newtonHessenberg p (n+1) i j) =
24       if i = Fin.last n then (n+1 : K)*b (n+1) else 0 := by
25   have hi : i.val+1 ≤ n+1 := i.isLt
26   simp only [newtonHessenberg, mul_sub, sum_sub_distrib]
27   rw [Fin.sum_univ_eq_sum_range (f := fun j => b j * (if j ≤ i.val then p (i.val+1-j) else 0)),
28     Fin.sum_univ_eq_sum_range (f := fun j => b j * (if j = i.val+1 then (i.val+1 : K) else 0))]
29   have hl : (∑ j ∈ range (n+1), b j * (if j ≤ i.val then p (i.val+1-j) else 0)) =
30       ∑ j ∈ range (i.val+1), b j*p (i.val+1-j) := by
31     symm
32     calc
33       _ = ∑ j ∈ range (i.val+1), b j * (if j ≤ i.val then p (i.val+1-j) else 0) := by
34         apply sum_congr rfl
35         intro j hj
36         simp only [show j ≤ i.val by have := mem_range.mp hj; omega, ite_true]
37       _ = _ := sum_subset (range_mono hi) (by
38         intro j _ hj
39         simp only [show ¬ j ≤ i.val by simpa using hj, ite_false, mul_zero])
40   rw [hl, hb]
41   by_cases hiLast : i = Fin.last n
42   · subst i
43     simp only [Fin.val_last, ite_true]
44     have hz : (∑ j ∈ range (n+1), b j * (if j = n+1 then (n+1 : K) else 0)) = 0 := by
45       apply sum_eq_zero
46       intro j hj
47       simp only [show j ≠ n+1 by have := mem_range.mp hj; omega, ite_false, mul_zero]
48     rw [hz, sub_zero]
49   · have hil : i.val+1 < n+1 := by
50       have hn : i.val ≠ n := by intro h; apply hiLast; exact Fin.ext h
51       omega
52     rw [sum_eq_single (i.val+1)]
53     · simp [hiLast, mul_comm]
54     · intro j _ hj
55       simp only [hj, ite_false, mul_zero]
56     · exact fun h => (h (mem_range.mpr hil)).elim
57 
58 /-- The exact determinant solution of Newton's recurrence. -/
59 theorem newtonHessenberg_det (p b : ℕ → K) (hb0 : b 0 = 1)
60     (hb : ∀ i, ∑ j ∈ range (i+1), b j*p (i+1-j) = (i+1 : K)*b (i+1))
61     (n : ℕ) : (newtonHessenberg p n).det = (n.factorial : K)*b n := by
62   classical
63   cases n with
64   | zero => simp [hb0, Matrix.det_isEmpty]
65   | succ n =>
66     let H := newtonHessenberg p (n+1)
67     let C := H.updateCol 0 (fun i => if i = Fin.last n then (n+1 : K)*b (n+1) else 0)
68     have hc : C.det = H.det := by
69       have h := Matrix.det_updateCol_sum H 0 (fun j => b j.val)
70       simp only [Fin.val_zero, hb0, smul_eq_mul, one_mul] at h
71       have he : (fun i => if i = Fin.last n then (n+1 : K)*b (n+1) else 0) =
72           (fun i => ∑ j, b j.val * H i j) := by
73         funext i
74         exact (newton_column_sum p b n hb i).symm
75       change (H.updateCol 0 _).det = H.det
76       rw [he]
77       exact h
78     have hm : C.submatrix (Fin.last n).succAbove Fin.succ =
79         fun i j : Fin n => newtonHessenberg p (n+1) i.castSucc j.succ := by
80       ext i j
81       simp [C, H, Matrix.submatrix, Matrix.updateCol, Fin.succ_ne_zero]
82     have ht : (C.submatrix (Fin.last n).succAbove Fin.succ).IsLowerTriangular := by
83       rw [hm]
84       intro i j hij
85       simp only [newtonHessenberg, Fin.val_castSucc, Fin.val_succ]
86       have hij' : i.val < j.val := hij
87       simp [show ¬ j.val+1 ≤ i.val by omega, show j.val ≠ i.val by omega]
88     have hd : (C.submatrix (Fin.last n).succAbove Fin.succ).det =
89         (-1 : K)^n * (n.factorial : K) := by
90       rw [Matrix.det_of_isLowerTriangular _ ht, hm]
91       simp only [newtonHessenberg, Fin.val_castSucc, Fin.val_succ,
92         Nat.not_succ_le_self, ite_false, ite_true, zero_sub]
93       rw [Finset.prod_neg]
94       simp only [card_univ, Fintype.card_fin]
95       congr 1
96       rw [Fin.prod_univ_eq_prod_range (f := fun j => (j : K)+1), Nat.factorial_eq_prod_range_add_one, Nat.cast_prod]
97       apply prod_congr rfl
98       intro j _
99       simp
100     rw [← hc, Matrix.det_succ_column_zero, sum_eq_single (Fin.last n)]
101     · simp only [C, Matrix.updateCol_self, Fin.val_last, ite_true]
102       rw [hd, Nat.factorial_succ, Nat.cast_mul, Nat.cast_add, Nat.cast_one]
103       have hs : (-1 : K)^n * (-1)^n = 1 := by rw [← mul_pow]; simp
104       calc
105         _ = ((-1 : K)^n * (-1)^n) * ((n+1 : K)*(n.factorial : K)*b (n+1)) := by ring
106         _ = _ := by rw [hs, one_mul]
107     · intro i _ hi
108       simp [C, hi]
109     · simp
110 
111 /-- The Newton determinant printed in Theorem 7.6. -/
112 theorem reciprocalOrdCoeff_newton_determinant (N n : ℕ) (r : K) :
113     reciprocalOrdCoeff N n r = (n.factorial : K)⁻¹ *
114       (newtonHessenberg (fun k => r*kappaCoeff N k) n).det := by
115   have hb0 : reciprocalOrdCoeff N 0 r = 1 := by simp [reciprocalOrdCoeff]
116   have hb (i : ℕ) : ∑ j ∈ range (i+1), reciprocalOrdCoeff N j r *
117       (r*kappaCoeff N (i+1-j)) = (i+1 : K)*reciprocalOrdCoeff N (i+1) r := by
118     have h := reciprocalOrdCoeff_newton N r i
119     rw [Finset.Nat.sum_antidiagonal_eq_sum_range_succ_mk] at h
120     rw [h, mul_sum]
121     apply sum_congr rfl
122     intro j hj
123     have he : i-j+1=i+1-j := by have := mem_range.mp hj; omega
124     rw [he]
125     ring
126   rw [newtonHessenberg_det _ _ hb0 hb]
127   have hf : (n.factorial : K) ≠ 0 := by exact_mod_cast n.factorial_ne_zero
128   field_simp
129 
130 end
131 end Komorebi.PaperIV

Back to source controls