Double-double: 31 digits of precision without leaving the FPU
If 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. With no heap allocation and no dependencies, this puts it perf-wise almost exactly halfway between a double and an arbitrary-precision library. This post explains how the format works and the math behind it. It also benchmarks the type against __float128 and MPFR, and shows where the trick runs out of steam.
1. The gap nobody fills#
Floating-point types come in many "flavors". For example, a float has ~7 decimal digits, a double ~15. An arbitrary-precision library gives you as many digits as you want but at a painful per-operation cost (~100x slower or more). Between "not quite enough precision" 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 was my motivation. It's the Mandelbrot set rendered using double precision (left) and double-double precision (right).
The Mandelbrot set is a famous fractal, full of endlessly repeating spirals and mini copies of itself, that I have explored before.
When zooming deeper and deeper, eventually two neighboring pixel positions are "rounded" to the same double value,
and the image stops being a picture of the fractal and starts being a picture of the number format.
Notice that the "ugly" blocks in Figure 1 are oriented differently in each image.
This is a consequence of how a floating point representation works.
The gap between neighboring numbers grows with the magnitude, so the larger of the two coordinates loses resolution first.
The left image is at (-0.74, 0.13) so the blocks are wider than tall, but the right image is at (-0.11, 0.92) where the ratio is reversed.
The canonical answer to "I need more precision than double" is to use a library for arbitrary-precision math, typically GMP or MPFR. But it is not a small leap to take, and you pay for it on every single operation:
- Every value is a heap block. Each value is a pointer to a data array (often called limbs), so creating a new number allocates. One has to be very careful to reuse values to avoid allocations on every operation.
- Every operation is a loop. 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.
So, to get a little more precision cheaply, the trick is:
Keep two doubles separate, but treat them as one number.
This gives you ~31 decimal digits with no heap allocations, no dependencies, and no native binary to ship.
The idea is not mine. Dekker described it in 1971, and the pieces it is built from are even older (see Further reading for details). I wanted to prove to myself (and to you) that this trick is still relevant today, after half a century of faster hardware and better arbitrary-precision libraries.
2. The double-double idea#
A floating-point type holds the same number of significant digits no matter how large or small the number is. The exponent then determines where the decimal point will be. Take these two numbers:
A = 111222333444B = 0.555666777888Each of them has 12 significant digits, so each fits in a double.
Their exact sum is 111222333444.555666777888 and has 24 digits, but a double can only hold around 15.
If you evaluate A + B, the digits that do not fit are rounded away, leaving you with just 111222333444.55566.
But what if we do not evaluate the addition? 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 the rest.
That is the whole representation:
A double-double value stores the unevaluated sum of two doubles.
x = xhi + xlowith the invariant that the high part xhi is exactly what xhi + xlo rounds to as a double.
The low part xlo carries the rest (the "error" of the high part after addition).
Every operation in the next section keeps this invariant true.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
public readonly partial struct DoubleDouble : IEquatable<DoubleDouble>, IComparable<DoubleDouble> {
private readonly double m_hi;
private readonly double m_lo;
/// <summary>
/// The leading component; equals the value rounded to nearest double.
/// </summary>
public double Hi => m_hi;
/// <summary>
/// The trailing component (the rounding error of <see cref="Hi"/>).
/// </summary>
public double Lo => m_lo;
}
128 bits, but not quadruple precision. The 128-bit type in IEEE 754 (the standard that also defines the double) is called binary128, better known as quadruple-precision.
It has a 113-bit significand and a 15-bit exponent (that reaches to ~1e4932).
A double-double also uses 128 bits of storage, but it is just two ordinary doubles next to each other, resulting in 106 combined significand bits and an unchanged exponent range (still overflowing at ~1e308).[1]
Figure 2: Number of significand bits of each format.3. Error-free transformations#
Everything rests on one insight:
When a double-arithmetic operation rounds, the amount it rounded away can still be recovered as another double.
When you add two doubles, the hardware computes the true sum, which is then rounded.
Let the rounded result be s = round(a + b).
The part that got rounded away is exactly e = (a + b) - s.
It may seem that this formula is self-defeating, it needs the exact sum but the exact sum is the one thing we cannot have.
And yet, it is possible to compute the error term e exactly.
And once we have it, we are done because a + b = s + e, and (s, e) is our double-double!
So how to compute e exactly? The algorithm is called Knuth's two-sum and you can see it in Code listing 2: six additions/subtractions, no branches, no bit twiddling, and no wider types needed.
1
2
3
4
5
6
7
8
private static (double Sum, double Error) twoSum(double a, double b) {
double sum = a + b;
double bKept = sum - a; // the part of b that made it into the sum
double aKept = sum - bKept; // and the part of a
double aLost = a - aKept;
double bLost = b - bKept;
return (sum, aLost + bLost);
}
Here is the algorithm on our example numbers A and B:
1
2
3
4
5
6
7
8
9
10
11
12
13
A = 111222333444
B = 0.555666777888
A + B, exactly = 111222333444.555666777888
sum = A + B = 111222333444.5556640625 <- what the FPU hands back
bKept = sum - A = 0.5556640625
aKept = sum - bKept = 111222333444
aLost = A - aKept = 0
bLost = B - bKept = 0.000002715388
error = aLost + bLost = 0.000002715388
sum + error = 111222333444.555666777888 <- exact, to the last digit
The double sum kept the leading 17 digits of the exact one and dropped the rest, exactly as expected.[2]
The line to look at is bLost. It is not an approximation of what got dropped, it is what got dropped.
In the extreme case, when the operands are too "far apart" in magnitude to overlap at all, the sum rounds to the larger one and the error is the smaller one.
Notice that A lost nothing in the example above (aLost is 0) and that is not a coincidence.
Whenever |a| >= |b| the sum keeps every bit of the larger operand, so aKept and aLost can be dropped from the function entirely.
That is Dekker's fast two-sum:
1
2
3
4
5
private static (double Sum, double Error) quickTwoSum(double a, double b) {
double sum = a + b;
double bKept = sum - a;
return (sum, b - bKept); // a kept everything, so only b's loss is left
}
Multiplication needs the same thing. The exact product of two 53-bit numbers needs 106 bits, so it splits into a high and a low part. Dekker's original method (Veltkamp splitting) took 17 operations. Once the hardware has an FMA (fused multiply-add), which every mainstream processor does, it takes just two:
1
2
3
4
5
private static (double Product, double Error) twoProduct(double a, double b) {
double product = a * b;
double error = Math.FusedMultiplyAdd(a, b, -product);
return (product, error);
}
The trick is that an FMA computes a*b + c with a single rounding at the very end.
Call it as fma(a, b, -product) and it evaluates a*b - product in one go.
The subtraction cancels everything product kept, leaving exactly the error that a * b threw away, small enough to be recovered exactly.
These three functions are called error-free transformations, and they are the entire foundation. They all assume round-to-nearest arithmetic and no overflow. Note that infinities and NaNs are not handled properly.
4. Building arithmetic operations for double-double#
With our twoSum, quickTwoSum, and twoProduct in hand, writing arithmetic operations for DoubleDouble is quite straightforward.
Just remember to use an error-free transformation wherever a rounding would lose digits we care about
and finish with a quickTwoSum so the result is a normalized pair.
Addition#
Add the two components pairwise, then incorporate lowSum and lowError one by one.
Addition of the lowError itself can round, so we need to use an extra quickTwoSum.
1
2
3
4
5
6
7
8
9
public DoubleDouble Add(DoubleDouble right) {
(double sum, double error) = twoSum(m_hi, right.m_hi);
(double lowSum, double lowError) = twoSum(m_lo, right.m_lo);
error += lowSum;
(sum, error) = quickTwoSum(sum, error);
error += lowError;
(sum, error) = quickTwoSum(sum, error);
return new DoubleDouble(sum, error);
}
There is a cheaper "sloppy" variant in circulation using only 11 operations instead of these 20,
which adds both low parts m_lo and right.m_lo in one step.
It is fine while the operands share a sign, but if they have opposite signs and cancel each other,
the relative error has no upper bound.
Multiplication#
Multiply the two high parts exactly with twoProduct, then add the three cross terms to the error.
1
2
3
4
5
6
public DoubleDouble Multiply(DoubleDouble right) {
(double product, double error) = twoProduct(m_hi, right.m_hi);
error += m_hi * right.m_lo + m_lo * right.m_hi + m_lo * right.m_lo;
(product, error) = quickTwoSum(product, error);
return new DoubleDouble(product, error);
}
The one exact product of the high halves carries the value.
The cross terms are corrections, and they do not need to be exact.
Each of them is smaller than the product it corrects by a factor of at least 2-53,
so whatever plain double arithmetic rounds away from them is at 2-106 and below,
under the precision floor of DoubleDouble.
Squaring can be optimized by using just two cross terms instead of three,
and multiplying by a plain double needs only one.
These are worth having as overloads for better performance.
Division#
There is no error-free transformation for division. Instead, take a double's worth of the quotient, subtract its contribution at full double-double width, and repeat as long division three times.
1
2
3
4
5
6
7
8
9
10
public DoubleDouble Divide(DoubleDouble right) {
double quotient1 = m_hi / right.m_hi;
DoubleDouble remainder = Subtract(right.Multiply(quotient1));
double quotient2 = remainder.m_hi / right.m_hi;
remainder = remainder.Subtract(right.Multiply(quotient2));
double quotient3 = remainder.m_hi / right.m_hi;
(double sum, double error) = quickTwoSum(quotient1, quotient2);
return new DoubleDouble(sum, error).Add(FromDouble(quotient3));
}
Square root uses the same idea. Take a double root, then one step of Newton refinement. Division and square root come out accurate to a couple of units in the last place but not correctly rounded (not the nearest representable pair). The error starts to matter only when the last bits have to be reproducible against a different implementation.
Comparison, for free#
The comparison is easy thanks to the invariant.
Note that if unnormalized pairs were allowed, (1, 0.75) and (1.75, 0) would be the same number and compare as different ones.
1
2
3
4
public int CompareTo(DoubleDouble other) {
int hiComparison = m_hi.CompareTo(other.m_hi);
return hiComparison != 0 ? hiComparison : m_lo.CompareTo(other.m_lo);
}
One pair, four shapes#
The invariant leaves the two 53-bit significand windows surprising freedom in their relative positions, and the pair behaves a little differently in each arrangement as shown in Figure 3.
Figure 3: The arrangements a normalized pair can take, and the one it cannot.- Single double: Any value that a
doublecan hold exactly is a pair withlo = 0. - Touching:
locontinues exactly wherehiends, giving 106 contiguous bits. For example,1 + 2-53becomes(1, 2-53), exactly. - Gapped: This is the interesting case, where
lostarts "further down", and every bit position in between is an implicit zero. This feels like a special case, a "double-double denormal", but it is not. The pair from our example happens to have a gap of 2 (Code listing 10). - Overlapping never exists. If the windows overlapped,
hiwould no longer be the value rounded to a double and the same number could be written as many different pairs. That is what the invariant forbids, and you cannot even build it. For example,FromTwoDoubles(1.0, 0.75)renormalizes and returns(1.75, 0).
The gap of 2 from our example can be seen if we write the digits out in binary.
1
2
3
4
5
6
7
A = ...0000100. <- 37 integer bits, the top one at 2^36
B = .1000111001000000001011011000...
sum = ...0000100.1000111001000000 <- 53-bit window, ends at 2^-16
error = 001011011000... <- first set bit at 2^-19
^^
the gap: two zero bits stored by neither half
In our case, the error term happens to begin with two zeros, and that is all a gap is.
Those zeros are real digits of the value, but double never spends storage on leading zeros.
Instead, its exponent simply points lower (0.001011 and 1.011 ∙ 2-3 are the same number).
The gapped shape has one catch.
It may seem like you are getting extra precision "for free", but the zeros may collapse on any operation and any extra precision will be lost.
Suppose this absurd case: 10300 + 10-300 is a perfectly legal pair with 1,940 zero bits between its islands, carried in a 16-byte DoubleDouble.
However, evaluate (10300 + 10-300) + 1 and the exact answer would need a third island in the middle, so the smallest one is dropped and the answer is simply 10300 + 1.
"But surely 10300 written in binary does not have loads of trailing zeros‽"
Right, it does not. Writing it exactly needs about 700 significant bits.
But the 1e300 in the pair is not that number, it is the nearest double, which stores a 53-bit significand times a power of two, with every position below those 53 bits implicitly zero.
Double-double promises ~106 significant bits at the leading edge of the value.
A wide-gapped pair can exceed that promise but only temporarily, as any arithmetic operation can erase it.
Printing runs into the same limit from the other side, which is why the type also carries a ToStringExact.
And if you are wondering why not 3 doubles, or 4? It works too! 3 doubles make a triple-double (~159 bits), 4 a quad-double (~212 bits). But each extra term adds the same 53 bits while renormalization keeps getting more expensive, so somewhere around the third or fourth double the precision stops being worth the cost.
5. What it actually costs#
Chained single-operation latency, i9-13900K, .NET 10, BenchmarkDotNet:
| Operation | double | DoubleDouble | ratio |
|---|---|---|---|
| add | 0.325 ns | 4.07 ns | 12.5x |
| multiply | 0.64 ns | 3.08 ns | 4.8x |
| square | 0.64 ns | 2.83 ns | 4.4x |
| multiply + add | 1.17 ns | 7.08 ns | 6.0x |
| divide | 2.20 ns | 24.09 ns | 10.9x |
Addition is the expensive one (12.5x) and multiplication the cheap one (4.8x), the exact opposite of every other numeric type. The FMA does the whole exact product in one instruction, while the exact sum needs two twoSums and two renormalizations.
And division, which has no error-free transformation at all, costs what long division costs.
Against a real arbitrary-precision library#
Now let's compare DoubleDouble against MPFR, an arbitrary-precision library.
The code I am comparing is a Mandelbrot set kernel that iteratively computes zn in the complex plane.
zn+1 = zn2 + cIf |zn| stays bounded forever, c belongs to the set and the pixel is set to black.
If it "escapes" to infinity, which provably happens once |z| > 2, the iteration count at that point determines the color.
One image totals roughly 390 million iterations.
MPFR goes through hand-written P/Invoke with no wrapper class, at 106 bits of precision, which matches double-double's.
| Kernel | iterations/second | vs double | vs DoubleDouble |
|---|---|---|---|
double | 8,700 M | 1x | |
DoubleDouble | 1,000 M | ~9x | 1x |
| MPFR @106, in place | 108 M | ~81x | ~9x |
| MPFR @106, fresh result per operation | 26 M | ~330x | ~38x |
The two MPFR styles run the identical operations, but the "fresh result" style does not reuse values and allocates per operation, the way many wrappers present it.
A double-double costs roughly 9x a double, and MPFR costs roughly 9x a double-double again. There is no reason its perf should end up in the geometric middle of the gap it fills, but in my case it did.
Read those as "about 9" and "about 80" rather than to the digit.[3]
The slowdown in the "fresh result" row is avoidable, but it is up to the code writer to keep the discipline and properly reuse values (e.g. never write a * b + c as an expression because each operation would allocate a new value).
DoubleDouble has no such style to get wrong.
A value is two doubles in a readonly struct, so zRe * zIm * 2.0 + cIm allocates nothing, and no later edit can regress it.
The 9x is the number for my particular case, and it comes from 7 arithmetic operations per iteration, so your mileage may vary.
Compared to arbitrary-precision floats, the DoubleDouble type wins on four counts:
- Inlineable. A double-double add or multiply is just 7-20 instructions, which will often be inlined, unlike MPFR, which is a function call (with a P/Invoke boundary on top in .NET).
- No loops, no branches. A fixed branch-less instruction sequence will make your CPU pipeline happy, unlike MPFR's data-dependent loops over limbs.
- Runs on the FPU. Native
doubleoperations do the heavy lifting. - It can vectorize. A structure-of-arrays rewrite over
Vector256<double>does four values at a time, and the same arithmetic ports to a GPU. While this is a rewrite rather than something the compiler/JIT does for you, such a rewrite would not be possible with an arbitrary-precision library.
Against the compiler's 128-bit float#
If you write C/C++ on GCC/Clang, a quadruple-precision __float128 is built in.
Declare it and every operator just works.
If it were fast, this article would be a lot shorter.
The catch is that no x86 or mainstream ARM processor implements __float128 in hardware, so the compiler emulates every operation.
Measured with the same chained loops and the same kernel, this time in C[4]:
| Kernel (C, one core) | iterations/second | vs double | vs DoubleDouble |
|---|---|---|---|
double | 531 M | 1x | |
long double (x87) | 441 M | 1.2x | |
DoubleDouble | 76 M | ~7x | 1x |
__float128 | 12.4 M | ~43x | ~6x |
Per operation __float128 is 2-4x behind double-double, but the entire Mandelbrot kernel runs 6x slower.
The long double row shows what hardware support is worth: nearly a double's speed, but only 64 significand bits.
The 6x slowdown does buy you some things: 7 more bits of precision, far more exponent range, correct rounding, and properly working infinities/NaNs.
And where __float128 is in hardware, IBM's POWER9 and later, it is the right answer.
On x86 the trade goes the other way, and in .NET the question never arises, as there is unfortunately no 128-bit float to reach for (yet?).
6. Showcase: a Mandelbrot zoom that double cannot reach#
Back to the zoom that started this. Here is what the renderer is asked to draw:
1
2
3
4
5
private const string CENTER_RE = "-0.74364389112485980556186";
private const string CENTER_IM = "0.13182591316154390843318";
private const double MARQUEE_WIDTH = 7.746e-15;
private const double ESCAPE_SQUARED = 65536.0 * 65536.0;
The kernel is deliberately the same code once per type: same operations, same order, same values. Code listing 12 shows the double-double kernel, and the rest are in Appendix B: the benchmarked kernels.
1
2
3
4
5
6
7
8
9
10
11
12
13
DoubleDouble zRe = DoubleDouble.Zero;
DoubleDouble zIm = DoubleDouble.Zero;
for (; n < maxIterations; n++) {
DoubleDouble zReSquared = zRe.Square();
DoubleDouble zImSquared = zIm.Square();
magnitudeSquared = zReSquared.Hi + zImSquared.Hi;
if (magnitudeSquared > ESCAPE_SQUARED) {
break;
}
zIm = zRe * zIm * 2.0 + cIm;
zRe = zReSquared - zImSquared + cRe;
}
Figure 4: A view 1e-3 across. At this depth both types render the same image.
double 
double-double
The pixel step here is about 1.5e-17, a seventh of the gap between neighboring doubles, so the blocks are about seven pixels wide, and every further zoom step widens them until the whole frame is a single one. With double-double the blocks are gone, and the detail holds down to a view width of about 1e-28.
Keep going, though, and the same thing happens to it, as Figure 6 shows.
double-double 
perturbation, for reference
7. Having 31 digits is not the same as using 31 digits#
So far, all we measured was how finely you can represent a point. But how much precision the computation is left with is a different question, and the Julia set makes this distinction very clear.
A Julia set uses Equation 4 again, but with the two roles swapped: c is one constant for the whole image, and the pixel coordinate is the starting z0 instead of zero.
Same arithmetic, same cost per iteration, and the image falls apart at a zoom the Mandelbrot renderer would not even notice.
I ran into this while browsing Julia sets in my own renderer. I often found that the most interesting structure is in the center of the image (the origin), so that's where I zoomed in. The glitches arrived so early that I first thought I had a bug in my code. They also looked very different. However, switching the kernel to double-double eliminated them completely, so the bug was not in the code, and I went looking for what had run out.
double 
double-double
The coordinate is not what ran out. The center is 9.7e-8 from the origin, and doubles are dense down there, with consecutive ones just 1.3e-23 apart.
The precision runs out in the first line of the loop, in the addition rather than the squaring.
Squaring makes the pixel's coordinate smaller, but a double keeps all the important digits.
The damage is done at the very next step, where that tiny result is added to c, a number around 0.75.
I've written an example out so you can see what survives:
1
2
3
4
5
6
7
8
9
10
11
12
13
z_0 = 0.000000012 the pixel
z_0^2 = 0.000000000000000144 exact, to a double's full 16 digits
c = 0.75
c + z_0^2 = 0.750000000000000144 what the sum should be
c + z_0^2 = 0.750000000000000111 the nearest double, 23% short
^^^
next pixel over, 3.1e-10 away:
z_0 = 0.0000000123125 a different coordinate
z_0^2 = 0.000000000000000152 and a different square
c + z_0^2 = 0.750000000000000111 onto the very same double
Near 0.75 the doubles are only 1.1e-16 apart, and this pixel is worth 1.44e-16, so the sum can record it only as one whole step, storing the pixel's contribution 23% short.
The neighbor's square is 5% larger but still rounds to that same number after the addition.
Two pixels that started 3.1e-10 apart are indistinguishable before the second iteration begins.
You can see that neighborhood in Figure 7, in the bottom right of the first picture where every pixel falls to the same z1 = c, shown as a large gray region.
The double-double next to it draws a spiral there.
That is also why the damage looks nothing like the earlier figures. Instead of axis-aligned boxes, you can see jagged tears coming from a rounding error amplified by the iteration dynamics of the system.
Around the origin of the Julia set, the usable zoom is the square root of what the same type manages on a Mandelbrot, half the digits. Move the same view away from the origin and the glitch vanishes.
8. What 31 digits buys, in real-life units#
Double-double renders at zooms in the range of 1029.
These numbers are quite hard to get a feel for, so here are some fun facts.
Put the whole Mandelbrot set on your phone screen, about 7 cm wide, some 1200 pixels across. Now pinch to zoom and keep going until neighboring pixels start rounding to the same coordinate, the way they do in Figure 5. Then ask how big the whole set has grown by the time that happens.
Say every pinch doubles the zoom and takes one second.
In double precision you are done in 44 seconds, with the screen down to 1.3e-13 of a unit.
The entire set drawn at that scale would be 1,200,000,000 km across, 8x the distance from the Earth to the Sun, reaching from here to past Jupiter.
In DoubleDouble you keep pinching for another 53 seconds, down to 1.5e-29 across.
The original set is now 12 observable universes, side by side!
And it goes far, far deeper still! Figure 8 shows the same pinching kept up for 22 minutes, to a view 2.5e-398 across, which needs the center represented to 400+ decimal digits. No fixed-width type reaches that, and no amount of gluing doubles together ever will. That one takes perturbation: one reference orbit at full precision, and a per-pixel delta at lower precision (down there, the delta is around 1e-401, far below the smallest number even a double can represent).
Figure 8: A view 2.5e-398 across, 1,322 doublings below the whole set, computed by perturbation.[5]
Placing its center takes about 400 decimal digits.[6] A double has 15 and a double-double has 31.9. Where the double-double trick stops working#
These are double-double's edges:
- Range is unchanged, and the usable range is smaller. Overflow is still at 1e308, but full precision only lasts down to about 2e-292. Below that the low component goes subnormal, the exact product stops being exact, and the pair degrades toward a plain double.
- Infinities do not survive. The error terms compute
inf - inf, so the first overflow or division by zero hands back NaN where a double would hand back an infinity. The Mandelbrot kernel above is safe only because it stops iteration long before anything can grow that large. - NaN sorts instead of poisoning. Comparison is lexicographic on the two halves, so a NaN sorts below
every real value rather than making every comparison false the way a double's does, and
EqualsandCompareTodisagree about whether two NaNs are the same. - Not correctly rounded. The guarantees are relative-error bounds (about
2-104for the product), not "the nearest representable 106-bit value". Good enough for everything in the target band, but do not claim more. - Non-associativity is worse. Without correct rounding, reordering an expression alters the result more than it would with a plain
double. This can complicate testing. - The compiler must not "help". Every error-free transformation depends on the operations happening
exactly as written, at exactly double precision. C and C++ with
-ffast-mathwill reassociatesum - aaway, and FMA contraction, which is on by default in GCC and Clang without any fast-math flag, rewrites Dekker's splitting steps into something that is no longer a split. .NET is a comfortable place to write this: there is no fast-math switch, RyuJIT does not reassociate floating-point expressions, and it never contracts a multiply and an add into an FMA on its own. - ~31 digits is the ceiling. If the requirement is "arbitrary" or "whatever the input needs", this is the wrong tool. In my own fractal renderer double-double is one tier of many, and it hands off at 1e-28.
10. Takeaway#
Here are all the formats discussed in this article:
| Format | Decimal digits | perf vs. double | Pros | Cons |
|---|---|---|---|---|
float | ~7 | ~1x | half the memory, twice the SIMD lanes | low precision |
double | ~15 | 1x | fast, precise | precision may run out, which is why you are here |
long double (x87) | ~19 | ~1.2x | extra precision for free | x86 and GCC/Clang only, no SIMD, only 4 more digits |
DoubleDouble | ~31 | ~9x | no heap, no dependencies, inlines, vectorizes | infinities break, not correctly rounded, ~200 lines you own |
__float128 | ~34 | ~40x | high precision and range up to 1e4932 | emulated in software on x86, limited support |
| MPFR | as many as you ask | ~80x and up | unlimited precision | heap per value, P/Invoke, allocation discipline needed |
I knew that double-double is faster than MPFR, but it surprised me that it happens to be right in the middle. That 9x speedup over MPFR is why I wrote the type in the first place. It kept my own explorer interactive at extended zoom levels and made high-resolution renders practical.
Now, let's be honest here. Double-double is a niche type.
Most code never needs more than a double.
It only helps when you need more than 15 digits but fewer than 31.
However, if you need your code to run on, say, a GPU or an embedded system, a double-double can suddenly be a great option. It is just arithmetic on two doubles, while MPFR needs heap allocations and a native library. In those places it is not competing with MPFR, because MPFR cannot run there at all.
You should consider two things before using it. First, having 31 digits is not the same as using 31 digits. A wider type hands you more digits to start with but does nothing about how fast you lose them. And second, there are edge cases. Infinities turning into NaNs, nothing being correctly rounded, and full precision stopping at about 2e-292.
I sat down to write a short note about gluing two doubles together, and it kept growing. Partly because I decided to measure MPFR and binary128 instead of guessing, and partly because the fractals I was using as a demo kept showing me things I had not planned to write about.
And finally, for deep Mandelbrot zooms, brute-force precision is not the state of the art. Perturbation theory enables very deep zooms without the need for using arbitrary precision for all computations, but more on that some other time. The fractal was the motivation here, not the point. The gap between doubles and arbitrary-precision libraries shows up in more places than just fractal renderers, and if you end up using the trick, let me know!
Figure 9: And a little something extra at the end.
A view 3.8e-810 across, 2,693 doublings below the whole set, computed by perturbation.
And we are still not at the bottom!11. Further reading#
- T. J. Dekker, "A Floating-Point Technique for Extending the Available Precision", Numerische Mathematik, 1971.
The original paper: the splitting trick, and the fast two-sum this article calls
quickTwoSum. - D. E. Knuth, The Art of Computer Programming, vol. 2: Seminumerical Algorithms, section 4.2.2.
The six-operation
twoSumand the proof that it is exact. - Y. Hida, X. S. Li and D. H. Bailey, "Algorithms for Quad-Double Precision Floating Point Arithmetic", 2001. The QD library paper - the reference double-double and quad-double implementation.
- M. Joldes, J.-M. Muller and V. Popescu, "Tight and Rigorous Error Bounds for Basic Building Blocks of Double-Word Arithmetic", ACM TOMS, 2017. The modern error analysis: tight, proven bounds for the addition, multiplication and division used here.
- J. R. Shewchuk, "Adaptive Precision Floating-Point Arithmetic and Fast Robust Geometric Predicates", 1997. Where the technique goes next: expansions of n doubles whose precision adapts to the input.
Footnotes
The naming gets more confusing with C's
long double, whose meaning depends on the platform and the compiler. On IBM POWER it is actually double-double, exactly this technique (GCC calls it__ibm128, and reserves__float128for the IEEE type). On x86 it is yet another type: GCC and Clang give you the 80-bit x87 extended format with a 64-bit significand, while MSVC makes it just another name fordouble. ↩Printed values are rounded for readability.
Bis really the double nearest 0.555666777888, which drags a tail of extra digits from the decimal-to-binary conversion intoB,bLost,errorand the final sum, and the algorithm carries that tail exactly like the digits shown. ↩Medians of three runs on an idle machine; the spread between runs is a few percent either way. The iteration counts, on the other hand, are identical across every run and every type. ↩
gcc 14.2,
-O2, FMA contraction off, one P-core of the same i9-13900K; double-double is a line-for-line C port of Appendix A, and reproduces the .NET per-operation numbers within a few percent. Single-threaded, so compare the ratios, not the rows, with the all-core table earlier. ↩Nobody had ever looked at this spot before, and without the coordinate written down nobody could find it again, including me. It is kind of like searching for one particular grain of sand in the entire observable universe. Even if you managed to find the right planet, the search would still be hopeless. And like space, the set is mostly empty, with everything worth looking at on the boundary. ↩
Written down here so it is not lost, 406 digits each. Real part:
-1.2546501186909057505269466250921981644070768100251708177516070303734046805262190991942188750230768678605788326475221138767425434553956510745751843244663895032007841539877137705604003669583875944006181706618704710525415050302476211752688927355924699367031357139525690997877831869648195623522823691538908055742989329820011761064778120010166893341652782500699958380682058425509820660627053365326418384253599781. Imaginary part:0.3819035642800001819461581237720802452234485990147800334446378567314132368999769468860848762828896508632064176850184911749715609976648123007802564910929725517436736759042462717235720004499322510654559848620036722449528324306932522044048459018966997955517496317108858725278768456014781914592217389991723653848006887239802789886568674879888329442270299742538205772252794459812326617996073246041743419429748869. The viewport is 1.906e-398 high, on a 4:3 frame. ↩
Appendix A: the whole type#
Everything discussed above, in two files.
First the arithmetic half, depending on nothing but System - this is the part to copy when you only need the math:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
// Double-double arithmetic in ~200 lines. This file depends on nothing but System, on
// purpose: it is meant to be readable top to bottom and droppable into any project.
// Printing and parsing need a big-integer detour, so they live in DoubleDoubleFormat.cs.
namespace DoubleDoubleSample;
/// <summary>
/// A value stored as the unevaluated sum of two doubles: ~106 significand bits
/// (~31 decimal digits) at a few times the cost of a double, with no heap allocation.
///
/// Invariant: Hi == RoundToNearest(Hi + Lo), i.e. |Lo| is at most half an ulp of Hi.
/// Every factory and operation below maintains it; comparisons and ToDouble rely on it.
///
/// Built on the Dekker/Knuth error-free transformations at the bottom of this file;
/// products use the hardware fused multiply-add. The exponent range is a plain double's
/// (this buys precision, not range), and non-finite values are not special-cased: the
/// first operation on an infinity yields NaN (the error terms compute inf - inf), so
/// test for escape before a value can blow up.
///
/// By Marek Fiser, from marekfiser.com/blog/double-double-arithmetic/ (CC BY 4.0).
/// </summary>
public readonly partial struct DoubleDouble : IEquatable<DoubleDouble>, IComparable<DoubleDouble> {
private readonly double m_hi;
private readonly double m_lo;
/// <summary>The leading component; equals the value rounded to nearest double.</summary>
public double Hi => m_hi;
/// <summary>The trailing component (the rounding error of <see cref="Hi"/>).</summary>
public double Lo => m_lo;
public bool IsZero => m_hi == 0.0;
/// <summary>-1, 0 or +1 (0 also for NaN).</summary>
public int Sign => m_hi > 0.0 ? 1 : m_hi < 0.0 ? -1 : 0;
public static DoubleDouble Zero { get; } = new DoubleDouble(0.0, 0.0);
public static DoubleDouble One { get; } = new DoubleDouble(1.0, 0.0);
private DoubleDouble(double hi, double lo) {
m_hi = hi;
m_lo = lo;
}
/// <summary>Exact, so the trailing component is zero.</summary>
public static DoubleDouble FromDouble(double value) {
return new DoubleDouble(value, 0.0);
}
/// <summary>
/// Creates the exact sum <paramref name="hi"/> + <paramref name="lo"/> of two arbitrary
/// doubles, renormalizing so the invariant holds.
/// </summary>
public static DoubleDouble FromTwoDoubles(double hi, double lo) {
(double sum, double error) = twoSum(hi, lo);
return new DoubleDouble(sum, error);
}
/// <summary>Relative error stays below 3 * 2^-106 for all inputs.</summary>
public DoubleDouble Add(DoubleDouble right) {
// Add the two components pairwise, then incorporate lowSum and lowError one by one.
// Adding lowError can itself round, which is what the second quickTwoSum cleans up.
(double sum, double error) = twoSum(m_hi, right.m_hi);
(double lowSum, double lowError) = twoSum(m_lo, right.m_lo);
error += lowSum;
(sum, error) = quickTwoSum(sum, error);
error += lowError;
(sum, error) = quickTwoSum(sum, error);
return new DoubleDouble(sum, error);
}
public DoubleDouble Subtract(DoubleDouble right) {
return Add(right.Negate());
}
public DoubleDouble Multiply(DoubleDouble right) {
// hi*hi exactly, then the three cross terms - each already below the error term's
// magnitude, so ordinary double addition is accurate enough for them.
(double product, double error) = twoProduct(m_hi, right.m_hi);
error += m_hi * right.m_lo + m_lo * right.m_hi + m_lo * right.m_lo;
(product, error) = quickTwoSum(product, error);
return new DoubleDouble(product, error);
}
/// <summary>Cheaper than the full product, with two cross terms fewer.</summary>
public DoubleDouble Multiply(double right) {
(double product, double error) = twoProduct(m_hi, right);
error += m_lo * right;
(product, error) = quickTwoSum(product, error);
return new DoubleDouble(product, error);
}
/// <summary>
/// Divides this value by <paramref name="right"/>. There is no error-free
/// transformation for division, so this is long division: take a double-precision digit
/// of the quotient, subtract its contribution exactly, repeat.
/// </summary>
public DoubleDouble Divide(DoubleDouble right) {
double quotient1 = m_hi / right.m_hi;
DoubleDouble remainder = Subtract(right.Multiply(quotient1));
double quotient2 = remainder.m_hi / right.m_hi;
remainder = remainder.Subtract(right.Multiply(quotient2));
double quotient3 = remainder.m_hi / right.m_hi;
(double sum, double error) = quickTwoSum(quotient1, quotient2);
return new DoubleDouble(sum, error).Add(FromDouble(quotient3));
}
/// <summary>Cheaper than the general product, with two cross terms instead of three.</summary>
public DoubleDouble Square() {
(double product, double error) = twoProduct(m_hi, m_hi);
error += 2.0 * m_hi * m_lo + m_lo * m_lo;
(product, error) = quickTwoSum(product, error);
return new DoubleDouble(product, error);
}
/// <summary>Exact, since negation never rounds.</summary>
public DoubleDouble Negate() {
return new DoubleDouble(-m_hi, -m_lo);
}
public DoubleDouble Abs() {
return Sign < 0 ? Negate() : this;
}
/// <summary>
/// The square root, by one Newton refinement of the double root: faithful to about an
/// ulp of the double-double, not correctly rounded. Zero stays zero; negative and NaN
/// inputs resolve through the double root, so non-finite values propagate as everywhere
/// else in this arithmetic.
/// </summary>
public DoubleDouble Sqrt() {
if (!(m_hi > 0.0)) { // deliberately not "m_hi <= 0.0": this form is also true for NaN
return FromDouble(Math.Sqrt(m_hi));
}
double approx = Math.Sqrt(m_hi);
DoubleDouble residual = Subtract(FromDouble(approx).Square());
double correction = residual.m_hi / (approx + approx);
(double hi, double lo) = quickTwoSum(approx, correction);
return new DoubleDouble(hi, lo);
}
/// <summary>The nearest double, which by the invariant is the high component itself.</summary>
public double ToDouble() {
return m_hi;
}
/// <summary>Compares values; thanks to the normalization invariant this is simply
/// lexicographic on (Hi, Lo).</summary>
public int CompareTo(DoubleDouble other) {
int hiComparison = m_hi.CompareTo(other.m_hi);
return hiComparison != 0 ? hiComparison : m_lo.CompareTo(other.m_lo);
}
public bool Equals(DoubleDouble other) {
return m_hi == other.m_hi && m_lo == other.m_lo;
}
public override bool Equals(object? obj) {
return obj is DoubleDouble other && Equals(other);
}
public override int GetHashCode() {
return HashCode.Combine(m_hi, m_lo);
}
public static DoubleDouble operator +(DoubleDouble left, DoubleDouble right) => left.Add(right);
public static DoubleDouble operator -(DoubleDouble left, DoubleDouble right) => left.Subtract(right);
public static DoubleDouble operator *(DoubleDouble left, DoubleDouble right) => left.Multiply(right);
public static DoubleDouble operator *(DoubleDouble left, double right) => left.Multiply(right);
public static DoubleDouble operator /(DoubleDouble left, DoubleDouble right) => left.Divide(right);
public static DoubleDouble operator -(DoubleDouble value) => value.Negate();
public static bool operator ==(DoubleDouble left, DoubleDouble right) => left.Equals(right);
public static bool operator !=(DoubleDouble left, DoubleDouble right) => !left.Equals(right);
public static bool operator <(DoubleDouble left, DoubleDouble right) => left.CompareTo(right) < 0;
public static bool operator >(DoubleDouble left, DoubleDouble right) => left.CompareTo(right) > 0;
public static bool operator <=(DoubleDouble left, DoubleDouble right) => left.CompareTo(right) <= 0;
public static bool operator >=(DoubleDouble left, DoubleDouble right) => left.CompareTo(right) >= 0;
/// <summary>Widening a double is exact, so it may happen implicitly.</summary>
public static implicit operator DoubleDouble(double value) => FromDouble(value);
/// <summary>Narrowing back to a double loses half the digits, so it must be asked for.</summary>
public static explicit operator double(DoubleDouble value) => value.ToDouble();
/// <summary>
/// Knuth's two-sum: sum + error == a + b <em>exactly</em>, for any two doubles. The
/// rounding error of a floating-point addition is itself a representable double, and
/// these six operations recover it. Branch-free; the intermediates may round, but the
/// steps are arranged so that whatever one loses another accounts for.
/// </summary>
private static (double Sum, double Error) twoSum(double a, double b) {
double sum = a + b;
double bKept = sum - a; // the part of b that made it into the sum
double aKept = sum - bKept; // and the part of a
double aLost = a - aKept;
double bLost = b - bKept;
return (sum, aLost + bLost);
}
/// <summary>Dekker's fast two-sum: the same guarantee in three operations, but only
/// valid when |a| >= |b| (or a == 0). Used where the caller knows the order.</summary>
private static (double Sum, double Error) quickTwoSum(double a, double b) {
double sum = a + b;
double bKept = sum - a;
return (sum, b - bKept); // a kept everything, so only b's loss is left
}
/// <summary>
/// Exact product: product + error == a * b exactly. The full product needs 106 bits;
/// the fused multiply-add computes a*b - product with a single rounding, which is
/// precisely the missing low half. This is what makes modern double-double
/// multiplication cheap (Dekker's original splitting trick needed 17 operations
/// instead of 2).
///
/// <para>Practically every mainstream CPU of the last decade has an FMA instruction;
/// on one that does not, Math.FusedMultiplyAdd computes the same answer in software,
/// so the type stays correct and merely stops being cheap.</para>
///
/// <para>The error is representable only while the product itself stays normal. Below
/// about 2e-292 the exact tail no longer fits a double and this stops being exact.</para>
/// </summary>
private static (double Product, double Error) twoProduct(double a, double b) {
double product = a * b;
double error = Math.FusedMultiplyAdd(a, b, -product);
return (product, error);
}
}
And the human half: exact printing and parsing, which cannot stay inside doubles and takes the BigInteger detour:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
// The half of DoubleDouble that has to leave the world of doubles. Arithmetic never needs
// this: only humans do. Printing 31 correct decimal digits means computing them exactly,
// which is a job for BigInteger (in the BCL, so still zero NuGet dependencies).
using System.Globalization;
using System.Numerics;
using System.Text;
namespace DoubleDoubleSample;
public readonly partial struct DoubleDouble {
/// <summary>Significant digits <see cref="ToString()"/> prints: enough to round-trip any
/// value whose two components are contiguous (see <see cref="Parse"/>).</summary>
private const int DEFAULT_DIGITS = 36;
/// <summary>Formats with <see cref="DEFAULT_DIGITS"/> significant digits.</summary>
public override string ToString() {
return ToString(DEFAULT_DIGITS);
}
/// <summary>Formats with the given number of significant decimal digits, correctly rounded.</summary>
public string ToString(int significantDigits) {
if (significantDigits < 1) {
throw new ArgumentOutOfRangeException(nameof(significantDigits));
}
if (!double.IsFinite(m_hi)) {
return m_hi.ToString(CultureInfo.InvariantCulture);
}
(BigInteger numerator, BigInteger denominator) = toRational(this);
return formatRational(numerator, denominator, significantDigits);
}
/// <summary>
/// Formats the exact value. A double-double is a sum of two binary fractions, so it
/// always has a finite decimal expansion - just a long one (hundreds of digits when the
/// components are far apart).
/// </summary>
public string ToStringExact() {
if (!double.IsFinite(m_hi)) {
return m_hi.ToString(CultureInfo.InvariantCulture);
}
(BigInteger numerator, BigInteger denominator) = toRational(this);
if (numerator.IsZero) {
return "0";
}
// denominator is a power of two, so multiplying by the matching power of five turns
// it into a power of ten and the digits fall out directly.
int twos = (int)(denominator.GetBitLength() - 1);
BigInteger scaled = BigInteger.Abs(numerator) * BigInteger.Pow(5, twos);
string digits = scaled.ToString(CultureInfo.InvariantCulture).PadLeft(twos + 1, '0');
string sign = numerator.Sign < 0 ? "-" : "";
if (twos == 0) {
return sign + digits;
}
string result = sign + digits[..^twos] + "." + digits[^twos..];
return result.TrimEnd('0').TrimEnd('.');
}
/// <summary>
/// Parses a decimal string to the nearest double-double (the value is built exactly as a
/// rational, then rounded once into each component). Round-trips with
/// <see cref="ToString()"/> for every value whose components are <em>contiguous</em> -
/// a hand-built pair such as (1.0, 1e-300) carries more information than 36 digits can.
/// </summary>
public static DoubleDouble Parse(string text) {
(BigInteger numerator, BigInteger denominator) = parseRational(text);
double hi = toNearestDouble(numerator, denominator);
if (!double.IsFinite(hi)) {
return FromDouble(hi);
}
// Subtract the leading component exactly, and round what is left into the trailing one.
(BigInteger hiNumerator, BigInteger hiDenominator) = toRational(FromDouble(hi));
BigInteger restNumerator = numerator * hiDenominator - hiNumerator * denominator;
BigInteger restDenominator = denominator * hiDenominator;
double lo = toNearestDouble(restNumerator, restDenominator);
return FromTwoDoubles(hi, lo);
}
/// <summary>The exact value of a double-double as a rational (the denominator is a power of two).</summary>
private static (BigInteger Numerator, BigInteger Denominator) toRational(DoubleDouble value) {
(BigInteger hiNumerator, BigInteger hiDenominator) = toRational(value.m_hi);
(BigInteger loNumerator, BigInteger loDenominator) = toRational(value.m_lo);
return (hiNumerator * loDenominator + loNumerator * hiDenominator, hiDenominator * loDenominator);
}
/// <summary>The exact value of a finite double as a rational, straight from its bits.</summary>
private static (BigInteger Numerator, BigInteger Denominator) toRational(double value) {
long bits = BitConverter.DoubleToInt64Bits(value);
int exponent = (int)((bits >> 52) & 0x7FF);
long mantissa = bits & 0xF_FFFF_FFFF_FFFFL;
if (exponent == 0) {
exponent = -1074; // subnormal: no implicit leading bit
} else {
mantissa |= 1L << 52;
exponent -= 1075;
}
BigInteger numerator = bits < 0 ? -mantissa : mantissa;
return exponent >= 0
? (numerator << exponent, BigInteger.One)
: (numerator, BigInteger.One << -exponent);
}
/// <summary>Rounds an exact rational to the nearest double, ties to even.</summary>
private static double toNearestDouble(BigInteger numerator, BigInteger denominator) {
if (numerator.IsZero) {
return 0.0;
}
int sign = numerator.Sign;
numerator = BigInteger.Abs(numerator);
// Scale so the quotient lands in [2^52, 2^54): one integer division then gives every
// bit of the significand plus the remainder needed to decide the rounding.
long exponent = numerator.GetBitLength() - denominator.GetBitLength() - 53;
if (exponent > 0) {
denominator <<= (int)exponent;
} else {
numerator <<= (int)-exponent;
}
BigInteger quotient = BigInteger.DivRem(numerator, denominator, out BigInteger remainder);
if (quotient.GetBitLength() > 53) {
// One bit too many: give it back to the exponent, adding the dropped bit into
// the remainder so the rounding decision below stays exact.
if (!quotient.IsEven) {
remainder += denominator;
}
quotient >>= 1;
denominator <<= 1;
exponent++;
}
int comparison = (remainder << 1).CompareTo(denominator);
if (comparison > 0 || (comparison == 0 && !quotient.IsEven)) {
quotient++;
if (quotient.GetBitLength() > 53) {
quotient >>= 1;
exponent++;
}
}
return sign * Math.ScaleB((double)quotient, (int)exponent);
}
/// <summary>Turns "-1.25e-7" and friends into an exact rational.</summary>
private static (BigInteger Numerator, BigInteger Denominator) parseRational(string text) {
text = text.Trim();
if (text.Length == 0) {
throw new FormatException("Empty input.");
}
int index = 0;
bool negative = text[index] == '-';
if (negative || text[index] == '+') {
index++;
}
BigInteger mantissa = BigInteger.Zero;
int fractionDigits = 0;
bool seenDot = false;
bool seenDigit = false;
for (; index < text.Length; index++) {
char c = text[index];
if (c == '.') {
if (seenDot) {
throw new FormatException($"Two decimal points in '{text}'.");
}
seenDot = true;
} else if (c is >= '0' and <= '9') {
mantissa = mantissa * 10 + (c - '0');
seenDigit = true;
if (seenDot) {
fractionDigits++;
}
} else if (c is 'e' or 'E') {
break;
} else {
throw new FormatException($"Unexpected character '{c}' in '{text}'.");
}
}
if (!seenDigit) {
throw new FormatException($"No digits in '{text}'.");
}
int exponent = -fractionDigits;
if (index < text.Length) {
exponent += int.Parse(text[(index + 1)..], CultureInfo.InvariantCulture);
}
if (negative) {
mantissa = -mantissa;
}
return exponent >= 0
? (mantissa * BigInteger.Pow(10, exponent), BigInteger.One)
: (mantissa, BigInteger.Pow(10, -exponent));
}
/// <summary>Rounds an exact rational to N significant digits and lays them out.</summary>
private static string formatRational(BigInteger numerator, BigInteger denominator, int digits) {
if (numerator.IsZero) {
return "0";
}
bool negative = numerator.Sign < 0;
numerator = BigInteger.Abs(numerator);
// Decimal exponent of the leading digit, from the bit lengths (log10(2) ~ 0.30103),
// then corrected by the one-step loop below when the estimate is off by one.
int exponent10 = (int)Math.Floor((numerator.GetBitLength() - denominator.GetBitLength()) * 0.30103);
string text;
while (true) {
int scale = digits - 1 - exponent10;
BigInteger scaledNumerator = numerator;
BigInteger scaledDenominator = denominator;
if (scale >= 0) {
scaledNumerator *= BigInteger.Pow(10, scale);
} else {
scaledDenominator *= BigInteger.Pow(10, -scale);
}
BigInteger rounded = roundToNearest(scaledNumerator, scaledDenominator);
text = rounded.ToString(CultureInfo.InvariantCulture);
if (text.Length == digits) {
break;
}
// Off by one either way (estimate too low, or rounding carried into a new digit).
exponent10 += text.Length - digits;
}
text = text.TrimEnd('0');
if (text.Length == 0) {
text = "0";
}
StringBuilder result = new(negative ? "-" : "");
if (exponent10 >= -5 && exponent10 < digits) {
// Positional, the range where it stays readable.
if (exponent10 >= 0) {
text = text.PadRight(exponent10 + 1, '0');
result.Append(text[..(exponent10 + 1)]);
if (text.Length > exponent10 + 1) {
result.Append('.').Append(text[(exponent10 + 1)..]);
}
} else {
result.Append("0.").Append('0', -exponent10 - 1).Append(text);
}
} else {
result.Append(text[0]);
if (text.Length > 1) {
result.Append('.').Append(text[1..]);
}
result.Append('e').Append(exponent10.ToString(CultureInfo.InvariantCulture));
}
return result.ToString();
}
/// <summary>Integer nearest to numerator/denominator, ties away from zero (both positive).</summary>
private static BigInteger roundToNearest(BigInteger numerator, BigInteger denominator) {
BigInteger quotient = BigInteger.DivRem(numerator, denominator, out BigInteger remainder);
return (remainder << 1) >= denominator ? quotient + 1 : quotient;
}
}
Appendix B: the benchmarked kernels#
The numbers in section 5: What it actually costs come from these benchmarks. The pixel loop and the threading are the same for all three and are left out.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
namespace DoubleDoubleSample.Mandelbrot;
/// <summary>
/// One escape-time kernel, three ways, plus the scratch the MPFR one needs. Everything else -
/// parsing the center, dividing up rows, timing - is in MandelbrotRenderer.cs, so the three can be
/// read against each other here without the driver in the way.
///
/// <para>All three run the same seven operations per iteration: two squares, one multiply, a
/// doubling, two adds and a subtract. The first two differ only in the declared type of the
/// coordinates, which is the point of the exercise. The third is the same arithmetic again as calls
/// into a library, and it is the shape of that third one, rather than the arithmetic in it, that
/// costs the order of magnitude.</para>
///
/// <para>Each does its escape test in whatever way is cheapest for the type it is written in. The
/// test's precision does not matter: past |z| > 2 the orbit diverges, and by the time it crosses a
/// bailout of 65536 it is squaring itself each iteration, so the last bits could only move the
/// escape by an iteration - measurably, they do not move it at all, and both MPFR variants below
/// count the same iterations to the digit as the double-double one. For double-double the test is
/// two field reads. For MPFR it is a native compare against a preallocated constant, because
/// getting a double out of an MPFR value costs a call, and that call would be per iteration.</para>
///
/// <para>A fourth variant, which allocates a fresh MPFR value per operation the way a normal
/// wrapper does, lives beside the driver in MandelbrotRenderer.cs. It is these same calls with the
/// lifetime bookkeeping spelled out, and it is four times slower for it.</para>
/// </summary>
public static unsafe partial class MandelbrotRenderer {
/// <summary>Bailout radius squared. Generous, so the smooth-iteration formula is accurate.</summary>
private const double ESCAPE_SQUARED = 65536.0 * 65536.0;
/// <summary>Plain double, for reference and for the pictures that fall apart.</summary>
private static double iterateDouble(
double cRe,
double cIm,
int maxIterations,
out int iterations
) {
double zRe = 0.0;
double zIm = 0.0;
double magnitudeSquared = 0.0;
int n = 0;
for (; n < maxIterations; n++) {
double zReSquared = zRe * zRe;
double zImSquared = zIm * zIm;
magnitudeSquared = zReSquared + zImSquared;
if (magnitudeSquared > ESCAPE_SQUARED) {
break;
}
zIm = 2.0 * zRe * zIm + cIm;
zRe = zReSquared - zImSquared + cRe;
}
iterations = n;
return magnitudeSquared;
}
/// <summary>
/// The same source, with <c>DoubleDouble</c> substituted for <c>double</c>. That is the whole
/// diff: operators, no allocation, no scratch, nothing to release. The escape test reads the
/// leading halves directly, which costs nothing at all.
/// </summary>
private static double iterateDoubleDouble(
DoubleDouble cRe,
DoubleDouble cIm,
int maxIterations,
out int iterations
) {
DoubleDouble zRe = DoubleDouble.Zero;
DoubleDouble zIm = DoubleDouble.Zero;
double magnitudeSquared = 0.0;
int n = 0;
for (; n < maxIterations; n++) {
DoubleDouble zReSquared = zRe.Square();
DoubleDouble zImSquared = zIm.Square();
magnitudeSquared = zReSquared.Hi + zImSquared.Hi;
if (magnitudeSquared > ESCAPE_SQUARED) {
break;
}
zIm = zRe * zIm * 2.0 + cIm;
zRe = zReSquared - zImSquared + cRe;
}
iterations = n;
return magnitudeSquared;
}
/// <summary>
/// The same recurrence with every double-double operation inlined, which lets the kernel do
/// three things a caller of the type cannot. It is the same arithmetic built from the same
/// error-free transformations, and it renders the same picture; what it drops is work that only
/// exists to hand a well-formed value back to a caller who, here, is the next line of this loop.
///
/// <list type="number">
/// <item><b>Two additions instead of three.</b> <c>zRe^2 - zIm^2 + cRe</c> is a three-term sum
/// and needs one renormalization, not the two it gets when spelled as two <c>Add</c> calls.
/// Addition is the expensive operation in this type, so this is most of the saving.</item>
/// <item><b>The squares are left unnormalized.</b> The closing <c>quickTwoSum</c> restores an
/// invariant for a recipient; here the recipients are an escape test that reads the leading
/// component and a subtraction that starts with an exact <c>twoSum</c> of leading components
/// anyway.</item>
/// <item><b>The doubling is componentwise.</b> Exact for a power of two, and it reuses the
/// <c>zRe + zRe</c> the square already needed.</item>
/// </list>
///
/// <para>The escape test is negated so a non-finite magnitude reads as escaped. This kernel
/// cannot produce one, but the cost of the safe spelling is zero and the failure mode it avoids
/// is silent: NaN fails every comparison, so a pixel would never escape and would paint as
/// interior.</para>
/// </summary>
private static double iterateDoubleDoubleFused(
DoubleDouble cRe,
DoubleDouble cIm,
int maxIterations,
out int iterations
) {
double cReHi = cRe.Hi, cReLo = cRe.Lo;
double cImHi = cIm.Hi, cImLo = cIm.Lo;
double zReHi = 0.0, zReLo = 0.0;
double zImHi = 0.0, zImLo = 0.0;
double magnitudeSquared = 0.0;
int n = 0;
for (; n < maxIterations; n++) {
// Both corrections of each square ride one FMA that never touches the leading product's
// dependency chain, so they compute alongside it rather than after it.
double doubleReLo = zReLo + zReLo;
double reSquare = zReHi * zReHi;
double reSquareError = Math.FusedMultiplyAdd(zReHi, zReHi, -reSquare)
+ Math.FusedMultiplyAdd(zReHi, doubleReLo, zReLo * zReLo);
double imSquare = zImHi * zImHi;
double imSquareError = Math.FusedMultiplyAdd(zImHi, zImHi, -imSquare)
+ Math.FusedMultiplyAdd(zImHi, zImLo + zImLo, zImLo * zImLo);
magnitudeSquared = reSquare + imSquare;
if (!(magnitudeSquared <= ESCAPE_SQUARED)) {
break;
}
// zIm = 2*zRe*zIm + cIm. The doubling is folded into the left operand, and the three
// cross terms are summed as a tree rather than chained through one accumulator.
double doubleReHi = zReHi + zReHi;
double cross = doubleReHi * zImHi;
double crossError = Math.FusedMultiplyAdd(doubleReHi, zImHi, -cross);
crossError = Math.FusedMultiplyAdd(doubleReLo, zImLo, crossError)
+ Math.FusedMultiplyAdd(doubleReHi, zImLo, doubleReLo * zImHi);
(double imSum, double imResidual) = twoSum(cross, cImHi);
(zImHi, zImLo) = quickTwoSum(imSum, imResidual + crossError + cImLo);
// zRe = zRe^2 - zIm^2 + cRe. The leading terms go through twoSum, so the cancellation
// this line is famous for is captured exactly; only the corrections are approximate.
(double difference, double differenceResidual) = twoSum(reSquare, -imSquare);
(double reSum, double reResidual) = twoSum(difference, cReHi);
double reTail = differenceResidual + reResidual + (reSquareError - imSquareError) + cReLo;
(zReHi, zReLo) = quickTwoSum(reSum, reTail);
}
iterations = n;
return magnitudeSquared;
}
/// <summary>Knuth's two-sum, as in the type. Private to the kernel because the type keeps its
/// own copy private; this file is the one place that needs it without a DoubleDouble around it.</summary>
private static (double Sum, double Error) twoSum(double a, double b) {
double sum = a + b;
double bKept = sum - a;
double aKept = sum - bKept;
return (sum, (a - aKept) + (b - bKept));
}
/// <summary>Dekker's fast two-sum, valid where the caller knows the first argument dominates.</summary>
private static (double Sum, double Error) quickTwoSum(double a, double b) {
double sum = a + b;
double bKept = sum - a;
return (sum, b - bKept);
}
/// <summary>
/// MPFR with destinations allocated once and written in place, which is how you would write this
/// in C and is MPFR at its best. Nothing here reaches the heap: at these precisions MPFR takes
/// its scratch space off the stack.
/// </summary>
private static double iterateInPlace(
MpfrScratch s,
int maxIterations,
out int iterations
) {
MpfrNative.SetDouble(s.ZRe, 0.0);
MpfrNative.SetDouble(s.ZIm, 0.0);
int n = 0;
for (; n < maxIterations; n++) {
MpfrNative.mpfr_sqr(s.ZReSquared, s.ZRe, MpfrNative.ROUND_NEAREST);
MpfrNative.mpfr_sqr(s.ZImSquared, s.ZIm, MpfrNative.ROUND_NEAREST);
MpfrNative.mpfr_add(s.Temp, s.ZReSquared, s.ZImSquared, MpfrNative.ROUND_NEAREST);
if (MpfrNative.mpfr_cmp(s.Temp, s.Escape) > 0) {
break;
}
MpfrNative.mpfr_mul(s.Temp, s.ZRe, s.ZIm, MpfrNative.ROUND_NEAREST);
MpfrNative.mpfr_add(s.Temp, s.Temp, s.Temp, MpfrNative.ROUND_NEAREST); // exact doubling
MpfrNative.mpfr_add(s.ZIm, s.Temp, s.CIm, MpfrNative.ROUND_NEAREST);
MpfrNative.mpfr_sub(s.Temp, s.ZReSquared, s.ZImSquared, MpfrNative.ROUND_NEAREST);
MpfrNative.mpfr_add(s.ZRe, s.Temp, s.CRe, MpfrNative.ROUND_NEAREST);
}
iterations = n;
// Converted once per pixel rather than once per iteration, and only the escape path uses it:
// a pixel that ran out of iterations is colored by its count alone. Temp still holds the
// magnitude there, because the loop broke before reaching the line that reuses it.
return MpfrNative.GetDouble(s.Temp);
}
/// <summary>
/// Per-worker MPFR values. MPFR values cannot be shared across threads for writing, so
/// <c>Parallel.For</c>'s localInit hands each worker its own set and localFinally releases them.
/// </summary>
private sealed class MpfrScratch : IDisposable {
public readonly void* CRe;
public readonly void* CIm;
public readonly void* ZRe;
public readonly void* ZIm;
public readonly void* ZReSquared;
public readonly void* ZImSquared;
public readonly void* Temp;
public readonly void* Temp2;
/// <summary>The bailout radius as an MPFR value, so the escape test never has to convert.</summary>
public readonly void* Escape;
public long Iterations;
public MpfrScratch(int precision) {
CRe = MpfrNative.Allocate(precision);
CIm = MpfrNative.Allocate(precision);
ZRe = MpfrNative.Allocate(precision);
ZIm = MpfrNative.Allocate(precision);
ZReSquared = MpfrNative.Allocate(precision);
ZImSquared = MpfrNative.Allocate(precision);
Temp = MpfrNative.Allocate(precision);
Temp2 = MpfrNative.Allocate(precision);
Escape = MpfrNative.Allocate(precision);
MpfrNative.SetDouble(Escape, ESCAPE_SQUARED);
}
/// <summary>Every value is live here: both kernels restore that before returning.</summary>
public void Dispose() {
MpfrNative.Free(CRe);
MpfrNative.Free(CIm);
MpfrNative.Free(ZRe);
MpfrNative.Free(ZIm);
MpfrNative.Free(ZReSquared);
MpfrNative.Free(ZImSquared);
MpfrNative.Free(Temp);
MpfrNative.Free(Temp2);
MpfrNative.Free(Escape);
}
}
}
The escape test of the double-double kernel differs slightly, but it's just a small optimization with no precision loss, because the orbit past |z| > 2 is diverging exponentially and the low part has no real effect on the number of iterations.
The fourth kernel (the last row of the kernel table) is the third one with each operation producing a new value instead of writing into the scratch.