Struct Mat3Fix 

#[repr(C)]
pub struct Mat3Fix {
    pub col0: Vec3Fix,
    pub col1: Vec3Fix,
    pub col2: Vec3Fix,
}

3x3 Matrix for inertia tensors and rotations

Fields

col0: Vec3Fix

Column 0

col1: Vec3Fix

Column 1

col2: Vec3Fix

Column 2

Implementations

impl Mat3Fix

pub const IDENTITY: Self

Identity matrix

pub const ZERO: Self

Zero matrix

pub const fn from_cols(col0: Vec3Fix, col1: Vec3Fix, col2: Vec3Fix) -> Self

Create from columns

pub const fn diagonal(x: Fix128, y: Fix128, z: Fix128) -> Self

Create diagonal matrix

pub fn mul_vec(self, v: Vec3Fix) -> Vec3Fix

Matrix-vector multiplication

pub const fn transpose(self) -> Self

Transpose

pub fn scale(self, s: Fix128) -> Self

Scale all elements

pub fn mul_mat(self, rhs: Self) -> Self

Matrix-matrix multiplication (self * rhs)

pub fn determinant(self) -> Fix128

Determinant

pub fn inverse(self) -> Option<Self>

Inverse matrix.

Returns None if the matrix is singular or if the inverse cannot be computed inside the Q64.64 range: a cofactor, the determinant, 1 / det (|det| ≤ 2⁻⁶³) or one of the entries cofactor / det does not fit. Earlier versions wrapped in the out-of-range cases and returned Some with flipped signs (e.g. diag(2²¹, 2²¹, 2²¹), whose determinant is 2⁶³).

The products inside each cofactor and inside the determinant are summed exactly (256-bit), so a single product beyond 2⁶³ that cancels out does not make the result None: in that case the wrapping arithmetic was already exact modulo 2¹²⁸, and the result is bit-identical to the unchecked formula, as it is for every input whose intermediates fit.

A matrix whose entries are all small (largest below 1/2) and whose determinant is below 2⁻³² is scaled up by a power of two first, so its determinant is not quantised away (diag(1e-5, 1e-5, 1e-5) had a relative error of 5e-5, and below 2⁻⁶³ the result was None although the inverse fits): A⁻¹ = 2ᵏ (2ᵏ A)⁻¹, both scalings exact. Every other matrix gets the unscaled result, bit for bit. k follows the largest entry, so the gain is for matrices whose entries are of similar size: one with entries many orders apart (diag(3·2⁻⁶⁴, 2⁻³³, 3·2⁻⁶⁴)) still has its determinant quantised and can be off by a large relative error (12.5 % there), where it used to be None.

pub fn max_abs_component(self) -> Fix128

Largest absolute component, as a scale for this matrix.

Useful for turning an intrinsic, scale-free criterion into an absolute one at the call site: a determinant scales as the cube of the matrix, so ‖F‖³ is the quantity a determinant floor has to be expressed against. See Self::polar_rotation.

pub fn polar_rotation( self, det_floor: Fix128, max_iters: u32, ) -> Result<Self, PolarError>

Rotation factor R of the polar decomposition F = R·U, where U is symmetric positive definite.

This is the piece a co-rotational formulation needs: the rotation a deformation gradient contains, separated from the stretch, so an element’s small-strain stiffness can be rotated into the deformed frame instead of reading a large rotation as a strain.

Method

Higham’s Newton iteration for the orthogonal polar factor,

R₀ = F,   R_{k+1} = ½ (R_k + R_k⁻ᵀ)

which needs only Self::inverse, Self::transpose, addition and a halving — no square root, no trigonometry, no eigensolver. Every operation is one Fix128 performs exactly, so two runs on different targets produce bit-identical rotations and nothing here can reach a platform libm. Nothing allocates, so this is available in no_std builds.

Convergence, measured

Quadratic near the answer — the number of correct bits doubles per step — and linear while a singular value is far from one, where the iteration halves the distance per step. Measured change per step for a general gradient:

4.17e-1 → 8.01e-2 → 3.20e-3 → 5.12e-6 → 1.31e-11 → 0

Measured step counts to the settling bound (tests/analytic_polar_decomposition.rs prints the first four; the rest were measured by raising max_iters until the call succeeded):

