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 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.

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.
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.

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.

I benchmarked all of that, and the gap turned out to be surprisingly well structured. 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. Let MPFR allocate a result per operation and the double-double-to-MPFR gap is closer to 40x instead of 9x.

The whole trick is:

Keep two doubles separate, but treat them as one number.

About 31 decimal digits, no heap, no dependencies, and no native binary to ship.

The complete type and the benchmarked kernels are in the two appendices.

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 only determines where the decimal point will be. Take these two numbers:

A = 111222333444``B = 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 about 15.
Evaluate A + B and the digits that do not fit are rounded away, leaving 111222333444.55566.

Now the key idea: 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 theunevaluated sum of two doubles.

x = xhi + xlowith the invariant that x is exactly what hix rounds to as a double.
The high part carries the value, the low part carries the rest (the "error" of the high part after addition).
Every operation in the next section is just keeping this invariant true.hi + xlo

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

public readonly partial struct DoubleDouble : IEquatable, IComparable {
private readonly double m_hi;
private readonly double m_lo;
///

/// The leading component; equals the value rounded to nearest double.
///
public double Hi => m_hi;
///
/// The trailing component (the rounding error of ).
///
public double Lo => m_lo;
}

One naming note before the math. The 128-bit type in IEEE 754 (the standard that also defines the double) is called binary128, better known as quadruple-precision.
It holds one contiguous 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 combined ~106 significand bits and unchanged exponent range (still overflowing at ~1e308).[1]

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.
Call the rounded result s = round(a + b).
The part that got rounded away is exactly (a + b) - s, a formula that looks self-defeating: it needs the exact sum, and the exact sum is the one thing we do not have.
And yet a handful of ordinary double operations recover the leftover exactly.
The algorithm is Knuth's two-sum, and Code listing 2 is the whole of it, six additions
and subtractions, with no branches, no bit twiddling and no wider type.

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 every line of that 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 hardware 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, and it always fits in one double.
It has to, because everything below the cut came from the smaller operand, and the smaller operand never had more than 53 bits.
In the extreme, when the operands are too far apart to overlap at all, the sum rounds to the larger one and the leftover is the smaller one, whole.

Which makes the last line the property the whole type is built on: (sum, error) is exactly A + B, to the last digit.
The exactness comes from the arrangement rather than from rounding errors cancelling out; only the first two lines round at all.
aKept and bKept split sum exactly between them, so the last four lines are exact and aLost + bLost is precisely what the first line threw away.
For the proofs, see Dekker 1971 and Knuth vol. 2 in Further reading.

Notice that A lost nothing in the example above, 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 double and a low double. Dekker's original method (Veltkamp splitting) took 17 operations. Once the hardware has an FMA, 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 a fused multiply-add 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 full. The subtraction cancels
everything product kept, leaving exactly the tail that a * b threw away, small enough to come back
without any rounding at all.

Here are all three side by side:

| Transformation | Operations | Precondition | Gives you |
|---|---|---|---|
| twoSum | 6 | none | exact sum of two doubles |
| quickTwoSum | 3 | larger operand first | exact sum of two doubles |
| twoProduct | 2 (one FMA) | product must not over- or underflow | exact product of two doubles |

These three are called error-free transformations, and they are the entire foundation. They all
assume round-to-nearest arithmetic and no overflow, which .NET gives you and offers no switch to take
away. Everything below is bookkeeping on top of them.

Building arithmetic operations for double-double

With twoSum, quickTwoSum, and twoProduct functions in hand, writing arithmetic operations for DoubleDouble is quite straightforward,
and one pattern repeats through all of them: use an error-free transformation wherever a rounding would lose
digits we care about, use plain double arithmetic wherever the digits at stake are already below
2, and finish with a -106quickTwoSum so the result is a normalized pair again.

Addition

Add the two components pairwise, then fold the errors back in.
The fold itself can round, which is what the second quickTwoSum cleans up.

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, 11 operations instead of these 20, which folds both low parts in one step. It is fine while the operands share a sign, but under cancellation it has no relative error bound at all.

Multiplication

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,
so whatever plain double arithmetic rounds away from them is at -532 and below,
under the precision floor of the type itself.
Knowing what you are allowed to be sloppy about is most of the art here.-106

