-- Finite dimensional real linear algebra
-- 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.FiniteDimensional
import Real.FiniteSum
import Real.Squares

-- Vectors are real sequences observed only below the dimension n. This
-- representation permits the same sequence to be reused at several sizes.
def dotProduct (u v : RealSequence, n : Nat) : Real :=
  realPartialSum(mulSequence(u, v), n)

theorem dotProduct.real (u v : Set, n : Nat)
  requires sequenceU : RealSequence(u)
  requires sequenceV : RealSequence(v)
  shows Real(dotProduct(u, v, n))
proof
  have product := mulSequence.sequence with sequenceU, sequenceV
  unfold dotProduct
  exact realPartialSum.real with product

theorem dotProduct.zero (u v : Set)
  shows dotProduct(u, v, 0) = (0:R)
proof
  unfold dotProduct
  exact realPartialSum.zero

theorem dotProduct.succ (u v : Set, n : Nat)
  requires sequenceU : RealSequence(u)
  requires sequenceV : RealSequence(v)
  shows dotProduct(u, v, succ(n)) =
    (dotProduct(u, v, n):R) + (u[n]:R) * (v[n]:R)
proof
  have un := RealSequence.value with sequenceU
  have vn := RealSequence.value with sequenceV
  have value := mulSequence.value with sequenceU, sequenceV, un, vn
  unfold dotProduct
  rw [realPartialSum.succ, value]
  unfold realAddValue
  unfold call2
  rfl

theorem dotProduct.symmetric (u v : Set, n : Nat)
  requires sequenceU : RealSequence(u)
  requires sequenceV : RealSequence(v)
  shows dotProduct(u, v, n) = dotProduct(v, u, n)
proof
  induction n
  | zero =>
    rw [dotProduct.zero, dotProduct.zero]
    rfl
  | succ(k) using ih =>
    have uk := RealSequence.value with sequenceU
    have vk := RealSequence.value with sequenceV
    rw [dotProduct.succ, dotProduct.succ, ih]
    ring

-- Only the first n coordinates are observed by an n-dimensional dot product.
theorem dotProduct.eq_prefix (u v w z : Set, n : Nat)
  requires sequenceU : RealSequence(u)
  requires sequenceV : RealSequence(v)
  requires sequenceW : RealSequence(w)
  requires sequenceZ : RealSequence(z)
  requires equalLeft : ∀ k ∈ ℕ, k < n → u[k] = w[k]
  requires equalRight : ∀ k ∈ ℕ, k < n → v[k] = z[k]
  shows dotProduct(u, v, n) = dotProduct(w, z, n)
proof
  have first := mulSequence.sequence with sequenceU, sequenceV
  have second := mulSequence.sequence with sequenceW, sequenceZ
  have equal : ∀ k ∈ ℕ, k < n →
      mulSequence(u, v)[k] = mulSequence(w, z)[k]
    intros k naturalK before
    have typedK := nat.omega.backward with naturalK
    have uk := RealSequence.value at u, k with sequenceU
    have vk := RealSequence.value at v, k with sequenceV
    have wk := RealSequence.value at w, k with sequenceW
    have zk := RealSequence.value at z, k with sequenceZ
    have firstValue := mulSequence.value with sequenceU, sequenceV, uk, vk
    have secondValue := mulSequence.value with sequenceW, sequenceZ, wk, zk
    have left := equalLeft with k, naturalK, before
    have right := equalRight with k, naturalK, before
    rw [firstValue, secondValue, left, right]
    rfl
  unfold dotProduct
  exact realPartialSum.eq_prefix with first, second, equal

theorem dotProduct.add_right (u v w : Set, n : Nat)
  requires sequenceU : RealSequence(u)
  requires sequenceV : RealSequence(v)
  requires sequenceW : RealSequence(w)
  shows dotProduct(u, addSequence(v, w), n) =
    (dotProduct(u, v, n):R) + (dotProduct(u, w, n):R)
proof
  induction n
  | zero =>
    rw [dotProduct.zero, dotProduct.zero, dotProduct.zero]
    ring
  | succ(k) using ih =>
    have combined := addSequence.sequence with sequenceV, sequenceW
    have uk := RealSequence.value with sequenceU
    have vk := RealSequence.value with sequenceV
    have wk := RealSequence.value with sequenceW
    have vkPlusWk := RealSequence.value with combined
    have value := addSequence.value with sequenceV, sequenceW, vk, wk
    rw [dotProduct.succ, dotProduct.succ, dotProduct.succ, ih, value]
    ring

theorem dotProduct.add_left (u v w : Set, n : Nat)
  requires sequenceU : RealSequence(u)
  requires sequenceV : RealSequence(v)
  requires sequenceW : RealSequence(w)
  shows dotProduct(addSequence(u, v), w, n) =
    (dotProduct(u, w, n):R) + (dotProduct(v, w, n):R)
