# This Source Code Form is subject to the terms of the Mozilla Public
# License, v. 2.0. If a copy of the MPL was not distributed with this
# file, You can obtain one at https://mozilla.org/MPL/2.0/.
#
# Copyright (c) 2026 Jan Hettenkofer
"""Small dense linear algebra written by hand, validated against numpy.
Scope at milestone M1 is deliberately narrow: an immutable 3-vector, an immutable
3x3 matrix, and a symmetric 3x3 eigendecomposition. The dense solve, the least-squares
fit, the matrix-free conjugate gradient and the geometric fits catalogued in PRD 11.13
are added by the milestones that need them; the Tier E sparse solver in particular is
**not** part of M1.
Why this is hand-written rather than a numpy call: the surface numpy would genuinely
replace is about ninety lines (PRD 11.13). numpy is feasible inside Fusion and runs at
full native speed, but costs 56 MB per platform pinned to the cp314 ABI to buy headroom
this workload does not need (ADR-013). The other half of that decision is not optional:
**every routine here is differentially tested against numpy over randomized input**
(PRD 13.3 T14, ``tests/unit/test_linalg.py``), and the numeric core runs with
``RuntimeWarning`` promoted to an error.
Units:
This module is unit-agnostic. Callers work in millimeters and radians (PRD 7.6);
nothing here inspects or converts units.
Hazards:
:func:`eigen_sym3` carries the one real bug the numpy oracle has already caught —
an overflow in the Jacobi rotation on a near-circular cross-section (PRD 11.13).
See its docstring before changing it.
"""
from __future__ import annotations
import math
from dataclasses import dataclass
from typing import Final
__all__ = ["Mat3", "Vec3", "eigen_sym3"]
#: Machine epsilon for IEEE-754 binary64, used to scale convergence thresholds.
_EPS: Final = 2.220446049250313e-16
#: Maximum Jacobi sweeps before :func:`eigen_sym3` gives up. A symmetric 3x3 needs
#: three or four sweeps in practice; 50 is a deliberately loose safety margin, an order
#: of magnitude above the observed worst case, and is never reached for a finite input.
_MAX_SWEEPS: Final = 50
[docs]
@dataclass(frozen=True, slots=True)
class Vec3:
"""An immutable 3-vector.
The core exchanges plain value objects across the adapter boundary rather than
``adsk`` handles (PRD 12.3, ADR-002), which is why this type exists even though
numpy would otherwise supply it (PRD 11.13).
Attributes:
x: First component. Unit-agnostic; millimeters in geometric use (PRD 7.6).
y: Second component.
z: Third component.
Example:
>>> Vec3(1.0, 2.0, 3.0) + Vec3(0.5, 0.5, 0.5)
Vec3(x=1.5, y=2.5, z=3.5)
"""
x: float
y: float
z: float
# -- construction ------------------------------------------------------------
[docs]
@classmethod
def zero(cls) -> Vec3:
"""Return the zero vector.
Returns:
``Vec3(0.0, 0.0, 0.0)``.
Example:
>>> Vec3.zero()
Vec3(x=0.0, y=0.0, z=0.0)
"""
return cls(0.0, 0.0, 0.0)
[docs]
@classmethod
def from_tuple(cls, values: tuple[float, float, float]) -> Vec3:
"""Build a vector from a 3-tuple.
Args:
values: ``(x, y, z)``.
Returns:
The equivalent :class:`Vec3`.
Example:
>>> Vec3.from_tuple((1.0, 0.0, 0.0))
Vec3(x=1.0, y=0.0, z=0.0)
"""
return cls(values[0], values[1], values[2])
[docs]
def as_tuple(self) -> tuple[float, float, float]:
"""Return the components as a plain tuple.
Returns:
``(x, y, z)``.
Example:
>>> Vec3(1.0, 2.0, 3.0).as_tuple()
(1.0, 2.0, 3.0)
"""
return (self.x, self.y, self.z)
# -- arithmetic --------------------------------------------------------------
def __add__(self, other: Vec3) -> Vec3:
return Vec3(self.x + other.x, self.y + other.y, self.z + other.z)
def __sub__(self, other: Vec3) -> Vec3:
return Vec3(self.x - other.x, self.y - other.y, self.z - other.z)
def __neg__(self) -> Vec3:
return Vec3(-self.x, -self.y, -self.z)
def __mul__(self, scalar: float) -> Vec3:
return Vec3(self.x * scalar, self.y * scalar, self.z * scalar)
def __rmul__(self, scalar: float) -> Vec3:
return self.__mul__(scalar)
def __truediv__(self, scalar: float) -> Vec3:
return Vec3(self.x / scalar, self.y / scalar, self.z / scalar)
[docs]
def scaled(self, scalar: float) -> Vec3:
"""Return this vector multiplied by ``scalar``.
Named alternative to ``v * scalar`` for call sites where the operator would be
ambiguous next to a dot or cross product.
Args:
scalar: The multiplier.
Returns:
``scalar * self``.
Example:
>>> Vec3(1.0, 2.0, 3.0).scaled(2.0)
Vec3(x=2.0, y=4.0, z=6.0)
"""
return self.__mul__(scalar)
# -- products and norms ------------------------------------------------------
[docs]
def dot(self, other: Vec3) -> float:
r"""Return the Euclidean inner product.
:math:`a \cdot b = a_x b_x + a_y b_y + a_z b_z`.
In plain language: how much of ``other`` points along ``self``, times the
lengths of both.
Args:
other: The second operand.
Returns:
The scalar product, in the square of the caller's length unit.
Example:
>>> Vec3(1.0, 2.0, 3.0).dot(Vec3(4.0, -5.0, 6.0))
12.0
"""
return self.x * other.x + self.y * other.y + self.z * other.z
[docs]
def cross(self, other: Vec3) -> Vec3:
r"""Return the right-handed cross product.
:math:`a \times b = (a_y b_z - a_z b_y,\; a_z b_x - a_x b_z,\;
a_x b_y - a_y b_x)`.
In plain language: a vector perpendicular to both inputs, whose length is the
area of the parallelogram they span, pointing the way a right-handed screw
turning ``self`` into ``other`` would advance.
Hazard:
Handedness (PRD 7.1). This is the right-handed convention and the template
coordinate convention of PRD 4 depends on it. A mirrored template looks
perfect and fits nothing.
Args:
other: The second operand.
Returns:
``self x other``.
Example:
>>> Vec3(1.0, 0.0, 0.0).cross(Vec3(0.0, 1.0, 0.0))
Vec3(x=0.0, y=0.0, z=1.0)
"""
return Vec3(
self.y * other.z - self.z * other.y,
self.z * other.x - self.x * other.z,
self.x * other.y - self.y * other.x,
)
[docs]
def norm_squared(self) -> float:
"""Return the squared Euclidean length.
Prefer this to ``norm() ** 2`` when comparing magnitudes: it avoids a square
root and the rounding it introduces.
Returns:
``self . self``.
Example:
>>> Vec3(3.0, 4.0, 0.0).norm_squared()
25.0
"""
return self.x * self.x + self.y * self.y + self.z * self.z
[docs]
def norm(self) -> float:
"""Return the Euclidean length.
Uses :func:`math.hypot`, which is correctly rounded and does not overflow for
inputs whose squares would (a plain ``sqrt(x*x + y*y + z*z)`` overflows above
about 1.3e154).
Returns:
The length, in the caller's length unit.
Example:
>>> Vec3(3.0, 4.0, 0.0).norm()
5.0
"""
return math.hypot(self.x, self.y, self.z)
[docs]
def normalized(self, *, tolerance: float = 1e-12) -> Vec3:
"""Return a unit vector in the same direction.
Args:
tolerance: Lengths at or below this are treated as zero and rejected.
Default 1e-12, chosen because the core works in millimeters
(PRD 7.6) and 1e-12 mm is a picometer — far below any length that
can arise from CAD geometry, so a vector this short is a degenerate
direction rather than a small one.
Returns:
``self / |self|``.
Raises:
ValueError: If ``|self| <= tolerance``; a direction cannot be recovered
from a zero-length vector and silently returning one would produce a
plausible but wrong frame (PRD 7.10).
Example:
>>> Vec3(0.0, 5.0, 0.0).normalized()
Vec3(x=0.0, y=1.0, z=0.0)
"""
length = self.norm()
if length <= tolerance:
msg = (
f"cannot normalize a vector of length {length!r} "
f"(tolerance {tolerance!r}); the direction is undefined"
)
raise ValueError(msg)
return Vec3(self.x / length, self.y / length, self.z / length)
[docs]
@dataclass(frozen=True, slots=True)
class Mat3:
"""An immutable 3x3 matrix, stored as three row vectors.
Row-major throughout: ``row0`` is the first row, so ``Mat3.apply`` computes
``M v`` with ``v`` treated as a column.
Attributes:
row0: First row.
row1: Second row.
row2: Third row.
Example:
>>> Mat3.identity().apply(Vec3(1.0, 2.0, 3.0))
Vec3(x=1.0, y=2.0, z=3.0)
"""
row0: Vec3
row1: Vec3
row2: Vec3
# -- construction ------------------------------------------------------------
[docs]
@classmethod
def identity(cls) -> Mat3:
"""Return the 3x3 identity matrix.
Returns:
``I``.
Example:
>>> Mat3.identity().determinant()
1.0
"""
return cls(
Vec3(1.0, 0.0, 0.0),
Vec3(0.0, 1.0, 0.0),
Vec3(0.0, 0.0, 1.0),
)
[docs]
@classmethod
def zero(cls) -> Mat3:
"""Return the 3x3 zero matrix.
Returns:
The all-zero matrix.
Example:
>>> Mat3.zero().determinant()
0.0
"""
return cls(Vec3.zero(), Vec3.zero(), Vec3.zero())
[docs]
@classmethod
def from_rows(cls, rows: tuple[tuple[float, float, float], ...]) -> Mat3:
"""Build a matrix from three row 3-tuples.
Args:
rows: Exactly three ``(a, b, c)`` tuples, in row order.
Returns:
The matrix whose ``i``-th row is ``rows[i]``.
Raises:
ValueError: If ``rows`` does not contain exactly three entries.
Example:
>>> diagonal = Mat3.from_rows(((1.0, 0.0, 0.0), (0.0, 2.0, 0.0), (0.0, 0.0, 3.0)))
>>> diagonal.determinant()
6.0
"""
if len(rows) != 3:
msg = f"a Mat3 needs exactly 3 rows, got {len(rows)}"
raise ValueError(msg)
return cls(
Vec3.from_tuple(rows[0]),
Vec3.from_tuple(rows[1]),
Vec3.from_tuple(rows[2]),
)
[docs]
@classmethod
def from_columns(cls, columns: tuple[Vec3, Vec3, Vec3]) -> Mat3:
"""Build a matrix whose columns are the given vectors.
This is the natural constructor for a frame or an eigenvector basis, where each
column is a direction.
Args:
columns: The three column vectors, in order.
Returns:
The matrix ``[c0 | c1 | c2]``.
Example:
>>> basis = (Vec3(0.0, 1.0, 0.0), Vec3(-1.0, 0.0, 0.0), Vec3(0.0, 0.0, 1.0))
>>> Mat3.from_columns(basis).apply(Vec3(1.0, 0.0, 0.0))
Vec3(x=0.0, y=1.0, z=0.0)
"""
c0, c1, c2 = columns
return cls(
Vec3(c0.x, c1.x, c2.x),
Vec3(c0.y, c1.y, c2.y),
Vec3(c0.z, c1.z, c2.z),
)
[docs]
def rows(self) -> tuple[Vec3, Vec3, Vec3]:
"""Return the three rows.
Returns:
``(row0, row1, row2)``.
Example:
>>> Mat3.identity().rows()[1]
Vec3(x=0.0, y=1.0, z=0.0)
"""
return (self.row0, self.row1, self.row2)
[docs]
def columns(self) -> tuple[Vec3, Vec3, Vec3]:
"""Return the three columns.
Returns:
``(col0, col1, col2)``.
Example:
>>> Mat3.identity().columns()[2]
Vec3(x=0.0, y=0.0, z=1.0)
"""
return (
Vec3(self.row0.x, self.row1.x, self.row2.x),
Vec3(self.row0.y, self.row1.y, self.row2.y),
Vec3(self.row0.z, self.row1.z, self.row2.z),
)
[docs]
def at(self, i: int, j: int) -> float:
"""Return the element in row ``i``, column ``j``.
Args:
i: Row index, 0-based.
j: Column index, 0-based.
Returns:
``M[i][j]``.
Raises:
IndexError: If either index is outside ``0..2``.
Example:
>>> Mat3.identity().at(2, 2)
1.0
"""
if not (0 <= i <= 2 and 0 <= j <= 2):
msg = f"Mat3 index out of range: ({i}, {j})"
raise IndexError(msg)
return self.rows()[i].as_tuple()[j]
# -- arithmetic --------------------------------------------------------------
def __add__(self, other: Mat3) -> Mat3:
return Mat3(
self.row0 + other.row0,
self.row1 + other.row1,
self.row2 + other.row2,
)
def __sub__(self, other: Mat3) -> Mat3:
return Mat3(
self.row0 - other.row0,
self.row1 - other.row1,
self.row2 - other.row2,
)
def __mul__(self, scalar: float) -> Mat3:
return Mat3(self.row0 * scalar, self.row1 * scalar, self.row2 * scalar)
def __rmul__(self, scalar: float) -> Mat3:
return self.__mul__(scalar)
[docs]
def apply(self, vector: Vec3) -> Vec3:
"""Return the matrix-vector product ``M v``.
Args:
vector: The column vector to transform.
Returns:
``M v``.
Example:
>>> Mat3.identity().apply(Vec3(1.0, 2.0, 3.0))
Vec3(x=1.0, y=2.0, z=3.0)
"""
return Vec3(
self.row0.dot(vector),
self.row1.dot(vector),
self.row2.dot(vector),
)
[docs]
def matmul(self, other: Mat3) -> Mat3:
"""Return the matrix product ``self @ other``.
Args:
other: The right-hand matrix.
Returns:
``self @ other``.
Example:
>>> Mat3.identity().matmul(Mat3.identity()) == Mat3.identity()
True
"""
cols = other.columns()
return Mat3(
Vec3(self.row0.dot(cols[0]), self.row0.dot(cols[1]), self.row0.dot(cols[2])),
Vec3(self.row1.dot(cols[0]), self.row1.dot(cols[1]), self.row1.dot(cols[2])),
Vec3(self.row2.dot(cols[0]), self.row2.dot(cols[1]), self.row2.dot(cols[2])),
)
def __matmul__(self, other: Mat3) -> Mat3:
return self.matmul(other)
[docs]
def transposed(self) -> Mat3:
"""Return the transpose.
Returns:
``M^T``.
Example:
>>> upper = Mat3.from_rows(((1.0, 2.0, 3.0), (0.0, 1.0, 0.0), (0.0, 0.0, 1.0)))
>>> upper.transposed().row0
Vec3(x=1.0, y=0.0, z=0.0)
"""
return Mat3.from_columns((self.row0, self.row1, self.row2))
[docs]
def trace(self) -> float:
"""Return the sum of the diagonal.
Returns:
``M[0][0] + M[1][1] + M[2][2]``.
Example:
>>> Mat3.identity().trace()
3.0
"""
return self.row0.x + self.row1.y + self.row2.z
[docs]
def determinant(self) -> float:
r"""Return the determinant.
:math:`\det M = r_0 \cdot (r_1 \times r_2)` for rows :math:`r_i`.
In plain language: the signed volume of the box spanned by the rows. Negative
means the matrix mirrors — a handedness hazard (PRD 7.1) wherever this matrix
is used as a frame.
Returns:
``det(M)``.
Example:
>>> scaling = Mat3.from_rows(((2.0, 0.0, 0.0), (0.0, 3.0, 0.0), (0.0, 0.0, 4.0)))
>>> scaling.determinant()
24.0
"""
return self.row0.dot(self.row1.cross(self.row2))
[docs]
def frobenius_norm(self) -> float:
"""Return the Frobenius norm, the square root of the sum of squared entries.
Used to scale convergence thresholds and test tolerances, since the absolute
error of an eigendecomposition scales with the size of the matrix.
Returns:
``||M||_F``.
Example:
>>> Mat3.identity().frobenius_norm() == math.sqrt(3.0)
True
"""
return math.sqrt(
self.row0.norm_squared() + self.row1.norm_squared() + self.row2.norm_squared()
)
[docs]
def is_symmetric(self, *, tolerance: float = 1e-12) -> bool:
"""Report whether the matrix equals its transpose to a relative tolerance.
Args:
tolerance: Relative tolerance. The comparison is
``|M[i][j] - M[j][i]| <= tolerance * max(1, ||M||_F)``. Default 1e-12,
about 4500 times machine epsilon: loose enough for a matrix assembled
as ``A^T A`` from real data, where the two triangles are computed by
different summation orders, and tight enough that a genuinely
asymmetric matrix is rejected.
Returns:
``True`` if symmetric within the tolerance.
Example:
>>> symmetric = Mat3.from_rows(((1.0, 2.0, 3.0), (2.0, 4.0, 5.0), (3.0, 5.0, 6.0)))
>>> symmetric.is_symmetric()
True
"""
scale = max(1.0, self.frobenius_norm())
limit = tolerance * scale
return (
abs(self.row0.y - self.row1.x) <= limit
and abs(self.row0.z - self.row2.x) <= limit
and abs(self.row1.z - self.row2.y) <= limit
)
[docs]
def eigen_sym3(
matrix: Mat3,
*,
symmetry_tolerance: float = 1e-12,
) -> tuple[tuple[float, float, float], tuple[Vec3, Vec3, Vec3]]:
r"""Return the eigenvalues and eigenvectors of a real symmetric 3x3 matrix.
Uses the cyclic Jacobi method: repeatedly apply a plane rotation that annihilates
the largest remaining off-diagonal entry, until the off-diagonal mass is negligible.
For a symmetric matrix the eigenvalues are real and the eigenvectors orthogonal, so
the accumulated rotations are the eigenvector basis.
The rotation angle for the ``(p, q)`` entry comes from
:math:`\theta = (a_{qq} - a_{pp}) / (2 a_{pq})`, with
:math:`t = \operatorname{sign}(\theta) / (|\theta| + \sqrt{\theta^2 + 1})`.
In plain language: find the three perpendicular directions along which the matrix
only stretches, and by how much.
Hazard:
**Overflow in the Jacobi rotation (PRD 11.13).** When the off-diagonal entry is
tiny relative to the diagonal gap — exactly what a **near-circular
cross-section** produces (PRD 11.10) — the textbook rotation formula overflows
silently and the rotation collapses to the identity, so the sweep stops making
progress. This routine is guarded in two independent places. The tangent comes
from ``_jacobi_tangent``, which takes a first-order branch before overflow
is possible, and the annihilated entry is written to zero explicitly rather
than left to the rotation, so even a lost rotation cannot stall the sweep. Do
not remove either guard. The bug this describes was caught by the numpy oracle
and by nothing else (PRD 13.3 T14); ``tests/unit/test_linalg.py`` holds the
regression.
Args:
matrix: A real symmetric 3x3 matrix. Only the upper triangle is read; the lower
triangle is checked against it and otherwise ignored. Unit-agnostic — for a
second-moment matrix of points in millimeters the eigenvalues are in mm^2.
symmetry_tolerance: Relative tolerance passed to
:meth:`Mat3.is_symmetric`. Default 1e-12; see that method for the
justification.
Returns:
A pair ``(eigenvalues, eigenvectors)``:
* ``eigenvalues`` — three floats in **descending** order.
* ``eigenvectors`` — three orthonormal :class:`Vec3`, ``eigenvectors[i]``
belonging to ``eigenvalues[i]``.
The basis is **right-handed**: ``eigenvectors[2] == eigenvectors[0].cross(
eigenvectors[1])`` and the determinant of the three as columns is ``+1``. This
is the convention PRD 4 fixes, and it is guaranteed here rather than left to
the caller because a mirrored frame is PRD 7.1's first hazard — a mirrored
template looks perfect and fits nothing, and the Jacobi accumulation on its own
yields a left-handed basis for about half of all inputs.
The result is also **deterministic**: the signs of ``eigenvectors[0]`` and
``eigenvectors[1]`` are pinned so that each one's largest-magnitude component
is positive, and the third follows from the first two. An eigenvector is only
defined up to sign, so pinning it costs nothing and it is what makes the export
golden files byte-comparable (ADR-012, PRD 13.3 T10).
Note that both guarantees are about the *basis*, not about which eigenvector is
which: the pairing with ``eigenvalues`` is exact either way.
Raises:
ValueError: If ``matrix`` is not symmetric within ``symmetry_tolerance``, or if
any entry is not finite. Both are programming errors upstream, and both
would otherwise yield a confident wrong answer.
ArithmeticError: If the sweeps do not converge. This cannot happen for a finite
symmetric 3x3 input and indicates a defect in this function.
Example:
>>> diagonal = Mat3.from_rows(((3.0, 0.0, 0.0), (0.0, 1.0, 0.0), (0.0, 0.0, 2.0)))
>>> values, vectors = eigen_sym3(diagonal)
>>> values
(3.0, 2.0, 1.0)
>>> vectors[0]
Vec3(x=1.0, y=0.0, z=0.0)
References:
Golub, G. H. and Van Loan, C. F., *Matrix Computations*, 4th ed., Johns Hopkins
University Press, 2013, section 8.5 (the cyclic and threshold Jacobi methods).
This is the reference for the method, the rotation formulas given above and the
stable form of the tangent in ``_jacobi_tangent``.
The "negligible after four sweeps" cutoff below is the threshold-Jacobi idea in
the specific form popularized by Press, W. H. et al., *Numerical Recipes*,
section 11.1, and is restated here from that description. **No code from that
book is reproduced or adapted** — it is not redistributable, so this matters.
The routine below is written for the 3x3 case directly: it accumulates no
``b``/``z`` correction vectors, has no ``tresh`` schedule, tests convergence
against a Frobenius-scaled threshold, and enumerates the three index pairs
rather than looping over ``n``.
"""
for row in matrix.rows():
for component in row.as_tuple():
if not math.isfinite(component):
msg = f"eigen_sym3 requires a finite matrix, got entry {component!r}"
raise ValueError(msg)
if not matrix.is_symmetric(tolerance=symmetry_tolerance):
msg = (
"eigen_sym3 requires a symmetric matrix; "
f"got rows {matrix.row0}, {matrix.row1}, {matrix.row2}"
)
raise ValueError(msg)
# Working copy of the upper triangle plus diagonal, and the accumulated rotations.
a = [
[matrix.row0.x, matrix.row0.y, matrix.row0.z],
[matrix.row0.y, matrix.row1.y, matrix.row1.z],
[matrix.row0.z, matrix.row1.z, matrix.row2.z],
]
v = [
[1.0, 0.0, 0.0],
[0.0, 1.0, 0.0],
[0.0, 0.0, 1.0],
]
# Convergence threshold: the off-diagonal mass is negligible once it falls below
# the rounding level of the diagonal. Scaling by the Frobenius norm makes the test
# invariant to the units the caller works in.
scale = matrix.frobenius_norm()
if scale == 0.0:
return (0.0, 0.0, 0.0), (
Vec3(1.0, 0.0, 0.0),
Vec3(0.0, 1.0, 0.0),
Vec3(0.0, 0.0, 1.0),
)
converged_below = _EPS * scale
# `a` is held fully symmetric throughout, so a[i][j] and a[j][i] always agree and
# either index may be read.
pairs = ((0, 1), (0, 2), (1, 2))
for sweep in range(_MAX_SWEEPS):
off_diagonal = abs(a[0][1]) + abs(a[0][2]) + abs(a[1][2])
if off_diagonal <= converged_below:
break
for p, q in pairs:
apq = a[p][q]
if apq == 0.0:
continue
# Threshold-Jacobi cutoff (see the References in the docstring): after four
# sweeps, an off-diagonal entry that is lost in the rounding of both its
# diagonals contributes nothing and is set to zero outright. This is what
# guarantees the loop terminates rather than chasing a residual that
# round-off keeps regenerating.
negligible = 100.0 * abs(apq)
if (
sweep >= 4
and abs(a[p][p]) + negligible == abs(a[p][p])
and abs(a[q][q]) + negligible == abs(a[q][q])
):
a[p][q] = 0.0
a[q][p] = 0.0
continue
# PRD 11.13 hazard — see _jacobi_tangent and the docstring above.
t = _jacobi_tangent(apq, a[q][q] - a[p][p])
c = 1.0 / math.sqrt(1.0 + t * t)
s = t * c
tau = s / (1.0 + c)
shift = t * apq
a[p][p] -= shift
a[q][q] += shift
a[p][q] = 0.0
a[q][p] = 0.0
# Rotate the one remaining off-diagonal entry: the index not in (p, q).
(r,) = {0, 1, 2} - {p, q}
arp = a[r][p]
arq = a[r][q]
new_rp = arp - s * (arq + tau * arp)
new_rq = arq + s * (arp - tau * arq)
a[r][p] = new_rp
a[p][r] = new_rp
a[r][q] = new_rq
a[q][r] = new_rq
for row_index in range(3):
vp = v[row_index][p]
vq = v[row_index][q]
v[row_index][p] = vp - s * (vq + tau * vp)
v[row_index][q] = vq + s * (vp - tau * vq)
else:
msg = (
f"eigen_sym3 did not converge in {_MAX_SWEEPS} sweeps; "
"this is a defect in eigen_sym3, not in the input"
)
raise ArithmeticError(msg)
eigenvalues = (a[0][0], a[1][1], a[2][2])
eigenvectors = tuple(Vec3(v[0][k], v[1][k], v[2][k]) for k in range(3))
order = sorted(range(3), key=lambda k: eigenvalues[k], reverse=True)
sorted_values = (
eigenvalues[order[0]],
eigenvalues[order[1]],
eigenvalues[order[2]],
)
# HAZARD PRD 7.1 (handedness). The Jacobi accumulation produces an orthonormal
# basis with no orientation guarantee, and fixing each vector's sign independently
# does not restore one: measured over 2000 random symmetric matrices, roughly half
# the results came out left-handed (det = -1). A caller that turns these axes into
# a cross-section frame (Tier C, PRD 6.4) or a Bishop-frame seed (Tier D, PRD 6.5)
# would then get a mirrored frame on half its inputs, and inconsistently between
# faces of one model -- a mirrored template looks perfect and fits nothing.
#
# So the sign rule is applied to the first two vectors only, and the third is taken
# as their cross product. That is still exactly +/- the third eigenvector, because
# the three are orthonormal and an eigenvector's sign is free, so the
# eigenvalue/eigenvector pairing is untouched. It is fully deterministic, which the
# byte-comparable exports require (ADR-012, PRD 13.3 T10), and it gives det = +1
# unconditionally.
first = _canonical_sign(eigenvectors[order[0]])
second = _canonical_sign(eigenvectors[order[1]])
# Renormalized because the cross product of two vectors that are orthonormal only
# to within round-off has a norm of 1 +/- a few ulps, and callers are promised unit
# vectors. The direction is unaffected.
third = first.cross(second).normalized()
return sorted_values, (first, second, third)
def _jacobi_tangent(off_diagonal: float, gap: float) -> float:
r"""Return the tangent of the Jacobi rotation that annihilates ``off_diagonal``.
The rotation angle :math:`\phi` for the ``(p, q)`` entry satisfies
:math:`\theta = (a_{qq} - a_{pp}) / (2 a_{pq})` and
:math:`t = \tan\phi = \operatorname{sign}(\theta) /
(|\theta| + \sqrt{\theta^2 + 1})`, the root of :math:`t^2 + 2\theta t - 1 = 0`
with the smaller magnitude. Choosing the smaller root keeps the rotation under 45
degrees, which is what makes cyclic Jacobi stable.
In plain language: how far to twist the ``(p, q)`` plane so that the coupling
between axes ``p`` and ``q`` disappears.
Hazard:
**Silent overflow (PRD 11.13).** When ``|gap|`` dwarfs ``|off_diagonal|``,
:math:`\theta` exceeds 1.3e154 and :math:`\theta^2` overflows to infinity.
Python floats overflow silently — no exception, no warning — so ``t`` collapses
to exactly ``0.0``: the rotation becomes the identity and the sweep makes no
progress. That input is not exotic; it is what a **near-circular cross-section**
produces, such as the 12.50 mm near-circular end of the reference model's
morphing tube (PRD 11.10). The
guard below switches to the first-order form :math:`t \approx a_{pq} / gap`
whenever ``100 * |off_diagonal|`` is lost in the rounding of ``|gap|``, which
happens at :math:`|\theta| \approx 2.3\mathrm{e}{17}` — safely below the
overflow threshold. ``tests/unit/test_linalg.py`` holds the regression.
Args:
off_diagonal: The entry ``a[p][q]`` to annihilate. Must be non-zero; a zero
entry needs no rotation and the caller skips it.
gap: The diagonal difference ``a[q][q] - a[p][p]``.
Returns:
The tangent ``t`` of the rotation angle, with ``|t| <= 1``.
Example:
A gap of 2 with an off-diagonal of 1e-200 puts theta at 1e200, where the
textbook formula overflows. The guard returns the first-order value instead:
>>> _jacobi_tangent(1e-200, 2.0)
5e-201
An equal-magnitude gap and off-diagonal give the full 45-degree rotation:
>>> _jacobi_tangent(1.0, 0.0)
1.0
"""
if abs(gap) + 100.0 * abs(off_diagonal) == abs(gap):
return off_diagonal / gap
theta = 0.5 * gap / off_diagonal
tangent = 1.0 / (abs(theta) + math.sqrt(1.0 + theta * theta))
return -tangent if theta < 0.0 else tangent
def _canonical_sign(vector: Vec3) -> Vec3:
"""Flip ``vector`` so its largest-magnitude component is positive.
An eigenvector is defined only up to sign. Pinning the sign makes
:func:`eigen_sym3` deterministic, which the byte-comparable exports require
(ADR-012, PRD 13.3 T10).
Args:
vector: The vector to canonicalize.
Returns:
Either ``vector`` or ``-vector``.
Example:
>>> _canonical_sign(Vec3(0.25, -1.0, 0.5))
Vec3(x=-0.25, y=1.0, z=-0.5)
"""
components = vector.as_tuple()
dominant = max(range(3), key=lambda k: abs(components[k]))
return -vector if components[dominant] < 0.0 else vector