1 Introduction↩︎

Multi-component multiple-precision arithmetic bundles several binary64 (or binary32) components and handles carries exactly through error-free transformations (EFTs); double-word (DW), triple-word (TW), and quad-word (QW) formats—implemented as DD/TD/QD and, with a single-precision base, DS/TS/QS—belong to this family. Because it reuses the existing floating-point units directly, it tends to achieve much higher performance than arbitrary-precision libraries. Its internals, however—in particular the renormalization stages—contain if-branches that depend on the magnitudes of intermediate values, and these branches have been a major performance obstacle for SIMD and GPU implementations, where they manifest as lane or warp divergence.

Zhang and Aiken proposed branch-free algorithms for DW/TW/QW addition and multiplication that eliminate these branches, together with machine-checked error bounds obtained with an SMT solver, and demonstrated high performance [1], [2]. We previously implemented the TW and QW variants, for which large speedups over the conventional algorithms can be expected, and confirmed that they indeed contribute to acceleration.

In this paper we propose, for further acceleration, fused multiply-add (FMA) operations \(z\leftarrow x\cdot y+c\) for DW/TW/QW (and for DS/TS/QS). The proposed FMA requires fewer operations than the composition of branch-free addition and multiplication (DW \(17\) / TW \(66\) / QW \(146\) flops versus \(29/96/209\) flops), its error bounds and the validity of every FastTwoSum precondition are machine-verified with FPANVerifier, and it is bitwise commutative under the exchange \(x\leftrightarrow y\). As a reference that measures the attainable accuracy, we also provide an Exact FMA variant with a result-relative error bound.

The paper is organized as follows. Section 2 presents the proposed algorithms and their machine-verified error bounds; Section 3 contrasts them with existing FMA research; Section 4 reports benchmark tests on two platforms, GB10 (Arm) and H100 (x86)—the acceleration of division and square root, and the basic linear-algebra kernels AXPY/GEMV/GEMM in serial, OpenMP-parallel, and GPU-parallel settings. Section 5 concludes with future work.

2 Branch-free FMA algorithms for DW, TW, and QW arithmetic↩︎

The proposed branch-free FMA computes \(z\leftarrow x\cdot y+c\) in a single fused stage for \(K\)-component operands \(x,y,c\) (\(K=2,3,4\) for DW/TW/QW), and has the following properties.

• Commutativity: it is commutative in \(x\) and \(y\); \(\mathrm{fma}(x,y,c)\) and \(\mathrm{fma}(y,x,c)\) agree bitwise for all six types (Section 2.10).

• Machine-verified error bounds: the bounds \(34u^2/184u^3/812u^4\cdot(|xy|+|c|)\) and the validity of every FastTwoSum precondition are machine-verified with FPANVerifier for all precisions \(p\) (Section 2.9).

• Fewer operations: it needs only \(17/66/146\) flops, fewer than the combination of branch-free addition and multiplication (\(29/96/209\) flops; Section 2.8).

• Error semantics: the proposed variant is bounded relative to \((|xy|+|c|)\), whereas the Exact variant is bounded relative to the result \(|xy+c|\) (Section 2.7).

2.1 Basic operations (error-free transformations)↩︎

The proposed FMA is built from only three error-free transformations (EFTs) listed in Table 1. Each decomposes its input into an unrounded two-term sum or product. FastTwoSum is 3 flops cheaper than TwoSum, but it is exact only when \(\mathrm{exponent}(a)\ge\mathrm{exponent}(b)\) (this precondition is the crux of the machine verification in this paper). TwoProd obtains the rounding error of a product in \(2\) flops using the hardware FMA.

2.2 Overall structure of the proposed FMA↩︎

For \(K\)-component (\(K=2,3,4\), i.e.DW/TW/QW) operands \(x,y,c\), the proposed FMA computes \(z=x\cdot y+c\) in three stages: (i) product expansion (TwoProd of the \(x_iy_j\)), (ii) weight-wise accumulation (terms arriving at the same power-of-two position are gathered by TwoSum chains), and (iii) renormalization (reshaping into \(K\) non-overlapping components). Terms of order \(u^{K+1}\) and below are truncated; the resulting error bound is machine-verified in Section 2.9. Algorithm 1 shows the DW case, Algorithm 3 the TW case, and Algorithm 5 the QW case.

2.3 DW-FMA (proposed)↩︎

The data flow of Algorithm 1 is shown in Figure 2. Each TwoProd/TwoSum/FastTwoSum is one gate; solid lines carry the high-order words and dashed lines the error (low-order) words. The truncated term \(x_1y_1\) is not drawn (\(O(u^2)\)). The output is made non-overlapping by the single trailing FastTwoSum gate (this is the source of the \(20\to17\) flops reduction).

2.4 TW-FMA (proposed)↩︎

For \(K\ge3\) the terms of the product expansion fall naturally into weight levels (\(u^0,u^1,u^2,\dots\)). The skeleton of the proposed FMA is a triangular structure that accumulates each level with a TwoSum chain while preserving errors, and drops each spilled error one level down; finally the level representatives \(B,A_1,A_2,\dots\) are made non-overlapping by a single renormalization chain. Figures 4 (TW) and 6 (QW) illustrate this. The key point is that the error words of TwoSum (dashed) always flow to the next lower level.

2.5 QW-FMA (proposed)↩︎

2.6 Choice of renormalization gates and provability↩︎

The key design principle of the proposed FMA is to choose each gate of the renormalization stage according to whether its FastTwoSum precondition can be machine-proved for all precisions \(p\). In TW/QW, the gates right after accumulation, \((B,m_1)\) and \((B,A_1)\), admit a proof of \(\mathrm{exp}(B)\ge\mathrm{exp}(\cdot)\) and can therefore be FastTwoSum, whereas the gates just before the output, \((w_0,w_1)\) and \((\rho,w_2)\), cannot be proved and are kept as TwoSum (Section 2.9). This placement minimizes flops while preserving the property that “FastTwoSum is used only where its precondition is proved.”

2.7 The Exact FMA (distillation variant): algorithm and error bound↩︎

As a counterpart of the proposed FMA we use a result-relative Exact FMA (distillation FMA). It is a reference for measuring the attainable accuracy and not the main proposal of this paper, but the contrast is essential for understanding the error semantics of the proposed FMA.

2.7.0.1 The only error source is the final truncation.

The decisive difference from the proposed FMA is where information is discarded.

• Proposed FMA: high-order products (\(i+j\ge K\)) are never formed, and the diagonal and error sums are rounded by plain additions. The discarded amount is governed by the magnitudes of \(|xy|\) and \(|c|\), so the bound is relative to \((|xy|+|c|)\) (\(34u^2/184u^3/812u^4\); Section 2.9).

• Exact FMA: step (1) of the expansion is exact. Since \(x\cdot y=\sum_{i,j}x_iy_j\) and each \(x_iy_j\) is decomposed exactly into \(P_{ij}+E_{ij}\) by TwoProd, the identity \(\textstyle\sum v = x\cdot y + c\) holds without any rounding error. Steps (2) and (3) use only TwoSum, which is an EFT, so the sweeps preserve the total of the multiset. Consequently the only error is the final truncation of \(v_K,\dots,v_{M-1}\): \[z-(x\cdot y+c) \;=\; -\sum_{i\ge K} v_i .\]

The sweeps arrange the expansion into a nearly non-overlapping form (\(|v_{i+1}|\lesssim u|v_i|\)), so the discarded tail is governed by \(v_0\), i.e.by the result, giving a bound of the form \[|z-(x\cdot y+c)| \;\le\; C_K\,u^K\,|z| \qquad(\text{\boldsymbol{result-relative}}).\] The constant \(C_K\) stems from the finiteness of the sweep count \(n\) (infinitely many sweeps would converge to Priest’s non-overlapping expansion with \(C_K\to O(1)\)).

2.7.0.2 Machine verification.

The verification scripts (dwfma_exact_sym.fpan etc.) bound each discarded wire after distillation with the result \(z_0\) as the anchor (in contrast to the .fpan files of the proposed FMA, which anchor to \(\max(|xy|,|c|)\)). For DW, bounds such as |P01| <= u^2 |P00| are proved, confirming that the sum of the discarded wires stays within \(O(u^2|z_0|)\). For TW/QW, whose expansions have 21/36 terms, tightening every wire has not yet completed within practical time (future work).

2.7.0.3 Measurements.

For the single operation \(x\cdot y+c\) (\(200{,}000\) trials each, against MPFR 600-bit), the relative errors are

For DD the maximum error stays within \(\varepsilon_2/2\), i.e.the result is effectively correctly rounded. The mean for TD/QD is also below \(0.2\varepsilon_K\), but the maximum reaches \(10^3\)–\(10^2\) times \(\varepsilon_K\). This is due to the finite sweep count (\(n=3,4\)): trials with strong cancellation do not reach a fully non-overlapping expansion and a tail remains (the constant \(C_K\) above). More sweeps would shrink it at the cost of speed—this is why this paper positions the Exact FMA as a reference for the attainable accuracy and proposes the cheap FMA as the practical replacement.

2.7.0.4 What result-relative means.

Under strong cancellation \(c\approx-xy\) we have \(|z|\ll|xy|+|c|\). There the \((|xy|+|c|)\)-relative proposed FMA may lose relative accuracy with respect to the result, whereas the Exact FMA is immune by definition. Measurements confirm the advantage of the Exact FMA under cancellation (errors \(1/4\)–\(1/8\) of BF in the division benchmark of Section 4.1). Its cost, however, is \(194/432/968\) flops, i.e.\(6.6\)–\(11\) times the \(17/66/146\) of the proposed FMA, and measured speeds stay at \(0.07\)–\(0.35\) of BF. It is useful only in local, accuracy-critical places.

2.8 Theoretical operation counts↩︎

Table 2 lists the operation counts per accumulated MAC \(z=a\cdot b+c\), counting a fused FMA instruction as \(1\), with two_sum\(=6\), quick_two_sum\(=3\), two_prod\(=2\), three_sum\(=18\), three_sum2\(=13\), and plain add/sub/mul/fma\(=1\).

3pt

The proposed FMA has the smallest operation count among the fused variants (DD 17 / TD 66 / QD 146 \(<\) BF 29 / 96 / 209), with FLOP ratios FMA/BF of \(1.71/1.45/1.43\). The measured GEMM speedups (\(1.50/1.32/1.46\) on AVX2, Table 11) follow this trend well (the largest speedup for DD, whose ratio is highest). The Exact FMA costs about \(11/6/6\) times the proposed FMA, matching its measured slowdown of 4–5\(\times\).

2.9 Machine verification of the error bounds↩︎

The error bounds of the proposed FMA and the validity of the renormalization gates were machine-verified with Zhang–Aiken’s FPANVerifier [2], a floating-point abstract interpreter driven by the SMT solver z3. The verification scripts (.fpan) treat the precision \(p\) symbolically, so each proof holds simultaneously for all formats (binary32/64/128, …). Every plain \(+\) and \(\times\) of the reference C implementation fma_ref.c (Appendix 7) corresponds one-to-one to a discarded wire (a single rounding error or a truncated term) in the .fpan file; moreover fma_ref.c \(\equiv\) the scalar reference \(\equiv\) all SIMD/GPU implementations agree bitwise, so the proofs carry over to the implementation code as they are.

2.9.1 Certified error bounds↩︎

The best provable bound of each discarded wire is found by binary search, and their sum gives the overall bound (Table 3).

2.9.2 Validity of FastTwoSum and minimization of the renormalization gates↩︎

