Emulating Double Precision in GLSL with Float-Float

Emulating 64-bit double precision on graphics hardware is essential when deploying compute tasks or high-zoom rendering (such as fractals) to platforms lacking native double support, including WebGL 1.0, WebGL 2.0, and many mobile GPUs. This article explains how double-single arithmetic (often called float-float or DSfun) emulates near-double precision in GLSL by storing numbers as an unevaluated sum of two standard 32-bit floats. We cover the underlying data representation, error-free transformations like Knuth's twoSum, and the fundamental algorithms for addition and multiplication.

The Float-Float Representation

In GLSL, an emulated double is represented as a vec2, where the components define a single real number \(x\):

\[x = x_{\text{hi}} + x_{\text{lo}}\]

Here, \(x_{\text{hi}}\) represents the most significant 24 bits of the mantissa, while \(x_{\text{lo}}\) captures the residual rounding error, effectively adding another 24 bits of precision. The condition \(|x_{\text{lo}}| \le \frac{1}{2}\text{ulp}(x_{\text{hi}})\) ensures that the two floats do not overlap in significant bits.

This approach provides approximately 46 to 48 bits of significand precision—closely approaching IEEE 754 double precision (53 bits). However, the exponent range remains restricted to that of 32-bit floats (\(2^{-126}\) to \(2^{127}\)), meaning this technique addresses precision loss, not underflow/overflow bounds.

Error-Free Transformations

The foundation of float-float arithmetic rests on error-free transformations (EFTs). These mathematical identities compute both the floating-point result of an operation and the exact rounding error discarded by the hardware.

QuickTwoSum

When it is known that \(|a| \ge |b|\), the error can be extracted using three floating-point additions/subtractions:

vec2 quickTwoSum(float a, float b) {
    float s = a + b;
    float e = b - (s - a);
    return vec2(s, e);
}

TwoSum (Knuth's Algorithm)

When the relative magnitudes of \(a\) and \(b\) are unknown, Knuth’s six-operation algorithm calculates the sum and exact rounding error:

vec2 twoSum(float a, float b) {
    float s = a + b;
    float v = s - a;
    float e = (a - (s - v)) + (b - v);
    return vec2(s, e);
}

Addition of Float-Float Pairs

Adding two emulated pairs, \(A = (a_{\text{hi}}, a_{\text{lo}})\) and \(B = (b_{\text{hi}}, b_{\text{lo}})\), requires accumulating the high components, the low components, and the rounding error produced when adding the high components:

vec2 ds_add(vec2 a, vec2 b) {
    vec2 s = twoSum(a.x, b.x);
    s.y += a.y + b.y;
    return quickTwoSum(s.x, s.y);
}

In this implementation, twoSum(a.x, b.x) splits the dominant sum into primary result s.x and truncation error s.y. The smaller terms (a.y, b.y, and s.y) are then combined into the low word before a final quickTwoSum normalizes the output.

Multiplication of Float-Float Pairs

Multiplication requires computing the cross-terms and the exact product of the high components. The product of two single-precision numbers \(a \times b\) produces up to 48 bits, which can be decomposed into an exact head and tail using twoProd.

If GLSL 4.0 or the fma() intrinsic is supported, twoProd is calculated using fused multiply-add:

vec2 twoProdFMA(float a, float b) {
    float p = a * b;
    float e = fma(a, b, -p);
    return vec2(p, e);
}

On hardware lacking fma(), Dekker's splitting algorithm is used. The input floats are split into two 12-bit halves to perform exact multiplication without hardware overflow:

vec2 split(float a) {
    const float c = 4097.0; // 2^12 + 1
    float p = a * c;
    float a_hi = p - (p - a);
    float a_lo = a - a_hi;
    return vec2(a_hi, a_lo);
}

vec2 twoProdDekker(float a, float b) {
    float p = a * b;
    vec2 a_split = split(a);
    vec2 b_split = split(b);
    float err = ((a_split.x * b_split.x - p) 
               + a_split.x * b_split.y 
               + a_split.y * b_split.x) 
               + a_split.y * b_split.y;
    return vec2(p, err);
}

Using twoProd, the full float-float multiplication (\(A \times B\)) is computed as:

vec2 ds_mul(vec2 a, vec2 b) {
    vec2 p = twoProdFMA(a.x, b.x); // or twoProdDekker
    p.y += a.x * b.y + a.y * b.x;
    return quickTwoSum(p.x, p.y);
}

Performance Characteristics

Float-float emulation carries an instruction overhead of approximately \(3\times\) to \(8\times\) compared to native single-precision math. However, because it relies strictly on primary 32-bit ALUs without requiring complex control flow or memory operations, it parallelizes cleanly across GPU warps/wavefronts and runs significantly faster than multi-pass rendering solutions or CPU readbacks.