From 02dc877361f05714ba269a1fd9ba495bf23c20f0 Mon Sep 17 00:00:00 2001 From: Mike Boyle Date: Sat, 20 Jun 2026 00:45:31 -0400 Subject: [PATCH 1/2] Add KAN decomposition function --- src/Lorentz.jl | 19 +++++++++++++++++++ src/Quaternionic.jl | 2 +- 2 files changed, 20 insertions(+), 1 deletion(-) diff --git a/src/Lorentz.jl b/src/Lorentz.jl index bc968e9..7557423 100644 --- a/src/Lorentz.jl +++ b/src/Lorentz.jl @@ -395,3 +395,22 @@ function vR(Λ::Lorentz{T}) where {T<:Real} v⃗ = v̂sinhη╱2 * (2coshη╱2 / (2coshη╱2^2 - 1)) return v⃗, R end + +# --------------------------------------------------------------------------- +# Iwasawa's KAN decomposition +# --------------------------------------------------------------------------- + +function KAN(Λ::Lorentz{T}) where {T<:Real} + 𝐭𝐳 = im * 𝐤 + u₊ = (1 + 𝐭𝐳) / T(2) + ℂℜΛu₊ = ℂreal(Λ * u₊) + Rₖ = rotor(ℂℜΛu₊) + RₐRₙ = conj(Rₖ) * Λ + coshφ╱2, sinhφ╱2 = real(RₐRₙ.w), imag(RₐRₙ.z) + Rₐ = coshφ╱2 + sinhφ╱2 * 𝐭𝐳 + Rₙraw = conj(Rₐ) * RₐRₙ + ξˣ╱sqrt2 = imag(Rₙraw.x) - real(Rₙraw.y) + ξʸ╱sqrt2 = imag(Rₙraw.y) + real(Rₙraw.x) + Rₙ = Lorentz{T}(1, (im*ξˣ╱sqrt2+ξʸ╱sqrt2)/2, (im*ξʸ╱sqrt2-ξˣ╱sqrt2)/2, 0) + return Rₖ, Rₐ, Rₙ +end diff --git a/src/Quaternionic.jl b/src/Quaternionic.jl index 677c80b..06bae6e 100644 --- a/src/Quaternionic.jl +++ b/src/Quaternionic.jl @@ -15,7 +15,7 @@ export QuatVec, quatvec, QuatVecF64, QuatVecF32, QuatVecF16 export components, basetype public value, iszerovalue public ℂconj, ℂreal, ℂimag, ℂreim -public RB, BR, Rv, vR +public RB, BR, Rv, vR, KAN export (⋅), (×), (×̂), normalize, norm export abs2vec, absvec export from_float_array, to_float_array, From 31ca9efd5f9d56b8db234087bc08cab2a0f06bda Mon Sep 17 00:00:00 2001 From: Mike Boyle Date: Sun, 21 Jun 2026 12:35:15 -0400 Subject: [PATCH 2/2] Add KAN decomposition, docs, and tests --- docs/src/spacetime_algebra.md | 188 ++++++++++++++++++++++++++++-- src/Lorentz.jl | 45 +++++-- test/lorentz_decomposition.jl | 213 +++++++++++++++++++++++++++++++++- 3 files changed, 425 insertions(+), 21 deletions(-) diff --git a/docs/src/spacetime_algebra.md b/docs/src/spacetime_algebra.md index 0572725..5fd99ed 100644 --- a/docs/src/spacetime_algebra.md +++ b/docs/src/spacetime_algebra.md @@ -27,13 +27,13 @@ this more general algebra. One crucial new feature is the *spinor norm*: for complex quaternions, the GA reverse gives ``𝐐𝐐̃ = w² + x² + y² + z²``, which is a -*"complex"* scalar, not the Euclidean norm ``|w|² + |x|² + |y|² + -|z|²`` which will always be real. The surprising point is that the -normalization condition ``𝐐𝐐̃ = 1`` is now *two real* conditions, -because this says that the "imaginary" part is zero. This reduces the -eight real degrees of freedom in a complex quaternion down to six, -which is the correct number for a Lorentz transformation in four -dimensions. +*complex* scalar, distinct from the Euclidean norm ``|w|² + |x|² + +|y|² + |z|²`` which will always be real. The surprising point is that +the normalization condition ``𝐐𝐐̃ = 1`` is now *two real* +conditions, because this says that the "imaginary" part is zero. This +reduces the eight real degrees of freedom in a complex quaternion down +to six, which is the correct number for a Lorentz transformation in +four dimensions. ## The spacetime algebra @@ -197,8 +197,6 @@ which is not equal to ``1`` in general. Phase factors therefore do *not* belong to the Lorentz rotor group, even though they are unit elements under the Euclidean norm (``|\exp(i\varphi)| = 1``). - - !!! note "Real vs. complex quaternions" For real quaternions, ``𝐐\,\widetilde{𝐐} = \sum r_i^2 = \sum |r_i|^2``, @@ -247,11 +245,177 @@ Lorentz transformations. and confirming that a positive rapidity boosts ``t \to \cosh\varphi\,t + \sinh\varphi\,x`` and ``x \to \sinh\varphi\,t + \cosh\varphi\,x``. + +## [Iwasawa's ``KA\,N`` decomposition](@id iwasawa-kan) + +The [Iwasawa +decomposition](https://en.wikipedia.org/wiki/Iwasawa_decomposition) is +a factorization of the Lorentz group or its double cover +``\mathrm{Spin}^+(3,1)`` (or any semisimple Lie group) into three +subgroups: + +```math +G = {KA\,N}. +``` + +This decomposition is most useful when we have a preferred time axis +and preferred spatial direction. We use those to pick out boosts +along that spatial direction, and null rotations that fix the +corresponding null vector. Specifically, we have the following +subgroups: + +- **``K``** is the maximal compact subgroup, which in the case of the + Lorentz group is isomorphic to the rotation group + ``\mathrm{SO}(3)``, and in the case of ``\mathrm{Spin}(3,1)`` is + isomorphic to ``\mathrm{Spin}(3)``. +- **``A``** is the abelian subgroup of boosts in a fixed direction. + We take this direction to be the z-axis, so that we have ``A = + \left\{ \exp\left[\tfrac{φₐ}{2} \, 𝐭𝐳 \right] \mid φₐ ∈ ℝ + \right\}``. +- **``N``** is the nilpotent subgroup of null rotations, which + consists of all Lorentz transformations that can be represented as + null rotations about a fixed null vector. We take this vector to be + the null vector ``\boldsymbol{ℓ}``, so that we have ``N = \left\{ + \exp\left[\tfrac{1}{2} \boldsymbol{ℓ ξ}\right] \mid \boldsymbol{ξ} = + ξˣ𝐱 + ξʸ𝐲 \right\}``, where ``\boldsymbol{ℓ} = (𝐭+𝐳)/√2``. + +Note that ``φₐ`` here is *not* the rapidity of the overall boost if we +factor a transformation just as a boost and a rotation. Rather, it is +the rapidity of the particular boost picked out as ``A`` by this +construction. ``N`` also contributes to the overall boost, so the +overall rapidity is not just ``φₐ``. + +One nice feature of this decomposition is that we can compute the +factorization of a given Lorentz transformation ``Λ`` fairly simply. +We first introduce the idempotent + +```math +u₊ = \frac{1}{2} (1 + 𝐭𝐳). +``` + +It is easy to verify ``u₊𝐭𝐳 = 𝐭𝐳u₊ = u₊`` and hence the idempotent +identity: ``u₊² = u₊``. Next, we write + +```math +Λ = Rₖ Rₐ Rₙ, +``` + +where ``Rₖ ∈ K``, ``Rₐ ∈ A``, and ``Rₙ ∈ N``. Starting with the +rightmost factor, we have + +```math +Rₙ = \exp\left[\tfrac{1}{2} \boldsymbol{ℓ ξ}\right] = 1 + \frac{1}{2} \boldsymbol{ℓ ξ}, +``` + +because ``\boldsymbol{ℓ}`` and ``\boldsymbol{ξ}`` anticommute, and +because ``\boldsymbol{ℓ}² = 0``, their product ``\boldsymbol{ℓ ξ}`` +itself is also nilpotent. (Remember, that's what the "N" stands for.) +We can simply compute the product ``Rₙ u₊`` by using the fact that +``\boldsymbol{ξ} u₊ = u₊ \boldsymbol{ξ}`` and expanding terms in our +basis, then find + +```math +Rₙ u₊ %&= \left(1 + \frac{1}{2} \boldsymbol{ℓ ξ}\right) u₊, \\ +% &= u₊ + \frac{1}{2} \boldsymbol{ℓ ξ} u₊, \\ +% &= u₊ + \frac{1}{2} \boldsymbol{ℓ} u₊ \boldsymbol{ξ}, \\ +% &= u₊ + \frac{1}{4\sqrt{2}} (𝐭+𝐳) (1+𝐭𝐳) \boldsymbol{ξ}, \\ +% &= u₊ + \frac{1}{4\sqrt{2}} (𝐭+𝐳+𝐭𝐭𝐳+𝐳𝐭𝐳) \boldsymbol{ξ}, \\ +% &= u₊ + \frac{1}{4\sqrt{2}} (𝐭+𝐳-𝐳-𝐭) \boldsymbol{ξ}, \\ += u₊. +``` + +Now, with + +```math +Rₐ = \exp\left[\tfrac{φₐ}{2} \, 𝐭𝐳 \right] += \cosh\left(\frac{φₐ}{2}\right) + \sinh\left(\frac{φₐ}{2}\right) 𝐭𝐳, +``` + +and the fact that ``𝐭𝐳u₊ = u₊``, we have + +```math +\begin{aligned} +Rₐ u₊ +&= \left[\cosh\left(\frac{φₐ}{2}\right) + \sinh\left(\frac{φₐ}{2}\right)\right] u₊ +&= e^{φₐ/2} u₊. +\end{aligned} +``` + +Putting these together, we have + +```math +Λ u₊ = e^{φₐ/2} Rₖ u₊. +``` + +Recall that ``e^{φₐ/2}`` is a strictly positive real number, and +``Rₖ`` is a pure rotation so it is "ℂ-real" — meaning that it is a +linear combination of only basis elements that do not have a factor of +``𝐭`` — whereas ``u₊`` is just a combination of 1 and a "ℂ-imaginary" +part. Therefore, we can take the real part of this expression to find + +```math +ℂ\Re\{Λ u₊\} = \frac{e^{φₐ/2}}{2} Rₖ. +``` + +This is just a strictly positive real number times ``Rₖ``, which we +know must have unit magnitude, so we can find ``Rₖ`` by normalizing +this real part: + +```math +Rₖ = \mathrm{normalize}\left(ℂ\Re\{Λ u₊\}\right). +``` + +We can continue solving for the remaining factors ``Rₐ`` and ``Rₙ``. +One option would be to take the norm of the expression above and solve +for ``φₐ``, but it requires no transcendental functions and is more +numerically stable to extract the components directly from the +remaining product. We can evaluate ``RₐRₙ = Rₖ^{-1} Λ``. But now, +let us examine the left-hand side analytically: + +```math +Rₐ Rₙ = \left[\cosh\left(\frac{φₐ}{2}\right) + \sinh\left(\frac{φₐ}{2}\right) 𝐭𝐳\right] + \left[1 + \frac{1}{2} \boldsymbol{ℓ ξ}\right] +``` + +The ``\boldsymbol{ℓ ξ}`` term is entirely within the span of ``\{𝐭𝐱, +𝐭𝐲, 𝐳𝐱, 𝐳𝐲\}``; therefore, it is not hard to show that the +scalar and ``𝐭𝐳`` terms in this product equal ``Rₐ``. As usual, we +write ``⟨·⟩₀`` to extract the scalar (grade‑0) part, and we have + +```math +Rₐ = \left⟨ Rₖ^{-1} Λ \right⟩₀ + \left⟨ Rₖ^{-1} Λ\, 𝐭𝐳 \right⟩₀\, 𝐭𝐳. +``` + +Finally, we can find the remaining factor by rearranging the original +factorization: + +```math +Rₙ = Rₐ^{-1} Rₖ^{-1} Λ. +``` + +This decomposition is implemented in the [`Quaternionic.KAN`](@ref) +function. + ## API reference -```@autodocs -Modules = [Quaternionic] -Pages = ["Lorentz.jl"] +```@meta +CurrentModule = Quaternionic +``` + +```@docs +Lorentz +Lorentz(::AbstractVector) +Boost +ga_components +RB +BR +Rv +vR +KAN +ℂreal +ℂimag +ℂreim +ℂconj ``` ## Further reading diff --git a/src/Lorentz.jl b/src/Lorentz.jl index 7557423..c20c360 100644 --- a/src/Lorentz.jl +++ b/src/Lorentz.jl @@ -268,7 +268,7 @@ This is distinct from the GA "real" part, given by `real(Λ)`, which is construc reverse, rather than the complex-conjugate. That will just be the scalar part of `Λ`, but may be complex; this function returns a full quaternion, with all components real. """ -ℂreal(Λ::Lorentz{T}) where {T<:Real} = Quaternion{T}( +ℂreal(Λ::AbstractQuaternion{Complex{T}}) where {T<:Real} = Quaternion{T}( real(Λ[1]), real(Λ[2]), real(Λ[3]), real(Λ[4]) ) @@ -285,7 +285,7 @@ the reverse, rather than the complex-conjugate. That will just be the quaternio part of `Λ`, but may be complex; this function returns a full quaternion, with all components real. """ -ℂimag(Λ::Lorentz{T}) where {T<:Real} = Quaternion{T}( +ℂimag(Λ::AbstractQuaternion{Complex{T}}) where {T<:Real} = Quaternion{T}( imag(Λ[1]), imag(Λ[2]), imag(Λ[3]), imag(Λ[4]) ) @@ -323,7 +323,8 @@ find `cosh(η/2)` as the magnitude of the complex-real part of `Λ`. Then, we f `sinh(η/2)*n̂` simply by multiplying the complex-imaginary part of `Λ` by `inv(R)`, and reconstruct `B` simply by adding the result to `cosh(η/2)`. -See also [`BR`](@ref). +See also [`BR`](@ref), [`vR`](@ref), and [`Rv`](@ref) for similar forms, and [`KAN`](@ref) +for the Iwasawa decomposition. """ function RB(Λ::Lorentz{T}) where {T<:Real} ℜΛ, ℑΛ = ℂreim(Λ) @@ -339,6 +340,9 @@ Polar decomposition `Λ = B * R`: pure boost `B` followed by pure rotation `R`. The algorithm here is the same as for [`RB`](@ref), but with the order of multiplication by `inv(R)` reversed. + +See also [`vR`](@ref) and [`Rv`](@ref) for similar forms, and [`KAN`](@ref) for the Iwasawa +decomposition. """ function BR(Λ::Lorentz{T}) where {T<:Real} ℜΛ, ℑΛ = ℂreim(Λ) @@ -370,7 +374,8 @@ where the last equality uses the double-angle identity ``\cosh(η) = 2\cosh^2(η The last form is made up of entirely of the scalar component of `B`, does not involve any cancellation or division by small numbers, and avoids computing the norm of the vector part. -See also [`vR`](@ref), [`RB`](@ref), and [`BR`](@ref). +See also [`vR`](@ref), [`RB`](@ref), and [`BR`](@ref) for similar forms, and [`KAN`](@ref) +for the Iwasawa decomposition. """ function Rv(Λ::Lorentz{T}) where {T<:Real} R, B = RB(Λ) @@ -386,7 +391,8 @@ end Return the boost velocity `v⃗` and pure rotation `R` such that `Λ = Boost(η, v̂) * R`. See [`Rv`](@ref) for more details. -See also [`BR`](@ref) and [`RB`](@ref). +See also [`BR`](@ref) and [`RB`](@ref) for similar forms, and [`KAN`](@ref) for the Iwasawa +decomposition. """ function vR(Λ::Lorentz{T}) where {T<:Real} B, R = BR(Λ) @@ -400,17 +406,40 @@ end # Iwasawa's KAN decomposition # --------------------------------------------------------------------------- +@doc raw""" + KAN(Λ::Lorentz{T}) → (Rₖ, Rₐ, Rₙ) + +Compute the Iwasawa ``KAN`` decomposition of the Lorentz transformation `Λ`, returning the +three factors such that `Λ = Rₖ * Rₐ * Rₙ`, where + +- `Rₖ ∈ K` is a pure rotation (the maximal compact factor), +- `Rₐ ∈ A` is a boost along the ``𝐳`` axis, `Rₐ = cosh(φₐ/2) + sinh(φₐ/2) 𝐭𝐳`, and +- `Rₙ ∈ N` is a null rotation fixing the null vector ``ℓ = (𝐭+𝐳)/√2``. + +See the [`KAN` decomposition section](@ref iwasawa-kan) of the spacetime-algebra +documentation for a full discussion and derivation of the method used in this function. + +See also [`BR`](@ref), [`RB`](@ref), [`vR`](@ref), and [`Rv`](@ref) for polar +decompositions. +""" function KAN(Λ::Lorentz{T}) where {T<:Real} 𝐭𝐳 = im * 𝐤 u₊ = (1 + 𝐭𝐳) / T(2) + # We readily find the rotation factor by projecting with the idempotent, taking the pure + # spatial part, and normalizing. ℂℜΛu₊ = ℂreal(Λ * u₊) Rₖ = rotor(ℂℜΛu₊) + # We avoid transcendental functions `ln`, `sinh`, and `cosh` by using the structure of + # the following product to read off just the relevant components. RₐRₙ = conj(Rₖ) * Λ coshφ╱2, sinhφ╱2 = real(RₐRₙ.w), imag(RₐRₙ.z) Rₐ = coshφ╱2 + sinhφ╱2 * 𝐭𝐳 + # Rather than returning the raw `Rₐ⁻¹ * RₐRₙ`, which may have accumulated roundoff + # errors, we project that result onto the null-rotation sector by averaging 𝐢 and 𝐣 + # components, and setting 𝟏 and 𝐤 components to their exact values. Rₙraw = conj(Rₐ) * RₐRₙ - ξˣ╱sqrt2 = imag(Rₙraw.x) - real(Rₙraw.y) - ξʸ╱sqrt2 = imag(Rₙraw.y) + real(Rₙraw.x) - Rₙ = Lorentz{T}(1, (im*ξˣ╱sqrt2+ξʸ╱sqrt2)/2, (im*ξʸ╱sqrt2-ξˣ╱sqrt2)/2, 0) + ξˣ╱2sqrt2 = (imag(Rₙraw.x) - real(Rₙraw.y)) / 2 + ξʸ╱2sqrt2 = (imag(Rₙraw.y) + real(Rₙraw.x)) / 2 + Rₙ = Lorentz{T}(1, ξʸ╱2sqrt2 + im * ξˣ╱2sqrt2, -ξˣ╱2sqrt2 + im * ξʸ╱2sqrt2, 0) return Rₖ, Rₐ, Rₙ end diff --git a/test/lorentz_decomposition.jl b/test/lorentz_decomposition.jl index 874b837..bb15b88 100644 --- a/test/lorentz_decomposition.jl +++ b/test/lorentz_decomposition.jl @@ -10,7 +10,7 @@ using Random: Xoshiro using Quaternionic using DoubleFloats: Double64 - import Quaternionic: RB, BR, ℂconj, ℂreal, ℂimag, ℂreim + import Quaternionic: RB, BR, KAN, ℂconj, ℂreal, ℂimag, ℂreim const n = 20 @@ -57,6 +57,35 @@ abs(imag(cs[1])) < atol && all(c -> abs(real(c)) < atol, cs[2:4]) end + + # ── KAN-specific helpers ─────────────────────────────────────────────────── + + # A null rotation in N (fixing ℓ = (𝐭+𝐳)/√2), built from screw parameters ξˣ, ξʸ + # using the same encoding as `KAN` returns for its Rₙ factor. + function null_rotation(::Type{T}, ξˣ, ξʸ) where {T} + a, b = T(ξˣ)/sqrt(T(2)), T(ξʸ)/sqrt(T(2)) + Lorentz{T}(1, (im*a + b)/2, (im*b - a)/2, 0) + end + + # A boost along the z-axis with rapidity φ — an element of A. + z_boost(::Type{T}, φ) where {T} = Boost(T(φ), QuatVec{T}(0, 0, 1)) + + # The null vector ℓ = (𝐭+𝐳)/√2 fixed by N, as a Minkowski 4-vector. + ℓ_vec(::Type{T}) where {T} = T[1, 0, 0, 1] / sqrt(T(2)) + + # Structural predicate for A: pure z-boost (w real, x=y=0, z pure imaginary). + # Accepts any AbstractQuaternion (KAN returns Rₐ as a bare `Quaternion`). + function is_z_boost(R; atol) + w, x, y, z = components(R) + abs(imag(w)) < atol && abs(x) < atol && abs(y) < atol && abs(real(z)) < atol + end + + # Structural predicate for N: null rotation (w=1, z=0, Re(x)=Im(y), Im(x)=−Re(y)). + function is_null_rotation(R; atol) + w, x, y, z = components(R) + abs(w - 1) < atol && abs(z) < atol && + abs(real(x) - imag(y)) < atol && abs(imag(x) + real(y)) < atol + end end # ── Special cases ───────────────────────────────────────────────────────────── @@ -381,3 +410,185 @@ end @test vr_error(Double64) ≤ 10eps(Double64) @test vr_error(BigFloat) ≤ 10eps(BigFloat) end + +# ── KAN (Iwasawa) decomposition ──────────────────────────────────────────────── +# +# KAN(Λ) → (Rₖ, Rₐ, Rₙ) with Λ = Rₖ·Rₐ·Rₙ, where Rₖ ∈ K is a rotation, Rₐ ∈ A is a +# boost along 𝐳, and Rₙ ∈ N is a null rotation fixing ℓ = (𝐭+𝐳)/√2. Errors scale +# with the overall transformation magnitude, captured by `maximum(abs, components(Λ))`. + +@testitem "KAN: reconstruction round-trip" tags=[:validation, :fast] setup=[LorentzDecompData] begin + import Quaternionic: KAN, components + using .LorentzDecompData: mixed, rand_rotations, rand_boosts, same_Λ + using DoubleFloats: Double64 + for T ∈ (Float32, Float64, Double64) + for Λ ∈ [mixed(T); rand_rotations(T); rand_boosts(T)] + Rₖ, Rₐ, Rₙ = KAN(Λ) + ε = 1024eps(T) * maximum(abs, components(Λ)) + @test same_Λ(Lorentz(Rₖ) * Lorentz(Rₐ) * Lorentz(Rₙ), Λ; atol=ε) + end + end +end + +@testitem "KAN: structural membership (K, A, N)" tags=[:validation, :fast] setup=[LorentzDecompData] begin + import Quaternionic: KAN, components + using .LorentzDecompData: mixed, rand_rotations, rand_boosts, is_real_rotor, is_z_boost, is_null_rotation + using DoubleFloats: Double64 + for T ∈ (Float32, Float64, Double64) + for Λ ∈ [mixed(T); rand_rotations(T); rand_boosts(T)] + Rₖ, Rₐ, Rₙ = KAN(Λ) + ε = 1024eps(T) * maximum(abs, components(Λ)) + @test is_real_rotor(Rₖ; atol=ε) # Rₖ ∈ K: pure rotation (no boost part) + @test is_z_boost(Rₐ; atol=ε) # Rₐ ∈ A: boost along 𝐳 only + @test is_null_rotation(Rₙ; atol=ε) # Rₙ ∈ N: null rotation, w=1, z=0 + end + end +end + +@testitem "KAN: fixed-point oracles" tags=[:validation, :fast] setup=[LorentzDecompData] begin + import Quaternionic: KAN, components + using .LorentzDecompData: mixed, rand_rotations, rand_boosts, ℓ_vec + using DoubleFloats: Double64 + # K fixes the time axis 𝐭; A (boost along 𝐳) fixes the transverse axes 𝐱, 𝐲; + # N (null rotation) fixes the null vector ℓ = (𝐭+𝐳)/√2. + for T ∈ (Float32, Float64, Double64) + 𝐭, 𝐱, 𝐲 = T[1,0,0,0], T[0,1,0,0], T[0,0,1,0] + ℓ = ℓ_vec(T) + for Λ ∈ [mixed(T); rand_rotations(T); rand_boosts(T)] + Rₖ, Rₐ, Rₙ = KAN(Λ) + ε = 1024eps(T) * maximum(abs, components(Λ)) + @test isapprox(Lorentz(Rₖ)(𝐭), 𝐭; atol=ε) + @test isapprox(Lorentz(Rₐ)(𝐱), 𝐱; atol=ε) + @test isapprox(Lorentz(Rₐ)(𝐲), 𝐲; atol=ε) + @test isapprox(Lorentz(Rₙ)(ℓ), ℓ; atol=ε) + end + end +end + +@testitem "KAN: special cases" tags=[:unit, :fast, :validation] setup=[LorentzDecompData] begin + import Quaternionic: KAN + using .LorentzDecompData: rand_rotations, null_rotation, z_boost, same_Λ + using DoubleFloats: Double64 + for T ∈ (Float32, Float64, Double64) + ε = 64eps(T) + 𝟙 = one(Lorentz{T}) + + # Identity and minus-identity → all three factors trivial + for Λ ∈ (𝟙, -𝟙) + Rₖ, Rₐ, Rₙ = KAN(Λ) + @test same_Λ(Lorentz(Rₖ), 𝟙; atol=ε) + @test same_Λ(Lorentz(Rₐ), 𝟙; atol=ε) + @test same_Λ(Lorentz(Rₙ), 𝟙; atol=ε) + @test same_Λ(Lorentz(Rₖ) * Lorentz(Rₐ) * Lorentz(Rₙ), Λ; atol=ε) + end + + # Pure rotation → Rₖ ≈ Λ, Rₐ ≈ Rₙ ≈ 𝟙 + for Λ ∈ rand_rotations(T) + Rₖ, Rₐ, Rₙ = KAN(Λ) + @test same_Λ(Lorentz(Rₖ), Λ; atol=ε) + @test same_Λ(Lorentz(Rₐ), 𝟙; atol=ε) + @test same_Λ(Lorentz(Rₙ), 𝟙; atol=ε) + end + + # Pure z-boost → Rₐ ≈ Λ, Rₖ ≈ Rₙ ≈ 𝟙 + for φ ∈ T[0.3, 1.0, 2.5] + Λ = z_boost(T, φ) + Rₖ, Rₐ, Rₙ = KAN(Λ) + @test same_Λ(Lorentz(Rₖ), 𝟙; atol=ε) + @test same_Λ(Lorentz(Rₐ), Λ; atol=ε) + @test same_Λ(Lorentz(Rₙ), 𝟙; atol=ε) + end + + # Pure null rotation → Rₙ ≈ Λ, Rₖ ≈ Rₐ ≈ 𝟙 + for (ξˣ, ξʸ) ∈ [(T(0.4), T(-0.7)), (T(1.5), T(0.0)), (T(0.0), T(2.0))] + Λ = null_rotation(T, ξˣ, ξʸ) + εₙ = 64eps(T) * maximum(abs, Quaternionic.components(Λ)) + Rₖ, Rₐ, Rₙ = KAN(Λ) + @test same_Λ(Lorentz(Rₖ), 𝟙; atol=εₙ) + @test same_Λ(Lorentz(Rₐ), 𝟙; atol=εₙ) + @test same_Λ(Lorentz(Rₙ), Λ; atol=εₙ) + end + end +end + +@testitem "KAN: factor recovery (uniqueness)" tags=[:validation, :fast] setup=[LorentzDecompData] begin + import Quaternionic: KAN, components + using .LorentzDecompData: rand_rotations, null_rotation, z_boost, same_Λ + using DoubleFloats: Double64 + # Build Λ from known canonical factors K₀·A₀·N₀ and verify KAN recovers each. + # The Iwasawa decomposition is unique, so recovery is an oracle test. + for T ∈ (Float32, Float64, Double64) + for K₀ ∈ rand_rotations(T) + for φ ∈ T[0.5, 2.0] + A₀ = z_boost(T, φ) + for (ξˣ, ξʸ) ∈ [(T(0.4), T(-0.7)), (T(1.3), T(0.9))] + N₀ = null_rotation(T, ξˣ, ξʸ) + Λ = K₀ * A₀ * N₀ + ε = 1024eps(T) * maximum(abs, components(Λ)) + Rₖ, Rₐ, Rₙ = KAN(Λ) + @test same_Λ(Lorentz(Rₖ), K₀; atol=ε) + @test same_Λ(Lorentz(Rₐ), A₀; atol=ε) + @test same_Λ(Lorentz(Rₙ), N₀; atol=ε) + end + end + end + end +end + +@testitem "KAN: idempotency (decomposing a factor)" tags=[:validation, :fast] setup=[LorentzDecompData] begin + import Quaternionic: KAN, components + using .LorentzDecompData: rand_rotations, null_rotation, z_boost, same_Λ + using DoubleFloats: Double64 + # Decomposing an element that already lives in a single subgroup returns it in + # the matching slot and the identity in the other two. + for T ∈ (Float32, Float64, Double64) + 𝟙 = one(Lorentz{T}) + ε = 64eps(T) + + for K₀ ∈ rand_rotations(T) + Rₖ, Rₐ, Rₙ = KAN(K₀) + @test same_Λ(Lorentz(Rₖ), K₀; atol=ε) + @test same_Λ(Lorentz(Rₐ), 𝟙; atol=ε) + @test same_Λ(Lorentz(Rₙ), 𝟙; atol=ε) + end + + for φ ∈ T[0.5, 2.0] + A₀ = z_boost(T, φ) + Rₖ, Rₐ, Rₙ = KAN(A₀) + @test same_Λ(Lorentz(Rₖ), 𝟙; atol=ε) + @test same_Λ(Lorentz(Rₐ), A₀; atol=ε) + @test same_Λ(Lorentz(Rₙ), 𝟙; atol=ε) + end + + for (ξˣ, ξʸ) ∈ [(T(0.4), T(-0.7)), (T(1.3), T(0.9))] + N₀ = null_rotation(T, ξˣ, ξʸ) + εₙ = 64eps(T) * maximum(abs, components(N₀)) + Rₖ, Rₐ, Rₙ = KAN(N₀) + @test same_Λ(Lorentz(Rₖ), 𝟙; atol=εₙ) + @test same_Λ(Lorentz(Rₐ), 𝟙; atol=εₙ) + @test same_Λ(Lorentz(Rₙ), N₀; atol=εₙ) + end + end +end + +@testitem "KAN: multi-precision convergence" tags=[:validation, :slow] begin + import Quaternionic: KAN, components + using DoubleFloats: Double64 + # Algebraically direct algorithm: reconstruction error stays at the machine-epsilon + # floor across precisions (no Float64-specific constants or thresholds). + function recon_error(T) + c = sin(T(13)/20) / √T(3) + K₀ = Lorentz(rotor(cos(T(13)/20), c, c, c)) + A₀ = Boost(T(8)/10, QuatVec{T}(0, 0, 1)) + a, b = T(4)/10/√T(2), T(-7)/10/√T(2) + N₀ = Lorentz{T}(1, (im*a + b)/2, (im*b - a)/2, 0) + Λ = K₀ * A₀ * N₀ + Rₖ, Rₐ, Rₙ = KAN(Λ) + recon = Lorentz(Rₖ) * Lorentz(Rₐ) * Lorentz(Rₙ) + maximum(abs, components(recon) - components(Λ)) + end + @test recon_error(Float32) ≤ 10eps(Float32) + @test recon_error(Float64) ≤ 10eps(Float64) + @test recon_error(Double64) ≤ 10eps(Double64) + @test recon_error(BigFloat) ≤ 10eps(BigFloat) +end