FPANVerifier’s validity condition for FastTwoSum is \(\mathrm{exp}(a)\ge\mathrm{exp}(b)\) (or one operand equal to \(0\)), not the textbook \(|a|\ge|b|\). Table 4 shows the per-gate test results. Maximizing the provable FastTwoSum gates yields the minimal flops of the proven variant for the given accumulation part: DW turns its trailing gate into FastTwoSum (\(20\to17\)), QW turns two unprovable FastTwoSum into TwoSum (\(143\to146\)), and TW stays at \(66\) after swapping gate roles. The candidates that would lower TW to \(63\) or QW to \(143\) (other FastTwoSum conversions) all fail in the SMT check (Section 5).

3pt

2.9.3 Bitwise agreement between the proven and naive variants (accuracy unchanged)↩︎

The proven variants (\(17/66/146\)) agree bitwise in IEEE binary32/64 with the naive variants that simply place FastTwoSum at the tail, because \(\mathrm{\small FastTwoSum}\equiv\mathrm{\small TwoSum}\) at every proven gate. We confirmed \(0\) mismatches over \(5\times10^6\) trials per type including deep cancellation (ab_provable.c). Hence all accuracy and cancellation tables in this paper are common to both variants; migrating changes only flops and wall-clock time.

2.9.4 Falsification test (MPFR)↩︎

Since the certified bounds are upper bounds, we independently check that the implementations never exceed them. Each \(x\cdot y+c\) is evaluated exactly with MPFR and the maximum of \(\eta=|z-(xy+c)|/(|xy|+|c|)\) is reported in units of \(u^K\). Inputs are non-overlapping expansions, the exponent of \(c\) is swept over about 125 binades relative to \(|xy|\), and a family of deep cancellations \(c\approx-xy\) is mixed in. Not only the scalar but also the AVX2/AVX-512 implementations are exercised directly. As Table 5 shows, there are no violations in 400,000 trials per type per backend. The observed constants (\(\sim5\)) are smaller than the certified ones (\(34/184/812\)) because the proofs are loosened by the “per-wire worst case plus max-anchoring” construction. A negative control that drops one cross term \(x_1y_0\) from DW yields an observed \(C_2=8.7\times10^{15}\) with 17,906 violations, confirming that the test does detect breakage.

2.9.5 A remark on the non-overlap of the QW output↩︎

Because the QW renormalization is a single pass, the leading pair \((z_0,z_1)\) may fail to be a non-overlapping expansion under extreme cancellation (\(c\approx-xy\)); even the weak s_dominates claim \(z_0\!\vartriangleright\!z_1\) is REFUTED by the machine check. The value nevertheless stays within the certified bound, and no accuracy degradation is observed in the recursive accumulation of GEMM/dot products (\(\mathrm{acc}\leftarrow\mathrm{fma}(a,b,\mathrm{acc})\)): the worst relative error equals that of renormalizing after every step. If the output is passed to a consumer that requires a non-overlapping expansion, one extra renormalization pass is needed.

2.10 Commutativity in \(x\leftrightarrow y\)↩︎

The existing branch-free implementation (Zhang–Aiken [1]) is bitwise invariant under operand exchange for both mul and add. Since the fused FMA \(z=x\cdot y+c\) is offered as a replacement for mul\(+\)add, the bitwise agreement of \(\mathrm{fma}(x,y,c)\) and \(\mathrm{fma}(y,x,c)\)—commutativity in \(x\) and \(y\)—is a property that should be preserved.

2.10.0.1 What breaks commutativity.

Exchanging \(x\leftrightarrow y\) transposes the indices of the product expansion, \(P_{ij}\to P_{ji}\) and \(E_{ij}\to E_{ji}\). two_prod and two_sum are EFTs whose outputs (sum and error) are determined by the values of the inputs alone, so they are commutative by themselves (\(\mathrm{fl}(a+b)=\mathrm{fl}(b+a)\); the error terms are exact and hence identical). What breaks commutativity is the association order of the plain additions, floating-point addition being commutative but non-associative. The first, naively written netlist had three associations that were not invariant under transposition: \[\begin{align} \text{TW:}\quad & \sigma = (P_{02} + P_{11}) + P_{20} \;\xrightarrow{\;x\leftrightarrow y\;}\;(P_{20} + P_{11}) + P_{02} \quad(\ne)\\ \text{QW:}\quad & D = ((P_{03} + P_{12}) + P_{21}) + P_{30}\\ & \text{\mathrm{\small TwoSum}{}}(A_2,E_{01})\;\text{immediately followed by}\;\text{\mathrm{\small TwoSum}{}}(A_2,E_{10})\\ & t_1 = E_{02}+E_{11},\quad t_2 = E_{20}+D \end{align}\] DW is automatically commutative because its cross terms involve only one addition, \(l = P_{01} + P_{10}\). Diagonal sums with three or more terms, where association matters, first appear at \(K\ge3\).

2.10.0.2 Symmetrization (operation count unchanged).

Re-associating so that the partners exchanged by transposition (\(P_{ij}\) with \(P_{ji}\), \(E_{ij}\) with \(E_{ji}\)) are paired first makes the whole netlist invariant under transposition, without adding a single operation: \[\begin{align} \text{TW:}\quad & \sigma = (P_{02} + P_{20}) + P_{11} && \text{(still 2 adds)}\\ \text{QW:}\quad & D = (P_{03} + P_{30}) + (P_{12} + P_{21}) && \text{(still 3 adds)}\\ & (\tilde{E},\,g_4) = \text{\mathrm{\small TwoSum}{}}(E_{01},E_{10}) && \text{(still two \mathrm{\small TwoSum}{}}\\ & (A_2,\,g_3) = \text{\mathrm{\small TwoSum}{}}(A_2,\tilde{E}) && \text{\;\;=12 flops)}\\ & t_1 = E_{02} + E_{20},\quad t_2 = E_{11} + D && \text{(still 2 adds)} \end{align}\] Thus the proposed FMA becomes commutative at unchanged DW \(17\) / TW \(66\) / QW \(146\) flops (likewise \(20/69/149\) for the div/sqrt-safe variant).

2.10.0.3 Verification.

After symmetrization, \(\mathrm{fma}(x,y,c)\) and \(\mathrm{fma}(y,x,c)\) agree bitwise for all six types (DD/TD/QD/DS/TS/QS) (\(200{,}000\) trials each, \(0\) mismatches; operands are the normalized non-overlapping components of Section 2.9). Before symmetrization only DW/DS were commutative; TW/TD disagreed in the trailing component in \(8.3\%\) of trials and QW/QD in \(14.7\%\) (the differences stayed inside the certified bounds, with maximum \(|z(x,y)-z(y,x)|/(|xy|+|c|)\) of TW \(5.4\times10^{-48}\) and QW \(6.7\times10^{-64}\)). Harmless for accuracy, but commutativity as bit reproducibility was lost.

2.10.0.4 Impact on the machine verification.

Since the netlist changes, all wire bounds and FastTwoSum gates were re-proved with FPANVerifier. The results improve: for TW individual constants merely permute and the total is unchanged (\(184u^3\)); for QW the bound shrinks from \(48u^4|xy| + 776u^4\max\) to \(44u^4|xy| + 768u^4\max\), i.e.the total from \(824u^4\) to \(812u^4\). No FastTwoSum gate reports an ERROR, and the REFUTED prove statements are the same as before symmetrization (the known issue on the QW output non-overlap discussed in Section 2.9). The bitwise agreement of the five backends (\(200{,}000\) trials each, \(0\) mismatches) was reconfirmed after symmetrization.

2.10.0.5 All numbers re-measured with the symmetric variant.

Following the symmetrization, all measurements in this paper—AVX2 / AVX-512 (Xeon Gold 6526Y) / Arm-NEON / Arm-SVE2 / single-precision base / division and square root / H100 and GB10 (CUDA)—were re-taken with the symmetric variant (the tables from here on). Every backend still agrees bitwise with the scalar reference (\(200{,}000\) trials each for AVX2/AVX-512, \(0\) mismatches for both FMA and Exact; device\(=\)host for H100 GEMM at \(200{,}000\) trials and div/sqrt at \(4096\) points). As symmetrization does not change the operation count (DW \(17\)/TW \(66\)/QW \(146\)), the speed ratios are essentially unchanged; indeed the GEMM ratios versus BF moved by at most \(\pm0.03\) (e.g.AVX2 QD GEMM \(\times 1.46\) \(\to\) \(\times 1.47\), NEON QD GEMM \(\times 1.44\) \(\to\) \(\times 1.47\); the three double-precision types on H100 and GB10 match exactly). The mean relative errors of TW/QW move slightly because the association order of the trailing component changes (e.g.the QW mean \(6.07\times10^{-65}\to6.16\times10^{-65}\)), while the maximum relative errors are unchanged.

2.10.0.6 Symmetrizing the Exact FMA.

The Exact FMA (distillation variant) is not commutative as-is either, because the sweep order is not invariant under transposition (before symmetrization we observed \(2\) mismatches in \(200{,}000\) trials for QS). Here the fix is to combine, before the sweep, each slot pair exchanged by transposition—\((P_{ij},P_{ji})\) and \((E_{ij},E_{ji})\)—with TwoSum. TwoSum is an EFT, so its outputs depend only on the input values; the input expansion becomes transposition-invariant, and the fixed-order sweep and the final truncation are then commutative as well. Each pair costs one TwoSum (\(6\) flops); DW has \(2\) pairs, TW \(6\), and QW \(12\), so the Exact FMA grows from \(182/396/896\) to \(194/432/968\) flops (the proposed FMA is unchanged). With this, both the proposed and the Exact FMA are commutative for all six types (\(200{,}000\) trials each, \(0\) mismatches).

3 Related work↩︎

Several lines of research exploit the FMA for multiple-precision and multi-component arithmetic, but what is called an “FMA” differs from work to work. This section takes up three representative threads—Boldo–Muller’s error of the scalar FMA, Ozaki–Koizumi’s DW-FMA for matrix multiplication, and Muller et al.’s FMA emulation—and positions the proposed FMA against them (Table 6).

2pt

3.1 Boldo–Muller: the exact error of a scalar FMA (ErrFma)↩︎

Boldo and Muller showed that the error of a single hardware FMA is always representable exactly by two floating-point numbers, and gave ErrFma (\(20\) operations), which computes \(r_1=\circ(ax+y)\) together with \(r_2,r_3\) such that \(ax+y=r_1+r_2+r_3\) exactly [3]. They further provided a Coq formalization, extended beyond radix 2 to general even radices including subnormals. An approximate variant, ErrFmaAppr, trades a little accuracy for only \(12\) operations.

The difference from our work lies in the level of the computation. ErrFma concerns the error of one FMA applied to scalars \(a,x,y\in\mathbb{F}\), not an FMA of multi-component numbers. In our proposed FMA, both inputs and output are \(K\)-component (DW/TW/QW), and the internal two_prod/two_sum use the hardware FMA. The two lines of work do not compete; they are complementary.

One point of contact is particularly important. The “fma-folded diagonal” listed as future work in Section 5 (folding \(a\cdot b+c\) with the single rounding of a fused operation would further reduce the counts to DW \(16\)/TW \(64\)/QW \(140\) flops) was not adopted in this paper precisely because the single rounding of a fused operation cannot be expressed as a composition of two_prod\(+\)two_sum, hence cannot be machine-verified with FPANVerifier. Boldo–Muller’s theorem supplies exactly this missing gate—the exact representation of the residual of a fused \(ax+y\)—and a verifier extended with an “FMA gate” based on their Coq theorem would bring the folded variants into the reach of machine verification. This is the most concrete bridge between the two works.

A caveat on operation counts: the \(20\) operations of ErrFma pay for the error of one scalar FMA, whereas the \(17\) flops of our DW-FMA pay for the whole \(K{=}2\) MAC—the two numbers measure different things. Boldo–Muller themselves state that no cheaper algorithm was found, and, like our claim that \(17/66/146\) is minimal within the given renormalization structure, both remain claims of “the current minimum.”

