-- Cramer's rule and the adjugate inverse over a bundled commutative ring:
-- for a solution x of the row system x · M = b, det(M) · x(i) is the
-- determinant of M with row i replaced by b, and when det(M) is a unit
-- the matrix det(M)⁻¹ · adj(M) is a right inverse of M; everything here
-- is ZF
-- 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.Algebra.Ring.MatrixCramer
import MatrixEntries

-- The scalar multiple a · M, entry by entry.
defgraph scalarMultiple (R n a M : Set) (p ∈ n × n) : Ring.carrier(R) :=
  call2(Ring.mul(R), a, M(p))
requires ring : Ring(R)
requires scalar : a ∈ Ring.carrier(R)
requires matrix : M : R ^ n × n
proof
  have function := Matrix.function with matrix
  have entry := FunctionBetween.apply_closed with function, member
  exact Ring.mul_mem with *

-- The row vector x · M of a vector x : n → R.
defgraph rowProduct (R : Set, n : Nat, x M : Set) (j ∈ n) : Ring.carrier(R) :=
  rowVectorEntry(n, Ring.carrier(R), Ring.add(R), Ring.zero(R), Ring.neg(R),
    Ring.mul(R), Ring.one(R), x, M, j)
requires comm : CommRing(R)
requires vector : x : n → Ring.carrier(R)
requires matrix : M : R ^ n × n
proof
  have parts := CommRing.parts with comm
  have square := Matrix.square with matrix
  have shape := squareMatrices.membership.forward with square
  have finite := Finite.natural at n
  have weights := rowVectorWeights.function at n, Ring.carrier(R), Ring.add(R),
    Ring.zero(R), Ring.neg(R), Ring.mul(R), Ring.one(R), x, M, j with *
  unfold rowVectorEntry
  exact finiteRingSum.closed with *

section Scalars
variable (R : Set) (n : Nat)
hypothesis comm : CommRing(R)

theorem scalarMultiple.entry (a : R) (M : R ^ n × n) (i j : Set)
  requires row : i ∈ n
  requires column : j ∈ n
  shows scalarMultiple(R, n, a, M)[i, j] = a * M[i, j]
proof
  have pair := product.pair at i, j, n, n with *
  exact scalarMultiple.value at R, n, a, M, (i, j) with *

instance scalarMultiple.matrix (a : R) (M : R ^ n × n)
  shows scalarMultiple(R, n, a, M) : R ^ n × n
proof
  have function := scalarMultiple.function at R, n, a, M with *
  exact Matrix.of_function with function

theorem rowProduct.entry (x : Set) (M : R ^ n × n) (j : Set)
  requires vector : x : n → R.carrier
  requires column : j ∈ n
  shows rowProduct(R, n, x, M)(j) = rowVectorEntry(n, R.carrier, R.add, R.zero,
    R.neg, R.mul, R.one, x, M, j)
proof
  exact rowProduct.value with *

end Scalars

-- The terms x(k) · M[k, j] · C(j), C the cofactor row, as an array over
-- n × n indexed by (j, k).
defunique cramerWeightValue (R : Set, n : Nat, x M i q : Set) (a ∈ Ring.carrier(R)) :=
  ∃ j ∈ n, ∃ k ∈ n, q = (j, k) ∧
    a = call2(Ring.mul(R), x(k), call2(Ring.mul(R), M[k, j],
      cofactorRow(n, Ring.carrier(R), Ring.add(R), Ring.zero(R), Ring.neg(R),
        Ring.mul(R), Ring.one(R), M, i)(j)))
requires comm : CommRing(R)
requires vector : x : n → Ring.carrier(R)
requires matrix : M : R ^ n × n
requires row : i ∈ n
requires pair : q ∈ n × n
proof
  obtain ⟨bound, j, hj, k, hk, shape⟩ := product.membership.forward at
    n, n, q with *
  have parts := CommRing.parts with comm
  have square := Matrix.square with matrix
  have shapeM := squareMatrices.membership.forward with square
  have finite := Finite.natural at n
  have cofactors := cofactorRow.function at n, Ring.carrier(R), Ring.add(R),
    Ring.zero(R), Ring.neg(R), Ring.mul(R), Ring.one(R), M, i with *
  have cofactor := FunctionBetween.apply_closed with cofactors, hj
  have coordinate := FunctionBetween.apply_closed with vector, hk
  have entry := Matrix.entry at M, R, n, k, j with *
  have ordinary := IsCommRing.ring with parts
  have inner := IsRing.mul_closed with ordinary, entry, cofactor
  have value := IsRing.mul_closed with ordinary, coordinate, inner
  refine ⟨call2(Ring.mul(R), x(k), call2(Ring.mul(R), M[k, j],
    cofactorRow(n, Ring.carrier(R), Ring.add(R), Ring.zero(R), Ring.neg(R),
      Ring.mul(R), Ring.one(R), M, i)(j))), value, j, hj, k, hk, shape, rfl⟩
