Feature: Hardening E4M3, E5M2, F16, & BF16 dot-products

[ashvardanian]

It's very easy to come up with a pair of E5M2 inputs that produce absurd numerical results, even when the inputs are upcast to BF16 and F32 is used for the actual product and accumulation. This is a fundamental and broad problem, but here is a real example from the CI on a tiny random-generated 8-dimensional input:

	i = 0	i = 1	i = 2	i = 3	i = 4	i = 5	i = 6	i = 7
aᵢ	0.00122	20480	−0.00122	1.5	−0.00586	−3072	−640	0.00146
bᵢ	−40	320	−1280	−7.63e⁻⁵	0	0.000427	10240	−4.58e⁻⁵
aᵢ·bᵢ	−0.04883	6553600	1.5625	−0.000114	0	−1.3125	−6553600	≈ 0

Why F32 accumulation fails here. The accurate sum of these 8 products is ≈ 0.201. After two vfmaq_f32 calls, the 4 accumulator lanes hold pairwise products: lanes 1 and 2 carry values around ±6.5 M. At that magnitude the F32 ULP is 0.5 — so the small meaningful terms (−0.049, 1.563, −1.313, −0.0001) are all below one ULP and get absorbed during pairwise reduction. The large terms then cancel exactly to zero, and the information is gone. Final F32 result: 0.0 instead of 0.201.

Sadly, this issue is not limited to E5M2. Pretty much all of E4M3, E5M2, F16, & BF16 types use similar multi-precision mechanisms — constrained by F32 precision.

That said, F32 dot-products use F64 numerics and aren't affected. F64 dot-products use "Dot2"-like schemes for stable multiplication and addition and aren't affected either.

There are a few ways to mitigate this:

1. F64 accumulators for all of E4M3, E5M2, F16, & BF16. Pretty much no modern hardware properly accelerates dot-products with 64-bit accumulators except for the smef64 capability level on Apple M4 chips. That said, even there — the upcast will be manual.
2. Separate accumulators for small & large magnitudes and positive & negative products. Such methods are often leveraged in traditional HPC environments. Still, they occupy too many registers for the "accumulator state", making them inapplicable to GEMM-like tiled kernels.
3. Exact integer arithmetic for E_x_M_y_ dot-products for low x and y. We already leverage the same approach for E2M3 via I8 and E3M2 via I16, leveraging trivial algebraic transforms.

For now, we'll just accept the catastrophic errors in some cases, assuming they don't affect applications like USearch, where many products with neighboring nodes will be computed before picking a new traversal direction. But if new functionality for cheaper widening and multi-precision numerics emerges, it should be considered for the next releases.

Can you contribute to the implementation?

• I can contribute

Is your feature request specific to a certain interface?

It applies to everything

Contact Details

No response

Is there an existing issue for this?

• I have searched the existing issues

Code of Conduct

• I agree to follow this project's Code of Conduct

[ashvardanian]

The last approach is the most promising, but still way too slow. Here's a draft of a possible solution for AVX-512:

__attribute__((target("avx512f,avx512bw,avx512vl")))
float dot_e4m3_exact_icelake(const unsigned char *a, const unsigned char *b, int n) {
    __m512i accumulator_even_i64x8 = _mm512_setzero_si512();
    __m512i accumulator_odd_i64x8 = _mm512_setzero_si512();

    __m512i const implicit_bit_i32x16 = _mm512_set1_epi32(8);
    __m512i const exponent_mask_i32x16 = _mm512_set1_epi32(0xF);
    __m512i const mantissa_mask_i32x16 = _mm512_set1_epi32(0x7);
    __m512i const sign_mask_i32x16 = _mm512_set1_epi32(0x80);
    __m512i const min_shift_i32x16 = _mm512_set1_epi32(1);
    __m512i const zero_i32x16 = _mm512_setzero_si512();

    for (int i = 0; i + 16 <= n; i += 16) {
        __m128i a_raw_u8x16 = _mm_loadu_si128((const __m128i *)(a + i));
        __m128i b_raw_u8x16 = _mm_loadu_si128((const __m128i *)(b + i));
        __m512i a_raw_i32x16 = _mm512_cvtepu8_epi32(a_raw_u8x16);
        __m512i b_raw_i32x16 = _mm512_cvtepu8_epi32(b_raw_u8x16);

        // Extract exponent and mantissa fields
        __m512i a_exponent_i32x16 = _mm512_and_si512(_mm512_srli_epi32(a_raw_i32x16, 3), exponent_mask_i32x16);
        __m512i a_mantissa_i32x16 = _mm512_and_si512(a_raw_i32x16, mantissa_mask_i32x16);
        __m512i b_exponent_i32x16 = _mm512_and_si512(_mm512_srli_epi32(b_raw_i32x16, 3), exponent_mask_i32x16);
        __m512i b_mantissa_i32x16 = _mm512_and_si512(b_raw_i32x16, mantissa_mask_i32x16);

        // Add implicit bit for normals (exp != 0): significand = M + 8
        __mmask16 a_normal_mask = _mm512_cmpneq_epi32_mask(a_exponent_i32x16, zero_i32x16);
        __mmask16 b_normal_mask = _mm512_cmpneq_epi32_mask(b_exponent_i32x16, zero_i32x16);
        __m512i a_significand_i32x16 = _mm512_mask_add_epi32(a_mantissa_i32x16, a_normal_mask, a_mantissa_i32x16, implicit_bit_i32x16);
        __m512i b_significand_i32x16 = _mm512_mask_add_epi32(b_mantissa_i32x16, b_normal_mask, b_mantissa_i32x16, implicit_bit_i32x16);

        // Shift: max(exp, 1) — normals shift by exp, subnormals by 1
        __m512i a_scaled_i32x16 = _mm512_sllv_epi32(a_significand_i32x16, _mm512_max_epu32(a_exponent_i32x16, min_shift_i32x16));
        __m512i b_scaled_i32x16 = _mm512_sllv_epi32(b_significand_i32x16, _mm512_max_epu32(b_exponent_i32x16, min_shift_i32x16));

        // Combined sign: negate b where (a ^ b) has sign bit set. a_scaled stays non-negative, b_scaled carries the product sign.
        __mmask16 negate_mask = _mm512_test_epi32_mask( _mm512_xor_si512(a_raw_i32x16, b_raw_i32x16), sign_mask_i32x16);
        b_scaled_i32x16 = _mm512_mask_sub_epi32(b_scaled_i32x16, negate_mask, zero_i32x16, b_scaled_i32x16);

        // VPMULDQ: signed i32 × i32 → i64, even/odd split for all 16 lanes
        accumulator_even_i64x8 = _mm512_add_epi64(accumulator_even_i64x8, _mm512_mul_epi32(a_scaled_i32x16, b_scaled_i32x16));
        accumulator_odd_i64x8 = _mm512_add_epi64( accumulator_odd_i64x8, _mm512_mul_epi32(_mm512_srli_epi64(a_scaled_i32x16, 32), _mm512_srli_epi64(b_scaled_i32x16, 32)));
    }

    __m512i total_i64x8 = _mm512_add_epi64(accumulator_even_i64x8, accumulator_odd_i64x8);
    int64_t full_sum = _mm512_reduce_add_epi64(total_i64x8);
    return (float)((double)full_sum * 9.5367431640625e-07); // 2^(-20)
}

Key idea. Each E4M3 value is decoded into an exact integer: significand << exponent. Normals become (8 + M) << E, subnormals become M << 1. The maximum per-element product is (15 << 15)² ≈ 2.4 × 10¹¹, which fits in 38 bits — well within i64. VPMULDQ gives us the exact i32 × i32 → i64 multiply, and the i64 accumulator never overflows for practical vector lengths (up to ~37M elements).

[ashvardanian]

I've explored several alternative exact E4M3 dot product strategies on a Sapphire Rapids CPU. All produce zero error vs. the F64 reference.

Strategy	Elem/iter	Accum state	cy/elem	GB/s
Decode each E4M3 → exact i32, variable-shift by exponent, multiply pairs via i32×i32→i64	16	2 × ZMM i64 (128 B)	0.53	7.9
Group by exponent class (7 bins), use byte LUT + u8×i8→i16 narrow multiply within each bin	64	7 × ZMM i32 (448 B)	0.44	9.5
Same idea but split exponent into nibbles (15 bins) — more bins, more overhead	64	15 × ZMM i32 (960 B)	0.82	5.1
Decompose into 4×4 exponent×mantissa groups, use i16×i16→i32 pair multiplies	64	7 × ZMM i32 (448 B)	0.86	4.9
Track F32 rounding error via 2Sum identity, accumulate correction separately	16	2 × ZMM f32 (128 B)	0.95	4.4
Naive F32 FMA — upcast to F32, multiply-add (not exact)	16	1 × ZMM f32 (64 B)	0.87	4.8
Hardware F8→F16 upcast via VCVTNE2PH, F32 FMA (not exact, Genoa only)	32	1 × ZMM f32 (64 B)	0.32	13.2

VPSLLVD + VPMULDQ (exact_icelake) — The conceptually simplest approach. Each E4M3 byte becomes an exact integer via variable shift; products are computed in i64. The 5-cycle VPMULDQ latency is the main bottleneck; 2× loop unrolling hides it via independent dependency chains (0.62 → 0.53 cy/elem). The tiny 128-byte accumulator state makes this the most GEMM-friendly exact kernel.

7-bin VPSHUFB (7bin_skylake) — Groups elements by E_hi = E >> 2 (4 exponent classes per input, 7 sum classes). Within each bin, all products share a common scale factor, so narrow VPSHUFB + VPMADDUBSW suffices. Fastest dot product throughput (0.44 cy/elem), but the 7-accumulator state is prohibitive for GEMM tiling: a 4×4 output tile alone requires 112 ZMM registers.

15-bin nibble decomposition (nibble_skylake) — Splits the exponent by E_hi = E >> 1, yielding 15 bins. More bins means more overhead with diminishing returns. Not competitive.

Compensated 2Sum (compensated_skylake) — Tracks both the main F32 sum and the rounding error using the 2Sum identity. Compact state (2 ZMM registers), but the extra additions in the error path add ~10% overhead over naive FMA, and it's still slower than the integer approaches.