3.2 Ozaki–Koizumi: a DW-FMA for matrix multiplication (FastFma_DW)↩︎

Ozaki and Koizumi gave the first error analysis of FastTwoFma, an approximate EFT for \(ab+c\) built from just two FMAs (\(x=\texttt{fma}(a,b,c)\), \(y=\texttt{fma}(a,b,c-x)\); Koizumi et al. [6]), and proposed the DW-FMA FastFma_DW (only \(6\) operations) together with the matrix multiplication Sesqui_Word_MM built on it [4]. On a Ryzen 9 7950X, it stays within \(\times 4.1\)–\(\times 4.7\) of OpenBLAS dgemm, faster than pair arithmetic (PA, \(\times 6.5\)–\(\times 7.4\)) and DW (\(\times 11\)–\(\times 13\)).

The decisive point is the applicability condition. The error of FastTwoFma stays at \(O(u^2)\) only when \(c-x\) is exactly representable; otherwise an \(O(u)|ab|\) term remains. In particular, with \(\mathrm{sign}(c)\ne\mathrm{sign}(ab)\) and \(|c|\approx|ab|\) (strong cancellation) the relative error can degrade badly. They prove sufficient conditions such as \(|c|\ge2|ab|\) and, for matrix multiplication, guarantee the condition at every \(k\) by adding a shift matrix \(E\) (whose entries \(e_{ij}\) are built from row/column norms with upward rounding).

Three differences from our work.

(1) Presence of preconditions. Their DW-FMA is a fused MAC with conditions on the operands, completed as a matrix multiplication algorithm together with the machinery (the precomputed shift matrix \(E\)) that enforces the condition. Our proposed FMA is unconditional for arbitrary operands, with every FastTwoSum precondition machine-proved for all precisions \(p\) (Section 2.9). It therefore drops into AXPY/GEMV/GEMM—and into division and square root—as a drop-in replacement of mul\(+\)add. This difference is not academic: as Section 4.1 shows, the inner MACs of division and square root receive non-non-overlapping residuals, and even in our work part of the gate placement proved for the standard MAC (the unproved FastTwoSum) had to be demoted to TwoSum (the div/sqrt-safe variant, DW \(20\)/TW \(69\)/QW \(149\)). Our finding that “a proven gate placement must be re-verified for each usage context” explains precisely why Ozaki–Koizumi engineered the condition with \(E\) specifically for matrix multiplication.

(2) Accuracy target. They deliberately choose an accuracy below DW (\(1.5\) words, “sesqui-word”) to buy speed. Indeed, at \(n{=}1000\) the maximum relative error of their method is \(5.5\times 10^{-22}\) versus PA \(2.4\times 10^{-24}\) and DW \(9.2\times 10^{-25}\), i.e.–3 digits worse [4]. Our proposed FMA requires the working accuracy of DW/TW/QW to be preserved: the bounds \(34u^2/184u^3/812u^4\) are machine-proved and the roundings agree bitwise across five backends.

(3) Component counts. They treat only \(K{=}2\) (DW). We treat \(K{=}2,3,4\) (DW/TW/QW) uniformly, and the benefit of removing renormalization stages grows with \(K\) (GEMM vs BF: QD \(\times 1.44\)–\(\times 1.61\); division: QD \(\times 2.03\)–\(\times 2.11\)).

3.2.0.1 Theoretical complement: tight bounds under parameterized overlap (Jeannerod–Joldeş–Louvet–Muller).

The condition required by Ozaki–Koizumi’s FastFma_DW—\(|c|\ge2|ab|\) guaranteeing the exactness of \(c-x\)—is no accidental trick but a structure common to fast FMA kernels of this kind. Jeannerod, Joldeş, Louvet, and Muller recently addressed it head-on, analyzing the error bounds of a family of fast FMA kernels that approximate \(ab+c\) in DW representation, under overlap-parameterized conditions \(|x_\ell|\le k\,\mathrm{ulp}(x_h)\) (including two- and three-parameter versions for multiple DW inputs) [5]. For the same dominance regime \(|c|\ge2|ab|\) as Ozaki–Koizumi and its DW forms (\(|c_h|\ge2|ab|\), \(|c_h|\ge2|ab_h|\), etc.; arising in polynomial evaluation of elementary functions and other cancellation-free settings), they exhibit explicit worst-case constants, together with examples attaining them, tightening the previous analyses (ARITH 2026).

The relation to our work is twofold.

(1) Confirmation that the guarantee is conditional. Their results quantitatively confirm the point of this section: fast FMAs of the FastFma_DW kind carry guarantees conditioned on dominance and overlap. The significance of our FMA being operand-unconditional (Section 2.9) stands out in this contrast.

(2) Implications for our own assumptions. Conversely, our proofs are not free of overlap assumptions either: the .fpan scripts assume strongly_dominates (non-overlap) on the input components, which is their parameter \(k\) taken at its strongest. Indeed, under extreme cancellation (\(c\approx-xy\)) the QW output may fail to be non-overlapping (Section 2.9), and feeding it back as input in a recursive dot-product accumulation violates the assume. We verified by measurement that the values stay within the certified bounds, but proving bounds under weakened \(k\) requires their parameterization. This is the most natural next step of our work, corresponding to relaxing the assume of FPANVerifier by the parameter \(k\) and re-proving.

In summary, the two approaches occupy different points on the accuracy–speed trade-off curve. If \(1.5\)-word accuracy suffices and the computation is confined to matrix multiplication, the \(6\)-operation FastFma_DW wins; if a general-purpose FMA that unconditionally preserves full DW/TW/QW accuracy is needed, our \(17/66/146\) is the answer. They are not exclusive: a natural combination is “FastFma_DW for GEMMs where the \(E\) shift can be applied, and the proposed FMA for division, square root, and wherever the condition cannot be manufactured.”

3.3 Muller et al.: FMA emulation and the building blocks of DW/TW↩︎

Graillat and Muller gave algorithms that emulate the FMA and the correctly rounded sum of three numbers (ADD3) using only round-to-nearest addition, multiplication, and comparisons, for platforms without a hardware FMA [7] (with earlier work by Boldo–Melquiond using rounding to odd [8]). This is exactly the inverse problem setting of ours: all five of our evaluation targets (AVX2 / AVX-512 / NEON / SVE2 / CUDA) have hardware FMA, and we actively exploit the fact that \(\texttt{two\_prod}(a,b)=(\circ(ab),\;\texttt{fma}(a,b,-\circ(ab)))\) costs 2 operations. Without hardware FMA, two_prod degrades to Dekker splitting (\(17\) operations), and the flop table of this paper (Table 2) itself ceases to be a valid cost model. Conversely, such platforms would need the Graillat–Muller emulation, whose cost easily exceeds the savings of the proposed FMA (\(29\to17\) flops, etc.).

The same school includes Joldeş–Muller–Popescu, who gave tight error bounds for the building blocks of DW arithmetic [9], its Coq formalization by Muller–Rideau [10], and triple-word algorithms by Fabiano–Muller–Picot [11]. These works built the very arena of DW/TW/QW on which we stand, but all of them target mul and add; none defines a fused FMA as a single operation. The contribution of this paper is to place a branch-free, machine-verified multi-component FMA in that gap and to fold the two stages of the existing mul\(+\)add into one, with real speedups.

4 Benchmark tests↩︎

We evaluate the proposed FMA on two platforms. The Arm platform is NVIDIA GB10 (Grace–Blackwell): the CPU is an Arm Cortex-X925 (3.9 GHz, NEON / SVE2 with VL\(=\)​128 bit, 20 big.LITTLE cores) and the GPU is Blackwell (sm_121, 48 SM, LPDDR5X). The x86 platform is an H100 node: the CPU is an Intel Xeon Gold 6526Y (AVX2 / AVX-512, 32 threads) and the GPU is an NVIDIA H100 NVL (sm_90, 132 SM, HBM3). Compilers are gcc 13.3 and nvcc; the accuracy reference is MPFR at 600 bits (DD\(\approx\)​106 / TD\(\approx\)​159 / QD\(\approx\)​212 bit). Section 4.1 evaluates division and square root, and Section 4.2 the basic linear-algebra kernels (AXPY/GEMV/GEMM; serial, OpenMP-parallel, and GPU-parallel).

4.0.0.1 Implementation.

The scalar two_sum / two_prod / fast_two_sum are mapped one-to-one to the BNCmatmul AVX2 primitives _bncavx2_dtwo_sum / _bncavx2_dtwo_prod / _bncavx2_dquick_two_sum (which use the hardware FMA), and plain \(+,\times\) become _mm256_add/mul_pd. Data are stored column-major in SoA form with one double* per component, so one instruction processes 4 output elements in parallel. Table 7 lists the four accumulation variants \(z \leftarrow a\cdot b + c\) compared throughout.

2pt

Measurement environment for x86: Intel (AVX2+FMA+AVX-512, 32 threads), gcc 13.3, -O3 -mavx2 -mfma -funroll-loops; accuracy reference MPFR 600-bit. DD\(\approx\)​106bit / TD\(\approx\)​159bit / QD\(\approx\)​212bit.

4.0.0.2 Verifying the translation.

For the single operation \(x\cdot y+c\), every backend implementation agrees bitwise with the scalar reference (DD/TD/QD \(\times\) {proposed, Exact}, 200,000 trials each, 0 mismatches). Vectorization does not change the rounding.

4.1 Accelerating division and square root↩︎

4.1.1 Motivation and algorithms↩︎

The inner loops of multi-component division and square root consist entirely of MACs of the form \(x\cdot y+c\). The proposed FMA should therefore improve not only stand-alone MACs and AXPY/GEMV/GEMM but also division and square root. This section verifies this for all six types DS/TS/QS/DD/TD/QD.

4.1.1.1 Division (long division).

Take \(K+1\) quotient digits, correcting the residual per digit: \[q_i=\frac{r_0}{b_0},\qquad r \leftarrow r - q_i\cdot b \quad (i=0,\dots,K-1),\qquad q_K=\frac{r_0}{b_0},\] finally distilling the \(K+1\) values \(q_0,\dots,q_K\) into \(K\) components (\(K\) bubble sweeps). The correction \(r\leftarrow r-q_i b\) is exactly a MAC; note that the multiplier \(q_i\) is a scalar.

4.1.1.2 Square root (Newton iteration for \(1/\sqrt{a}\)).

\[r_0=\frac{1}{\sqrt{a_0}},\qquad h=\frac{a}{2},\qquad e \leftarrow \tfrac12 - h\cdot r^2,\qquad r \leftarrow r + r\cdot e,\qquad \sqrt{a}\simeq a\cdot r .\] Both updates of \(e\) and \(r\) are MACs. Since the iteration starts at about \(p-1\) bits and doubles per step, we use 2 iterations for \(K=2,3\) and 3 for \(K=4\) (2 iterations at \(K=4\) would give only \(208<212\) bits for QD and \(92<96\) bits for QS). The iteration count is common to the three methods and does not affect the speed ratios.

4.1.1.3 Three methods.

The MACs above are implemented in three ways.

• BF (baseline): two stages of branch-free mul\(+\)add (Zhang–Aiken); the product is rounded back to \(K\) components before the addition.

• FMA (proposed): the branch-free fused FMA with a single distillation.

• Exact FMA: the result-relative distillation FMA.

The square \(r\cdot r\) and the final \(a\cdot r\) use the same mul in all three methods (a real library always has the cheaper mul). Only the MAC part differs, and that difference is exactly the contribution of the proposed FMA.

4.1.1.4 Zero-folded variants mul_d / fma_d.

