From 278ef74c6a8ded33ef6af34e4262a2e7a7501a22 Mon Sep 17 00:00:00 2001
From: David Widmann <devmotion@users.noreply.github.com>
Date: Mon, 2 Sep 2024 13:31:29 +0200
Subject: [PATCH] Remove code for unsupported Julia versions (#155)

* Remove code for unsupported Julia versions

* Partially revert change

* Slightly improve `_expm1(::Float16)`

* Another fix

* Stay in log space

* Another typo...

* Add `_log2mexp`

* Just fall back to Float32

* Always use `LogExpFunctions.log2mexp`
---
 Project.toml        |  2 +-
 src/distrs/tdist.jl | 21 ++++++---------------
 test/misc.jl        |  7 +------
 3 files changed, 8 insertions(+), 22 deletions(-)

diff --git a/Project.toml b/Project.toml
index 203d730..8240a2e 100644
--- a/Project.toml
+++ b/Project.toml
@@ -25,7 +25,7 @@ ChainRulesCore = "1"
 HypergeometricFunctions = "0.3.10"
 InverseFunctions = "0.1"
 IrrationalConstants = "0.1, 0.2"
-LogExpFunctions = "0.3.2"
+LogExpFunctions = "0.3.26"
 Reexport = "1"
 Rmath = "0.6.1, 0.7"
 SpecialFunctions = "1.8.4, 2.1.4"
diff --git a/src/distrs/tdist.jl b/src/distrs/tdist.jl
index e2673af..ce2cfd1 100644
--- a/src/distrs/tdist.jl
+++ b/src/distrs/tdist.jl
@@ -48,25 +48,16 @@ end
 tdistinvcdf(ν::Real, p::Real) = tdistinvcdf(map(float, promote(ν, p))...)
 
 tdistinvccdf(ν::Real, p::Real) = -tdistinvcdf(ν, p)
-
-if VERSION < v"1.7.0-DEV.1172"
-    function _expm1(x::Float16)
-        if -0.2 < x < 0.1
-            return @evalpoly(x, Float16(0), Float16(1), Float16(1/2), Float16(1/6), Float16(1/24), Float16(1/120))
-        else
-            return exp(x) - 1
-        end
-    end
-end
-_expm1(x::Number) = expm1(x)
-
 function tdistinvlogcdf(ν::T, logp::T) where T<:Real
     if isinf(ν)
         return norminvlogcdf(logp)
-    elseif logp < -log(2)
-        return -sqrt(fdistinvlogccdf(one(ν), ν, logp + logtwo))
     else
-        return sqrt(fdistinvccdf(one(ν), ν, -2*_expm1(logp)))
+        logq = logp + logtwo
+        if logq < 0
+            return -sqrt(fdistinvlogccdf(one(ν), ν, logq))
+        else
+            return sqrt(fdistinvlogccdf(one(ν), ν, log2mexp(logq)))
+        end
     end
 end
 tdistinvlogcdf(ν::Real, logp::Real) = tdistinvlogcdf(map(float, promote(ν, logp))...)
diff --git a/test/misc.jl b/test/misc.jl
index fb1828f..5faf2a5 100644
--- a/test/misc.jl
+++ b/test/misc.jl
@@ -43,12 +43,7 @@ end
             #  B₁(a, b) = B(a, b)
             a = rand(eltya)
             b = rand(eltyb)
-            # Older SpecialFunctions + Julia versions require slightly larger tolerance
-            if VERSION < v"1.3" && promote_type(eltya, eltyb) === Float64
-                @test logmvbeta(1, a, b) ≈ logbeta(a, b) rtol = 10 * sqrt(eps(Float64))
-            else
-                @test logmvbeta(1, a, b) ≈ logbeta(a, b)
-            end
+            @test logmvbeta(1, a, b) ≈ logbeta(a, b)
         end
     end