Get started#

This page takes you from installing ExactRef to a test that checks an implementation and says what went wrong when it fails.

Install#

pip install exactref

ExactRef needs Python 3.10 or newer. This also installs python-flint, its only dependency, which supplies proven bounds for functions such as exp and exact integer matrix products. NumPy and PyTorch are optional: ExactRef reads their arrays and tensors when you pass them, and needs neither.

Your first reference#

Code
import exactref as er

# A 2x2 and a 2x1 float32 matrix
A = er.asarray([[1.0, 2.0], [3.0, 4.0]], "float32")
B = er.asarray([[0.1], [0.2]], "float32")

# The reference for A @ B
r = er.matmul(A, B)
# The exact product, rounded to float32
print("rounded:", r.rounded.tolist())
# rounded is an ExactRef array, with a dtype and a shape
print("dtype:", r.rounded.dtype.name)
print("shape:", r.rounded.shape)
# Per element, the most a correct float32 output can differ from rounded
print("tolerance:", r.tolerance().tolist())
Output
rounded: [[0.5], [1.100000023841858]]
dtype: float32
shape: (2, 1)
tolerance: [[5.960464477539063e-08], [1.1920928955078125e-07]]

er.asarray rounds 0.1 and 0.2 to the nearest float32 values. r.rounded is then the exact product of those float32 values, rounded once to float32: the best answer float32 can hold. r.tolerance() says, per element, how far a correct float32 implementation can land from it.

Every operation works the same way. er.sum adds over the dimensions you name:

Code
# The sums down A's columns (over axis 0), rounded to float32
print("column sums:", er.sum(A, axis=0).rounded.tolist())
Output
column sums: [4.0, 6.0]

Bring your data#

Operations accept arrays from anywhere, and er.asarray imports them explicitly:

  • Python values and nested lists. Without a dtype, floats become float64, ints int64, complex numbers complex128 and bools bool. With one, each value is rounded once to that format.

  • NumPy arrays and anything else with the buffer protocol, imported exactly in their own dtype.

  • PyTorch tensors (and JAX arrays) on the CPU, imported exactly through DLPack, bfloat16 and float8 included. A tensor that requires grad is detached first.

Code
import numpy as np
import torch

# A NumPy float32 array, imported exactly as float32
print("numpy:", er.asarray(np.array([0.1, 0.2], dtype=np.float32)).tolist())

# A PyTorch bfloat16 tensor, imported exactly as bfloat16
x = er.asarray(torch.tensor([1.0, 0.1, 3.3], dtype=torch.bfloat16))
print("dtype:", x.dtype.name)
# torch.tensor rounded 0.1 and 3.3 to the nearest bfloat16 values
print("values:", x.tolist())
Output
numpy: [0.10000000149011612, 0.20000000298023224]
dtype: bfloat16
values: [1.0, 0.10009765625, 3.296875]

A tensor on a GPU has to be copied to the CPU first: er.asarray(t.cpu()). Passing it directly raises ValueError: ... is not CPU memory; copy the tensor to the CPU first.

To convert while importing, give a dtype, and rounding if you want other than round-to-nearest-even:

Code
# 0.1 rounded to the nearest bfloat16 value, the default
nearest = er.asarray([0.1], "bfloat16")
print("nearest:", nearest.tolist())

# 0.1 rounded toward zero, to the bfloat16 value just below it
toward_zero = er.asarray([0.1], "bfloat16", rounding=er.Rounding.TOWARD_ZERO)
print("toward zero:", toward_zero.tolist())
Output
nearest: [0.10009765625]
toward zero: [0.099609375]

Check an implementation#

Run the operation on the same inputs your implementation got, then ask the result whether each output element is one a correct implementation can produce. Here the implementation is a plain Python loop:

Code
# Two lists of Python floats
xs, ys = [0.3, 0.6, 0.9], [3.0, 2.0, 1.0]

# The implementation: a dot product, rounding each product and sum to float64
total = 0.0
for x, y in zip(xs, ys):
    total += x * y

print("loop:", total)

