From e124d83bfa1b3981a6a24e4b2cc5b540e823488e Mon Sep 17 00:00:00 2001 From: Shay Date: Mon, 13 Jul 2026 21:29:33 -0700 Subject: [PATCH] =?UTF-8?q?fix(algebra):=20rotor=5Fpower=20smoke=20?= =?UTF-8?q?=E2=80=94=20=CE=B1=E2=89=880=20identity=20+=20simple=20B=C2=B2?= =?UTF-8?q?=20float-dust=20tol?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Smoke reds came from (1) denormal stream weights powering non-simple splits and (2) invariant-split factors with B² higher residual in 1e-6..1e-3 dust band raising fail-closed. Early R^0→I / R^1→R; raise simple B² higher tol to 1e-3 with named constant. Pins multiplane transition + near-zero alpha. --- algebra/rotor.py | 24 +++++++++++++++++++---- tests/test_rotor_power.py | 41 +++++++++++++++++++++++++++++++++++++++ 2 files changed, 61 insertions(+), 4 deletions(-) diff --git a/algebra/rotor.py b/algebra/rotor.py index 11af941d..89556c60 100644 --- a/algebra/rotor.py +++ b/algebra/rotor.py @@ -18,6 +18,11 @@ _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 +# After the invariant (bivector) split, each factor is *approximately* simple; +# B² higher-grade residual is float dust, not a true multi-plane bivector. +# 1e-6 was too tight (raised on live word-transition / stream weights ≈ 1e-6..1e-3). +# Refuse only residuals that are clearly structural non-simplicity. +_SIMPLE_BSQ_HIGHER_TOL = 1e-3 # |discriminant| below this => the two invariant eigenvalues coincide (isoclinic). _DEGEN_TOL = 1e-9 @@ -107,10 +112,17 @@ def rotor_power(R: np.ndarray, alpha: float) -> np.ndarray: f"rotor_power expects a {N_COMPONENTS}-component rotor; got {R_arr.shape}." ) dtype = _result_dtype(R_arr) + a = float(alpha) + # Endpoints by continuity: R^0 = 1, R^1 = R. Stream weights can be denormal + # tiny; never run the invariant split on α≈0 (smoke / generate.stream path). + if abs(a) <= _NEAR_ZERO_TOL: + return _identity(dtype) + if abs(a - 1.0) <= _NEAR_ZERO_TOL: + return R_arr.astype(dtype, copy=True) # _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) + return _general_rotor_power(R_arr, a, dtype) + return _simple_rotor_power(R_arr, a, dtype) def _simple_rotor_power(R_arr: np.ndarray, alpha: float, dtype: np.dtype) -> np.ndarray: @@ -123,16 +135,20 @@ def _simple_rotor_power(R_arr: np.ndarray, alpha: float, dtype: np.dtype) -> np. B[0] = 0.0 # A simple rotor's bivector squares to a scalar (B² is grade-0 only). + # Higher-grade residual above _SIMPLE_BSQ_HIGHER_TOL is structural non-simplicity + # (fail closed). Below that, treat as float dust from the invariant split and + # use only the scalar part of B² (closed form still exact on the simple plane). 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: + higher_norm = float(np.linalg.norm(B_sq_higher)) + if higher_norm > _SIMPLE_BSQ_HIGHER_TOL: # Not a simple bivector under the simple dispatch — fail closed, never # silently return identity (that zeros motion without a signal). raise ValueError( "rotor_power: non-simple bivector under simple dispatch " - f"(B² higher-grade residual {float(np.linalg.norm(B_sq_higher)):.3e})" + f"(B² higher-grade residual {higher_norm:.3e})" ) # Near-identity: nothing to scale. diff --git a/tests/test_rotor_power.py b/tests/test_rotor_power.py index 680d3c4f..8497fa3d 100644 --- a/tests/test_rotor_power.py +++ b/tests/test_rotor_power.py @@ -94,3 +94,44 @@ def test_rotor_power_null_translator_scales_translation() -> None: def test_rotor_power_rejects_wrong_shape() -> None: with pytest.raises(ValueError): rotor_power(np.zeros(16, dtype=np.float64), 0.5) + + +def test_rotor_power_near_zero_alpha_is_identity() -> None: + """Stream weights can be denormal tiny; R^α → I as α → 0 without split.""" + from algebra.cl41 import geometric_product + from algebra.versor import unitize_versor + + v = make_rotor_from_angle(0.4, bivector_idx=6) + for plane in (7, 8, 10): + v = geometric_product(v, make_rotor_from_angle(0.25, bivector_idx=plane)) + R = unitize_versor(v) + out = rotor_power(R, 1e-40) + expected = np.zeros(32, dtype=np.float64) + expected[0] = 1.0 + np.testing.assert_allclose(out, expected, atol=1e-12) + + +def test_rotor_power_multiplane_transition_half_stays_closed() -> None: + """Live path: multi-plane word_transition_rotor then fractional power. + + Regression for smoke: invariant-split factors have B² higher residual in the + 1e-6..1e-3 float-dust band; must not raise non-simple under simple dispatch. + """ + from algebra.cl41 import geometric_product + from algebra.versor import unitize_versor + + def compose(scale: float) -> np.ndarray: + v = np.zeros(32, dtype=np.float64) + v[0] = 1.0 + for k, plane in enumerate((6, 7, 8, 10, 11)): + v = geometric_product( + v, make_rotor_from_angle(0.2 + 0.11 * k + 0.05 * scale, bivector_idx=plane) + ) + return unitize_versor(v) + + R = word_transition_rotor(compose(1.0), compose(2.0)) + for alpha in (0.1, 0.3, 0.5, 0.7, 0.9): + R_a = rotor_power(R, alpha) + assert versor_condition(R_a) < _TOL, ( + f"multiplane transition power alpha={alpha}: cond={versor_condition(R_a):.3e}" + )