GNU Linux and the Extended Precision on x86 Processors

The Goal of the Extended Precision

Before listing the drawbacks of the extended precision, one should note that the (optional) extended precision was introduced in the IEEE-754 standard essentially to allow an accurate double-precision implementation of the elementary functions of the mathematical library (exponential, logarithm, trigonometric and hyperbolic functions, etc), thus for very particular cases.

The file /usr/include/fpu_control.h under Linux/x86 says:

#define _FPU_EXTENDED 0x300     /* libm requires double extended precision.  */

but tests (with mpcheck in particular) show that the math library does not seem to need extended precision: double precision seems to be sufficient. Or does anyone have a counter-example? Update: As shown in glibc bug 706, the pow function is inaccurate in some cases, and even more when the processor is configured in double precision.

The Double Rounding

Even if in a programming language, a type corresponding to the double precision is used (for instance, the type double on most C implementations), a processor computing in extended precision internally will not necessarily give the same results as with a processor computing in double precision. For instance, let us consider the two floating-point double-precision numbers x = 2⁵³ + 2 and y = 1 − 2⁻¹⁶, and their sum in the rounding-to-nearest mode z = x + y. On a system working in double precision, the result is the rounding of 2⁵³ + 3 − 2⁻¹⁶, that is, 2⁵³ + 2 (closer to the exact result than the following machine number 2⁵³ + 4 is). On a system working in extended precision, the result will be first rounded in extended precision, that is, 2⁵³ + 3, then rounded in double precision (when stored to memory, for instance), that is, 2⁵³ + 4 (by using the even rounding rule, since 2⁵³ + 3 is the middle of two consecutive machine numbers). Therefore the obtained final result is different.

Note: this phenomenon, called double rounding, occurs in the rounding-to-nearest mode only, but this is the rounding mode that is used in general.

Portability Problems

First the extended precision leads to portability problems, even if much care is taken at the type level, in particular because the support for extended precision is optional, and in practice, many processors do not support it. Also note that the extended precision is not completely normalized. Even if obtaining exactly the same result on two different machines is not regarded as important and the requirements are linked only to the accuracy of the results, the extended precision is useless, since a portable program must assume the worst case.

If the portability is not important and the software needs to compute in extended precision on some given machine (by using the long double type in C, for instance), the software can still configure the processor in rounding to extended precision (in a non-portable way unfortunately, but anyway, the portability has already been sacrificed).

Concerning the portability, it is thus preferable to have everywhere a processor whose default behavior is to round to the double precision, without any major drawback for the programs needing extended precision.

Accuracy Problems

The goal of the extended precision is to allow to have more accurate basic operations in order to obtain more accurate final results. But most programs have not been written specially for the extended precision, and in general, one gets no guarantee. For instance, many programs written in C use the type double, and even if the intermediate computations are carried out in extended precision, some results may automatically be converted to double precision, and these conversions cannot be controlled.

The extended precision can be needed by some software. But in this case, it is important to be able to prove the error bounds and not to rely on assumptions which could be wrong. In particular, using the extended precision with the type double will not allow the user to obtain better certified error bounds; on the contrary, the bounds may be worse because of the double-rounding problem. The only solution is to write non-portable code (as told above), whether the default precision is the double one or the extended one.

For the other software, the choice of the default precision has little importance, but rounding in double precision is preferable concerning the proofs related to the error bounds.

Problems with Some Algorithms

Some algorithms are based on the correct-rounding property and may become completely wrong if a double rounding can occur. Concerning this point, my paper The Euclidean division implemented with a floating-point division and a floor gives an example particularly useful for programs written in ECMAScript (often referred to as Javascript) or using XPath.

Moreover, a same floating-point expression evaluated twice can yield different results and make a comparison fail. It is particularly true with GCC, which does not conform to some points of the ISO C standard. Several people encountered these problems; see the GCC bug 323 and the duplicates (which are in fact not always real duplicates, since the problem is not always the same, even though the first cause is the extended precision), and my program tst-ieee754.c. Note that GCC bug 323 has been fixed in GCC 4.5, as long as the right compiler options are used, like asking to conform to the C standard (see also the corresponding patch): GCC now supports handling floating-point excess precision arising from use of the x87 floating-point unit in a way that conforms to ISO C99. This is enabled with -fexcess-precision=standard and with standards conformance options such as -std=c99, and may be disabled using -fexcess-precision=fast.