Squaring is the same with two cross terms instead of three,
and multiplying by a plain double needs only one - worth having as overloads, because kernels use them constantly (say, a pixel step times a loop index).

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. Long division, three digits deep.

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 is the same shape: take the double root, then one Newton refinement. Both come out accurate to a couple of units in the last place rather than correctly rounded - not the nearest representable pair, and not always within one of it. The error starts to matter only when the last bits have to be reproducible against a different implementation.

Comparison, for free

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);
}

This is plain lexicographic ordering; the invariant makes it correct.
Since the low half is never more than half the gap to the next double, it can never grow large enough to overturn a decision the high halves already made.
It is also one more reason every operation ends by renormalizing.
If unnormalized pairs were allowed, (1, 0.75) and (1.75, 0) would be the same number and compare as different ones.

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.

Single double:Any value one double holds exactly is the pair with- lo = 0; that is all- FromDouble(1.5)does.
Touching:- locontinues exactly where- hiends, giving 106 contiguous bits, the "31 digits" from the title. The smallest example is- 1 + 2, which becomes the pair- -53- (1, 2.- -53)
Gapped:This is the interesting case, where- lostarts further down, and every bit position in between is a 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:- FromTwoDoubles(1.0, 0.75)renormalizes on entry and returns- (1.75, 0).

The example's gap of 2 can be seen if we write the digits out in binary.
A owns all the integer bits, so the fractional bits of the exact sum are exactly B's own digits,
and the pair splits that digit stream in two. sum's 53-bit window (which started way up at 2)
runs out 16 places after the binary point, and 36error carries the rest:

1
2
3
4
5

B = .1000111001000000001011011000... sum = .1000111001000000 window ends at 2^-16 error = 001011011000... first set bit at 2^-19 ^^ the gap: two zero bits stored by neither half

And in our case, error's share of the stream happens to begin with two zeros, and that is all a gap is.
Those zeros are real digits of the value, but a float never spends storage on leading zeros.
Instead, its exponent simply points lower, the way 0.001011 and 1.011 â 2 are the same number.-3

The gapped shape has one catch, though.
The gap itself costs nothing and loses nothing; the pair holds exactly two "islands" of bits.
Suppose this absurd case: 10 is a perfectly legal pair with 1,940 zero bits between its islands, carried in 16 bytes.
However, evaluate 300 + 10-300(10 and the exact answer would need a third island in the middle, so the smallest one is dropped and the answer is simply 300 + 10-300) + 110.300 + 1

"But surely 10 written in binary does not have loads of trailing zerosâ½"
Right, it does not. Writing it exactly needs about 700 significant bits.
But the 3001e300 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 snap it back.
Printing runs into the same limit from the other side, which is why the type also carries a ToStringExact.

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

The comparison that matters is against MPFR, made in a real kernel rather than a loop that repeats one operation.

The kernel is the whole definition of the Mandelbrot set, which is a small thing to write down.
Take a point c in the complex plane, start at z, and iterate:0 = 0

zn+1 = zn2 + cIf the orbit stays bounded forever, c belongs to the set and the pixel is black.
If it escapes, and it provably has once |z| > 2, the iteration count on the way out is the color.
Every kernel number below is the cost of running those seven real operations, once per iteration, in one type or another.
Operands change every time, the escape test is included, and one frame is roughly 390 million iterations with every core busy.
MPFR goes through hand-written P/Invoke with no wrapper class in the way, at 106 bits, which is double-double's significand exactly.

| 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 |
| MPFR @212, in place | 57 M | ~150x | ~17x |
| MPFR @212, fresh result per operation | 22 M | ~400x | ~45x |

The two MPFR styles run the identical seven operations. The only difference is where results go: into destinations allocated once and reused, the way you would write it in C, or into a fresh value per operation, the way many wrappers present it.

A double-double costs roughly 9x a double, MPFR at the same 31 digits costs roughly 9x a double-double again, and those two steps multiply into the ~81x. There is no reason the option in the middle should land in the geometric middle of the gap it fills, but this one does. Take the API at face value and let it allocate, and your perf penalty shoots up to ~38x instead.

Read those as "about nine" and "about forty" rather than to the digit.[3]

