Search code examples
javacfloating-pointexp

Java Implementation of Math.exp()

Math.exp() delegate the computation to StrictMath, which in turns delegate it to the fdlibm library in C. At least this is how it looks from what is mentioned in https://docs.oracle.com/javase/8/docs/api/java/lang/StrictMath.html.

However, this C implementation (https://www.netlib.org/fdlibm/e_exp.c) appears to give different results compared to what Java returns: I'm evaluating exp(-0.271153846153846134203746487401076592504978179931640625).

The C implementation gives: 0.7624991798148184063421695100259967148303985595703125-> 0x1.86664ae1101cap-1

while Java gives: 0.76249917981481829531986704751034267246723175048828125-> 0x1.86664ae1101c9p-1 Java code is just:

import java.math.BigDecimal;
double expon2 = Math.exp(-0.271153846153846134203746487401076592504978179931640625);
System.out.println(new BigDecimal(expon2).toPlainString());

C code is https://www.netlib.org/cgi-bin/netlibfiles.pl?filename=/fdlibm/e_exp.c.

How is this possible? I would expect a perfect match, instead.

Solution

  • You're not using the fdlibm library properly. __ieee754_exp is what they call an "elementary function", in other words it's an internal function that is used by the real exp function. Also, the library is multi-platform and does a lot of bit-twiddling and is very sensitive to CPU endianness, so it needs to be configured with the appropriate macros.

    To use the library properly, you need to build it according to the instructions in the readme, and then include the built library in your application. Do not call __ieee754_exp directly.

    Update

    The question has been updated after it had been pointed out that we were using the wrong endianness for fdlibm. The question now is: why Java gives an answer differing from C by just one ulp.

    The answer: the C fdlibm library and Java's StrictMath class actually give the same answer (which is not surprising, since StrictMath uses the same implementation as fdlibm):

    jshell> Double.toHexString(StrictMath.exp(-0.271153846153846134203746487401076592504978179931640625))
$6 ==> "0x1.86664ae1101cap-1"

    On the other hand, Java's Math class uses an intrinsic implementation and gives the same answer as the gcc library in C:

    jshell> Double.toHexString(Math.exp(-0.271153846153846134203746487401076592504978179931640625))
$7 ==> "0x1.86664ae1101c9p-1"

    In C:

    printf("%a", exp(-0.271153846153846134203746487401076592504978179931640625));
Output:
0x1.86664ae1101c9p-1

    The intrinsic implementation is typically the one that is microcoded into the processor, and gives a very slightly different approximation by the very last bit differing by 1.