Recent Progress of the Java Grande
Numerics Working Group

Subject: Feedback for the Numerics Working Group

Hi Ron,

The following gives my perspective on the Java Grande Numerics Working Group and its relationship and relevance to Sun in Sun's role as the steward of Java development.

As you know, I'm a Distinguished Engineer at Sun, one of the members of the original Java project, the author of The Java Virtual Machine Specification, and currently one of the architects of the Java 2 platform. I work closely with the other architects of the Java technologies, such as James Gosling, and while I can't speak for them in detail I can say that the opinions I express below are not out of line with theirs.

During the development of the Java 2 platform version 1.2 (formerly known as JDK 1.2), I was handed responsiblity for creating licensee and industry consensus around changes to Java's primitive floating-point arithmetic targeted at improving performance on Intel CPUs. This was an extremely difficult task because it required careful attention to the balance between many factors, and was being done in an extremely charged political environment. We who were responsible did not have a broad background in numerics, and had not been successful finding help within Sun. We understood that there was high risk: Java had taken a rather different approach to floating point than many other languages, and the wrong decision in 1.2 could throw away many of the advantages of that approach.

Our best attempts led to a public proposal that we considered a bad compromise and were not happy with, but were resigned to. At this point the Numerics Working Group wrote a counterproposal to Sun's public proposal that gave new technical insight on how to balance the demands of performance on Intel without throwing away the important part of reproducibility. The counterproposal was both very sensitive to the spirit of Java and satisfactory as a solution for the performance problem. When we saw the new proposal we revived efforts to reach a better answer.

We were subsequently aided by email and phone calls with a number of members on the Numerics Working Group (you, Roldan, and Cleve). Joe Darcy, one of the authors of the counterproposal, helped us to understand the proposal generally, and specifically to evaluate the effects of modifications we required of the counterproposal. All of you helped us understand how the Numerics Working Group represented a number of rather different fields with different perspectives on numerics. This helped us gain confidence that the new solution reflected a consensus that would satisfy a broad range of Java users.

We are sure that we ended up with a better answer for 1.2, and arrived at it through more complete consideration of real issues, because of the efforts of the Numerics Working Group. We have internally resolved to consult with the Group on future numeric issues, and to do so earlier in the process. Our attendance at the 3/11 Numerics Working Group meeting at Sun and all the work we accomplished there is evidence of this resolution. We think that the Group continues to show great sensitivity to the needs of the Java platform and the difficulty of introducing change while preserving compatibility and robustness. We look forward to continuing this relationship into the future.

-- Tim

The current JVM 1.2 spec supports two double formats: IEEE format with 53 bit mantissa, 12 bit exponent; and a default format, that has exactly the same mantissa, but may have a larger exponent. Internal values (not visible to the user) can be held in the default format, and converted to the IEEE format when assigned to program variables. If a method is declared strictfp, then the use of large exponents is disallowed. A similar design is used for floats. This design enables Intel to take (partial) advantage of its 80 bit floating point registers, using the mode where only 53 bit mantissas are used in floating point operations. It does not allow Power PC to take advantage of the fused multiply add, and does not allow the use of full sized mantissas in the 80 bit Intel floating point registers.

The JGNWG has discussed several proposals to remediate this problem. To be adopted, a proposal should

1. require limited implementation effort, and negligible implementation effort on platforms that are not affected by current Java floating point restrictions. ,
2. introduce only limited platform dependence (similar to the platform dependence due to the default/strictfp design),
3. not worsen accuracy, and
4. improve performance of codes without required significant, platform-dependent, code rewrite.

The proposal outlined below seems to achieve these four goals. An "extended" floating point mode is added, where larger mantissas are allowed, in the same circumstances where larger exponents are allowed in JVM 1.2: internal values can be stored with more digits of precision than required by IEEE format, but need be converted back to IEEE format when assigned to program variables. This solution would allow (i) the use of fused multiply adds, where the intermediary output from the multiplication is not rounded, (ii) the promotion of floats to doubles, on machines where double arithmetic is faster, and (iii) the use of full 80 bit precision on Intel machines – at least for internal values.

We can have two versions of this design. V1: The use of the "extended" format would be allowed only for methods (and classes) that were explicitly tagged with the keyword extendedfp in the source. The byte file will contain an extendedfp flag. Otherwise, the rules of JVM 1.2 apply. We then have two floating-point modes: strictfp and extendedfp (in addition to the default mode). V2: "extended" format can be used in default mode – no new keyword.

Note that the user can always prohibit the use of extended format by explicitly naming each internal value (replace x=a*b+c with x1=a*b; x=x1+c), by using strictfp, or, if V1 is adopted, just by avoiding the use of the extendedfp attribute. Note, too, that the implementation effort is minimal, for implementations that do not take advantage of the extended format: the Java compilers need generate a new flag, but this flag may be ignored by Java interpreters and JITs.

A downside of this design is that it is unlikely that Sun will add new keywords and new floating-point modes. Thus, either (i) we agree upfront what extendedfp entails (does it allow usage of associativity, as proposed with associativefp?), or (ii) we agree that extendedfp may add new relaxations, over time, or (iii) we agree that none will be added. (i) is unlikely, due to short time.

