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 bc968e9..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(Λ)
@@ -395,3 +401,45 @@ 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
+# ---------------------------------------------------------------------------
+
+@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ₙ
+    ξˣ╱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/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,
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