uniqueness
  intros u hu first v hv second
  obtain ⟨j, hj, k, hk, shapeU, valueU⟩ := first
  obtain ⟨a, ha, b, hb, shapeV, valueV⟩ := second
  have equal : (j, k) = (a, b)
    cc
  have coordinates := orderedPair.injective at j, k, a, b with equal
  cc

defgraph cramerWeights (R : Set, n : Nat, x M i : Set) (q ∈ n × n) : Ring.carrier(R) :=
  cramerWeightValue(R, n, x, M, i, q)
requires comm : CommRing(R)
requires vector : x : n → Ring.carrier(R)
requires matrix : M : R ^ n × n
requires row : i ∈ n
proof
  have selected := cramerWeightValue.selected at R, n, x, M, i, q with *
  exact selected.left

theorem cramerWeights.pair (R : Set, n : Nat, x M i j k : Set)
  requires comm : CommRing(R)
  requires vector : x : n → Ring.carrier(R)
  requires matrix : M : R ^ n × n
  requires row : i ∈ n
  requires hj : j ∈ n
  requires hk : k ∈ n
  shows cramerWeights(R, n, x, M, i)((j, k)) =
    x(k) * (M[k, j] * cofactorRow(n, R.carrier, R.add, R.zero, R.neg, R.mul,
      R.one, M, i)(j))
proof
  have pair := product.pair at j, k, n, n with *
  have selected := cramerWeightValue.selected at R, n, x, M, i, (j, k) with *
  obtain ⟨a, ha, b, hb, shape, result⟩ := selected.right
  have coordinates := orderedPair.injective with shape
  have value := cramerWeights.value at R, n, x, M, i, (j, k) with *
  cc

section Cramer
variable (R : Set) (n : Nat)
hypothesis comm : CommRing(R)

-- Cramer's rule: a solution x of x · M = b has det(M) · x(i) equal to the
-- determinant of the matrix N obtained from M by replacing row i by b.
theorem Matrix.cramer (M N : R ^ n × n) (x b i : Set)
  requires vector : x : n → R.carrier
  requires solution : ∀ j ∈ n, rowProduct(R, n, x, M)(j) = b(j)
  requires row : i ∈ n
  requires replaced : ∀ j ∈ n, N[i, j] = b(j)
  requires outside : ∀ k ∈ n, k ≠ i → ∀ j ∈ n, N[k, j] = M[k, j]
  shows det(M) * x(i) = det(N)
