Improve support for Float16

[devmotion]

I noticed that I had a local branch with some improvements for Float16 (IIRC these issues popped up in StatsFuns).

Ref: JuliaLang/julia#42223

[tpapp]

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

[devmotion]

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.

[devmotion]

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).

[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.

[devmotion]

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.

[oscardssmith]

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

[devmotion]

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?

[devmotion]

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

[oscardssmith]

so...

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

[devmotion]

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?

[oscardssmith]

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

[oscardssmith]

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

[oscardssmith]

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

[devmotion]

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).

[devmotion]

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

[tpapp]

lgtm