# onemath docs, all pages (OneMath 1.0.0)

Every page of https://recursiveai.org/onemath/docs/ and https://recursiveai.org/onemath/cookbook/ as markdown, getting started first. Index: https://recursiveai.org/onemath/llms.txt


<!-- https://recursiveai.org/onemath/docs/getting-started/ -->

# Getting started with OneMath

Install OneMath for .NET 10 and use its correctly rounded functions, vectors and matrices, exact predicates and fixed-point bridge, with results that are the same bits on every platform.

## Install

OneMath targets `net10.0`, so you need the .NET 10 SDK. This guide uses 1.0.0.

```shell
dotnet new console -o OneMathStart
cd OneMathStart
dotnet add package RecursiveAI.OneMath --version 1.0.0
```

The package has no dependency beyond the .NET base class library and works under NativeAOT, trimming and browser WebAssembly.

## A first program

```csharp
using RecursiveAI.OneMath;
using RecursiveAI.OneMath.Predicates;

// Elementary functions: correctly rounded, the same bits everywhere.
double s = Math64.Sin(1e300);
float e = Math32.Exp(1.5f);
Console.WriteLine($"Math64.Sin(1e300) = {s:R}  bits {BitConverter.DoubleToInt64Bits(s):X16}");
Console.WriteLine($"Math32.Exp(1.5f)  = {e:R}");

// Vectors, quaternions and matrices, in float (F) and double (D).
var q = QuaternionD.CreateFromAxisAngle(new Vector3D(0, 1, 0), Math.PI / 2);
Matrix4x4D m = Matrix4x4D.CreateFromQuaternion(q) * Matrix4x4D.CreateTranslation(new Vector3D(10, 0, 0));
Vector3D p = m.TransformPoint(new Vector3D(1, 0, 0));
Console.WriteLine($"Rotated and moved: {p}");

// Exact predicates: the sign of the exact determinant.
int side = ExactPredicates.Orient2D(0, 0, 1, 0, 0.5, 1e-10);
Console.WriteLine($"Orient2D: {side} (+1 left, -1 right, 0 on the line)");

// The Q31.32 fixed-point bridge.
long raw = FixedBridge.FromDouble(1.5);
Console.WriteLine($"FixedBridge.FromDouble(1.5) = {raw} = 0x{raw:X}, back: {FixedBridge.ToDouble(raw)}");
```

```text
Math64.Sin(1e300) = -0.8178819121159085  bits BFEA2C16B010E385
Math32.Exp(1.5f)  = 4.481689
Rotated and moved: (2.4424906541753444E-15, 0, -11)
Orient2D: 1 (+1 left, -1 right, 0 on the line)
FixedBridge.FromDouble(1.5) = 6442450944 = 0x180000000, back: 1.5
```

What each part is:

