Guarantees and limits#

What ExactRef promises, what its tolerances cover and what they do not, and the practical limits of using it.

References#

Every output is the exact result rounded once to the output format, or the call raises er.ReferenceUnavailableError. ExactRef never returns an answer it has not certified.

Results with transcendental values, such as exp in attention, are refined until the rounding is certain, which happens before the precision cap in practice. The one documented exception is attention's backward: grad_key, grad_value, and grad_query summed over broadcast copies of a query row, add terms from rows with different softmax denominators. They are found exact only when every term is rational, so such a sum that is exactly a rounding boundary, or exactly zero by cancellation, reaches the cap and raises.

What a tolerance covers#

tolerance() covers every implementation that loses precision only at the operation's required rounding steps, rounding to nearest there:

  • one rounding in the compute format at each required step,

  • the additions in normal accumulation (index order, or 2 or 4 SIMD-style lanes), or with standard="any_order" in any order that does not overflow,

  • library functions such as exp within ulps units in the last place (default 2),

  • and the final rounding to the output format.

The rounding argument of an operation sets the rounding of the reference and of that final step, not of the steps inside the computation.

An output within tolerance() is no worse than normal accumulation can be. Beyond it but within the any-order tolerance, it is correct but worse than normal accumulation can be. Beyond the any-order tolerance, no implementation that rounds only at the required steps, without overflowing, could produce it. Where an overflow is possible, may_be() says which infinities or NaN a correct implementation may produce, erring on the side of yes.

What a tolerance does not cover#

  • Roundings an implementation adds beyond the required ones, such as FlashAttention's cast of the weights to bfloat16, or an intermediate result stored in the dtype. Such an implementation can exceed tolerance(). Algorithms and chains budget these.

  • Materialized attention probabilities: PyTorch's math path normalizes and rounds P before the product with V. No tolerance budgets that.

  • Every bug. A wrong formula can still produce an output within the bounds, for example when its effect is smaller than the rounding allowed on those inputs. A tolerance says whether rounding can explain a difference.

Matrix exponential#

matrix_exp has no rounding step that every algorithm shares, so its tolerance is a backward error instead: a correct result is exp(A + dA), rounded to the compute format and then to the output format, with

||dA||_inf <= ulps * (n + 16) * u * max(||A||_inf, 1)

for n x n matrices, where u is the compute format's unit roundoff and ulps defaults to 1. The n and the 16 count the roundings of an algorithm's inner products and its other steps (short sums, scalings, the linear solve), and the floor of 1 covers the rounding of the identity part of exp(A). matrix_exp_frechet and matrix_exp's backward use 2n + 16 for their 2n x 2n block matrix. There is no normal answer.

The allowance is an assumption, measured on random matrices (n from 2 to 64, ||A||_inf from 0.03 to 32): Higham's 2005 scaling and squaring with its degree-13 Pade approximant, and PyTorch's matrix_exp, stay within half of it. Some results fall outside it: PyTorch's own at small norms (float64 near 0.03, float32 up to 1), and SciPy 1.18.1's expm for some float64 matrices from ||A||_inf 4. There the reference matches mpmath at high precision.

Formats without infinity#

float8_e4m3fn and the fnuz formats have no infinity. A finite result past the largest value rounds within the format's values: to the largest value when rounding to nearest, toward zero or toward it, and to NaN when rounding outward, where the format has no value. PyTorch's finite float8_e4m3fn results agree. PyTorch also turns infinities into the largest value, which overflow="saturate" reproduces; its conversions into the fnuz formats give NaN instead, which overflow="nan" reproduces. exact keeps the exact value of every element, so behind a NaN it shows whether the true value was a finite number too large for the format, an infinity, or NaN.

Special values#

Infinities, NaN and signed zeros follow IEEE 754 applied to the exact expression. Where mathematics leaves the answer open, ExactRef follows PyTorch; where mathematics gives an answer, ExactRef gives it, even when PyTorch differs. Each operation's reference page states its rules.

Threads#

python-flint's working precision is a process-wide setting, so do not call ExactRef from several threads at once. Separate processes are fine, including pytest-xdist's workers. ExactRef sets FLINT's own thread count (up to 8), also a process-wide setting; no result depends on it.

The floating-point environment#

ExactRef computes in the default floating-point environment. When the calling thread flushes subnormals, for example after torch.set_flush_denormal(True), each call turns that off and restores it before returning. Where the platform offers no way to do that, the call raises ReferenceUnavailableError.

Proofs#

Every guarantee and tolerance is proven in plain-text Markdown files in the repository:

ExactRef's answers are also tested against independent references: MPFR and MPC (through gmpy2), SymPy, NumPy, PyTorch and the Python standard library.