CODEWITHSUNDEEP
javacmathfloating-pointieee-754
Martijn Courteaux
Martijn Courteaux

Reputation: 68887

Java/C: OpenJDK native tanh() implementation wrong?

I was digging through some of the Java Math functions native C source code. Especially tanh(), as I was curious to see how they implemented that one. However, what I found surprised me:

double tanh(double x) {
    ...
    if (ix < 0x40360000) {          /* |x|<22 */
        if (ix<0x3c800000)          /* |x|<2**-55 */
            return x*(one+x);       /* tanh(small) = small */
    ...
}

As the comment indicates, the taylor series of tanh(x) around 0, starts with:

tanh(x) = x - x^3/3 + ...

Then why does it look like they implemented it as:

tanh(x) = x * (1 + x)
        = x + x^2

Which is clearly not the correct expansion, and even a worse approximation than just using tanh(x) = x (which would be faster), as indicated by this plot:

enter image description here

(The bold line is the one indicated on top. The other gray one is the log(abs(x(1+x) - tanh(x))). The sigmoid is of course the tanh(x) itself.)

So, is this a bug in the implementation, or is this a hack to fix some problem (like numerical issues, which I can't really think of)? Note that I expect the outcome of both approaches to be exactly the same as there are not enough mantisse bits to actually perform the addition 1 + x, for x < 2^(-55).

EDIT: I will include a link to the version of the code at the time of writing, for future reference, as this might get fixed.

Upvotes: 8

Views: 251

Answers (1)

John Bollinger
John Bollinger

Reputation: 181104

Under the conditions in which that code is executed, and supposing that IEEE-754 double-precision floating point representations and arithmetic are in use, 1.0 + x will always evaluate to 1.0, so x * (1.0 + x) will always evaluate to x. The only externally (to the function) observable effect of performing the computation as is done instead of just returning x would be to set the IEEE "inexact" status flag.

Although I know no way to query the FP status flags from Java, other native code could conceivably query them. More likely than not, however, the practical reason for the implementation is given by by these remarks in the Javadocs for java.StrictMath:

To help ensure portability of Java programs, the definitions of some of the numeric functions in this package require that they produce the same results as certain published algorithms. These algorithms are available from the well-known network library netlib as the package "Freely Distributable Math Library," fdlibm. These algorithms, which are written in the C programming language, are then to be understood as executed with all floating-point operations following the rules of Java floating-point arithmetic.

The Java math library is defined with respect to fdlibm version 5.3. Where fdlibm provides more than one definition for a function (such as acos), use the "IEEE 754 core function" version (residing in a file whose name begins with the letter e). The methods which require fdlibm semantics are sin, cos, tan, asin, acos, atan, exp, log, log10, cbrt, atan2, pow, sinh, cosh, tanh, hypot, expm1, and log1p.

(Emphasis added.) You will note in the C source code an #include "fdlibm.h" that seems to tie it to the Javadoc comments.

Upvotes: 6

Related Questions