Problems Related to Standards Conformance

Finally, some standards such as Java and XPath (and also XSLT, as it is partly based on XPath) require IEEE double precision with correct rounding, in order to obtain completely reproducible results on different machines or software. Default rounding in extended precision makes their implementation non portable (as explained above); and in practice, most implementations do not change the internal precision of the processor, thus are buggy (see below).

Note: Java also has a single-precision type (float), but the choice between the rounding in double precision and the rounding in extended precision does not solve the rounding problems concerning this type; in both cases, without special code, the implementation would be incorrect on this point. However one may assume that no developers concerned in the numerical behavior of their program use the float type for the calculations, especially as this would not be faster.

XSLT Processors

Test the XPath arithmetic with this XSLT stylesheet, to be applied to any XML file, on itself for instance. The following XSLT processors give an incorrect result (message Not a conforming XPath implementation.) under Linux/x86:

• LibXSLT (xsltproc command). Last tested version: 1.1.7, under Debian/unstable. Bug reports: bug 123297 on the GNOME BTS, bug 206549 on the Debian BTS (closed), bug 247465 on the Debian BTS.

• Sablotron (sabcmd command). Last tested version: 1.0, under Debian/unstable. Bug report: bug 257843 on the Debian BTS.

• Xalan-C++ (xalan command). Last tested version: 1.8.0 under Debian/unstable.

Java Virtual Machines (JVMs)

Test the arithmetic of your JVM with this Java class (Java source). Correct result:

$ java test
z = 9.007199254740994E15
d = 0.0

Incorrect result, when the processor rounds in extended precision:

$ java test
z = 9.007199254740996E15
d = 2.0

The GNU interpreter gij (last tested version: 4.3.3) gives this incorrect result under Linux/x86. Bug 255525 on the Debian BTS. Bug 16122 on GCC Bugzilla.

This bug was fixed in the following JVMs:

• SableVM. Fixed in the 1.1.7 release (it was bug 1 on the SableVM BTS). Fixed in the Debian package 1.1.6-4 (see bug 255524 on the Debian BTS).

• JamVM. Fixed in the 1.2.1 release (see bug 260410 on the Debian BTS).

• Kaffe Virtual Machine. Fixed in the 1.1.6-2 Debian package (see bug 255502 on the Debian BTS).

• CACAO. Fixed in 0.94-2 Debian package (version 0.94-1 was incorrect, but not previous versions, see bug 350729 on the Debian BTS).

Note: Sun's and IBM's JVMs always give a correct result.

For more information about Java numerics: Improving Java for Numerical Computation (in particular, see Section Use of floating-point hardware) on the JavaNumerics web site, and Updates to the Java Language Specification for JDK Release 1.2 Floating Point.

The problem of the single precision remains (is it often used?). Developers of JIT compilers may be interested in the article Optimizing precision overhead for x86 processors by Takeshi Ogasawara, Hideaki Komatsu and Toshio Nakatani, in Software — Practice and Experience, 34(9):875–893, July 2004.

Javascript (ECMAScript)

Test of the Javascript arithmetic of your browser: . The number 0 is the correct result. If you get 2, this probably means that the Javascript implementation of your browser uses extended precision, and this is incorrect (for more information, see Section 8.5 of the ECMAScript Language Specification). If you get nothing, this means that your browser does not support Javascript or Javascript has been disabled.

The following web browsers give the incorrect result 2 under Linux/x86:

• Mozilla (last tested version: 1.8b Gecko/20050114), but see below concerning Firefox. Bug 309797 on the Debian BTS.

• Opera (last tested version: 11.00).

The bug has apparently been fixed in Firefox. Version 1.5 Gecko/20060110 gave an incorrect result (idem until Gecko 1.9.1), but recent versions such as Firefox 3.5.6 and later seem to behave correctly. Bug 264912 on Bugzilla (now closed).

PHP

A critical bug has been found in PHP, but it is related to the wider exponent range (and GCC bug 323), rather than the extended precision (the code in question runs in the double-precision mode). See below.