Double-double: 31 digits of precision without leaving the FPU

Posted 16 Aug 2026By Marek Fiser25 min read

If you ever need more precision than what 15 decimal digits of the double format can offer, and arbitrary-precision arithmetic is too slow, too heavy, or simply not available, there is a neat trick: glue two doubles together and treat them as one number. This gives you ~31 decimal digits for roughly 9x the cost of a plain double in a real kernel (4-12x per isolated operation). With no heap allocation and no dependencies, this lands it almost exactly halfway between a double and an arbitrary-precision library perf-wise. This post explains the error-free transformations that make it work, measures it against MPFR, and shows where the trick runs out of steam.

The gap nobody fills

Floating-point precision comes in many sizes. For example, a double gives you about 15 decimal digits, essentially 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 is that gap in one picture: the same view rendered twice, blocky on the left where a plain double has run out, sharp on the right where a double-double has not.

The Mandelbrot set is a famous fractal with these mesmerizing patterns that I have explored before. A deep zoom runs out of precision in the most literal way possible: eventually two neighboring pixels land on the same double, and the image stops being a picture of the fractal and starts being a picture of the number format, as you can see in Figure 1.

Figure 1: Two deep locations rendered in double precision: whole blocks of pixels share one coordinate and come out the same color. The gap between neighboring doubles grows with the number's size, so the real axis is the coarse one on the left (re -0.74, im 0.13, wide blocks) and the imaginary axis is the coarse one on the right (re -0.11, im 0.92, tall blocks).

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 genuinely 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 that double-double => MPFR gap is closer to 40x instead of 9x.

The insight needed to unlock double-double precision 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 the number is. The exponent decides where the digits sit, it does not change how many there are. So take these two numbers:

(1)A = 111222333444

(2)B = 0.555666777888

Each of them has 12 significant digits, so each fits in a double with room to spare. Their exact sum, 111222333444.555666777888, has 24 digits, and a double holds about 15: evaluate A + B and the tail of B is quietly rounded away to 111222333444.55566.

But what if we simply 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 digits the first one has no room for. Nothing was rounded, because nothing was added.

That is the whole representation:

A double-double value stores the unevaluated sum of two doubles.

(3)x = x_hi + x_lo

with the invariant that x_hi is exactly what x_hi + x_lo rounds to as a double. The high part carries the value, the low part carries the error of the high part. Keeping that invariant true is what every operation in the next section is really doing.

Code listing 1: The whole double-double representation. Everything else in this post is arithmetic on these two fields.
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;

}

One naming note before the math. The 128-bit type in IEEE 754, the standard that also defines the double you already use, is binary128, usually called quadruple precision: one contiguous 113-bit significand and an exponent that reaches ~1e4932.

A double-double is also 128 bits of storage, but it is two ordinary doubles: ~106 significand bits and a normal double's exponent range, still overflowing at ~1e308. It buys significant digits, not range.

The two names collide in the wild. On IBM POWER, long double is double-double, exactly this technique (GCC calls it __ibm128, and reserves __float128 for the IEEE type). On x86, long double is a third thing again, and which thing depends on the compiler: 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. So "128-bit floating point" means quadruple precision, double-double, or neither, depending on who is speaking.

Figure 2: Where the significand bits sit in each format.

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 and then has to round it. 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.

Code listing 2: Knuth's two-sum: sum + error == a + b exactly, for any two doubles.
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:

Code listing 3: Knuth's two-sum on the two example numbers.

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.^[1] 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: 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 not an approximation of A + B, it is A + B, to the last digit. The arrangement is what earns that, and it is not the rounding 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. That is a theorem, not an accident; 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, and it is what knowing the order buys you:

Code listing 4: Dekker's fast two-sum: the same guarantee as Knuth's two-sum in half the operations, valid only when |a| >= |b|.
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:

Code listing 5: Exact product via fused multiply-add. This one instruction is why modern double-double is cheap.
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^-106, and finish with a quickTwoSum 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.

Code listing 6: Addition. Relative error stays below 3 * 2^-106 for all inputs.
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

Code listing 7: Multiplication: one exact product, plus the three cross terms.
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 sits at 2^-106 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.

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.

Code listing 8: Division by iterated remainder.
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

Code listing 9: Comparison is trivial: the normalization invariant pays for itself here.
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, and the invariant is what 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 where they sit relative to each other, 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 one double holds exactly is the pair with lo = 0; that is all FromDouble(1.5) does.
Touching: lo continues exactly where hi ends: 106 contiguous bits, the "31 digits" from the title. The smallest example is 1 + 2^-53, which becomes the pair (1, 2^-53).
Gapped: This is the interesting case, where lo starts 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, hi would no longer be the value rounded to a double and the same number could be written as many different pairs. That is precisely 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 error carries the rest:

Code listing 10: B's binary digits split between the two halves of the pair. A owns all the integer bits, so the fraction of the exact sum is B itself.

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^-3 are the same number.

