Error-Free Transformations and Double-Double Arithmetic

The previous articles treat rounding as a loss. It is better understood as a partition: a correctly rounded operation returns the leading part of the exact result and discards the trailing part, and the trailing part is itself a floating-point number. If both parts are kept, the operation has lost nothing. The algorithms that compute the second part are called error-free transformations, and the arithmetic built on them is the standard way to obtain more than the hardware’s precision without a wider hardware format.

The residual of an operation is representable

The foundation is a theorem about rounding rather than about any particular algorithm. If s = fl(a + b) and no overflow occurs, then a + b − s is exactly representable in the same format. The reason is a counting argument on the significands: the error of a correctly rounded sum is the tail that was shifted out during alignment, and that tail is a multiple of the smallest spacing involved and shorter than a full significand.

The same holds for multiplication, with a condition: a*b − fl(a*b) is exactly representable provided the product does not underflow. In the subnormal range the fixed spacing can make the residual unrepresentable, and the guarantee evaporates.

Given representability, an error-free transformation is a fixed sequence of ordinary floating-point operations that returns the pair. It costs a handful of operations and requires no access to the rounding mode, no wider type, and no library support.

TwoSum and the cost of not knowing the order

Knuth’s TwoSum computes s = fl(a + b) and the exact error e with six operations:

/* s + e == a + b exactly, for any a and b, provided no overflow. */
static void two_sum(double a, double b, double *s, double *e)
{
    *s = a + b;
    double bb = *s - a;              /* b as seen through the rounded sum */
    double err = (a - (*s - bb)) + (b - bb);
    *e = err;
}

Each subtraction in the sequence is exact, and the composition reconstructs the part of b that the addition discarded. The algorithm makes no assumption about the relative magnitudes of the operands, which is why it costs six operations rather than three.

Dekker’s FastTwoSum is the three-operation version, and it requires |a| ≥ |b|:

static void fast_two_sum(double a, double b, double *s, double *e)
{
    *s = a + b;
    *e = b - (*s - a);   /* exact only when |a| >= |b| */
}

When the magnitude condition holds, the intermediate *s - a recovers b’s contribution exactly, and one subtraction gives the residual. When it does not, the recovered quantity is the wrong one and the transformation silently returns a residual for a different pair. The condition is not a performance hint; it is the correctness requirement, and code that sorts operands to satisfy it pays a comparison per operation.

The difference in cost matters in a loop. Six operations per term is a factor of six over a plain addition, and both variants are serially dependent: each iteration’s sum feeds the next. A dot product computed with TwoSum for the additions and TwoProduct for the products costs roughly twenty operations per term, which is why compensated summation with a single correction term is preferred when the full exactness of an expansion is not required.

Splitting, and the product without FMA

Multiplication is harder because the residual is not obtainable from the product alone. Dekker’s approach splits each operand into a high part with half the significand and a low part that holds the rest, multiplies the four pieces, and reassembles.

The split uses a multiplier of 2^⌈p/2⌉ + 1: 2^27 + 1 = 134,217,729 for binary64, 2^12 + 1 = 4097 for binary32. Multiplying by that constant forces the significand to break at bit 26, and two subtractions separate the pieces.

#include <math.h>

#define SPLIT 134217729.0   /* 2^27 + 1 for binary64 */

static void split(double a, double *hi, double *lo)
{
    double t = SPLIT * a;
    *hi = t - (t - a);
    *lo = a - *hi;
}

/* p + e == a * b exactly, without a fused multiply-add. */
static void two_product_dekker(double a, double b, double *p, double *e)
{
    double ah, al, bh, bl;
    split(a, &ah, &al);
    split(b, &bh, &bl);
    *p = a * b;
    *e = ((ah * bh - *p) + ah * bl + al * bh) + al * bl;
}

Seventeen operations, versus two for the fused version shown in Fused Multiply-Add, Single Rounding, and the Cost of Contraction. The split has a range limit as well: SPLIT * a overflows when a approaches 2^997 in binary64, so the transformation is not total. Values near the top of the exponent range need a scaled variant.

