-- Rotations about the z-axis without trigonometry: planeTurn(a, b) is the
-- matrix with rows (a, −b, 0), (b, a, 0), (0, 0, 1), turn(t) is
-- planeTurn((1 − t²)/(1 + t²), 2t/(1 + t²)), a rotation for every real t;
-- the plane turns multiply like complex numbers, t ↦ turn(t) is injective,
-- and powers of 3 × 3 matrices
-- 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.Geometry.Absorption.Turn
import Triples
import Real.Squares

-- Cancelling a nonzero real factor.
theorem realMulValue.right_cancel_of_nonzero (X Y D : Real)
  requires nonzero : D ≠ 0
  requires same : X * D = Y * D
  shows X = Y
proof
  have inverse := realInvValue.mul_right with nonzero
  have inverseReal := realInvValue.real at D
  have leftUnit : X * (D * realInvValue(D)) = X
    rw inverse
    ring
  have rightUnit : Y * (D * realInvValue(D)) = Y
    rw inverse
    ring
  calc
    X = X * (D * realInvValue(D)) := by rw leftUnit
    _ = (X * D) * realInvValue(D) := by ring
    _ = (Y * D) * realInvValue(D) := by rw same
    _ = Y * (D * realInvValue(D)) := by ring
    _ = Y := rightUnit

-- Powers of a 3 × 3 real matrix in the matrix ring.
def matrixPow3 (M : Set, n : Nat) : Set :=
  monoidPow(Ring.mul(matrixRing(realRing, 3)), Ring.one(matrixRing(realRing, 3)), M, n)

notation infixr 80 "^" => matrixPow3 for Matrix3, Nat

theorem matrixPow3.zero (M : Matrix3)
  shows M ^ 0 = Ring.one(matrixRing(realRing, 3))
proof
  unfold matrixPow3
  exact monoidPow.zero

theorem matrixPow3.succ (M : Matrix3) (n : Nat)
  shows M ^ (n + 1) = matrixMul3(M, M ^ n)
proof
  have step := monoidPow.succ at Ring.mul(matrixRing(realRing, 3)),
    Ring.one(matrixRing(realRing, 3)), M, n
  have succForm : n + 1 = succ(n) := add.right_one at n
  unfold matrixPow3
  rw matrixMul3.graph
  simp [call2]

instance matrixPow3.matrix3 (M : Matrix3) (n : Nat)
  shows Matrix3(M ^ n)
proof
  induction n
  | zero =>
    rw matrixPow3.zero
    exact Matrix3.one
  | succ(k) using ih =>
    have succForm : k + 1 = succ(k) := add.right_one at k
    have step := matrixPow3.succ at M, k with *
    rw succForm at step
    rw step
    exact Matrix3.mul with ih

-- The plane turn by the point (a, b) and the rational turn by t.
def planeTurn (a b : Real) : Set :=
  matrix3of(a, -b, (0:R), b, a, (0:R), (0:R), (0:R), (1:R))

theorem planeTurn.graph (a b : Real)
  shows planeTurn(a, b) = matrix3of(a, -b, (0:R), b, a, (0:R), (0:R), (0:R), (1:R))
proof
  simp [planeTurn]

instance planeTurn.matrix (a b : Real)
  shows planeTurn(a, b) : realRing ^ 3 × 3
proof
  rw planeTurn.graph
  exact matrix3of.matrix with *

theorem planeTurn.entry00 (a b : Real)
  shows planeTurn(a, b)[0, 0] = a
proof
  rw [planeTurn.graph, matrix3of.entry00]
  rfl

theorem planeTurn.entry10 (a b : Real)
  shows planeTurn(a, b)[1, 0] = b
proof
  rw [planeTurn.graph, matrix3of.entry10]
  rfl

-- Plane turns multiply as complex numbers do.
theorem planeTurn.mul (a b c d : Real)
  shows matrixMul3(planeTurn(a, b), planeTurn(c, d)) = planeTurn(a * c - b * d, a * d + b * c)
proof
  rw [planeTurn.graph, planeTurn.graph, planeTurn.graph, matrix3of.mul]
  apply matrix3of.congr
  . ring
  . ring
  . ring
  . ring
  . ring
  . ring
  . ring
  . ring
  . ring

theorem planeTurn.one
  shows planeTurn((1:R), (0:R)) = Ring.one(matrixRing(realRing, 3))
proof
  have negZero : -(0:R) = (0:R) := realNegValue.zero
  rw [planeTurn.graph, negZero]
  exact matrix3of.one

