It is somewhat surprising that the cexpexp tests fail because of this. Maybe the precision in Float16 is not enough for that function?

Yes, I assume that's likely the problem. Maybe we should add a specialization that performs calculations with Float32 and converts the result to Float16.

Tests with Float32 seem to be affected from the accuracy problems as well. I strongly suspect that the initial culprit is JuliaLang/julia#37440 as tests pass with Julia 1.6 but fail with Julia >= 1.7.0-beta1 (cc @oscardssmith).

I'm not fully sure what the accuracy problem would be. expm1(Float16) and exp(Float16) return the correctly rounded result for all inputs.

The accuracy problem seems to be caused by expm1 (as this one is affected by the linked PR). I assume it is only apparent in the tests of cexpexp since the first exp seems to create inputs in a problematic range. These tests (comparisons with highly accurate BigFloat computations) pass on Julia < 1.7.0-beta1, so to me there seems to be a clear regression.

Iterating over all Float16s, expm1(x) returns identical results to Float16(expm1(big(x))

The tests clearly show that cexpexp tests for Float32 and Float16 pass with old Julia versions but fail with recent ones. Doesn't this clearly indicate that a regression/bug was introduced at some point?

Simple example:

On Julia >= 1.7.0-beta1:

julia> Float32(expm1(-exp(big(Float32(0.1)))))
-0.6688457f0

julia> expm1(-exp(Float32(0.1)))
-0.3376915f0

On Julia < 1.7.0-beta1:

julia> Float32(expm1(-exp(big(Float32(0.1)))))
-0.6688457f0

julia> expm1(-exp(Float32(0.1)))
-0.6688458f0

so...

1. that doesn't use Float16.
2. I can't reproduce that. (at least on 1.9.2)

As mentioned above, Float32 is failing as well, so it's not a Float16-exclusive problem. That's why the problem can't be fixed by computing the Float16 result with Float32 (as tried above) and why I was pointing to the Float32-expm1 PR. The same problem can be observed with Float16 as well though:

On Julia >= 1.7.0-beta1:

julia> Float16(expm1(-exp(big(Float16(0.1)))))
Float16(-0.669)

julia> expm1(-exp(Float16(0.1)))
Float16(-0.338)

On Julia < 1.7.0-beta1:

julia> Float16(expm1(-exp(big(Float16(0.1)))))
Float16(-0.669)

julia> Base.expm1(x::Float16) = Float16(expm1(Float32(x))) # not defined on Julia <= 1.6.7

julia> expm1(-exp(Float16(0.1)))
Float16(-0.669)

I can't reproduce that. (at least on 1.9.2)

I can, both the Float32 and the Float16 issue:

julia> Float16(expm1(-exp(big(Float16(0.1)))))
Float16(-0.669)

julia> expm1(-exp(Float16(0.1)))
Float16(-0.338)

julia> Float32(expm1(-exp(big(Float32(0.1)))))
-0.6688457f0

julia> expm1(-exp(Float32(0.1)))
-0.3376915f0

julia> versioninfo()
Julia Version 1.9.2
Commit e4ee485e909 (2023-07-05 09:39 UTC)
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 8 × 11th Gen Intel(R) Core(TM) i7-1165G7 @ 2.80GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-14.0.6 (ORCJIT, tigerlake)
  Threads: 8 on 8 virtual cores
Environment:
  JULIA_PKG_USE_CLI_GIT = true
  JULIA_NUM_THREADS = auto

I assume that maybe you use a different OS and/or CPU?

Something is very seriously wrong. Luckily I have exactly the same CPU as you on my laptop and I can reproduce this. To make things very weird, this doesn't show up in the debugger

Thank you very much for finding this. This appears to be a very bad bug in LLVM.

never mind. It was a bug in expm1 JuliaLang/julia#50989. I can't believe no one noticed this for 2 years...

I disabled tests of cexpexp with Float32 and Float16 for the affected Julia versions (on Julia 1.0 and Julia nightly tests are run, only on Julia 1.9.2 they are omitted).

Fine with the PR and the disabled tests @tpapp? 🙂