-- Matrices as set-theoretic functions on pairs of natural indices.
-- Part of Cumulus, the Grace17 mathematics library.
-- Copyright (c) 2026 Grant Passmore.
-- Released under the Apache License 2.0. See LICENSE in the repository root.
-- Authors: Grant Passmore

module Grace.Cumulus.LinearAlgebra.Matrix
import FiniteDimensional

-- Dimensions are cutoffs on a total array, so one matrix can be reused at
-- several finite sizes without changing its underlying set-theoretic graph.
predicate RealMatrix (M : Set) := FunctionBetween(M, ℕ × ℕ, ℝ)

-- Coordinates are uniquely determined on the Cartesian product.
-- Extracting them needs only ZF and supports matrix graphs over index pairs.
defunique matrixFirstIndex (p : Set) (i ∈ ℕ) :=
  ∃ j ∈ ℕ, p = (i, j)
requires pair : p ∈ ℕ × ℕ
proof
  have shape := product.membership.forward at ℕ, ℕ, p with pair
  rcases shape with ⟨bound, i, naturalI, j, naturalJ, equation⟩
  witness i
  refine ⟨naturalI, j, naturalJ, equation⟩

uniqueness
  intros a aMember first c cMember second
  obtain ⟨b, bMember, firstPair⟩ := first
  obtain ⟨d, dMember, secondPair⟩ := second
  have same : (a, b) = (c, d)
    cc
  have coordinates := orderedPair.injective at a, b, c, d with same
  tauto

defunique matrixSecondIndex (p : Set) (j ∈ ℕ) :=
  ∃ i ∈ ℕ, p = (i, j)
requires pair : p ∈ ℕ × ℕ
proof
  have shape := product.membership.forward at ℕ, ℕ, p with pair
  rcases shape with ⟨bound, i, naturalI, j, naturalJ, equation⟩
  witness j
  refine ⟨naturalJ, i, naturalI, equation⟩

uniqueness
  intros b bMember first d dMember second
  obtain ⟨a, aMember, firstPair⟩ := first
  obtain ⟨c, cMember, secondPair⟩ := second
  have same : (a, b) = (c, d)
    cc
  have coordinates := orderedPair.injective at a, b, c, d with same
  tauto

theorem matrixFirstIndex.pair (i j : Nat)
  shows matrixFirstIndex((i, j)) = i
proof
  have naturalI : i ∈ ℕ by simp [Nat]
  have naturalJ : j ∈ ℕ by simp [Nat]
  have pair := product.pair with naturalI, naturalJ
  obtain ⟨natural, k, naturalK, shape⟩ := matrixFirstIndex.selected with pair
  obtain ⟨first, second⟩ := orderedPair.injective with shape
  rw ←first
  rfl

theorem matrixSecondIndex.pair (i j : Nat)
  shows matrixSecondIndex((i, j)) = j
proof
  have naturalI : i ∈ ℕ by simp [Nat]
  have naturalJ : j ∈ ℕ by simp [Nat]
  have pair := product.pair with naturalI, naturalJ
  obtain ⟨natural, k, naturalK, shape⟩ := matrixSecondIndex.selected with pair
  obtain ⟨first, second⟩ := orderedPair.injective with shape
  rw ←second
  rfl

theorem RealMatrix.entry (M : Set, i j : Nat)
  requires matrix : RealMatrix(M)
  shows Real(M[i, j])
proof
  unfold RealMatrix at matrix
  have iNatural : i ∈ ℕ by simp [Nat]
  have jNatural : j ∈ ℕ by simp [Nat]
  have pair := product.pair with iNatural, jNatural
  exact FunctionBetween.apply_closed with matrix, pair

defgraph matrixTranspose (M : Set) (p ∈ ℕ × ℕ) : ℝ :=
  M[matrixSecondIndex(p), matrixFirstIndex(p)]
requires matrix : RealMatrix(M)
proof
  have first := matrixFirstIndex.selected with member
  have second := matrixSecondIndex.selected with member
  rcases first with ⟨naturalFirst, j, naturalJ, shapeFirst⟩
  rcases second with ⟨naturalSecond, i, naturalI, shapeSecond⟩
  have typedFirst := nat.omega.backward with naturalFirst
  have typedSecond := nat.omega.backward with naturalSecond
  exact RealMatrix.entry at
    M, matrixSecondIndex(p), matrixFirstIndex(p) with matrix

theorem matrixTranspose.matrix (M : Set)
  requires matrix : RealMatrix(M)
  shows RealMatrix(matrixTranspose(M))
proof
  unfold RealMatrix
  exact matrixTranspose.function with matrix

theorem matrixTranspose.entry (M : Set, i j : Nat)
  requires matrix : RealMatrix(M)
  shows matrixTranspose(M)[i, j] = M[j, i]
proof
  have naturalI : i ∈ ℕ by simp [Nat]
  have naturalJ : j ∈ ℕ by simp [Nat]
  have pair := product.pair with naturalI, naturalJ
  have value := matrixTranspose.value with matrix
  rw [value, matrixFirstIndex.pair, matrixSecondIndex.pair]
  rfl

defgraph matrixAdd (A B : Set) (p ∈ ℕ × ℕ) : ℝ :=
  call2(realAdd, A[matrixFirstIndex(p), matrixSecondIndex(p)],
    B[matrixFirstIndex(p), matrixSecondIndex(p)])
requires matrixA : RealMatrix(A)
requires matrixB : RealMatrix(B)
proof
  have first := matrixFirstIndex.selected with member
  have second := matrixSecondIndex.selected with member
  rcases first with ⟨naturalFirst, j, naturalJ, shapeFirst⟩
  rcases second with ⟨naturalSecond, i, naturalI, shapeSecond⟩
  have typedFirst := nat.omega.backward with naturalFirst
  have typedSecond := nat.omega.backward with naturalSecond
  have left := RealMatrix.entry at A, matrixFirstIndex(p), matrixSecondIndex(p) with matrixA
  have right := RealMatrix.entry at B, matrixFirstIndex(p), matrixSecondIndex(p) with matrixB
  exact realAddValue.real with left, right

