Grace/Cumulus/Geometry/Absorption/Turn.g17

Open in Workbench Download source

Imports

Triples · Real.Squares

Declarations

Source

1-- Rotations about the z-axis without trigonometry: planeTurn(a, b) is the2-- matrix with rows (a, −b, 0), (b, a, 0), (0, 0, 1), turn(t) is3-- planeTurn((1 − t²)/(1 + t²), 2t/(1 + t²)), a rotation for every real t;4-- the plane turns multiply like complex numbers, t ↦ turn(t) is injective,5-- and powers of 3 × 3 matrices6-- Part of Cumulus, the Grace17 mathematics library.7-- Copyright (c) 2026 Grant Passmore.8-- Released under the Apache License 2.0. See LICENSE in the repository root.9-- Authors: Grant Passmore1011module Grace.Cumulus.Geometry.Absorption.Turn12import Triples13import Real.Squares1415-- Cancelling a nonzero real factor.16theorem realMulValue.right_cancel_of_nonzero (X Y D : Real)17  requires nonzero : D ≠ 018  requires same : X * D = Y * D19  shows X = Y20proof21  have inverse := realInvValue.mul_right with nonzero22  have inverseReal := realInvValue.real at D23  have leftUnit : X * (D * realInvValue(D)) = X24    rw inverse25    ring26  have rightUnit : Y * (D * realInvValue(D)) = Y27    rw inverse28    ring29  calc30    X = X * (D * realInvValue(D)) := by rw leftUnit31    _ = (X * D) * realInvValue(D) := by ring32    _ = (Y * D) * realInvValue(D) := by rw same33    _ = Y * (D * realInvValue(D)) := by ring34    _ = Y := rightUnit3536-- Powers of a 3 × 3 real matrix in the matrix ring.37def matrixPow3 (M : Set, n : Nat) : Set :=38  monoidPow(Ring.mul(matrixRing(realRing, 3)), Ring.one(matrixRing(realRing, 3)), M, n)3940notation infixr 80 "^" => matrixPow3 for Matrix3, Nat4142theorem matrixPow3.zero (M : Matrix3)43  shows M ^ 0 = Ring.one(matrixRing(realRing, 3))44proof45  unfold matrixPow346  exact monoidPow.zero4748theorem matrixPow3.succ (M : Matrix3) (n : Nat)49  shows M ^ (n + 1) = matrixMul3(M, M ^ n)50proof51  have step := monoidPow.succ at Ring.mul(matrixRing(realRing, 3)),52    Ring.one(matrixRing(realRing, 3)), M, n53  have succForm : n + 1 = succ(n) := add.right_one at n54  unfold matrixPow355  rw matrixMul3.graph56  simp [call2]5758instance matrixPow3.matrix3 (M : Matrix3) (n : Nat)59  shows Matrix3(M ^ n)60proof61  induction n62  | zero =>63    rw matrixPow3.zero64    exact Matrix3.one65  | succ(k) using ih =>66    have succForm : k + 1 = succ(k) := add.right_one at k67    have step := matrixPow3.succ at M, k with *68    rw succForm at step69    rw step70    exact Matrix3.mul with ih7172-- The plane turn by the point (a, b) and the rational turn by t.73def planeTurn (a b : Real) : Set :=74  matrix3of(a, -b, (0:R), b, a, (0:R), (0:R), (0:R), (1:R))7576theorem planeTurn.graph (a b : Real)77  shows planeTurn(a, b) = matrix3of(a, -b, (0:R), b, a, (0:R), (0:R), (0:R), (1:R))78proof79  simp [planeTurn]8081instance planeTurn.matrix (a b : Real)82  shows planeTurn(a, b) : realRing ^ 3 × 383proof84  rw planeTurn.graph85  exact matrix3of.matrix with *8687theorem planeTurn.entry00 (a b : Real)88  shows planeTurn(a, b)[0, 0] = a89proof90  rw [planeTurn.graph, matrix3of.entry00]91  rfl9293theorem planeTurn.entry10 (a b : Real)94  shows planeTurn(a, b)[1, 0] = b95proof96  rw [planeTurn.graph, matrix3of.entry10]97  rfl9899-- Plane turns multiply as complex numbers do.100theorem planeTurn.mul (a b c d : Real)101  shows matrixMul3(planeTurn(a, b), planeTurn(c, d)) = planeTurn(a * c - b * d, a * d + b * c)102proof103  rw [planeTurn.graph, planeTurn.graph, planeTurn.graph, matrix3of.mul]104  apply matrix3of.congr105  . ring106  . ring107  . ring108  . ring109  . ring110  . ring111  . ring112  . ring113  . ring114115theorem planeTurn.one116  shows planeTurn((1:R), (0:R)) = Ring.one(matrixRing(realRing, 3))117proof118  have negZero : -(0:R) = (0:R) := realNegValue.zero119  rw [planeTurn.graph, negZero]120  exact matrix3of.one121122def turnCos (t : Real) : Real := ((1:R) - t * t) / ((1:R) + t * t)123def turnSin (t : Real) : Real := (t + t) / ((1:R) + t * t)124def turn (t : Real) : Set := planeTurn(turnCos(t), turnSin(t))125126theorem turnCos.graph (t : Real)127  shows turnCos(t) = ((1:R) - t * t) / ((1:R) + t * t)128proof129  simp [turnCos]130131theorem turnSin.graph (t : Real)132  shows turnSin(t) = (t + t) / ((1:R) + t * t)133proof134  simp [turnSin]135136theorem turn.graph (t : Real)137  shows turn(t) = planeTurn(turnCos(t), turnSin(t))138proof139  simp [turn]140141instance turn.matrix (t : Real)142  shows turn(t) : realRing ^ 3 × 3143proof144  rw turn.graph145  exact planeTurn.matrix with *146147theorem turnDenominator.positive (t : Real)148  shows (0:R) < (1:R) + t * t149proof150  have square := realMulValue.square_nonnegative at t151  have one := realZero.lt_one152  linarith153154theorem turnDenominator.nonzero (t : Real)155  shows (1:R) + t * t ≠ 0156proof157  have positive := turnDenominator.positive at t158  intro zero159  rw zero at positive160  linarith161162theorem turnCos.mul_denominator (t : Real)163  shows turnCos(t) * ((1:R) + t * t) = (1:R) - t * t164proof165  have nonzero := turnDenominator.nonzero at t166  rw turnCos.graph167  exact realDiv.mul_right with nonzero168169theorem turnSin.mul_denominator (t : Real)170  shows turnSin(t) * ((1:R) + t * t) = t + t171proof172  have nonzero := turnDenominator.nonzero at t173  rw turnSin.graph174  exact realDiv.mul_right with nonzero175176-- The point (turnCos(t), turnSin(t)) lies on the unit circle.177theorem turn.circle (t : Real)178  shows turnCos(t) * turnCos(t) + turnSin(t) * turnSin(t) = 1179proof180  have nonzero := turnDenominator.nonzero at t181  have cosRule := turnCos.mul_denominator at t182  have sinRule := turnSin.mul_denominator at t183  have regroup : (turnCos(t) * turnCos(t) + turnSin(t) * turnSin(t)) *184      ((1:R) + t * t) * ((1:R) + t * t) =185      (turnCos(t) * ((1:R) + t * t)) * (turnCos(t) * ((1:R) + t * t)) +186      (turnSin(t) * ((1:R) + t * t)) * (turnSin(t) * ((1:R) + t * t))187    ring188  have identity : ((1:R) - t * t) * ((1:R) - t * t) + (t + t) * (t + t) =189      (1:R) * ((1:R) + t * t) * ((1:R) + t * t)190    ring191  have expanded : (turnCos(t) * turnCos(t) + turnSin(t) * turnSin(t)) *192      ((1:R) + t * t) * ((1:R) + t * t) = (1:R) * ((1:R) + t * t) * ((1:R) + t * t)193    rw [regroup, cosRule, sinRule, identity]194    rfl195  have once := realMulValue.right_cancel_of_nonzero with nonzero, expanded196  exact realMulValue.right_cancel_of_nonzero with nonzero, once197198-- Recovering t from the turn: t = turnSin(t) / (1 + turnCos(t)).199theorem turnCos.one_add_mul_denominator (t : Real)200  shows ((1:R) + turnCos(t)) * ((1:R) + t * t) = (1:R) + (1:R)201proof202  have cosRule := turnCos.mul_denominator at t203  have expand : ((1:R) + turnCos(t)) * ((1:R) + t * t) =204      ((1:R) + t * t) + turnCos(t) * ((1:R) + t * t)205    ring206  rw [expand, cosRule]207  ring208209theorem turnCos.one_add_nonzero (t : Real)210  shows (1:R) + turnCos(t) ≠ 0211proof212  have rule := turnCos.one_add_mul_denominator at t213  have one := realZero.lt_one214  intro zero215  rw zero at rule216  have vanish : (0:R) * ((1:R) + t * t) = 0 := realMulValue.zero_left at (1:R) + t * t217  rw vanish at rule218  linarith219220theorem turn.recover (t : Real)221  shows turnSin(t) / ((1:R) + turnCos(t)) = t222proof223  have nonzero := turnDenominator.nonzero at t224  have oneAddNonzero := turnCos.one_add_nonzero at t225  have oneAddRule := turnCos.one_add_mul_denominator at t226  have sinRule := turnSin.mul_denominator at t227  have regroup : (t * ((1:R) + turnCos(t))) * ((1:R) + t * t) =228      t * (((1:R) + turnCos(t)) * ((1:R) + t * t))229    ring230  have same : (t * ((1:R) + turnCos(t))) * ((1:R) + t * t) =231      turnSin(t) * ((1:R) + t * t)232    rw [regroup, oneAddRule, sinRule]233    ring234  have product := realMulValue.right_cancel_of_nonzero with nonzero, same235  have selfQuotient := realDiv.self with oneAddNonzero236  rw [← product, realDiv.mul_numerator_left, selfQuotient]237  ring238239theorem turn.entry00 (t : Real)240  shows turn(t)[0, 0] = turnCos(t)241proof242  rw turn.graph243  exact planeTurn.entry00 with *244245theorem turn.entry10 (t : Real)246  shows turn(t)[1, 0] = turnSin(t)247proof248  rw turn.graph249  exact planeTurn.entry10 with *250251-- The parametrization t ↦ turn(t) is injective: t is recovered from the252-- point (turnCos(t), turnSin(t)).253theorem turn.injective_of_coordinates (s t : Real)254  requires cosEqual : turnCos(s) = turnCos(t)255  requires sinEqual : turnSin(s) = turnSin(t)256  shows s = t257proof258  have recoverS := turn.recover at s259  have recoverT := turn.recover at t260  calc261    s = turnSin(s) / ((1:R) + turnCos(s)) := by rw recoverS262    _ = turnSin(t) / ((1:R) + turnCos(t)) := by rw [cosEqual, sinEqual]263    _ = t := recoverT264265theorem turn.injective (s t : Real)266  requires same : turn(s) = turn(t)267  shows s = t268proof269  have cosS := turn.entry00 at s270  have cosT := turn.entry00 at t271  have sinS := turn.entry10 at s272  have sinT := turn.entry10 at t273  have cosEqual : turnCos(s) = turnCos(t)274    calc275      turnCos(s) = turn(s)[0, 0] := by rw cosS276      _ = turn(t)[0, 0] := by rw same277      _ = turnCos(t) := cosT278  have sinEqual : turnSin(s) = turnSin(t)279    calc280      turnSin(s) = turn(s)[1, 0] := by rw sinS281      _ = turn(t)[1, 0] := by rw same282      _ = turnSin(t) := sinT283  exact turn.injective_of_coordinates with cosEqual, sinEqual284285-- The transpose of a 3 × 3 matrix, and orthogonality: M · Mᵀ = 1.286def transpose3 (M : Set) : Set :=287  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])288289theorem transpose3.graph (M : Set)290  shows transpose3(M) =291    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])292proof293  simp [transpose3]294295predicate Orthogonal3 (M : Set) :=296  Matrix(M, realRing, 3, 3) ∧297  matrixMul3(M, transpose3(M)) = Ring.one(matrixRing(realRing, 3))298299theorem planeTurn.transpose (a b : Real)300  shows transpose3(planeTurn(a, b)) = planeTurn(a, -b)301proof302  rw [transpose3.graph, planeTurn.graph, planeTurn.graph]303  matrix3of_entries304  apply matrix3of.congr305  . rfl306  . ring307  . rfl308  . rfl309  . rfl310  . rfl311  . rfl312  . rfl313  . rfl314315theorem planeTurn.orthogonal (a b : Real)316  requires circle : a * a + b * b = 1317  shows Orthogonal3(planeTurn(a, b))318proof319  have matrix := planeTurn.matrix at a, b320  have product := planeTurn.mul at a, b, a, -b321  have first : a * a - b * (-b) = (1:R)322    rw ← circle323    ring324  have second : a * (-b) + b * a = (0:R)325    ring326  unfold Orthogonal3327  constructor328  . exact matrix329  . rw [planeTurn.transpose, product, first, second]330    exact planeTurn.one331332-- turn(t) is orthogonal.333theorem turn.orthogonal (t : Real)334  shows Orthogonal3(turn(t))335proof336  have circle := turn.circle at t337  rw turn.graph338  exact planeTurn.orthogonal with circle339