An important design goal of Java has been machine independence. A single program should run on a wide range of computers, and get the same result everywhere. But, as Java matures and becomes widely used for a broad range of applications, the difficulty of achieving this goal in practice is becoming apparent.

The machine independence goal conflicts with execution speed. Keeping intermediate quantities in extended precision floating point registers increases speed, but may produce different results and lead to different exceptions on different machines.

Tradeoffs between machine independence and speed become particularly important in the design of a library for the elementary math functions, square root, exponential, logarithm, power, and a half dozen trig functions.

These fundamental tradeoff considerations are confounded in at least one currently available Java math library by serious design flaws. The exponential, sine and cosine functions in Microsoft's math library are so inaccurate for large arguments that the library is unsuitable for general-purpose use.

However, very few people are paying any attention to this specification. There are two widely used Java environments for Microsoft operating systems on Intel PCs, one from Microsoft and one from Sun. Neither one of them uses currently FDLIM; both use the "libm" from Visual C++. As a consequence, Java numerical results on the PC may not be the same as those obtained on a Sun SPARC, and, worse yet, may be inaccurate or completely wrong.

FDLIBM deserves to be better known and more widely used. It is a very valuable asset. Authored by K. C. Ng, David Hough and possibly others at Sun, it is essentially a publicly available version of Sun's elementary function library. See http://www.netlib.org/fdlibm.