theorem matrixAdd.matrix (A B : Set)
  requires matrixA : RealMatrix(A)
  requires matrixB : RealMatrix(B)
  shows RealMatrix(matrixAdd(A, B))
proof
  unfold RealMatrix
  exact matrixAdd.function with matrixA, matrixB

theorem matrixAdd.entry (A B : Set, i j : Nat)
  requires matrixA : RealMatrix(A)
  requires matrixB : RealMatrix(B)
  shows matrixAdd(A, B)[i, j] = (A[i, j]:R) + (B[i, j]:R)
proof
  have naturalI : i ∈ ℕ by simp [Nat]
  have naturalJ : j ∈ ℕ by simp [Nat]
  have pair := product.pair with naturalI, naturalJ
  have value := matrixAdd.value with matrixA, matrixB
  rw [value, matrixFirstIndex.pair, matrixSecondIndex.pair]
  rfl

defgraph matrixScale (C : Real, M : Set) (p ∈ ℕ × ℕ) : ℝ :=
  realMul(C, M[matrixFirstIndex(p), matrixSecondIndex(p)])
requires matrix : RealMatrix(M)
proof
  have first := matrixFirstIndex.selected with member
  have second := matrixSecondIndex.selected with member
  rcases first with ⟨naturalFirst, j, naturalJ, shapeFirst⟩
  rcases second with ⟨naturalSecond, i, naturalI, shapeSecond⟩
  have typedFirst := nat.omega.backward with naturalFirst
  have typedSecond := nat.omega.backward with naturalSecond
  have entry := RealMatrix.entry at M, matrixFirstIndex(p), matrixSecondIndex(p) with matrix
  have scaled := realMulValue.real at C,
    M[matrixFirstIndex(p), matrixSecondIndex(p)] with assumption, entry
  exact scaled

theorem matrixScale.matrix (C : Real, M : Set)
  requires matrix : RealMatrix(M)
  shows RealMatrix(matrixScale(C, M))
proof
  unfold RealMatrix
  exact matrixScale.function with matrix

theorem matrixScale.entry (C : Real, M : Set, i j : Nat)
  requires matrix : RealMatrix(M)
  shows matrixScale(C, M)[i, j] = C * (M[i, j]:R)
proof
  have naturalI : i ∈ ℕ by simp [Nat]
  have naturalJ : j ∈ ℕ by simp [Nat]
  have pair := product.pair with naturalI, naturalJ
  have value := matrixScale.value with matrix
  rw [value, matrixFirstIndex.pair, matrixSecondIndex.pair]
  rfl

-- Equality of indices uniquely determines the diagonal coefficient in ZF.
defunique diagonalCoefficient (i j : Nat) (X ∈ ℝ) :=
  (i = j → X = (1:R)) ∧ (i ≠ j → X = (0:R))
proof
  have alternatives : i = j ∨ i ≠ j by tauto
  cases alternatives with equal | distinct
  . witness (1:R)
    constructor
    . exact realOne.real
    . constructor
      . intro same
        rfl
      . intro impossible
        exfalso
        exact impossible with equal
  . witness (0:R)
    constructor
    . exact realZero.real
    . constructor
      . intro impossible
        exfalso
        exact distinct with impossible
      . intro different
        rfl

uniqueness
  intros u uMember first v vMember second
  have alternatives : i = j ∨ ¬(i = j) by tauto
  cases alternatives with present | absent
  . have uValue : u = (1:R)
      tauto
    have vValue : v = (1:R) by tauto
    cc
  . have uValue : u = (0:R)
      tauto
    have vValue : v = (0:R) by tauto
    cc

theorem diagonalCoefficient.same (i : Nat)
  shows diagonalCoefficient(i, i) = (1:R)
proof
  have selected := diagonalCoefficient.selected at i, i with *
  rcases selected with ⟨real, diagonal, offDiagonal⟩
  have self : i = i by rfl
  exact diagonal with self

theorem diagonalCoefficient.distinct (i j : Nat)
  requires different : i ≠ j
  shows diagonalCoefficient(i, j) = (0:R)
proof
  have selected := diagonalCoefficient.selected at i, j with *
  rcases selected with ⟨real, diagonal, offDiagonal⟩
  exact offDiagonal with different

defgraph identityMatrix (n : Nat) (p ∈ ℕ × ℕ) : ℝ :=
  diagonalCoefficient(matrixFirstIndex(p), matrixSecondIndex(p))
proof
  have first := matrixFirstIndex.selected with member
  have second := matrixSecondIndex.selected with member
  rcases first with ⟨naturalFirst, j, naturalJ, shapeFirst⟩
  rcases second with ⟨naturalSecond, i, naturalI, shapeSecond⟩
  have typedFirst := nat.omega.backward with naturalFirst
  have typedSecond := nat.omega.backward with naturalSecond
  obtain ⟨real, diagonal, offDiagonal⟩ := diagonalCoefficient.selected with typedFirst, typedSecond
  exact real

theorem identityMatrix.matrix (n : Nat)
  shows RealMatrix(identityMatrix(n))
proof
  unfold RealMatrix
  exact identityMatrix.function with assumption

theorem identityMatrix.entry (n i j : Nat)
  shows identityMatrix(n)[i, j] = diagonalCoefficient(i, j)
proof
  have naturalI : i ∈ ℕ by simp [Nat]
  have naturalJ : j ∈ ℕ by simp [Nat]
  have pair := product.pair with naturalI, naturalJ
  have value := identityMatrix.value at n, (i, j)
  rw [value, matrixFirstIndex.pair, matrixSecondIndex.pair]
  rfl

theorem identityMatrix.diagonal (n i : Nat)
  shows identityMatrix(n)[i, i] = (1:R)
proof
  rw identityMatrix.entry
  exact diagonalCoefficient.same with assumption

theorem identityMatrix.off_diagonal (n i j : Nat)
  requires different : i ≠ j
  shows identityMatrix(n)[i, j] = (0:R)
proof
  rw identityMatrix.entry
  exact diagonalCoefficient.distinct with *

theorem matrixAdd.commutative_entry (A B : Set, i j : Nat)
  requires matrixA : RealMatrix(A)
  requires matrixB : RealMatrix(B)
  shows matrixAdd(A, B)[i, j] = matrixAdd(B, A)[i, j]