F	steps
simple shear γ = 0.4	6
diag(1.5, 0.8, 1.1)	7
rotation ∘ diag(1.3, 0.9, 1.05)	7
a general full matrix	7
diag(1, 1, 0.1)	9
diag(1e3, 1e3, 1e2)	9
diag(4, 1, 0.25)	10
diag(1e-3, 1e-3, 1e-4)	10
diag(1, 1, 0.01)	13
diag(1, 1, 1e-4)	19

So max_iters = 32 covers everything measured with margin, and 16 is not enough — the flattest gradient the intrinsic floor still accepts needs 19.

The two uniformly scaled rows are the point of the normalisation below: diag(1e3, 1e3, 1e2) and diag(1e-3, 1e-3, 1e-4) are the same shape six decades apart and both finish in 9 or 10 steps, so the budget a caller has to pass does not depend on the units its field is expressed in. What the count does depend on is the spread: diag(1, 1, 1e-4) is slow because one singular value is four decades from the others, and no uniform scaling can fix that.

An exactly representable rotation is returned unchanged: its largest component lies in [1/√3, 1] so the divisor is exactly one, and for orthogonal F, F⁻ᵀ = F, so the first step is ½(F + F) with no rounding.

Why the iteration stops on a change of a few bits, not on no change

⚠️ There is in general no matrix this iteration leaves unchanged. Fix128 multiplication truncates rather than rounding to nearest, so the rounded map carries a systematic bias of a fraction of a unit in the last place per step, and the iterate keeps creeping after it has converged to everything the representation can hold. Measured on a rotation about x built from 3/5 and 4/5 (whose determinant is 1 + 819·2⁻⁶⁴, not exactly one), tracked for 24 steps:

step  1: col1.y lo=…730314   col1.z lo=…641784
step  4: col1.y lo=…730312   col1.z lo=…641785
step  9: col1.y lo=…730310   col1.z lo=…641786
step 14: col1.y lo=…730309   col1.z lo=…641787   ← no fixed point, no cycle

So next == r is not a reachable stopping condition, and neither is a period-two cycle: both were tried and both burned the whole budget, which the caller then received as a failure to decompose. The drift is 1 unit in the last place per one or two steps; the iteration stops at a change of 4, which is above the drift with margin and far below anything a caller of a rotation matrix can observe.

Rounding to nearest would not remove the need for this rule — it would make the bias unbiased rather than systematic, and the map would still be a discretisation with no guaranteed fixed point.

The bound is scale free because the iterate is normalised: ‖r‖ ≈ 1, so “four units in the last place” is a relative quantity, not the absolute threshold that always breaks at some other problem size.

Why det_floor is a parameter and not a constant

A matrix does not know the scale of the field it came from, so a floor written inside this function would be an absolute threshold — the shape that always breaks at some problem size. The policy belongs to the caller: a finite-element assembly derives it from the element volume, and a caller that only wants the intrinsic numerical criterion passes

det_floor = ε · ‖F‖³      (Self::max_abs_component gives ‖F‖)

which is scale free, because both sides scale as the cube of the matrix.

Returns None

• When det F ≤ 0, whatever det_floor is given. A non-positive determinant means the deformation has turned the element inside out, and the nearest orthogonal matrix is then a reflection with det = −1. Returning it would silently flip the sign of every stress component, so this rejection is not delegated to the caller’s floor: a caller that passes a negative det_floor still cannot get a reflection back.
• When det F ≤ det_floor.
• When max_iters is exhausted, or an intermediate matrix is exactly singular.

Errors

One of the four PolarError variants, which are kept separate because the caller’s next action differs for each: PolarError::Inverted and PolarError::Degenerate are about the data, PolarError::NotConverged is about the budget, and PolarError::SingularStep is about the path.

Trait Implementations

impl Clone for Mat3Fix

fn clone(&self) -> Self

Returns a duplicate of the value.

fn clone_from(&mut self, source: &Self)

Performs copy-assignment from source.

impl Copy for Mat3Fix

impl Debug for Mat3Fix

fn fmt(&self, f: &mut Formatter<'_>) -> Result

Formats the value using the given formatter.

impl Default for Mat3Fix

fn default() -> Self

Returns the “default value” for a type.

impl Eq for Mat3Fix

impl Mul for Mat3Fix

type Output = Mat3Fix

The resulting type after applying the * operator.

fn mul(self, rhs: Self) -> Self

Performs the * operation.

impl PartialEq for Mat3Fix

fn eq(&self, other: &Self) -> bool

Equality operator ==.

fn ne(&self, other: &Rhs) -> bool

Inequality operator !=.

impl StructuralPartialEq for Mat3Fix