def turnCos (t : Real) : Real := ((1:R) - t * t) / ((1:R) + t * t)
def turnSin (t : Real) : Real := (t + t) / ((1:R) + t * t)
def turn (t : Real) : Set := planeTurn(turnCos(t), turnSin(t))

theorem turnCos.graph (t : Real)
  shows turnCos(t) = ((1:R) - t * t) / ((1:R) + t * t)
proof
  simp [turnCos]

theorem turnSin.graph (t : Real)
  shows turnSin(t) = (t + t) / ((1:R) + t * t)
proof
  simp [turnSin]

theorem turn.graph (t : Real)
  shows turn(t) = planeTurn(turnCos(t), turnSin(t))
proof
  simp [turn]

instance turn.matrix (t : Real)
  shows turn(t) : realRing ^ 3 × 3
proof
  rw turn.graph
  exact planeTurn.matrix with *

theorem turnDenominator.positive (t : Real)
  shows (0:R) < (1:R) + t * t
proof
  have square := realMulValue.square_nonnegative at t
  have one := realZero.lt_one
  linarith

theorem turnDenominator.nonzero (t : Real)
  shows (1:R) + t * t ≠ 0
proof
  have positive := turnDenominator.positive at t
  intro zero
  rw zero at positive
  linarith

theorem turnCos.mul_denominator (t : Real)
  shows turnCos(t) * ((1:R) + t * t) = (1:R) - t * t
proof
  have nonzero := turnDenominator.nonzero at t
  rw turnCos.graph
  exact realDiv.mul_right with nonzero

theorem turnSin.mul_denominator (t : Real)
  shows turnSin(t) * ((1:R) + t * t) = t + t
proof
  have nonzero := turnDenominator.nonzero at t
  rw turnSin.graph
  exact realDiv.mul_right with nonzero

-- The point (turnCos(t), turnSin(t)) lies on the unit circle.
theorem turn.circle (t : Real)
  shows turnCos(t) * turnCos(t) + turnSin(t) * turnSin(t) = 1
proof
  have nonzero := turnDenominator.nonzero at t
  have cosRule := turnCos.mul_denominator at t
  have sinRule := turnSin.mul_denominator at t
  have regroup : (turnCos(t) * turnCos(t) + turnSin(t) * turnSin(t)) *
      ((1:R) + t * t) * ((1:R) + t * t) =
      (turnCos(t) * ((1:R) + t * t)) * (turnCos(t) * ((1:R) + t * t)) +
      (turnSin(t) * ((1:R) + t * t)) * (turnSin(t) * ((1:R) + t * t))
    ring
  have identity : ((1:R) - t * t) * ((1:R) - t * t) + (t + t) * (t + t) =
      (1:R) * ((1:R) + t * t) * ((1:R) + t * t)
    ring
  have expanded : (turnCos(t) * turnCos(t) + turnSin(t) * turnSin(t)) *
      ((1:R) + t * t) * ((1:R) + t * t) = (1:R) * ((1:R) + t * t) * ((1:R) + t * t)
    rw [regroup, cosRule, sinRule, identity]
    rfl
  have once := realMulValue.right_cancel_of_nonzero with nonzero, expanded
  exact realMulValue.right_cancel_of_nonzero with nonzero, once

-- Recovering t from the turn: t = turnSin(t) / (1 + turnCos(t)).
theorem turnCos.one_add_mul_denominator (t : Real)
  shows ((1:R) + turnCos(t)) * ((1:R) + t * t) = (1:R) + (1:R)
proof
  have cosRule := turnCos.mul_denominator at t
  have expand : ((1:R) + turnCos(t)) * ((1:R) + t * t) =
      ((1:R) + t * t) + turnCos(t) * ((1:R) + t * t)
    ring
  rw [expand, cosRule]
  ring

theorem turnCos.one_add_nonzero (t : Real)
  shows (1:R) + turnCos(t) ≠ 0
proof
  have rule := turnCos.one_add_mul_denominator at t
  have one := realZero.lt_one
  intro zero
  rw zero at rule
  have vanish : (0:R) * ((1:R) + t * t) = 0 := realMulValue.zero_left at (1:R) + t * t
  rw vanish at rule
  linarith

theorem turn.recover (t : Real)
  shows turnSin(t) / ((1:R) + turnCos(t)) = t