proof
  have typedA := RealMatrix.entry at A, i, j with matrixA
  have typedB := RealMatrix.entry at B, i, j with matrixB
  rw [matrixAdd.entry, matrixAdd.entry]
  ring

theorem matrixAdd.associative_entry (A B C : Set, i j : Nat)
  requires matrixA : RealMatrix(A)
  requires matrixB : RealMatrix(B)
  requires matrixC : RealMatrix(C)
  shows matrixAdd(matrixAdd(A, B), C)[i, j] =
    matrixAdd(A, matrixAdd(B, C))[i, j]
proof
  have typedAB := matrixAdd.matrix with matrixA, matrixB
  have typedBC := matrixAdd.matrix with matrixB, matrixC
  have typedA := RealMatrix.entry at A, i, j with matrixA
  have typedB := RealMatrix.entry at B, i, j with matrixB
  have typedC := RealMatrix.entry at C, i, j with matrixC
  rw [matrixAdd.entry, matrixAdd.entry, matrixAdd.entry, matrixAdd.entry]
  ring

theorem matrixScale.add_entry (A B : Set, C : Real, i j : Nat)
  requires matrixA : RealMatrix(A)
  requires matrixB : RealMatrix(B)
  shows matrixScale(C, matrixAdd(A, B))[i, j] =
    matrixAdd(matrixScale(C, A), matrixScale(C, B))[i, j]
proof
  have combined := matrixAdd.matrix with matrixA, matrixB
  have scaledA := matrixScale.matrix with matrixA
  have scaledB := matrixScale.matrix with matrixB
  have typedA := RealMatrix.entry at A, i, j with matrixA
  have typedB := RealMatrix.entry at B, i, j with matrixB
  rw [matrixScale.entry, matrixAdd.entry, matrixAdd.entry, matrixScale.entry, matrixScale.entry]
  ring

theorem matrixTranspose.add_entry (A B : Set, i j : Nat)
  requires matrixA : RealMatrix(A)
  requires matrixB : RealMatrix(B)
  shows matrixTranspose(matrixAdd(A, B))[i, j] =
    matrixAdd(matrixTranspose(A), matrixTranspose(B))[i, j]
proof
  have combined := matrixAdd.matrix with matrixA, matrixB
  have transposedA := matrixTranspose.matrix with matrixA
  have transposedB := matrixTranspose.matrix with matrixB
  rw [matrixTranspose.entry, matrixAdd.entry, matrixAdd.entry, matrixTranspose.entry, matrixTranspose.entry]
  rfl

theorem matrixTranspose.scale_entry (C : Real, M : Set, i j : Nat)
  requires matrix : RealMatrix(M)
  shows matrixTranspose(matrixScale(C, M))[i, j] =
    matrixScale(C, matrixTranspose(M))[i, j]
proof
  have scaled := matrixScale.matrix with matrix
  have transposed := matrixTranspose.matrix with matrix
  rw [matrixTranspose.entry, matrixScale.entry, matrixScale.entry, matrixTranspose.entry]
  rfl

defgraph matrixRow (M : Set, i : Nat) (j ∈ ℕ) : ℝ := M[i, j]
requires matrix : RealMatrix(M)
proof
  have typedJ := nat.omega.backward with member
  exact RealMatrix.entry at M, i, j with matrix

theorem matrixRow.sequence (M : Set, i : Nat)
  requires matrix : RealMatrix(M)
  shows RealSequence(matrixRow(M, i))
proof
  unfold RealSequence
  exact matrixRow.function with matrix

defgraph matrixColumn (M : Set, j : Nat) (i ∈ ℕ) : ℝ := M[i, j]
requires matrix : RealMatrix(M)
proof
  have typedI := nat.omega.backward with member
  exact RealMatrix.entry at M, i, j with matrix

theorem matrixColumn.sequence (M : Set, j : Nat)
  requires matrix : RealMatrix(M)
  shows RealSequence(matrixColumn(M, j))
proof
  unfold RealSequence
  exact matrixColumn.function with matrix

theorem matrixTranspose.involutive_entry (M : Set, i j : Nat)
  requires matrix : RealMatrix(M)
  shows matrixTranspose(matrixTranspose(M))[i, j] = M[i, j]
proof
  have transposed := matrixTranspose.matrix with matrix
  have outer := matrixTranspose.entry at matrixTranspose(M), i, j with transposed
  have inner := matrixTranspose.entry at M, j, i with matrix
  rwa outer

theorem matrixTranspose.row_column (M : Set, i j : Nat)
  requires matrix : RealMatrix(M)
  shows matrixRow(matrixTranspose(M), i)[j] = matrixColumn(M, i)[j]
proof
  have transposed := matrixTranspose.matrix with matrix
  have entry := matrixTranspose.entry at M, i, j with matrix
  rwa [matrixRow.value, matrixColumn.value]

defgraph matrixProduct (A B : Set, n : Nat) (p ∈ ℕ × ℕ) : ℝ :=
  dotProduct(matrixRow(A, matrixFirstIndex(p)),
    matrixColumn(B, matrixSecondIndex(p)), n)
requires matrixA : RealMatrix(A)
requires matrixB : RealMatrix(B)
proof
  have first := matrixFirstIndex.selected with member
  have second := matrixSecondIndex.selected with member
  rcases first with ⟨naturalFirst, j, naturalJ, shapeFirst⟩
  rcases second with ⟨naturalSecond, i, naturalI, shapeSecond⟩
  have typedFirst := nat.omega.backward with naturalFirst
  have typedSecond := nat.omega.backward with naturalSecond
  have row := matrixRow.sequence at A, matrixFirstIndex(p) with matrixA
  have column := matrixColumn.sequence at B, matrixSecondIndex(p) with matrixB
  have value := dotProduct.real at matrixRow(A, matrixFirstIndex(p)),
    matrixColumn(B, matrixSecondIndex(p)), n with row, column
  exact value

theorem matrixProduct.matrix (A B : Set, n : Nat)
  requires matrixA : RealMatrix(A)
  requires matrixB : RealMatrix(B)
  shows RealMatrix(matrixProduct(A, B, n))
