Double-double: 31 digits of precision without leaving the FPU | Marek's blog

Skip to main content

If you ever need more precision than what 15 decimal digits of the double format can offer, and arbitrary-precision arithmetic is too slow, too heavy, or simply not available, there is a neat trick: glue two doubles together and treat them as one number. This gives you ~31 decimal digits for roughly 9x the cost of a plain double in a real kernel (4-12x per isolated operation). With no heap allocation and no dependencies, this lands it almost exactly halfway between a double and an arbitrary-precision library perf-wise. This post explains the error-free transformations that make it work, measures it against MPFR, and shows where the trick runs out of steam.<br>The gap nobody fills<br>Floating-point precision comes in many sizes.<br>For example, a double gives you about 15 decimal digits, essentially for free.<br>An arbitrary-precision library gives you as many digits as you want, at a painful per-operation cost.<br>Between "not quite enough" and "orders of magnitude slower" there is a gap and I fell into it while zooming deep into the Mandelbrot set.<br>The image at the top of this page is that gap in one picture: the same view rendered twice, blocky on the left where a plain double has run out, sharp on the right where a double-double has not.<br>The Mandelbrot set is a famous fractal with these mesmerizing patterns that I have explored before.<br>A deep zoom runs out of precision in the most literal way possible: eventually two neighboring pixels land on the same double, and the image stops being a picture of the fractal and starts being a picture of the number format, as you can see in Figure 1.

Figure 1: Two deep locations rendered in double precision: whole blocks of pixels share one coordinate and come out the same color. The gap between neighboring doubles grows with the number's size, so the real axis is the coarse one on the left (re -0.74, im 0.13, wide blocks) and the imaginary axis is the coarse one on the right (re -0.11, im 0.92, tall blocks).<br>The canonical answer to "I need more precision than double" is to use a library for arbitrary-precision math, typically GMP or MPFR.<br>That is genuinely the right answer when the precision you need is open-ended.<br>But it is not a small leap to take, and you pay for it on every single operation:<br>Every value is a heap block. A value is a pointer to a data array (limbs), so bringing one into existence allocates. An API where each operation returns a new value, which is what many wrappers offer, allocates on every operation.<br>Every operation is a loop. Add or multiply, everything walks the limb array: loads, stores, and data-dependent branches.<br>In .NET it is worse still. GMP and MPFR are native C libraries, so you are looking at P/Invoke, a native binary per platform to ship, and marshalling on the boundary.<br>And then there is the license: GMP and MPFR are LGPL, which static linking, closed platforms, or company policy can turn into a hard no, regardless of benchmark results.<br>I benchmarked all of that, and the gap turned out to be surprisingly well structured.<br>In a real kernel at matched precision, a double-double costs roughly 9x a plain double, and MPFR carrying those same 31 digits costs roughly 9x a double-double again, which is about 81x the double it started from.<br>Let MPFR allocate a result per operation and that double-double => MPFR gap is closer to 40x instead of 9x.<br>The insight needed to unlock double-double precision is:<br>Keep two doubles separate, but treat them as one number.<br>About 31 decimal digits, no heap, no dependencies, and no native binary to ship.<br>The complete type and the benchmarked kernels are in the two appendices.<br>The double-double idea<br>A floating-point type holds the same number of significant digits no matter how large the number is.<br>The exponent decides where the digits sit, it does not change how many there are.<br>So take these two numbers:<br>(1)A = 111222333444<br>(2)B = 0.555666777888<br>Each of them has 12 significant digits, so each fits in a double with room to spare.<br>Their exact sum, 111222333444.555666777888, has 24 digits, and a double holds about 15: evaluate A + B and the tail of B is quietly rounded away to 111222333444.55566.<br>But what if we simply do not evaluate the addition?<br>If we keep A and B side by side and agree to treat the pair as one number, then one of them carries the leading digits and the other carries digits the first one has no room for.<br>Nothing was rounded, because nothing was added.<br>That is the whole representation:<br>A double-double value stores the unevaluated sum of two doubles .<br>(3)x = xhi + xlo<br>with the invariant that xhi is exactly what xhi + xlo rounds to as a double.<br>The high part carries the value, the low part carries the error of the high part.<br>Keeping that invariant true is what every operation in the next section is really doing.<br>Code listing 1: The whole double-double representation. Everything else in this post is arithmetic...