In the correction stage of division the multiplier \(q_i\) is a scalar, so we use dedicated variants with \(y_1=\dots=y_{K-1}=0\) folded in. Since \(\mathrm{two\_prod}(x,0)=(0,0)\) and \(\mathrm{two\_sum}(a,0)=(a,0)\) are both exact, the folding does not change any rounding: indeed fma_d\((x,q,c)\) agrees bitwise with fma\((x,\{q,0,\dots\},c)\) (6 types \(\times\) 20,000 trials, mismatch\(=0\)). The BF side likewise uses mul_d, so the comparison is fair, and this matches the construction that BNCmatmul’s _bncneon_rqd_div etc.actually use (mul_d\(+\)sub).

4.1.2 Implementation↩︎

Division and square root are written as backend-independent C++ templates in a single divsqrt.hpp, expanded from the same source into four backends: Scalar / NEON (float64x2_t, float32x4_t) / SVE2 (fixed-length svfloat64_t, svfloat32_t) / CUDA (scalar per thread). Everything reduces to IEEE-754 basic operations and true hardware FMA, so the four agree bitwise by construction. We confirmed by measurement that each NEON / SVE2 / CUDA(device) kernel agrees bitwise with the scalar run of the same source (\(4096\)–\(8192\) points each, mismatch\(=0\)); the NEON and SVE2 relative-error tables agree in all 36 rows. The host requires -ffp-contract=off and nvcc requires --fmad=false (removing them breaks two_sum/two_prod).

4.1.2.1 Context dependence of the gate placement, and the machine-verified div/sqrt-safe variant.

Gate changes go in two directions. Demoting an unproved FastTwoSum to TwoSum (safe) never loses accuracy, since TwoSum is unconditionally exact, and makes the proof possible (\(+3\) flops/gate); promoting TwoSum to FastTwoSum (dangerous) is exact only in contexts where the precondition holds. The standard MAC can prove the FastTwoSum preconditions at \((B,m_1)\)/\((B,A_1)\), but the FMAs inside div/sqrt receive intermediate values that are not non-overlapping (iteration residuals), where the precondition fails and this optimization degrades accuracy. For div/sqrt we therefore use a “safe” variant that merely demotes the unproved FastTwoSum to TwoSum (DW \(20\) / TW \(69\) / QW \(149\) flops). It agrees bitwise with the original in relative error (all 12 H100 cases) and is machine-verifiable (DW: dwfma.fpan, TW: twfma_fix.fpan, QW: qwfma_safe.fpan; FastTwoSum ERROR \(0\) in all). Since fma_d (\(y=\{y,0,\dots\}\)) agrees bitwise with fma\((x,\{y,0,\dots\},c)\) (\(3\times10^6\) trials/type, \(0\) mismatches), it inherits the proof of the full fma (\(y_1{=}y_2{=}0\) satisfies the dominance assumptions). A machine-proved gate placement must be re-verified for each usage context; this is why the minimal proven flops differ between the standard MAC (\(17/66/146\)) and div/sqrt (\(20/69/149\)).

4.1.3 Dynamic operation counts↩︎

Table 8 shows the number of operations actually executed (add/sub/mul/neg/div/sqrt/fma counted as 1 each), obtained at run time with a counting backend—the numbers of the very code the benchmarks run. They depend only on the component count \(K\), not on the component type (float/double).

3pt

4.1.4 Speed↩︎

Table 9 shows the speed ratios of the proposed FMA against the BF baseline. CPU runs use 1 core (taskset) with \(n=131072\) elements; GPU runs use \(n=4{,}194{,}304\) elements. The GPUs are NVIDIA GB10 (sm_121, 48 SM, LPDDR5X) and NVIDIA H100 NVL (sm_90, 132 SM, HBM3), evaluated with the same binary and arguments.

3pt

4.1.4.1 Why the CPU exceeds the operation-count ratio.

QD division measures \(\times 2.03\)–\(\times 2.11\), well above the operation-count ratio of \(1.51\). In the BF path, two serial chains of sequential quick_two_sum—the normalization of mul_d and that of add—are folded into one by the proposed FMA. Besides the reduction in operations, the dependence chain shortens and ILP improves; the effect grows with \(K\).

4.1.4.2 H100: all six types drift to the bandwidth-bound side.

Running the same source and arguments on the H100 NVL (sm_90, 132 SM, HBM3), the elementwise kernels show FMA/BF of \(\times 1.02\)–\(\times 1.37\) for division and \(\times 1.11\)–\(\times 1.25\) for square root. Unlike GB10, even the three double-precision types do not become compute-bound, because the H100’s FP64 is strong (half of FP32) and memory saturates first. Effective bandwidths for DIV-FMA are DD \(3146\) / TD \(2600\) / QD \(2116\) GB/s and DS \(2767\) / TS \(2414\) / QS \(2005\) GB/s, i.e.\(51\)–\(81\%\) of the HBM3 peak (about \(3900\) GB/s). The most saturated case, DD division (\(81\%\)), gives FMA/BF \(=1.02\); the least saturated, QD division (\(54\%\)), gives \(1.35\)—saturation and speedup anti-correlate cleanly. The compute-only chains give division \(\times 1.25\)–\(\times 1.56\) and square root \(\times 1.22\)–\(\times 1.28\), matching GB10 and the theoretical operation-count ratios. The operation-count advantage is thus invariant across GPU generations and memory technologies; only the degree to which it appears in wall-clock time is set by the compute/bandwidth balance of the kernel. In absolute terms the H100 is 9–53 times faster than GB10 (e.g.QD division FMA \(10.11 \to 0.190\) ms).

4.1.4.3 Why GB10’s three single-precision types sit at \(\times 1.00\).

The elementwise kernels merely read and write \(3K\) component arrays, so the arithmetic intensity is low. The measured effective bandwidths are DS \(243\) / TS \(238\) / QS \(238\) GB/s, pinned at the LPDDR5X peak of GB10 (about \(273\) GB/s): bandwidth-bound. The three double-precision types become compute-bound because GB10’s FP64 is weak (DD \(130\) / TD \(61\) / QD \(40\) GB/s), so the operation-count reduction shows up directly in wall-clock time—the same picture as the GEMMs of Section 4.2.8. Extracting pure arithmetic with register-resident compute-only chains (chdiv / chsqrt) that never touch global memory inside the loop, all six types match the theoretical operation-count ratios (Table 9): the operation-count advantage of division and square root is real on GPUs as well.

4.1.4.4 Exact FMA.

For division, Exact/BF is NEON \(\times 0.28\)–\(\times 0.35\), SVE2 \(\times 0.13\)–\(\times 0.27\), GB10 \(\times 0.20\)–\(\times 0.34\), H100 \(\times 0.22\)–\(\times 0.57\) (the symmetric Exact is slightly slower than before symmetrization due to the added operations). The large drop on SVE2 occurs because the 36 fixed-length vectors (QD) demanded by the distillation sweep do not fit in registers and spill—the same phenomenon as the MAC in Section 4.2.6. The proposed FMA and BF carry no expansion array and are unaffected.

4.1.5 Relative error↩︎

Table 10 shows the relative errors against MPFR 600-bit (\(n=8192\); max / mean). \(\varepsilon_K=2^{-Kp}\) is the rounding unit of each type (DS \(3.55\times 10^{-15}\), TS \(2.12\times 10^{-22}\), QS \(1.26\times 10^{-29}\), DD \(1.23\times 10^{-32}\), TD \(1.37\times 10^{-48}\), QD \(1.52\times 10^{-64}\)).

3pt

4.1.5.1 Division.

In units of \(\varepsilon_K\), the maximum errors are \(1.3\)–\(2.3\) for BF, \(1.5\)–\(4.9\) for the proposed FMA, and \(0.2\)–\(0.5\) for the Exact FMA. The proposed FMA is slightly worse than BF (\(1\)–\(2\times\), up to \(3\times\) for TD): BF renormalizes once when rounding the product back to \(K\) components in the correction stage, whereas the proposed FMA discards more of the low-order cross terms. Still, all stay within a few \(\varepsilon_K\) and are practically equivalent. The Exact FMA has \(1/4\)–\(1/8\) of BF’s error (closest to correct rounding) at \(1/3\)–\(1/8\) of the speed.

4.1.5.2 Square root.

The errors of the three methods coincide almost exactly (for DS and DD the maximum errors are bitwise identical). The Newton iteration converges far below \(\varepsilon_K\), so the final error is determined solely by the rounding of the last multiplication \(a\cdot r\), common to all three (\(2.5\)–\(3.9\) times \(\varepsilon_K\)). For square root the proposed FMA is therefore \(\times 1.51\)–\(\times 1.61\) faster without losing a single bit of accuracy—a pure gain. For the same reason there is little point in building square root on the Exact FMA (the final multiplication dominates the error).

4.1.6 Summary↩︎

The proposed FMA accelerates not only the MAC-style kernels (AXPY/GEMV/GEMM) but also division and square root. Division is \(\times 1.54\)–\(\times 2.11\) faster on CPUs and \(\times 1.25\)–\(\times 1.61\) on GPUs (when compute-bound) with accuracy on par with BF; square root is \(\times 1.51\)–\(\times 1.61\) faster with accuracy essentially identical to BF. The gain is largest for large \(K\) (QD/QS), since the number of renormalization stages the proposed FMA removes grows with \(K\).

4.2 Basic linear algebra: AXPY, GEMV, GEMM↩︎

We evaluate accuracy and speed of the three kernels AXPY (\(y\leftarrow\alpha x+y\), \(n=10^6\)), GEMV (\(n=2048\)), and GEMM (\(n=512\)).

4.2.1 Accuracy (relative error against MPFR)↩︎

Table ¿tbl:tab:acc? shows the maximum and mean relative errors. All four variants maintain the working accuracy (DD\(\sim\)​31, TD\(\sim\)​46, QD\(\sim\)​61 digits). The proposed FMA is on par with the existing variants (Q/BF), with the worst case at most about one digit worse (due to the \((|xy|+|c|)\)-relative semantics). The Exact FMA, being result-relative, is the most accurate.

4pt

4.2.1.1 Cancellation (ill-conditioned dot products).

Even at condition number \(2\times10^{13}\) the proposed FMA degrades numerically almost identically to the existing variants (DD: Q \(1.3\times 10^{-20}\) \(=\) FMA \(1.3\times 10^{-20}\)). Replacing with the proposed FMA causes no loss of accuracy.

4.2.2 Serial performance (x86: AVX2)↩︎

Table 11 shows the measured time per call of the proven variant (\(17/66/146\) flops) on AVX2 (Xeon Gold 6526Y, best of 3). The fair baseline is BF.

3pt

Branch freedom (BF) already accounts for most of the existing speedup (Q\(\to\)BF gives e.g.\(\times 4.0\) for TD AXPY); the proposed FMA adds on top of it. In the compute-bound GEMM it reaches DD \(\times 1.49\) / TD \(\times 1.32\) / QD \(\times 1.47\), with DD fastest, reflecting the flops of the proven variant (\(17/66/146\)). Compared with the naive variant (\(20/143\)), DW gains from the FastTwoSum-ized tail (GEMM \(\times 1.32\) \(\to\) \(\times 1.49\)) and QW loses slightly from the two FastTwoSum demoted to TwoSum (\(\times 1.50\) \(\to\) \(\times 1.47\)); TW is unchanged. In the bandwidth-bound AXPY/GEMV the differences are small, and QD GEMV ties BF (\(\times 1.03\)).

4.2.3 Serial performance (x86: AVX-512)↩︎

On the same Xeon Gold 6526Y we prepared an AVX-512 implementation, doubling the SIMD width from 4 to 8. The translation is mechanical: __m256d\(\to\)__m512d, _mm256_*\(\to\)_mm512_*, EFT primitives _bncavx2_*\(\to\)_bncavx512_*, loop stride i+=8, and alignment 32\(\to\)​64 bytes; the algorithms are otherwise identical. Compilation: gcc -O3 -mavx2 -mfma -mavx512f -mavx512dq -mavx512vl -mavx512cd -funroll-loops.

