Double-double arithmetic: 31 digits of precision from two doubles
If you ever need more precision than what 15 decimal digits of the double format can offer, there is a neat trick: glue two doubles together and treat them as one number. This gives you ~31 decimal digits for roughly 9x the cost of a plain double. With no heap allocation and no dependencies, this puts it perf-wise almost exactly halfway between a double and an arbitrary-precision library. This post explains how the format works and the math behind it. It also benchmarks the type against __float128 and MPFR, and shows where the trick runs out of steam.
1. The gap that double-double fills#
Floating-point types come in many "flavors". For example, a float has ~7 decimal digits, a double ~15. An arbitrary-precision library gives you as many digits as you want but at a painful per-operation cost (~100x slower or more). Between "not quite enough precision" and "orders of magnitude slower" there is a gap, and I fell into it implementing an interactive Mandelbrot set renderer. The image at the top of this page was my motivation. It's the Mandelbrot set rendered using double precision (left) and double-double precision (right).
The Mandelbrot set is a famous fractal, full of endlessly repeating spirals and mini copies of itself, that I have explored before. When zooming deeper and deeper, eventually the image starts looking like a map of the number format, not a fractal.
Notice that the flat-colored blocks in Figure 1 are oriented differently in each image.
This is a consequence of how a floating point representation works.
The gap between neighboring numbers grows with the magnitude, so the larger of the two coordinates loses resolution first.
The left image is at (-0.74, 0.13) so the blocks are wider than tall, but the right image is at (-0.11, 0.92) where the ratio is reversed.
The canonical answer to "I need more precision than double" is to use a library for arbitrary-precision math, typically GMP or MPFR. But it comes with some caveats:
- Each value is a pointer to a heap-allocated data array (called limbs), so creating a new number allocates. One has to be very careful to reuse values to avoid allocations on every operation.
- Most operations loop over the limb array, which means memory traffic, data-dependent branches, etc.
- In .NET there are even more hurdles. GMP and MPFR are native C libraries, so you are looking at P/Invoke, a native binary per platform to ship, and marshalling on the boundary.
So, to get a little more precision cheaply, the trick is:
Keep two doubles separate, but treat them as one number.
This gives you ~31 decimal digits with no heap allocations and no dependencies.
The idea is not mine. Dekker described it in 1971, and the pieces it is built from are even older (see Further reading for details). I wanted to prove to myself (and to you) that this trick is still relevant today, after half a century of faster hardware and better arbitrary-precision libraries.
2. The double-double idea#
A floating-point type holds the same number of significant digits no matter how large or small the number is. The exponent then determines where the decimal point will be. Take these two numbers:
A = 111222333444B = 0.555666777888Each of them has 12 significant digits, so each fits in a double.
Their exact sum is 111222333444.555666777888 and has 24 digits, but a double can only hold around 15.
If you evaluate A + B, the digits that do not fit are rounded away, leaving you with just 111222333444.55566.
But what if we do not evaluate the addition? If we keep A and B side by side and agree to treat the pair as one number, then one of them carries the leading digits and the other carries the rest.
A double-double value stores the unevaluated sum of two doubles:
x = xhi + xlowith the invariant that the high part xhi is exactly what xhi + xlo rounds to as a double.
The low part xlo carries the rest (the "error" of the high part after addition).
Every operation in the next section keeps this invariant true.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
public readonly partial struct DoubleDouble : IEquatable<DoubleDouble>, IComparable<DoubleDouble> {
private readonly double m_hi;
private readonly double m_lo;
/// <summary>
/// The leading component; equals the value rounded to nearest double.
/// </summary>
public double Hi => m_hi;
/// <summary>
/// The trailing component (the rounding error of <see cref="Hi"/>).
/// </summary>
public double Lo => m_lo;
}
A double-double has 128 bits, but it is not quadruple precision.
The 128-bit type in IEEE 754 (the standard that also defines the double) is called binary128, better known as quadruple-precision.
It has a 113-bit significand and a 15-bit exponent (that reaches to ~1e4932).
A double-double also uses 128 bits of storage, but it is just two ordinary doubles next to each other, resulting in 106 combined significand bits and an unchanged exponent range (still overflowing at ~1e308).[1]
Figure 2: Number of significand bits of each format.3. Error-free transformations#
It all builds on one observation:
When a double-arithmetic operation rounds, the amount it rounded away can still be recovered as another double.
When you add two doubles, the hardware computes the true sum, which is then rounded.
Let the rounded result be s = round(a + b).
The part that got rounded away is exactly e = (a + b) - s.
And once we have it, we are done because a + b = s + e, and (s, e) is our double-double!
But how to get e if we are unable to compute a + b exactly? The algorithm is called Knuth's two-sum and you can see it in Code listing 2: six additions/subtractions with no bit twiddling and no wider types needed.
1
2
3
4
5
6
7
8
private static (double Sum, double Error) twoSum(double a, double b) {
double sum = a + b;
double bKept = sum - a; // The part of b that made it into the sum.
double aKept = sum - bKept; // The part of a that made it into the sum.
double aLost = a - aKept;
double bLost = b - bKept;
return (sum, aLost + bLost);
}
Here is the algorithm on our example numbers A and B:
1
2
3
4
5
6
7
8
9
10
11
12
13
A = 111222333444
B = 0.555666777888
A + B, exactly = 111222333444.555666777888
sum = A + B = 111222333444.5556640625 <- what the FPU returns
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, the rest was rounded away, as expected.[2]
But the bLost is exactly what was lost!
In the extreme case, when the operands are too "far apart" in magnitude to overlap at all, the sum rounds to the larger of A or B and the error is the smaller one.
Notice that A lost nothing in the example above (aLost is 0) and that is not a coincidence.
Whenever |a| >= |b| the sum keeps every bit of the larger operand, so aKept and aLost can be omitted 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 lost nothing, so only b's loss remains.
}
In multiplication we need the exact product of two 53-bit numbers, which will need 106 bits. So we split it into a high and a low part. Dekker's original method (Veltkamp splitting) took 17 operations. Once the hardware has an FMA (fused multiply-add), which every mainstream processor does, it takes just two!
1
2
3
4
5
private static (double Product, double Error) twoProduct(double a, double b) {
double product = a * b;
double error = Math.FusedMultiplyAdd(a, b, -product);
return (product, error);
}
The trick is that an FMA computes a*b + c with a single rounding at the very end.
So we can call it as fma(a, b, -product) which evaluates a*b - product and the subtraction cancels everything product kept, leaving exactly the error that the initial a * b rounded away.
These three functions, twoSum, quickTwoSum, and twoProduct, are called error-free transformations, and they are the foundation of double-double arithmetic.
They all assume round-to-nearest arithmetic and no overflow.
Note that they do not handle infinities and NaNs properly.
4. Building arithmetic operations for double-double#
Now that we have twoSum, quickTwoSum, and twoProduct, writing arithmetic operations for DoubleDouble is quite straightforward.
We have to keep in mind to use an error-free transformation wherever a rounding would lose digits,
and finish with a quickTwoSum so the result is properly normalized.
Addition#
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);
}
Add the two components pairwise, then incorporate lowSum and lowError one by one.
Addition of the lowError itself can round, so we need to use an extra quickTwoSum.
There is a cheaper "sloppy" variant in circulation using only 11 operations instead of these 20,
which adds both low parts m_lo and right.m_lo in one step.
It is fine while the operands share a sign, but if they have opposite signs and cancel each other,
the relative error has no upper bound.
Multiplication#
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);
}
Multiply the two high parts exactly with twoProduct, then add the three cross terms to the error.
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 what plain double arithmetic rounds away from them is around 2-106 of the product, the last bit the pair can hold.
Because of this rounding, the multiplication is not correctly rounded.
Measured against an exact product on 20,000 random pairs, it returns the nearest representable pair about two thirds of the time, and is within a few units of the last bit otherwise.
Squaring can be optimized by using just 2 cross terms instead of 3,
and multiplying by a plain double needs only one.
These are worth having as overloads for better performance.
Division#
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));
}
There is no error-free transformation for division.
Instead, we take a double's worth of the quotient, subtract its contribution at double-double precision,
and repeat as long division 3 times.
Square root uses the same idea. Take a double root, then one step of Newton refinement. Division and square root come out accurate to a couple of units in the last place but not correctly rounded (not the nearest representable pair). The error starts to matter only when the last bits have to be reproducible against a different implementation.
Comparison#
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);
}
The comparison is easy thanks to our double-double invariant.
Note that if unnormalized pairs were allowed, (1, 0.75) and (1.75, 0) would be the same number and compare as different ones.
Three possible shapes of double-double#
Until now you probably imagined that double-doubles are just two doubles tightly "glued" together into 106 contiguous bits. However, it does not have to be that way as shown in Figure 3.
Figure 3: The three possible shapes of a double-double pair, and the overlapping arrangement that breaks the invariant.- Single double: Any value that a
doublecan hold exactly is a pair withlo = 0. - Touching:
locontinues exactly wherehiends, giving 106 contiguous bits. For example,1 + 2-53becomes(1, 2-53), exactly. - Gapped: This is the interesting case, where
lostarts "further down", and every bit position in between is an implicit zero. This may feel 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 invariant would be broken. For example,FromTwoDoubles(1.0, 0.75)renormalizes and returns(1.75, 0).
The gap of 2 from our example can be seen if we write the digits out in binary.
1
2
3
4
5
6
7
A = ...0000100. <- 37 integer bits, the top one at 2^36
B = .1000111001000000001011011000...
sum = ...0000100.1000111001000000 <- 53-bit window, ends at 2^-16
error = 001011011000... <- first set bit at 2^-19
^^
the gap: two zero bits stored by neither half
In our case, the error term happens to begin with two zeros, so the gap naturally emerges.
Those zeros are real digits of the value, but a double does not store leading zeros.
Instead, its exponent points lower (0.001011 and 1.011 ∙ 2-3 are the same number).
The third shape with a gap has one catch.
It may seem like you are getting extra precision "for free", and you are, but the zeros in the gap may collapse on any operation and the extra precision can be lost.
Suppose this absurd case: 10300 + 10-300 is a perfectly legal pair with 1,940 zero bits in the gap.
However, evaluate (10300 + 10-300) + 1 and the exact answer collapses to simply 10300 + 1, losing the 10-300 part entirely.
"But surely 10300 written in binary does not have loads of trailing zeros‽"
Right, it does not. Writing it exactly needs about 700 significant bits.
But the 1e300 in the example 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.
Now you may be wondering why not glue together 3 doubles, or 4?
It works too! 3 doubles make a triple-double (~159 bits), 4 a quad-double (~212 bits, the QD library by Hida, Li and Bailey in Further reading). But each extra term adds another 53 bits while renormalization keeps getting more expensive, so somewhere around the fourth double the precision stops being worth the cost.
5. Double-double operations benchmark#
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 what we are usually used to (thanks to the FMA). And division is the most expensive one for both types.
Comparison with MPFR#
Now let's compare DoubleDouble against MPFR, an arbitrary-precision library.
The code I am comparing is a Mandelbrot set kernel that iteratively computes zn in the complex plane.
zn+1 = zn2 + cIf |zn| stays bounded forever, c belongs to the set and the pixel is set to black.
If it "escapes" to infinity, which provably happens once |z| > 2, the iteration count at that point determines the color.
One image totals roughly 390 million iterations.
MPFR goes through hand-written P/Invoke with no wrapper class, at 106 bits of precision, which matches double-double's.
| Kernel | iterations/second | vs double | vs DoubleDouble |
|---|---|---|---|
double | 8,700 M | 1x | |
DoubleDouble | 1,000 M | ~9x | 1x |
| MPFR @106, in place | 108 M | ~81x | ~9x |
| MPFR @106, fresh result per operation | 26 M | ~330x | ~38x |
The two MPFR styles run the identical operations, but the "fresh result" does not reuse values and allocates per operation, the way many wrappers present it.
As you can see, double-double costs roughly 9x a double, and MPFR costs roughly 9x a double-double again.[3]
The slowdown in the "fresh result" row is avoidable, but it is up to the code writer to keep the discipline and properly reuse values (e.g. never write a * b + c as an expression because each operation would allocate a new value).
A DoubleDouble has no heap allocations, so you can write zRe * zIm * 2.0 + cIm as usual and it allocates nothing.
The 9x is the number for my particular case, and it comes from 7 arithmetic operations per iteration, so your mileage with double-doubles may vary.
Compared to arbitrary-precision floats, the DoubleDouble has some advantages:
- A double-double add or multiply is just 7-20 instructions, which will often be inlined, unlike MPFR, which is a function call (with a P/Invoke boundary on top in .NET).
- The instruction sequence is fixed and branch-less, which will make your CPU pipeline happy, unlike MPFR's data-dependent loops over arrays.
- It can be vectorized and ported. It is possible to rewrite double-double operations using
Vector256<double>that process 4 values at a time, and the same arithmetic ports to a GPU. While this is a rewrite rather than something the compiler/JIT does for you, such a rewrite would not be possible with an arbitrary-precision library.
Comparison with __float128#
If you write C/C++ on GCC/Clang, a quadruple-precision __float128 is built in and it has all the arithmetic operations available.
The catch is that no x86 or mainstream ARM processor implements __float128 in hardware, so the compiler emulates every operation.
Measured with the same chained loops and the same kernel, this time in C[4]:
| Kernel (C, one core) | iterations/second | vs double | vs DoubleDouble |
|---|---|---|---|
double | 531 M | 1x | |
long double (x87) | 441 M | 1.2x | |
DoubleDouble | 76 M | ~7x | 1x |
__float128 | 12.4 M | ~43x | ~6x |
Per operation __float128 is 2-4x slower than double-double (chained loops, not shown), but the entire Mandelbrot kernel runs 6x slower.
The long double row shows an additional data point for an FPU-accelerated type but with only 64 significand bits.
The 6x slowdown of __float128 does buy you some things: 7 more bits of precision, far more exponent range, correct rounding, and properly working infinities/NaNs.
On platforms where __float128 is hardware accelerated, e.g. IBM's POWER9 and later, use it!
On x86 it is not worth it.
In .NET the question never arises since there is unfortunately no 128-bit float type.
6. Showcase: a deep Mandelbrot zoom#
Now let's take a look at some Mandelbrot renders that inspired this work. 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;
Code listing 12 shows the double-double kernel in C#, and the rest are in Appendix B: the benchmarked kernels.
1
2
3
4
5
6
7
8
9
10
11
12
13
DoubleDouble zRe = DoubleDouble.Zero;
DoubleDouble zIm = DoubleDouble.Zero;
for (; n < maxIterations; n++) {
DoubleDouble zReSquared = zRe.Square();
DoubleDouble zImSquared = zIm.Square();
magnitudeSquared = zReSquared.Hi + zImSquared.Hi;
if (magnitudeSquared > ESCAPE_SQUARED) {
break;
}
zIm = zRe * zIm * 2.0 + cIm;
zRe = zReSquared - zImSquared + cRe;
}
Figure 4: A view 1e-3 across. At this depth both types render the same image.
double 
double-double
The pixel step here is about 1.5e-17, a seventh of the gap between neighboring doubles, so the blocks are about seven pixels wide, and every further zoom step widens them until the whole frame is a single one. With double-double the blocks are gone, and the detail holds down to a view width of about 1e-28.
Keep going, though, and the same thing happens to it, as Figure 6 shows.
double-double 
perturbation, for reference
7. The Julia set problem#
So far, all we measured was how finely you can represent a point. But how much precision the computation is left with is a different question, and the Julia set made it very clear to me.
A Julia set uses Equation 4 again, but with the two roles swapped: c is one constant for the whole image, and the pixel coordinate is the starting z0 instead of zero.
So the arithmetic and the cost per iteration are the same, but the image falls apart earlier than expected.
I ran into this while browsing Julia sets in my own renderer. I often found that the most interesting structure is in the center of the image (the origin), so that's where I zoomed in. The glitches arrived so early that I first thought I had a bug in my code. They also looked very different. However, switching the kernel to double-double eliminated them completely, so the bug was not in the code, and I went looking for what had run out.
double 
double-double
The issue here is not with the "density" of doubles like in our earlier example.
The precision runs out in the first iteration of the loop, in the addition rather than the squaring.
Squaring makes the pixel's coordinate smaller, but a double keeps all the important digits.
The digits are lost at the very next step, where that tiny squared result is added to c, a number around 0.75.
Here is an example:
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
From the example you can see that two pixels that started 3.1e-10 apart are indistinguishable before the second iteration begins.
You can see that neighborhood in Figure 7, in the bottom right of the first picture where every pixel falls to the same z1 = c, shown as a large gray region.
The double-double next to it draws a spiral there.
Also notice that the glitches look different. Instead of axis-aligned boxes you can see jagged tears coming from a rounding error amplified by the iteration dynamics of the system.
Around the origin of the Julia set, the usable zoom is the square root of what the same type manages on a Mandelbrot, half the digits. But as you move away from the origin, the glitch disappears.
8. Zoom depth in real-life units#
Double-double renders at zooms in the range of 1029.
These numbers are quite hard to get a feel for, so here are some fun facts.
Put the whole Mandelbrot set on your phone screen, about 7 cm wide, some 1200 pixels across. Now pinch to zoom and keep going until neighboring pixels start rounding to the same coordinate, the way they do in Figure 5. Then ask how big the whole set has grown by the time that happens.
Say every pinch doubles the zoom and takes one second.
In double precision you are done in 44 seconds, with the screen down to 1.3e-13 of a unit.
The entire set drawn at that scale would be 1,200,000,000 km across, 8x the distance from the Earth to the Sun, reaching from here to past Jupiter.
In DoubleDouble you keep pinching for another 53 seconds, down to 1.5e-29 across.
The original set is now 12 observable universes wide!
And it goes far, far deeper! 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.
Figure 8: 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.[5]9. Limitations of double-double#
Double-double has the following limitations:
- The representable range is unchanged, and the usable range is slightly 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 are not properly supported. The error terms compute
inf - inf, so the first overflow or division by zero returns NaN where a double would result in an infinity. The Mandelbrot kernel above is safe only because it stops iteration long before anything can grow that large. - NaN breaks sorting. Comparison is lexicographic on the two halves, so a NaN sorts below
every real value rather than making every comparison false the way a double's does, and
EqualsandCompareTodisagree about whether two NaNs are the same. - The operations are not correctly rounded. The guarantees are only about relative errors (e.g.
2-104for the product). Without correct rounding, reordering an expression (e.g. using associativity rules) alters the result more than it would with a plaindouble. This can complicate testing. - The compiler must not "help". Every error-free transformation depends on the operations happening
exactly as written, at exactly double precision. For example, C and C++ with
-ffast-mathwill reassociatesum - aaway. Fortunately, .NET (RyuJIT) is not doing any optimizations on floating-point expressions.
10. Conclusion#
Here are all the formats discussed in this article:
| Format | Decimal digits | perf vs. double | Pros | Cons |
|---|---|---|---|---|
float | ~7 | ~1x | half the memory, twice the SIMD lanes | low precision |
double | ~15 | 1x | fast, precise | precision may run out |
long double (x87) | ~19 | ~1.2x | extra precision for free | x86 and GCC/Clang only, no SIMD, only 4 more digits |
DoubleDouble | ~31 | ~9x | no heap, no dependencies, inlines, vectorizes | infinities break, not correctly rounded, ~200 lines you own |
__float128 | ~34 | ~40x | high precision and range up to 1e4932 | emulated in software on x86, limited support |
| MPFR | as many as you ask | ~80x and up | unlimited precision | heap per value, P/Invoke, allocation discipline needed |
I knew that double-double is faster than MPFR, but it surprised me that it happens to be right in the middle. The speedup over MPFR is why I wrote the type in the first place. It kept my own explorer interactive at extended zoom levels and made high-resolution renders practical.
Double-double is a very niche type, but despite that, I hope that I convinced you that it can be useful. For example, if you need your code to run in an environment without access to arbitrary-precision libraries like a GPU or a web browser, or when just a little more precision is sufficient and performance is critical.
And finally I should mention that for deep Mandelbrot zooms, brute-force precision is not the state of the art. Perturbation theory enables very deep zooms without the need for using arbitrary precision for all computations, but more on that some other time. The fractal was the motivation here, not the point. If you know about neat uses for double-double, or end up using them yourself, let me know!
Figure 9: And a little something extra at the end.
A view 3.9e-811 across, 2,693 doublings below the whole set, computed by perturbation.
And we are still not at the bottom!11. Further reading#
- T. J. Dekker, "A Floating-Point Technique for Extending the Available Precision", Numerische Mathematik, 1971. The original paper. Both the splitting trick and the fast two-sum come from here.
- D. E. Knuth, The Art of Computer Programming, vol. 2: Seminumerical Algorithms, section 4.2.2.
The six-operation
twoSumand the proof that it is exact. - Y. Hida, X. S. Li and D. H. Bailey, "Algorithms for Quad-Double Precision Floating Point Arithmetic", 2001. The paper behind the QD library, the reference double-double and quad-double implementation.
- M. Joldes, J.-M. Muller and V. Popescu, "Tight and Rigorous Error Bounds for Basic Building Blocks of Double-Word Arithmetic", ACM TOMS, 2017. Proven error bounds for the addition, multiplication and division used here.
- J. R. Shewchuk, "Adaptive Precision Floating-Point Arithmetic and Fast Robust Geometric Predicates", 1997. Extends the technique to expansions of n doubles whose precision adapts to the input.
Footnotes
The naming gets more confusing with C's
long double, whose meaning depends on the platform and the compiler. On IBM POWER it is actually double-double (GCC calls it__ibm128, and reserves__float128for the IEEE type). On x86 it is yet another type: GCC and Clang give you the 80-bit x87 extended format with a 64-bit significand, while MSVC makes it just another name fordouble. ↩Printed values are rounded for readability. The numbers in decimal form would have extra digits from the decimal-to-binary conversion, but the algorithm carries them exactly as shown anyway. ↩
Medians of three runs on an idle machine, and the spread between runs is a few percent. The resulting iteration counts were the same across all runs. ↩
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. ↩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. So here is the exact coordinate (406 digits each). Real part:
-1.2546501186909057505269466250921981644070768100251708177516070303734046805262190991942188750230768678605788326475221138767425434553956510745751843244663895032007841539877137705604003669583875944006181706618704710525415050302476211752688927355924699367031357139525690997877831869648195623522823691538908055742989329820011761064778120010166893341652782500699958380682058425509820660627053365326418384253599781. Imaginary part:0.3819035642800001819461581237720802452234485990147800334446378567314132368999769468860848762828896508632064176850184911749715609976648123007802564910929725517436736759042462717235720004499322510654559848620036722449528324306932522044048459018966997955517496317108858725278768456014781914592217389991723653848006887239802789886568674879888329442270299742538205772252794459812326617996073246041743419429748869. The viewport is 1.906e-398 high. ↩
Appendix A: the whole type#
The complete type is in two files.
The first file contains the arithmetic and depends only on System.
This is the part to copy if 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
// Double-double arithmetic in about 200 lines. This file depends only on System so that it
// can be copied into any project. Printing and parsing use BigInteger and are in
// DoubleDoubleFormat.cs.
namespace DoubleDoubleSample;
/// <summary>
/// A value stored as the unevaluated sum of two doubles: ~106 significand bits
/// (~31 decimal digits) at a few times the cost of a double, with no heap allocation.
///
/// Invariant: Hi == RoundToNearest(Hi + Lo), i.e. |Lo| is at most half an ulp of Hi.
/// Every factory and operation below maintains it; comparisons and ToDouble rely on it.
///
/// Operations are not correctly rounded but stay within a few units of 2^-106 of the
/// exact result. The exponent range is the same as a plain double's, and infinities are
/// not special-cased: the first operation on one returns NaN (the error terms compute
/// inf - inf).
///
/// By Marek Fiser, from marekfiser.com/blog/double-double-arithmetic/ (CC BY 4.0).
/// </summary>
public readonly partial struct DoubleDouble : IEquatable<DoubleDouble>, IComparable<DoubleDouble> {
private readonly double m_hi;
private readonly double m_lo;
/// <summary>The leading component; equals the value rounded to nearest double.</summary>
public double Hi => m_hi;
/// <summary>The trailing component (the rounding error of <see cref="Hi"/>).</summary>
public double Lo => m_lo;
public bool IsZero => m_hi == 0.0;
/// <summary>-1, 0 or +1 (0 also for NaN).</summary>
public int Sign => m_hi > 0.0 ? 1 : m_hi < 0.0 ? -1 : 0;
public static DoubleDouble Zero { get; } = new DoubleDouble(0.0, 0.0);
public static DoubleDouble One { get; } = new DoubleDouble(1.0, 0.0);
private DoubleDouble(double hi, double lo) {
m_hi = hi;
m_lo = lo;
}
/// <summary>Exact conversion from a double.</summary>
public static DoubleDouble FromDouble(double value) {
return new DoubleDouble(value, 0.0);
}
/// <summary>The exact sum <paramref name="hi"/> + <paramref name="lo"/> of two arbitrary doubles.</summary>
public static DoubleDouble FromTwoDoubles(double hi, double lo) {
// Renormalizes so that the invariant holds.
(double sum, double error) = twoSum(hi, lo);
return new DoubleDouble(sum, error);
}
public DoubleDouble Add(DoubleDouble right) {
// Add the two components pairwise, then incorporate lowSum and lowError one by one.
// Adding lowError can round, so a second quickTwoSum renormalizes the result.
(double sum, double error) = twoSum(m_hi, right.m_hi);
(double lowSum, double lowError) = twoSum(m_lo, right.m_lo);
error += lowSum;
(sum, error) = quickTwoSum(sum, error);
error += lowError;
(sum, error) = quickTwoSum(sum, error);
return new DoubleDouble(sum, error);
}
public DoubleDouble Subtract(DoubleDouble right) {
return Add(right.Negate());
}
public DoubleDouble Multiply(DoubleDouble right) {
// hi*hi exactly, then the three cross terms. Each cross term is at least 2^-53 smaller
// than the product, so plain double arithmetic 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. Faster than the full double-double 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>Slower and slightly less accurate than the other operations.</summary>
public DoubleDouble Divide(DoubleDouble right) {
// Three steps of long division.
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>Slightly faster 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);
}
public DoubleDouble Negate() {
// Negation never rounds, so this is exact.
return new DoubleDouble(-m_hi, -m_lo);
}
public DoubleDouble Abs() {
return Sign < 0 ? Negate() : this;
}
/// <summary>Slightly less accurate than the other operations. Negative and NaN inputs
/// return NaN, as Math.Sqrt does.</summary>
public DoubleDouble Sqrt() {
if (!(m_hi > 0.0)) { // Not "m_hi <= 0.0", because this form is also true for NaN.
return FromDouble(Math.Sqrt(m_hi));
}
// One Newton step from the double square root.
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);
}
public double ToDouble() {
return m_hi;
}
public int CompareTo(DoubleDouble other) {
// Correct only for normalized pairs.
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;
// Conversion from a double is exact, so it is implicit. Conversion to a double drops the
// low component, so it is explicit.
public static implicit operator DoubleDouble(double value) => FromDouble(value);
public static explicit operator double(DoubleDouble value) => value.ToDouble();
/// <summary>Knuth's two-sum: sum + error == a + b exactly, for any two doubles.</summary>
private static (double Sum, double Error) twoSum(double a, double b) {
// The rounding error of a floating-point addition is itself representable as a double.
// The intermediate results below may round, but their rounding errors cancel out.
double sum = a + b;
double bKept = sum - a; // The part of b that made it into the sum.
double aKept = sum - bKept; // The part of a that made it into the sum.
double aLost = a - aKept;
double bLost = b - bKept;
return (sum, aLost + bLost);
}
/// <summary>Dekker's fast two-sum: the same guarantee as twoSum, valid only when
/// |a| >= |b| (or a == 0).</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 lost nothing, so only b's loss remains.
}
/// <summary>Exact product: product + error == a * b exactly, as long as the product is a
/// normal double. Below about 2e-292 the error no longer fits a double.</summary>
private static (double Product, double Error) twoProduct(double a, double b) {
// The fused multiply-add computes a*b - product with a single rounding, which is the
// low half of the 106-bit product. On a CPU without an FMA instruction,
// Math.FusedMultiplyAdd computes the same result in software, much slower.
double product = a * b;
double error = Math.FusedMultiplyAdd(a, b, -product);
return (product, error);
}
}
The second file contains exact printing and parsing, which are implemented with BigInteger:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
// Printing and parsing of DoubleDouble. Printing 31 correct decimal digits requires
// computing them exactly, which is done with BigInteger (part of the BCL, so still no
// 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, but it can be hundreds of digits long 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";
}
// The denominator is a power of two, so multiplying by the matching power of five
// turns it into a power of ten, which gives the decimal digits 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 contiguous;
/// a hand-built pair such as (1.0, 1e-300) holds more information than 36 digits can
/// represent.
/// </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 the remainder 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, 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 that the quotient is in [2^52, 2^54). One integer division then gives all
// bits of the significand plus the remainder needed for 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. Move it to the exponent and add the dropped bit to the
// remainder so that 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>Parses a decimal string such as "-1.25e-7" 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 formats them.</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);
// Estimate of the decimal exponent of the leading digit from the bit lengths
// (log10(2) ~ 0.30103). The loop below corrects it 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 notation for the range where it is readable.
if (exponent10 >= 0) {
text = text.PadRight(exponent10 + 1, '0');
result.Append(text[..(exponent10 + 1)]);
if (text.Length > exponent10 + 1) {
result.Append('.').Append(text[(exponent10 + 1)..]);
}
} else {
result.Append("0.").Append('0', -exponent10 - 1).Append(text);
}
} else {
result.Append(text[0]);
if (text.Length > 1) {
result.Append('.').Append(text[1..]);
}
result.Append('e').Append(exponent10.ToString(CultureInfo.InvariantCulture));
}
return result.ToString();
}
/// <summary>Integer nearest to numerator/denominator, ties away from zero (both positive).</summary>
private static BigInteger roundToNearest(BigInteger numerator, BigInteger denominator) {
BigInteger quotient = BigInteger.DivRem(numerator, denominator, out BigInteger remainder);
return (remainder << 1) >= denominator ? quotient + 1 : quotient;
}
}
Appendix B: the benchmarked kernels#
The numbers in section 5: Double-double operations benchmark come from these benchmarks. The pixel loop and the threading are the same for all three kernels and are not shown.
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
namespace DoubleDoubleSample.Mandelbrot;
/// <summary>
/// The Mandelbrot escape-time kernel in three versions: double, DoubleDouble and MPFR, plus the
/// scratch values the MPFR version needs. Parsing the center, dividing up rows and timing are in
/// MandelbrotRenderer.cs.
///
/// All three run the same seven operations per iteration and count the same iterations for
/// every pixel; each uses the cheapest escape test for its type. A fourth variant, which
/// allocates a fresh MPFR value per operation the way a typical wrapper does, is also in
/// MandelbrotRenderer.cs.
/// </summary>
public static unsafe partial class MandelbrotRenderer {
/// <summary>Bailout radius squared. Large, so that the smooth-iteration formula is accurate.</summary>
private const double ESCAPE_SQUARED = 65536.0 * 65536.0;
/// <summary>Plain double, for reference.</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 code with <c>DoubleDouble</c> substituted for <c>double</c>. The escape test reads
/// the leading components directly.
/// </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, the way it would be written in C.
/// At these precisions MPFR takes its temporary space from the stack, so the loop does not
/// allocate.
/// </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;
// Only escaped pixels use the magnitude, and for those Temp still holds it because
// the loop broke before 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 gives 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 does not need a conversion.</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>Frees all values. Both kernels leave every value allocated when they return.</summary>
public void Dispose() {
MpfrNative.Free(CRe);
MpfrNative.Free(CIm);
MpfrNative.Free(ZRe);
MpfrNative.Free(ZIm);
MpfrNative.Free(ZReSquared);
MpfrNative.Free(ZImSquared);
MpfrNative.Free(Temp);
MpfrNative.Free(Temp2);
MpfrNative.Free(Escape);
}
}
}
The escape test of the double-double kernel differs slightly, but it's just a small optimization with no precision loss, because the orbit past |z| > 2 is diverging exponentially and the low part has no real effect on the number of iterations.
The fourth kernel (the last row of the kernel table) is the third one with each operation producing a new value instead of writing into the scratch.