Welcome to Fractal Forums

Double-double-precision is a number format that uses two double-precision numbers to represent a single number. (Technically, this technique could use any multiple of floating point numbers, even of different kind)

The first double of the number is called the high part, the other the low part. The high part encodes the main portion of the number, eg. 3.141592653589793, and the low part encodes the exact residue of (π - hi), which is

   π = 3.14159265358979323846264338327950288...
 - hi = 3.14159265358979311599796346854419
 -----------------------------------------
   lo = 0.00000000000000012246467991473531

As two double precision numbers are used together, the combination does not have the same precision as an IEEE 128-bit floating-point number, even though it occupies the same space. This is because they store two signs, two exponents and two mantissas.

Number of bits:

Type	Sign	Exponent	significand	total
Half precision (IEEE 754r)	1	5	10	16
Single	1	8	23	32
Double	1	11	52	64
Double-double	2	22	104	128
Quad	1	15	112	128

Even though the double-double has less precision than the quad-precision, it is able to store numbers Quad-precision can't. Consider this number:
1.00000000000000000000000000000000001
which in Quad-precision would just be 1.0, but can be expressed in double-double precision by the tuple (1.0, 1.0^-45})

---

Arithmetic with Double-double precision is not straight-forward, but is it nonetheless a lot faster than pure software emulation (eg. arbitrary-precision).

Consider these two basic functions in Java:

public DoubleDouble add(DoubleDouble y)
{
   double a, b, c, d, e, f;
   e = this.hi + y.hi;
   d = this.hi - e;
   a = this.lo + y.lo;
   f = this.lo - a;
   d = ((this.hi - (d + e)) + (d + y.hi)) + a;
   b = e + d;
   c = ((this.lo - (f + a)) + (f + y.lo)) + (d + (e - b));
   a = b + c;
   return new DoubleDouble(a, c + (b - a));
}

public DoubleDouble mul(DoubleDouble y)
{
   double a, b, c, d, e;
   a = 0x08000001 * this.hi;
   a += this.hi - a;
   b = this.hi - a;
   c = 0x08000001 * y.hi;
   c += y.hi - c;
   d = y.hi - c;
   e = this.hi * y.hi;
   c = (((a * c - e) + (a * d + b * c)) + b * d) + (this.lo * y.hi + this.hi * y.lo);
   a = e + c;
   return new DoubleDouble(a, c + (e - a));
}

In particular, after some (sub-)operations (bold lines) and at the very end, the numbers must be 'normalized'. This is the process of trying to store the complete result (hi+lo) in the hi part, and then subtracting the hi part from the original number, storing the remainder in the lo part. This introduces a lot of odd structures in the code using parenthesis. When these structures are simplified according to mathematical rules, the functions will actually break. Beware of this when choosing optimization flags.

I used the (somewhat buggy) ''QD'' library in C++ from the site http://crd.lbl.gov/~dhbailey/mpdist/index.html (http://crd.lbl.gov/~dhbailey/mpdist/index.html) as a reference, learned the underlying technique, and fixed many bugs and improved efficiency.

ps. I think Windows Calculator (WinXP and newer) uses Double-double precision. The numbers I get for every operation I can think off, have the same amount of digits as my Double-double implementation generates, which is a rather odd coincidence if it uses another kind of emulation.

---

I've attached my up-to-date Java library, but for those who think about giving the original C version a shot; It's not recommended. Just take a look at the included the improvement notes I recorded as I was porting it.

Usage example

(Curse Java for not supporting operator overloading)


Full example attached, but requires my Digisoft library. Major rewrite necessary to make it standard-Java.

The MPFR library (multiple precision floating-point reliable library, see http://www.mpfr.org) does exactly that... ;)

I see nowhere on that site that it uses floats as internal representation, and furthermore, because their precision is user-specifyable at run time and x-double-precision is not (by nature; doubling precision quadruples the algorithm complexity) this is probably just an arbitrary-precision instance.

Have you tried to port the quad-double arithmetic too? I also implemented double-double arithmetic in my own Java fractal generator, based on the referenced "QD" library, but quad-double arithmetic seems far more complex.

I have not taken the time to fully understand the technique. But should it not be possible to arrive at quad precision by a recursive application of the method? That is, quad precision numbers are simply pairs of doubled doubles, carrying out the above formulas in double double arithmetic.

MPFR is indeed a different kind of beast (but for some hardware platforms with very slow integer multiply (SPARC) it did use floating point internally).

Neither of those methods is the clear winner. It depends on specific strengths and weaknesses of a particular piece of hardware.

I tried to implement quad-double precision using double-double numbers as hi and lo parts, but I got stuck at implementing the multiplication. There is an integer constant in the following line:

a = 0x08000001 * this.hi;

Does anyone know where this mysterious 0x08000001 comes from?

I apologize for resurrecting this old topic. On the one hand, I found out that a "quad-double precision" technique exists. It is not simply a recursive application of the double-double formulas, but based on the same principles. A paper with more details is here: http://crd.lbl.gov/~dhbailey/dhbpapers/arith15.pdf

On the other hand, I have an obscure question about the theory behind this method. I want to apply the technique to pairs of single precision floating point variables. Most of the program does not require more than single precision, but there are a few isolated spots where higher accuracy would be quite useful. I want to port that program to GPUs later on, so I don't want to use double precision variables (too few GPU models have double precision hardware).

The context is OpenCL, where addition, multiplication and fused multiply add are defined to be exact as mandated in the IEEE 754 floating point standard. But division in OpenCL is only accurate to 2.5 ulp. I am unsure if the "doubled" division algorithm can be made to work when the underlying divide operation does not adhere to IEEE 754 rules.

The doubled divide seems to be somewhat self-correcting, but I am just not sure. Does someone perhaps have a pointer to a mathematical analysis of the technique?

Since this topic has been revived I'd like to ask a couple of questions.
I have never actually tried extending precision using either my own code or any of the available libraries but have noticed that probably the fastest way is extending the standard fpu formats (as in double+double) rather than full arbitrary precision.

If extending (fpu) double to double+double then I take it that it must be ensured that the fpu values are always plain double and not in extended precision mode ?
Also can anyone tell me the derivation of the division algorithm for double+double ? (I did work out one myself but would be interested in knowing the "correct" method, mine seems a little long-winded)

Given the relative speeds of floating point hardware and integer ALUs, the double-double method should win for "small" extended precisions (two to four words of machine precision, I guess). On GPUs, there is really no alternative. But beware if the GPU adheres to IEEE 754 floating point rules.

Exact IEEE-754 rounding seems to be required to make the double-double trick work. The paper I mentioned above has some references, but I found them either hard to locate on the web, or hard to digest. :)

OK.

BTW my (hypothetical) method for division is based on:

        (a+b) / (c+d)
 ((a+b)*(c-d)) / ((c+d)*(c-d))
  ((a+b)*(c-d)) / (c^2-d^2)

Assuming a and c in a+b and c+d are normalised then if d is (at least) 2^n smaller than c then d^2 is 2^2n smaller so drops out of the calculations (assuming c uses n bits of precision for the positive mantissa).

     ((a+b)*(c-d))/c^2
 (a+b)/c - ((a+b)*d)/c^2