- **`Math64` and `Math32`** have sin, cos, tan, sincos, asin, acos, atan, atan2, exp, exp2, expm1, log, log2, log10, log1p, pow, cbrt, hypot, sinh, cosh and tanh for `double` and `float`. Each is correctly rounded: the result is the representable number nearest to the exact value, even for an argument like 10^300. See [elementary functions](https://recursiveai.org/onemath/docs/functions/).
- **The types** are `Vector2/3/4`, `Matrix2x2/3x3/4x4`, `Quaternion`, `Affine3`, `RigidTransform` and `Bounds2/3`, each with an `F` (float) and a `D` (double) form. Matrices use **column vectors**: `R * T` applies `T` first. That is the opposite of System.Numerics; see [Port from System.Numerics](https://recursiveai.org/onemath/cookbook/port-from-system-numerics/).
- **`ExactPredicates`** (namespace `RecursiveAI.OneMath.Predicates`) returns the sign of the exact determinant, never a rounding accident: see [Exact geometric predicates](https://recursiveai.org/onemath/cookbook/exact-predicates/).
- **`FixedBridge`** converts between Q31.32 fixed-point raw values and `double`, exactly within ±2^21: see [The fixed-point bridge](https://recursiveai.org/onemath/cookbook/fixed-point-bridge/).

The rotated point shows double rounding at work: `2.44E-15` is the cosine of a rounded π/2, not zero. What OneMath promises is that it is the same `2.44E-15` everywhere.

## The determinism promise

A result depends only on its inputs, not on the OS, the CPU and its SIMD or fused multiply-add support, the runtime version, the JIT tier, NativeAOT or WebAssembly, or the thread count. OneMath never calls `Math.Sin` and friends and never uses fused multiply-add.

It also holds across versions: no 1.x release changes the result of any existing function, type or predicate. A deterministic host can update within 1.x without breaking replays, saved games or lockstep peers. [Check that two machines agree](https://recursiveai.org/onemath/cookbook/check-machines-agree/) shows how to test it on your own targets.

## Next

- [Exact geometric predicates](https://recursiveai.org/onemath/cookbook/exact-predicates/)
- [Port from System.Numerics](https://recursiveai.org/onemath/cookbook/port-from-system-numerics/)
- [Inertia tensors and principal axes](https://recursiveai.org/onemath/cookbook/inertia-tensors/)
- [Evaluation order](https://recursiveai.org/onemath/docs/evaluation-order/): the order of every operation, so another implementation can reproduce the bits.


<!-- https://recursiveai.org/onemath/docs/overview/ -->

# onemath

Deterministic math for .NET (`RecursiveAI.OneMath`): vectors, matrices, quaternions, elementary functions and exact geometric predicates that give the same bits on every OS, CPU, runtime, WebAssembly mode and thread count. Only IEEE-exact `+ - * /` and `sqrt` decide a result; every other function is implemented here. No dependencies beyond the .NET base class library; NativeAOT and trimming compatible.

**Status: 1.0.0.** The public API is frozen under semantic versioning. Every result that exists in `0.1.0-preview.5` is bit-identical in 1.0.0; the additions are listed in the changelog. See [CHANGELOG.md](https://recursiveai.org/onemath/docs/changelog/).

```
dotnet add package RecursiveAI.OneMath
```

```csharp
using RecursiveAI.OneMath;
using RecursiveAI.OneMath.Predicates;

double s = Math64.Sin(1e300);                                  // correctly rounded, same on every platform
int side = ExactPredicates.Orient2D(0, 0, 1, 0, 0.5, 1e-10);   // sign of the exact determinant: +1, -1 or 0
long raw = FixedBridge.FromDouble(1.5);                        // Q31.32 fixed point, exact within +-2^21
```

## What is in it

- **Types:** `Vector2/3/4`, `Matrix2x2/3x3/4x4`, `Quaternion`, `Affine3`, `RigidTransform`, `Bounds2/3`, each in float (`F`) and double (`D`) form. The evaluation order of every operation is written down in [docs/evaluation-order.md](https://recursiveai.org/onemath/docs/evaluation-order/), so another implementation can reproduce the bits.
- **Elementary functions:** `Math32` and `Math64` provide sin, cos, tan, sincos, asin, acos, atan, atan2, exp, exp2, expm1, log, log2, log10, log1p, pow, cbrt, hypot, sinh, cosh and tanh. They are **correctly rounded**: the result is the float or double nearest to the exact value. How this is checked, and the error bounds of each internal stage, are in [docs/functions.md](https://recursiveai.org/onemath/docs/functions/).
- **Exact predicates:** `Orient2D`, `Orient3D`, `InCircle`, `InSphere` and segment classification, adaptive arithmetic whose fast filter always returns the sign of the exact value ([docs/predicates.md](https://recursiveai.org/onemath/docs/predicates/)). Coordinates must be finite and either zero or of magnitude in [2^-140, 2^200]; **out-of-range input may throw** `ArgumentOutOfRangeException`.
- **Fixed-point bridge:** conversion between a deterministic host's Q31.32 integers and double, exact for magnitudes below 2^21 and rejecting anything beyond ([docs/fixed-bridge.md](https://recursiveai.org/onemath/docs/fixed-bridge/)).
- **Physics math:** symmetric 3x3 eigen decomposition and inertia tensors ([docs/physics.md](https://recursiveai.org/onemath/docs/physics/)).

## Matrix convention

Matrices use **column vectors**: `m * v` transforms the column `v`, the translation of a 4x4 sits in column 4 (`M14`, `M24`, `M34`), and `a * b` applies `b` first. **System.Numerics is the opposite**: row vectors, translation in `M41..M43`, and `a * b` applies `a` first, although the field names `Mrc` are the same. Porting by copying fields, or keeping operand order, silently gives a different transform. Translate then rotate:

```csharp
// System.Numerics (row vectors): T first, then R
var sn = System.Numerics.Matrix4x4.CreateTranslation(10, 0, 0) * System.Numerics.Matrix4x4.CreateRotationY(angle);
// OneMath (column vectors): the same transform is R * T
var m = Matrix4x4F.CreateRotationY(angle) * Matrix4x4F.CreateTranslation(new Vector3F(10, 0, 0));
// Convert instead of copying fields: both directions transpose
Matrix4x4F same = Matrix4x4F.FromSystemNumerics(sn);   // same transform as m
System.Numerics.Matrix4x4 back = m.ToSystemNumerics();
```

Vectors and quaternions convert with `FromSystemNumerics` / `ToSystemNumerics` too (quaternions keep their components and their Hamilton product). Converting from double to float is only available as the explicit `ToSystemNumerics()` and rounds to nearest even.

## The determinism promise

A result depends only on its inputs: not on the OS, the CPU or its features (SIMD or scalar, with or without FMA hardware), the runtime version, the JIT tier, NativeAOT or WebAssembly, or the thread count. The library never calls `Math.Sin` and friends and never uses fused multiply-add; a banned-symbol analyzer enforces this at build time.

**Across versions.** Within a major version, every release gives the same bits for the same inputs: no 1.x minor or patch release changes the result of any existing function, type or predicate. A change that would alter an existing result, even a more accurate one, either ships as a new member that leaves the old one as it was, or waits for the next major version and is listed in the changelog. Members added in a minor release then keep their results for the rest of that major version. So a deterministic host can update within 1.x without breaking replays, saved games or lockstep peers that run a different 1.x release.

The promise is tested by a hash matrix: [tools/HashMatrix](https://github.com/recursive-ai-oss/onemath-csharp/tree/1986f32ecb9f5f3001c0121484f104c489720334/tools/HashMatrix) runs fixed scenarios and writes a SHA-256 manifest, and the manifests must be equal across x64 JIT, scalar (`DOTNET_EnableHWIntrinsic=0`), `DOTNET_EnableAVX2=0`, one thread and all cores, NativeAOT, browser-wasm (interpreted and AOT-compiled) and osx-arm64. See [docs/hash-matrix.md](https://recursiveai.org/onemath/docs/hash-matrix/). Across versions, every line of the previous release's manifest must appear unchanged in the next release's; a minor release may only add lines.

## Speed

Correct rounding costs time. Informal single-thread measurements on a shared development machine (ours divided by the platform's `Math`; not a gate measurement): sin 2.8x, cos 4.3x, tan 4.9x, atan2 5.5x, exp 12x, log 15x, pow 16x, cosh 21x. The full table and how to reproduce it are in [docs/functions.md](https://recursiveai.org/onemath/docs/functions/).

## Platforms

`net10.0` on Linux (x64, arm64), macOS (arm64), Windows (x64) and browser-wasm. CI builds and tests on the four native targets, publishes the NativeAOT smoke app with warnings as errors, and compares the hash-matrix manifest on WebAssembly.

## Layout

- `src/OneMath/` the library
- `tests/OneMath.Tests/` xUnit v3 tests; `tests/AotSmoke/` and `tests/WasmSmoke/` publish-and-run checks
- `tools/HashMatrix/` scenario runner; `tools/hash-matrix.sh` runs it across platforms
- `bench/Micro/` BenchmarkDotNet; `docs/` normative notes

## Build

```
dotnet build -c Release
dotnet test -c Release
dotnet build -c Checked     # Release optimisations plus the CHECKED define
```

Build output goes to `$ONEMATH_ARTIFACTS` (or `~/.cache/onemath-artifacts`), never into the tree.

## Licence

Licensed under either of [MIT](https://github.com/recursive-ai-oss/onemath-csharp/tree/1986f32ecb9f5f3001c0121484f104c489720334/LICENSE-MIT) or [Apache-2.0](https://github.com/recursive-ai-oss/onemath-csharp/tree/1986f32ecb9f5f3001c0121484f104c489720334/LICENSE-APACHE), at your option. Copyright Recursive AI, LLC.


<!-- https://recursiveai.org/onemath/docs/functions/ -->

# Elementary functions

`Math64` (double) and `Math32` (float32) in `RecursiveAI.OneMath`: `Sin Cos Tan SinCos Asin Acos Atan Atan2 Exp Exp2 Expm1 Log Log2 Log10 Log1p Pow Cbrt Hypot Sinh Cosh Tanh`.

## Contract

- **Every function is correctly rounded for every input**: the result is the IEEE 754 round-to-nearest, ties-to-even value of the exact mathematical result, with gradual underflow and overflow to infinity. Results are therefore bit-identical on every platform and equal to any other correctly rounded implementation. Nothing in the library calls `Math.*`/`MathF.*` for a function result: only `+ - * / sqrt` on float32/double (IEEE-exact), integer arithmetic, and `System.Numerics.BigInteger` (integer arithmetic) decide a bit.
- Special values follow IEEE 754-2019 clause 9.2 and C99 Annex F (`pow(x, ±0) = pow(1, y) = 1` even for NaN, `atan2(±0, -0) = ±π`, `hypot(±∞, NaN) = +∞`, `expm1(-0) = -0`, ...). Every NaN result is the canonical `double.NaN` / `float.NaN` regardless of input payload.
- Float32 inputs are converted to double (exact) and run through the same code as double; only the final rounding differs.
- Domain errors return NaN (no exceptions). Large arguments are fine: `Sin(1e300)` is exact.

## How a result is produced: three stages and a rounding test

Each call produces a value `v` with a known relative error bound `δ`, then asks Ziv's question (Ziv 1991; Muller, *Elementary Functions*, ch. 10): do all reals in `v(1 ± δ)` round to the same output? If yes, that output is correct. If not, the next, more accurate stage runs.

1. **Fast stage** (`Fast.cs`, tables in `FastTables.cs`): table-driven, plain double arithmetic plus error-free transformations; relative error ≤ 2^-62 .. 2^-65 depending on the function (bounds below). Decides 98–99.8% of double results (measured per function below) and all but a few in 10^6 of float32 results.
2. **Double-double stage** (`Cores.cs`, `DD.cs`): relative error ≈ 2^-100 (measured: `pow` 2^-96, `exp`/`exp2` 2^-106, the rest 2^-101 to 2^-105). Taylor-type series with exact-rational coefficients computed at start-up by double-double division. Decides all but about 1 double in 2^31 (and effectively no float32).
3. **Exact stage** (`Exact/`): about 500 bits using `BigInteger`, then one rounding. Used when the double-double test fails, when a double result is subnormal or at the overflow edge, and for `hypot`/`cbrt` ties (exact integer roots with a sticky bit). This stage shares no code with the other two, and it is what the tests treat as ground truth, alongside an independent mpmath (Python) run for the known-answer file.

Rounding test, double result from a double-double `(hi, lo)` with bound `δ`: `err = |hi| δ`; `a = hi + (lo + err)`, `b = hi + (lo - err)`; the result is decided iff `a == b` (monotonicity of rounding), and equals `a`. The result must be a normal double, otherwise stage 3 runs. Float32: `f = (float) hi`, `d = (hi - f) + lo`, compared with half the distance to the neighbouring float in the direction of `d`, with margin `err + 2.3e-16 |d|`; inconclusive or overflowing results go to the next stage.

The declared δ for stages 1 and 2 are the *inputs to the test*; they are validated, not proved, as follows: `FastStageTests` compares stage 1 with stage 2 on ≥ 4·10^6 inputs per function and requires every sample's error to be at most a quarter of its δ (measured margins in the table below); `Double_random` compares stage 2 with stage 3 on 2·10^7 inputs per function. A violated δ could only produce a wrong rounding when the exact result lies within the violation of a rounding boundary, so the exhaustive float32 run and the exact-checked double samples are the backstop.

## Argument handling

- **Representation.** `DD` is an unevaluated sum `Hi + Lo`, |Lo| ≤ ulp(Hi)/2. `TwoSum` is Knuth's, `TwoProd` is Dekker/Veltkamp splitting with constant 2^27+1 (no fused multiply-add; inputs below 2^996 in magnitude — callers scale first). `+`: accurate variant (relative error < 2^-103); `*`: < 2^-103; `/`: three quotient digits (< 2^-102); `sqrt`: one Karp–Markstein step.
- **Results as `m · 2^k`.** Cores return a double-double mantissa and a separate integer exponent, so the rounding test runs on a value near 1 and the scaling by `2^k` is exact. A double result that would be subnormal or overflow goes to stage 3.
- **Tiny arguments.** `sin tan asin atan sinh tanh x → x` and `cos, cosh → 1` for |x| < 2^-27; `expm1, log1p x → x` for |x| < 2^-55 (the first dropped term is at most 2^-56 of x, below the smallest rounding half-gap 2^-54 of x; this is the argument, and the known-answer and exhaustive float32 runs cross-check it). `tanh(x) = ±1` for |x| ≥ 22.
- **Thresholds** that return `±∞`/`0` directly are far outside the range where the result could still be finite (`exp`: x > 710 or x < -746; `exp2`: |x| > 1100; `sinh/cosh`: |x| > 711); anything between a threshold and the true boundary is decided by the stages.

## Methods per function

Constants (π/2 and ln 2 parts, 1/ln 2, 1/ln 10, π, 2/π bits for reduction, the stage-1 tables) are generated by `tools/GenConstants/gen_constants.py` and `gen_tables.py` with mpmath at 600–1700 bits and written as exact bit patterns (`Constants.cs`) or shortest round-trip decimal literals (`FastTables.cs`; the C# compiler parses them exactly). No table is computed at run time from a platform function. Series coefficients `1/n!`, `1/(2k+1)` are computed at start-up as `DD` quotients (`+ - * /` only, so identical everywhere).

| Function | Stage 1 (fast) | Stage 2 (double-double) |
|---|---|---|
| **exp, exp2** | `n = round(t·64/ln2)`; `r = t − n·ln2/64` with ln2/64 in three parts (two 34-bit parts, so `n·part` is exact; `r` as a double-double renormalised); `2^(n/64) = 2^(n>>6)·T[n&63]`, `T` = 64 double-double table entries; `exp(r) = 1 + r + r²/2 + … + r⁶/720` (r ≤ ln2/128; the dropped `r⁷/5040 < 2^-64`); product by `T` with the high×high product formed exactly (`TwoProd`). exp2: `r = (x − n/64)·ln2` (exact subtraction). | ln2 in three parts (40, 40, 53 bits) reduces to |r| ≤ 0.35; scale by 2^-8; Taylor series for `expm1` (5 terms in double-double, 5 in double); then 8 times `e ← 2e + e²` (relative error preserved in `expm1` form). exp2: `r = (x − round x)·ln2` in double-double. |
| **expm1** | |x| ≥ 0.5: stage-1 `exp` minus 1 (≤ 1 bit of cancellation); smaller: stage 2. | As exp but not adding 1 when k = 0, so tiny results keep full relative accuracy; exact `2^k(1+e) − 1` for |k| ≤ 600. |
| **log, log2, log10** | Not within 1/16 of 1 (those go to stage 2). `x = 2^e·m`, `m ∈ [1,2)`; `j` = top 7 mantissa bits; `inv[j] ≈ 1/(1 + (j+½)/128)`; `r = m·inv[j] − 1` exactly (TwoProd), |r| ≤ 2^-8; `log m = r − r²/2 + … + r⁷/7 + L[j]`, `L[j] = −log(inv[j])` as a double-double (table of 128); plus `e·ln2` (three parts). log2/log10 multiply by `1/ln2`, `1/ln10` as double-doubles. | `m ∈ [√½, √2)`, `z = (m−1)/(m+1)`, `log m = 2 atanh z` as a series in `z²` (11 double-double terms, 11 double terms), plus `e·ln2`. Subnormal inputs scaled by 2^54 first. |
| **log1p** | stage 2 only | `|x| < 0.2`: `2 atanh(x/(2+x))` (relative accuracy for any small x); otherwise `1+x` as an exact double-double and the log path above. |
| **sin, cos, tan, sincos** | |x| ≤ 2^19 only. Cody–Waite reduction with π/2 in parts of 33, 33, 33 and 53 bits (`k·part` exact for k < 2^20), `r` as a double-double; `a = |r|`, `j = round(a·128)` (≤ 101), `b = a − j/128` (exact), `sin a = S_j cos b + C_j sin b`, `cos a = C_j cos b − S_j sin b` with tables of `sin(j/128)`, `cos(j/128)` as double-doubles; `sin b`, `cos b − 1` from short polynomials (`b ≤ 2^-8`), the cross term of `r`'s low part included; the large product `C_j·b` is formed exactly. tan = sin/cos with a corrected double division, `−cos/sin` for odd quadrants. | Payne–Hanek reduction (below) for any |x| ≥ π/4, then Taylor series for sin and cos of `r` (8 + 7 double-double terms, 7 + 7 double terms); tan = quotient. |
| **asin, acos, atan, atan2** | `atan` of a non-negative double-double `t`: `t > 1 → 1/t` and `π/2 −`; `j = round(u·64)`; `b = (u − j/64)/(1 + u·j/64)` as a double-double; `atan u = A_j + b − b³/3 + b⁵/5 − b⁷/7 + b⁹/9` (|b| ≤ 2^-7), `A_j = atan(j/64)` double-doubles (table of 65). asin: `atan(x/√(1−x²))` with `1−x²` from an exact `x²`; acos: `2·atan(√((1−x)/(1+x)))` (no cancellation near ±1); atan2: ratio of the scaled magnitudes (smaller over larger), quadrant by `π − a`, `π/2 − a`. | Same identities with `u ← u/(1+√(1+u²))` applied twice, then the Taylor series (12 double-double + 11 double terms), result times 4. |
| **sinh, cosh, tanh** | `u = e^|x|` from stage 1, then `(u − 1/u)/2` (|x| ≥ 0.5), `(u + 1/u)/2`, `(u² − 1)/(u² + 1)` (|x| ≥ 0.5); |x| < 35 (tanh: 22); beyond 35, `e^|x|/2` is the answer to 2^-101. Smaller |x| go to stage 2. | `E = expm1(x)`: `sinh = (E + E/(1+E))/2`, `cosh = (u + 1/u)/2`, `tanh = E₂/(E₂ + 2)`, `E₂ = expm1(2x)`; no cancellation anywhere. |
| **pow** | `|x|^y = exp(y·log|x|)`: stage-1 log as a double-double (x not within 1/16 of 1), times `y` (double-double × double), stage-1 exp (needs |y log x| < 64); error bound `δ = 2^-61 + |y log x|·2^-62`. | Same with the stage-2 log and exp (|y log x| up to 746; δ = 2^-80). Beyond: ±∞/±0. `y = 1, 2, -1, ½` use `x`, `x·x`, `1/x`, `sqrt x`, which IEEE arithmetic already rounds correctly. Sign for negative bases: negative iff `y` is an odd integer. |
| **cbrt** | stage 2 only | `x = m·2^(3q)`, `m ∈ [1,8)`; initial `1 + (m−1)/3.5`; six Newton steps in double; one Newton step in double-double. |
| **hypot** | stage 2 only | Scale so the larger magnitude is in [1,2), `a² + b²` from two exact `TwoProd`s added in double-double, double-double `sqrt`. |

**Payne–Hanek reduction** (`Cores.Reduce`, Payne & Hanek 1983): `|x| = m·2^e`, `m` an integer below 2^53; the binary expansion of 2/π (24 words of 64 bits, generated, `Consts.TwoOverPi`) starting at word `j0 = max(0, (e−2)/64)` for five words is multiplied by `m` into a 6-limb integer; the 2 integer bits give the quadrant, the next 192 bits the fraction; a fraction ≥ ½ is negated (two's complement) and the quadrant incremented; the fraction is converted to a double-double (six 32-bit pieces added smallest first) and multiplied by π/2 (double-double). Truncation error is below 2^-192 of a quadrant, so `r` keeps more than 2^-100 relative accuracy even when `r ≈ 2^-62` (the closest a double gets to a multiple of π/2). The x = `6381956970095103·2^797` case (Muller's worst case for reduction) is in the known-answer file.

**Exact stage algorithms** (`Exact.cs`; value = `M·2^E`, `BigInteger` mantissa truncated to 512 bits after each multiplication, division and square root): π by Machin (`16 atan(1/5) − 4 atan(1/239)`) at 1900 bits, ln 2 = `2 atanh(1/3)`, ln 10 = log 10. `exp`: `k = round(x/ln2)`, reduce, `expm1` series on `r/2^24` and square 24 times. `log`: `x = 2^e m`, eight square roots, `atanh` series. `sin/cos`: `k` by exact integer division by π/2 (1900 bits), series on `|r| ≤ π/4`. `atan`: `1/x` for x > 1, six halvings `u ← u/(1+√(1+u²))`, series. `asin(x) = atan(x/√(1−x²))`, `acos(x) = 2 atan √((1−x)/(1+x))`, `sinh/cosh` from `e^x` (|x| ≥ 1) or `expm1`, `tanh = E₂/(E₂+2)`, `pow = exp(y log x)` with exact power-of-two bases (`2^(jy)` is exact when `jy` is an integer, which matters on subnormal ties), `exp2` of an integer is an exact power of two. `cbrt` and `hypot`: integer cube/square root of the exact integer with ≥ 160 bits and a sticky bit, so exact cases and ties round correctly. Final rounding (`ToDouble`, `ToSingle`) is integer shifting with round-half-even, producing subnormals and overflow correctly. The ≈ 2^-500 truncation error is far below the worst cases published for double (≈ 2^-120 for exp, log, trig; Lefèvre & Muller 2001, Lefèvre 2005 for pow); the only inputs where a value lands exactly on a rounding boundary are those with exact results, which the integer algorithms and the exact power-of-two rules handle.

## Correct rounding: how it is established

- **float32 (2^32 inputs, all 17 unary functions): exhaustive.** Every bit pattern is run; the platform's `Math.*` (double, accurate to ≲ 2^-50) screens each result, and every input whose float rounding the screen cannot settle (within 2^-40 of a rounding boundary — about 5·10^4 inputs per function) or where the result differs from the screened one is checked against the exact stage. Counts of non-correctly-rounded results: see the table. `atan2`, `pow`, `hypot` (two arguments): 4·10^7 random pairs each through the same procedure.
- **double: argument for correctness plus sampling.** Stage 2 decides a result only when its interval `v(1 ± 2^-85)` rounds uniquely, and the measured stage-2 error is ≤ 2^-96, so a double is wrong only if stage 2's error exceeded its declared bound by 2^11 — a claim the tests probe on 2·10^7 inputs per function (including ≈ 2·10^5 compared with the exact stage per function) and 1.6·10^4 per-function error measurements against the exact stage. Everything the tests could not settle goes to stage 3, which uses 500 bits. That is the standard argument of CRlibm/CORE-MATH-style two-step schemes; the unproved part is only that the polynomial/table error bound declared for stages 1 and 2 holds on every input, which is measured, not proved. Known hard cases cited: the argument-reduction worst case above; exact results (`exp2` integers, `log2` of powers of two, `cbrt`, `hypot`, `pow` with exact results) are in the known-answer file and special-case tests, including exact ties on subnormal boundaries (`exp2(-1075)`, `pow(2, -1075)`, `pow(4, -537)`). We did not reproduce the published worst-case tables for double (their inputs are not in the sources we may use); the exact stage is available as the arbiter on any input (`Exact.Eval`, internal; the tests call it).

## Hash matrix

Scenarios `functions/<fn>/f64|f32` (20 functions + `sincos`): 2000 pseudo-random inputs per function and precision plus 19 special values, generated by integer arithmetic only (SplitMix64, exponent/mantissa assembled from bits), every result's bit pattern hashed. The manifest of the fast + double-double + exact design is identical to the manifest of the double-double + exact design (the fast stage adds speed, never a different result).

## Test and tooling map

| What | Where |
|---|---|
| Known answers (3200, both precisions, mpmath at 700 bits) | `testdata/functions/kat.txt`, generator `tools/GenConstants/gen_kat.py`, test `FunctionTests.Known_answers_are_bit_exact` |
| Fast/double-double stages agree with the exact stage | `FunctionTests.Fast_stage_agrees_with_exact_stage` |
| Special values against IEEE | `SpecialValueTests` |
| Stage-1 error bound | `FastStageTests` |
| float32 exhaustive, binary float32 random, double random | `Sweeps` (off unless `ONEMATH_SWEEP=float|double`) |
| Speed against `Math`/`MathF` | `bench/Micro/FunctionsBench.cs` |

## Results (this build; every row is reproducible with the commands in the last section)

**float32, all 2^32 inputs, one run per function (nothing differed from the exact stage):**

| Function | inputs checked against the exact stage after screening | not correctly rounded |
|---|---|---|
| sin / cos / tan | 62,974 / 64,990 / 56,584 | 0 / 0 / 0 |
| asin / acos / atan | 5,332 / 11,939 / 21,954 | 0 / 0 / 0 |
| exp / exp2 / expm1 | 13,614 / 14,038 / 13,079 | 0 / 0 / 0 |
| log / log2 / log10 / log1p | 50,630 / 55,097 / 51,331 / 36,070 | 0 / 0 / 0 / 0 |
| cbrt / sinh / cosh / tanh | 103,318 / 7,916 / 6,324 / 6,720 | 0 / 0 / 0 / 0 |
| atan2 / pow / hypot, 4·10^7 random pairs each (hypot: 1.6·10^7 finite pairs) | 876 / 432 / 4 | 0 / 0 / 0 |

NaN, ±0, ±∞ and out-of-domain inputs (up to half of all bit patterns for the inverse trigonometric and logarithmic functions) are compared with the platform's result, including sign of zero (`expm1(-0)` keeps its sign here; the runtime's series screen does not, which is a screen artifact).

**double, 2·10^7 inputs per function** (generators in `Table.cs`: uniform and log-uniform over the function's range, near ±1 for asin/acos, integer and fractional exponents for pow): 2·10^5 per function compared with the exact stage: **0 differences for all 20 functions**. The rate at which the fast stage could not decide and the slow stage ran (per 5·10^4 probed inputs): sin 108, cos 122, tan 120, asin 115, acos 160, atan 124, atan2 161, exp 909, exp2 902, expm1 281, log 512, log2 205, log10 504, pow 1200, sinh 331, cosh 389, tanh 442 — 0.2% to 2.4%. Stage-2 error was measured against the exact stage on 1.6·10^4 inputs per function (before the fast stage existed, so the double-double stage was the one measured) and never exceeded 2^-96 against a declared 2^-85 (pow: 2^-80).

**Fast-stage error bound (stage 1 against stage 2, ≥ 4·10^6 inputs per function, `FastStageTests`)** — worst relative error seen and the bound handed to the rounding test:

| Function | worst seen | bound used |
|---|---|---|
| sin, cos, tan | 2^-67.5 | 2^-62 |
| asin, acos, atan, atan2 | 2^-66.4 | 2^-62 |
| exp, exp2 | 2^-64.7 | 2^-61 |
| expm1 | 2^-63.6 | 2^-60 |
| log, log2, log10 | 2^-63.1 | 2^-60 |
| sinh, cosh, tanh | 2^-63.8 | 2^-60 |
| pow | 2^-62.1 | 2^-61 + \|y log x\|·2^-62 (at least 3 bits below its bound on every sample) |

**Speed, indicative only** (one shared development machine under load, single thread, 10^5 inputs, best of 5; not a gate measurement): ours / `Math`: sin 2.8×, cos 4.3×, tan 4.9×, asin 9×, acos 11×, atan 12×, atan2 5.5×, exp 12×, exp2 7.7×, log 15×, log2 14×, log10 12×, pow 16×, cbrt 10×, hypot 4.8×, sinh 13×, cosh 21×, tanh 12×. Before the fast stage (double-double stage only) the ratios were 11× to 73×. Where the time goes is the double-double glue (a few dozen flops per call beyond the table-and-polynomial core) and the by-value `R` record.

**Hash matrix** (`functions/*`, 42 scenarios, manifest `12f01fe38319903cec79bea6cdb88a6a04447cec418d13b7aec71d6dd5b660d6`): equal on x64 default, `DOTNET_EnableHWIntrinsic=0`, `DOTNET_EnableAVX2=0` with 4 threads, `TieredCompilation=0`, quick-jit-for-loops off, NativeAOT x64, osx-arm64 and browser-wasm (Node). The same manifest was produced before and after the fast stage was added.

## Reproducing

```
# known answers (needs mpmath):   python3 tools/GenConstants/gen_kat.py > testdata/functions/kat.txt
# tables and constants:           python3 tools/GenConstants/gen_tables.py > src/OneMath/Functions/FastTables.cs
#                                 python3 tools/GenConstants/gen_constants.py > src/OneMath/Functions/Constants.cs
dotnet test tests/OneMath.Tests -c Release                       # unit tests, error bounds, known answers
ONEMATH_SWEEP=float  ONEMATH_FN=exp dotnet test --no-build -c Release tests/OneMath.Tests --filter-method '*Float32_exhaustive_unary'   # one function, ~2-6 min on 8 threads
ONEMATH_SWEEP=float  ONEMATH_N=40000000 dotnet test ... --filter-method '*Float32_random_binary'
ONEMATH_SWEEP=double ONEMATH_N=20000000 ONEMATH_EXACT_EVERY=100 dotnet test ... --filter-method '*Double_random'
dotnet run -c Release --project bench/Micro -- --filter '*FunctionsBench*'   # speed against Math/MathF, ~20-25 min (add --job short for ~8 min)
dotnet run -c Release --project tools/HashMatrix -- --filter functions/       # hash-matrix manifest
```


<!-- https://recursiveai.org/onemath/docs/predicates/ -->

# Exact geometric predicates

Normative description of `RecursiveAI.OneMath.Predicates`. A port that follows this text produces the same signs on every input, and the same expansion components in the internal engine.

Sources: Shewchuk, *Adaptive Precision Floating-Point Arithmetic and Fast Robust Geometric Predicates* (Discrete & Computational Geometry 18, 1997); Dekker (1971), Knuth (TAOCP vol. 2) and Priest (1991) for the error-free transformations.

## Arithmetic

Only IEEE 754 binary64 `+`, `-`, `*` with round-to-nearest-even and `abs`. No fused multiply-add anywhere; products are made exact by Dekker splitting. Float32 inputs convert to double exactly, so every predicate accepts them unchanged.

## Public API

| Call | Returns the sign of |
|---|---|
| `Orient2D(a, b, c)` | `(ax-cx)(by-cy) - (ay-cy)(bx-cx)`. +1 = counter-clockwise (c left of a→b, y up). |
| `Orient3D(a, b, c, d)` | `det[a-d; b-d; c-d]`. +1 = d below the plane abc, where a, b, c run counter-clockwise seen from above (Shewchuk's convention: a point above the plane gives -1). |
| `InCircle(a, b, c, d)` | `det[(a-d, \|a-d\|²); (b-d, \|b-d\|²); (c-d, \|c-d\|²)]`. +1 = d inside the circle when a, b, c are counter-clockwise. |
| `InSphere(a, b, c, d, e)` | `det[(a-e, \|a-e\|²); …; (d-e, \|d-e\|²)]` (4×4, rows a, b, c, d). +1 = e inside the sphere when `Orient3D(a, b, c, d) > 0`. |
| `SegmentPredicates.Classify(a, b, c, d)` | relation of closed segments ab and cd: `Disjoint`, `Proper`, `Touching`, `CollinearOverlap`. |

The result is `-1`, `0` or `+1`, and equals the sign of the determinant of the *exact real values* of the input doubles.

## Input requirements and range

- Every coordinate must be finite. NaN or ±∞ throws `ArgumentOutOfRangeException`.
- Exactness is guaranteed when each coordinate is 0 or has magnitude in [2^-140, 2^200]. The determinants have degree at most 5 (in-sphere), so every product then stays within [2^-960, 2^1005]: nothing overflows and the error term of no product underflows.
- **Out-of-range input may throw.** Outside that range a call that needs the exact path throws `ArgumentOutOfRangeException`; its `ParamName` is the coordinate parameter of the method that was called (`cx`, `dz`, ...). A call the filter decides by margin returns that sign, which is still correct, so whether an out-of-range call throws depends on the values; do not rely on either outcome and keep inputs inside the range. `SegmentPredicates.Classify` always checks the range first, so it throws for every out-of-range coordinate.
- A later release may extend the exact path (for example by power-of-two pre-scaling) so that out-of-range finite input stops throwing; removing an exception is compatible, so callers must not depend on it being thrown.
- The coordinate rule of the family (|x| < 2^21, resolution 2^-32) is far inside the range.

## Stage A: the floating-point filter

Each predicate first evaluates its determinant in double with a fixed operation order and also accumulates the *permanent*: the same expression with every subtraction replaced by addition of absolute values. With ε = 2^-53 the result is accepted when `|det| > c·permanent + 2^-1000`, with

| predicate | c |
|---|---|
| orient2d | (3 + 16ε)ε |
| orient3d | (7 + 56ε)ε |
| incircle | (10 + 96ε)ε |
| insphere | (16 + 224ε)ε |

These are Shewchuk's published stage-A bounds for exactly this evaluation order (the orders are written down in the code comments and must not be reassociated). `2^-1000` is an absolute term that has no effect inside the supported range. NaN and ∞ never pass the comparison, so they reach the exact path and are rejected there.

Evaluation orders (translated differences first):

- orient2d: `detleft = (ax-cx)*(by-cy)`, `detright = (ay-cy)*(bx-cx)`, `det = detleft - detright`; permanent `|detleft| + |detright|`.
- orient3d: `det = adz*(bdx*cdy - cdx*bdy) + bdz*(cdx*ady - adx*cdy) + cdz*(adx*bdy - bdx*ady)` with `adx = ax-dx` etc., left-to-right; permanent with absolute values.
- incircle: lifts `alift = adx*adx + ady*ady`, then `det = alift*(bdx*cdy - cdx*bdy) + blift*(cdx*ady - adx*cdy) + clift*(adx*bdy - bdx*ady)`; permanent likewise.
- insphere: 2×2 minors `xi*yj - xj*yi` of the translated (x, y) columns for the pairs ab, bc, cd, da, ac, bd; 3×3 minors `abc = aez*bc - bez*ac + cez*ab`, `bcd = bez*cd - cez*bd + dez*bc`, `cda = cez*da + dez*ac + aez*cd`, `dab = dez*ab + aez*bd + bez*da`; lifts `x²+y²+z²` added left to right; `det = (dlift*abc - clift*dab) + (blift*cda - alift*bcd)`; permanent mirrors it with absolute values.

## Stage D: exact evaluation

When the filter cannot decide, the determinant is evaluated exactly with expansion arithmetic and its sign read from the largest component. (The intermediate adaptive stages of the paper are not used: the exact stage is entered directly. Its sign is the same by construction; only the cost of near-degenerate inputs differs.)

- Translated coordinate differences are exact two-component expansions (`TwoDiff`), except orient2d which uses the translation-free form `ax(by-cy) + bx(cy-ay) + cx(ay-by)`, equal to the filter's determinant as a real number.
- 3×3 determinants expand along the first row of 2×2 minors; the in-sphere 4×4 uses the cofactor form of the filter's last line, with the four 3×3 minors computed on the translated rows in the order named above. Lifts are exact sums of exact squares.
- Products of expansions are sums, in increasing component order of the right factor, of `Scale` results.

## Expansion arithmetic (internal engine)

`ExpansionArithmetic` is internal: it is the engine behind `ExactPredicates` and `SegmentPredicates`, not public API (its preconditions are unchecked, and its output-size rules are not a promise). The text below is normative for the predicates' results.

An expansion is a list of doubles, sorted by increasing magnitude, nonoverlapping (the highest set bit of each component is below the lowest set bit of the next), with no zero component. The value is the exact sum; the sign is the sign of the last component; zero is the empty list.

- `TwoSum(a,b)`: `x = a+b; bv = x-a; av = x-bv; y = (a-av)+(b-bv)`. `TwoDiff` likewise with `bv = a-x; av = x+bv; y = (a-av)+(bv-b)`.
- `Split(a)`: `c = (2^27+1)·a; abig = c-a; hi = c-abig; lo = a-hi`. `TwoProduct(a,b)`: `x = a·b`, then with both operands split, `y = alo·blo - (((x - ahi·bhi) - alo·bhi) - ahi·blo)`. Requires |a|,|b| < 2^995 and the product's error term not to underflow.
- `Sum(e,f)`: copy e, then for each component of f run Grow-Expansion (Shewchuk Alg. 3): `Q = b; for each h_i: (Q, hh) = TwoSum(Q, h_i); keep hh if nonzero`; append Q if nonzero. `Difference` adds the negated components.
- `Scale(e,b)` (Shewchuk Alg. 6): `(Q, h) = TwoProduct(e_1, b)`; for i ≥ 2: `(T, t) = TwoProduct(e_i, b)`; `(S, h) = TwoSum(Q, t)`; `(Q, h') = FastTwoSum(T, S)`; keep nonzero h, h'; append Q.
- `Product(e,f)`: the sum over components f_j (increasing) of `Scale(e, f_j)`, accumulated with `Sum`.
- Output sizes: `Sum`/`Difference` ≤ |e|+|f|; `Scale` ≤ 2|e|; `Product` ≤ 2|e||f| (the caller's buffer needs that capacity).

Scratch space for intermediates is a per-thread bump allocator of 2^17 doubles allocated on first exact evaluation. Which thread ran a call never affects its result.

## Segment classification

With `o1 = Orient2D(a,b,c)`, `o2 = Orient2D(a,b,d)`, `o3 = Orient2D(c,d,a)`, `o4 = Orient2D(c,d,b)`:

1. A segment with equal endpoints is a point; point–point is `Touching` if equal else `Disjoint`; point–segment is `Touching` if collinear and inside the segment's bounding box, else `Disjoint`.
2. `Disjoint` if `o1·o2 > 0` or `o3·o4 > 0`.
3. If `o1 = o2 = 0` (collinear): order each segment's endpoints lexicographically (x, then y, the order along the line), take the larger lower end and smaller upper end: `CollinearOverlap` if lower < upper, `Touching` if equal, else `Disjoint`.
4. `Proper` if `o1·o2 < 0` and `o3·o4 < 0`; otherwise `Touching`.

## Checks

`tests/OneMath.Tests/PredicateTests.cs` compares every predicate with exact BigInteger evaluation (inputs scaled by 2^1074 so each double is an integer) on random, small-integer, wide-exponent, large-offset, and near-degenerate sets (collinear, coplanar, cocircular, cospherical points nudged by 0–2 ulps); the case count comes from `ONEMATH_PREDICATE_CASES` (default 100 000 per predicate). The input generator is `tools/HashMatrix/Scenarios/Predicates/PredicateInputs.cs`; the hash-matrix scenarios `predicates/*` hash the signs over fixed sets. Benchmarks: `bench/Micro/PredicateBench.cs`.


<!-- https://recursiveai.org/onemath/docs/evaluation-order/ -->

# Evaluation order of the vector, matrix and quaternion types

Normative: a port must evaluate in this order to produce the same bits.

## Rules

- Only IEEE 754 binary32/binary64 `+ - * /` and square root. Each operation is rounded to nearest, ties to even, in the type's own precision (binary32 for the `F` types, binary64 for the `D` types). No fused multiply-add, no reciprocal or square-root estimate instructions, no extended precision.
- Every `a*b + c*d` rounds each product, then rounds the sum. Sums of three or more terms accumulate left to right: `((t1 + t2) + t3) + t4`. Subtractions in the formulas below are written in the order given.
- Division is a true division (`x / s`), never a multiplication by a precomputed reciprocal.
- Square root is the correctly rounded square root. For binary32 it is computed as the binary64 square root of the widened value, rounded once to binary32; this equals the correctly rounded binary32 square root (double rounding is innocuous when the wide format has at least 2p + 2 significand bits; Figueroa 1995).
- Scaling by 0.5 or 2 is exact and written as a multiplication.
- Min and Max return the second operand when the comparison `a < b` (or `a > b`) is false, so a NaN in either operand gives the second operand.
- **NaN:** whether an operation that produces a NaN yields a positive or negative NaN, and which payload, is decided by the CPU (x64 and arm64 differ) and is not part of the contract. Everything else, including infinities, signed zeros and subnormals, is bit-exact. The hash-matrix runner writes every NaN as one canonical NaN.

## Vectors

- `Dot(a, b)` = `((a.X*b.X + a.Y*b.Y) + a.Z*b.Z) + a.W*b.W` (as many terms as the dimension).
- `LengthSquared` = the dot product with itself; `Length` = sqrt of that. `Normalized` divides each component by `Length` and returns zero when `Length == 0`.
- `Cross` (3D) = `(a.Y*b.Z - a.Z*b.Y, a.Z*b.X - a.X*b.Z, a.X*b.Y - a.Y*b.X)`; 2D cross = `a.X*b.Y - a.Y*b.X`.
- `Lerp(a, b, t)` = `a + (b - a) * t`, component-wise.

## Matrices

Row-major fields `Mrc`; column-vector convention (`M * v`). Translation sits in column 4 (`M14`, `M24`, `M34`) and `a * b` applies `b` first. This is the transpose of System.Numerics' row-vector convention, which keeps the same field names; `FromSystemNumerics` / `ToSystemNumerics` on `Matrix4x4` transpose. The product of System.Numerics matrices `T * R` (translate, then rotate) corresponds to `R * T` here.

- `CreateRotationX/Y/Z(radians)` with `(s, c) = SinCos(radians)`: X is rows `(1,0,0) (0,c,-s) (0,s,c)`, Y is `(c,0,s) (0,1,0) (-s,0,c)`, Z is `(c,-s,0) (s,c,0) (0,0,1)`. The 4x4 forms embed these with zero translation (`CreateFromMatrix3x3`: entries copied, `0 0 0 1` bottom row). `CreateFromQuaternion(q)` is `q.ToMatrix()`.

- Product entry (r, c) = `((a.Mr1*b.M1c + a.Mr2*b.M2c) + a.Mr3*b.M3c) + a.Mr4*b.M4c`; matrix-vector likewise.
- 2x2 determinant `M11*M22 - M12*M21`; inverse entries `M22/det, -M12/det, -M21/det, M11/det`.
- 3x3 cofactors `c11 = M22*M33 - M23*M32`, `c12 = M23*M31 - M21*M33`, `c13 = M21*M32 - M22*M31`, `c21 = M13*M32 - M12*M33`, `c22 = M11*M33 - M13*M31`, `c23 = M12*M31 - M11*M32`, `c31 = M12*M23 - M13*M22`, `c32 = M13*M21 - M11*M23`, `c33 = M11*M22 - M12*M21`; determinant `(M11*c11 + M12*c12) + M13*c13`; inverse = `[c11 c21 c31; c12 c22 c32; c13 c23 c33] / det` (each entry one division).
- 4x4: minors of rows 1-2 `s0 = M11*M22 - M21*M12`, `s1 = M11*M23 - M21*M13`, `s2 = M11*M24 - M21*M14`, `s3 = M12*M23 - M22*M13`, `s4 = M12*M24 - M22*M14`, `s5 = M13*M24 - M23*M14`; of rows 3-4 `c5 = M33*M44 - M43*M34`, `c4 = M32*M44 - M42*M34`, `c3 = M32*M43 - M42*M33`, `c2 = M31*M44 - M41*M34`, `c1 = M31*M43 - M41*M33`, `c0 = M31*M42 - M41*M32`. Determinant `s0*c5 - s1*c4 + s2*c3 + s3*c2 - s4*c1 + s5*c0` summed left to right. Inverse entries are the three-term sums in `Matrix4x4.TryInvert` (each term a product of one matrix entry and one minor, summed left to right, then divided by the determinant).
- `TryInvert` fails (returns false, zero matrix) when the determinant is zero or not finite. There is no conditioning test.

## Quaternions

- Hamilton product `a * b`: `X = ((a.W*b.X + a.X*b.W) + a.Y*b.Z) - a.Z*b.Y`, `Y = ((a.W*b.Y - a.X*b.Z) + a.Y*b.W) + a.Z*b.X`, `Z = ((a.W*b.Z + a.X*b.Y) - a.Y*b.X) + a.Z*b.W`, `W = ((a.W*b.W - a.X*b.X) - a.Y*b.Y) - a.Z*b.Z`.
- `Rotate(v)`: `t = 2 * cross(q.xyz, v)`, result `(v + t*W) + cross(q.xyz, t)`.
- `ToMatrix`: `xx = X*X` etc., diagonal `1 - 2*(yy + zz)`, `1 - 2*(xx + zz)`, `1 - 2*(xx + yy)`, off-diagonals `2*(xy - wz)`, `2*(xz + wy)`, `2*(xy + wz)`, `2*(yz - wx)`, `2*(xz - wy)`, `2*(yz + wx)`.
- `CreateFromRotationMatrix`: trace method; branch order trace > 0, then `M11 > M22 && M11 > M33`, then `M22 > M33`, else M33.
- `NLerp`: negate `b` when `Dot(a, b) < 0`; then component-wise `Lerp`; then `Normalized` (identity for a zero result).
- `CreateFromAxisSinCos(axis, sinHalf, cosHalf)` only multiplies: `(axis.X*sinHalf, axis.Y*sinHalf, axis.Z*sinHalf, cosHalf)`. `CreateFromAxisAngle(axis, radians)` is exactly `CreateFromAxisSinCos(axis, SinCos(radians / 2))` with the library's own `SinCos` (`Math64` for double, `Math32` for float).
- `Slerp(a, b, t)`: `d = Dot(a, b)`; if `d < 0` negate `b` and `d`; if `d > 0.9995` return `Lerp(a, b, t).Normalized()`; otherwise `theta = Acos(d)`, `s = Sin(theta)`, `wa = Sin((1 - t) * theta) / s`, `wb = Sin(t * theta) / s`, and each component is `a.c * wa + b.c * wb` (two products, one sum; not renormalized).

## Transforms and bounds

- `Affine3`: `TransformPoint(p) = Linear*p + Translation`; composition `a*b` has linear part `a.Linear*b.Linear` and translation `a.Linear*b.Translation + a.Translation`; the inverse is `(L^-1, -(L^-1 * t))`.
- `RigidTransform`: same shape with the quaternion in place of the matrix; `Inverse` uses the conjugate.
- Bounds: `Center = Min + (Max - Min) * 0.5`; `Volume = (sx*sy)*sz`; `SurfaceArea = 2 * ((sx*sy + sy*sz) + sz*sx)`; empty boxes have `Min = +inf, Max = -inf`.

## Naming and conversions

- Factories that build a value from parts or another concept are `Create...` (`CreateTranslation`, `CreateScale`, `CreateRotationX`, `CreateFromQuaternion`, `CreateFromRotationMatrix`, `CreateFromPoints`). `From...` is only for pure representation changes (`FromRows`, `FromColumns`, `FromSystemNumerics`). Instance methods that return a new value use participles (`Normalized`, `Transposed`, `Inverted`, `Conjugated`); `TryInvert(out)` is the fallible form. Row and column parameters of `FromRows` / `FromColumns` are 0-based (`row0`, `column0`, ...) like `Row(int)` and `Column(int)`.
- `FromSystemNumerics` / `ToSystemNumerics` exist for vectors, quaternions (same components, same Hamilton product) and 4x4 matrices (transposed). Float to float is a plain copy; double to float rounds to nearest even and is only available as the explicit method `ToSystemNumerics()`.
- JSON (`System.Text.Json`): every public value type serializes its public fields (or properties for the physics result types) by name and reads them back through its component constructor, bit for bit with default options and with source-generated metadata. Derived members (`QuaternionD.Vector`, `Vector4D.XYZ`, `Bounds.IsEmpty`/`Size`/`Extents`/`Center`, `PrincipalAxes.Orientation`) are not written. JSON has no NaN or infinity: such a value (for example `Bounds.Empty`) cannot be written with default options.


<!-- https://recursiveai.org/onemath/docs/fixed-bridge/ -->

# The Q31.32 fixed-point bridge

A deterministic game host keeps world coordinates in Q31.32: a signed 64-bit integer (the raw value) holding the value times 2^32. `FixedBridge` converts the raw `long` to and from floating point. It has no dependency on any fixed-point type.

## Range

| | Bound |
|---|---|
| Accepted value | finite, \|v\| < 2^21 = 2,097,152 |
| Accepted raw | \|raw\| <= 2^53 - 1 |

Anything else (including NaN and infinities) throws `ArgumentOutOfRangeException` whose message states the offending value and the bound (`TryToDouble` and `TryFromDouble` return false instead).

The limit is the exactness limit: 32 fractional bits plus 21 integer bits is 53 bits, the binary64 significand. Inside the range, every raw value is a binary64 number.

## Conversions

| Function | Result | Rounding |
|---|---|---|
| raw to double | `raw * 2^-32` | none: exact |
| raw to float32 | the exact value rounded once to binary32 | nearest, ties to even |
| double to raw | `round(v * 2^32)` | nearest integer, ties to even; `-0` gives 0 |
| float32 to raw | widen to double (exact), then as double to raw | as above |

- raw to float32 goes through the exact double, so there is exactly one rounding. Converting a `long` straight to `float` can round twice on some platforms.
- A double of magnitude 2^20 or more is always an exact multiple of 2^-32, so double to raw rounds only below that; the largest accepted double, 2^21 - 2^-32, maps to raw 2^53 - 1.
- Scaling by 2^32 and 2^-32 is exact (powers of two, no overflow or underflow inside the range), so raw to double and back is the identity on every raw value in range.
- The span forms convert in index order, throw on the first out-of-range element (earlier elements are already written) and require the destination to be at least as long as the source.

## Why float32 is not exact

A value like 1 + 2^-32 needs 33 significand bits, so no float32 holds it. Positions are therefore carried as double at the public boundary and as tile-local float32 (relative to a double origin) inside; a tile of extent 2^t resolves 2^(t-24).


<!-- https://recursiveai.org/onemath/docs/physics/ -->

# Physics math: symmetric 3x3 eigen decomposition and inertia tensors

Normative: a port must evaluate in this order to produce the same bits. Both precisions (`F` = binary32, `D` = binary64) use the same code with the type swapped. Everything follows the rules of `evaluation-order.md`: only `+ - * /` and the correctly rounded square root, each rounded in the type's precision, no fused operations, products and sums in the written order, sums left to right.

Types: `SymmetricEigenF/D`, `EigenDecompositionF/D`, `InertiaF/D`, `MassPropertiesF/D`, `PrincipalAxesF/D`, namespace `RecursiveAI.OneMath.Physics`.

## Input checks

- Any NaN or infinite entry, mass, length or offset throws `ArgumentOutOfRangeException` naming the parameter and the field (`matrix.M13 is NaN; a finite value is expected.`). A negative mass or length throws the same type.
- A tensor must be symmetric: all nine entries finite and each mirrored pair within `64 * eps * max|entry|` (`eps` = 2^-52 for D, 2^-23 for F), else `ArgumentException` (`not symmetric`). Functions read the **upper triangle** only (M11 M12 M13 M22 M23 M33); every tensor they return is exactly symmetric.
- `Combine` of an empty span, or of parts whose total mass is zero or not finite, throws `ArgumentException`. Errors for a part name it (`parts[1].Mass`).
- Rotation matrices and quaternions are checked for finiteness only; they are assumed to be rotations and unit quaternions.

## Symmetric eigen decomposition

Cyclic Jacobi (Golub and Van Loan, *Matrix Computations*, section 8.5) with Rutishauser's stable rotation.

1. **Scaling.** `m` = largest absolute upper-triangle entry. If `m == 0` the result is eigenvalues 0 and the identity frame. Otherwise `e` = the binary exponent of `m` (2^e <= m < 2^(e+1), subnormals included), `e1 = e >> 1` (arithmetic shift), `e2 = e - e1`. Every entry becomes `(x * 2^-e1) * 2^-e2` (exact except for entries that underflow, which are negligible against `m`). The scaled matrix has largest entry in [1, 2). The eigenvalues are scaled back by `(x * 2^e1) * 2^e2`. A consequence: scaling the input by a power of two scales the eigenvalues by it and leaves the eigenvectors bit-identical (while nothing leaves the normal range).
2. **Sweeps.** At most 16 sweeps; a sweep rotates the pairs (p,q) = (0,1), (0,2), (1,2) in that order; iteration stops after the first sweep in which no pair is rotated. The eigenvector matrix starts as the identity and is updated in the same rotations.
3. **One pair.** With `apq` the pivot, `g = 2*|apq|`: if `apq == 0` or both `|app| + g == |app|` and `|aqq| + g == |aqq|`, set `apq = aqp = 0` and do not rotate. Otherwise `d = aqq - app`. If `|d| > B * g` (B = 2^60 for D, 2^30 for F) then `t = apq / d`; else `theta = d / (apq*2)`, `t = 1 / (|theta| + sqrt(theta*theta + 1))`, negated when `theta < 0`. Then `c = 1 / sqrt(t*t + 1)`, `s = t*c`, `tau = s / (1 + c)`, `shift = t*apq`, `app = app - shift`, `aqq = aqq + shift`, `apq = aqp = 0`. For the third index `r`: `arp' = arp - s*(arq + tau*arp)`, `arq' = arq + s*(arp - tau*arq)`, mirrored to the symmetric positions. Each row `k` of V: `vkp' = vkp - s*(vkq + tau*vkp)`, `vkq' = vkq + s*(vkp - tau*vkq)`.
4. **Output order.** Eigenvalues ascending, by a stable insertion sort on the indices (0,1,2) with comparison `lambda[order[j]] > lambda[key]`; equal values keep the positions the iteration left them (a diagonal input keeps axis order).
5. **Signs.** For the first two eigenvectors the component of largest magnitude is made positive (the lowest index wins a tie; a flip is `0 - x`, so zeros stay +0). The third eigenvector is flipped if `dot(cross(c0, c1), c2) < 0`, so the frame is right-handed (determinant +1). The third vector therefore need not have a positive largest component.

### Accuracy

Stated bound, with `u` = 2^-53 (D) or 2^-24 (F) and `A` the input with `V` its eigenvector matrix and `L` the eigenvalues:

- `max |V L V^T - A| <= 16 u max|A|`
- `max |V^T V - I| <= 10 u` and `|det V - 1| <= 10 u`

It holds when the eigenvalues are in the normal range; a result that lands in the subnormal range is quantized by the format itself (a few subnormal steps are added to the first bound). By Weyl's inequality the first bound also bounds each eigenvalue's absolute error by about `3 * 16 u max|A|`. Eigenvalues that overflow the type become infinite (only possible when the largest entry is within a factor of about 3 of the type's maximum).

Measured over 10^6 random matrices per precision (random, wide dynamic range, rank one and two, repeated and nearly repeated eigenvalues, near-singular spreads, diagonal with repeats, small integers, binary scales across the whole exponent range) against a double-double oracle: D worst 12.4 u reconstruction and 7.5 u orthonormality; F worst 10.4 u and 6.2 u.

## Inertia

Tensors are row-major `Matrix3x3`; offsets, centres and tensors are in one frame.

- **Rotate** `R I R^T`: `S` = the symmetric matrix from the upper triangle, `M = R * S` (the matrix product of `evaluation-order.md`), entry `(i, j)` for `i <= j` = `Dot(row i of M, row j of R)`; the lower triangle mirrors the upper. A quaternion first becomes `q.ToMatrix()` (unit quaternion assumed).
- **ShiftToPoint** (centre of mass to a point at `d` = point minus centre; `m` the mass): `xx = dy*dy + dz*dz`, `yy = dx*dx + dz*dz`, `zz = dx*dx + dy*dy`, `mx = m*dx`, `my = m*dy`; `I11 + m*xx`, `I22 + m*yy`, `I33 + m*zz`, `I12 - mx*dy`, `I13 - mx*dz`, `I23 - my*dz`. **ShiftToCenter** subtracts the same terms on the diagonal and adds them off it. A round trip is exact only up to rounding.
- **Combine** parts (mass `m_i`, centre `c_i`, tensor about `c_i`): in span order, `M = ((0 + m_0) + m_1) + ...`, `S = Σ m_i * c_i` per component the same way, `C = S / M` per component; each tensor entry of the result is the left-to-right sum, from 0, of `ShiftToPoint(I_i, m_i, c_i - C)`. The result is the tensor about `C` with mass `M`.
- **PrincipalAxes**: the eigen decomposition above; `Moments` = eigenvalues, `Rotation` = eigenvector matrix (`I = Rotation * diag(Moments) * Rotation^T`), `Orientation` = `Quaternion.CreateFromRotationMatrix(Rotation)`.
- **Solids** (axes along the shape; `h` half extents): box diagonal `(m*(hy*hy + hz*hz))/3`, `(m*(hx*hx + hz*hz))/3`, `(m*(hx*hx + hy*hy))/3`; sphere `((2*m)*(r*r))/5`; cylinder along Y, `rr = r*r`, `hq = hh*hh` for half height `hh`: Ixx = Izz = `(m*(3*rr + 4*hq))/12`, Iyy = `(m*rr)/2`.

## Tests and hash-matrix scenarios

`tests/OneMath.Tests/Physics/` holds known-answer tests (diagonal, 2- and 3-fold repeats, rank-deficient, tridiagonal, extreme scales, power-of-two commutation, hostile input), the 10^6-case oracle sweeps per precision (`ONEMATH_PHYSICS_CASES` overrides the count) and the inertia identities. The scenarios `physics/{f32,f64}/{eigen,inertia}` are in the hash matrix.


<!-- https://recursiveai.org/onemath/docs/hash-matrix/ -->

# Hash matrix

A scenario is a deterministic function that writes bytes to a `ScenarioContext`; its SHA-256 goes in a manifest. Two platforms agree when their manifests are identical.

## Adding a scenario

Create a file under `tools/HashMatrix/Scenarios/<area>/` (the `HashMatrix` project and the smoke apps compile everything beneath it) with a module initializer:

```csharp
namespace RecursiveAI.OneMath.HashMatrix;

internal static class MyScenarios
{
    [System.Runtime.CompilerServices.ModuleInitializer]
    internal static void Init() => Scenarios.Register("area/name", ctx => { ctx.Write(1.0); });
}
```

No shared file changes. Names are lower-case, slash separated and unique; the manifest orders them by ordinal comparison. Writers (`Write(int|long|float|double|bytes)`) are little-endian by bit pattern; every NaN is written as one canonical NaN because arithmetic NaN sign and payload are chosen by the CPU. `ctx.Threads` is the worker count a parallel scenario should use; the output must not depend on it (split by index, combine in index order). Use a seeded generator, never `Random` or time.

## Running

`HashMatrix [--threads N] [--filter PREFIX] [--out FILE] [--list]` prints `name sha256` per scenario and a final `manifest sha256` over those lines. `tools/hash-matrix.sh` runs every column (see its header) and compares them.

Columns: x64 JIT; x64 JIT with `DOTNET_EnableHWIntrinsic=0`; with `DOTNET_EnableAVX2=0`; 1 thread and all cores; x64 NativeAOT (`tests/AotSmoke`); browser-wasm interpreted and AOT-compiled (`tests/WasmSmoke`, run under Node); osx-arm64.


<!-- https://recursiveai.org/onemath/docs/changelog/ -->

# Changelog

All notable changes to this project are documented here. The format follows [Keep a Changelog](https://keepachangelog.com/en/1.1.0/), and the project uses [Semantic Versioning](https://semver.org/) with pre-release suffixes for previews.

## [1.0.0]

The first stable release. The public API is frozen under semantic versioning. Results are bit-identical to 0.1.0-preview.5: there is no numeric change.

From this release on, results are also bit-identical across every release of a major version: no 1.x minor or patch release changes an existing result (see "The determinism promise" in the README).

### Changed (breaking since 0.1.0-preview.5)

One verb per concept: `Create...` builds a value from parts or another concept, `From...` is only a representation change, and instance methods that return a new value are participles. Numeric results are unchanged; only names move.

| Before | Now |
|---|---|
| `Matrix2x2/3x3/4x4 D/F .Transpose()` | `.Transposed()` |
| `QuaternionD/F .Conjugate()` | `.Conjugated()` |
| `RigidTransformD/F .Inverse()` | `.Inverted()` |
| `QuaternionD/F .FromAxisSinCos(axis, sin, cos)` | `.CreateFromAxisSinCos(axis, sin, cos)` |
| `QuaternionD/F .FromMatrix(m)` | `.CreateFromRotationMatrix(m)` |
| `Affine3D/F .FromTranslation(v)` | `.CreateTranslation(v)` |
| `Bounds2/3 D/F .FromPoints(points)` | `.CreateFromPoints(points)` |
| `Bounds2/3 D/F .FromCenterExtents(c, e)` | `.CreateFromCenterExtents(c, e)` |
| `FixedBridge.MaxMagnitude` (`const`) | `FixedBridge.ExclusiveMagnitudeBound` (static property; the bound is exclusive) |
| `PrincipalAxesD/F(moments, rotation, orientation)` | `PrincipalAxesD/F(moments, rotation)`; `Orientation` is derived from `Rotation`, so the two cannot disagree |
| `MassPropertiesD/F(mass, ...)` accepted a negative mass | throws `ArgumentOutOfRangeException` (parameter `mass`) |

Also: the parameters of `FromRows` / `FromColumns` are now `row0..` / `column0..` (0-based, like `Row(int)` and `Column(int)`); `ParamName` of every argument exception names a parameter of the public method that was called (predicates: `ax`..`ez`; `InertiaD/F.PrincipalAxes`: `tensor`); all public value types serialize with System.Text.Json (public fields and component constructors, derived members ignored; bit-exact round trip with default options and with source-generated metadata); the package project URL points at the repository; the matrix convention (column vectors, the opposite of System.Numerics) and the predicate input range (out-of-range input may throw) are documented in the README, the package readme and the type and method docs.

### Added

- `QuaternionD/F`: `CreateFromAxisAngle` (exactly `CreateFromAxisSinCos(axis, SinCos(radians / 2))`) and `Slerp` (shortest arc; evaluation order in `docs/evaluation-order.md`).
- `Matrix3x3D/F`: `CreateRotationX/Y/Z`, `CreateFromQuaternion`. `Matrix4x4D/F`: `CreateRotationX/Y/Z`, `CreateFromQuaternion`, `CreateFromMatrix3x3`, `CreateFromAffine`.
- `Affine3D/F`: `CreateScale`, `CreateFromLinear`, `CreateFromQuaternion`, `CreateRotationX/Y/Z`.
- `FromSystemNumerics` / `ToSystemNumerics` on `Vector2/3/4`, `Quaternion` and `Matrix4x4` (D and F). The matrix conversions transpose; double to float rounds to nearest even.
- Hash-matrix scenarios `types/f32/rotation`, `types/f64/rotation` and `types/json` (a source-generated JSON round trip that also runs in the NativeAOT and WebAssembly smoke apps). Every earlier manifest line is unchanged.

### Removed

- `ExpansionArithmetic` is no longer public. Its preconditions are not checked (invalid input gives silently wrong results), and publishing it would freeze Shewchuk's output-size rules. It remains the internal engine behind `ExactPredicates` and `SegmentPredicates`; exposing a checked form later would be an additive change.

## [0.1.0-preview.5]

The first version published to nuget.org. It is a preview: public API may change before 1.0.

### Added

- Vectors (2, 3, 4), matrices (2x2, 3x3, 4x4), quaternions, affine and rigid transforms and axis-aligned bounds, in float and double, with a documented evaluation order (`docs/evaluation-order.md`).
- Correctly rounded elementary functions for float and double (`Math32`, `Math64`): sin, cos, tan, asin, acos, atan, atan2, exp, exp2, expm1, log, log2, log10, log1p, pow, cbrt, hypot, sinh, cosh, tanh, sincos. Float results were checked exhaustively; double results against an exact big-number stage. They are slower than the platform's `Math` functions (about 3x to 20x single-threaded in an informal measurement; see `docs/functions.md`), in exchange for identical bits everywhere.
- Exact geometric predicates (`Orient2D`, `Orient3D`, `InCircle`, `InSphere`, segment classification) with a fast filter that always agrees with the exact sign.
- The Q31.32 fixed-point bridge (`FixedBridge`): exact conversion for magnitudes below 2^21, rejecting beyond.
- Symmetric 3x3 eigen decomposition and inertia-tensor math (`RecursiveAI.OneMath.Physics`).
- NativeAOT and trimming compatible; no dependencies beyond the .NET base class library.
- A hash-matrix runner that compares results across JIT, scalar, NativeAOT and WebAssembly builds, and across x64 and arm64.

### Known limits

- Preview status: no API stability promise yet.
- The correctly rounded functions favour exactness over speed; faster opt-in variants are not part of this release.

[1.0.0]: https://github.com/recursive-ai-oss/onemath-csharp/releases/tag/v1.0.0
[0.1.0-preview.5]: https://github.com/recursive-ai-oss/onemath-csharp/releases/tag/v0.1.0-preview.5


<!-- https://recursiveai.org/onemath/cookbook/exact-predicates/ -->

# Exact geometric predicates

Decide orientation, in-circle and segment intersection with OneMath’s exact predicates, and see how often the naive floating-point determinant gets the sign wrong.

The classic test, after Shewchuk: a line through (12, 12) and (24, 24), and a grid of 256 × 256 points around (0.5, 0.5), one unit in the last place apart. Each of them lies on, above or below the line; the naive formula is asked which.

```csharp
using RecursiveAI.OneMath.Predicates;

// Shewchuk's example: is p left of, right of, or on the line through a and b?
// p sweeps a 256 x 256 grid of the doubles nearest (0.5, 0.5), one ulp apart.
double ax = 12, ay = 12, bx = 24, by = 24;
double ulp = Math.Pow(2, -53);
int wrong = 0;
for (int i = 0; i < 256; i++)
for (int j = 0; j < 256; j++)
{
    double px = 0.5 + i * ulp, py = 0.5 + j * ulp;
    double naive = (ax - px) * (by - py) - (ay - py) * (bx - px);
    if (Math.Sign(naive) != ExactPredicates.Orient2D(ax, ay, bx, by, px, py)) wrong++;
}
Console.WriteLine($"naive determinant has the wrong sign for {wrong} of {256 * 256} points");

double qx = 0.5 + 3 * ulp, qy = 0.5 + 1 * ulp;
Console.WriteLine($"one of them: naive {Math.Sign((ax - qx) * (by - qy) - (ay - qy) * (bx - qx))}, exact {ExactPredicates.Orient2D(ax, ay, bx, by, qx, qy)}");

// In-circle test for Delaunay triangulations: +1 inside, -1 outside, 0 on the circle.
Console.WriteLine($"InCircle (1,1): {ExactPredicates.InCircle(0, 0, 1, 0, 0, 1, 1, 1)}");
Console.WriteLine($"InCircle (0.5,0.5): {ExactPredicates.InCircle(0, 0, 1, 0, 0, 1, 0.5, 0.5)}");

// Segment classification.
Console.WriteLine(SegmentPredicates.Classify(0, 0, 2, 2, 0, 2, 2, 0));   // cross in the middle
Console.WriteLine(SegmentPredicates.Classify(0, 0, 2, 2, 2, 2, 3, 0));   // share an endpoint
Console.WriteLine(SegmentPredicates.Classify(0, 0, 2, 2, 1, 1, 3, 3));   // overlap on a line
Console.WriteLine(SegmentPredicates.Classify(0, 0, 1, 0, 0, 1, 1, 1));   // parallel, apart
```

```text
naive determinant has the wrong sign for 11972 of 65536 points
one of them: naive 0, exact -1
InCircle (1,1): 0
InCircle (0.5,0.5): 1
Proper
Touching
CollinearOverlap
Disjoint
```

## What they return

- `Orient2D(ax, ay, bx, by, cx, cy)`: +1 if c is to the left of the directed line a→b (counterclockwise), −1 to the right, 0 on it. `Orient3D` does the same for a point and a plane.
- `InCircle(a, b, c, d)`: +1 if d is inside the circle through a, b and c (taken counterclockwise), −1 outside, 0 on it. `InSphere` is the 3D form. Here (1, 1) is exactly on the circle through (0, 0), (1, 0) and (0, 1).
- `SegmentPredicates.Classify` returns a `SegmentRelation`: `Disjoint`, `Proper` (they cross at one interior point), `Touching` (they meet at an endpoint) or `CollinearOverlap`.

The arithmetic is adaptive: a fast floating-point filter answers when it can prove the sign, and exact expansion arithmetic finishes the rest, so the common case stays cheap.

## Good to know

- Coordinates must be finite and either zero or of magnitude in [2^-140, 2^200]. Out-of-range input may throw `ArgumentOutOfRangeException`.
- Exact signs are what keep triangulations, convex hulls and polygon clipping consistent: two calls about the same points can never contradict each other.

The algorithm and its error bounds are in [exact predicates](https://recursiveai.org/onemath/docs/predicates/).


<!-- https://recursiveai.org/onemath/cookbook/port-from-system-numerics/ -->

# Port from System.Numerics

Move matrix code from System.Numerics to OneMath without silently transposing your transforms: the conventions are opposite, and the conversions handle it.

```csharp
using RecursiveAI.OneMath;

float angle = MathF.PI / 2;

// System.Numerics (row vectors): T * R translates first, then rotates.
var sn = System.Numerics.Matrix4x4.CreateTranslation(10, 0, 0) * System.Numerics.Matrix4x4.CreateRotationY(angle);
// OneMath (column vectors): the same transform is R * T.
var m = Matrix4x4F.CreateRotationY(angle) * Matrix4x4F.CreateTranslation(new Vector3F(10, 0, 0));

var p = new Vector3F(1, 0, 0);
Console.WriteLine($"System.Numerics: {System.Numerics.Vector3.Transform(p.ToSystemNumerics(), sn)}");
Console.WriteLine($"OneMath:         {m.TransformPoint(p)}");

// Where the translation lives.
Console.WriteLine($"System.Numerics translation M41..M43: {sn.M41}, {sn.M42}, {sn.M43}");
Console.WriteLine($"OneMath translation M14, M24, M34:    {m.M14}, {m.M24}, {m.M34}");

// Convert, don't copy fields: both directions transpose.
Matrix4x4F same = Matrix4x4F.FromSystemNumerics(sn);
Console.WriteLine($"FromSystemNumerics(sn) == m: {same == m}");
Console.WriteLine($"Round trip equal: {Matrix4x4F.FromSystemNumerics(m.ToSystemNumerics()) == m}");

// Copying fields one to one gives the transpose: a different transform.
var copied = new Matrix4x4F(sn.M11, sn.M12, sn.M13, sn.M14, sn.M21, sn.M22, sn.M23, sn.M24,
                            sn.M31, sn.M32, sn.M33, sn.M34, sn.M41, sn.M42, sn.M43, sn.M44);
Console.WriteLine($"Field copy:      {copied.TransformPoint(p)}");
```

```text
System.Numerics: <-4.808253E-07, 0, -11>
OneMath:         (-4.808253E-07, 0, -11)
System.Numerics translation M41..M43: -4.371139E-07, 0, -10
OneMath translation M14, M24, M34:    -4.371139E-07, 0, -10
FromSystemNumerics(sn) == m: True
Round trip equal: True
Field copy:      (-4.371139E-08, 0, 1)
```

## The rule

- **OneMath uses column vectors**: `m * v` transforms the column `v`, a 4x4’s translation sits in column 4 (`M14`, `M24`, `M34`), and `a * b` applies `b` first.
- **System.Numerics uses row vectors**: translation in `M41`…`M43`, and `a * b` applies `a` first. The field names `Mrc` are the same, which is what makes a field-by-field copy look right and transform wrongly.
- So reverse the operand order of every product, and convert with `FromSystemNumerics` and `ToSystemNumerics`, which transpose. Both directions round-trip exactly.

Vectors and quaternions convert with the same pair of methods; quaternions keep their components and their Hamilton product. Converting from double to float is only the explicit `ToSystemNumerics()` and rounds to nearest even.

## Good to know

Both libraries give `-11` here. The difference is the promise: OneMath’s results are the same bits on every OS, CPU and runtime, and System.Numerics makes no such promise. Keep System.Numerics for rendering if you like, and use OneMath for whatever the simulation, a replay or a lockstep peer has to agree on.


<!-- https://recursiveai.org/onemath/cookbook/fixed-point-bridge/ -->

# The fixed-point bridge

Convert between Q31.32 fixed-point world coordinates and `double` with OneMath’s `FixedBridge`, exactly, with one documented rounding where one is needed and a clear error outside the range.

A deterministic game host often keeps positions as Q31.32: a 64-bit integer holding the value times 2^32. `FixedBridge` converts the raw `long` to and from floating point without depending on any fixed-point type.

```csharp
using RecursiveAI.OneMath;

// A deterministic host stores positions as Q31.32 raw longs (value * 2^32).
long x = FixedBridge.FromDouble(1234.5678);
Console.WriteLine($"raw {x}, back to double {FixedBridge.ToDouble(x):R}, as float {FixedBridge.ToSingle(x):R}");

// raw -> double -> raw is the identity for every raw value in range.
long[] raws = [1, -1, (1L << 53) - 1, 5_000_000_000];
double[] values = new double[raws.Length];
long[] again = new long[raws.Length];
FixedBridge.ToDouble(raws, values);
FixedBridge.FromDouble(values, again);
Console.WriteLine($"identity: {raws.SequenceEqual(again)}");

// Rounding: nearest, ties to even.
Console.WriteLine($"2^-33 -> {FixedBridge.FromDouble(Math.Pow(2, -33))}, 3*2^-33 -> {FixedBridge.FromDouble(3 * Math.Pow(2, -33))}");

// The range is |v| < 2^21.
Console.WriteLine($"bound {FixedBridge.ExclusiveMagnitudeBound}, in range 2e6: {FixedBridge.IsInRange(2e6)}, 3e6: {FixedBridge.IsInRange(3e6)}");
Console.WriteLine($"TryFromDouble(3e6): {FixedBridge.TryFromDouble(3e6, out _)}");
try { FixedBridge.FromDouble(double.NaN); }
catch (ArgumentOutOfRangeException e) { Console.WriteLine(e.Message); }
```

```text
raw 5302428325695, back to double 1234.567800000077, as float 1234.5677
identity: True
2^-33 -> 0, 3*2^-33 -> 2
bound 2097152, in range 2e6: True, 3e6: False
TryFromDouble(3e6): False
The value NaN is outside the exact range: it must be finite with magnitude below 2^21 = 2097152. (Parameter 'value')
Actual value was NaN.
```

## The rules

| Conversion | Result |
|---|---|
| raw to `double` (`ToDouble`) | `raw × 2^-32`, exact |
| raw to `float` (`ToSingle`) | the exact value rounded once to float, nearest, ties to even |
| `double` to raw (`FromDouble`) | `round(v × 2^32)`, nearest integer, ties to even; −0 gives 0 |
| `float` to raw (`FromSingle`) | widened to `double` (exact), then as above |

- The range is a finite value with |v| < 2^21 = 2,097,152, or a raw value with |raw| ≤ 2^53 − 1: 21 integer bits and 32 fractional bits are the 53 bits of a double’s significand, so every raw value in range is a double, and raw → double → raw is the identity.
- Outside it, `FromDouble` and `ToDouble` throw `ArgumentOutOfRangeException` naming the value and the bound; `TryFromDouble` and `TryToDouble` return false instead.
- The span forms convert in index order and throw on the first out-of-range element.

`2^-33` is exactly half a raw unit, so it rounds to the even neighbor, 0; `3 × 2^-33` is one and a half, and rounds to 2.

More in [the fixed-point bridge](https://recursiveai.org/onemath/docs/fixed-bridge/), including why float32 can’t hold every Q31.32 value.


<!-- https://recursiveai.org/onemath/cookbook/inertia-tensors/ -->

# Inertia tensors and principal axes

Build the mass properties of a body from parts with OneMath: solid shapes, a combined center of mass and inertia tensor, and its principal axes, with the same bits on every platform.

A body made of a 2 kg box and a 1 kg ball, offset from each other:

```csharp
using RecursiveAI.OneMath;
using RecursiveAI.OneMath.Physics;

// A body made of two parts: a 2 kg box and a 1 kg sphere offset along x.
var box = new MassPropertiesD(2, new Vector3D(0, 0, 0), InertiaD.SolidBox(2, new Vector3D(1, 0.5, 0.25)));
var ball = new MassPropertiesD(1, new Vector3D(1.5, 0.5, 0), InertiaD.SolidSphere(1, 0.4));

MassPropertiesD body = InertiaD.Combine([box, ball]);
Console.WriteLine($"mass {body.Mass}, center {body.Center}");
Console.WriteLine($"inertia about the center:\n{body.Inertia}");

PrincipalAxesD axes = InertiaD.PrincipalAxes(body.Inertia);
Console.WriteLine($"principal moments {axes.Moments}");
Console.WriteLine($"orientation {axes.Orientation}");
```

```text
mass 3, center (0.5, 0.16666666666666666, 0)
inertia about the center:
[0.43900000000000006, -0.5, 0; -0.5, 2.2723333333333335, 0; 0, 0, 2.564]
principal moments (0.3115029928215277, 2.399830340511806, 2.564)
orientation (0, 0, 0.12451268641694911, 0.9922180158217418)
```

## The pieces

- `InertiaD.SolidBox(mass, halfExtents)`, `SolidSphere(mass, radius)` and `SolidCylinder` (axis along Y, from its mass, radius and half height) give the tensor of a solid about its own center, axes along the shape. `InertiaF` is the float form.
- `MassPropertiesD(mass, center, inertia)` is one part; `InertiaD.Combine(parts)` returns the total mass, the center of mass and the tensor about that center (parallel-axis shifts included). `ShiftToPoint`, `ShiftToCenter` and `Rotate` are there for doing it by hand.
- `InertiaD.PrincipalAxes(tensor)` diagonalizes it: `Moments` are the principal moments in ascending order, `Rotation` the axes as a matrix and `Orientation` the same as a quaternion, with a right-handed frame and fixed signs.

The decomposition is cyclic Jacobi with a written-down order of operations, so the principal frame comes out the same, bit for bit, on every platform. NaN, infinity, a negative mass or a non-symmetric tensor throw an exception naming the field. Details in [inertia and eigen decomposition](https://recursiveai.org/onemath/docs/physics/).


<!-- https://recursiveai.org/onemath/cookbook/check-machines-agree/ -->

# Check that two machines agree

Hash a fixed OneMath computation and compare the hash across machines, runtimes and CPU features: the cheapest way to confirm a deterministic simulation will stay in sync.

```csharp
using System.Security.Cryptography;
using RecursiveAI.OneMath;

// A fixed computation whose hash two machines (or a replay and the original) can compare.
var buffer = new byte[8 * 100_000];
for (int i = 0; i < 100_000; i++)
{
    double t = i * 0.001;
    double v = Math64.Sin(t) * Math64.Exp(-t / 50) + Math64.Atan2(t, 3) * Math64.Log1p(t);
    BitConverter.TryWriteBytes(buffer.AsSpan(8 * i), v);
}
Console.WriteLine($"sha256 {Convert.ToHexString(SHA256.HashData(buffer))[..32]}");
```

Built once in Release and run with the CPU’s vector instructions turned off step by step, and with the JIT’s tiering off:

```text
$ dotnet run -c Release
sha256 D72748EDFF396FE462CF2D674C8D342C
$ DOTNET_EnableAVX512F=0 dotnet run -c Release
sha256 D72748EDFF396FE462CF2D674C8D342C
$ DOTNET_EnableAVX2=0 dotnet run -c Release
sha256 D72748EDFF396FE462CF2D674C8D342C
$ DOTNET_EnableHWIntrinsic=0 dotnet run -c Release
sha256 D72748EDFF396FE462CF2D674C8D342C
$ DOTNET_TieredCompilation=0 dotnet run -c Release
sha256 D72748EDFF396FE462CF2D674C8D342C
```

## Using it

- Run the same program on each target you ship (x64 and arm64, NativeAOT, WebAssembly) and compare one line. If a hash ever differs, the cause is in your own code: a call to `Math` or `MathF`, a `float` sum in a different order, or a dependence on thread scheduling.
- For a game, hash the simulation state every N frames and exchange it with peers; a mismatch then points at the first frame that diverged.

OneMath runs its own version of this, the hash matrix: fixed scenarios whose SHA-256 manifest must be equal on x64 JIT, scalar, AVX2 off, one thread and all cores, NativeAOT, browser WebAssembly (interpreted and AOT-compiled) and macOS arm64. Each 1.x release must also reproduce every line of the previous release’s manifest. See [hash matrix](https://recursiveai.org/onemath/docs/hash-matrix/).