proof
  have nonzero := turnDenominator.nonzero at t
  have oneAddNonzero := turnCos.one_add_nonzero at t
  have oneAddRule := turnCos.one_add_mul_denominator at t
  have sinRule := turnSin.mul_denominator at t
  have regroup : (t * ((1:R) + turnCos(t))) * ((1:R) + t * t) =
      t * (((1:R) + turnCos(t)) * ((1:R) + t * t))
    ring
  have same : (t * ((1:R) + turnCos(t))) * ((1:R) + t * t) =
      turnSin(t) * ((1:R) + t * t)
    rw [regroup, oneAddRule, sinRule]
    ring
  have product := realMulValue.right_cancel_of_nonzero with nonzero, same
  have selfQuotient := realDiv.self with oneAddNonzero
  rw [← product, realDiv.mul_numerator_left, selfQuotient]
  ring

theorem turn.entry00 (t : Real)
  shows turn(t)[0, 0] = turnCos(t)
proof
  rw turn.graph
  exact planeTurn.entry00 with *

theorem turn.entry10 (t : Real)
  shows turn(t)[1, 0] = turnSin(t)
proof
  rw turn.graph
  exact planeTurn.entry10 with *

-- The parametrization t ↦ turn(t) is injective: t is recovered from the
-- point (turnCos(t), turnSin(t)).
theorem turn.injective_of_coordinates (s t : Real)
  requires cosEqual : turnCos(s) = turnCos(t)
  requires sinEqual : turnSin(s) = turnSin(t)
  shows s = t
proof
  have recoverS := turn.recover at s
  have recoverT := turn.recover at t
  calc
    s = turnSin(s) / ((1:R) + turnCos(s)) := by rw recoverS
    _ = turnSin(t) / ((1:R) + turnCos(t)) := by rw [cosEqual, sinEqual]
    _ = t := recoverT

theorem turn.injective (s t : Real)
  requires same : turn(s) = turn(t)
  shows s = t
proof
  have cosS := turn.entry00 at s
  have cosT := turn.entry00 at t
  have sinS := turn.entry10 at s
  have sinT := turn.entry10 at t
  have cosEqual : turnCos(s) = turnCos(t)
    calc
      turnCos(s) = turn(s)[0, 0] := by rw cosS
      _ = turn(t)[0, 0] := by rw same
      _ = turnCos(t) := cosT
  have sinEqual : turnSin(s) = turnSin(t)
    calc
      turnSin(s) = turn(s)[1, 0] := by rw sinS
      _ = turn(t)[1, 0] := by rw same
      _ = turnSin(t) := sinT
  exact turn.injective_of_coordinates with cosEqual, sinEqual

-- The transpose of a 3 × 3 matrix, and orthogonality: M · Mᵀ = 1.
def transpose3 (M : Set) : Set :=
  matrix3of(M[0, 0], M[1, 0], M[2, 0], M[0, 1], M[1, 1], M[2, 1], M[0, 2], M[1, 2], M[2, 2])

theorem transpose3.graph (M : Set)
  shows transpose3(M) =
    matrix3of(M[0, 0], M[1, 0], M[2, 0], M[0, 1], M[1, 1], M[2, 1], M[0, 2], M[1, 2], M[2, 2])
proof
  simp [transpose3]

predicate Orthogonal3 (M : Set) :=
  Matrix(M, realRing, 3, 3) ∧
  matrixMul3(M, transpose3(M)) = Ring.one(matrixRing(realRing, 3))

theorem planeTurn.transpose (a b : Real)
  shows transpose3(planeTurn(a, b)) = planeTurn(a, -b)
proof
  rw [transpose3.graph, planeTurn.graph, planeTurn.graph]
  matrix3of_entries
  apply matrix3of.congr
  . rfl
  . ring
  . rfl
  . rfl
  . rfl
  . rfl
  . rfl
  . rfl
  . rfl

theorem planeTurn.orthogonal (a b : Real)
  requires circle : a * a + b * b = 1
  shows Orthogonal3(planeTurn(a, b))
proof
  have matrix := planeTurn.matrix at a, b
  have product := planeTurn.mul at a, b, a, -b
  have first : a * a - b * (-b) = (1:R)
    rw ← circle
    ring
  have second : a * (-b) + b * a = (0:R)
    ring
  unfold Orthogonal3
  constructor
  . exact matrix
  . rw [planeTurn.transpose, product, first, second]
    exact planeTurn.one

-- turn(t) is orthogonal.
theorem turn.orthogonal (t : Real)
  shows Orthogonal3(turn(t))
proof
  have circle := turn.circle at t
  rw turn.graph
  exact planeTurn.orthogonal with circle