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
