Numerics contract

Numerics contract#

When aether computes a sum, a dot product or a division, which exact floating-point operations does it perform, and in what order?

This page states them. Floating-point addition is not associative, so a different order is a different (merely close) answer. Every target that computes against aether, whether the C++ templates or a generated language core, follows the rules below.

Folds#

A fold is any reduction over the components of an item: sum, dot, squaredNorm, norm, maxNorm. It is sequential and right-nested, with the highest index outermost and no zero seed:

e[N-1] + (e[N-2] + (... + e[0]))

so a 3-vector dot is (a2*b2) + ((a1*b1) + (a0*b0)). The static forms (aether/expr/Reduce.h) and the rolled-loop forms in aether::rt (aether/expr/Rt.h) use this same order. A matrix-vector row, a matrix product entry and the Frobenius matDot (row-major) are folds of a_k*b_k + acc in the same order. Folds run in the working carrier of the element type; for float and double that is the type itself.

There is no pairwise or blocked fold inside an item. The one tree in aether is the 32-slot pair tree of the compensated reductions below.

Division by a scalar#

expr / s computes r = 1/s once and multiplies every component by it (operator/ in aether/expr/Operations.h). It is not a correctly rounded division per component: x * (1/s) can differ from x / s in the last bit, and a tiny s whose reciprocal overflows gives inf * x.

Per-lane linear algebra#

aether::linalg (aether/linalg/) factors and solves one small matrix per lane. A factorisation is not a fold of an item, so it states its own order:

  • Every sum runs in ascending index. LU reduces each entry by one finished dot product (acc += a_ik*a_kj, then a_ij -= acc), the order LAPACK’s getrf kernels use; the triangular solves and Cholesky subtract term by term (s -= l*x).

  • A division by a diagonal entry is a correctly rounded division, except LU’s multipliers, which multiply by one reciprocal of the pivot.

  • Nothing is compensated inside a factorisation.

  • The pivot is the first entry of largest magnitude in its column; a zero pivot is the singular test. Cholesky tests !(d > 0).

  • cholesky and eigh read 0.5*a_ij + 0.5*a_ji, so a symmetric input is read unchanged. eigh runs a fixed number of Jacobi sweeps (eighSweeps) in fixed order.

The same bits follow wherever the translation unit does not contract products, as above; on the device the values agree to a few ulps and info agrees exactly.

Other fixed spellings#

  • neg is x * -1.

  • norm is sqrt of the squared norm. There is no hypot.

  • Nothing promotes. A float expression stays in float and uses the f-suffixed libm; a double expression stays in double.

Compensated sums#

Plain folds are not compensated. Compensation is used where many contributions land in one slot or one value is reduced across lanes (aether/accum/atomic.h):

  • Two-sum. twoSumErr_(a, b, s) returns the exact residual (a + b) - s of a rounded sum, for any order and signs of a and b. Its inputs pass an optimisation barrier, so a product feeding them is rounded first and never fused into the subtractions.

  • Accumulation plane. AccumPlane with the default CompensatedAtomic policy adds each term to a value slot with an atomic add and adds the two-sum residual to a companion slot. The delivered value is value + comp. The components of a term are folded in ascending component order, host and device alike (see Accumulation).

  • Neumaier step and pair combine. neumaierAdd_ adds one term to a (hi, comp) pair. pairCombine_ merges two pairs: a two-sum of the hi parts, with the comp parts added; swapping its operands gives the same bits.

  • The 32-slot tree. reduce32Host_ combines 32 per-lane pairs with pairCombine_ in rounds of xor 16, 8, 4, 2, 1. The device butterfly32_ leaves each reduced column bit-equal to it; lanes past a ragged end carry +0 pairs.

Non-finite values#

  • The compensated helpers drop the correction once the high part is not finite, so the pair then follows IEEE through hi: a lone +Inf gives +Inf, +Inf with -Inf gives NaN, and NaN stays NaN.

  • isFinite of an expression is true only if every component is finite.

  • maxNorm folds a > b ? a : b over the absolute values, so which operand wins against a NaN follows that expression’s operand order.

  • math::fmax and fmin have C fmax/fmin semantics: a NaN operand is ignored when the other is not NaN.

  • math::sign of a NaN is 0.

Where bits are guaranteed#

The guarantee is checked, not assumed. Against a reference set of inputs and expected results (below), the plain-C core gives the same bits as the C++ templates on the host, under two conditions:

  • the C translation unit is compiled with -ffp-contract=off. GNU-mode GCC defaults to -ffp-contract=fast, which fuses a*b - c*d (the shape of cross, quat_mul and every row dot) into an FMA and changes the low bits;

  • the templates take their default std:: host route: no AETHER_HOST_VECTOR_MATH (which reroutes some functions to faithfully-rounded packet implementations) and no fast-math.

One result is outside the guarantee: the sign of a zero from min or max of +0 and -0. C leaves the sign of fmin(+0, -0) and fmax(+0, -0) unspecified, and compilers differ, so the check accepts either zero for those two kinds. A NaN is accepted for a NaN; every other result is a byte comparison.

Device code is not part of the bit guarantee. nvcc contracts a product and a sum into an FMA per statement, which the host build with -ffp-contract=off does not. A core whose target has no float64, or different fused multiply-add rules, states a ULP band against the reference set instead of bit identity.

The key cases for the plain-C core:

Operation

Generated arithmetic

sum, dot

right-nested, no seed: (a2*b2) + ((a1*b1) + (a0*b0))

norm

sqrt(ssq), where ssq is the right-nested sum of squares

div by scalar

r = 1/s once, then r * x_k per component

neg

x * -1

quaternion

a vec4 laid out (w, x, y, z); a 3-vector is a vec3

Language cores#

A target whose language cannot compile aether’s C++23 templates needs a language core: the contract above written out in that language. A core is generated, so one source states each rule, and it is checked against the same reference set as every other core. Plain-C reference core is the first.

The reference set is the committed corpus under tests/c_core/ (the umbrella aether_c_corpus.h and the md5-fenced parts in corpus/). The CCore gtest suite (tests/test_CCore.cpp, tests/test_CCore.cu) checks every generated function against the template expression it mirrors. A future core (C++17 for older toolchains and MSL, OpenCL C) adds its generator next to c_core and runs the same check.