News: Visit the official fractalforums.com Youtube Channel  ## The All New FractalForums is now in Public Beta Testing! Visit FractalForums.org and check it out!

 Pages:    Go Down       Author Topic: [Java] Double-double library for 128-bit precision.  (Read 14979 times) Description: Higher precision than FPU/SSE2, faster than arbitrary precision. 0 Members and 1 Guest are viewing this topic.
Zom-B
Guest « on: November 17, 2009, 01:31:21 PM »

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:

{
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 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. DoubleDouble.txt (44.44 KB - downloaded 909 times.) DoubleDoubleImprove.txt (2.88 KB - downloaded 591 times.) « Last Edit: November 17, 2009, 01:41:00 PM by Zom-B » Logged
Zom-B
Guest « Reply #1 on: November 17, 2009, 02:00:14 PM »

Usage example

Code:
int iter;
// init: z = c
zx = cx.clone();
zy = cy.clone();
for (iter = 0; iter < maxiter; iter++)
{
// Calculate complex z' = z^2 + c in real arithmetic:
// x' = x*x - y*y + cx
// y' = x*y + y*x + cy
// and for bailout |z|:
// x*x + y*y
// Note that x*x and y*y occur twice. Precalculate the
// squares and reverse the order of iteration and bailout.
// An incidental benefit is that there is no conflict
// between old and new z components anymore when y' is
// calculated before x'.

// xx = x * x
xx = zx.sqr();
// yy = y * y
yy = zy.sqr();

// test (xx + yy > bailout)
// Actually only a few digits accuracy needed, so the .hi
// component of the DoubleDouble result suffices.
{
break;
}

// y' = y * x * 2 + cy
zy.mulSelf(zx);
// x' = xx - yy + cx
xx.subSelf(yy);
}

Full example attached, but requires my Digisoft library. Major rewrite necessary to make it standard-Java.
 « Last Edit: November 17, 2009, 02:04:35 PM by Zom-B, Reason: attachment » Logged
Axel
Guest « Reply #2 on: December 16, 2009, 11:40:02 PM »

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,
The MPFR library (multiple precision floating-point reliable library, see http://www.mpfr.org) does exactly that...  Logged
Zom-B
Guest « Reply #3 on: December 24, 2009, 12:59:09 PM »

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,
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. Logged
Botond Kósa
Fractal Lover  Posts: 233   « Reply #4 on: March 31, 2010, 02:45:43 PM »

I used the (somewhat buggy) ''QD'' library in C++ from the site http://crd.lbl.gov/~dhbailey/mpdist/index.html as a reference

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

Check out my Mandelbrot set explorer:
http://web.t-online.hu/kbotond/mandelmachine/
hobold
Fractal Bachius Posts: 573 « Reply #5 on: March 31, 2010, 04:15:53 PM »

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. Logged
Botond Kósa
Fractal Lover  Posts: 233   « Reply #6 on: March 31, 2010, 04:48:41 PM »

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.

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? Logged

Check out my Mandelbrot set explorer:
http://web.t-online.hu/kbotond/mandelmachine/
hobold
Fractal Bachius Posts: 573 « Reply #7 on: December 19, 2010, 01:59:01 PM »

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? Logged
David Makin
Global Moderator
Fractal Senior      Posts: 2286    « Reply #8 on: December 19, 2010, 05:45:35 PM »

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) Logged

The meaning and purpose of life is to give life purpose and meaning.

http://www.fractalgallery.co.uk/
"Makin' Magic Music" on Jango
hobold
Fractal Bachius Posts: 573 « Reply #9 on: December 19, 2010, 07:50:57 PM »

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.  Logged
David Makin
Global Moderator
Fractal Senior      Posts: 2286    « Reply #10 on: December 20, 2010, 04:03:48 AM »

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 Logged

The meaning and purpose of life is to give life purpose and meaning.

http://www.fractalgallery.co.uk/
"Makin' Magic Music" on Jango
 Pages:    Go Down
 Related Topics Subject Started by Replies Views Last post  double VS long double performance in x86 (Core2Duo) Programming aluminumstudios 14 9185 April 15, 2010, 11:45:20 PM by Botond Kósa  Ultra Fractal double precision Programming cabrif 11 2490 January 20, 2012, 11:48:35 PM by David Makin  Is there any way to get double precision camera data? Fragmentarium laser blaster 1 1043 April 21, 2013, 09:56:09 PM by Syntopia  End Of Double Precision Gestaltlupe Gallery trafassel 1 898 November 12, 2013, 04:50:56 PM by eiffie  Playing Double Precision Mandelball Mandelbulber Gallery mclarekin 0 679 July 02, 2015, 01:28:22 PM by mclarekin