Per-operation timings and whole-kernel throughput disagree, as they should. A loop repeating one operation on fixed operands is the friendliest case there is, not the one you ship. It is not a threading artifact either, running the same frame on one thread and on all of them speeds both types up by about the same factor.

The allocation is avoidable, and that is the catch

The 38x row is avoidable, but the discipline is manual and permanent: keep destinations alive across iterations, never write a * b + c as an expression because each operation would allocate a new value.
People do exactly that, and it works. It is also most of what you were hoping a number type would spare you.

Double-double 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.

What a fixed-width value type buys

The 9x is this kernel's number, not a constant. It comes from seven arithmetic operations per iteration on operands that are already in registers. Shift the mix toward division, where double-double is at its weakest, or toward memory, where neither type is doing the work, and it moves. What does not move is where the difference comes from:

There is no call.A double-double add or multiply is 7 to 20 instructions the JIT inlines into your loop (a divide is about eighty, and still inlined); an MPFR operation is an opaque function call (in .NET across the P/Invoke boundary, about 1.6 ns before any arithmetic starts).
No loops, no branches.A fixed instruction sequence the CPU pipelines, against a data-dependent walk over limbs, then a normalize, then a round.
It rides the FPU.Native- doubleadds and multiplies, retired one per cycle with several in flight, against integer paths plus software rounding.
Precision is decided at compile time.No precision field, no rounding policy, none of the bookkeeping a library spends much of its time on.
It can vectorize.Nothing branches, so a structure-of-arrays rewrite over- Vector256<double>does four values at a time and the same arithmetic ports to a GPU. That is a rewrite, not something the JIT will do to the type above, but arbitrary precision has no such rewrite available at all.

You get 2x the precision for about 9x the cost. The trick even composes: three doubles make a triple-double (~159 bits), four 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 a real library becomes the cheaper way to get digits.

Against the compiler's 128-bit float

If you write C or C++ on GCC or Clang, the quadruple precision from the naming note is built in, spelled __float128. Declare it and every operator works, 113 significand bits, no library, no license.
If it were fast, this article would be a lot shorter.

The catch is that no x86 or mainstream ARM processor implements binary128 in hardware, so the compiler lowers every operation into a call to a software routine (__addtf3 and friends).
That gives up 3 of the 5 advantages above: it is a call, the arithmetic leaves the FPU, and nothing vectorizes.
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 isolated operation it is 2x to 4x behind double-double, with one exception: division, double-double's weakest operation and soft-float's least bad one, goes to __float128 by 1.6x.
The long double row shows what hardware support is worth: nearly a double's speed, but only 64 significand bits, widening a double rather than doubling it.

The 6x perf price does buy you some things: 7 more bits, far more exponent range, correct rounding, and properly working infinities/NaNs.
And where binary128 is in hardware, IBM's POWER9 and later, __float128 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?).

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, with each type using the cheapest spelling of each operation it has. All of them are in section Appendix B: the benchmarked kernels, side by side; here is the double-double one:

1
2
3
4
5
6
7
8
9
10
11
12

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;
}

The pixelated areas in the double render are not a bug. Each one is a group of neighboring pixels whose coordinates rounded to exactly the same double, iterated to exactly the same result (color). 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.

Having 31 digits is not the same as using 31 digits

Everything so far measures one thing: how finely you can represent a point. How much precision the computation then needs is a different question, and the Julia set makes the gap between them 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 z instead of zero.
Same arithmetic, same cost per iteration, and it falls apart at a zoom the Mandelbrot renderer would not notice.0

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. 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.

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, which costs nothing; a double keeps all its digits no matter how small the number gets.
The damage is the very next step, where that tiny result is added to c, a number around 0.75.
Write one iteration out and you can count 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 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 again and still rounds to that same single step, which is why the two lines end identically.
Two pixels that started 3.1e-10 apart are indistinguishable before the second iteration begins.
Closer to the origin there is nothing left to lose at all. A pixel 5e-9 out has z, under half of that step, so 02 = 2.5 â 10-17c + z_0^2 returns 0.75 exactly and the pixel contributes nothing at all.
You can see that neighborhood in Figure 7, in the bottom right of the first picture, where every pixel falls to the same z. The double-double beside it draws a spiral there.1 = c