proof
  unfold RealMatrix
  exact matrixProduct.function with matrixA, matrixB

theorem matrixProduct.entry (A B : Set, n i j : Nat)
  requires matrixA : RealMatrix(A)
  requires matrixB : RealMatrix(B)
  shows matrixProduct(A, B, n)[i, j] =
    dotProduct(matrixRow(A, i), matrixColumn(B, j), n)
proof
  have naturalI : i ∈ ℕ by simp [Nat]
  have naturalJ : j ∈ ℕ by simp [Nat]
  have pair := product.pair with naturalI, naturalJ
  have value := matrixProduct.value at A, B, n, (i, j) with matrixA, matrixB
  rw [value, matrixFirstIndex.pair, matrixSecondIndex.pair]
  rfl

theorem matrixProduct.zero (A B : Set, i j : Nat)
  requires matrixA : RealMatrix(A)
  requires matrixB : RealMatrix(B)
  shows matrixProduct(A, B, 0)[i, j] = (0:R)
proof
  have value := matrixProduct.entry at A, B, 0, i, j with matrixA, matrixB
  rw value
  exact dotProduct.zero

theorem matrixProduct.succ (A B : Set, n i j : Nat)
  requires matrixA : RealMatrix(A)
  requires matrixB : RealMatrix(B)
  shows matrixProduct(A, B, succ(n))[i, j] =
    (matrixProduct(A, B, n)[i, j]:R) + (A[i, n]:R) * (B[n, j]:R)
proof
  have row := matrixRow.sequence at A, i with matrixA
  have column := matrixColumn.sequence at B, j with matrixB
  have step := dotProduct.succ at matrixRow(A, i), matrixColumn(B, j), n with row, column
  rw [matrixRow.value, matrixColumn.value] at step
  have next := matrixProduct.entry at A, B, succ(n), i, j with matrixA, matrixB
  have previous := matrixProduct.entry at A, B, n, i, j with matrixA, matrixB
  rwa [next, previous]

theorem matrixProduct.transpose_entry (A B : Set, n i j : Nat)
  requires matrixA : RealMatrix(A)
  requires matrixB : RealMatrix(B)
  shows matrixTranspose(matrixProduct(A, B, n))[i, j] =
    matrixProduct(matrixTranspose(B), matrixTranspose(A), n)[i, j]
proof
  induction n
  | zero =>
    have transposedA := matrixTranspose.matrix with matrixA
    have transposedB := matrixTranspose.matrix with matrixB
    have product := matrixProduct.matrix at A, B, 0 with matrixA, matrixB
    have flip := matrixTranspose.entry at matrixProduct(A, B, 0), i, j with product
    have left := matrixProduct.zero at A, B, j, i with matrixA, matrixB
    have right := matrixProduct.zero at matrixTranspose(B), matrixTranspose(A), i, j with
      transposedB, transposedA
    rw [flip, left, right]
    rfl
  | succ(k) using ih =>
    have transposedA := matrixTranspose.matrix with matrixA
    have transposedB := matrixTranspose.matrix with matrixB
    have product := matrixProduct.matrix at A, B, succ(k) with matrixA, matrixB
    have flip := matrixTranspose.entry at matrixProduct(A, B, succ(k)), i, j with product
    have left := matrixProduct.succ at A, B, k, j, i with matrixA, matrixB
    have right := matrixProduct.succ at
      matrixTranspose(B), matrixTranspose(A), k, i, j with transposedB, transposedA
    have oldProduct := matrixProduct.matrix at A, B, k with matrixA, matrixB
    have oldFlip := matrixTranspose.entry at matrixProduct(A, B, k), i, j with oldProduct
    have Bentry := matrixTranspose.entry at B, i, k with matrixB
    have Aentry := matrixTranspose.entry at A, k, j with matrixA
    have realA := RealMatrix.entry at A, j, k with matrixA
    have realB := RealMatrix.entry at B, k, i with matrixB
    have realOld := RealMatrix.entry at matrixProduct(A, B, k), j, i with oldProduct
    have reverseProduct := matrixProduct.matrix at matrixTranspose(B),
      matrixTranspose(A), k with transposedB, transposedA
    have realReverse := RealMatrix.entry at
      matrixProduct(matrixTranspose(B), matrixTranspose(A), k), i, j with reverseProduct
    rw [flip, left, right, Bentry, Aentry, ←oldFlip, ih]
    ring

theorem matrixProduct.add_right_entry (A B C : Set, i j n : Nat)
  requires matrixA : RealMatrix(A)
  requires matrixB : RealMatrix(B)
  requires matrixC : RealMatrix(C)
  shows matrixProduct(A, matrixAdd(B, C), n)[i, j] =
    matrixAdd(matrixProduct(A, B, n), matrixProduct(A, C, n))[i, j]
proof
  induction n
  | zero =>
    have combined := matrixAdd.matrix with matrixB, matrixC
    have productB := matrixProduct.matrix at A, B, 0 with matrixA, matrixB
    have productC := matrixProduct.matrix at A, C, 0 with matrixA, matrixC
    have left := matrixProduct.zero at A, matrixAdd(B, C), i, j with matrixA, combined
    have rightB := matrixProduct.zero at A, B, i, j with matrixA, matrixB
    have rightC := matrixProduct.zero at A, C, i, j with matrixA, matrixC
    rw [left, matrixAdd.entry, rightB, rightC]
    ring
  | succ(k) using ih =>
    have combined := matrixAdd.matrix with matrixB, matrixC
    have nextB := matrixProduct.matrix at A, B, succ(k) with matrixA, matrixB
    have nextC := matrixProduct.matrix at A, C, succ(k) with matrixA, matrixC
    have oldB := matrixProduct.matrix at A, B, k with matrixA, matrixB
    have oldC := matrixProduct.matrix at A, C, k with matrixA, matrixC
    have oldCombined := matrixProduct.matrix at A, matrixAdd(B, C), k with matrixA, combined
    have left := matrixProduct.succ at A, matrixAdd(B, C), k, i, j with matrixA, combined
    have rightB := matrixProduct.succ at A, B, k, i, j with matrixA, matrixB
    have rightC := matrixProduct.succ at A, C, k, i, j with matrixA, matrixC
    have entry := matrixAdd.entry at B, C, k, j with matrixB, matrixC
    have previous := matrixAdd.entry at matrixProduct(A, B, k),
      matrixProduct(A, C, k), i, j with oldB, oldC
    have realA := RealMatrix.entry at A, i, k with matrixA
    have realB := RealMatrix.entry at B, k, j with matrixB
    have realC := RealMatrix.entry at C, k, j with matrixC
    have realOldB := RealMatrix.entry at matrixProduct(A, B, k), i, j with oldB
    have realOldC := RealMatrix.entry at matrixProduct(A, C, k), i, j with oldC
    have realOldCombined := RealMatrix.entry at matrixProduct(A, matrixAdd(B, C), k),
      i, j with oldCombined
    rw [left, entry, matrixAdd.entry, rightB, rightC]
    rw previous at ih
    rw ih
    ring