4.2.3.1 Verifying the translation.

For the single operation \(x\cdot y+c\), the AVX-512 implementation also agrees bitwise with the scalar reference (DD/TD/QD \(\times\) {proposed, Exact}, 200,000 trials each, 0 mismatches). Moreover all 18 max/mean relative-error values of BF, proposed, and Exact, and all 9 rows of the cancellation test, coincide exactly with AVX2 (Table ¿tbl:tab:acc?). Rounding is invariant under the SIMD width.

3pt

The proposed FMA achieves FMA/BF \(\times 1.03\)–\(\times 1.77\) in every case on AVX-512 (Table 12), most prominently QD AXPY (\(\times 1.77\)) and GEMM (DD \(\times 1.39\) / TD \(\times 1.43\) / QD \(\times 1.48\)). QD GEMV, which tied on AVX2, is clearly faster at \(\times 1.39\). The effect of doubling the SIMD width itself approaches the width ratio \(\times 2\) for the compute-bound GEMM (proposed FMA: TD \(\times 1.73\) / QD \(\times 1.68\)). Being branch-free, the proposed FMA scales cleanly with lane width.

4.2.3.2 A pitfall of the default implementation (Q).

The QD FMA/Q of Table 12 (\(\times 21.4\)) means less than it seems. In BNCmatmul, _bncavx512_rqd_{mul,add}_sloppy has its #if 0 inverted relative to the AVX2 version, so the SIMD body is disabled and the path that processes the 8 lanes individually with scalar rqd_mul is active. The AVX-512 QD-Q is thus not vectorized and is 4.5 times slower than the AVX2 QD-Q (GEMM \(9.36\) s vs \(2.08\) s). Also, TD-Q on AVX-512 defaults to _rtd_mult, which the AVX2 header explicitly disables as “fast on Xeon but include bugs (drop accuracy),” making the TD-Q accuracy about one digit worse than AVX2. Both are pre-existing issues of the library, independent of this proposal, but they reinforce that the fair baseline is consistently BF. In addition, _bncavx512_vseb initializes only j[0..3] of the per-lane counters int j[8], and the default TD Q segfaults (an incomplete port from the 4-lane version); our measurements were taken after fixing this.

4.2.3.3 OpenMP.

The FMA/BF ratio is stable with respect to thread count, most prominently in the compute-bound GEMM (QD GEMM: \(\times 1.73\)/\(\times 1.74\)/\(\times 1.60\)/\(\times 1.61\) at 1/4/16/32 threads; TD GEMM: \(\times 1.41\)/\(\times 1.41\)/\(\times 1.37\)/\(\times 1.37\)). The parallel scaling of TD GEMM(\(n{=}512\)) from 1 to 32 threads is BF \(\times 27.8\) / FMA \(\times 27.1\). The accuracy defect of the existing _bncavx512_rdd_fma (\(2.2\times 10^{-13}\) in GEMM) reproduces as on AVX2/NEON.

4.2.4 Serial performance (Arm: NEON)↩︎

Note (this and the next two subsections): the Arm/NEON/SVE2 implementations have migrated to the proven variant (DW 17/TW 66/QW 146), and the timings below are re-measured with the proven variant on a Cortex-X925. Accuracy and bitwise agreement are unchanged (the proven variant agrees bitwise with the naive one, so the relative errors of Table ¿tbl:tab:acc? are unaffected).

The AVX2 implementation above was ported to Arm NEON. The scalar EFT primitives map one-to-one to _bncneon_dtwo_sum / _dtwo_prod / _dquick_two_sum (using the hardware FMA), and plain \(+,\times\) become vaddq_f64 / vmulq_f64. Since NEON’s float64x2_t holds 2 double-precision elements, the SoA row parallelism drops from \(4\) to \(2\) and the loop stride becomes i+=2; the algorithm is otherwise identical to AVX2. Environment: Arm Cortex-X925 (3.9 GHz, NEON/asimd), gcc -O3 -march=native -DBNC_ENABLE_NEON -funroll-loops, accuracy reference MPFR 600-bit.

4.2.4.1 Verifying the translation.

For the single operation \(x\cdot y + c\), the NEON implementation also agrees bitwise with the scalar reference (DD/TD/QD \(\times\) {proposed, Exact}, 200,000 trials each, 0 mismatches). Maximum relative errors and cancellation behavior are numerically identical to x86 (Table ¿tbl:tab:acc?).

3pt

The proposed FMA shows FMA/BF \(\times 1.26\)–\(\times 2.01\) in every case on NEON (Table 13), most prominently AXPY (TD \(\times 2.01\), QD \(\times 1.81\)); GEMM gives DD \(\times 1.40\) / TD \(\times 1.42\) / QD \(\times 1.47\). QD GEMV, the only tie on x86, is clearly faster at \(\times 1.36\). The Exact FMA runs at \(\times 0.20\)–\(\times 0.37\) of BF, much as on x86. With 2 lanes the absolute times are roughly 1.0–1.7 times x86, but the FMA/BF ratios are preserved at x86 level or better.

4.2.4.2 OpenMP.

The FMA/BF ratio is stable in thread count, especially for the compute-bound GEMM (TD GEMM: \(\times 1.60\)/\(\times 1.57\)/\(\times 1.45\)/\(\times 1.42\) at 1/4/8/16 threads; QD GEMM: \(\times 1.43\)/\(\times 1.43\)/\(\times 1.43\)/\(\times 1.42\)). TD GEMM(\(n{=}512\)) scales almost linearly from 1 to 16 threads (BF \(\times 15.8\) / FMA \(\times 14.1\); 20 big.LITTLE cores). The bandwidth-bound AXPY/GEMV converge to ratios near 1 at high thread counts. The accuracy defect of the existing _bncneon_rdd_fma (\(2.2\times 10^{-13}\) in GEMM) reproduces as on x86.

4.2.5 Single-precision base DS/TS/QS (Arm: NEON)↩︎

We extended the proposed FMA to DS/TS/QS (double-single, triple-single, quad-single), whose components are single-precision float. The algorithms have the same structure as DD/TD/QD, with double\(\to\)float, fma\(\to\)fmaf, NEON float64x2_t (2 lanes) \(\to\) float32x4_t (4 lanes), and EFT primitives such as _bncneon_ftwo_sum. The NEON implementation agrees bitwise with the scalar reference (DS/TS/QS \(\times\) , 200,000 trials each, 0 mismatches) and maintains the working accuracy (DS\(\approx\)​14, TS\(\approx\)​21, QS\(\approx\)​28 digits).

4pt

As Table ¿tbl:tab:acc95single? shows, all three variants maintain the working accuracy. The Exact FMA gives the best value in 7 of 9 rows for the mean and 6 of 9 rows for the maximum relative error, improving the DS maximum by nearly 3 digits (AXPY: \(3.1\times 10^{-12}\) \(\to\) \(1.7\times 10^{-15}\)). The rows where it is not best (Q for TS GEMV/GEMM and QS AXPY, BF for QS GEMM) are cases where the maximum error is determined not by the rounding of individual fused operations but by cancellation in the accumulated sum, which the result-relative semantics of the Exact FMA cannot improve. The proposed FMA is on par with the existing variants (Q/BF); the worst-case maximum error can be up to \(\times 7.9\) worse than BF (QS GEMM), which stems from the \((|xy|+|c|)\)-relative semantics as in the double base, while the mean errors stay at BF level (a drop-in replacement without loss of accuracy).

3pt

The proposed FMA is \(\times 1.37\)–\(\times 1.83\) versus BF in every case with the single-precision base as well (Table 14), mirroring the double base. Under OpenMP the GEMM advantage is stable in thread count (TS GEMM \(\times 1.40\)–\(\times 1.63\), QS GEMM \(\times 1.39\)–\(\times 1.41\); TS GEMM scales from 1 to 16 threads at BF \(\times 15.9\) / FMA \(\times 13.6\)). Note that BNCmatmul’s neon/_bncneon_f.h shares the include guard __BNCNEON_DS_H with the full DS version neon/_bncneon_ds.h, so the DS sloppy/bf/fma variants are shadowed when included via bncneon.h (needs fixing); our evaluation includes the full DS version first to work around this.

4.2.6 Serial performance (Arm: SVE2)↩︎

The NEON suite was ported to Arm SVE2 (the SIMD extension of Armv9-A) and re-run on the same machine (Cortex-X925). The VL-agnostic primitives in BNCmatmul’s include/sve2/, _bncsve2_dtwo_sum / _dtwo_prod / _dquick_two_sum (predicated with pg, using the hardware FMA), are mapped one-to-one, with \(+,\times\) as svadd_f64_x / svmul_f64_x. Row loops are predicated with svwhilelt_b64 and need no scalar tail.

4.2.6.1 The sizeless-type constraint and fixed-length SVE types.

SVE’s svfloat64_t is a sizeless type: it cannot be an array element or a struct member and has no sizeof. The Exact FMA performs distillation sweeps over an indexed expansion \(v[0..35]\) and cannot be written with plain svfloat64_t. Our implementation therefore builds arrays with the GCC/Clang fixed-length SVE types __attribute__((arm_sve_vector_bits(N))) (requiring -msve-vector-bits=N) and calls the VL-agnostic library primitives through the implicit conversion between fixed-length and sizeless types. A macro bridges to the library’s per-limb-variable signatures (f(pg,&r0,&r1,a0,a1,b0,b1)). At run time we check that svcntd() matches the build-time lane count.

4.2.6.2 Correctness and accuracy.

The SVE2 implementation also agrees bitwise with the scalar reference (DD/TD/QD \(\times\) , 200,000 trials each, 0 mismatches). Together with AVX2 / AVX-512 / NEON / CUDA, rounding is invariant across all five implementations. Relative errors and cancellation behavior equal x86/NEON (QD shows relative error \(0\) in all tested conditions).

The Cortex-X925 has VL\(=128\) bit (svcntd()\(=2\)), i.e.the same double-precision lane count as NEON, so this subsection is an ISA comparison at equal lane width.

3pt

The proposed FMA achieves FMA/BF \(\times 1.29\)–\(\times 1.88\) in every case (Table 15), mirroring NEON. In particular, QD GEMM reaches its best ratio \(\times 1.62\) on SVE2 (\(\times 1.47\) on NEON), because BF is slower on SVE2 (\(1.45\) s vs NEON \(1.30\) s) while the proposed FMA is equal (\(0.884\) s vs \(0.884\) s)—a consequence of register count (z0–z31) and latency-hiding differences. At equal lane width the absolute times of the proposed FMA are about equal to NEON (\(\pm15\%\)): the advantage of the proposed FMA does not depend on the ISA. Two independent runs reproduced the ratios within \(\pm0.06\).

Under OpenMP the advantage in the compute-bound TD/QD GEMM is invariant in thread count (TD GEMM \(\times 1.47\)–\(\times 1.64\), QD GEMM \(\times 1.48\)–\(\times 1.50\); TD GEMM scales from 1 to 16 threads at BF \(\times 15.8\) / FMA \(\times 14.1\)). The bandwidth-bound GEMV and the low-intensity DD GEMM saturate at BF level with ratios near \(1\) at high thread counts.

4.2.6.3 Stack spills of the Exact FMA.