(At the MathWorks, we use FDLIBM as the underlying library for MATLAB on all platforms. We did this primarily because we got tired of working around the bugs, poor algorithms and inconsistent handling of exceptional cases that we found in the libm's provided by other hardware and operating systems vendors.)

FDLIBM has been translated to Java by John Brophy of Visual Numerics in Houston. This is useful in situations where the native code generated by the C version of FDLIBM is inappropriate. Brophy's source code is publicly available, http://www.vni.com/corner/garage/grande/index.html.

Let F(x) be any of the elementary math functions under consideration. How can we assess the accuracy of a routine that computes F(x)? One possibility is to regard any floating point input argument, x, as a precisely known real number, call upon some oracle to evaluate F(x) exactly, and then round the result to the nearest floating point number. There is only one rounding error in this entire theoretical computation. This is what we mean by "the correctly rounded result" or "getting the right answer".

This ultimate accuracy may not be justified by itself, but the implied consequences are very important. If a particular routine is known to fit this model, then it will always get the same results on different machines. If a particular function F(x) has monotonicity properties, such as x <= y implies exp(x) <= exp(y), or range properties, such as -1 <= cos(x) <= 1, then the computed function will inherit these properties.

Several years ago, a widely used Unix math library had cos(0) > 1. The value was an accurate approximation to the exact value; it differed from 1 by only one bit in the last place. But it violated a fundamental property of the cosine function. This kind of violation is not possible when computations produce correctly rounded results.

It is often argued that the input x is not known exactly, so there is no reason to evaluate F(x) exactly. This same argument could be (and actually has been) used to justify sloppy rounding with basic arithmetic operations like additional and multiplication. Again, the goal is not the accuracy itself, but the desirable mathematical consequences that correct rounding brings.

Whether or not this idealized model is justified in any particular context can be discussed in that context. Whether or not such a model can be achieved with reasonable cost in time and program size is a key question in the current discussion.

Ulps stands for "units in the last place". This is a strong measure of error in which the error in y is expressed relative to the spacing of the floating-point numbers near y. For all double precision numbers between 1 and 2, an ulp is 2^-52. For numbers between 2^k and 2^k+1, an ulp is 2^k-52.

The distance from any real number to the nearest floating-point number is less than or equal to the half the spacing between the nearby numbers. So the correctly rounded results in our idealized model have errors <= 0.5 ulps. Documentation for various elementary function libraries, including FDLIBM, often says that the error is less than one ulp. That's good, but not perfect. Half an ulp is the ultimate goal.

The documentation for FDLIBM itself says:

     *   ------------------------------------------
     *   |  Use the hardware sqrt if you have one |
     *   ------------------------------------------

It may be true that the FDLIM software square root and various hardware square roots all produce the same results. But, if that's not the case, what should Java use? For square root, the choice is clear: the hardware square root should be producing the correctly rounded answer, so use it. But, for other functions, the resolution is not so clear.

Let pi be the transcendental mathematical quantity usually denoted by a Greek letter and let PI be Java.lang.Math.PI, the floating point number closest to pi. What is the correct value for Java.lang.Math.sin(PI)? The answer is not zero, because PI is not equal to pi. Here are two answers, one obtained from the FSIN instruction on an Intel Pentium and the other obtained from FDLIBM:

     1.224606353822377e-16
     1.224646799147353e-16

These two values agree to only five significant figures. If the FDLIBM result is regarded as the correct answer, then the error in the FSIN result is 1.64e+11 ulps, or 164 gigaulps. It can be argued that both values are so small that either one is an acceptable result, but the mere fact that more than one answer is possible violates both the machine independence objective and the idealized model.

The first stage in the computation of sin(x) is argument reduction. If abs(x) is greater than pi/2, then an integer multiple of some approximation to pi is subtracted from x to bring it into the desired range. Different hardware and software manufacturers use different approximations to pi in this argument reduction. The FSIN instruction on Intel's Pentium has a 66-bit approximation to pi available in silicon for this purpose. Before the Pentium, earlier Intel 80x87 chips used less accurate approximations to pi. Documentation for new Cyrix and AMD chips indicate they may use more than 66 bits.

FDLIBM has a 1144-bit approximation to pi. This allows argument reduction equivalent to an infinitely precise pi for any input value less than 2¹⁰²⁴, the largest possible double precision floating point number. (See the papers by K. C. Ng available from http://www.validgh.com/.)

So, for sin(PI), the FDLIBM sin function gives the "correct" result required by our idealized model, while the hardware instructions on various chips give different, less correct, results. On the Pentium, the Java Development Kits from both Microsoft and Sun appear to use the less accurate FSIN instruction.

Using these different values for pi for the argument reduction in sin(), cos() and tan() is the most serious deterrent to machine independent results. When measured if terms of relative error, or ulps, the discrepancies can be huge.

Moreover, it is going to be difficult to get everyone to agree on one particular value for pi. The infinitely precise scheme in FDLIBM has the best overall mathematical properties, but its advantages may not be strong enough to convince PC uses to abandon their builtin functions.

Because the Pentium FSIN instruction only "knows" 66 bits of pi, it can not do accurate reduction of large arguments. When invoked with an argument larger than 2⁶³, the instruction sets an error bit in the floating-point status word and returns the argument unchanged. Systems, like Java, that use the FSIN instruction without checking the error bit, find that for large x, sin(x) is computed really fast, but is equal to x. This is unacceptable for general purpose technical computing.

Like sin(x), computation of exp(x) begins with an argument reduction. The role of pi is played by ln2, the natural logarithm of 2. For any x, let n = round(x/ln2). Then

The reduced argument is in a range where exp() can be approximated quickly and accurately by a modest number of terms in a power series. The final scaling by 2ⁿ is accomplished by simply adding n to the floating-point exponent.

The integer n has been chosen so that severe subtractive cancellation occurs in the computation of the reduced argument, x - n*ln2. If the reduced argument is to have full relative precision, then ln2 must be represented to, and the computation must be carried out in, higher precision. It appears that the exp function in Microsoft's libm, which is used by Visual C and by both Microsoft's and Sun's JDKs, fails to do this. For example,

     java.lang.Math.exp(100) = 2.688117141816161e+043
                                              ^^

whereas the exact answer is

                    exp(100) = 2.688117141816135448... * 10^43

                                              ^^^^^

More extensive experiments show that as x increases, the maximum error in exp(x) doubles whenever x/ln2 cross a power of 2. This is strong evidence that the argument reduction is not being done carefully.

Is a perfectly rounded elementary math library possible and practical? The difficulty is that nobody knows how much computation is required to achieve perfect rounding in all cases. If F(x) is a transcendental function, then, except for special arguments like x = 0 and x = 1, the value of F(x) is irrational. So, it can't be a floating-point number, and it can't fall exactly halfway between two floating-point numbers. But it can fall arbitrarily close to halfway between two numbers. A potentially very large number of extra digits is required to decide which way to round. This is the Table Maker's Dilemma.

The exact value of x is 7543731635572464 x 2^-53, or, in hexadecimal, x = 1.ACCFBE46B4EF0 x 2^-1. The decimal representation of exp(x) begins with

The reason why we are interested in this number is revealed by the hexadecimal expansion of exp(x)

The first blank is at the point where we have to decide how to round. Should we round down to a floating-point number that ends in EE70, or up to EE71? The binary expansion after the rounding point has a zero, followed by 54 ones, then another zero. We need to go into the third double precision number before we can make the correct decision. If we decide to accept

the error is 0.4999999999999999818 ulps. If, instead, we were to round up to EE71 the error would be 0.5000000000000000182 ulps. We have found our correctly rounded result, but it wasn't easy to get it.

Lefevre and Muller have done an exhaustive search to establish that this is the worst case for the exponential function with arguments between 1/2 and 1. But no one knows how much extra precision is needed to do perfect rounding for all elementary functions and all double precision floating point arguments.

Two research groups, one at Sun in Menlo Park and one at IBM in Haifa, are currently working on elementary function libraries that deliver correctly rounded results. As far as we know, their efforts have not yet been reviewed by others. We do not know how the execution time compares with machine instructions or FDLIBM. And, the libraries may use fairly large tables, so their memory requirements may be significant.

Proving that such a library is perfect by exhaustive search is out of the question. A rigorous proof will have to involve careful mathematical analysis of the underlying algorithms and code.

Here are my recommendations for the vendors of Java libraries and the maintainers of the Java specifications.

• A. Ziv, Fast evaluation of elementary mathematical functions with correctly rounded last bit, ACM Transactions on Mathematical Software 17 (1991), pp. 410-423. http://www.acm.org/pubs/citations/journals/toms/1991-17-3/p410-ziv/