theorem matrixProduct.add_left_entry (A B C : Set, i j n : Nat)
  requires matrixA : RealMatrix(A)
  requires matrixB : RealMatrix(B)
  requires matrixC : RealMatrix(C)
  shows matrixProduct(matrixAdd(A, B), C, n)[i, j] =
    matrixAdd(matrixProduct(A, C, n), matrixProduct(B, C, n))[i, j]
proof
  induction n
  | zero =>
    have combined := matrixAdd.matrix with matrixA, matrixB
    have productA := matrixProduct.matrix at A, C, 0 with matrixA, matrixC
    have productB := matrixProduct.matrix at B, C, 0 with matrixB, matrixC
    have left := matrixProduct.zero at matrixAdd(A, B), C, i, j with combined, matrixC
    have rightA := matrixProduct.zero at A, C, i, j with matrixA, matrixC
    have rightB := matrixProduct.zero at B, C, i, j with matrixB, matrixC
    rw [left, matrixAdd.entry, rightA, rightB]
    ring
  | succ(k) using ih =>
    have combined := matrixAdd.matrix with matrixA, matrixB
    have nextA := matrixProduct.matrix at A, C, succ(k) with matrixA, matrixC
    have nextB := matrixProduct.matrix at B, C, succ(k) with matrixB, matrixC
    have oldA := matrixProduct.matrix at A, C, k with matrixA, matrixC
    have oldB := matrixProduct.matrix at B, C, k with matrixB, matrixC
    have left := matrixProduct.succ at matrixAdd(A, B), C, k, i, j with combined, matrixC
    have rightA := matrixProduct.succ at A, C, k, i, j with matrixA, matrixC
    have rightB := matrixProduct.succ at B, C, k, i, j with matrixB, matrixC
    have entry := matrixAdd.entry at A, B, i, k with matrixA, matrixB
    have previous := matrixAdd.entry at matrixProduct(A, C, k),
      matrixProduct(B, C, k), i, j with oldA, oldB
    have realA := RealMatrix.entry at A, i, k with matrixA
    have realB := RealMatrix.entry at B, i, k with matrixB
    have realC := RealMatrix.entry at C, k, j with matrixC
    have realOldA := RealMatrix.entry at matrixProduct(A, C, k), i, j with oldA
    have realOldB := RealMatrix.entry at matrixProduct(B, C, k), i, j with oldB
    rw [left, entry, matrixAdd.entry, rightA, rightB]
    rw previous at ih
    rw ih
    ring

-- A matrix with n columns acts on a vector by the dot product with each row.
defgraph matrixVectorProduct (M u : Set, n : Nat) (i ∈ ℕ) : ℝ :=
  dotProduct(matrixRow(M, i), u, n)
requires matrix : RealMatrix(M)
requires vector : RealSequence(u)
proof
  have typedI := nat.omega.backward with member
  have row := matrixRow.sequence at M, i with matrix
  exact dotProduct.real at matrixRow(M, i), u, n with row, vector

theorem matrixVectorProduct.sequence (M u : Set, n : Nat)
  requires matrix : RealMatrix(M)
  requires vector : RealSequence(u)
  shows RealSequence(matrixVectorProduct(M, u, n))
proof
  unfold RealSequence
  exact matrixVectorProduct.function with matrix, vector

theorem matrixProduct.column_action (A B : Set, n i j : Nat)
  requires matrixA : RealMatrix(A)
  requires matrixB : RealMatrix(B)
  shows matrixProduct(A, B, n)[i, j] =
    matrixVectorProduct(A, matrixColumn(B, j), n)[i]
proof
  have column := matrixColumn.sequence at B, j with matrixB
  have product := matrixProduct.entry at A, B, n, i, j with matrixA, matrixB
  rw [product, matrixVectorProduct.value]
  rfl

theorem matrixVectorProduct.zero (M u : Set, i : Nat)
  requires matrix : RealMatrix(M)
  requires vector : RealSequence(u)
  shows matrixVectorProduct(M, u, 0)[i] = (0:R)
proof
  rw matrixVectorProduct.value
  exact dotProduct.zero

theorem matrixVectorProduct.succ (M u : Set, n i : Nat)
  requires matrix : RealMatrix(M)
  requires vector : RealSequence(u)
  shows matrixVectorProduct(M, u, succ(n))[i] =
    (matrixVectorProduct(M, u, n)[i]:R) + (M[i, n]:R) * (u[n]:R)
proof
  have row := matrixRow.sequence at M, i with matrix
  have step := dotProduct.succ at matrixRow(M, i), u, n with row, vector
  rw matrixRow.value at step
  rwa [matrixVectorProduct.value, matrixVectorProduct.value]

theorem matrixVectorProduct.add (M u v : Set, n i : Nat)
  requires matrix : RealMatrix(M)
  requires vectorU : RealSequence(u)
  requires vectorV : RealSequence(v)
  shows matrixVectorProduct(M, addSequence(u, v), n)[i] =
    (matrixVectorProduct(M, u, n)[i]:R) +
    (matrixVectorProduct(M, v, n)[i]:R)
proof
  have row := matrixRow.sequence at M, i with matrix
  have combined := addSequence.sequence with vectorU, vectorV
  have distributive := dotProduct.add_right at matrixRow(M, i), u, v, n with
    row, vectorU, vectorV
  rwa [matrixVectorProduct.value, matrixVectorProduct.value, matrixVectorProduct.value]

