JDK 9 RFR of JDK-8136874: Bug in port of fdlibm pow to Java

[joe darcy]

Hello,

After a good amount of staring at the code, I was able to figure out the 
cause of the previously reported problem in the initial port of FDLIBM 
pow to Java:

http://mail.openjdk.java.net/pipermail/core-libs-dev/2015-September/035276.html

In brief, I had replaced some magnitude comparisons based directly on 
the high-order 32-bits of a floating-point value with comparisons of the 
64-bit floating-point value instead. In some cases, the trailing 32-bit 
need to be accounted for as well.

For the problematic cases in this bug:

Original fdlibm code where iy and ix and the high-order 32-bits of the 
absolute values of y and x, respectively:

  201     /* |y| is huge */
  202         if(iy>0x41e00000) { /* if |y| > 2**31 */
[snip]       // Not relevant here
  207         /* over/underflow if x is not close to one */
  208             if(ix<0x3fefffff) return (hy<0)? s*huge*huge:s*tiny*tiny;
  209             if(ix>0x3ff00000) return (hy>0)? s*huge*huge:s*tiny*tiny;

Code in corrected port:

  354             // |y| is huge
  355             if (y_abs > 0x1.00000_ffff_ffffp31) { // if |y| > ~2**31
[snip]          // Not relevant here
  360                 // Over/underflow if x is not close to one
  361                 if (x_abs < 0x1.fffff_0000_0000p-1) // |x| < 
~0.9999995231628418
  362                     return (y < 0.0) ? s * INFINITY : s * 0.0;
  363                 if (x_abs > 0x1.00000_ffff_ffffp0) // |x| > ~1.0
  364                     return (y > 0.0) ? s * INFINITY : s * 0.0;

In the end, only two lines of code were meaningfully changed, the ones 
with a greater-than comparison:

              // |y| is huge
-            if (y_abs > 0x1.0p31) { // if |y| > 2**31
+            if (y_abs > 0x1.00000_ffff_ffffp31) { // if |y| > ~2**31
                  final double INV_LN2   =  0x1.7154_7652_b82fep0; //  
1.44269504088896338700e+00 = 1/ln2
                  final double INV_LN2_H =  0x1.715476p0; //  
1.44269502162933349609e+00 = 24 bits of 1/ln2
                  final double INV_LN2_L =  0x1.4ae0_bf85_ddf44p-26; //  
1.92596299112661746887e-08 = 1/ln2 tail

                  // Over/underflow if x is not close to one
-                if (x_abs < 0x1.fffffp-1) // |x| < 0.9999995231628418
+                if (x_abs < 0x1.fffff_0000_0000p-1) // |x| < 
~0.9999995231628418
                      return (y < 0.0) ? s * INFINITY : s * 0.0;
-                if (x_abs > 1.0)         // |x| > 1.0
+                if (x_abs > 0x1.00000_ffff_ffffp0)         // |x| > ~1.0
                      return (y > 0.0) ? s * INFINITY : s * 0.0;
                  /*
                   * now |1-x| is tiny <= 2**-20, sufficient to compute
                   * log(x) by x - x^2/2 + x^3/3 - x^4/4
                   */

In the hex floating-point notation, the first five digits of the 
significant align with the high-order 32-bits. For example, the 
floating-point values with the high-order bits corresponding to 
0x1.00000p31 range over

     0x1.00000_0000_0000p31
     0x1.00000_0000_0001p31
      ...
     0x1.00000_ffff_fffep31
     0x1.00000_ffff_ffffp31

Therefore, when the original sources used

     ix>0x3ff00000

this is equivalent to

     x_abs > 0x1.00000_ffff_ffffp31

since that is the largest value whose high-order 32 bits are 0x3ff00000.

Full webrev of the changes up at

     JDK-8136874: Bug in port of fdlibm pow to Java - Java Bug System
     http://cr.openjdk.java.net/~darcy/8136874.0/

The bulk of the changeset is a new battery of tests including the 
failing values reported earlier.

I checked the other replacements of ix in the pow code and they seem 
okay, but I wouldn't mind another pair of eyes looking over it; webrev 
of the initial port is at http://cr.openjdk.java.net/~darcy/8134795.6/.

I'll review the port of hypot for analogous problems.

Thanks again to Jeff Hain for the reporting the incorrect behavior of 
the port.

Cheers,

-Joe