-- Explicit solutions of two real linear equations.
-- 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.TwoByTwo
import Real.Inverse

namespace LinearSystem2

-- These total expressions are solutions only when the determinant is nonzero.
def solutionX (a b c d u v : Real) : Real :=
  realInvValue(a * d - b * c) * (d * u - b * v)

def solutionY (a b c d u v : Real) : Real :=
  realInvValue(a * d - b * c) * (a * v - c * u)

theorem solutionX.real (a b c d u v : Real)
  shows Real(solutionX(a, b, c, d, u, v))
proof
  unfold solutionX
  typecheck

theorem solutionY.real (a b c d u v : Real)
  shows Real(solutionY(a, b, c, d, u, v))
proof
  unfold solutionY
  typecheck

theorem correct (a b c d u v : Real)
  requires nonsingular : a * d - b * c ≠ 0
  shows a * solutionX(a, b, c, d, u, v) + b * solutionY(a, b, c, d, u, v) = u ∧
    c * solutionX(a, b, c, d, u, v) + d * solutionY(a, b, c, d, u, v) = v
proof
  have inverse := realInvValue.mul_left with *
  constructor
  . unfold solutionX
    unfold solutionY
    linear_combination u * inverse
  . unfold solutionX
    unfold solutionY
    linear_combination v * inverse

-- Uniqueness applies directly to any two solutions, without unfolding Cramer formulas.
theorem unique (a b c d u v x y s t : Real)
  requires nonsingular : a * d - b * c ≠ 0
  requires first : a * x + b * y = u
  requires second : c * x + d * y = v
  requires otherFirst : a * s + b * t = u
  requires otherSecond : c * s + d * t = v
  shows x = s ∧ y = t
proof
  have inverse := realInvValue.mul_left with *
  constructor
  . linear_combination (realInvValue(a * d - b * c) * d) * first - (realInvValue(a * d - b * c) * b) * second - (realInvValue(a * d - b * c) * d) * otherFirst + (realInvValue(a * d - b * c) * b) * otherSecond - (x - s) * inverse
  . linear_combination (realInvValue(a * d - b * c) * a) * second - (realInvValue(a * d - b * c) * c) * first - (realInvValue(a * d - b * c) * a) * otherSecond + (realInvValue(a * d - b * c) * c) * otherFirst - (y - t) * inverse

end LinearSystem2