The gapped shape is the one place that needs a warning. The gap itself costs nothing and loses nothing: the pair holds exactly two "islands" of bits. Suppose this absurd case: 10³⁰⁰ + 10^-300 is a perfectly legal pair with 1,940 zero bits between its islands, carried in 16 bytes. However, evaluate (10³⁰⁰ + 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 10³⁰⁰ + 1, the lowest island is gone.

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

What double-double promises is ~106 significant bits at the leading edge of the value; a wide gapped pair exceeds that promise only until the next operation snaps it back down. Printing hits the same wall from the other side, which is why the type also carries a ToStringExact.

What it actually costs

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

Two things stand out. 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 anyone actually wants is against MPFR, and the honest place to make it is a real kernel rather than a loop that repeats one operation.

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

(4)z_n+1 = z_n² + c

If 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 that change every time, the escape test included, roughly 390 million iterations, 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. No law says the option in the middle lands in the geometric middle of the gap it fills. 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.^[2]

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 only by you, on every line, forever: 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. Lean 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 double adds 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, no "how wide should this result be", 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 exactly 2x the precision for about 9x the cost, and not a digit more. The trick even composes: three doubles make a triple-double (~159 bits), four a quad-double (~212 bits). But each extra term buys the same 53 bits while the renormalization bill keeps growing, so the economics quietly run out, and somewhere around the third or fourth a real library becomes the cheaper way to buy digits.

Showcase: a Mandelbrot zoom that double cannot reach

Back to the zoom that started this. Here is what the renderer is asked to draw:

Code listing 11: A deep-zoom center in the seahorse valley, and a view width that is past double's floor.
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:

Code listing 12: The double-double inner loop. Note the escape test: it decides a branch, so a double's accuracy is plenty for it.
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

Figure 5: The 7.7e-15 view. The double render quantizes into blocks; the double-double render is unbothered.

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. The double-double precision resolves this issue and recovers all the beautiful details 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

Figure 6: A nearby center at 4.9e-31 across, sixteen orders of magnitude deeper. Now it is the double-double precision that has run out. The reference on the right is the same view computed by perturbation, which is what is really there.

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.

double
double-double

Figure 7: One Julia view just about 3e-7 across. The Mandelbrot at the top of this article is more than seven orders of magnitude deeper and comes out clean.

The coordinate is not what ran out. The center sits 9.7e-8 from the origin, and doubles are dense down there: consecutive ones are 1.3e-23 apart.

What runs out is the first line of the loop, and it is the addition that does it, not 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:

Code listing 13: One pixel's first iteration, 1.2e-8 from the origin, and its neighbor one pixel step away. Both squares land on the same double once c is added.

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: the pixel's contribution is stored 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₀² = 2.5 ∙ 10^-17, under half of that step, so c + z_0^2 returns 0.75 exactly and the pixel contributes nothing whatsoever. You can see that neighborhood in Figure 7, in the bottom right of the first picture: every pixel in that area falls to the same z₁ = c. The double-double beside it draws a spiral there.

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 is what sets the rate. Zooming toward the origin, the pixel's signal in z₁ shrinks as |z₀|² while the view width only shrinks as |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 tax vanishes completely, matching the Mandelbrot depth for depth, which is the tell that this is a fact about small values of z₀ rather than about Julia sets.

Where the double-double trick stops working

Every trick has edges. These are double-double's:

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 lands below every real value rather than making every comparison false the way a double's does, and Equals and CompareTo disagree 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 2^-104 for the product), not "the nearest representable 106-bit value". Good enough for everything in the target band, but do not claim more.
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-math will cheerfully reassociate sum - a into oblivion, 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 and no amount of cleverness changes that. In my own fractal renderer double-double is exactly one tier of three, and it hands off at 1e-28.

What 31 digits buys, in real-life units

Double-double is able to render at zooms in the range of 10²⁹. 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 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 by zooming that deep, you have only scratched the surface. 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).

A view 2.5e-398 across, 1,322 doublings below the whole set, computed by perturbation.
Placing its center takes about 400 decimal digits. A double has 15 and a double-double has 31.

Figure 8: A view 2.5e-398 across, 1,322 doublings below the whole set, computed by perturbation.^[3] Placing its center takes about 400 decimal digits.^[4] A double has 15 and a double-double has 31.

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. Rendering 9x faster can turn an interactive fractal explorer from unusable into usable, and it lets you produce much higher-resolution images in reasonable time.

You cannot hold it wrong. It is two doubles in a readonly struct: no destination to preallocate, no temporaries to hoist, no allocation discipline to defend in code review. No native binary per platform, no 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 burns the same 14 digits either way: a double walks in with just under 16 and comes out with 2, a double-double walks 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 honest 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

Printed values are rounded for readability: B is really the double nearest 0.555666777888, which drags a tail of extra digits from the decimal-to-binary conversion into B, bLost, error and 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. ↩
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:

Code listing 14: DoubleDouble.cs, complete.
// 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, 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).
/// </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;

	/// <summary>True when the value is exactly zero.</summary>
	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>Creates the exact value of a double; 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>Adds two values. Relative error stays below 3 * 2^-106 for all inputs.</summary>
	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);
	}

	/// <summary>Subtracts <paramref name="right"/> from this value.</summary>
	public DoubleDouble Subtract(DoubleDouble right) {
		return Add(right.Negate());
	}

	/// <summary>Multiplies two values.</summary>
	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>Multiplies by a plain double - two cross terms fewer than the full product.</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>Returns this value squared (cheaper than the general product).</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>Returns the negated value (exact - negation never rounds).</summary>
	public DoubleDouble Negate() {
		return new DoubleDouble(-m_hi, -m_lo);
	}

	/// <summary>Returns the absolute value.</summary>
	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>Rounds to the nearest double.</summary>
	public double ToDouble() {
		return m_hi + m_lo;
	}

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

Code listing 15: DoubleDoubleFormat.cs, complete.
// 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, 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);
	}

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

Code listing 16: Kernels.cs, complete: one kernel, three ways.
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>
	/// 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 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.

blog

C#Performance Fractals

This post is licensed under CC BY 4.0 by the author.