theorem matrixVectorProduct.add_matrix (A B u : Set, i n : Nat)
  requires matrixA : RealMatrix(A)
  requires matrixB : RealMatrix(B)
  requires vector : RealSequence(u)
  shows matrixVectorProduct(matrixAdd(A, B), u, n)[i] =
    (matrixVectorProduct(A, u, n)[i]:R) + (matrixVectorProduct(B, u, n)[i]:R)
proof
  induction n
  | zero =>
    have combined := matrixAdd.matrix with matrixA, matrixB
    have left := matrixVectorProduct.zero with combined, vector
    have rightA := matrixVectorProduct.zero with matrixA, vector
    have rightB := matrixVectorProduct.zero with matrixB, vector
    rw [left, rightA, rightB]
    ring
  | succ(k) using ih =>
    have combined := matrixAdd.matrix with matrixA, matrixB
    have left := matrixVectorProduct.succ at matrixAdd(A, B), u, k, i with combined, vector
    have rightA := matrixVectorProduct.succ at A, u, k, i with matrixA, vector
    have rightB := matrixVectorProduct.succ at B, u, k, i with matrixB, vector
    have entry := matrixAdd.entry at A, B, i, k with matrixA, matrixB
    have typedA := RealMatrix.entry at A, i, k with matrixA
    have typedB := RealMatrix.entry at B, i, k with matrixB
    have typedU := RealSequence.value at u, k with vector
    have seqA := matrixVectorProduct.sequence at A, u, k with matrixA, vector
    have seqB := matrixVectorProduct.sequence at B, u, k with matrixB, vector
    have oldA := RealSequence.value at matrixVectorProduct(A, u, k), i with seqA
    have oldB := RealSequence.value at matrixVectorProduct(B, u, k), i with seqB
    rw [left, rightA, rightB, entry, ih]
    ring

theorem matrixVectorProduct.scale_matrix (M u : Set, C : Real, i n : Nat)
  requires matrix : RealMatrix(M)
  requires vector : RealSequence(u)
  shows matrixVectorProduct(matrixScale(C, M), u, n)[i] =
    C * (matrixVectorProduct(M, u, n)[i]:R)
proof
  induction n
  | zero =>
    have scaled := matrixScale.matrix with matrix
    have left := matrixVectorProduct.zero with scaled, vector
    have right := matrixVectorProduct.zero with matrix, vector
    rw [left, right]
    ring
  | succ(k) using ih =>
    have scaled := matrixScale.matrix with matrix
    have left := matrixVectorProduct.succ at matrixScale(C, M), u, k, i with scaled, vector
    have right := matrixVectorProduct.succ at M, u, k, i with matrix, vector
    have entry := matrixScale.entry at C, M, i, k with matrix
    have typedM := RealMatrix.entry at M, i, k with matrix
    have typedU := RealSequence.value at u, k with vector
    have sequence := matrixVectorProduct.sequence at M, u, k with matrix, vector
    have old := RealSequence.value at matrixVectorProduct(M, u, k), i with sequence
    rw [left, right, entry, ih]
    ring

theorem matrixVectorProduct.scale (M u : Set, C : Real, n i : Nat)
  requires matrix : RealMatrix(M)
  requires vector : RealSequence(u)
  shows matrixVectorProduct(M, mulSequence(constantSequence(C), u), n)[i] =
    C * (matrixVectorProduct(M, u, n)[i]:R)
proof
  have row := matrixRow.sequence at M, i with matrix
  have constants := constantSequence.sequence at C
  have scaled := mulSequence.sequence with constants, vector
  have distributive := dotProduct.scale_right at matrixRow(M, i), u, C, n with row, vector
  rwa [matrixVectorProduct.value, matrixVectorProduct.value]

theorem matrixVectorProduct.eq_prefix (M u v : Set, n i : Nat)
  requires matrix : RealMatrix(M)
  requires vectorU : RealSequence(u)
  requires vectorV : RealSequence(v)
  requires equal : ∀ k ∈ ℕ, k < n → u[k] = v[k]
  shows matrixVectorProduct(M, u, n)[i] = matrixVectorProduct(M, v, n)[i]
proof
  have row := matrixRow.sequence at M, i with matrix
  have sameRow : ∀ k ∈ ℕ, k < n →
      matrixRow(M, i)[k] = matrixRow(M, i)[k]
    intros k naturalK before
    rfl
  have prefix := dotProduct.eq_prefix with row, vectorU, row, vectorV, sameRow, equal
  rwa [matrixVectorProduct.value, matrixVectorProduct.value]

theorem matrixVectorProduct.equal_rows (M N u : Set, n i : Nat)
  requires matrixM : RealMatrix(M)
  requires matrixN : RealMatrix(N)
  requires vector : RealSequence(u)
  requires equal : ∀ j ∈ ℕ, j < n → M[i, j] = N[i, j]
  shows matrixVectorProduct(M, u, n)[i] = matrixVectorProduct(N, u, n)[i]
proof
  have rowM := matrixRow.sequence at M, i with matrixM
  have rowN := matrixRow.sequence at N, i with matrixN
  have sameRow : ∀ k ∈ ℕ, k < n →
      matrixRow(M, i)[k] = matrixRow(N, i)[k]
    intros k naturalK before
    have entry := equal with k, naturalK, before
    rwa [matrixRow.value, matrixRow.value]
  have sameVector : ∀ k ∈ ℕ, k < n → u[k] = u[k]
    intros k naturalK before
    rfl
  have prefix := dotProduct.eq_prefix with rowM, vector, rowN, vector, sameRow, sameVector
  rwa [matrixVectorProduct.value, matrixVectorProduct.value]

-- Before its diagonal coordinate, an identity row contributes only zeros.
theorem identityMatrix.action_before_aux (m : Nat)
  shows ∀ n ∈ ℕ, ∀ i ∈ ℕ, ∀ u, RealSequence(u) → m ≤ i →
    matrixVectorProduct(identityMatrix(n), u, m)[i] = (0:R)
