Improve support for Float16 (#74)

* Improve support for Float16

* Base `cexpexp(::Float16)` on `Float32` computations

* Test only non-broken Julia versions

* Update Project.toml

Project.toml

@@ -1,7 +1,7 @@
 name = "LogExpFunctions"
 uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
 authors = ["StatsFun.jl contributors, Tamas K. Papp <[email protected]>"]
-version = "0.3.25"
+version = "0.3.26"
 
 [deps]
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"

src/LogExpFunctions.jl

@@ -10,6 +10,16 @@ export xlogx, xlogy, xlog1py, xexpx, xexpy, logistic, logit, log1psq, log1pexp,
     softplus, invsoftplus, log1pmx, logmxp1, logaddexp, logsubexp, logsumexp, logsumexp!, softmax,
     softmax!, logcosh, logabssinh, cloglog, cexpexp
 
+# expm1(::Float16) is not defined in older Julia versions,
+# hence for better Float16 support we use an internal function instead
+# https://github.com/JuliaLang/julia/pull/40867
+if VERSION < v"1.7.0-DEV.1172"
+    _expm1(x) = expm1(x)
+    _expm1(x::Float16) = Float16(expm1(Float32(x)))
+else
+    const _expm1 = expm1
+end
+
 include("basicfuns.jl")
 include("logsumexp.jl")
 

src/basicfuns.jl

@@ -237,7 +237,7 @@ function log1mexp(x::Real)
     if x < oftype(float(x), IrrationalConstants.loghalf)
         return log1p(-exp(x))
     else
-        return log(-expm1(x))
+        return log(-_expm1(x))
     end
 end
 
@@ -246,15 +246,15 @@ $(SIGNATURES)
 
 Return `log(2 - exp(x))` evaluated as `log1p(-expm1(x))`
 """
-log2mexp(x::Real) = log1p(-expm1(x))
+log2mexp(x::Real) = log1p(-_expm1(x))
 
 """
 $(SIGNATURES)
 
 Return `log(exp(x) - 1)` or the “invsoftplus” function.  It is the inverse of
 [`log1pexp`](@ref) (aka “softplus”).
 """
-logexpm1(x::Real) = x <= 18.0 ? log(expm1(x)) : x <= 33.3 ? x - exp(-x) : oftype(exp(-x), x)
+logexpm1(x::Real) = x <= 18.0 ? log(_expm1(x)) : x <= 33.3 ? x - exp(-x) : oftype(exp(-x), x)
 logexpm1(x::Float32) = x <= 9f0 ? log(expm1(x)) : x <= 16f0 ? x - exp(-x) : oftype(exp(-x), x)
 
 const softplus = log1pexp
@@ -440,4 +440,4 @@ $(SIGNATURES)
 
 Compute the complementary double exponential, `1 - exp(-exp(x))`.
 """
-cexpexp(x) = -expm1(-exp(x))
+cexpexp(x) = -_expm1(-exp(x))

test/basicfuns.jl

@@ -159,23 +159,28 @@ end
 end
 
 @testset "log1mexp" begin
-    @test log1mexp(-1.0)  ≈ log1p(- exp(-1.0))
-    @test log1mexp(-10.0) ≈ log1p(- exp(-10.0))
+    for T in (Float64, Float32, Float16)
+        @test @inferred(log1mexp(-T(1))) isa T
+        @test log1mexp(-T(1))  ≈ log1p(- exp(-T(1)))
+        @test log1mexp(-T(10)) ≈ log1p(- exp(-T(10)))
+    end
 end
 
 @testset "log2mexp" begin
-    @test log2mexp(0.0)  ≈ 0.0
-    @test log2mexp(-1.0) ≈ log(2.0 - exp(-1.0))
+    for T in (Float64, Float32, Float16)
+        @test @inferred(log2mexp(T(0))) isa T
+        @test iszero(log2mexp(T(0)))
+        @test log2mexp(-T(1)) ≈ log(2 - exp(-T(1)))
+    end
 end
 
 @testset "logexpm1" begin
-    @test logexpm1(2.0)            ≈  log(exp(2.0) - 1.0)
-    @test logexpm1(log1pexp(2.0))  ≈  2.0
-    @test logexpm1(log1pexp(-2.0)) ≈ -2.0
-
-    @test logexpm1(2f0)            ≈  log(exp(2f0) - 1f0)
-    @test logexpm1(log1pexp(2f0))  ≈  2f0
-    @test logexpm1(log1pexp(-2f0)) ≈ -2f0
+    for T in (Float64, Float32, Float16)
+        @test @inferred(logexpm1(T(2))) isa T
+        @test logexpm1(T(2))            ≈  log(exp(T(2)) - 1)
+        @test logexpm1(log1pexp(T(2)))  ≈  T(2)
+        @test logexpm1(log1pexp(-T(2))) ≈ -T(2)
+    end
 end
 
 @testset "log1pmx" begin
@@ -433,9 +438,17 @@ end
     cloglog_big(x::T) where {T} = T(log(-log(1 - BigFloat(x))))
     cexpexp_big(x::T) where {T} = 1 - exp(-exp(BigFloat(x)))
 
-    for x in 0.1:0.1:0.9
-        @test cloglog(x) ≈ cloglog_big(x)
-        @test cexpexp(x) ≈ cexpexp_big(x)
+    for T in (Float64, Float32, Float16)
+        @test @inferred(cloglog(T(1//2))) isa T
+        @test @inferred(cexpexp(T(0))) isa T
+        for x in 0.1:0.1:0.9
+            @test cloglog(T(x)) ≈ cloglog_big(T(x))
+            # Julia bug for Float32 and Float16 initially introduced in https://github.com/JuliaLang/julia/pull/37440
+            # and fixed in https://github.com/JuliaLang/julia/pull/50989
+            if T === Float64 || VERSION < v"1.7.0-DEV.887" || VERSION >= v"1.11.0-DEV.310"
+                @test cexpexp(T(x)) ≈ cexpexp_big(T(x))
+            end
+        end
     end
     for _ in 1:10
         randf = rand(Float64)