Neighboring pixels are now separated by the size of a rounding error rather than the size of a coordinate, and a chaotic map spends the rest of the run amplifying both. Across a row of the double render, 99% of the pixels come out with the wrong escape count.

That is also why the damage looks nothing like the earlier figures. There the input grid was quantized, so it failed as a lattice aligned to the pixel axes; here the grid is fine and the noise is inside the dynamics, so the image tears along the set's own structure instead.

The square sets the rate.
Zooming toward the origin, the pixel's signal in z shrinks as 1|z while the view width only shrinks as 0|2|z, so 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 penalty vanishes completely, matching the Mandelbrot depth for depth, which shows that this is a fact about small values of 0|z rather than about Julia sets.0

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 its bailout fires 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- Equalsand- CompareTodisagree about whether two NaNs are the same. Given the bullet above, that NaN is reachable.
Not correctly rounded.The guarantees are relative-error bounds (about- 2for the product), not "the nearest representable 106-bit value". Good enough for everything in the target band, but do not claim more.- -104
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 reassociate- sum - 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.
Non-associativity is worse than a double's, so tests have to compare against an exact oracle rather than against remembered literals.
~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 three, and it hands off at 1e-28.

What 31 digits buys, in real-life units

Double-double renders at zooms in the range of 10.
These numbers are quite hard to get a feel for, so here are some fun facts.29

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 landing on 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).

Takeaway

I sat down to write a short note about gluing two doubles together, and it kept growing. Partly because I decided to measure MPFR instead of guessing at it, and partly because the fractals I was using as a demo kept showing me things I had not planned to write about.

Three things are worth taking away even if you skipped everything above.

The performance gap is split almost exactly in half. A double-double costs about 9x a double, and an arbitrary-precision float (MPFR here) at the same 31 digits costs about 9x a double-double again.
That second 9x is why I wrote the type at all.
It kept my own explorer interactive at these depths and made high-resolution renders practical (until I implemented perturbation-based algorithms).

You cannot hold it wrong. It is two doubles in a readonly struct: nothing to preallocate for better performance, no native binary or P/Invoke, no precision field to keep making decisions about.
It inlines, and it vectorizes if you want it to.

Having 31 digits is not the same as using 31 digits. The type decides how many digits you start with. Your operations and numbers decide how many useful digits are still left at the end.
In the Julia section a single addition costs the same 14 digits either way: a double goes in with just under 16 and comes out with 2, a double-double goes in with just under 32 and comes out with 18.
A wider type hands you more digits to start with. It does nothing about how fast you lose them.

The decision rule I ended up with:

Need up to ~31 digits, in a hot loop, and want it fast without maintaining allocation discipline by hand? Double-double: two hundred lines of arithmetic and no dependencies.
Need arbitrary or input-dependent precision, or simply more than ~31 digits? Use a real library, and budget for roughly 80x a plain double at those same 31 digits, rising as you ask for more, and 300x or worse if you let the API allocate a result per operation.
Need range rather than precision? Different trick entirely. I have a type for that one too, and it may get its own article.

And one note on the motivating example: for deep Mandelbrot zooms, brute-force precision is not the state of the art. Perturbation goes where no fixed-width type can, and I have implemented that too, but it deserves an article of its own. The fractal was the motivation here, not the point. The gap between doubles and arbitrary-precision libraries shows up in far more places than fractal renderers.

Footnotes

The naming gets even muddier 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 for- double. â©
Printed values are rounded for readability. - Bis really the double- nearest0.555666777888, which drags a tail of extra digits from the decimal-to-binary conversion into- B,- 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
233
234
235
236

// 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;
///

