Finding a bug in FDLIBM pow
By Darcy-Oracle on Feb 12, 2010
Writing up a piece of old work for some more
Friday fun, an example of
testing where the failures are likely to be led to my independent discovery of a bug in the FDLIBM
pow function, one of only two bugs fixed in
Even back when this bug was fixed for Java some time ago
the FDLIBM library was well-established, widely used in the Java platform and elsewhere, and already thoroughly tested so I was quite proud my tests found a new problem. The next most recent change to the
pow implementation was eleven years prior to the fix in 5.3.
The specification for
Math.pow is involved, with over two dozen special cases listed. When setting out to write tests for this method, I re-expressed the specification in a tabular form to understand what was going on. After a few iterations reminiscent of tweaking a Karnaugh map, the table below was the result.
|x||–∞||–∞ < y < 1||–1||–1 < y < 0||–0.0||+0.0||0 < y < 1||1||1 < y < +∞||+∞||NaN|
|–∞ < y < –1||+0.0||f3(x, y)||f3(x, y)||+∞|
|–1 < y < 0||+∞||+0.0|
|0 < y < 1||+∞||x||+0.0|
|1 < y < +∞||+0.0||x||+∞|
f1(y) = isOddInt(y) ? –∞ : +∞;
f2(y) = isOddInt(y) ? –0.0 : +0.0;
f3(x, y) = isEvenInt(y) ? |x|y : (isOddInt(y) ? –|x|y : NaN);
† Defined to be +1.0 in C99, see §F.9.4.4 of the C99 specification. Large magnitude finite floating-point numbers are all even integers (since the precision of a typical floating-point format is much less than its exponent range, a large number will be an integer times the base raised to a power). Therefore, by the reasoning of the C99 committee,
pow(-1.0, ∞)was like
pow(-1.0, Unknown large even integer)so the result was defined to be
The range of arguments in each row and column are partitioned into eleven categories, ten categories of finite values together with NaN (Not a Number). Some combination of x and y arguments are covered by multiple clauses of the specification.
A few helper functions are defined to simplify the presentation. As noted in the table, a cross-platform wrinkle is that the C99 specification, which came out after Java was first released, defined certain special cases differently than in FDLIBM and Java's
A regression test based on this tabular representation of
pow special cases is
jdk/test/java/lang/Math/PowTests.java. The test makes sure each interesting combination in the table is probed at least once. Some combinations receive multiple probes.
When an entry represents a range, the exact endpoints of the range are tested; in addition, other interesting interior points are tested too. For example, for the range 1 < x< +∞ the individual points tested are:
+1.0000000000000002, // nextAfter(+1.0, +oo) +1.0000000000000004, +2.0, +Math.E, +3.0, +Math.PI, -(double)Integer.MIN_VALUE - 1.0, -(double)Integer.MIN_VALUE, -(double)Integer.MIN_VALUE + 1.0, double)Integer.MAX_VALUE + 4.0, (double) ((1L<<53)-1L), (double) ((1L<<53)), (double) ((1L<<53)+2L), -(double)Long.MIN_VALUE, Double.MAX_VALUE,
Besides the endpoints, the interesting interior points include points worth checking because of transitions either in the IEEE 754
double format or a 2's complement integer format.
Inputs that used to fail under this testing include a range of severities, from the almost always numerical benign error of returning a wrongly signed zero, to returning a zero when the result should be finite nonzero result, to returning infinity for a finite result, to even returning a wrongly singed infinity!
Selected Failing InputsFailure for StrictMath.pow(double, double): For inputs -0.5 (-0x1.0p-1) and 9.007199254740991E15 (0x1.fffffffffffffp52) expected -0.0 (-0x0.0p0) got 0.0 (0x0.0p0). Failure for StrictMath.pow(double, double): For inputs -0.9999999999999999 (-0x1.fffffffffffffp-1) and 9.007199254740991E15 (0x1.fffffffffffffp52) expected -0.36787944117144233 (-0x1.78b56362cef38p-2) got -0.0 (-0x0.0p0). Failure for StrictMath.pow(double, double): For inputs -1.0000000000000004 (-0x1.0000000000002p0) and 9.007199254740994E15 (0x1.0000000000001p53) expected 54.598150033144236 (0x1.b4c902e273a58p5) got 0.0 (0x0.0p0). Failure for StrictMath.pow(double, double): For inputs -0.9999999999999998 (-0x1.ffffffffffffep-1) and 9.007199254740992E15 (0x1.0p53) expected 0.13533528323661267 (0x1.152aaa3bf81cbp-3) got 0.0 (0x0.0p0). Failure for StrictMath.pow(double, double): For inputs -0.9999999999999998 (-0x1.ffffffffffffep-1) and -9.007199254740991E15 (-0x1.fffffffffffffp52) expected -7.38905609893065 (-0x1.d8e64b8d4ddaep2) got -Infinity (-Infinity). Failure for StrictMath.pow(double, double): For inputs -3.0 (-0x1.8p1) and 9.007199254740991E15 (0x1.fffffffffffffp52) expected -Infinity (-Infinity) got Infinity (Infinity).
The code changes to address the bug were fairly simple; corrections were made to extracting components of the floating-point inputs and sign information was propagated properly.
Even expertly written software can have errors and even long-used software can have unexpected problems. Estimating how often this bug in FDLIBM caused an issue is difficult, while the errors could be egregious, the needed inputs to elicit the problem were arguably unusual (even though perfectly valid mathematically). Thorough testing is key aspect of assuring the quality of numerical software, it is also helpful for end-users to be able to examine the output of their programs to help notice problems.