proof
  induction m
  | zero =>
    intros n naturalN i naturalI u sequence outside
    have typedN := nat.omega.backward with naturalN
    have typedI := nat.omega.backward with naturalI
    have matrix := identityMatrix.matrix with typedN
    exact matrixVectorProduct.zero at identityMatrix(n), u, i with matrix, sequence
  | succ(k) using ih =>
    intros n naturalN i naturalI u sequence outside
    have typedN := nat.omega.backward with naturalN
    have typedI := nat.omega.backward with naturalI
    have smaller : k < i by exact outside
    have before := Less.implies_LessEqual at k, i with typedI, smaller
    have different : i ≠ k
      intro equal
      rw equal at smaller
      have impossible := Less.irreflexive at k
      exact impossible with smaller
    have matrix := identityMatrix.matrix with typedN
    have prefix := ih with n, naturalN, i, naturalI, u, sequence, before
    have zeroEntry := identityMatrix.off_diagonal with typedN, typedI, different
    have term := RealSequence.value at u, k with sequence
    have products := matrixVectorProduct.sequence at identityMatrix(n), u, k with matrix, sequence
    have old := RealSequence.value at matrixVectorProduct(identityMatrix(n), u, k), i with products
    rw [matrixVectorProduct.succ, prefix, zeroEntry]
    ring

theorem identityMatrix.action_before (n m i : Nat, u : Set)
  requires vector : RealSequence(u)
  requires outside : m ≤ i
  shows matrixVectorProduct(identityMatrix(n), u, m)[i] = (0:R)
proof
  exact identityMatrix.action_before_aux at m with
    n, nat.omega.forward at n with assumption,
    i, nat.omega.forward at i with assumption, u, vector, outside

theorem identityMatrix.action_at_aux (m : Nat)
  shows ∀ n ∈ ℕ, ∀ i ∈ ℕ, ∀ u, RealSequence(u) → i < m →
    matrixVectorProduct(identityMatrix(n), u, m)[i] = u[i]
proof
  induction m
  | zero =>
    intros n naturalN i naturalI u sequence before
    have typedI := nat.omega.backward with naturalI
    have weak := Grace.Core.Nat.Order.LessEqual.zero with typedI
    have typedZero : Nat(0) by typecheck
    have cycle := Less.transitive_LessEqual with typedI, typedZero, typedI, before, weak
    have impossible := Less.irreflexive with typedI
    exfalso
    exact impossible with cycle
  | succ(k) using ih =>
    intros n naturalN i naturalI u sequence before
    have typedN := nat.omega.backward with naturalN
    have typedI := nat.omega.backward with naturalI
    have matrix := identityMatrix.matrix with typedN
    have term := RealSequence.value at u, k with sequence
    have products := matrixVectorProduct.sequence at identityMatrix(n), u, k with matrix, sequence
    have old := RealSequence.value at matrixVectorProduct(identityMatrix(n), u, k), i with products
    have member := membership.of_Less at i, succ(k) with
      typedI, Naturals.constructor.succ at k with assumption, before
    have alternatives := membership.successor_forward{n := k} with member
    cases alternatives with earlier | equal
    . have smaller := Less.of_membership at i, k with typedI, earlier
      have prefix := ih with n, naturalN, i, naturalI, u, sequence, smaller
      have different : i ≠ k
        intro same
        rw same at smaller
        have impossible := Less.irreflexive at k
        exact impossible with smaller
      have zeroEntry := identityMatrix.off_diagonal with typedN, typedI, different
      rw [matrixVectorProduct.succ, prefix, zeroEntry]
      ring
    . rw equal
      have reflexive := LessEqual.reflexive at k
      have prefix := identityMatrix.action_before at n, k, k, u with
        sequence, reflexive
      have diagonal := identityMatrix.diagonal at n, k with typedN
      rw [matrixVectorProduct.succ, prefix, diagonal]
      ring

theorem identityMatrix.action (n i : Nat, u : Set)
  requires vector : RealSequence(u)
  requires inRange : i < n
  shows matrixVectorProduct(identityMatrix(n), u, n)[i] = u[i]
proof
  exact identityMatrix.action_at_aux at n with
    n, nat.omega.forward at n with assumption,
    i, nat.omega.forward at i with assumption, u, vector, inRange

theorem matrixProduct.identity_left (A : Set, n i j : Nat)
  requires matrix : RealMatrix(A)
  requires inRange : i < n
  shows matrixProduct(identityMatrix(n), A, n)[i, j] = A[i, j]
proof
  have identity := identityMatrix.matrix at n
  have column := matrixColumn.sequence at A, j with matrix
  have action := identityMatrix.action with column, inRange
  have product := matrixProduct.column_action at identityMatrix(n), A, n, i, j with
    identity, matrix
  rw [product, action, matrixColumn.value]
  rfl

theorem diagonalCoefficient.not_equal_symmetric (i j : Nat)
  requires different : i ≠ j
  shows j ≠ i
proof
  simp

theorem diagonalCoefficient.symmetric (i j : Nat)
  shows diagonalCoefficient(i, j) = diagonalCoefficient(j, i)
proof
  have alternatives : i = j ∨ i ≠ j by tauto
  cases alternatives with equal | different
  . rw equal
    rfl
  . have reverse := diagonalCoefficient.not_equal_symmetric with different
    have first := diagonalCoefficient.distinct with different
    have second := diagonalCoefficient.distinct with reverse
    rw [first, second]
    rfl

theorem identityMatrix.symmetric_entry (n i j : Nat)
  shows identityMatrix(n)[i, j] = identityMatrix(n)[j, i]
proof
  rw [identityMatrix.entry, identityMatrix.entry]
  exact diagonalCoefficient.symmetric with *

theorem matrixProduct.identity_right (A : Set, n i j : Nat)
  requires matrix : RealMatrix(A)
  requires inRange : j < n
  shows matrixProduct(A, identityMatrix(n), n)[i, j] = A[i, j]