Only the Exact FMA is much slower on SVE2 than on NEON (QD GEMM: NEON \(6.09\) s \(\to\) SVE2 \(15.21\) s, \(\times 2.5\)). The cause is the stack spill of the distillation expansion \(v[0..35]\): the fixed-length SVE array \(v[36]\) does not fit in the 32 z registers, and the compiler spills with predicated st1d stores. Disassembly counts of vector loads/stores: SVE2 gemm_f_qd (proposed FMA, no array) \(4+4=8\); SVE2 gemm_x_qd (Exact FMA) \(19+50=69\); NEON gemm_x_qd \(6+23=29\). The spill volume is proportional to the expansion length and matches the measured degradation (DD \(\times 1.1\) \(\to\) TD \(\times 2.1\) \(\to\) QD \(\times 2.5\)). The proposed FMA, the subject of this paper, uses no expansion array and is unaffected. Using the Exact FMA on SVE2 would require manually unrolling the expansion into per-limb variables.

4.2.7 OpenMP parallelization↩︎

GEMM is parallelized over columns and GEMV/AXPY over row blocks with #pragma omp parallel for. The FMA/BF ratio of the proposed FMA is essentially invariant in thread count (Table 16); the speedup is orthogonal to parallelization. Parallel scaling itself is good: BF GEMM(\(n{=}512\)) scales about 50\(\times\) at 32 threads. The best effective throughput is obtained by “proposed FMA \(\times\) parallelization.”

3pt

4.2.8 GPU parallel (CUDA: H100 and GB10, all six types)↩︎

We implemented the proposed branch-free FMA and the Exact FMA for all six types in CUDA and evaluated them on an NVIDIA H100 NVL (Hopper, sm_90, 132 SM) and an NVIDIA GB10 (Grace-Blackwell, sm_121, 48 SM). A C++ template unifies the component type \(R\in\{\texttt{float},\texttt{double}\}\) and the component count \(K\in\{2,3,4\}\); the device EFTs use the true hardware FMA. With --fmad=false (and host-side -ffp-contract=off) forbidding automatic FMA contraction of \(a\cdot b+c\) and thereby guaranteeing the EFTs, the device FMA agrees bitwise with the CPU/NEON scalar reference for all six types (200,000 trials each, 0 mismatches).

4.2.8.1 Correcting the baseline.

As in the CPU/AVX-512/NEON evaluations, the fair baseline is the branch-free BF (Zhang–Aiken mul_bf+add_bf), implemented on the device. An earlier version of this paper reported \(\times 1.8\)–\(\times 12.6\) against compose (fully distilled expansion mul+add); these figures were inflated and are retracted. compose costs \(248/810/1856\) flops at \(K{=}2/3/4\), i.e.more than the Exact FMA (symmetric, \(194/432/968\)), and measures slower than Exact. It is not a realistic mul+add implementation and was an inappropriate baseline; below it appears only as the reference column compose\(^{*}\).

2pt

2pt

4.2.8.2 Machines where FMA/BF is \(\times 1.00\) and machines where it is \(\times 1.4\)–\(\times 1.7\).

With the fair BF baseline, the GEMM FMA/BF is \(\times 0.84\)–\(\times 1.02\) on the H100—the \(\times 1.3\)–\(\times 1.9\) advantage seen on CPU/NEON does not appear—whereas the three double-precision types on GB10 reach \(\times 1.40\)–\(\times 1.71\) (Table 18), matching the theoretical FLOP ratios of each type (DW \(1.71\) / TW \(1.45\) / QW \(1.43\); Table 2). In particular DD: BF \(29\) flops versus the proposed \(17\) flops gives ratio \(1.71\), matching the measured \(\times 1.71\). What separates the two regimes is whether the kernel is compute-bound or bandwidth-bound (our GEMM/AXPY are naive, without shared-memory tiling). GB10 has very low FP64 throughput (DD GEMM at \(11.2\) GMAC/s, about \(1/9\) of the H100’s \(106\) GMAC/s), so its three double-precision types become compute-bound and the FLOP reduction shows directly in wall-clock time. On the same GB10 the three single-precision types are bandwidth-bound with FMA/BF \(\approx 1.01\), and on the H100, whose FP64 is strong (\(1/2\) of FP32), all six types fall to the bandwidth/latency-bound side.

To confirm this interpretation independently, a register-resident compute-only MAC chain without memory traffic gives FMA/BF of \(\times 1.21\)–\(\times 1.62\) on H100 (DS \(1.58\) / TS \(1.46\) / QS \(1.26\) / DD \(1.62\) / TD \(1.42\) / QD \(1.21\)) and \(\times 1.21\)–\(\times 1.67\) on GB10 (DS \(1.54\) / TS \(1.47\) / QS \(1.24\) / DD \(1.67\) / TD \(1.42\) / QD \(1.21\))—the two machines agree. DW (theory \(1.71\)) and TW (\(1.45\)) reproduce the theoretical ratios well, while QW stays at \(1.21\)–\(1.24\), short of the theoretical \(1.43\): with its long dependence chain and high register pressure, the chain-shaped measurement leans latency-bound for QW, and the memory-carrying GEMM (\(\times 1.40\)) is actually closer to theory. In other words, the operation-count advantage is real on GPUs, and it appears in wall-clock time whenever the kernel is compute-bound. Extracting it on the bandwidth-bound side requires making the kernel compute-bound by tiling and register blocking—future work.

Against the Exact FMA, the proposed FMA is \(\times 1.30\)–\(\times 4.08\) faster on H100 and \(\times 1.34\)–\(\times 9.09\) on GB10, consistent with the operation-count ratios (the symmetric Exact costs more, so the ratios grew slightly). Throughput of the proposed FMA in GEMM: H100 single base DS \(209\) / TS \(137\) / QS \(103\) GMAC/s and double base DD \(106\) / TD \(59\) / QD \(52\) GMAC/s; GB10 DS \(101\) / TS \(68\) / QS \(51\) GMAC/s and DD \(11.2\) / TD \(2.9\) / QD \(1.3\) GMAC/s. In conclusion, the proposed FMA agrees bitwise across CPU/NEON/SVE2/GPU, has the smallest operation count, and its wall-clock advantage on GPUs is obtained whenever the kernel is compute-bound.

5 Conclusion and future work↩︎

• A machine-verified minimal construction: the proposed FMA is the minimal branch-free fused FMA whose error bounds (\(34u^2\)/\(184u^3\)/\(812u^4\)) and all FastTwoSum preconditions are machine-verified with FPANVerifier (DW \(17\) / TW \(66\) / QW \(146\) flops, Section 2.9). Choosing each renormalization gate by “is the FastTwoSum precondition provable for all precisions \(p\)” preserves the property that FastTwoSum is used only with proven preconditions. The MPFR falsification test (400,000 trials per type per backend, 0 violations) supports the bounds, and a negative control confirms the test’s detection power.

• It is commutative in \(x\leftrightarrow y\) (bitwise agreement for all six types, 200,000 trials each, 0 mismatches; Section 2.10). two_prod/two_sum are EFTs and hence commutative from the outset; only the association order of the diagonal sums breaks commutativity. Re-associating so that transposition partners pair first restores commutativity without adding a single operation, and as a by-product the certified QW bound shrinks from \(824u^4\) to \(812u^4\). The commutativity of the existing BF mul/add is thus not sacrificed by the fused FMA.

• The proven variant agrees bitwise with the naive one (\(5\times10^6\) trials per type, 0 mismatches), so the relative errors are unchanged; migration changes only flops (DW \(20\to17\), QW \(143\to146\), TW unchanged).

• The proposed FMA is faster than branch-free BF (AVX2 GEMM: DD \(\times 1.49\) / TD \(\times 1.32\) / QD \(\times 1.47\)), reflecting the proven-variant operation counts (\(17/66/146\)). DD gains most from the FastTwoSum-ized tail; QW loses slightly from the two FastTwoSum gates demoted to TwoSum.

• Most of the existing speedup is due to Zhang–Aiken’s BF; the proposed FMA stacks on top of it. In environments with only the old Bailey variant (Q), the apparent gain is \(\times 2.8\)–\(\times 6.6\).

• The Exact FMA is useful where accuracy is paramount (about one digit better, result-relative) but runs 2.5–5.0 times slower than BF (its operation count grew from \(182/396/896\) to \(194/432/968\) with the symmetric variant).

• All four variants maintain the working accuracy, and under cancellation the proposed FMA behaves numerically almost identically to the existing ones (Q/BF).

• Doubling the SIMD width to AVX-512 leaves rounding bitwise identical; the proposed FMA is \(\times 1.03\)–\(\times 1.77\) versus BF and \(\times 1.68\)–\(\times 1.73\) faster than its own AVX2 GEMM (Section 4.2.3). Being branch-free, it scales cleanly with lane width.

• Porting the same evaluation to Arm/NEON leaves rounding bitwise identical, with FMA/BF \(\times 1.26\)–\(\times 2.01\) (Section 4.2.4); the advantage is consistent across x86 and Arm.

• Extending to the single-precision base DS/TS/QS keeps the working accuracy, with FMA/BF \(\times 1.37\)–\(\times 1.83\) (Section 4.2.5).

• Porting to Arm/SVE2 leaves rounding bitwise identical, with FMA/BF \(\times 1.29\)–\(\times 1.88\) (Section 4.2.6). The test machine has VL\(=128\) bit, i.e.the same lane width as NEON, and the absolute times match within \(\pm15\%\), showing that the advantage is ISA-independent; QD GEMM attains its best \(\times 1.62\) on SVE2. The sizeless-type constraint of SVE causes stack spills of the Exact FMA’s expansion array (\(\times 2.5\) slower than NEON); the proposed FMA uses no expansion array and is unaffected.

• The CUDA implementations of all six types agree bitwise with CPU/NEON/SVE2 (Section 4.2.8). With the fair BF baseline, the FMA/BF of the naive (untiled) GEMM kernel is decided by whether the kernel is compute- or bandwidth-bound: it stays at \(\times 0.84\)–\(\times 1.02\) on the bandwidth-bound H100 (all six types) and GB10 (three single types), but reaches \(\times 1.40\)–\(\times 1.71\) for the compute-bound double-precision types on GB10, matching the theoretical FLOP ratios (DW \(1.71\) / TW \(1.45\) / QW \(1.43\)). Register-resident compute-only MAC chains give \(\times 1.21\)–\(\times 1.67\) on both machines, so the operation-count advantage is real on GPUs; extracting it on the bandwidth-bound side requires making kernels compute-bound.

• The benefit of the proposed FMA is not limited to MAC kernels. Division (long division) and square root (Newton iteration for \(1/\sqrt{a}\)) have inner loops made entirely of \(x\cdot y+c\) MACs and accelerate for all six types (Section 4.1): division at NEON \(\times 1.54\)–\(\times 2.08\) / SVE2 \(\times 1.54\)–\(\times 2.11\) / GB10 compute-only \(\times 1.25\)–\(\times 1.61\); square root at NEON \(\times 1.51\)–\(\times 1.58\) / SVE2 \(\times 1.59\)–\(\times 1.61\) / GB10 compute-only \(\times 1.25\)–\(\times 1.29\). QD division exceeds its operation-count ratio \(1.51\), reaching \(\times 2.03\)–\(\times 2.11\), because the two serial normalization chains of BF fold into one and the dependence chain shortens.

• For square root the proposed FMA loses no accuracy at all: the Newton iteration converges far below \(\varepsilon_K\), so the final error is dominated by the last multiplication \(a\cdot r\) common to all three methods, and the maximum errors of BF and the proposed FMA agree bitwise for DS/DD. A pure gain of \(\times 1.51\)–\(\times 1.61\) at zero accuracy cost; for the same reason there is little point in an Exact-FMA square root. For division the proposed FMA’s error is \(1\)–\(3\) times BF’s, all within a few \(\varepsilon_K\) and practically equivalent (the Exact FMA has \(1/4\)–\(1/8\) of BF’s error at \(1/3\)–\(1/8\) of the speed).

