8

The following C program produces different results on my Mac and on Linux. I'm suprised because I assumed that the implementation of libm is somehow standardized

#include<math.h>
#include<stdio.h>

int main()
{
  double x18=-6.899495205106946e+01;
  double x19 = exp(-x18);
  printf("x19     = %.15e\n", x19);
  printf("x19 hex = %llx\n", *((unsigned long long *)(&x19)));
}

The output on Mac is

x19     = 9.207186811339878e+29
x19 hex = 46273e0149095886

and on Linux 

x19     = 9.207186811339876e+29
x19 hex = 46273e0149095885

Both were compiled without any optimization flags as follows:

gcc -lm ....

I know that I never should compare floats to be excatly the same.

This issue came up during debugging, regrettably the algorithm using this calculation proofs to be numerically unstable and this slight difference leads to significant deviations in the final result. But this is a different problem.

I'm just surprised that such basic operations as exp are not standardized as I can expect for the basic algebraic operations specified by IEEE 754.

Are there any assumptions about precision I can rely on for different implementations of libm for different machines or for different versions ?

Because of the discussion below I used mpmath to compute the value with higher than machine precision and I get with two more figures the result 9.2071868113398768244, so for both of my results the last figure is already wrong. The result on linux can be explained by down rounding this value, the Mac result is also off if the computer uses up rounding.

rocksportrocker
  • 7,251
  • 2
  • 31
  • 48
  • 1
    Well... the standard operations are listed https://en.wikipedia.org/wiki/IEEE_754-1985#Standard_operations . What are you asking though? – Eugene Sh. Jun 26 '17 at 18:08
  • Sorry, added extra sentence. – rocksportrocker Jun 26 '17 at 18:17
  • 1
    BTW, you can print "exact" floating values in hex by using the `%a` format specifier. – Jens Gustedt Jun 26 '17 at 18:38
  • Ah, thanks ! Never seen `%a` before. – rocksportrocker Jun 26 '17 at 18:45
  • 1
    Also see Jonathan Wakely's [Why `` is more complicated than you might think](https://developers.redhat.com/blog/2016/02/29/why-cstdlib-is-more-complicated-than-you-might-think/) from the Red Hat blogs. Wakely is one of GCC's C++ standard library maintainers. I think `` vs `` is a much more interesting case study. – jww Jun 26 '17 at 18:53
  • 1
    I do think you should specify what compiler, compiler flags and CPU architecture was being used. The most simple math.h functions end up being replaced by intrinsics by the compiler in many settings, mapping directly to CPU instructions. This could mean doing 80-bit x87 operations internally that are then rounded to double, for example. – cnettel Jun 26 '17 at 19:13
  • Mac uses *clang*, Linux (probably) *gcc*. I really don't think this is a question of `libm`, but rather the compiler (most FP ops should end up as intrinsics). – tofro Jun 26 '17 at 19:24
  • https://gist.github.com/uweschmitt/e56f124784fcf4cbe428ccb56ceda234 shows the output of both calls of `gcc -v ....`. I apologize for the plenty of output. – rocksportrocker Jun 26 '17 at 20:17
  • I calculated `exp(68.99495205106946)` with PARI/GP and some online calculators, and always got the result (rounded to 15 decimal places) `9.207186811339866e+29`. Which would mean that both your results are wrong. – Martin R Jun 27 '17 at 09:18
  • @MartinR I get different results using `mpmath` for multi precison math !? – rocksportrocker Jun 28 '17 at 14:13

2 Answers2

2

The C99 Specification states (other version should be similar):

J.3 Implementation-defined behavior

1 A conforming implementation is required to document its choice of behavior in each of the areas listed in this subclause. The following are implementation-defined:

...

J.3.6 Floating point

— The accuracy of the floating-point operations and of the library functions in <math.h> and <complex.h> that return floating-point results (5.2.4.2.2).

Meaning GNU libm and BSD libm are free to have different levels of accuracy. Likely what is happening, is that the BSD implemention on OSX rounds to the nearest (unit in the last place) ULP, and the GNU implementation truncates down to the next ULP.

dlasalle
  • 1,615
  • 3
  • 19
  • 25
  • 1
    And the IEEE-754 standard only requires correctly rounded (error <= 0.5 ULP) for addition, subtraction, multiplication, division, and square root. Other operations are not required to be correctly rounded. The two results are within 1 ULP so the error, in this case, is <= 1 ULP. – casevh Jun 26 '17 at 23:03
  • I extended my question with a more precise result using `mpmath` indicating that the error is above 1 ULP on mac. – rocksportrocker Jun 28 '17 at 14:25
  • @rocksportrocker Please see my answer below. It has more details. – casevh Jun 28 '17 at 17:57
1

IEEE-754 behavior is specified at the binary level. Using a Linux, I get identical values for Python's native math library, mpmath, and MPFR (via gmpy2). However, conversion to decimal varies between the three methods.

>>> import mpmath, gmpy2
>>> import mpmath, gmpy2, math
>>> x18=68.99495205106946
>>> x19=math.exp(x18)
>>> mp18=mpmath.mpf("68.99495205106946")
>>> mp19=mpmath.exp(mp18)
>>> gp18=gmpy2.mpfr("68.99495205106946")
>>> gp19=gmpy2.exp(gp18)
>>> x18 == mp18
True
>>> x18 == gp18
True
>>> x19 == mp19
True
>>> x19 == gp19
True
>>> print(x18, mp18, gp18)
68.99495205106946 68.9949520510695 68.994952051069461
>>> print(x19, mp19, gp19)
9.207186811339876e+29 9.20718681133988e+29 9.2071868113398761e+29

After conversion to Python's arbitrary precision integer form, all three results also show as exact.

>>> hex(int(x19))
'0xb9f00a484ac42800000000000'
>>> hex(int(mp19))
'0xb9f00a484ac42800000000000'
>>> hex(int(gp19))
'0xb9f00a484ac42800000000000'

So (at least one) Linux math library, mpmath, and gmpy2.mpfr agree.

Disclaimer: I maintain gmpy2 and have contributed to mpmath in the past.

casevh
  • 11,093
  • 1
  • 24
  • 35
  • It is worth noting that IEEE-754 is not part of the C standard (or C++), nor is it required by Linux (given the restricted set of hardware OSX can legally run on it's probably safe to say it's implementation will use IEEE-754). – dlasalle Jun 28 '17 at 18:28
  • @dlasalle True. While IEEE-754 is not required, GNU libc tries to comply with the IEEE-754-2008 standard and many|most Linux systems will therefore be compliant. But it is not a safe assumption. If your application relies on reproducible results, including transcendental functions, then you should use a separate library, for example MPFR. – casevh Jun 28 '17 at 19:21