diff --git a/algebra/rotor.py b/algebra/rotor.py index 208e0932..3b39075e 100644 --- a/algebra/rotor.py +++ b/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 (_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) + # _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 = _0``, ``H = _2``, ``W = _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 == _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. diff --git a/tests/test_rotor_power_general.py b/tests/test_rotor_power_general.py new file mode 100644 index 00000000..8454e3a3 --- /dev/null +++ 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