What a double-double pair represents

A double-double stores a value as hi + lo, where hi is a binary64 and lo is the residual of some earlier operation. The pair carries about 106 bits of significand — roughly 31 or 32 decimal digits — and its arithmetic is built from the transformations.

Addition of two pairs uses TwoSum on the high parts, adds both low parts, and renormalises so that the high part is the larger magnitude:

typedef struct { double hi, lo; } dd;

static dd dd_add(dd a, dd b)
{
    double s, e;
    two_sum(a.hi, b.hi, &s, &e);
    e += a.lo + b.lo;
    dd r;
    r.hi = s + e;          /* renormalise: hi is the rounded sum */
    r.lo = e - (r.hi - s); /* and lo is what that rounding dropped */
    return r;
}

The renormalisation step is what keeps the representation canonical: after it, |lo| ≤ ulp(hi)/2, so the pair is a non-overlapping expansion of two components. Without it, the components drift and the effective precision degrades.

Multiplication follows the same pattern with TwoProduct and four products: hi*hi exactly, plus the three cross terms, each accumulated with TwoSum.

Three properties are worth stating precisely because they are commonly overstated.

The precision is not uniform. A double-double has 106 bits of significand near its high part, but the spacing between adjacent representable pairs is not constant, and there are gaps in the representation that a true 106-bit binary format would not have. The effective precision is often quoted as “about 106 bits”, which is an average, not a guarantee for every value.

The arithmetic is not correctly rounded. A double-double operation has a relative error of order u² — about 1.2e-32 — rather than being the correctly rounded result of the exact operation. Obtaining a correctly rounded 106-bit result requires an additional step, and libraries that need it implement it separately.

The operations are not associative and not closed under all inputs. Cancellation between the high parts beyond the capacity of the low part loses the guarantee, and mixing double-double values with plain doubles requires converting, which is exact only in one direction.

Where the transformations break

Four failure modes account for most of the ways this machinery goes wrong.

The first is compiler optimisation. Every transformation in this article is a sequence of operations whose exact-arithmetic value is zero or a simple expression; a compiler entitled to reassociate or to simplify algebraic identities will delete the steps that carry the information. The TwoSum residual is (a − (s − bb)) + (b − bb), which simplifies to zero under real arithmetic. Compiling with flags that permit reassociation destroys it, and the resulting code is a plain, slightly more expensive sum.

The second is the rounding mode. The derivations assume round-to-nearest. Under a directed mode the residual is still exact in most cases, but the renormalisation steps that rely on the tie behaviour of nearest rounding no longer produce a canonical pair, and the accumulated error bounds change.

The third is underflow. A product that falls into the subnormal range has a residual that may not be representable, and a sum whose result is subnormal has an error that is a multiple of the subnormal spacing but may still be representable — the guarantee is weaker and the bounds acquire absolute terms. Code that runs at the bottom of the exponent range needs a scaled formulation.

The fourth is the assumption that hi dominates. The renormalisation in dd_add assumes the incoming pairs are already normalised. A pair constructed by hand, or one that has passed through an operation that failed to renormalise, can have |lo| > |hi|/2, and every subsequent operation then produces a result whose stated error bound does not apply.

What the extra precision is for

The practical value is not in replacing binary128. It is in making a small number of critical operations exact while the bulk of the computation stays in binary64: a dot product that must not lose the cancellation, a residual computation in iterative refinement, an accumulation of many small terms where the condition number is large. The Dot2 algorithm of Ogita, Rump, and Oishi is the canonical example — a dot product built from one TwoProduct and one TwoSum per term, whose error bound does not grow with n at first order, as a naive accumulation’s does.

That is the same guarantee an exact accumulator provides, at a fraction of the cost, and it is the reason the transformations survive in libraries despite the compiler hazards. They are also the reference implementation against which correctly rounded conversions are tested, since an exact computation in double-double is a convenient oracle for the last bit of a decimal parse. The limits of what any of this can guarantee once the compiler is allowed to reassociate are the subject of Reassociation, Excess Precision, and the Limits of Bitwise Reproducibility.