Compensated Arithmetic in Core

Recent changes have exposed the raw compensated arithmetic primitives like two_sum and two_prod in Thermite core, whereas they were previously only in the thermite-compensated crate. Why? Turns out it's pretty essential to keep the basic compensated ops around at the very low level.

# What is Compensated Arithmetic?

Compensated arithmetic primitives, more formally known as Error-Free Transforms, are a way of recovering the error lost by rounding when applying basic operators such as addition, subtraction, multiplication, and division. Almost everything else falls out from that.

a + b in an f64 isn't a + b. The exact sum can need far more bits than the mantissa has (1e300 + 1e-300 would take about two thousand), so the hardware rounds it to the nearest representable value and drops the rest. What gets dropped is the interesting part, because for addition it's always exactly representable in the same format. Knuth's 2Sum recovers it in six adds, no branches, no extra precision.

\begin{aligned}
s &= \mathrm{fl}(a + b) \\
e &= (a + b) - s \\
s + e &= a + b
\end{aligned}

$s$ is the result you would have gotten anyway, and $e$ is exactly what the rounding lost, so the pair $(s, e)$ holds the exact sum in two words. Nothing is lost, the value is just re-encoded across two floats.

Multiplication works out the same way. The exact product of two $p$-bit values needs $2p$ bits, and the residual left after rounding fits in $p$ of them, so it's recoverable as well. With a true fused multiply-add that's two instructions:

let p = a * b;
let e = a.mul_sub(b, p); // fma(a, b, -p)

Without FMA, you need Dekker's algorithm, which cuts both operands into half-width pieces with Veltkamp splitting and sums the four cross products in the right order. That's 17 operations instead of 2, so many algorithms prefer other methods when native FMA isn't available.

Everything else is built on those two. Double-double arithmetic (Compensated<f64>, around 106 bits of significand) is a two_sum/two_prod pair plus renormalization. Kahan and Neumaier summation are two_sum with the error fed back into the next iteration. A correctly-rounded mul_add on hardware without an FMA instruction is Boldo and Melquiond's construction, built out of two_prod and two_sum. Interval arithmetic needs the residual to know which way to nudge an endpoint. It's pretty handy all things considered.

One catch is that this buys precision, not range. The exponent is untouched, so anything that overflows or flushes to zero is still gone, and Veltkamp splitting has its own threshold ($2^{115}$ for f32) past which operands have to be rebalanced before they can be split at all.

# algebraic_scalar

Thermite has the algebraic_scalar Cargo feature to enable the "Algebraic" float ops (algebraic_* on f32/f64), which allow LLVM to casually reassociate operations to improve performance or improve autovectorization opportunities. With this feature enabled, all basic operations (+,-,*,/) are set to algebraic ops at the lowest library level, so everything inherits those everywhere. This is honestly great for most scalar code. If you are not using or cannot use SIMD for one reason or another, making the scalar fallbacks as fast as possible is ideal.

Summation with a division in it is the clearest example:

// NOTE: "Scalars" in Thermite are still wrapped in `Vector`
let mut sum: Vector<f64> = Vector::splat(0.0);
for &x in xs {
    sum += x / Vector::splat(3.0);
}

Strict, at -O -C target-cpu=x86-64-v3, the inner loop is unrolled 8x and every copy looks like this:

...
vmovsd  (%rax), %xmm2
vdivsd  %xmm1, %xmm2, %xmm2   # xmm1 = 3.0
vaddsd  %xmm2, %xmm0, %xmm0
# repeated 8 times

Enabling the algebraic_scalar feature gives this instead:

vbroadcastsd __real@3fd5555555555555(%rip), %ymm1   # 1/3
vfmadd231pd  (%rcx,%r10,8), %ymm1, %ymm0
vfmadd231pd  32(%rcx,%r10,8), %ymm1, %ymm2
vfmadd231pd  64(%rcx,%r10,8), %ymm1, %ymm3
vfmadd231pd  96(%rcx,%r10,8), %ymm1, %ymm4

