Finding a bug in FDLIBM pow
By darcy 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
FDLIBM 5.3.
Even back when this bug was fixed for Java some time ago
(5033578),
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} | y | ||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|
x | –∞ | –∞ < y < 1 | –1 | –1 < y < 0 | –0.0 | +0.0 | 0 < y < 1 | 1 | 1 < y < +∞ | +∞ | NaN |
–∞ | +0.0 | f2(y) | 1.0 | f1(y) | +∞ | NaN | |||||
–∞ < y < –1 | +0.0 | f3(x, y) | f3(x, y) | +∞ | |||||||
–1 | NaN^{†} | NaN^{†} | |||||||||
–1 < y < 0 | +∞ | +0.0 | |||||||||
–0.0 | +∞ | f1(y) | f2(y) | +0.0 | |||||||
+0.0 | +∞ | +0.0 | |||||||||
0 < y < 1 | +∞ | x | +0.0 | ||||||||
1 | NaN^{†} | 1.0 | NaN^{†} | ||||||||
1 < y < +∞ | +0.0 | x | +∞ | ||||||||
+∞ | +0.0 | +∞ | |||||||||
NaN | NaN | NaN |
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 likepow(-1.0, Unknown large even integer)
so the result was defined to be1.0
instead ofNaN
.
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 Math.pow
.
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 Inputs
Failure 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.
You very definitely get to keep your nerd card.
You can have mine too and adjust the photo should you ever need to escape to somewhere without a von Neumann machine for a break;
Rgds
Damon
Posted by Damon Hart-Davis on February 12, 2010 at 05:46 AM PST #
Is there a typo in the table at column 1, row 4? The expression:
-Inf < y < 1
should be:
-Inf < y < -1
The symmetry of the table caught my eye. The small part that is asymmetric
Posted by Anonymous on February 12, 2010 at 06:19 PM PST #
sorry for the hanging sentence, hit the post button accidentally...
The symmetry is sort of a divided into two bi-diagonal matrices with one asymmetric exception for an x value of 1 given positive y values. Anyway, the table is a nice visualisation aid.
Posted by Anonymous on February 12, 2010 at 06:24 PM PST #
@Anonymous,
Thanks for the typo correction; I've updated the table accordingly.
Posted by Joe Darcy on February 13, 2010 at 03:00 AM PST #