• In division and square root as well, GB10’s three single types are bandwidth-bound (effective \(238\)–\(243\) GB/s against the LPDDR5X peak of about \(273\) GB/s) with elementwise FMA/BF at \(\times 1.00\), while the register-resident compute-only chains match the theoretical operation-count ratios for all six types—reproducing the GEMM conclusion that the advantage appears in wall-clock time when compute-bound.

• The figure of “\(\times 1.8\)–\(\times 12.6\) over non-fused compose on GPUs” reported in an earlier version of this paper is retracted: the compose baseline (fully distilled mul+add, \(248/810/1856\) flops) is more expensive than the Exact FMA and measured slower; it was an inappropriate baseline.

5.0.0.1 Future work.

Within the current renormalization structure, \(17/66/146\) is the minimum of the proven variant (the FastTwoSum candidates that would lower TW to \(63\) or QW to \(143\) all fail the SMT check; Section 2.9, Table 4). This is, however, a minimum within a given structure, not a global one. Directions for improvement include (i) redesigning the renormalization structure (a Priest-style renormalization, in which FastTwoSum chains are exact on magnitude-sorted expansions, could make the whole renormalization FastTwoSum and might lower QW from \(146\) to \(140\) or below); (ii) provably minimal three_sum/vec_sum accumulation (lemma search with FPANVerifier’s conjecturer); and (iii) fma-folded diagonals (DW \(16\)/TW \(64\)/QW \(140\), \(1\)–\(3\) flops fewer, but the single rounding of a fused operation cannot be expressed as two_prod\(+\)two_sum and hence cannot be machine-verified). All ask “can we go faster while staying machine-verifiable?”—the next stage of this research.

Acknowledgments↩︎

This work was supported by JSPS KAKENHI Grant Number 26K14846.

6 Appendix: a GCC 13.3 SVE2 miscompilation encountered in this work↩︎

When symmetrizing the Exact FMA added 12 TwoSum gates, we encountered a phenomenon in which the proposed FMA on SVE2/QD disagreed with the scalar reference in \(12{,}238/200{,}000\) trials. Isolation showed that this is not an error of either algorithm but a code-generation bug of GCC 13.3, based on the following evidence.

• The output of the Exact FMA itself is correct (its own mismatches are \(0\)). What gets corrupted are the values of the proposed FMA computed and stored before the call.

• Without calling the Exact FMA, the proposed FMA has \(0\) mismatches.

• At -O0 / -O1 and under AddressSanitizer there are \(0\) mismatches; the issue reproduces only at -O2 / -O3, and ASan detects no out-of-bounds access.

• The same algorithms do not reproduce it on NEON / AVX2 / AVX-512 / CUDA (SVE2-specific).

That is, inlining a spill-heavy function that uses fixed-length SVE vectors (-msve-vector-bits=128) corrupts values of the caller. The trigger was the register pressure crossing a threshold when the 12 TwoSum gates were added to the QW Exact’s \(v[36]\) (as noted in Section 4.2.6, QW Exact already spilled with \(50\) st1d stores and is the most spill-heavy SVE2 function).

As a workaround we attach __attribute__((noinline)) to the three SVE2 Exact-FMA functions. The Exact FMA is a reference implementation for accuracy bounds, and one call’s overhead is negligible. With the attribute, all types return to \(0\) mismatches at -O2/-O3. The bug was latent in the pre-symmetrization code as well, merely not manifest because the register pressure had not crossed the threshold. Implementations using fixed-length SVE types should beware of this kind of corruption.

7 Appendix: reference implementation source code↩︎

