Merge pull request 'feat(algebra): exact fractional powers of non-simple rotors (invariant split)' (#23) from feat/rotor-power-nonsimple-bivector into main

algebra/rotor.py

@@ -8,13 +8,18 @@ it describes a transformation being applied, not a property of the vocabulary.
 
 import numpy as np
 
-from .cl41 import N_COMPONENTS, geometric_product, reverse
+from .cl41 import N_COMPONENTS, geometric_product, grade_project, reverse, scalar_part
 from .versor import unitize_versor, versor_condition
 
 _TRANSITION_CONDITION_TOL = 1e-4
 _NEAR_ZERO_TOL = 1e-12
 _SAME_POINT_TOL = 1e-6
 _STRICT_RESIDUE_TOL = 1e-2
+# A rotor is SIMPLE iff its grade-4 part vanishes (<R>_4 == 0 <=> R = R1 with a
+# single invariant plane). Above this, the rotor needs the invariant split.
+_SIMPLE_GRADE4_TOL = 1e-10
+# |discriminant| below this => the two invariant eigenvalues coincide (isoclinic).
+_DEGEN_TOL = 1e-9
 
 
 def _identity(dtype: np.dtype) -> np.ndarray:
@@ -75,39 +80,55 @@ def make_rotor_from_angle(angle: float, bivector_idx: int = 6) -> np.ndarray:
 def rotor_power(R: np.ndarray, alpha: float) -> np.ndarray:
     """Return R^alpha — the rotor on the manifold path from identity to R by alpha.
 
-    For a simple unit rotor decomposed as ``R = a + B`` (scalar + bivector):
+    EXACT for ANY closed unit rotor in Cl(4,1), simple or not. A general rotor
+    factors (invariant / bivector decomposition) into two commuting SIMPLE
+    rotors ``R = R1 R2`` with distinct invariant planes; then, because they
+    commute, ``R^α = R1^α R2^α`` and each factor uses the simple closed form
+    below. The isoclinic case (coincident invariant planes) has its own closed
+    form. There is no iteration, no approximation, and no external library —
+    the split is built from the Cl(4,1) geometric product alone.
+
+    Simple factor ``R_i = a + B`` (scalar + simple bivector):
 
     - rotation plane (``B² < 0``):  ``R^α = cos(α·θ/2) + (sin(α·θ/2)/|B|) · B``
       where ``θ/2 = atan2(|B|, a)``.
     - boost plane (``B² > 0``):     ``R^α = cosh(α·η/2) + (sinh(α·η/2)/|B|) · B``
       where ``η/2 = atanh(|B|/a)``.
 
-    This is the proper slerp on the rotor manifold: it stays on the manifold
-    by construction, so ``versor_condition(rotor_power(R, α)) < 1e-6`` for any
-    α whenever ``R`` is itself a closed unit rotor.
-
-    Falls back to the identity rotor when ``R`` is not a closed scalar+bivector
-    rotor (e.g. carries higher-grade components or a non-simple bivector) so
-    callers never receive a manifold-violating output.
+    The result stays on the rotor manifold by construction, so
+    ``versor_condition(rotor_power(R, α)) < 1e-6`` for any α whenever ``R`` is a
+    closed unit rotor. (Historically this returned the *identity* for non-simple
+    rotors — an approximation where exactness was available, which silently
+    collapsed geodesic interpolation to a no-op. That corner is now closed.)
     """
     R_arr = np.asarray(R, dtype=np.float64)
     if R_arr.shape != (N_COMPONENTS,):
         raise ValueError(
             f"rotor_power expects a {N_COMPONENTS}-component rotor; got {R_arr.shape}."
         )
-
     dtype = _result_dtype(R_arr)
+    # <R>_4 == 0  <=>  R is a single simple rotor. Otherwise take the split path.
+    if float(np.linalg.norm(grade_project(R_arr, 4))) >= _SIMPLE_GRADE4_TOL:
+        return _general_rotor_power(R_arr, alpha, dtype)
+    return _simple_rotor_power(R_arr, alpha, dtype)
+
+
+def _simple_rotor_power(R_arr: np.ndarray, alpha: float, dtype: np.dtype) -> np.ndarray:
+    """R^alpha for a SIMPLE rotor (scalar + one simple bivector). Exact closed form.
+
+    Behaviour is unchanged from the original ``rotor_power`` on simple inputs.
+    """
     a = float(R_arr[0])
     B = R_arr.copy()
     B[0] = 0.0
 
-    # Quick guard: bivector must be a simple bivector (B² is grade-0 only).
+    # A simple rotor's bivector squares to a scalar (B² is grade-0 only).
     B_sq_full = geometric_product(B, B).astype(np.float64)
     bsq_scalar = float(B_sq_full[0])
     B_sq_higher = B_sq_full.copy()
     B_sq_higher[0] = 0.0
     if float(np.linalg.norm(B_sq_higher)) > 1e-6:
-        # Non-simple bivector — return identity to avoid drift.
+        # Not a simple bivector — should not reach here via the public dispatch.
         return _identity(dtype)
 
     # Near-identity: nothing to scale.
@@ -144,6 +165,92 @@ def rotor_power(R: np.ndarray, alpha: float) -> np.ndarray:
     return result.astype(dtype, copy=False)
 
 
+def _isoclinic_power_coeffs(x: float, alpha: float) -> tuple[float, float, float]:
+    """Power coefficients ``(A, f, c)`` for one of two identical (isoclinic) simple
+    factors with ``c² = x``: ``R_i^α = A + f · G_i``. Handles rotation, boost, and
+    the null limit uniformly.
+    """
+    gsq = x - 1.0
+    c = float(np.sqrt(max(x, 0.0)))
+    if gsq < -1e-15:                        # rotation: c = cos(theta)
+        theta = float(np.arccos(min(1.0, max(-1.0, c))))
+        slin = float(np.sin(theta))
+        A = float(np.cos(alpha * theta))
+        f = float(np.sin(alpha * theta) / slin) if slin > 1e-300 else float(alpha)
+    elif gsq > 1e-15:                       # boost: c = cosh(eta)
+        eta = float(np.arccosh(max(1.0, c)))
+        slin = float(np.sinh(eta))
+        A = float(np.cosh(alpha * eta))
+        f = float(np.sinh(alpha * eta) / slin) if slin > 1e-300 else float(alpha)
+    else:                                   # null / parabolic limit
+        A, f = 1.0, float(alpha)
+    return A, f, c
+
+
+def _split_commuting_simple(
+    P: float, H: np.ndarray, W: np.ndarray, h0: float, disc: float
+) -> tuple[np.ndarray, np.ndarray]:
+    """Invariant decomposition of a non-simple rotor into two commuting SIMPLE
+    unit rotors ``R = R1 R2`` (distinct-eigenvalue branch).
+
+    With ``P = <R>_0``, ``H = <R>_2``, ``W = <R>_4``: the squared scalars of the
+    two simple factors are ``x_i = c_i²`` — the roots of ``t² − (2P²−h0) t + P²``
+    — and each simple bivector ``G_i`` is recovered by the linear system in
+    ``{H, HW}``. Returns ``(R1, R2)`` as 32-component rotors.
+    """
+    b = 2.0 * P * P - h0
+    sq = float(np.sqrt(disc))
+    x1 = 0.5 * (b + sq)
+    x2 = 0.5 * (b - sq)
+    c1 = float(np.sqrt(max(x1, 0.0)))
+    c2 = float(np.sqrt(max(x2, 0.0)))
+    if P < 0.0:
+        c2 = -c2                            # fix product sign so c1·c2 == <R>_0
+    g1sq = x1 - 1.0
+    g2sq = x2 - 1.0
+    HW = grade_project(geometric_product(H, W), 2).astype(np.float64)
+    det = c2 * c2 * g1sq - c1 * c1 * g2sq
+    if abs(det) < _NEAR_ZERO_TOL:
+        raise ValueError(
+            "rotor_power: singular invariant split (unexpected for distinct eigenvalues)"
+        )
+    G1 = (c2 * g1sq * H - c1 * HW) / det
+    G2 = (c2 * HW - c1 * g2sq * H) / det
+    R1 = G1.copy()
+    R1[0] = c1
+    R2 = G2.copy()
+    R2[0] = c2
+    return R1, R2
+
+
+def _general_rotor_power(R_arr: np.ndarray, alpha: float, dtype: np.dtype) -> np.ndarray:
+    """R^alpha for a NON-simple rotor via the invariant (bivector) decomposition."""
+    P = float(R_arr[0])
+    H = grade_project(R_arr, 2).astype(np.float64)
+    W = grade_project(R_arr, 4).astype(np.float64)
+    h0 = float(scalar_part(geometric_product(H, H)))
+    b = 2.0 * P * P - h0
+    disc = b * b - 4.0 * P * P
+    if disc <= _DEGEN_TOL:
+        # Isoclinic: coincident invariant planes (x1 == x2 == b/2). The result
+        # depends only on the symmetric functions H and W, so no per-plane split
+        # is needed: R^α = A² + (A·f/c)·H + f²·W.
+        A, f, c = _isoclinic_power_coeffs(0.5 * b, alpha)
+        if c < _NEAR_ZERO_TOL:
+            raise ValueError(
+                "rotor_power: isoclinic rotor at theta~pi/2 has no principal power"
+            )
+        out = (A * f / c) * H + (f * f) * W
+        out[0] += A * A
+        return out.astype(dtype, copy=False)
+    R1, R2 = _split_commuting_simple(P, H, W, h0, disc)
+    out = geometric_product(
+        _simple_rotor_power(R1, alpha, np.dtype(np.float64)),
+        _simple_rotor_power(R2, alpha, np.dtype(np.float64)),
+    )
+    return out.astype(dtype, copy=False)
+
+
 def word_transition_rotor(A: np.ndarray, B: np.ndarray) -> np.ndarray:
     """
     Compute the closed transition operator from source versor A to target B.

tests/test_rotor_power_general.py

@@ -0,0 +1,135 @@
+"""Exact fractional powers of GENERAL (non-simple) rotors in Cl(4,1).
+
+`rotor_power` previously returned the identity for any non-simple rotor — an
+approximation where exactness was available (Pillar II), which silently
+collapsed geodesic interpolation (slerp, supervised blend) to a no-op. It now
+uses the invariant (bivector) decomposition: a general rotor factors into two
+commuting simple rotors, R = R1 R2, so R^a = R1^a R2^a exactly, with a closed
+form for the isoclinic case. First-principles, no library (Pillar III).
+
+These tests pin the group-theoretic identities to machine precision on rotors
+that exercise every plane type (Euclidean rotations, e5 boosts, mixed, and the
+isoclinic degenerate case) — the regimes the pre-existing `test_rotor_power.py`
+(simple rotors only) never covered.
+"""
+
+from __future__ import annotations
+
+import numpy as np
+import pytest
+
+from algebra.cl41 import N_COMPONENTS, geometric_product, grade_project
+from algebra.rotor import make_rotor_from_angle, rotor_power
+from algebra.versor import versor_condition
+
+_ROT = 1e-8          # power-identity tolerance
+_CLOSE = 1e-6        # the versor_condition invariant
+
+
+def _identity() -> np.ndarray:
+    v = np.zeros(N_COMPONENTS, dtype=np.float64)
+    v[0] = 1.0
+    return v
+
+
+def _rotor(seed: int, nplanes: int, lo: float = -1.2, hi: float = 1.2) -> np.ndarray:
+    """A reproducible composed rotor over `nplanes` distinct planes (grades e1..e5)."""
+    rng = np.random.default_rng(seed)
+    v = _identity()
+    planes = rng.choice(range(6, 16), size=nplanes, replace=False)
+    for idx in planes:
+        v = geometric_product(v, make_rotor_from_angle(float(rng.uniform(lo, hi)), bivector_idx=int(idx)))
+    return np.asarray(v, dtype=np.float64)
+
+
+def _is_non_simple(R: np.ndarray) -> bool:
+    return float(np.linalg.norm(grade_project(R, 4))) > 1e-9
+
+
+# A spread of seeds; most compose into non-simple rotors.
+_SEEDS = list(range(40))
+
+
+@pytest.mark.parametrize("seed", _SEEDS)
+def test_power_one_reconstructs_rotor(seed):
+    R = _rotor(seed, nplanes=3)
+    assert np.linalg.norm(rotor_power(R, 1.0) - R) < _ROT
+
+
+@pytest.mark.parametrize("seed", _SEEDS)
+def test_half_power_squares_to_rotor(seed):
+    R = _rotor(seed, nplanes=3)
+    half = rotor_power(R, 0.5)
+    assert np.linalg.norm(geometric_product(half, half) - R) < _ROT
+
+
+@pytest.mark.parametrize("seed", _SEEDS)
+def test_group_law_additive_exponents(seed):
+    R = _rotor(seed, nplanes=4)
+    lhs = geometric_product(rotor_power(R, 0.3), rotor_power(R, 0.45))
+    rhs = rotor_power(R, 0.75)
+    assert np.linalg.norm(lhs - rhs) < _ROT
+
+
+@pytest.mark.parametrize("seed", _SEEDS)
+def test_zero_power_is_identity(seed):
+    R = _rotor(seed, nplanes=3)
+    assert np.linalg.norm(rotor_power(R, 0.0) - _identity()) < 1e-12
+
+
+@pytest.mark.parametrize("seed", _SEEDS)
+@pytest.mark.parametrize("alpha", [0.0, 0.25, 0.5, 0.75, 1.0, 1.5])
+def test_closure_preserved(seed, alpha):
+    R = _rotor(seed, nplanes=3)
+    assert versor_condition(rotor_power(R, alpha)) < _CLOSE
+
+
+def test_covers_non_simple_rotors():
+    """Guard: the seed pool actually exercises the non-simple path (else the suite
+    would be vacuous, à la the fidelity finding it fixes)."""
+    n = sum(_is_non_simple(_rotor(s, 3)) for s in _SEEDS)
+    assert n >= len(_SEEDS) // 2
+
+
+@pytest.mark.parametrize("seed", _SEEDS[:20])
+def test_non_simple_power_is_not_identity(seed):
+    """The exact regression for the bug: a non-simple rotor's interior power must
+    MOVE off the identity (the old code returned identity → no-op geodesic)."""
+    R = _rotor(seed, nplanes=3)
+    if not _is_non_simple(R):
+        pytest.skip("simple rotor")
+    mid = rotor_power(R, 0.5)
+    assert np.linalg.norm(mid - _identity()) > 1e-3
+
+
+def test_isoclinic_degenerate_is_exact():
+    """Coincident invariant planes (same-angle disjoint Euclidean rotations) take the
+    closed-form isoclinic branch and reconstruct exactly."""
+    for (p, q, th) in [(6, 13, 0.7), (6, 13, 1.2), (7, 15, 0.4)]:
+        R = geometric_product(
+            make_rotor_from_angle(th, bivector_idx=p),
+            make_rotor_from_angle(th, bivector_idx=q),
+        )
+        assert _is_non_simple(R)
+        h = rotor_power(R, 0.5)
+        assert np.linalg.norm(geometric_product(h, h) - R) < 1e-10
+        assert np.linalg.norm(rotor_power(R, 1.0) - R) < 1e-10
+
+
+@pytest.mark.parametrize("seed", _SEEDS[:15])
+def test_determinism_replay(seed):
+    """Same input → byte-identical output (replay determinism, Pillar II)."""
+    R = _rotor(seed, nplanes=3)
+    a = rotor_power(R, 0.37)
+    b = rotor_power(R, 0.37)
+    assert np.array_equal(a, b)
+
+
+def test_backward_compat_simple_rotor_matches_analytic():
+    """On a SIMPLE rotor, R^a is the analytic half-angle interpolation, unchanged."""
+    theta = 0.9
+    R = make_rotor_from_angle(theta, bivector_idx=6)
+    for alpha in (0.0, 0.25, 0.5, 0.75, 1.0):
+        got = rotor_power(R, alpha)
+        expect = make_rotor_from_angle(alpha * theta, bivector_idx=6)
+        assert np.linalg.norm(got - expect) < 1e-12

tests/test_third_door_blueprint_fidelity.py

@@ -15,9 +15,9 @@ citing the blueprint section + audit finding. When an operator is implemented to
 spec, its xfail flips to xpass (strict) and CI forces the marker's removal.
 
 Empirical findings (2026-07-11 audit, reproduced deterministically below):
-  #1  supervised_blend is a no-op for interior alpha on composed versors
-      (rotor_power returns identity for the non-simple transition rotor);
-      versor_condition stays ~1e-16 and masks it.
+  #1  [RESOLVED by #23] supervised_blend was a no-op for interior alpha on
+      composed versors (rotor_power returned identity for the non-simple
+      transition rotor). Exact fractional powers now implemented.
   #2  cartan_iwasawa_factorize raises "factor R not closed" on composed
       conformal versors (~45% of the time), instead of the "mathematically
       exact, guaranteed" decomposition the Super-Blueprint §2.2 specifies.
@@ -58,34 +58,6 @@ def _composed_versor(planes: tuple[int, ...], seed: float) -> np.ndarray:
 
 
 # --- Finding #1: supervised_blend geodesic ----------------------------------
-
-def test_supervised_blend_currently_degenerates_on_composed_versors():
-    """CHARACTERIZATION of finding #1 — locks the current (defective) behaviour.
-
-    On a composed pair, the "geodesic" collapses onto the source for every
-    interior alpha while closure stays green. This test PASSES today. When
-    supervised_blend is fixed to spec it SHOULD start failing — at which point
-    delete this test and drop the xfail on the spec test below.
-    """
-    a = _composed_versor((6, 7, 8, 10, 11), seed=1.0)
-    b = _composed_versor((6, 7, 8, 10, 11), seed=2.0)
-    m = GoldTetherMonitor()
-    for alpha in (0.25, 0.5, 0.75):
-        mid = m.supervised_blend(a, b, alpha)
-        assert np.allclose(mid, a, atol=1e-12), "no-op degeneration expected (finding #1)"
-
-
-@pytest.mark.xfail(
-    reason=(
-        "Finding #1 / R&D-Revised §2.3: supervised_blend applies rotor_power to the "
-        "whole (non-simple) transition rotor B*rev(A); rotor_power returns identity "
-        "for non-simple bivectors, so the blend is a no-op for interior alpha and only "
-        "reaches the target via the explicit alpha>=1 shortcut. versor_condition passes "
-        "and cannot detect it. Spec requires a genuine BCH-free geodesic (Cartan-Iwasawa "
-        "factor-wise slerp) that interpolates strictly between the endpoints."
-    ),
-    strict=True,
-)
 def test_supervised_blend_should_interpolate_composed_versors():
     a = _composed_versor((6, 7, 8, 10, 11), seed=1.0)
     b = _composed_versor((6, 7, 8, 10, 11), seed=2.0)