proof
  have parts := CommRing.parts with comm
  have squareM := Matrix.square at M, R, n with assumption
  have squareN := Matrix.square at N, R, n with assumption
  have shapeM := squareMatrices.membership.forward with squareM
  have shapeN := squareMatrices.membership.forward with squareN
  have finite := Finite.natural at n
  have adjugateShape := adjugate.matrix at n, R.carrier, R.add, R.zero, R.neg,
    R.mul, R.one, M with *
  have determinantIn := Matrix.det_mem at R, n, M with *
  have coordinate := FunctionBetween.apply_closed with vector, row
  have cofactors := cofactorRow.function at n, R.carrier, R.add, R.zero, R.neg,
    R.mul, R.one, M, i with *
  have array := cramerWeights.function at R, n, x, M, i with *
  have transposed := transposedMatrix.function at n, n, R.carrier,
    cramerWeights(R, n, x, M, i) with *
  have rowsL := finiteRowSums.function at n, n, R.carrier, R.add, R.zero, R.neg,
    R.mul, R.one, cramerWeights(R, n, x, M, i) with *
  have rowsR := finiteRowSums.function at n, n, R.carrier, R.add, R.zero, R.neg,
    R.mul, R.one, transposedMatrix(n, n, R.carrier, cramerWeights(R, n, x, M, i)) with *
  have pairWeights := cofactorPairWeights.function at n, R.carrier, R.add,
    R.zero, R.neg, R.mul, R.one, M, N, i, i with *
  have selector := coefficientSingle.function at n, R.carrier, R.zero,
    det(M) * x(i), i with *
  -- Row expansion of det(N) along the replaced row, as a double sum.
  have expansion := determinant.row_expansion at R.carrier, R.add, R.zero, R.neg,
    R.mul, R.one, n, M, N, i with *
  have leftRows : ∀ j ∈ n, cofactorPairWeights(n, Ring.carrier(R), R.add, R.zero, R.neg, R.mul, R.one, M, N, i, i)(j) =
      finiteRowSums(n, n, Ring.carrier(R), R.add, R.zero, R.neg, R.mul, R.one, cramerWeights(R, n, x, M, i))(j)
    intros j hj
    have cofactor := FunctionBetween.apply_closed with cofactors, hj
    have rowWeights := finiteRowWeights.function at n, n, R.carrier,
      cramerWeights(R, n, x, M, i), j with *
    have inner := rowVectorWeights.function at n, R.carrier, R.add, R.zero,
      R.neg, R.mul, R.one, x, M, j with *
    have scaled : ∀ k ∈ n, finiteRowWeights(n, n, Ring.carrier(R), cramerWeights(R, n, x, M, i), j)(k) =
        cofactorRow(n, Ring.carrier(R), R.add, R.zero, R.neg, R.mul, R.one, M, i)(j) * rowVectorWeights(n, Ring.carrier(R), R.add, R.zero, R.neg, R.mul, R.one, x, M, j)(k)
      intros k hk
      have entry := Matrix.entry at M, R, n, k, j with *
      have coordinateK := FunctionBetween.apply_closed with vector, hk
      have rowWeight := finiteRowWeights.value at n, n, R.carrier,
        cramerWeights(R, n, x, M, i), j, k with *
      have cell := cramerWeights.pair at R, n, x, M, i, j, k with *
      have weight := rowVectorWeights.value at n, R.carrier, R.add, R.zero,
        R.neg, R.mul, R.one, x, M, j, k with *
      rw [rowWeight, cell, weight]
      ring
    have scale := finiteRingSum.scale_eq at n, R.carrier, R.add, R.zero, R.neg,
      R.mul, R.one, rowVectorWeights(n, R.carrier, R.add, R.zero, R.neg, R.mul, R.one, x, M, j),
      finiteRowWeights(n, n, R.carrier, cramerWeights(R, n, x, M, i), j),
      cofactorRow(n, R.carrier, R.add, R.zero, R.neg, R.mul, R.one, M, i)(j) with *
    have rowSum := finiteRowSums.value at n, n, R.carrier, R.add, R.zero, R.neg,
      R.mul, R.one, cramerWeights(R, n, x, M, i), j with *
    have weight := cofactorPairWeights.value at n, R.carrier, R.add, R.zero,
      R.neg, R.mul, R.one, M, N, i, i, j with *
    have entryN := replaced with j, hj
    have solved := solution with j, hj
    have product := rowProduct.entry at R, n, x, M, j with *
    have closed := finiteRingSum.closed at n, R.carrier, R.add, R.zero, R.neg,
      R.mul, R.one, rowVectorWeights(n, R.carrier, R.add, R.zero, R.neg, R.mul, R.one, x, M, j) with *
    unfold rowVectorEntry at product
    rw [weight, rowSum, scale, entryN, ← solved, product]
    ring
  -- After interchanging, the inner sums are the entries of M · adj(M).
  have rightRows : ∀ k ∈ n, finiteRowSums(n, n, Ring.carrier(R), R.add, R.zero, R.neg, R.mul, R.one, transposedMatrix(n, n, Ring.carrier(R), cramerWeights(R, n, x, M, i)))(k) =
      coefficientSingle(n, Ring.carrier(R), R.zero, det(M) * x(i), i)(k)
    intros k hk
    have coordinateK := FunctionBetween.apply_closed with vector, hk
    have rowWeights := finiteRowWeights.function at n, n, R.carrier,
      transposedMatrix(n, n, R.carrier, cramerWeights(R, n, x, M, i)), k with *
    have inner := matrixProductWeights.function at n, R.carrier, R.add, R.zero,
      R.neg, R.mul, R.one, M, adj(M), k, i with *
    have scaled : ∀ j ∈ n, finiteRowWeights(n, n, Ring.carrier(R), transposedMatrix(n, n, Ring.carrier(R), cramerWeights(R, n, x, M, i)), k)(j) =
        x(k) * matrixProductWeights(n, Ring.carrier(R), R.add, R.zero, R.neg, R.mul, R.one, M, adj(M), k, i)(j)
      intros j hj
      have rowWeight := finiteRowWeights.value at n, n, R.carrier,
        transposedMatrix(n, n, R.carrier, cramerWeights(R, n, x, M, i)), k, j with *
      have swap := transposedMatrix.pair at n, n, R.carrier,
        cramerWeights(R, n, x, M, i), j, k with *
      have cell := cramerWeights.pair at R, n, x, M, i, j, k with *
      have weight := matrixProductWeights.value at n, R.carrier, R.add, R.zero,
        R.neg, R.mul, R.one, M, adj(M), k, i, j with *
      have entry := Matrix.adj_entry at R, n, M, j, i with *
      rw [rowWeight, swap, cell, weight, entry]
      rfl
    have scale := finiteRingSum.scale_eq at n, R.carrier, R.add, R.zero, R.neg,
      R.mul, R.one, matrixProductWeights(n, R.carrier, R.add, R.zero, R.neg, R.mul, R.one, M, adj(M), k, i),
      finiteRowWeights(n, n, R.carrier, transposedMatrix(n, n, R.carrier, cramerWeights(R, n, x, M, i)), k),
      x(k) with *
    have rowSum := finiteRowSums.value at n, n, R.carrier, R.add, R.zero, R.neg,
      R.mul, R.one, transposedMatrix(n, n, R.carrier, cramerWeights(R, n, x, M, i)), k with *
    have identity := adjugate.right_identity at R.carrier, R.add, R.zero, R.neg,
      R.mul, R.one, n, M, k, i with *
    have zeroIn : R[0] ∈ R := Ring.zero_mem with *
    unfold matrixProductEntry at identity
    rw [rowSum, scale, identity]
    have alternatives : k = i ∨ k ≠ i by tauto
    cases alternatives with same | different
    . have delta := singleValue.same at R.carrier, R.zero, det(M), i, k with *
      have selected := coefficientSingle.same at n, R.carrier, R.zero, det(M) * x(i), i with *
      rw [delta, same, selected]
      ring
    . have reverse : i ≠ k
        intro equal
        apply different
        cc
      have delta := singleValue.other at R.carrier, R.zero, det(M), i, k with *
      have selected := coefficientSingle.other at n, R.carrier, R.zero, det(M) * x(i), i, k with *
      rw [delta, selected]
      ring
  have sameL := finiteRingSum.congr with leftRows
  have sameR := finiteRingSum.congr with rightRows
  have interchange := finiteDoubleSum.interchange at n, n, R.carrier, R.add, R.zero,
    R.neg, R.mul, R.one, cramerWeights(R, n, x, M, i) with *
  have single := finiteRingSum.single at n, R.carrier, R.add, R.zero, R.neg,
    R.mul, R.one, det(M) * x(i), i with *
  unfold finiteDoubleSum at interchange
  unfold cofactorPairing at expansion
  cc

