GoKawiilIf you ever need more precision than what 15 decimal digits of the double format can offer, 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 puts 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.
The gap nobody fills#
Floating-point types come in many "flavors". For example, a float has ~7 decimal digits, a double ~15, both basically for free. An arbitrary-precision library gives you as many digits as you want, at a painful per-operation cost. 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. The image at the top of this page shows this gap: the same view rendered twice, blocky on the left where a double has run out of precision, sharp on the right rendered with double-double.
The Mandelbrot set is a famous fractal, full of endlessly repeating spirals and mini copies of itself, that I have explored before. A deep zoom literally runs out of precision. Eventually, two neighboring pixel positions are "rounded" to 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 thus are colored the same. The gap between neighboring doubles grows with the number's magnitude, 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).
The canonical answer to "I need more precision than double" is to use a library for arbitrary-precision math, typically GMP or MPFR. That is the right answer when the precision you need is open-ended. But it is not a small leap to take, and you pay for it on every single operation:
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.
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. Every operation is a loop. Add or multiply, everything walks the limb array: loads, stores, and data-dependent branches.
Add or multiply, everything walks the limb array: loads, stores, and data-dependent branches. 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.
... continue reading