The division is gone entirely, replaced by a multiply against a rounded $1/3$, which is less precise but "close enough" for most use cases. The add contracted into that multiply, and the serial accumulator turned into four independent ones at four doubles apiece, so each iteration retires 16 elements in 4 instructions where the strict version needed 32. With a runtime divisor it does the same thing and just hoists a single vdivsd out of the loop. This is a great trade for most code.

However, taking a look at two_sum:

fn two_sum(a: Self, b: Self) -> (Self, Self) {
    let s = a + b;
    let v = s - a;
    let e = (a - (s - v)) + (b - v);
    (s, e)
}

As real-number algebra, with $s = a + b$ and $v = s - a$, reassociation collapses $v = (a + b) - a = b$, and with it the whole error term becomes zero:

\begin{aligned}
e &= \bigl(a - (s - v)\bigr) + (b - v) \\
  &= (a - a) + (b - b) = 0
\end{aligned}

LLVM is free to do this in multiple places, including common Cody-Waite style range reductions that occur in almost everything.

For example, average-precision sin/cos reduces its argument modulo $\pi$ by splitting the constant across several floats, $\pi \approx \pi_A + \pi_B + \pi_C$, where $\pi_A$ has enough low mantissa bits zeroed that $n \pi_A$ is exact for every $n$ in range, and each following term picks up the bits the previous one couldn't hold:

let n = (x * FRAC_1_PI).round();
let r = ((x - n * PI_A) - n * PI_B) - n * PI_C;

The subtractions are deliberately sequenced so that each one cancels against a quantity of its own magnitude, keeping $r$ accurate to full precision even when $x$ is large and $n \pi$ is nowhere near representable. Algebraically, though:

\begin{aligned}
r &= x - n(\pi_A + \pi_B + \pi_C) \\
  &= x - n\pi
\end{aligned}

so LLVM is free to fold the three constants into one rounded $\pi$ and emit a single $x - n\pi$. The extra $\approx 50$ bits of the constant vanish, and the reduction degrades from full precision to whatever $x - n\pi$ happens to give. Not good.

exp and ln have similar range reduction. algebraic_scalar would silently break that when using the math library with scalars. The effects of it spread quite far.

For that reason, thermite-compensated previously would emit a compile_error! if used with the algebraic_scalar feature enabled. It would have been simply broken.

# The Solution is more than the two_sum of its parts

Fixing all of these weak points was simply a matter of moving two_sum/two_diff/two_prod/two_quot from thermite-compensated and putting them on FloatRegister/FloatVectorWithBits.

That way they can bypass the algebraic math intrinsics in the scalar backend implementation, using strict operators for compensated arithmetic, while also allowing anyone to use them anywhere that FloatVectorWithBits is available. That's actually an important distinction, because if it were only FloatVector, two_sum would have to be available on Complex or Dual composite types that implement FloatVector. Compensated arithmetic on those only makes sense as an implementation detail, not on the public interface.

Furthermore, FMA, even mul_adde (opportunistic FMA: fused where the hardware has it, separate mul+add where not) and family, inherently use strict operators as well so they cannot be reassociated.

Although not taking advantage of algebraic reassociation, LLVM is smart enough to see that a.two_sum(b).0 (accessing only the sum, ignoring the error) can ignore the error calculation entirely, returning a + b, but with a strict operator. In a way it provides an escape hatch for algorithms to ensure that the sum is exactly that, not mixed in with nearby values. Same with the others.

This pattern and just using two_sum/two_prod deep in the math implementations helped improve or stabilize precision across a wide variety of functions, and fixed a few latent bugs that would have crippled the emulated FMA paths. Range reduction, accumulators, etc., are now fixed and stable under the feature.

Originally, this change was just because I wanted thermite-compensated to compile regardless of what crate features users used. However, I've since found quite a few places in Thermite's math libraries where either strict operators or compensated arithmetic significantly improve precision, nearly for free in many cases, and the autovectorizer can do an even better job everywhere if that feature is enabled. Win-win.