support fp128 and maybe quad-double

[cgbaker]

Extended precision types are occasionally useful in prototyping codes and developing numerical intuition of an algorithm. I am thinking here of double-double/quad-double from the QD package [1]. These are statically sized with reasonable performance. The downside is supporting them via LAPACK/BLAS and other libraries; this would definitely be a value-added component of the Julia environment.

[1] http://crd-legacy.lbl.gov/~dhbailey/mpdist/

[JeffBezanson]

I agree this would be good to have, especially since most of our library routines "just work" when given new data types.
Do you think it's worth exposing 80-bit floats on x86, or does software double-double/quad-double provide a better tradeoff? I'm also a little concerned that those packages seem to require disabling extended precision on x86.

[cgbaker]

80-bit is better than nothing, but I'm a big fan of dd/qd. There's something to be said about 60 decimal digits of accuracy when you're doing a convergence study of some iterative method. The disabling of extended precision on x86 is a concern, and a hassle. You could theoretical turn this on and off before operations on arrays of dd/qd objects.

[JeffBezanson]

Hmm, that is pretty disappointing for composability, and writing generic code. It might be easiest to disable it globally once those libraries are included. I wonder if it is possible to develop x86-friendly versions of these libraries; I know I probably couldn't do it though.

[JeffreySarnoff]

Extended precision that 'just works' at the same precision on different systems offers more valuable value.

When run on x86, double-double and quad-double disable extended-precision before each call to a sensitive function and re-enable it after each such call. Using dd/qd does not involve additional management of that state. Using Float80 is OK; considering it an extended precision solution for x86s -- less so.

[JeffBezanson]

Ok, that is useful feedback. True that Float80 offers relatively little in exchange for using 2x the memory. Sounds like we can just write up an interface to this library any time then, taking care to reproduce the license notice for commercial users. I like that their BSD license is provided as a Word document. Somebody clearly has a sense of humor.

[StefanKarpinski]

The big problem I can see here after downloading the library and playing around with it for a while is that it's written in C++ and doesn't seem to expose basic things like addition via a C-callable interface. It has the ability to be called from Fortran, however, so that gives hope, although I can't figure out where signatures for e.g. addition are in the header files.

I'm also slightly worried about printing these things. We have a "honest printing" policy, in that we use the double-conversion library to print the minimal number of digits necessary to exactly reconstruct each floating-point value. I wonder what guarantees the printing of dd and qd values makes? I'm not super keen on re-implementing the grisu algorithm in Julia, although that may eventually be necessary.

[cgbaker]

I think you are looking for these:
http://qd.sourcearchive.com/documentation/2.3.4/c__dd_8h-source.html
http://qd.sourcearchive.com/documentation/2.3.4/c__qd_8h-source.html

And note c_dd_swrite() in c_dd.cc:

void c_dd_swrite(const double *a, int precision, char *s, int len) {
  dd_real(a).write(s, len, precision);
}

[StefanKarpinski]

@cgbaker: I saw those, but they seem to work on doubles, rather than dd_real and qd_real structs. I see that you can print these things within the libraries, but the grisu algorithm is quite new (2010), and is, as far as I know, the only printing algorithm that guarantees accurate minimal printing of floats as decimals and doesn't have to use slow arbitrary precision numbers to do it.

[cgbaker]

@StefanKarpinski A dd_real in the C interface is a double[2], a qd_real is a double[4]. The double* in those methods is, in fact, qd objects. See here, for example:
http://qd.sourcearchive.com/documentation/2.3.4/c__test_8c-source.html

(learning this as I go ;)

[JeffBezanson]

We don't need minimal printing (at least at first) for dd and qd to be useful. In fact in many cases computing with dd/qd in memory and ultimately producing a double for output probably makes sense. Certainly useful for sending minimal data to be plotted in javascript, for example.

[StefanKarpinski]

@cgbaker: Ah, I see! I did not at all figure that out from looking at the code :-)

@JeffBezanson: That's absolutely true.

So, all-in-all, this looks very usable from Julia.

[nolta]

Another option is libquadmath, which is part of gcc >= 4.6. It includes math and string functions.

[StefanKarpinski]

@cgbaker: does libquadmath also provide all the functionality you'd need? Is there some reason for preferring one library over the other? All other things being equal, I'd be inclined to go with the one that's part of gcc. Even better, of course, would be if there was dd and qd support in LLVM...

[cgbaker]

I wasn't aware of libquadmath. It looks interesting. It seems only to provide 128-bit, comparable to double-double, so QD would seem to have some advantage their. However, I don't have a good sense for how widely used either of these extended precision types would be. I only use them for exploration, never yet in any production problems.

[StefanKarpinski]

Doubling the number of bits already seems like a lot of extra precision. Are there situations where quad doubles allow detecting problems that weren't revealed by double doubles? Although, if we're going to support this, we may as well support both. I wonder if libquadmath is faster or just has a more GPL-ish license or what the deal is.