# 1-D by 1-D: a dot product
r = er.matmul(xs, ys)
# The exact answer rounded once to float64, and its tolerance
print("rounded:", r.rounded.tolist())
print("tolerance:", r.tolerance().tolist())
# Can a correct implementation give the loop's output?
print("covered:", r.covers(total).all())
Output
loop: 2.9999999999999996
rounded: 3.0
tolerance: 8.881784197001252e-16
covered: True

The loop is not exact, but it is correct. Rounding each product and each partial sum to float64 keeps every correct implementation's output within this range, and the loop's output is inside it:

Code
# finite_range() returns the low and high ends as two arrays
print("finite range:", [end.tolist() for end in r.finite_range()])
Output
finite range: [2.999999999999999, 3.000000000000001]

covers() checks each finite element against finite_range(), the lowest and highest finite outputs a correct implementation can produce (at most tolerance() away from rounded). An infinity or NaN passes only where r.may_be("+inf"), r.may_be("-inf") or r.may_be("nan") says a correct implementation can produce it, for example near the top of the format's range. Reading complicated tolerances shows these cases.

The output you pass must have the result's dtype and shape. A PyTorch tensor from the implementation already does. Python floats are float64, so for a float32 result, convert them first:

Code
# The float32 product from the start of this page
r32 = er.matmul(A, B)
# Python floats, converted to float32 to match the result's dtype
out32 = er.asarray([[0.5], [1.100000023841858]], "float32")
print("covered:", r32.covers(out32).all())
Output
covered: True

Passing the Python floats directly raises ValueError: output must have the result's dtype float32 and shape (2, 1), got float64 (2, 1).

Say how the implementation computes#