The reference implementations of the proposed FMA and the Exact FMA (scalar C; the baseline of the bitwise-agreement checks of all backends) are listed below. Both must be compiled with -ffp-contract=off (or #pragma STDC FP_CONTRACT OFF).

7.1 Proposed FMA (fma_ref.c)↩︎

/********************************************************************************/
/* fma_ref.c : Branch-free fused multiply-add for DW/TW/QW arithmetic           */
/*             (scalar reference, machine-verified with FPANVerifier)           */
/* Copyright (C) 2026 Tomonori Kouya                                            */
/*                                                                              */
/* This program is free software: you can redistribute it and/or modify it      */
/* under the terms of the GNU Lesser General Public License as published by the */
/* Free Software Foundation, either version 3 of the License or any later       */
/* version.                                                                     */
/*                                                                              */
/* This program is distributed in the hope that it will be useful, but WITHOUT  */
/* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or        */
/* FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public License */
/* for more details.                                                            */
/*                                                                              */
/* You should have received a copy of the GNU Lesser General Public License     */
/* along with this program.  If not, see <http://www.gnu.org/licenses/>.        */
/*                                                                              */
/********************************************************************************/
/*
 * Reference-implementation notes:
 * Each routine matches fpan/{dw,tw,qw}fma.fpan exactly, operation by
 * operation and in the same order ("certified variants").  Every plain
 * add / plain mul corresponds to a discarded wire (one rounding error
 * or a truncated term) in the .fpan file.  The faster fma()-folded
 * variants (folding the diagonal terms into one value by an fma chain)
 * save DW:-1, TW:-2, QW:-3 flops but cannot be modeled by FPANVerifier
 * and would require manual analysis.
 * Compile with: -O2 -march=native  (FP_CONTRACT off recommended:
 *   #pragma STDC FP_CONTRACT OFF or -ffp-contract=off; if the compiler
 *   contracts a*b+c into fma on its own, the code no longer matches
 *   the verified netlist).
 */
#pragma STDC FP_CONTRACT OFF
#include <math.h>

typedef struct { double hi, lo; }         dw_t;
typedef struct { double x0, x1, x2; }     tw_t;
typedef struct { double x0, x1, x2, x3; } qw_t;

static inline void two_sum(double a, double b, double *s, double *e) {
    double t;
    *s = a + b;
    t  = *s - a;
    *e = (a - (*s - t)) + (b - t);
}
static inline void fast_two_sum(double a, double b, double *s, double *e) {
    *s = a + b;
    *e = b - (*s - a);
}
static inline void two_prod(double a, double b, double *p, double *e) {
    *p = a * b;
    *e = fma(a, b, -*p);
}

/* ---------------------------------------------------------------
 * DW-FMA: z = x*y + c   (17 flops, matches dwfma_f2s.fpan; trailing FastTwoSum)
 * Error: |z-(xy+c)| <= 34 u^2 (|xy|+|c|)  (FPANVerifier-certified, u=2^-p)
 * --------------------------------------------------------------- */
dw_t dw_fma(dw_t x, dw_t y, dw_t c) {
    double P00, E00, P01, P10, l, v, w, s, t, tp;
    dw_t z;
    two_prod(x.hi, y.hi, &P00, &E00);       /* 2 */
    P01 = x.hi * y.lo;                      /* 1  (E01 discarded) */
    P10 = x.lo * y.hi;                      /* 1  (E10 discarded) */
    l   = P01 + P10;                        /* 1  (eps discarded) */
    v   = E00 + c.lo;                       /* 1  (delta5 discarded) */
    w   = v + l;                            /* 1  (delta6 discarded) */
    two_sum(P00, c.hi, &s, &t);             /* 6 */
    tp  = t + w;                            /* 1  (delta7 discarded) */
    fast_two_sum(s, tp, &z.hi, &z.lo);      /* 3  trailing FastTwoSum:
                                                  exp(s)>=exp(tp) machine-proved by FPANVerifier
                                                  (dwfma_f2s.fpan). 20->17 flops */
    return z;
}

/* ---------------------------------------------------------------
 * TW-FMA: z = x*y + c   (66 flops, matches twfma_fix2.fpan; gates swapped)
 * Error: |z-(xy+c)| <= 184 u^3 (|xy|+|c|)  (FPANVerifier-certified)
 * --------------------------------------------------------------- */
tw_t tw_fma(tw_t x, tw_t y, tw_t c) {
    double P00,E00,P01,E01,P10,E10,P02,P11,P20,sg,G;
    double A,q1,q2,q3,B,r,m1,m2,w0,w1,w2,rho;
    tw_t z;
    two_prod(x.x0, y.x0, &P00, &E00);                 /* 2 */
    two_prod(x.x0, y.x1, &P01, &E01);                 /* 2 */
    two_prod(x.x1, y.x0, &P10, &E10);                 /* 2 */
    P02 = x.x0 * y.x2;                                /* 1 */
    P11 = x.x1 * y.x1;                                /* 1 */
    P20 = x.x2 * y.x0;                                /* 1 */
    sg  = (P02 + P20) + P11;                          /* 2 diagonal i+j=2 */
    G   = ((E01 + E10) + sg) + c.x2;                  /* 3 */
    two_sum(P01, P10, &A, &q1);                       /* 6 */
    two_sum(A,   E00, &A, &q2);                       /* 6 */
    two_sum(A,   c.x1,&A, &q3);                       /* 6 */
    G   = G + ((q1 + q2) + q3);                       /* 3 */
    two_sum(P00, c.x0, &B, &r);                       /* 6 */
    two_sum(r, A, &m1, &m2);                          /* 6 */
    m2  = m2 + G;                                     /* 1 */
    fast_two_sum(B,  m1, &w0, &w1);                   /* 3 renormalize (twfma_fix2: provable) */
    two_sum(w1, m2, &w1, &w2);                        /* 6 */
    two_sum(w0, w1, &z.x0, &rho);                     /* 6 (old FastTwoSum unproved -> TwoSum) */
    fast_two_sum(rho, w2, &z.x1, &z.x2);              /* 3 */
    return z;
}

/* ---------------------------------------------------------------
 * QW-FMA: z = x*y + c   (146 flops, matches qwfma_fix3.fpan; two unproved F2S demoted to TwoSum)
 * Error: |z-(xy+c)| <= 812 u^4 (|xy|+|c|)  (FPANVerifier-certified, symmetrized)
 * Note: the 9 sequential two_sum of L2 form a long dependence chain.
 *       For more ILP the pairing may be tree-shaped (flops unchanged; re-verify).
 * --------------------------------------------------------------- */
qw_t qw_fma(qw_t x, qw_t y, qw_t c) {
    double P00,E00,P01,E01,P10,E10,P02,E02,P11,E11,P20,E20;
    double P03,P12,P21,P30,D,B,r;
    double A1,f1,f2,f3,f4;
    double A2,g1,g2,g3,g4,g5,g6,g7,g8,g9;
    double A3,t1,t2,t3,t4;
    double w0,w1,w2,w3,rho,sig;
    qw_t z;
    two_prod(x.x0, y.x0, &P00, &E00);                 /* 2 */
    two_prod(x.x0, y.x1, &P01, &E01);                 /* 2 */
    two_prod(x.x1, y.x0, &P10, &E10);                 /* 2 */
    two_prod(x.x0, y.x2, &P02, &E02);                 /* 2 */
    two_prod(x.x1, y.x1, &P11, &E11);                 /* 2 */
    two_prod(x.x2, y.x0, &P20, &E20);                 /* 2 */
    P03 = x.x0 * y.x3;                                /* 1 */
    P12 = x.x1 * y.x2;                                /* 1 */
    P21 = x.x2 * y.x1;                                /* 1 */
    P30 = x.x3 * y.x0;                                /* 1 */
    D   = (P03 + P30) + (P12 + P21);                  /* 3 diagonal i+j=3 */
    two_sum(P00, c.x0, &B, &r);                       /* 6  L0 */
    /* L1 */
    two_sum(P01, P10,  &A1, &f1);                     /* 6 */
    two_sum(A1,  E00,  &A1, &f2);                     /* 6 */
    two_sum(A1,  c.x1, &A1, &f3);                     /* 6 */
    two_sum(A1,  r,    &A1, &f4);                     /* 6 */
    /* L2 */
    two_sum(P02, P20,  &A2, &g1);                     /* 6 */
    two_sum(A2,  P11,  &A2, &g2);                     /* 6 */
    two_sum(E01, E10, &E01, &g4);                     /* 6 */   /* SYM: pair E01,E10 first (commutative) */
    two_sum(A2, E01, &A2, &g3);                     /* 6 */
    two_sum(A2,  c.x2, &A2, &g5);                     /* 6 */
    two_sum(A2,  f1,   &A2, &g6);                     /* 6 */
    two_sum(A2,  f2,   &A2, &g7);                     /* 6 */
    two_sum(A2,  f3,   &A2, &g8);                     /* 6 */
    two_sum(A2,  f4,   &A2, &g9);                     /* 6 */
    /* L3: rounded sums (same order as the discarded wires of qwfma.fpan) */
    t1 = E02 + E20;                                   /* 1 */
    t2 = E11 + D;                                     /* 1 */
    t3 = (t1 + t2) + c.x3;                            /* 2 */
    t1 = g1 + g2;                                     /* 1 */
    t2 = g3 + g4;                                     /* 1 */
    t1 = t1 + t2;                                     /* 1 */
    t2 = g6 + g7;                                     /* 1 */
    t4 = g8 + g9;                                     /* 1 */
    t2 = t2 + t4;                                     /* 1 */
    t1 = t1 + t2;                                     /* 1 */
    t1 = t1 + g5;                                     /* 1 */
    A3 = t3 + t1;                                     /* 1  13 adds in total */
    /* renormalization (qwfma_fix3: only (B,A1) and (sig,w3) admit provable FastTwoSum) */
    fast_two_sum(B,  A1, &w0, &w1);                   /* 3 */
    two_sum(w1, A2, &w1, &w2);                        /* 6 */
    two_sum(w2, A3, &w2, &w3);                        /* 6 */
    two_sum(w0, w1, &z.x0, &rho);                     /* 6 (old FastTwoSum unproved -> TwoSum) */
    two_sum(rho, w2, &z.x1, &sig);                    /* 6 (old FastTwoSum unproved -> TwoSum) */
    fast_two_sum(sig, w3, &z.x2, &z.x3);              /* 3 */
    return z;
}

7.2 Exact FMA (fma_exact_ref.c)↩︎

/********************************************************************************/
/* fma_exact_ref.c : Exact (result-relative) fused multiply-add                 */
/*                   for DW/TW/QW arithmetic (distillation-based reference)     */
/* Copyright (C) 2026 Tomonori Kouya                                            */
/*                                                                              */
/* This program is free software: you can redistribute it and/or modify it      */
/* under the terms of the GNU Lesser General Public License as published by the */
/* Free Software Foundation, either version 3 of the License or any later       */
/* version.                                                                     */
/*                                                                              */
/* This program is distributed in the hope that it will be useful, but WITHOUT  */
/* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or        */
/* FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public License */
/* for more details.                                                            */
/*                                                                              */
/* You should have received a copy of the GNU Lesser General Public License     */
/* along with this program.  If not, see <http://www.gnu.org/licenses/>.        */
/*                                                                              */
/********************************************************************************/
/*
 * Reference-implementation notes:
 * Semantics: |z - (x*y+c)| <= K u^n |x*y+c|  (result-relative, robust
 * to cancellation) -- a different semantics from the (|xy|+|c|)-relative
 * cheap variant (fma_ref.c).  All n^2 TwoProd terms are kept, a
 * rounding-free TwoSum distillation sweep is applied, and only the top
 * n words are output.  The gate order matches
 * fpan/{dw,tw,qw}fma_exact.fpan one-to-one (generated together by
 * gen_exact_fma.py).  Compile with: -O2 -ffp-contract=off (mandatory).
 */
#pragma STDC FP_CONTRACT OFF
#include <math.h>

typedef struct { double hi, lo; }         dw_t;
typedef struct { double x0, x1, x2; }     tw_t;
typedef struct { double x0, x1, x2, x3; } qw_t;

static inline void two_sum(double a, double b, double *s, double *e) {
    double t;
    *s = a + b;
    t  = *s - a;
    *e = (a - (*s - t)) + (b - t);
}
static inline void two_prod(double a, double b, double *p, double *e) {
    *p = a * b;
    *e = fma(a, b, -*p);
}

/* dw_fma_exact: z = fl_2(x*y+c), result-relative
 * semantics. 4 TwoProd + 29 TwoSum =
 * 182 flops. Gate order matches
 * dwfma_exact.fpan exactly. */
dw_t dw_fma_exact(dw_t x, dw_t y, dw_t c) {
    double v[10], e;
    dw_t z;
    v[1] = c.hi;          /* c0 */
    v[4] = c.lo;          /* c1 */
    two_prod(x.hi, y.hi, &v[0], &v[5]); /* P00,E00 */
    two_prod(x.hi, y.lo, &v[2], &v[7]); /* P01,E01 */
    two_prod(x.lo, y.hi, &v[3], &v[8]); /* P10,E10 */
    two_prod(x.lo, y.lo, &v[6], &v[9]); /* P11,E11 */
    /* SYM: TwoSum each transposed pair (P_ij,P_ji) / (E_ij,E_ji) before the sweep,
       so the expansion depends only on the multiset -> Exact FMA is commutative. */
    two_sum(v[2], v[3], &v[2], &v[3]);
    two_sum(v[7], v[8], &v[7], &v[8]);
    for (int p = 0; p < 3; p++)       /* distillation sweeps */
        for (int i = 10 - 2; i >= 0; i--)
            two_sum(v[i], v[i+1], &v[i], &v[i+1]);
    for (int i = 2 - 1; i >= 0; i--) /* top refresh */
        two_sum(v[i], v[i+1], &v[i], &v[i+1]);
    z.hi = v[0];
    z.lo = v[1];
    (void)e;
    return z;   /* v[2]..v[9] discarded */
}
/* tw_fma_exact: z = fl_3(x*y+c), result-relative
 * semantics. 9 TwoProd + 63 TwoSum =
 * 396 flops. Gate order matches
 * twfma_exact.fpan exactly. */
tw_t tw_fma_exact(tw_t x, tw_t y, tw_t c) {
    double v[21], e;
    tw_t z;
    v[1] = c.x0;          /* c0 */
    v[4] = c.x1;          /* c1 */
    v[9] = c.x2;          /* c2 */
    two_prod(x.x0, y.x0, &v[0], &v[5]); /* P00,E00 */
    two_prod(x.x0, y.x1, &v[2], &v[10]); /* P01,E01 */
    two_prod(x.x0, y.x2, &v[6], &v[14]); /* P02,E02 */
    two_prod(x.x1, y.x0, &v[3], &v[11]); /* P10,E10 */
    two_prod(x.x1, y.x1, &v[7], &v[15]); /* P11,E11 */
    two_prod(x.x1, y.x2, &v[12], &v[18]); /* P12,E12 */
    two_prod(x.x2, y.x0, &v[8], &v[16]); /* P20,E20 */
    two_prod(x.x2, y.x1, &v[13], &v[19]); /* P21,E21 */
    two_prod(x.x2, y.x2, &v[17], &v[20]); /* P22,E22 */
    /* SYM: TwoSum each transposed pair (P_ij,P_ji) / (E_ij,E_ji) before the sweep,
       so the expansion depends only on the multiset -> Exact FMA is commutative. */
    two_sum(v[2], v[3], &v[2], &v[3]);
    two_sum(v[6], v[8], &v[6], &v[8]);
    two_sum(v[10], v[11], &v[10], &v[11]);
    two_sum(v[12], v[13], &v[12], &v[13]);
    two_sum(v[14], v[16], &v[14], &v[16]);
    two_sum(v[18], v[19], &v[18], &v[19]);
    for (int p = 0; p < 3; p++)       /* distillation sweeps */
        for (int i = 21 - 2; i >= 0; i--)
            two_sum(v[i], v[i+1], &v[i], &v[i+1]);
    for (int i = 3 - 1; i >= 0; i--) /* top refresh */
        two_sum(v[i], v[i+1], &v[i], &v[i+1]);
    z.x0 = v[0];
    z.x1 = v[1];
    z.x2 = v[2];
    (void)e;
    return z;   /* v[3]..v[20] discarded */
}
/* qw_fma_exact: z = fl_4(x*y+c), result-relative
 * semantics. 16 TwoProd + 144 TwoSum =
 * 896 flops. Gate order matches
 * qwfma_exact.fpan exactly. */
qw_t qw_fma_exact(qw_t x, qw_t y, qw_t c) {
    double v[36], e;
    qw_t z;
    v[1] = c.x0;          /* c0 */
    v[4] = c.x1;          /* c1 */
    v[9] = c.x2;          /* c2 */
    v[16] = c.x3;          /* c3 */
    two_prod(x.x0, y.x0, &v[0], &v[5]); /* P00,E00 */
    two_prod(x.x0, y.x1, &v[2], &v[10]); /* P01,E01 */
    two_prod(x.x0, y.x2, &v[6], &v[17]); /* P02,E02 */
    two_prod(x.x0, y.x3, &v[12], &v[23]); /* P03,E03 */
    two_prod(x.x1, y.x0, &v[3], &v[11]); /* P10,E10 */
    two_prod(x.x1, y.x1, &v[7], &v[18]); /* P11,E11 */
    two_prod(x.x1, y.x2, &v[13], &v[24]); /* P12,E12 */
    two_prod(x.x1, y.x3, &v[20], &v[29]); /* P13,E13 */
    two_prod(x.x2, y.x0, &v[8], &v[19]); /* P20,E20 */
    two_prod(x.x2, y.x1, &v[14], &v[25]); /* P21,E21 */
    two_prod(x.x2, y.x2, &v[21], &v[30]); /* P22,E22 */
    two_prod(x.x2, y.x3, &v[27], &v[33]); /* P23,E23 */
    two_prod(x.x3, y.x0, &v[15], &v[26]); /* P30,E30 */
    two_prod(x.x3, y.x1, &v[22], &v[31]); /* P31,E31 */
    two_prod(x.x3, y.x2, &v[28], &v[34]); /* P32,E32 */
    two_prod(x.x3, y.x3, &v[32], &v[35]); /* P33,E33 */
    /* SYM: TwoSum each transposed pair (P_ij,P_ji) / (E_ij,E_ji) before the sweep,
       so the expansion depends only on the multiset -> Exact FMA is commutative. */
    two_sum(v[2], v[3], &v[2], &v[3]);
    two_sum(v[6], v[8], &v[6], &v[8]);
    two_sum(v[10], v[11], &v[10], &v[11]);
    two_sum(v[12], v[15], &v[12], &v[15]);
    two_sum(v[13], v[14], &v[13], &v[14]);
    two_sum(v[17], v[19], &v[17], &v[19]);
    two_sum(v[20], v[22], &v[20], &v[22]);
    two_sum(v[23], v[26], &v[23], &v[26]);
    two_sum(v[24], v[25], &v[24], &v[25]);
    two_sum(v[27], v[28], &v[27], &v[28]);
    two_sum(v[29], v[31], &v[29], &v[31]);
    two_sum(v[33], v[34], &v[33], &v[34]);
    for (int p = 0; p < 4; p++)       /* distillation sweeps */
        for (int i = 36 - 2; i >= 0; i--)
            two_sum(v[i], v[i+1], &v[i], &v[i+1]);
    for (int i = 4 - 1; i >= 0; i--) /* top refresh */
        two_sum(v[i], v[i+1], &v[i], &v[i+1]);
    z.x0 = v[0];
    z.x1 = v[1];
    z.x2 = v[2];
    z.x3 = v[3];
    (void)e;
    return z;   /* v[4]..v[35] discarded */
}

References↩︎