proof
  have combined := addSequence.sequence with sequenceU, sequenceV
  have symmetryCombined := dotProduct.symmetric with combined, sequenceW
  have distributive := dotProduct.add_right with sequenceW, sequenceU, sequenceV
  have symmetryU := dotProduct.symmetric with sequenceW, sequenceU
  have symmetryV := dotProduct.symmetric with sequenceW, sequenceV
  rw [symmetryCombined, distributive, symmetryU, symmetryV]
  rfl

theorem dotProduct.scale_right (u v : Set, C : Real, n : Nat)
  requires sequenceU : RealSequence(u)
  requires sequenceV : RealSequence(v)
  shows dotProduct(u, mulSequence(constantSequence(C), v), n) =
    C * (dotProduct(u, v, n):R)
proof
  induction n
  | zero =>
    rw [dotProduct.zero, dotProduct.zero]
    ring
  | succ(k) using ih =>
    have constants := constantSequence.sequence at C
    have scaled := mulSequence.sequence with constants, sequenceV
    have uk := RealSequence.value with sequenceU
    have vk := RealSequence.value with sequenceV
    have constantK := RealSequence.value with constants
    have scaledK := RealSequence.value with scaled
    have value := mulSequence.value with constants, sequenceV, constantK, vk
    have constant := constantSequence.value at C, k with *
    rw [dotProduct.succ, dotProduct.succ, ih, value, constant]
    ring

theorem dotProduct.scale_left (u v : Set, C : Real, n : Nat)
  requires sequenceU : RealSequence(u)
  requires sequenceV : RealSequence(v)
  shows dotProduct(mulSequence(constantSequence(C), u), v, n) =
    C * (dotProduct(u, v, n):R)
proof
  have constants := constantSequence.sequence at C
  have scaled := mulSequence.sequence with constants, sequenceU
  have symmetricScaled := dotProduct.symmetric with scaled, sequenceV
  have right := dotProduct.scale_right with sequenceV, sequenceU
  have symmetric := dotProduct.symmetric with sequenceV, sequenceU
  rw [symmetricScaled, right, symmetric]
  rfl

theorem dotProduct.sub_left (u v w : Set, n : Nat)
  requires sequenceU : RealSequence(u)
  requires sequenceV : RealSequence(v)
  requires sequenceW : RealSequence(w)
  shows dotProduct(subSequence(u, v), w, n) =
    (dotProduct(u, w, n):R) - (dotProduct(v, w, n):R)
proof
  induction n
  | zero =>
    simp [dotProduct.zero]
    ring
  | succ(k) using ih =>
    have difference := subSequence.sequence at u, v with *
    have uk := RealSequence.value at u, k with *
    have vk := RealSequence.value at v, k with *
    have wk := RealSequence.value at w, k with *
    have dk := RealSequence.value at subSequence(u, v), k with *
    simp [dotProduct.succ, subSequence.value]
    ring

def squaredNorm (u : RealSequence, n : Nat) : Real := dotProduct(u, u, n)

predicate Orthogonal (u v : Set, n : Nat) := dotProduct(u, v, n) = (0:R)

theorem Orthogonal.symmetric (u v : Set, n : Nat)
  requires sequenceU : RealSequence(u)
  requires sequenceV : RealSequence(v)
  requires orthogonal : Orthogonal(u, v, n)
  shows Orthogonal(v, u, n)
proof
  unfold Orthogonal at orthogonal
  unfold Orthogonal
  have symmetry := dotProduct.symmetric with sequenceU, sequenceV
  rwa symmetry at orthogonal

theorem squaredNorm.real (u : Set, n : Nat)
  requires sequence : RealSequence(u)
  shows Real(squaredNorm(u, n))
proof
  unfold squaredNorm
  exact dotProduct.real with sequence, sequence

theorem squaredNorm.succ (u : Set, n : Nat)
  requires sequence : RealSequence(u)
  shows squaredNorm(u, succ(n)) =
    (squaredNorm(u, n):R) + (u[n]:R) * (u[n]:R)
proof
  unfold squaredNorm
  exact dotProduct.succ with sequence, sequence

theorem squaredNorm.nonnegative (u : Set, n : Nat)
  requires sequence : RealSequence(u)
  shows (0:R) ≤ (squaredNorm(u, n):R)
proof
  induction n
  | zero =>
    unfold squaredNorm
    rw dotProduct.zero
    exact realLe.reflexive with realZero.real
  | succ(k) using ih =>
    have uk := RealSequence.value with sequence
    have prefix := squaredNorm.real with sequence
    have term := realMulValue.real with uk, uk
    have nonnegative := realMulValue.square_nonnegative with uk
    have added := realLe.add with realZero.real, prefix, realZero.real, term, ih, nonnegative
    have zero := realAddValue.zero_right with realZero.real
    rw zero at added
    unfold squaredNorm at added
    unfold squaredNorm
    rwa dotProduct.succ

theorem squaredNorm.scale (u : Set, C : Real, n : Nat)
  requires sequence : RealSequence(u)
  shows squaredNorm(mulSequence(constantSequence(C), u), n) =
    (C * C) * (squaredNorm(u, n):R)