By default, tolerance() and covers() describe the most common kind of implementation: one that keeps its products and partial sums in the usual compute format, adds them by normal accumulation, and keeps subnormals (the tiny values below a format's smallest normal number, which some hardware modes flush to zero):

Assumption

Default

Change it with

the compute format

float32 for inputs narrower than float32 (bfloat16, float16, float8), else the input format

compute="bfloat16", or any other format

the order of additions

index order, or 2 or 4 interleaved running sums as SIMD code keeps them

standard="any_order": any order, any summation tree, with or without fused multiply-add

subnormals

kept

flush_to_zero=True, as torch.set_flush_denormal(True) does

When your implementation differs, say so. Here is a bfloat16 matrix product computed two ways: by PyTorch's own @, and by an implementation that keeps its running sums in bfloat16.

Code
import random

# Seeded, so every run uses the same values
rng = random.Random(0)

# A rows x cols bfloat16 tensor of random values
def randn(rows, cols):
    # Normally distributed, mean 0 and standard deviation 1
    values = [[rng.gauss(0, 1) for _ in range(cols)] for _ in range(rows)]
    return torch.tensor(values, dtype=torch.bfloat16)

# A random 64x256 and 256x32 bfloat16 matrix
A, B = randn(64, 256), randn(256, 32)

# The reference for A @ B
r = er.matmul(A, B)
# PyTorch's own @, judged under the default float32 compute format
print("A @ B covered:", r.covers(A @ B).all())

# An implementation that keeps its running sums in bfloat16
acc = torch.zeros(64, 32, dtype=torch.bfloat16)
for k in range(256):
    # Add A[i, k] * B[k, j] to every element (i, j), rounding to bfloat16
    acc = acc + A[:, k:k + 1] * B[k:k + 1, :]

# Judged under the default float32 compute format
print("bfloat16 sums covered:", r.covers(acc).all())
# Judged as computing in bfloat16
print("bfloat16 sums covered, bfloat16 compute:", r.covers(acc, "bfloat16").all())
Output
A @ B covered: True
bfloat16 sums covered: False
bfloat16 sums covered, bfloat16 compute: True

Under the default float32 compute format, the sums are so accurate that 2,004 of these 2,048 outputs have a tolerance of 0: every implementation that adds in float32 the normal way lands exactly on the reference. Sums kept in bfloat16 drift much further, and tolerance("bfloat16") reaches 100.0 for the worst element.

On the operation itself, dtype= sets the output format and rounding= the rounding of the reference and of the final step into the output format. The tolerance assumes round-to-nearest for the roundings inside the computation.

tolerance(), finite_range(), may_be() and covers() all take the same options, so you can look at the numbers behind a verdict.

Put it in a test#

In a test, assert r.covers(out).all() is enough to pass or fail. To see what failed, copy a helper like this one into your tests. It lists the elements no correct implementation produces:

Code
import exactref as er


def check(result, output, **options):
    """Fail, listing the elements of output no correct implementation produces.

    For a result with one real output; options are covers()'s, such as
    compute="bfloat16" or standard="any_order".
    """
    # Per element, whether a correct implementation can produce output
    ok = result.covers(output, **options)
    if ok.all():
        return

    # The shape, then the output, reference and tolerance as flat row-by-row lists
    shape = result.rounded.shape
    out = er.asarray(output).scalars()
    ref = result.rounded.scalars()
    tol = result.tolerance(**options).scalars()
    # The lowest and highest finite outputs a correct implementation can produce
    low, high = (end.scalars() for end in result.finite_range(**options))
    # The flat positions of the elements that failed
    bad = [i for i, good in enumerate(ok.scalars()) if not good]

    # The report: how many elements failed, then the first five of them
    lines = [f"{len(bad)} of {len(ref)} elements no correct implementation"
             " produces:"]
    for i in bad[:5]:
        # Its index, the output, and the reference with its tolerance and range
        lines.append(f"  {unravel(i, shape)}: got {out[i]}, reference {ref[i]} "
                     f"+/- {tol[i]:.3g}, range [{low[i]}, {high[i]}]")

    raise AssertionError("\n".join(lines))


def unravel(i, shape):
    """The index of flat position i in an array of this shape."""
    index = []
    for n in reversed(shape):
        # k is the index on this axis; i carries on to the axes before it
        i, k = divmod(i, n)
        index.append(k)

    # The indices were found last axis first
    return tuple(reversed(index))

A test then reads like this:

Code
def test_matmul_bfloat16():
    # Check PyTorch's A @ B against ExactRef's reference for it
    check(er.matmul(A, B), A @ B)

On the implementation that keeps its sums in bfloat16, the helper reports:

Code
# check() raises AssertionError with its report as the message
try:
    check(r, acc)
except AssertionError as error:
    print(error)
Output
1847 of 2048 elements no correct implementation produces:
  (0, 0): got -13.125, reference -13.0 +/- 0, range [-13.0, -13.0]
  (0, 1): got 2.375, reference 2.109375 +/- 0, range [2.109375, 2.109375]
  (0, 2): got 39.25, reference 39.75 +/- 0, range [39.75, 39.75]
  (0, 3): got -14.5625, reference -15.0 +/- 0, range [-15.0, -15.0]
  (0, 4): got 13.25, reference 13.0 +/- 0, range [13.0, 13.0]

pytest and pytest-xdist both work: pytest-xdist runs tests in separate processes, which is fine. Do not call ExactRef from several threads of one process at once, because python-flint's working precision is shared by the whole process.

When a check fails#

  1. Look at the failing elements. The helper above lists them with the reference, the tolerance and the range a correct implementation can reach.

  2. Check against the any-order standard, covers(out, standard="any_order"). If that passes, the implementation is correct but adds in an order that is worse than ordinary kernels can be.

  3. Check what you told ExactRef. Does the implementation compute in bfloat16 or float16 rather than float32? Does it flush subnormals? Does your GPU use TF32 for float32 matrix products? (Round the inputs to TF32 first and pass those.)

  4. Look for roundings the operation does not require. An implementation that stores an intermediate result, or a fused kernel that rounds its weights to bfloat16, rounds more often than the definition. Pass the earlier result in as an input (Chains and any_of), or use the matching algorithm, so the tolerance budgets those roundings too. PyTorch equivalents lists where PyTorch does this.

  5. Compare with r.normal(), what an ordinary implementation that adds in index order computes, emulated exactly.

If you are sure the implementation is right and the reference is wrong, open an issue with the inputs and your exactref version.

If an operation raises ReferenceUnavailableError, ExactRef could not certify the exact answer for those inputs, and it says so rather than guess.

Where to go next#