Multiword Arithmetic and Parallel Computing

Abstract: In many applications, the precision by the available hardware arithmetic is insufficient to guarantee accurate results. Multiword arithmetic is a special type of multiprecision arithmetic where a multiple double is an unevaluated sum of 64-bit doubles, or where a multiple integer is an unevaluated sum of 64-bit integers. Parallel computing is applied to compensate for the cost overhead of multiword arithmetic. This type of arithmetic exploits naturally the optimized hardware, allows for efficient type conversions, memory layouts, all favorable for parallel computing. For example, storing a multiword in registers rather than arrays is beneficial to parallel computing by tasking and acceleration by graphics processing units. Code for multiword arithmetic is available in the software PHCpack, written mainly in Ada, publicly available at github, and as an Alire crate, released under the GNU GPL v3.0 license.

• The paper introduces PHCpack implementations of multiword doubles and integers, using error-free transformations and deferred renormalization to support robust high-precision computation.
• The paper restructures double-double inner products as double-precision convolution and matrix multiplication, achieving vectorized speedups of up to 6.2× over ordinary multiword methods.
• The paper shows that 96-thread execution can compute hexa-double inner products in 293 milliseconds, but its method currently depends on aligned exponents and bounded addition counts.

This paper by Jan Verschelde (2606.20863) describes recent additions to PHCpack implementing multiword arithmetic — multiple doubles (unevaluated sums of nonoverlapping 64-bit doubles) and multiple integers (unevaluated sums of 64-bit integers) — together with a vectorized inner product algorithm whose cost overhead is offset by parallelism. The work builds on the classical error-free transformation tradition originating with Dekker's techniques, with prior software including QDlib and CAMPARY, and extends the author's earlier contributions on parallel higher-precision computation in Ada.

Motivation and cost overhead

The motivating application is robustness in the sense of Shewchuk: algorithms should not fail under small perturbations of degenerate inputs, and multiprecision arithmetic is one route to robustness. A specific driver is power series computation for polynomial homotopies, where errors in leading coefficients propagate to trailing coefficients, so leading coefficients must be computed beyond double precision. The paper also frames this code as computational preparation for GPU tensor cores via the Ozaki scheme.

The central obstacle is overhead. The paper tabulates operation counts showing that an $m$-double addition requires $m=2$: 20 doubles ops for add, 23 for mul, 70 for div; at $m=16$ these grow to 925, 11,499, and 33,041 respectively. These steep factors motivate both algorithmic restructuring and parallel execution as compensation.

Error-free summation and vectorized inner products

The core algorithmic idea assumes all doubles share the same exponent, reducing the problem to exact summation of 52-bit fractions. Doubles are split into quarters; if the number of additions stays below a threshold, enough zero bits remain at the tails of the operands to represent partial sums exactly, yielding an error-free accumulation. The paper notes this assumption explicitly — it holds only when exponents are aligned and addition counts are bounded, which constrains generality of the method.

For inner products $\sum_k x_k \star y_k$ over length-$n$ vectors of double doubles, each component $(x^{\mathrm{hi}}_k, x^{\mathrm{lo}}_k)$ is split into eight quarters $x_{k,0},\ldots,x_{k,7}$. The convolution structure

$$s_i = \sum_{k=1}^n \sum_{j=0}^i x_{k,j}\, y_{k,i-j}$$

is then evaluated entirely in double precision, with the eight partial sums $s_0,\ldots,s_7$ combined at the end in double double arithmetic. The key consequence is that multiword inner products reduce to matrix multiplications in plain double precision — precisely the workload GPUs and tensor cores accelerate best. A further efficiency observation is that postponing renormalization of the multiple doubles benefits performance.

The test data uses random doubles whose 52-bit fractions follow a pattern of four 12-bit blocks separated by leading ones, so that all quarters have fixed exponents (e.g., 0, −13, −26, −39). This is a deliberate experimental design choice that guarantees the exponent-alignment assumption; results on arbitrary inputs are not claimed.

Computational results

Benchmarks compute 1,024 inner products of length 6,144 on an Intel Xeon 5318Y Ice Lake-SP, compiled with GNAT 12.2.0 at -O3:

Precision	Ordinary time	Overhead vs. 1d	Vectorized time	Vectorized speedup
1d	12 ms	—	30 ms	0.4x
2d	158 ms	13x	69 ms	2.3x
4d	1 s 977 ms	12x	318 ms	6.2x
8d	6 s 428 ms	3.3x	1 s 520 ms	4.2x
16d	40 s 780 ms	6.3x	9 s 491 ms	4.3x

Vectorization yields speedups of roughly 2–6x over ordinary multiword inner products, with quad double benefiting most (6.2x). Notably, single-double vectorized code is slower than the scalar version (0.4x), indicating the splitting overhead only pays off once precision exceeds hardware precision. In hexa double precision, 1,024 inner products complete in about 9 seconds wall clock.

Multithreading amplifies the effect: with 96 threads across two 24-core Xeons, each thread handling one inner product, wall clock drops to 293 milliseconds. The paper highlights a "quality up" result: hexa double (16d) computation with 96 threads runs faster than quad double (4d) on a single thread (293 ms vs. 318 ms) — quadrupling the effective precision at no wall-clock cost.

Limitations and open questions

Several caveats bear directly on the reported numbers. The error-free summation argument depends on equal exponents across operands and on bounded addition counts leaving sufficient zero bits; the benchmark generator enforces fixed per-quarter exponents, so behavior on general data with wide dynamic range is not established. The multithreaded comparison also involves different thread counts on either side of the quality-up claim, so the result reflects throughput scaling rather than a like-for-like kernel speedup. Finally, the reduction to double-precision matrix multiplication is presented as preparation for tensor-core acceleration via the Ozaki scheme; actual GPU implementation and its achieved accuracy and throughput remain open questions this paper does not answer.

Conclusion

The paper demonstrates that restructuring multiword inner products as convolutions evaluated in double precision, combined with deferred renormalization and task-level parallelism, substantially offsets the cost of extended precision: up to 4.3x vectorized speedup and effective quadrupling of precision at constant wall-clock time using 96 threads. The implementation is available in PHCpack under GPL v3.0, distributed via GitHub and as an Alire crate.