end Cramer

section Inverse
variable (R : Set) (n : Nat)
hypothesis comm : CommRing(R)

-- The adjugate identity in the matrix ring.
theorem Matrix.mul_adj (M : R ^ n × n)
  shows M * adj(M) = scalarMultiple(R, n, det(M), matrixRing(R, n)[1])
proof
  have parts := CommRing.parts with comm
  have square := Matrix.square at M, R, n with assumption
  have shape := squareMatrices.membership.forward with square
  have finite := Finite.natural at n
  have determinantIn := Matrix.det_mem at R, n, M with *
  have adjugate := Matrix.adj_matrix at R, n, M with *
  have one := Matrix.of_mem at matrixRing(R, n)[1], R, n with (Ring.one_mem with *)
  have product := Matrix.of_mem at M * adj(M), R, n with (Ring.mul_mem with *)
  have scaled := scalarMultiple.matrix at R, n, det(M), matrixRing(R, n)[1] with *
  apply Matrix.ext with product, scaled
  intros i hi j hj
  have entry := matrixRing.mul_entry at R, n, M, adj(M), i, j with *
  have identity := adjugate.right_identity at R.carrier, R.add, R.zero, R.neg,
    R.mul, R.one, n, M, i, j with *
  have scale := scalarMultiple.entry at R, n, det(M), matrixRing(R, n)[1], i, j with *
  have unit := matrixRing.one_entry at R, n, i, j with *
  have zeroIn : R[0] ∈ R := Ring.zero_mem with *
  have oneIn : R[1] ∈ R := Ring.one_mem with *
  rw [entry, identity, scale, unit]
  have alternatives : i = j ∨ i ≠ j by tauto
  cases alternatives with same | different
  . have deltaD := singleValue.same at R.carrier, R.zero, det(M), j, i with *
    have deltaO := singleValue.same at R.carrier, R.zero, R.one, j, i with *
    rw [deltaD, deltaO]
    ring
  . have deltaD := singleValue.other at R.carrier, R.zero, det(M), j, i with *
    have deltaO := singleValue.other at R.carrier, R.zero, R.one, j, i with *
    rw [deltaD, deltaO]
    ring