/// 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, and comparisons 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).
///
public readonly partial struct DoubleDouble : IEquatable, IComparable {
private readonly double m_hi;
private readonly double m_lo;
/// The leading component; equals the value rounded to nearest double.
public double Hi => m_hi;
/// The trailing component (the rounding error of ).
public double Lo => m_lo;
/// True when the value is exactly zero.
public bool IsZero => m_hi == 0.0;
/// -1, 0 or +1 (0 also for NaN).
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;
}
/// Creates the exact value of a double; the trailing component is zero.
public static DoubleDouble FromDouble(double value) {
return new DoubleDouble(value, 0.0);
}
///
/// Creates the exact sum + of two arbitrary
/// doubles, renormalizing so the invariant holds.
///
public static DoubleDouble FromTwoDoubles(double hi, double lo) {
(double sum, double error) = twoSum(hi, lo);
return new DoubleDouble(sum, error);
}
/// Adds two values. Relative error stays below 3 * 2^-106 for all inputs.
public DoubleDouble Add(DoubleDouble right) {
// Add the two components pairwise, then fold the errors back in twice: the first
// fold 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);
}
/// Subtracts from this value.
public DoubleDouble Subtract(DoubleDouble right) {
return Add(right.Negate());
}
/// Multiplies two values.
public DoubleDouble Multiply(DoubleDouble right) {
// hihi 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);
}
/// Multiplies by a plain double - two cross terms fewer than the full product.
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);
}
///
/// Divides this value by . 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.
///
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));
}
/// Returns this value squared (cheaper than the general product).
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);
}
/// Returns the negated value (exact - negation never rounds).
public DoubleDouble Negate() {
return new DoubleDouble(-m_hi, -m_lo);
}
/// Returns the absolute value.
public DoubleDouble Abs() {
return Sign < 0 ? Negate() : this;
}
///
/// 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.
///
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);
}
/// Rounds to the nearest double.
public double ToDouble() {
return m_hi + m_lo;
}
/// Compares values; thanks to the normalization invariant this is simply
/// lexicographic on (Hi, Lo).
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;
/// Widening a double is exact, so it may happen implicitly.
public static implicit operator DoubleDouble(double value) => FromDouble(value);
/// Narrowing back to a double loses half the digits, so it must be asked for.
public static explicit operator double(DoubleDouble value) => value.ToDouble();
///
/// Knuth's two-sum: sum + error == a + b exactly, 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.
///
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);
}
/// 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.
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
}
///
/// Exact product: product + error == a * b exactly. The full product needs 106 bits;
/// the fused multiply-add computes ab - 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).
///
/// 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.
///
/// 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.
///
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:

// 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 {
///

Significant digits prints: enough to round-trip any
/// value whose two components are contiguous (see ).
private const int DEFAULT_DIGITS = 36;
/// Formats with significant digits.
public override string ToString() {
return ToString(DEFAULT_DIGITS);
}
/// Formats with the given number of significant decimal digits, correctly rounded.
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);
}
///
/// 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).
///
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('.');
}
///
/// 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
/// for every value whose components are contiguous -
/// a hand-built pair such as (1.0, 1e-300) carries more information than 36 digits can.
///
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);
}
/// The exact value of a double-double as a rational (the denominator is a power of two).
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);
}
/// The exact value of a finite double as a rational, straight from its bits.
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);
}
/// Rounds an exact rational to the nearest double, ties to even.
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, folding 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);
}
/// Turns "-1.25e-7" and friends into an exact rational.
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));
}
/// Rounds an exact rational to N significant digits and lays them out.
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 What it actually costs come from these, so here they are in full rather than as a claim. The pixel loop and the threading are the same for all three and are left out; what follows is the arithmetic, plus the per-worker scratch the MPFR kernel needs and the other two do not.

Read them against each other.

namespace DoubleDoubleSample.Mandelbrot;
///

/// 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.
///
/// 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.
///
/// 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.
///
/// 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.
///
public static unsafe partial class MandelbrotRenderer {
/// Bailout radius squared. Generous, so the smooth-iteration formula is accurate.
private const double ESCAPE_SQUARED = 65536.0 * 65536.0;
/// Plain double, for reference and for the pictures that fall apart.
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;
}
///
/// The same source, with DoubleDouble substituted for double. 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.
///
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;
}
///
/// 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.
///
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);
}
///
/// Per-worker MPFR values. MPFR values cannot be shared across threads for writing, so
/// Parallel.For's localInit hands each worker its own set and localFinally releases them.
///
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;
/// The bailout radius as an MPFR value, so the escape test never has to convert.
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);
}
/// Every value is live here: both kernels restore that before returning.
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 does differ slightly, but it's just a small optimization resulting in no precision loss.
The orbit past |z| > 2 is diverging exponentially, so the low part of double-double has no real effect on the number of iterations.

The fourth kernel, the MPFR one with a fresh value per operation, is very similar to the third, but with each operation producing a new value. That's the 38x row.