Theme/Product.lean — FrcCore

1import FrcCore.Theme.Cayley2/-!3# FrcCore.Theme.Product — the determinant of a product (the operator theme)45`det (A · B) = det A · det B` for square matrices over `K`, by elementary row operations. An operation — the addition of6a multiple of one row to another — is a datum, `Op`; a list of them applies to a matrix, `apply`, leaves the determinant7unchanged (`det_apply`, from `Theme/Det`'s `det_addRow`) and commutes with the right multiplication by any matrix8(`apply_mul`: the rows of `M · B` are the rows of `M` times `B`, a linear map of the row). Gauss–Jordan elimination9(`gj`) reads its operations off `A` — no swaps: a pivot below the current row is added to it — and ends either in a10diagonal matrix with nonzero diagonal or in a matrix with a zero row (`gj_spec`); applied to `A · B` the same operations11give `D · B`, whose rows are the scaled rows of `B`, or a matrix with a zero row. The two cases give the theorem12(`det_mul`). The identities rest on `Theme/Dirac`'s `Mat.mul` (`dot`, a left fold, is `Theme/Elim`'s `dotL`). No axioms.13-/1415namespace FRC.Dirac.Det1617open FRC FRC.Extension FRC.Dirac FRC.Dirac.Hermitian18open FRC.Dirac.Mat (nth)19open FRC.Shell (Frame)20open FRC.Shell.Frame (IsSquare)2122variable {p : Nat} [Pos p] {ν : Shell p}2324/-! ## The product, row by row -/2526/-- The row vector times the matrix: the `i`-th row of `M · B` is `rowMul (row M i) B`. -/27def rowMul (r : List (Ext p ν)) (B : Mat p ν) : List (Ext p ν) := (Mat.transpose B).map fun c => Mat.dot r c2829theorem mul_eq_map (M B : Mat p ν) : M * B = M.map fun r => rowMul r B := rfl3031theorem length_mul (M B : Mat p ν) : (M * B).length = M.length := length_map _ M3233theorem row_mul (M B : Mat p ν) (i : Nat) (hi : i < M.length) : row (M * B) i = rowMul (row M i) B :=34  nth_map_lt (fun r => rowMul r B) [] [] M i hi3536theorem foldl_add_eq (l : List (Ext p ν)) : ∀ a : Ext p ν, l.foldl (· + ·) a = a + l.foldl (· + ·) 0 := by37  induction l with38  | nil => intro a; exact (add_zero a).symm39  | cons x l ih =>40    intro a41    show l.foldl (· + ·) (a + x) = a + l.foldl (· + ·) (0 + x)42    rw [ih (a + x), ih (0 + x), zero_add, add_assoc]4344theorem dot_eq_dotL : ∀ r c : List (Ext p ν), Mat.dot r c = dotL r c45  | [], _ => rfl46  | _ :: _, [] => rfl47  | a :: r, x :: c => by48    show (List.zipWith (· * ·) r c).foldl (· + ·) (0 + a * x) = a * x + dotL r c49    rw [foldl_add_eq, zero_add]; exact congrArg (a * x + ·) (dot_eq_dotL r c)5051theorem zipWith_map_same {α β : Type} (f g : α → β) (h : β → β → β) : ∀ l : List α,52    List.zipWith h (l.map f) (l.map g) = l.map fun x => h (f x) (g x)53  | [] => rfl54  | _ :: l => congrArg (_ :: ·) (zipWith_map_same f g h l)5556theorem rowMul_radd (r s : List (Ext p ν)) (B : Mat p ν) (h : r.length = s.length) :57    rowMul (radd r s) B = radd (rowMul r B) (rowMul s B) := by58  show (Mat.transpose B).map (fun c => Mat.dot (radd r s) c) =59    radd ((Mat.transpose B).map fun c => Mat.dot r c) ((Mat.transpose B).map fun c => Mat.dot s c)60  rw [← zipWith_add_eq_radd ((Mat.transpose B).map fun c => Mat.dot r c), zipWith_map_same]61  exact map_congr' (fun c => by rw [dot_eq_dotL, dot_eq_dotL, dot_eq_dotL, dotL_radd r s c h]) _6263theorem rowMul_rsmul (k : Ext p ν) (r : List (Ext p ν)) (B : Mat p ν) : rowMul (rsmul k r) B = rsmul k (rowMul r B) := by64  show (Mat.transpose B).map (fun c => Mat.dot (rsmul k r) c) = ((Mat.transpose B).map _).map (k * ·)65  rw [map_map']66  exact map_congr' (fun c => by rw [dot_eq_dotL, dot_eq_dotL, dotL_rsmul]) _6768theorem length_rowMul (r : List (Ext p ν)) (B : Mat p ν) : (rowMul r B).length = (Mat.transpose B).length := length_map _ _6970theorem length_transpose (B : Mat p ν) : (Mat.transpose B).length = (B.headD []).length := by71  show ((List.range _).map _).length = _; rw [length_map, range_eq_rangeS, length_rangeS]7273theorem headD_eq_row (B : Mat p ν) : B.headD [] = row B 0 := by74  cases B with75  | nil => rfl76  | cons r _ => rfl7778/-- A square `B` of size `n > 0` has `n` columns, so `M · B` has rows of length `n`. -/79theorem length_transpose_square (B : Mat p ν) (hB : Square B) (h0 : 0 < B.length) : (Mat.transpose B).length = B.length := by80  rw [length_transpose, headD_eq_row]; exact hB 0 h08182theorem square_mul (M B : Mat p ν) (hB : Square B) (hMB : M.length = B.length) (h0 : 0 < B.length) : Square (M * B) := fun i hi => by83  rw [length_mul] at hi ⊢84  rw [row_mul M B i hi, length_rowMul, length_transpose_square B hB h0, hMB]8586/-! ## Elementary row operations as data -/8788/-- The elementary row operation `addRow i k c`: `c` times row `k` is added to row `i` (`i ≠ k`). -/89inductive Op (p : Nat) [Pos p] (ν : Shell p) where90  | addRow (i k : Nat) (c : Ext p ν)9192def apply1 : Op p ν → Mat p ν → Mat p ν93  | .addRow i k c, M => setAt M i (radd (row M i) (rsmul c (row M k)))9495def apply : List (Op p ν) → Mat p ν → Mat p ν96  | [], M => M97  | o :: os, M => apply os (apply1 o M)9899/-- The operation stays inside an `n × n` matrix. -/100def Valid (n : Nat) : Op p ν → Prop101  | .addRow i k _ => i < n ∧ k < n ∧ i ≠ k102103def ValidL (n : Nat) : List (Op p ν) → Prop104  | [] => True105  | o :: os => Valid n o ∧ ValidL n os106107theorem validL_append (n : Nat) : ∀ (l m : List (Op p ν)), ValidL n l → ValidL n m → ValidL n (l ++ m)108  | [], _, _, hm => hm109  | _ :: l, m, ⟨ho, hl⟩, hm => ⟨ho, validL_append n l m hl hm⟩110111theorem apply_append : ∀ (l m : List (Op p ν)) (M : Mat p ν), apply (l ++ m) M = apply m (apply l M)112  | [], _, _ => rfl113  | o :: l, m, M => apply_append l m (apply1 o M)114115theorem length_apply1 : ∀ (o : Op p ν) (M : Mat p ν), (apply1 o M).length = M.length116  | .addRow i _ _, M => length_setAt M i _117118theorem length_apply : ∀ (l : List (Op p ν)) (M : Mat p ν), (apply l M).length = M.length119  | [], _ => rfl120  | o :: l, M => by show (apply l (apply1 o M)).length = M.length; rw [length_apply l, length_apply1]121122theorem square_apply1 : ∀ (o : Op p ν) (M : Mat p ν), Valid M.length o → Square M → Square (apply1 o M)123  | .addRow i k _, M, ⟨hi, hk, _⟩, hS => square_setAt M hS i _ (by rw [length_radd _ _ (by rw [length_rsmul, hS i hi, hS k hk]), hS i hi])124125theorem square_apply : ∀ (l : List (Op p ν)) (M : Mat p ν), ValidL M.length l → Square M → Square (apply l M)126  | [], _, _, hS => hS127  | o :: l, M, ⟨ho, hl⟩, hS => by128    show Square (apply l (apply1 o M))129    exact square_apply l _ (by rw [length_apply1]; exact hl) (square_apply1 o M ho hS)130131/-! ## The determinant under the operations -/132133theorem det_apply1 : ∀ (o : Op p ν) (M : Mat p ν), Valid M.length o → Square M → det (apply1 o M) = det M134  | .addRow i k c, M, ⟨hi, hk, hik⟩, hS => det_addRow M hS i k hik hi hk c135136theorem det_apply : ∀ (l : List (Op p ν)) (M : Mat p ν), ValidL M.length l → Square M → det (apply l M) = det M137  | [], _, _, _ => rfl138  | o :: l, M, ⟨ho, hl⟩, hS => by139    show det (apply l (apply1 o M)) = det M140    rw [det_apply l _ (by rw [length_apply1]; exact hl) (square_apply1 o M ho hS), det_apply1 o M ho hS]141142/-! ## The operations commute with the right multiplication -/143144theorem addRow_mul (M B : Mat p ν) (hS : Square M) (i k : Nat) (hi : i < M.length) (hk : k < M.length) (c : Ext p ν) :145    setAt M i (radd (row M i) (rsmul c (row M k))) * B = setAt (M * B) i (radd (row (M * B) i) (rsmul c (row (M * B) k))) := by146  show (setAt M i (radd (row M i) (rsmul c (row M k)))).map (rowMul · B) = _147  rw [row_mul M B i hi, row_mul M B k hk, map_setAt, rowMul_radd _ _ _ (by rw [length_rsmul, hS i hi, hS k hk]), rowMul_rsmul]148  rfl149150theorem apply1_mul : ∀ (o : Op p ν) (M B : Mat p ν), Valid M.length o → Square M → apply1 o M * B = apply1 o (M * B)151  | .addRow i k c, M, B, ⟨hi, hk, _⟩, hS => addRow_mul M B hS i k hi hk c152153theorem apply_mul : ∀ (l : List (Op p ν)) (M B : Mat p ν), ValidL M.length l → Square M → apply l M * B = apply l (M * B)154  | [], _, _, _, _ => rfl155  | o :: l, M, B, ⟨ho, hl⟩, hS => by156    show apply l (apply1 o M) * B = apply l (apply1 o (M * B))157    rw [apply_mul l _ B (by rw [length_apply1]; exact hl) (square_apply1 o M ho hS), apply1_mul o M B ho hS]158159160/-! ## Gauss–Jordan elimination, its operations read off the matrix161162No swaps: when the pivot of column `j` lies in a row `k` below `r`, row `r` gains row `k` (one `addRow r k 1`, which makes163`M[r][j] = M[k][j] ≠ 0` since the rows between have a zero there), then every other row loses its multiple of row `r`. -/164165/-- The entry `M[i][c]`. -/166def E (M : Mat p ν) (i c : Nat) : Ext p ν := nth 0 (row M i) c167168/-- The least row `k ≥ i` among the next `f` rows with `M[k][j] ≠ 0`. -/169def pivotCol (M : Mat p ν) (j : Nat) : Nat → Nat → Option Nat170  | 0, _ => none171  | f + 1, i => if E M i j = 0 then pivotCol M j f (i + 1) else some i172173theorem pivotCol_none (M : Mat p ν) (j : Nat) : ∀ (f i : Nat), pivotCol M j f i = none → ∀ k, i ≤ k → k < i + f → E M k j = 0174  | 0, _, _, _, hk, hk' => absurd (Nat.lt_of_le_of_lt hk hk') (Nat.lt_irrefl _)175  | f + 1, i, h, k, hk, hk' => by176    unfold pivotCol at h177    match decEq (E M i j) 0 with178    | .isTrue e =>179      rw [ite_eq_left e] at h180      match Nat.decEq i k with181      | .isTrue e' => subst e'; exact e182      | .isFalse e' => exact pivotCol_none M j f (i + 1) h k (Nat.lt_of_le_of_ne hk e') (Nat.add_right_comm i 1 f ▸ hk')183    | .isFalse e => rw [ite_eq_right e] at h; exact nomatch h184185theorem pivotCol_some (M : Mat p ν) (j : Nat) : ∀ (f i k : Nat), pivotCol M j f i = some k →186    i ≤ k ∧ k < i + f ∧ E M k j ≠ 0 ∧ ∀ i', i ≤ i' → i' < k → E M i' j = 0187  | 0, _, _, h => nomatch h188  | f + 1, i, k, h => by189    unfold pivotCol at h190    match decEq (E M i j) 0 with191    | .isTrue e =>192      rw [ite_eq_left e] at h193      obtain ⟨h1, h2, h3, h4⟩ := pivotCol_some M j f (i + 1) k h194      refine ⟨Nat.le_of_succ_le h1, Nat.add_right_comm i 1 f ▸ h2, h3, fun i' hi' hik => ?_⟩195      match Nat.decEq i i' with196      | .isTrue e' => subst e'; exact e197      | .isFalse e' => exact h4 i' (Nat.lt_of_le_of_ne hi' e') hik198    | .isFalse e =>199      rw [ite_eq_right e] at h200      cases h201      exact ⟨Nat.le_refl i, Nat.lt_add_of_pos_right (Nat.succ_pos f), e, fun i' hi' hik => absurd (Nat.lt_of_le_of_lt hi' hik) (Nat.lt_irrefl _)⟩202203/-- The entries after a row addition: row `i` gains `c` times row `k`, the others stay. -/204theorem nth_radd : ∀ (r s : List (Ext p ν)) (c : Nat), r.length = s.length → nth 0 (radd r s) c = nth 0 r c + nth 0 s c205  | [], [], _, _ => (add_zero 0).symm206  | [], _ :: _, _, h => absurd h (Nat.ne_of_lt (Nat.succ_pos _))207  | _ :: _, [], _, h => absurd h (Nat.ne_of_gt (Nat.succ_pos _))208  | _ :: _, _ :: _, 0, _ => rfl209  | _ :: r, _ :: s, c + 1, h => nth_radd r s c (Nat.succ.inj h)210211theorem E_addRow (M : Mat p ν) (hS : Square M) (i k : Nat) (hi : i < M.length) (hk : k < M.length) (c : Ext p ν) (i' c' : Nat) :212    E (apply1 (.addRow i k c) M) i' c' = if i' = i then E M i c' + c * E M k c' else E M i' c' := by213  show nth 0 (row (setAt M i (radd (row M i) (rsmul c (row M k)))) i') c' = _214  match Nat.decEq i' i with215  | .isTrue e => rw [e, ite_eq_left rfl, row_setAt_eq _ _ _ hi, nth_radd _ _ _ (by rw [length_rsmul, hS i hi, hS k hk]), nth_rsmul]; rfl216  | .isFalse e => rw [ite_eq_right e, row_setAt_ne _ _ _ _ (fun h => e h.symm)]; rfl217218/-- The clearing of column `j` below and above the pivot row `r`, from the matrix `M1` whose pivot is `a` with `a · ainv = 1`:219the operations for the rows `i < m`, `i ≠ r`, each reading its coefficient off `M1`. -/220def clearOps (M1 : Mat p ν) (r j : Nat) (ainv : Ext p ν) : Nat → List (Op p ν)221  | 0 => []222  | i + 1 => clearOps M1 r j ainv i ++ (if i = r then [] else [Op.addRow i r (-(E M1 i j) * ainv)])223224theorem validL_clearOps (M1 : Mat p ν) (r j : Nat) (ainv : Ext p ν) (n : Nat) (hr : r < n) : ∀ m, m ≤ n → ValidL n (clearOps M1 r j ainv m)225  | 0, _ => trivial226  | m + 1, hm => by227    show ValidL n (clearOps M1 r j ainv m ++ _)228    refine validL_append n _ _ (validL_clearOps M1 r j ainv n hr m (Nat.le_of_succ_le hm)) ?_229    match Nat.decEq m r with230    | .isTrue e => rw [ite_eq_left e]; trivial231    | .isFalse e => rw [ite_eq_right e]; exact ⟨⟨Nat.lt_of_succ_le hm, hr, e⟩, trivial⟩232233theorem square_clearOps (M1 : Mat p ν) (hS : Square M1) (r j : Nat) (hr : r < M1.length) (ainv : Ext p ν) (m : Nat) (hm : m ≤ M1.length) :234    Square (apply (clearOps M1 r j ainv m) M1) :=235  square_apply _ _ (validL_clearOps M1 r j ainv M1.length hr m hm) hS236237/-- The rows the clearing leaves alone: row `r` and the rows `i ≥ m`. -/238theorem E_clearOps_fixed (M1 : Mat p ν) (hS : Square M1) (r j : Nat) (hr : r < M1.length) (ainv : Ext p ν) :239    ∀ m, m ≤ M1.length → ∀ i c, (i = r ∨ m ≤ i) → E (apply (clearOps M1 r j ainv m) M1) i c = E M1 i c240  | 0, _, _, _, _ => rfl241  | m + 1, hm, i, c, hi => by242    show E (apply (clearOps M1 r j ainv m ++ _) M1) i c = _243    rw [apply_append]244    have ih := E_clearOps_fixed M1 hS r j hr ainv m (Nat.le_of_succ_le hm)245    have hi' : i = r ∨ m ≤ i := match hi with | Or.inl h => Or.inl h | Or.inr h => Or.inr (Nat.le_of_succ_le h)246    match Nat.decEq m r with247    | .isTrue e => rw [ite_eq_left e]; exact ih i c hi'248    | .isFalse e =>249      rw [ite_eq_right e]250      show E (apply1 (.addRow m r (-(E M1 m j) * ainv)) (apply (clearOps M1 r j ainv m) M1)) i c = _251      rw [E_addRow _ (square_clearOps M1 hS r j hr ainv m (Nat.le_of_succ_le hm)) m r252            (by rw [length_apply]; exact Nat.lt_of_succ_le hm) (by rw [length_apply]; exact hr)]253      have him : i ≠ m := fun h => match hi with254        | Or.inl h' => e (h.symm.trans h')255        | Or.inr h' => Nat.lt_irrefl _ (Nat.lt_of_lt_of_le (Nat.lt_succ_self m) (h ▸ h'))256      rw [ite_eq_right him]; exact ih i c hi'257258/-- The cleared rows `i < m`, `i ≠ r`: `E M1 i c − E M1 i j · ainv · E M1 r c`. -/259theorem E_clearOps_cleared (M1 : Mat p ν) (hS : Square M1) (r j : Nat) (hr : r < M1.length) (ainv : Ext p ν) :260    ∀ m, m ≤ M1.length → ∀ i c, i < m → i ≠ r →261      E (apply (clearOps M1 r j ainv m) M1) i c = E M1 i c + (-(E M1 i j) * ainv) * E M1 r c262  | 0, _, i, _, hi, _ => absurd hi (Nat.not_lt_zero i)263  | m + 1, hm, i, c, hi, hir => by264    show E (apply (clearOps M1 r j ainv m ++ _) M1) i c = _265    rw [apply_append]266    match Nat.decEq m r with267    | .isTrue e =>268      rw [ite_eq_left e]269      exact E_clearOps_cleared M1 hS r j hr ainv m (Nat.le_of_succ_le hm) i c270        (Nat.lt_of_le_of_ne (Nat.le_of_lt_succ hi) (fun h => hir (h.trans e))) hir271    | .isFalse e =>272      rw [ite_eq_right e]273      show E (apply1 (.addRow m r (-(E M1 m j) * ainv)) (apply (clearOps M1 r j ainv m) M1)) i c = _274      rw [E_addRow _ (square_clearOps M1 hS r j hr ainv m (Nat.le_of_succ_le hm)) m r275            (by rw [length_apply]; exact Nat.lt_of_succ_le hm) (by rw [length_apply]; exact hr)]276      match Nat.decEq i m with277      | .isTrue e' =>278        rw [ite_eq_left e', e']279        rw [E_clearOps_fixed M1 hS r j hr ainv m (Nat.le_of_succ_le hm) m c (Or.inr (Nat.le_refl m)),280          E_clearOps_fixed M1 hS r j hr ainv m (Nat.le_of_succ_le hm) r c (Or.inl rfl)]281      | .isFalse e' =>282        rw [ite_eq_right e']283        exact E_clearOps_cleared M1 hS r j hr ainv m (Nat.le_of_succ_le hm) i c (Nat.lt_of_le_of_ne (Nat.le_of_lt_succ hi) e') hir284285286/-- The pivot fix: when the pivot of column `j` among the rows `r, …, n−1` is the row `k ≠ r`, row `r` gains row `k`. -/287def fixOps (n r j : Nat) (M : Mat p ν) : List (Op p ν) :=288  match pivotCol M j (n - r) r with289  | none => []290  | some k => if k = r then [] else [Op.addRow r k 1]291292/-- The operations of column `j`: the fix, then the clearing of the column outside the pivot row. -/293def colOps (n r j : Nat) (M : Mat p ν) : List (Op p ν) :=294  match pivotCol M j (n - r) r with295  | none => []296  | some _ =>297    let M1 := apply (fixOps n r j M) M298    fixOps n r j M ++ clearOps M1 r j (inv (E M1 r j)) n299300/-- The rows `r, …, n−1` are zero in the columns `c < j`. -/301def InvA (n r j : Nat) (M : Mat p ν) : Prop := ∀ i c, r ≤ i → i < n → c < j → E M i c = 0302303/-- The columns `c < j` carry a nonzero diagonal entry and nothing else. -/304def InvC (n j : Nat) (M : Mat p ν) : Prop := ∀ c, c < j → E M c c ≠ 0 ∧ ∀ i, i < n → i ≠ c → E M i c = 0305306theorem add_sub_of_le' {r n : Nat} (h : r ≤ n) : r + (n - r) = n := FRC.Nat.add_sub_of_le h307308/-- One column with a pivot in row `k`: the entries of the matrix after the fix. -/309theorem E_fix (n r j : Nat) (M : Mat p ν) (hM : M.length = n) (hS : Square M) (k : Nat) (hk : pivotCol M j (n - r) r = some k)310    (hkn : k < n) (hrk : r ≤ k) (i c : Nat) :311    E (apply (fixOps n r j M) M) i c = if i = r ∧ k ≠ r then E M r c + E M k c else E M i c := by312  unfold fixOps; rw [hk]313  show E (apply (if k = r then [] else [Op.addRow r k 1]) M) i c = _314  match Nat.decEq k r with315  | .isTrue e =>316    rw [ite_eq_left e]; show E M i c = _317    rw [ite_eq_right (fun h => h.2 e)]318  | .isFalse e =>319    rw [ite_eq_right e]320    show E (apply1 (.addRow r k 1) M) i c = _321    rw [E_addRow M hS r k (hM ▸ Nat.lt_of_le_of_lt hrk hkn) (hM ▸ hkn) 1 i c, Ext.one_mul]322    match Nat.decEq i r with323    | .isTrue e' => rw [ite_eq_left e', ite_eq_left ⟨e', e⟩]324    | .isFalse e' => rw [ite_eq_right e', ite_eq_right (fun h => e' h.1)]325326theorem validL_fixOps (n r j : Nat) (M : Mat p ν) (hr : r ≤ n) : ValidL n (fixOps n r j M) := by327  unfold fixOps328  match h : pivotCol M j (n - r) r with329  | none => exact trivial330  | some k =>331    obtain ⟨h1, h2, _, _⟩ := pivotCol_some M j _ r k h332    rw [add_sub_of_le' hr] at h2333    show ValidL n (if k = r then [] else [Op.addRow r k 1])334    match Nat.decEq k r with335    | .isTrue e => rw [ite_eq_left e]; exact trivial336    | .isFalse e => rw [ite_eq_right e]; exact ⟨⟨Nat.lt_of_le_of_lt h1 h2, h2, fun h' => e h'.symm⟩, trivial⟩337338/-- The column step. From a square `n × n` matrix with `InvA n r j` (and `InvC n j` when `r = j`), the column's339operations are valid, the new rank `r'` is `r + 1` or `r`, and the invariants advance to column `j + 1`. -/340theorem colOps_spec {κ : Nat} {g : Shell p} (F : Frame p κ g) (hν : ¬IsSquare ν) (n r j : Nat) (M : Mat p ν)341    (hM : M.length = n) (hS : Square M) (hr : r ≤ n) (hA : InvA n r j M) (hC : r = j → InvC n j M) :342    ValidL n (colOps n r j M) ∧343    ((pivotCol M j (n - r) r = none ∧ InvA n r (j + 1) (apply (colOps n r j M) M)) ∨344     (∃ k, pivotCol M j (n - r) r = some k ∧ r < n ∧ InvA n (r + 1) (j + 1) (apply (colOps n r j M) M) ∧345       (r = j → InvC n (j + 1) (apply (colOps n r j M) M)))) := by346  unfold colOps347  match hk : pivotCol M j (n - r) r with348  | none =>349    refine ⟨trivial, Or.inl ⟨rfl, fun i c hi hin hc => ?_⟩⟩350    show E M i c = 0351    match Nat.lt_or_ge c j with352    | Or.inl hcj => exact hA i c hi hin hcj353    | Or.inr hcj =>354      have e : c = j := Nat.le_antisymm (Nat.le_of_lt_succ hc) hcj355      rw [e]; exact pivotCol_none M j _ r hk i hi (by rw [add_sub_of_le' hr]; exact hin)356  | some k =>357    obtain ⟨hrk, hkn, hkj, hbelow⟩ := pivotCol_some M j _ r k hk358    rw [add_sub_of_le' hr] at hkn359    have hrn : r < n := Nat.lt_of_le_of_lt hrk hkn360    -- the fixed matrix and its pivot361    have hM1len : (apply (fixOps n r j M) M).length = n := by rw [length_apply, hM]362    have hM1sq : Square (apply (fixOps n r j M) M) := square_apply _ _ (hM ▸ validL_fixOps n r j M hr) hS363    have hE1 := E_fix n r j M hM hS k hk hkn hrk364    have ha : E (apply (fixOps n r j M) M) r j ≠ 0 := by365      rw [hE1]366      match Nat.decEq k r with367      | .isTrue e => rw [ite_eq_right (fun h => h.2 e), ← e]; exact hkj368      | .isFalse e => rw [ite_eq_left ⟨rfl, e⟩, hbelow r (Nat.le_refl r) (Nat.lt_of_le_of_ne hrk (fun h => e h.symm)), zero_add]; exact hkj369    -- the entries after the clearing370    have hfix : ∀ i c, E (apply (fixOps n r j M) M) i c = if i = r ∧ k ≠ r then E M r c + E M k c else E M i c := hE1371    have hrn1 : r < (apply (fixOps n r j M) M).length := by rw [hM1len]; exact hrn372    have hcl_fixed := E_clearOps_fixed _ hM1sq r j hrn1 (inv (E (apply (fixOps n r j M) M) r j)) n (Nat.le_of_eq hM1len.symm)373    have hcl := E_clearOps_cleared _ hM1sq r j hrn1 (inv (E (apply (fixOps n r j M) M) r j)) n (Nat.le_of_eq hM1len.symm)374    rw [apply_append]375    -- row r of the fixed matrix is zero in the processed columns376    have hr0 : ∀ c, c < j → E (apply (fixOps n r j M) M) r c = 0 := fun c hc => by377      rw [hfix]378      match Nat.decEq k r with379      | .isTrue e => rw [ite_eq_right (fun h => h.2 e)]; exact hA r c (Nat.le_refl r) hrn hc380      | .isFalse e => rw [ite_eq_left ⟨rfl, e⟩, hA r c (Nat.le_refl r) hrn hc, hA k c hrk hkn hc, add_zero]381    have hinv : E (apply (fixOps n r j M) M) r j * inv (E (apply (fixOps n r j M) M) r j) = 1 := mul_inv F hν ha382    -- every row other than r is zero in column j after the clearing383    have hcolj : ∀ i, i < n → i ≠ r → E (apply (clearOps _ r j _ n) (apply (fixOps n r j M) M)) i j = 0 := fun i hin hir => by384      rw [hcl i j hin hir, ← neg_mul, ← neg_mul, Ext.mul_assoc, Ext.mul_comm (inv _), hinv, Ext.mul_one, add_neg]385    refine ⟨validL_append n _ _ (validL_fixOps n r j M hr) (validL_clearOps _ r j _ n hrn n (Nat.le_refl n)),386      Or.inr ⟨k, rfl, hrn, fun i c hi hin hc => ?_, fun hrj c hc => ?_⟩⟩387    · -- InvA n (r + 1) (j + 1)388      have hir : i ≠ r := fun h => Nat.lt_irrefl r (h ▸ hi)389      match Nat.lt_or_ge c j with390      | Or.inl hcj =>391        rw [hcl i c hin hir, hr0 c hcj, mul_zero, add_zero, hfix, ite_eq_right (fun h => hir h.1)]392        exact hA i c (Nat.le_of_succ_le hi) hin hcj393      | Or.inr hcj =>394        have e : c = j := Nat.le_antisymm (Nat.le_of_lt_succ hc) hcj395        rw [e]; exact hcolj i hin hir396    · -- InvC n (j + 1) when r = j397      subst hrj398      match Nat.lt_or_ge c r with399      | Or.inl hcr =>400        obtain ⟨hdiag, hoff⟩ := hC rfl c hcr401        have hcr' : c ≠ r := Nat.ne_of_lt hcr402        have hr0c : E (apply (fixOps n r r M) M) r c = 0 := hr0 c hcr403        have hfixc : ∀ i, i ≠ r → E (apply (fixOps n r r M) M) i c = E M i c := fun i hir => by404          rw [hfix, ite_eq_right (fun h => hir h.1)]405        refine ⟨?_, fun i hin hic => ?_⟩406        · rw [hcl c c (Nat.lt_trans hcr hrn) hcr', hr0c, mul_zero, add_zero, hfixc c hcr']; exact hdiag407        · match Nat.decEq i r with408          | .isTrue e => rw [e, hcl_fixed r c (Or.inl rfl)]; exact hr0c409          | .isFalse e => rw [hcl i c hin e, hr0c, mul_zero, add_zero, hfixc i e]; exact hoff i hin hic410      | Or.inr hcr =>411        have e : c = r := Nat.le_antisymm (Nat.le_of_lt_succ hc) hcr412        subst e413        exact ⟨by rw [hcl_fixed c c (Or.inl rfl)]; exact ha, fun i hin hic => hcolj i hin hic⟩414415/-- The rank after column `j`: one more when the column has a pivot. -/416def rankStep (n r j : Nat) (M : Mat p ν) : Nat :=417  match pivotCol M j (n - r) r with418  | none => r419  | some _ => r + 1420421theorem rankStep_none (n r j : Nat) (M : Mat p ν) (h : pivotCol M j (n - r) r = none) : rankStep n r j M = r := by422  unfold rankStep; rw [h]423424theorem rankStep_some (n r j : Nat) (M : Mat p ν) (k : Nat) (h : pivotCol M j (n - r) r = some k) : rankStep n r j M = r + 1 := by425  unfold rankStep; rw [h]426427/-- Gauss–Jordan over `f` more columns from column `j` with rank `r`: the operations and the final rank. -/428def gj (n : Nat) : Nat → Nat → Nat → Mat p ν → List (Op p ν) × Nat429  | 0, r, _, _ => ([], r)430  | f + 1, r, j, M =>431    (colOps n r j M ++ (gj n f (rankStep n r j M) (j + 1) (apply (colOps n r j M) M)).1,432     (gj n f (rankStep n r j M) (j + 1) (apply (colOps n r j M) M)).2)433434theorem gj_spec {κ : Nat} {g : Shell p} (F : Frame p κ g) (hν : ¬IsSquare ν) (n : Nat) :435    ∀ (f r j : Nat) (M : Mat p ν), M.length = n → Square M → r ≤ j → j + f = n → InvA n r j M → (r = j → InvC n j M) →436      ValidL n (gj n f r j M).1 ∧ (gj n f r j M).2 ≤ n ∧437      ((gj n f r j M).2 = n → InvC n n (apply (gj n f r j M).1 M)) ∧438      ((gj n f r j M).2 < n → InvA n (gj n f r j M).2 n (apply (gj n f r j M).1 M))439  | 0, r, j, M, _, _, hrj, hj, hA, hC => by440    rw [Nat.add_zero] at hj; subst hj441    exact ⟨trivial, hrj, fun h => hC h, fun _ => hA⟩442  | f + 1, r, j, M, hM, hS, hrj, hj, hA, hC => by443    have hr : r ≤ n := Nat.le_trans hrj (hj ▸ Nat.le_add_right j (f + 1))444    obtain ⟨hv, hcase⟩ := colOps_spec F hν n r j M hM hS hr hA hC445    have hj' : j + 1 + f = n := by rw [Nat.add_right_comm]; exact hj446    show ValidL n (colOps n r j M ++ (gj n f (rankStep n r j M) (j + 1) (apply (colOps n r j M) M)).1) ∧447      (gj n f (rankStep n r j M) (j + 1) (apply (colOps n r j M) M)).2 ≤ n ∧448      ((gj n f (rankStep n r j M) (j + 1) (apply (colOps n r j M) M)).2 = n →449        InvC n n (apply (colOps n r j M ++ (gj n f (rankStep n r j M) (j + 1) (apply (colOps n r j M) M)).1) M)) ∧450      ((gj n f (rankStep n r j M) (j + 1) (apply (colOps n r j M) M)).2 < n →451        InvA n (gj n f (rankStep n r j M) (j + 1) (apply (colOps n r j M) M)).2 n452          (apply (colOps n r j M ++ (gj n f (rankStep n r j M) (j + 1) (apply (colOps n r j M) M)).1) M))453    have hM' : (apply (colOps n r j M) M).length = n := by rw [length_apply, hM]454    have hS' : Square (apply (colOps n r j M) M) := square_apply _ _ (hM ▸ hv) hS455    match hcase with456    | Or.inl ⟨hnone, hA'⟩ =>457      rw [rankStep_none n r j M hnone]458      obtain ⟨v2, hle, hC2, hA2⟩ := gj_spec F hν n f r (j + 1) _ hM' hS' (Nat.le_succ_of_le hrj) hj' hA'459        (fun h => absurd h (Nat.ne_of_lt (Nat.lt_succ_of_le hrj)))460      refine ⟨validL_append n _ _ hv v2, hle, ?_, ?_⟩461      · intro h; rw [apply_append]; exact hC2 h462      · intro h; rw [apply_append]; exact hA2 h463    | Or.inr ⟨k, hsome, _, hA', hC'⟩ =>464      rw [rankStep_some n r j M k hsome]465      obtain ⟨v2, hle, hC2, hA2⟩ := gj_spec F hν n f (r + 1) (j + 1) _ hM' hS' (Nat.succ_le_succ hrj) hj' hA'466        (fun h => hC' (Nat.succ.inj h))467      refine ⟨validL_append n _ _ hv v2, hle, ?_, ?_⟩468      · intro h; rw [apply_append]; exact hC2 h469      · intro h; rw [apply_append]; exact hA2 h470471472/-! ## From the end state to the theorem -/473474theorem nth_unitRow (n i c : Nat) (hc : c < n) : nth 0 (unitRow (ν := ν) n i) c = if i = c then 1 else 0 := by475  show nth 0 ((rangeS n).map fun j => if i = j then (1 : Ext p ν) else 0) c = _476  rw [nth_map_lt (fun j => if i = j then (1 : Ext p ν) else 0) 0 0 _ c (by rw [length_rangeS]; exact hc), nth_rangeS n c hc]477478theorem rows_ext (M N : Mat p ν) (hl : M.length = N.length) (h : ∀ i, i < M.length → row M i = row N i) : M = N :=479  ext_nth [] M N hl h480481/-- A diagonal matrix (`InvC n n`): its `i`-th row is the `i`-th diagonal entry times the unit row. -/482theorem row_of_invC (D : Mat p ν) (hD : D.length = n) (hS : Square D) (hC : InvC n n D) (i : Nat) (hi : i < n) :483    row D i = rsmul (E D i i) (unitRow n i) := by484  refine ext_nth 0 _ _ (by rw [length_rsmul, length_unitRow, hS i (hD ▸ hi), hD]) fun c hc => ?_485  rw [hS i (hD ▸ hi), hD] at hc486  rw [nth_rsmul, nth_unitRow n i c hc]487  match Nat.decEq i c with488  | .isTrue e => rw [ite_eq_left e, Ext.mul_one, e]; rfl489  | .isFalse e => rw [ite_eq_right e, mul_zero]; exact (hC c hc).2 i hi e490491/-- The rows `i < k` of `N` scaled by `d i`. -/492def scaleRows (d : Nat → Ext p ν) (N : Mat p ν) : Nat → Mat p ν493  | 0 => N494  | k + 1 => setAt (scaleRows d N k) k (rsmul (d k) (row N k))495496def prodUpto (d : Nat → Ext p ν) : Nat → Ext p ν497  | 0 => 1498  | k + 1 => prodUpto d k * d k499500theorem length_scaleRows (d : Nat → Ext p ν) (N : Mat p ν) : ∀ k, (scaleRows d N k).length = N.length501  | 0 => rfl502  | k + 1 => by show (setAt (scaleRows d N k) k _).length = _; rw [length_setAt, length_scaleRows d N k]503504theorem row_scaleRows_ge (d : Nat → Ext p ν) (N : Mat p ν) : ∀ k i, k ≤ i → row (scaleRows d N k) i = row N i505  | 0, _, _ => rfl506  | k + 1, i, h => by507    show row (setAt (scaleRows d N k) k _) i = _508    rw [row_setAt_ne _ _ _ _ (Nat.ne_of_lt h), row_scaleRows_ge d N k i (Nat.le_of_succ_le h)]509510theorem row_scaleRows_lt (d : Nat → Ext p ν) (N : Mat p ν) : ∀ k i, i < k → i < N.length → row (scaleRows d N k) i = rsmul (d i) (row N i)511  | 0, i, h, _ => absurd h (Nat.not_lt_zero i)512  | k + 1, i, h, hi => by513    show row (setAt (scaleRows d N k) k _) i = _514    match Nat.decEq i k with515    | .isTrue e => rw [e, row_setAt_eq _ _ _ (by rw [length_scaleRows]; exact e ▸ hi)]516    | .isFalse e => rw [row_setAt_ne _ _ _ _ (fun h' => e h'.symm), row_scaleRows_lt d N k i (Nat.lt_of_le_of_ne (Nat.le_of_lt_succ h) e) hi]517518theorem det_scaleRows (d : Nat → Ext p ν) (N : Mat p ν) : ∀ k, k ≤ N.length → det (scaleRows d N k) = prodUpto d k * det N519  | 0, _ => (Ext.one_mul _).symm520  | k + 1, hk => by521    show det (setAt (scaleRows d N k) k (rsmul (d k) (row N k))) = prodUpto d k * d k * det N522    have hlen : (scaleRows d N k).length = N.length - 1 + 1 := by523      rw [length_scaleRows]; exact length_form_one (Nat.lt_of_lt_of_le (Nat.succ_pos k) hk)524    have hk' : k < N.length - 1 + 1 := by rw [← length_form_one (Nat.lt_of_lt_of_le (Nat.succ_pos k) hk)]; exact hk525    rw [det_setAt_smul (N.length - 1) _ k (d k) (row N k) hlen hk', ← row_scaleRows_ge d N k k (Nat.le_refl k),526      setAt_row_self _ _ (by rw [length_scaleRows]; exact hk), det_scaleRows d N k (Nat.le_of_succ_le hk),527      ← Ext.mul_assoc, Ext.mul_comm (d k)]528529theorem map_nth_rangeS : ∀ r : List (Ext p ν), (rangeS r.length).map (fun j => nth 0 r j) = r530  | [] => rfl531  | x :: r => by532    show (rangeS (r.length + 1)).map (fun j => nth 0 (x :: r) j) = x :: r533    rw [rangeS_succ_eq_map, List.map, map_map']534    exact congrArg (x :: ·) (map_nth_rangeS r)535536/-- The unit row times `B` is the corresponding row of `B`. -/537theorem rowMul_unitRow (B : Mat p ν) (hB : Square B) (hBn : B.length = n) (h0 : 0 < n) (i : Nat) (hi : i < n) :538    rowMul (unitRow n i) B = row B i := by539  show ((List.range (B.headD []).length).map (Mat.col B)).map (fun c => Mat.dot (unitRow n i) c) = row B i540  rw [headD_eq_row, hB 0 (hBn ▸ h0), hBn, range_eq_rangeS, map_map']541  have hr : (row B i).length = n := by rw [hB i (hBn ▸ hi), hBn]542  have e : ∀ j, j < n → Mat.dot (unitRow n i) (Mat.col B j) = nth 0 (row B i) j := fun j _ => by543    rw [dot_eq_dotL, dotL_unitRow n i _ hi (by show n ≤ (B.map _).length; rw [length_map, hBn]; exact Nat.le_refl n)]544    exact nth_map_lt (fun r => nth 0 r j) [] 0 B i (hBn ▸ hi)545  calc (rangeS n).map (fun j => Mat.dot (unitRow n i) (Mat.col B j))546      = (rangeS n).map (fun j => nth 0 (row B i) j) := map_congr_rangeS n _ _ e547    _ = row B i := by rw [← hr]; exact map_nth_rangeS (row B i)548where549  map_congr_rangeS : ∀ (n : Nat) (f g : Nat → Ext p ν), (∀ j, j < n → f j = g j) → (rangeS n).map f = (rangeS n).map g550    | 0, _, _, _ => rfl551    | n + 1, f, g, h => by552      show (rangeS n ++ [n]).map f = (rangeS n ++ [n]).map g553      rw [map_append_single, map_append_single, map_congr_rangeS n f g (fun j hj => h j (Nat.lt_succ_of_lt hj)), h n (Nat.lt_succ_self n)]554555theorem nth_beyond {α : Type} (d : α) : ∀ (l : List α) (i : Nat), l.length ≤ i → nth d l i = d556  | [], _, _ => rfl557  | _ :: _, 0, h => absurd h (Nat.not_succ_le_zero _)558  | _ :: l, i + 1, h => nth_beyond d l i (Nat.le_of_succ_le_succ h)559560/-- A matrix with a zero row has determinant zero. -/561theorem det_zero_row (M : Mat p ν) (i : Nat) (hi : i < M.length) (h : ∀ c, nth 0 (row M i) c = 0) : det M = 0 := by562  have e : row M i = rsmul 0 (row M i) :=563    ext_nth 0 _ _ (by rw [length_rsmul]) fun c _ => by rw [nth_rsmul, zero_mul, h c]564  have hlen : M.length = M.length - 1 + 1 := length_form_one (Nat.lt_of_le_of_lt (Nat.zero_le i) hi)565  calc det M = det (setAt M i (rsmul 0 (row M i))) := by rw [← e, setAt_row_self M i hi]566    _ = 0 * det (setAt M i (row M i)) := det_setAt_smul (M.length - 1) M i 0 (row M i) hlen (by rw [← hlen]; exact hi)567    _ = 0 := zero_mul _568569theorem nth_rowMul_zero (r : List (Ext p ν)) (B : Mat p ν) (h : ∀ c, nth 0 r c = 0) (c : Nat) : nth 0 (rowMul r B) c = 0 := by570  have e : r = rsmul 0 r := ext_nth 0 _ _ (by rw [length_rsmul]) fun c _ => by rw [nth_rsmul, zero_mul, h c]571  rw [e, rowMul_rsmul, nth_rsmul, zero_mul]572573/-- `det (A · B) = det A · det B` for square matrices of the same size over `K`. -/574theorem det_mul {κ : Nat} {g : Shell p} (F : Frame p κ g) (hν : ¬IsSquare ν) (A B : Mat p ν) (hA : Square A) (hB : Square B)575    (hAB : A.length = B.length) : det (A * B) = det A * det B := by576  match hn : A.length with577  | 0 =>578    have hA0 : A = [] := match A, hn with | [], _ => rfl579    have hB0 : B = [] := match B, hAB.symm.trans hn with | [], _ => rfl580    rw [hA0, hB0]; show (1 : Ext p ν) = 1 * 1; exact (Ext.one_mul 1).symm581  | n + 1 =>582    have hn' : A.length = n + 1 := hn583    have hBn : B.length = n + 1 := hAB.symm.trans hn584    have h0 : 0 < n + 1 := Nat.succ_pos n585    obtain ⟨hv, hle, hdiag, hzero⟩ := gj_spec F hν (n + 1) (n + 1) 0 0 A hn' hA (Nat.le_refl 0) (Nat.zero_add _)586      (fun _ c _ _ hc => absurd hc (Nat.not_lt_zero c)) (fun _ c hc => absurd hc (Nat.not_lt_zero c))587    -- the end state D and the transported product588    have hDlen : (apply (gj (n + 1) (n + 1) 0 0 A).1 A).length = n + 1 := by rw [length_apply, hn']589    have hDsq : Square (apply (gj (n + 1) (n + 1) 0 0 A).1 A) := square_apply _ _ (by rw [hn']; exact hv) hA590    have hdetD : det (apply (gj (n + 1) (n + 1) 0 0 A).1 A) = det A := det_apply _ _ (by rw [hn']; exact hv) hA591    have hABsq : Square (A * B) := square_mul A B hB hAB (by rw [hBn]; exact h0)592    have hABlen : (A * B).length = n + 1 := by rw [length_mul, hn']593    have hdetAB : det (apply (gj (n + 1) (n + 1) 0 0 A).1 A * B) = det (A * B) := by594      rw [apply_mul _ _ _ (by rw [hn']; exact hv) hA]; exact det_apply _ _ (by rw [hABlen]; exact hv) hABsq595    match Nat.lt_or_ge (gj (n + 1) (n + 1) 0 0 A).2 (n + 1) with596    | Or.inr hge =>597      have hC := hdiag (Nat.le_antisymm hle hge)598      -- D = scaleRows d (ident) and D · B = scaleRows d B, with d the diagonal599      have hrow := row_of_invC _ hDlen hDsq hC600      have eD : apply (gj (n + 1) (n + 1) 0 0 A).1 A = scaleRows (fun i => E (apply (gj (n + 1) (n + 1) 0 0 A).1 A) i i) (Mat.ident (n + 1)) (n + 1) := by601        refine rows_ext _ _ (by rw [length_scaleRows, length_ident, hDlen]) fun i hi => ?_602        rw [hDlen] at hi603        rw [hrow i hi, row_scaleRows_lt _ _ _ i hi (by rw [length_ident]; exact hi), row_ident _ _ hi]604      have eDB : apply (gj (n + 1) (n + 1) 0 0 A).1 A * B = scaleRows (fun i => E (apply (gj (n + 1) (n + 1) 0 0 A).1 A) i i) B (n + 1) := by605        refine rows_ext _ _ (by rw [length_scaleRows, length_mul, hDlen, hBn]) fun i hi => ?_606        rw [length_mul, hDlen] at hi607        rw [row_mul _ _ i (by rw [hDlen]; exact hi), hrow i hi, rowMul_rsmul, rowMul_unitRow B hB hBn h0 i hi,608          row_scaleRows_lt _ _ _ i hi (by rw [hBn]; exact hi)]609      have h1 : det (apply (gj (n + 1) (n + 1) 0 0 A).1 A * B) =610          prodUpto (fun i => E (apply (gj (n + 1) (n + 1) 0 0 A).1 A) i i) (n + 1) * det B := by611        rw [eDB]; exact det_scaleRows _ _ _ (Nat.le_of_eq hBn.symm)612      have h2 : det (apply (gj (n + 1) (n + 1) 0 0 A).1 A) = prodUpto (fun i => E (apply (gj (n + 1) (n + 1) 0 0 A).1 A) i i) (n + 1) := by613        have := congrArg det eD614        rw [det_scaleRows _ _ _ (Nat.le_of_eq (length_ident (n + 1)).symm), det_ident, Ext.mul_one] at this615        exact this616      calc det (A * B) = det (apply (gj (n + 1) (n + 1) 0 0 A).1 A * B) := hdetAB.symm617        _ = prodUpto (fun i => E (apply (gj (n + 1) (n + 1) 0 0 A).1 A) i i) (n + 1) * det B := h1618        _ = det (apply (gj (n + 1) (n + 1) 0 0 A).1 A) * det B := by rw [h2]619        _ = det A * det B := by rw [hdetD]620    | Or.inl hlt =>621      have hAz := hzero hlt622      -- row n of D is zero: D and D · B are singular623      have hrow0 : ∀ c, nth 0 (row (apply (gj (n + 1) (n + 1) 0 0 A).1 A) n) c = 0 := fun c => by624        match Nat.lt_or_ge c (n + 1) with625        | Or.inl hc => exact hAz n c (Nat.le_of_lt_succ hlt) (Nat.lt_succ_self n) hc626        | Or.inr hc => exact nth_beyond _ _ _ (by rw [hDsq n (by rw [hDlen]; exact Nat.lt_succ_self n), hDlen]; exact hc)627      have hD0 : det (apply (gj (n + 1) (n + 1) 0 0 A).1 A) = 0 := det_zero_row _ n (by rw [hDlen]; exact Nat.lt_succ_self n) hrow0628      have hDB0 : det (apply (gj (n + 1) (n + 1) 0 0 A).1 A * B) = 0 :=629        det_zero_row _ n (by rw [length_mul, hDlen]; exact Nat.lt_succ_self n)630          (fun c => by rw [row_mul _ _ n (by rw [hDlen]; exact Nat.lt_succ_self n)]; exact nth_rowMul_zero _ B hrow0 c)631      rw [← hdetAB, hDB0, ← hdetD, hD0, zero_mul]632633end FRC.Dirac.Det634