exp(1) != e

[StefanKarpinski]

See comment in 77a8888. Regardless of syntax, it seems problematic that exp(1) does not quite yield Euler's constant to all digits:

julia> exp(1) == e
false

julia> exp(1)
2.7182818284590455

julia> float(e)
2.718281828459045

julia> bits(exp(1))
"0100000000000101101111110000101010001011000101000101011101101010"

julia> bits(float(e))
"0100000000000101101111110000101010001011000101000101011101101001"

[StefanKarpinski]

Since exp(1) just calls openlibm's exp function, this appears to be an openlibm issue. cc: @ViralBShah

[ViralBShah]

This seems to be an artefact of openlibm. I verify that this is not an issue with OS X libm or in octave (which presumably uses the os x libm). It is worth checking if the FreeBSD upstream has fixed this in msun, or else we may have to do it ourselves.

[StefanKarpinski]

It's worth noting that MPFR does not have this problem:

julia> exp(big(1))
2.718281828459045235360287471352662497757247093699959574966967627724076630353555e+00 with 256 bits of precision

julia> big(e)
2.718281828459045235360287471352662497757247093699959574966967627724076630353555e+00 with 256 bits of precision

julia> exp(big(1)) == e
true

I know it's a totally different can of worms, but thought it was worth checking.

[JeffBezanson]

Both are less than 1 ulp from the true value, but e is closer.

[StefanKarpinski]

Even if both values were acceptable, we still need to be consistent. Currently we have x^1 == x for any other x but e^1 != e. Since float(big(e)) == e we really have no choice but to go with the correctly rounded value that the constant has, so openlibm needs to change here.

[ViralBShah]

I was going to also say that it was less than 1 ulp, but figured it is better to match other libm implementations and be as accurate as possible.

[JeffBezanson]

Looking at the code, there is something intriguing about the constants:
https://github.com/JuliaLang/openlibm/blob/master/src/e_exp.c#L95

The coefficients are supposed to be 1/6, -1/360, etc. but instead seem to be "fuzzed" a bit. The code specifically states that these constants are intended and even provides them in hex to be extra clear, but no reason is given.

Just for fun I tried replacing them with the "correct" fractions (determined using rationalize), and I got an answer for exp(1) off by 1 ulp in the other direction, and worse overall. I wonder what is going on there.

[mschauer]

They say it: "we use a special Remez algorithm on [0,0.34658] to generate a polynomial of degree 5 to approximate R" [which is degree infinity].

exp(x = 1) is not a special value here, it first gets range reduced to
r = x-log(2) = 0.30... so that it lies in the [0,0.34658] interval where the approximation is valid. If one adapts P1...P2 to exp(1) = e, then the statement "according to an error analysis, the error is always less than 1 ulp" possibly becomes wrong

[mschauer]

Btw, under which circumstances does this difference exist, I don't see this here

Julia Version 0.3.0-prerelease+261
Commit 97b5983 (2013-11-30 12:55 UTC)
Platform Info:
  System: Linux (i686-linux-gnu)
  WORD_SIZE: 32
  BLAS: libopenblas (DYNAMIC_ARCH NO_AFFINITY)
  LAPACK: libopenblas
  LIBM: libopenlibm

julia> bits(exp(1))
"0100000000000101101111110000101010001011000101000101011101101001"

julia> bits(float(e))
"0100000000000101101111110000101010001011000101000101011101101001"

[JeffBezanson]

Ah, I see, some iterations are done on the coefficients.

glibc does not have this problem. In some cases we might be picking up the function in the libm linked to llvm or elsewhere; it can be hard to avoid.

I believe we are already using the newest version of the msun code. I've tried a few other implementations and have yet to find one that is BSD/MIT licensed and avoids this problem. We might have to reproduce their derivation of the series and add a term.

[StefanKarpinski]

Although a more general fix is certainly desirable, we could always just special-case the argument value 1.0.

On Dec 1, 2013, at 2:24 AM, Jeff Bezanson notifications@github.com wrote:

Looking at the code, there is something intriguing about the constants:
https://github.com/JuliaLang/openlibm/blob/master/src/e_exp.c#L95

The coefficients are supposed to be 1/6, -1/360, etc. but instead seem to be "fuzzed" a bit. The code specifically states that these constants are intended and even provides them in hex to be extra clear, but no reason is given.

Just for fun I tried replacing them with the "correct" fractions (determined using rationalize), and I got an answer for exp(1) off by 1 ulp in the other direction, and worse overall. I wonder what is going on there.

—
Reply to this email directly or view it on GitHub.

[JeffBezanson]

I thought about special-casing maybe integer arguments, but I don't think it's worth it since we are also off for nearby arguments such as 1+eps(1.0).

According to sage the next coefficient in the taylor series is -1.0568380277374986e-9 (for the t^6 term). Now I just need to rerun the Remez optimization.

[JeffBezanson]

I tried a couple programs that compute Remez polynomials, without much success. I was not able to exactly reproduce the coefficients given in the code, but the numbers I got were very close. I generated a couple higher-order polynomials, but none of them were any better. For example

    P1 =  1.66666666666666666e-1,
    P2 = -2.77777777777777784e-3,
    P3 =  6.61375661375056689e-5,
    P4 = -1.65343916232876607e-6,
    P5 =  4.17532165128832177e-8,
    P6 = -1.05908122239406444e-9,
    P7 =  2.45907466209372507e-11;

[JeffBezanson]

For reference, I tried the remez function in chebfun, which gave warnings on the relevant function r*(exp(r)+1)/(exp(r)-1). I couldn't get it to work but then again I don't really know anything about function approximation. I also tried a remez routine in the matlab file exchange, but it hit lots of singular matrices and didn't really work. This thing (http://lolengine.net/wiki/doc/maths/remez) actually seemed to work pretty well and gave me the above result, but still no dice.

[ViralBShah]

Can @alanedelman help us here?