proof
  have constants := constantSequence.sequence at C
  have scaled := mulSequence.sequence with constants, sequence
  have first := dotProduct.scale_left with sequence, scaled
  have second := dotProduct.scale_right with sequence, sequence
  simp [squaredNorm]
  ring

theorem squaredNorm.add (u v : Set, n : Nat)
  requires sequenceU : RealSequence(u)
  requires sequenceV : RealSequence(v)
  shows squaredNorm(addSequence(u, v), n) =
    (squaredNorm(u, n):R) +
      (2:R) * (dotProduct(u, v, n):R) +
        (squaredNorm(v, n):R)
proof
  have combined := addSequence.sequence with sequenceU, sequenceV
  have first := dotProduct.add_left with sequenceU, sequenceV, combined
  have second := dotProduct.add_right with sequenceU, sequenceU, sequenceV
  have third := dotProduct.add_right with sequenceV, sequenceU, sequenceV
  have symmetry := dotProduct.symmetric with sequenceV, sequenceU
  simp [squaredNorm]
  ring

theorem squaredNorm.sub (u v : Set, n : Nat)
  requires sequenceU : RealSequence(u)
  requires sequenceV : RealSequence(v)
  shows squaredNorm(subSequence(u, v), n) =
    ((squaredNorm(u, n):R) -
      (2:R) * (dotProduct(u, v, n):R)) +
        (squaredNorm(v, n):R)
proof
  induction n
  | zero =>
    unfold squaredNorm
    rw [dotProduct.zero, dotProduct.zero, dotProduct.zero, dotProduct.zero]
    ring
  | succ(k) using ih =>
    have difference := subSequence.sequence with sequenceU, sequenceV
    have uk := RealSequence.value with sequenceU
    have vk := RealSequence.value with sequenceV
    have differenceK := RealSequence.value with difference
    have value := subSequence.value with sequenceU, sequenceV, uk, vk
    rw [squaredNorm.succ, squaredNorm.succ, squaredNorm.succ, dotProduct.succ, value, ih]
    ring

theorem squaredNorm.parallelogram (u v : Set, n : Nat)
  requires sequenceU : RealSequence(u)
  requires sequenceV : RealSequence(v)
  shows (squaredNorm(addSequence(u, v), n):R) +
      (squaredNorm(subSequence(u, v), n):R) =
    (2:R) * (squaredNorm(u, n):R) +
      (2:R) * (squaredNorm(v, n):R)
proof
  have combined := addSequence.sequence with sequenceU, sequenceV
  have difference := subSequence.sequence with sequenceU, sequenceV
  have normU := squaredNorm.real with sequenceU
  have normV := squaredNorm.real with sequenceV
  have cross := dotProduct.real with sequenceU, sequenceV
  have plus := squaredNorm.add with sequenceU, sequenceV
  have minus := squaredNorm.sub with sequenceU, sequenceV
  have algebra : ∀ A ∈ ℝ, ∀ B ∈ ℝ, ∀ C ∈ ℝ,
      ((A:R) + B + C) + ((A - B) + C) =
        (2:R) * A + (2:R) * C
    intros A realA B realB C realC
    ring
  have memberNormU : squaredNorm(u, n) ∈ ℝ by exact normU
  have memberNormV : squaredNorm(v, n) ∈ ℝ by exact normV
  have memberCross : (2:R) * dotProduct(u, v, n) ∈ ℝ by typecheck
  have combinedAlgebra := algebra with squaredNorm(u, n), memberNormU,
    (2:R) * dotProduct(u, v, n), memberCross,
    squaredNorm(v, n), memberNormV
  calc
    (squaredNorm(addSequence(u, v), n):R) +
        (squaredNorm(subSequence(u, v), n):R) =
      ((squaredNorm(u, n):R) +
        (2:R) * (dotProduct(u, v, n):R) +
          (squaredNorm(v, n):R)) +
      (((squaredNorm(u, n):R) -
        (2:R) * (dotProduct(u, v, n):R)) +
          (squaredNorm(v, n):R)) := by cc
    _ = (2:R) * (squaredNorm(u, n):R) +
          (2:R) * (squaredNorm(v, n):R) := combinedAlgebra

theorem squaredNorm.pythagorean (u v : Set, n : Nat)
  requires sequenceU : RealSequence(u)
  requires sequenceV : RealSequence(v)
  requires orthogonal : Orthogonal(u, v, n)
  shows squaredNorm(addSequence(u, v), n) =
    (squaredNorm(u, n):R) + (squaredNorm(v, n):R)
proof
  have expansion := squaredNorm.add with sequenceU, sequenceV
  unfold Orthogonal at orthogonal
  rw orthogonal at expansion
  calc
    squaredNorm(addSequence(u, v), n) =
        (squaredNorm(u, n):R) + (2:R) * (0:R) +
          (squaredNorm(v, n):R) := expansion
    _ = (squaredNorm(u, n):R) + (squaredNorm(v, n):R) := by ring