Algorithm

Double-word and multi-double arithmetic

Cost: O(N^2) double operations for a product of N-term expansions, O(N^2) for a long division with tapered remainders; about N^2 fused multiply-adds dominate.

A floating-point expansion represents a real number as an unevaluated sum x = x_0 + x_1 + \cdots + x_{N-1} of machine doubles, each much smaller than the previous one. With N = 2 (double-word or double-double arithmetic) this gives about 106 bits of precision, with N = 4 (quad-double) about 212 bits, using only hardware floating-point operations. The idea goes back to Møller [Mol1965] and Dekker [Dek1971]; general expansions are due to Priest [Pri1991] and Shewchuk [She1997], the classical quad-double library is that of Hida, Li and Bailey [HLB2001], and tight error bounds for double-word operations were proved by Joldes, Muller and Popescu [JMP2017]. A comprehensive treatment is in [MBD2018, chapters 4 and 14].

FLINT implements expansions with N = 1, \ldots, 4 in the dfloat module. Its kernels are designed for ball arithmetic: every operation computes the exact terms of the result with error-free transformations, keeps N of them, and returns the magnitude of everything it dropped, which becomes the radius of a ball.

Throughout, doubles are IEEE 754 binary64 numbers with round-to-nearest (\operatorname{RN}), u = 2^{-53} is the unit roundoff, and \operatorname{ulp}(x) is the weight of the last bit of the double x. For any real t such that \operatorname{RN}(t) is a normal double, \operatorname{RN}(t) = t(1 + \delta) with |\delta| \le u.

Error-free transformations#

The sum and the product of two doubles are not doubles in general, but their rounding errors are, and can be computed exactly.

Lemma 1 (Sterbenz). If a and b are doubles with b/2 \le a \le 2b, then a - b is a double, so it is computed exactly [Ste1974].

Algorithm 1: TwoSum and FastTwoSum