proof
  have identity := identityMatrix.matrix at n
  have rowA := matrixRow.sequence at A, i with matrix
  have rowI := matrixRow.sequence at identityMatrix(n), j with identity
  have colI := matrixColumn.sequence at identityMatrix(n), j with identity
  have equal : ∀ k ∈ ℕ, k < n →
      matrixColumn(identityMatrix(n), j)[k] = matrixRow(identityMatrix(n), j)[k]
    intros k naturalK before
    have typedK := nat.omega.backward with naturalK
    have symmetric := identityMatrix.symmetric_entry at n, k, j with
      assumption, typedK, assumption
    rwa [matrixColumn.value, matrixRow.value]
  have sameRow : ∀ k ∈ ℕ, k < n →
      matrixRow(A, i)[k] = matrixRow(A, i)[k]
    intros k naturalK before
    rfl
  have equalDot := dotProduct.eq_prefix with colI, rowA, rowI, rowA, equal, sameRow
  have swap := dotProduct.symmetric at matrixRow(A, i),
    matrixColumn(identityMatrix(n), j), n with rowA, colI
  have product := matrixProduct.entry at A, identityMatrix(n), n, i, j with matrix, identity
  have action := identityMatrix.action with rowA, inRange
  rw [product, swap, equalDot, ←matrixVectorProduct.value, action, matrixRow.value]
  rfl

-- The two orders for summing a finite rectangle of a total matrix.
defgraph matrixRowTotals (M : Set, m : Nat) (i ∈ ℕ) : ℝ :=
  realPartialSum(matrixRow(M, i), m)
requires matrix : RealMatrix(M)
proof
  have typedI := nat.omega.backward with member
  have row := matrixRow.sequence at M, i with matrix
  exact realPartialSum.real at matrixRow(M, i), m with row

theorem matrixRowTotals.sequence (M : Set, m : Nat)
  requires matrix : RealMatrix(M)
  shows RealSequence(matrixRowTotals(M, m))
proof
  unfold RealSequence
  exact matrixRowTotals.function with matrix

defgraph matrixColumnTotals (M : Set, n : Nat) (j ∈ ℕ) : ℝ :=
  realPartialSum(matrixColumn(M, j), n)
requires matrix : RealMatrix(M)
proof
  have typedJ := nat.omega.backward with member
  have column := matrixColumn.sequence at M, j with matrix
  exact realPartialSum.real at matrixColumn(M, j), n with column

theorem matrixColumnTotals.sequence (M : Set, n : Nat)
  requires matrix : RealMatrix(M)
  shows RealSequence(matrixColumnTotals(M, n))
proof
  unfold RealSequence
  exact matrixColumnTotals.function with matrix

theorem matrixColumnTotals.succ (M : Set, n j : Nat)
  requires matrix : RealMatrix(M)
  shows matrixColumnTotals(M, succ(n))[j] =
    (matrixColumnTotals(M, n)[j]:R) + (M[n, j]:R)
proof
  have column := matrixColumn.sequence at M, j with matrix
  have next := matrixColumnTotals.value at M, succ(n), j with matrix
  have previous := matrixColumnTotals.value at M, n, j with matrix
  have entry := matrixColumn.value at M, j, n with matrix
  rw [next, previous, realPartialSum.succ, entry]
  unfold realAddValue
  unfold call2
  rfl

theorem matrixColumnTotals.empty_sum (M : Set, m : Nat)
  requires matrix : RealMatrix(M)
  shows realPartialSum(matrixColumnTotals(M, 0), m) = (0:R)
proof
  induction m
  | zero =>
    rw realPartialSum.zero
    rfl
  | succ(k) using ih =>
    have value := matrixColumnTotals.value at M, 0, k with matrix
    have zero := realPartialSum.zero at matrixColumn(M, k)
    have sequence := matrixColumnTotals.sequence with matrix
    have term := RealSequence.value with sequence
    rw [realPartialSum.succ, ih, value, zero]
    have unit := realAddValue.zero_right with realZero.real
    unfold realAddValue
    unfold call2 at unit
    exact unit

theorem matrixTotals.interchange_aux (n : Nat)
  shows ∀ M, RealMatrix(M) → ∀ m ∈ ℕ,
    realPartialSum(matrixRowTotals(M, m), n) =
      realPartialSum(matrixColumnTotals(M, n), m)
proof
  induction n
  | zero =>
    intros M matrix m naturalM
    have typedM := nat.omega.backward with naturalM
    rw realPartialSum.zero
    have empty := matrixColumnTotals.empty_sum with matrix
    rw empty
    rfl
  | succ(k) using ih =>
    intros M matrix m naturalM
    have typedM := nat.omega.backward with naturalM
    have row := matrixRow.sequence at M, k with matrix
    have earlier := matrixColumnTotals.sequence at M, k with matrix
    have later := matrixColumnTotals.sequence at M, succ(k) with matrix
    have combined := addSequence.sequence with earlier, row
    have equal : ∀ j ∈ ℕ, j < m →
        matrixColumnTotals(M, succ(k))[j] =
          addSequence(matrixColumnTotals(M, k), matrixRow(M, k))[j]
      intros j naturalJ before
      have typedJ := nat.omega.backward with naturalJ
      have step := matrixColumnTotals.succ at M, k, j with matrix
      have left := RealSequence.value at matrixColumnTotals(M, k), j with earlier
      have right := RealSequence.value at matrixRow(M, k), j with row
      have value := addSequence.value with earlier, row, left, right
      have entry := matrixRow.value at M, k, j with matrix
      rw [step, value, entry]
      rfl
    have prefix := realPartialSum.eq_prefix with later, combined, equal
    have sumAdd := realPartialSum.add at matrixColumnTotals(M, k),
      matrixRow(M, k), m with earlier, row
    have inductionStep := ih with M, matrix, m, naturalM
    have rowTotal := matrixRowTotals.value at M, m, k with matrix
    rw [realPartialSum.succ, rowTotal, inductionStep, prefix, sumAdd]
    unfold realAddValue
    unfold call2
    rfl

theorem matrixTotals.interchange (M : Set, n m : Nat)
  requires matrix : RealMatrix(M)
  shows realPartialSum(matrixRowTotals(M, m), n) =
    realPartialSum(matrixColumnTotals(M, n), m)
proof
  exact matrixTotals.interchange_aux at n with M, matrix,
    m, nat.omega.forward at m with assumption