-- When det(M) has an inverse u, the matrix u · adj(M) is a right inverse.
theorem Matrix.mul_inverse (M : R ^ n × n) (u : R)
  requires inverse : u * det(M) = R[1]
  shows M * scalarMultiple(R, n, u, adj(M)) = matrixRing(R, n)[1]
proof
  have parts := CommRing.parts with comm
  have square := Matrix.square at M, R, n with assumption
  have shape := squareMatrices.membership.forward with square
  have finite := Finite.natural at n
  have determinantIn := Matrix.det_mem at R, n, M with *
  have adjugate := Matrix.adj_matrix at R, n, M with *
  have adjugateShape := adjugate.matrix at n, R.carrier, R.add, R.zero, R.neg,
    R.mul, R.one, M with *
  have scaled := scalarMultiple.matrix at R, n, u, adj(M) with *
  have scaledSquare := Matrix.square at scalarMultiple(R, n, u, adj(M)), R, n with *
  have scaledShape := squareMatrices.membership.forward with scaledSquare
  have one := Matrix.of_mem at matrixRing(R, n)[1], R, n with (Ring.one_mem with *)
  have product := Matrix.of_mem at M * scalarMultiple(R, n, u, adj(M)), R, n with (Ring.mul_mem with *)
  have uIn : u ∈ R by assumption
  apply Matrix.ext with product, one
  intros i hi j hj
  have entry := matrixRing.mul_entry at R, n, M, scalarMultiple(R, n, u, adj(M)), i, j with *
  have weightsA := matrixProductWeights.function at n, R.carrier, R.add, R.zero,
    R.neg, R.mul, R.one, M, adj(M), i, j with *
  have weightsS := matrixProductWeights.function at n, R.carrier, R.add, R.zero,
    R.neg, R.mul, R.one, M, scalarMultiple(R, n, u, adj(M)), i, j with *
  have scaledWeights : ∀ k ∈ n, matrixProductWeights(n, Ring.carrier(R), R.add, R.zero, R.neg, R.mul, R.one, M, scalarMultiple(R, n, u, adj(M)), i, j)(k) =
      u * matrixProductWeights(n, Ring.carrier(R), R.add, R.zero, R.neg, R.mul, R.one, M, adj(M), i, j)(k)
    intros k hk
    have entryM := Matrix.entry at M, R, n, i, k with *
    have entryA := Matrix.entry at adj(M), R, n, k, j with *
    have weightA := matrixProductWeights.value at n, R.carrier, R.add, R.zero,
      R.neg, R.mul, R.one, M, adj(M), i, j, k with *
    have weightS := matrixProductWeights.value at n, R.carrier, R.add, R.zero,
      R.neg, R.mul, R.one, M, scalarMultiple(R, n, u, adj(M)), i, j, k with *
    have scale := scalarMultiple.entry at R, n, u, adj(M), k, j with *
    rw [weightA, weightS, scale]
    ring
  have scale := finiteRingSum.scale_eq at n, R.carrier, R.add, R.zero, R.neg,
    R.mul, R.one, matrixProductWeights(n, R.carrier, R.add, R.zero, R.neg, R.mul, R.one, M, adj(M), i, j),
    matrixProductWeights(n, R.carrier, R.add, R.zero, R.neg, R.mul, R.one, M, scalarMultiple(R, n, u, adj(M)), i, j), u with *
  have identity := adjugate.right_identity at R.carrier, R.add, R.zero, R.neg,
    R.mul, R.one, n, M, i, j with *
  have unit := matrixRing.one_entry at R, n, i, j with *
  have zeroIn : R[0] ∈ R := Ring.zero_mem with *
  have oneIn : R[1] ∈ R := Ring.one_mem with *
  unfold matrixProductEntry at identity
  unfold matrixProductEntry at entry
  rw [entry, unit, scale, identity]
  have alternatives : i = j ∨ i ≠ j by tauto
  cases alternatives with same | different
  . have deltaD := singleValue.same at R.carrier, R.zero, det(M), j, i with *
    have deltaO := singleValue.same at R.carrier, R.zero, R.one, j, i with *
    rw [deltaD, deltaO]
    exact inverse
  . have deltaD := singleValue.other at R.carrier, R.zero, det(M), j, i with *
    have deltaO := singleValue.other at R.carrier, R.zero, R.one, j, i with *
    rw [deltaD, deltaO]
    ring

end Inverse