TwoSum(a, b) [Knu1997; Mol1965]: s \gets \operatorname{RN}(a + b), b' \gets \operatorname{RN}(s - a), e \gets \operatorname{RN}(\operatorname{RN}(a - \operatorname{RN}(s - b')) + \operatorname{RN}(b - b')).

FastTwoSum(a, b) [Dek1971]: s \gets \operatorname{RN}(a + b), e \gets \operatorname{RN}(b - \operatorname{RN}(s - a)).

Theorem 2. If a + b does not overflow, TwoSum returns s = \operatorname{RN}(a + b) and e with s + e = a + b exactly, using six operations and no branch. FastTwoSum gives the same result in three operations if the exponent of a is at least that of b (in particular if |a| \ge |b|) or a = 0. In both cases |e| \le \operatorname{ulp}(s)/2.

Proof. For FastTwoSum with |a| \ge |b|, we first show that z = s - a is computed exactly. If a and b have the same sign, then |a| \le |s| \le 2|a| with s of the sign of a, so this follows from Lemma 1. If they have opposite signs and |b| \ge |a|/2, then a + b is exact by the same lemma, so s = a + b and z = b. If they have opposite signs and |b| < |a|/2, then |a|/2 \le |s| \le |a| (rounding is monotone and |a|/2 is a double), and again Lemma 1 applies. Now b - z = (a + b) - s is the rounding error of a nearest-rounded sum of two doubles, which is itself a double (its magnitude is at most \operatorname{ulp}(s)/2 and it is a multiple of the ulp of the operand with the smaller exponent), so it is also computed exactly. TwoSum reconstructs the same quantities without knowing which operand is larger; see [MBD2018, §4.3] for the complete proof. Sums are exact in the subnormal range, so no underflow condition is needed.

Algorithm 2: TwoProd

TwoProd(a, b): p \gets \operatorname{RN}(ab), e \gets \operatorname{RN}(ab - p) computed with one fused multiply-add.

Theorem 3. If ab does not overflow and |ab| \ge 2^{-969} or ab = 0, TwoProd returns p + e = ab exactly.

Proof. The exact product ab has at most 106 significant bits and is a multiple of \operatorname{ulp}(a)\operatorname{ulp}(b). Its rounding error ab - p is at most \operatorname{ulp}(p)/2 in magnitude and is a multiple of the same quantity, so it fits in 53 bits; the condition on |ab| ensures that this multiple is not below the subnormal range. An FMA rounds the exact ab - p once, which is exact for a double.

Without an FMA, TwoProd can be computed with Dekker's product, which splits each operand into two 26-bit halves (Veltkamp splitting) and costs 17 operations [Dek1971]; _d_vec_dot_thrice uses this form. In C, all of these transformations break if the compiler contracts a separate multiplication and addition into an FMA, or evaluates intermediate results in extended precision (x87). The dfloat module is therefore compiled with -ffp-contract=off and with SSE2 arithmetic, and dfloat_is_supported reports whether these requirements hold.

Expansions and renormalisation#

An expansion (x_0, \ldots, x_{N-1}) is nonoverlapping if |x_{i+1}| \le \operatorname{ulp}(x_i)/2 for all i with x_i \ne 0, and canonical if in addition zeros occur only at the end. A canonical expansion is obtained from a real number x by rounding term by term, x_0 = \operatorname{RN}(x), x_1 = \operatorname{RN}(x - x_0), and so on; it is unique up to ties.

Arithmetic produces a list of terms, ordered by decreasing magnitude but overlapping, which must be converted to N nonoverlapping components.

Algorithm 3: VecSum (bottom-up renormalisation)

Input: doubles t_0, \ldots, t_{m-1}, roughly decreasing in magnitude.

Output: t_0', \ldots, t_{m-1}' with the same exact sum.

For i = m - 2, m - 3, \ldots, 0: (t_i, t_{i+1}) \gets \text{TwoSum}(t_i, t_{i+1}).

Since each step is error-free, the sum is preserved exactly, and after the pass each t_{i+1} is the rounding error of a sum whose rounded value is t_i, so |t_{i+1}| \le \operatorname{ulp}(t_i)/2. The result is not necessarily canonical: a cancellation in the sum that produced t_{i+1} can leave t_{i+2} larger than t_{i+1}. FLINT calls this weak normalisation. It is all that the error analysis needs, and it is what the dfloat kernels produce; a second, top-down pass of TwoSums (VecSumErr) followed by a check makes the expansion canonical (dN_canonicalise, or every result in a context with DFLOAT_STRONG). For N = 2 weak normalisation is already canonical. The general renormalisation _dfloat_renorm of an arbitrary number of terms adds a branchy pass (VecSumErrBranch of [JMPT2016]) that removes zero components before nonzero ones.

Addition and multiplication#

The kernels organise the exact terms of a result by order: a term of order k has magnitude about u^k times that of the head.

Algorithm 4: Addition of N-term expansions

For i = 0, \ldots, N-1: (s_i, e_i) \gets \text{TwoSum}(x_i, y_i).
S_0 \gets s_0. For k = 1, \ldots, N-1: accumulate S_k \gets s_k + e_{k-1} + (\text{residuals of order } k-1) with TwoSums, keeping their residuals as terms of order k.
Renormalise S_0, \ldots, S_{N-1} with Algorithm 3.
The dropped part is e_{N-1} plus the residuals of order N-1.

Algorithm 5: Multiplication of N-term expansions

For i + j < N: (h_{ij}, l_{ij}) \gets \text{TwoProd}(x_i, y_j).
S_0 \gets h_{00}. For k = 1, \ldots, N-1: accumulate with TwoSums S_k \gets \sum_{i+j=k} h_{ij} + \sum_{i+j=k-1} l_{ij} + (\text{residuals of order } k-1).
Renormalise S_0, \ldots, S_{N-1} with Algorithm 3.
The dropped part consists of the residuals of order N - 1, the low parts l_{ij} with i + j = N - 1, and the products x_i y_j with i + j \ge N, which are not computed exactly and are bounded by |\operatorname{RN}(x_i y_j)|.

The products need N(N+1)/2 TwoProds, each a multiplication and an FMA. Squaring uses the symmetry h_{ij} = h_{ji} (doubling is exact). The product with a single double needs N TwoProds.

Theorem 4. Let \mathit{res} be the renormalised result and E the sum of the absolute values of the dropped terms of Algorithm 4 or Algorithm 5. If no product of nonzero components has magnitude below 2^{-969} and nothing overflows, then |x \circ y - \mathit{res}| \le E, where \circ is the exact operation.

Proof. Every step is an error-free transformation by Theorem 2 and Theorem 3, so x \circ y equals the sum of the kept and the dropped terms exactly, except for the products of order \ge N, whose exact values satisfy |x_i y_j| \le |\operatorname{RN}(x_i y_j)|(1 + u); the factor 1 + u and the rounding errors in forming E itself are absorbed by the slack factor of Lemma 5.

For the plain formats, E is not computed, and the error is of the order of u^N times the magnitude of the result, with small constants [HLB2001; JMP2017]. In the ball formats, E is added to the radius. When the inputs are not exact, tracking the last-order residuals exactly costs as much as the product itself for N = 4, so the ball products then accumulate the last order with plain additions and bound their rounding errors by u times the accumulated magnitudes; this loses only a constant factor in a radius that is nonzero anyway, and keeps exact products of exact balls exact.

Radius arithmetic#

The radius of a ball is a double computed in round-to-nearest arithmetic from nonnegative quantities, and is made an upper bound by a single multiplication at the end.

Lemma 5. Let r be a nonnegative real number obtained from nonnegative doubles by at most 2^7 additions and multiplications, and \tilde r the result of performing the same operations in round-to-nearest arithmetic, where no operation has an inexact subnormal result. Then r \le \operatorname{RN}(\tilde r (1 + 2^{-45})).

Proof. Each operation on nonnegative numbers satisfies t \le \operatorname{RN}(t)/(1 - u), so by induction r \le \tilde r (1 - u)^{-2^7} \le \tilde r \exp(2^7 u (1 + u)) \le \tilde r (1 + 2^{-46} + 2^{-90}). On the other hand \operatorname{RN}(\tilde r (1 + 2^{-45})) \ge \tilde r (1 + 2^{-45})(1 - u) \ge \tilde r (1 + 2^{-45} - 2^{-52}), which is larger. Sums of doubles in the subnormal range are exact, and products and quotients that could underflow are bumped up by 2^{-1070} (_dfloat_rad_mul), which preserves the inequality.

The ball operations then use the standard propagation bounds. For the product, |x'y' - xy| \le |x| r_y + |y| r_x + r_x r_y for |x' - x| \le r_x and |y' - y| \le r_y, with |x| bounded by the sum of the absolute values of the components. For the quotient:

Proposition 6. Let q be any number and |y| > r_y. For all x' and y' with |x' - x| \le r_x and |y' - y| \le r_y,

\left|\frac{x'}{y'} - q\right| \le \frac{|x - qy| + r_x + |q| r_y}{|y| - r_y}.

Proof. x'/y' - q = (x' - qy')/y', with |x' - qy'| \le |x - qy| + |x' - x| + |q||y' - y| and |y'| \ge |y| - r_y.

The residual x - qy for the computed quotient q is enclosed with the tracked kernels, so that an exact quotient of exact balls gets radius zero. For the reciprocal square root, with e = 1 - x r^2 for the computed r, the series (1 - e)^{-1/2} = \sum_k \binom{2k}{k} (e/4)^k has positive coefficients and \sum_{k \ge 2} \binom{2k}{k} 4^{-k} |e|^{k-2} \le 1/2 for |e| \le 1/4, so that |x^{-1/2} - r| \le |r|(|e|/2)(1 + |e|).

Division and square root#

Algorithm 6: Long division with tapered remainders

Input: expansions x, y with y_0 \ne 0.

Output: q = q_0 + \cdots + q_{N-1} \approx x/y.

r \gets x, v \gets \operatorname{RN}(1/y_0).
For k = 0, \ldots, N-1:
1. q_k \gets \operatorname{RN}(r_0 v) (for the last digit, \operatorname{RN}(r_0/y_0)).
2. If k < N - 1: (p, e) \gets \text{TwoProd}(q_k, y_0), d_0 \gets r_0 - p (exact), and r \gets (d_0 - e) + (r_1 + \cdots) - q_k (y_1 + \cdots), kept with N - k - 1 components.
Renormalise q_0, \ldots, q_{N-1}.

The digit q_k is within a relative 2u + u^2 of r_0/y_0, so that p = \operatorname{RN}(q_k y_0) is within a factor 2 of r_0 and d_0 = r_0 - p is exact by Lemma 1: the cancellation of the heads loses no information, and the next remainder absorbs the error of the digit. After step k the remainder is about u^{k+1}|x| and is only needed to absolute precision u^N |x|, so it is kept with fewer components (tapered precision); the remainder update is a three-term sum of expansions of the size of the new remainder. A Newton–Karp–Markstein division [KM1997] (a reciprocal at half precision, one product and a correction) turned out to be slower at every N, because the steps of the long division are cheap products of a double with the divisor and the divisions are off the critical path.

The square root is the Newton iteration s \gets s + (x - s^2)/(2s) from s = \operatorname{RN}(\sqrt{x_0}), doubling the number of correct components per step. The residual x - s^2 is computed as (x_0 - s_0^2) - (s - s_0)(s + s_0) + (x_1 + \cdots) with the head exact (TwoProd and Lemma 1), and the correction is a long division producing only the digits it needs. The reciprocal square root iterates r \gets r + r(1 - x r^2)/2 in the same way.

Use in FLINT#

dfloat: the plain formats dN use Algorithm 4, Algorithm 5, Algorithm 6 and the Newton iterations; the ball formats dNb add the dropped terms of Theorem 4 to the radius with Lemma 5, and use Proposition 6 for dNb_div and the bound above for dNb_rsqrt. The kernels are written once (src/dfloat/kernels.inc) and instantiated for doubles and for four-lane SIMD vectors, which the vector functions _dN_vec_mul, _dN_vec_dot, etc. use. The error-free transformations are exposed as the macros DFLOAT_TWO_SUM, DFLOAT_FAST_TWO_SUM and DFLOAT_TWO_PROD.
The elementary functions of dfloat are built on these kernels; see Elementary functions in double-word to quad-double precision.
_dfloat_platt_smk computes the moments of Platt's multi-evaluation in plain double-double arithmetic with the algorithms DWTimesDW1 and AccurateDWPlusDW of [JMP2017], using their proven error bounds a priori.
_d_vec_dot_thrice in d_vec is the compensated dot product of Ogita, Rump and Oishi [ORO2005], built on TwoSum and Dekker's product.

void d2_mul(d2_t res, d2_srcptr x, d2_srcptr y)

void d2b_mul(d2b_t res, d2b_srcptr x, d2b_srcptr y)

void d2_div(d2_t res, d2_srcptr x, d2_srcptr y)

void _d2_mul_tracked(double * res, double * err, const double * x, const double * y)

void _dfloat_renorm(double * res, int n, double * err, double * t, int len)

void d2_canonicalise(d2_t res, d2_srcptr x)

References#

[Mol1965] O. Møller. Quasi double-precision in floating point addition. BIT 5 (1965), 37–50.
[Dek1971] T. J. Dekker. A floating-point technique for extending the available precision. Numerische Mathematik 18 (1971), 224–242. https://doi.org/10.1007/BF01397083
[Pri1991] D. M. Priest. Algorithms for arbitrary precision floating point arithmetic. Proceedings of the 10th IEEE Symposium on Computer Arithmetic (ARITH-10), 1991, 132–143.
[She1997] J. R. Shewchuk. Adaptive precision floating-point arithmetic and fast robust geometric predicates. Discrete & Computational Geometry 18 (1997), 305–363. https://doi.org/10.1007/PL00009321
[HLB2001] Y. Hida, X. S. Li and D. H. Bailey. Algorithms for quad-double precision floating point arithmetic. Proceedings of the 15th IEEE Symposium on Computer Arithmetic (ARITH-15), 2001, 155–162.
[JMP2017] M. Joldes, J.-M. Muller and V. Popescu. Tight and rigorous error bounds for basic building blocks of double-word arithmetic. ACM Transactions on Mathematical Software 44:2 (2017), 15:1–15:27. https://doi.org/10.1145/3121432
[MBD2018] J.-M. Muller, N. Brunie, F. de Dinechin, C.-P. Jeannerod, M. Joldes, V. Lefèvre, G. Melquiond, N. Revol and S. Torres, Handbook of Floating-Point Arithmetic, second edition, Birkhäuser (2018), https://doi.org/10.1007/978-3-319-76526-6
[Ste1974] P. H. Sterbenz. Floating-Point Computation. Prentice-Hall, 1974.
[Knu1997] Knuth, D. E. The Art of Computer Programming, volume 2: Seminumerical algorithms, 1997
[JMPT2016] M. Joldes, O. Marty, J.-M. Muller and V. Popescu. Arithmetic algorithms for extended precision using floating-point expansions. IEEE Transactions on Computers 65:4 (2016), 1197–1210. https://doi.org/10.1109/TC.2015.2441714
[KM1997] A. H. Karp and P. Markstein, "High-precision division and square root", ACM Transactions on Mathematical Software 23:4 (1997), 561–589
[ORO2005] T. Ogita, S. M. Rump and S. Oishi. Accurate sum and dot product. SIAM Journal on Scientific Computing 26:6 (2005), 1955–1988. https://doi.org/10.1137/030601818

Cite this page

The FLINT team, “Double-word and multi-double arithmetic,” FLINT: Fast Library for Number Theory, documentation, version 3.7.0-dev, 2026. https://flintlib.org/doc/algorithms/floating-point/double-word.html (generated 2026-10-20).

@misc{flint-doc-algorithms-floating-point-double-word,
  author       = {{The FLINT team}},
  title        = {{Double-word and multi-double arithmetic}},
  howpublished = {FLINT: Fast Library for Number Theory, documentation, version 3.7.0-dev},
  year         = {2026},
  url          = {https://flintlib.org/doc/algorithms/floating-point/double-word.html},
  note         = {Generated 2026-10-20}
}

Edit this page · PDF version · Built with fdoc from f12f9f517b