From 8b896cf3e0a604ad73530569e1b3f008966bb5f8 Mon Sep 17 00:00:00 2001
From: Simone Carlo Surace <51025924+simsurace@users.noreply.github.com>
Date: Wed, 27 Apr 2022 17:56:36 +0200
Subject: [PATCH 1/4] Add fallbacks for `log1pmx` and `logmxp1`

This looks like it would make sense, given that the other functions dispatch on `Real` argument, and it also fixes #44.
---
 src/basicfuns.jl | 8 ++++++++
 1 file changed, 8 insertions(+)

diff --git a/src/basicfuns.jl b/src/basicfuns.jl
index e30809c2..cc247319 100644
--- a/src/basicfuns.jl
+++ b/src/basicfuns.jl
@@ -246,6 +246,7 @@ $(SIGNATURES)
 Return `log(1 + x) - x`.
 
 Use naive calculation or range reduction outside kernel range.  Accurate ~2ulps for all `x`.
+This will fall back to the naive calculation for argument types different from `Float64`.
 """
 function log1pmx(x::Float64)
     if !(-0.7 < x < 0.9)
@@ -267,10 +268,14 @@ function log1pmx(x::Float64)
     end
 end
 
+# Naive fallback
+log1pmx(x::Real) = log1p(x) - x
+
 """
 $(SIGNATURES)
 
 Return `log(x) - x + 1` carefully evaluated.
+This will fall back to the naive calculation for argument types different from `Float64`.
 """
 function logmxp1(x::Float64)
     if x <= 0.3
@@ -286,6 +291,9 @@ function logmxp1(x::Float64)
     end
 end
 
+# Naive fallback
+logmxp1(x::Real) = (log(x) + one(x)) - x
+
 # The kernel of log1pmx
 # Accuracy within ~2ulps for -0.227 < x < 0.315
 function _log1pmx_ker(x::Float64)

From 5e6ae6ff88a9992e1a1cb570f6005e22862bc2b3 Mon Sep 17 00:00:00 2001
From: Simone Carlo Surace <51025924+simsurace@users.noreply.github.com>
Date: Fri, 29 Apr 2022 19:26:42 +0200
Subject: [PATCH 2/4] Use less naive heuristic

Co-authored-by: David Widmann <devmotion@users.noreply.github.com>
---
 src/basicfuns.jl | 10 +++++++++-
 1 file changed, 9 insertions(+), 1 deletion(-)

diff --git a/src/basicfuns.jl b/src/basicfuns.jl
index cc247319..5b989c37 100644
--- a/src/basicfuns.jl
+++ b/src/basicfuns.jl
@@ -292,7 +292,15 @@ function logmxp1(x::Float64)
 end
 
 # Naive fallback
-logmxp1(x::Real) = (log(x) + one(x)) - x
+function logmxp1(x::Real)
+    one_x = one(x)
+    if 2 * x < one_x
+        # for small values of `x` the other branch returns non-finite values
+        return (log(x) + one_x) - x
+    else
+        return log1pmx(x - one_x)
+    end
+end
 
 # The kernel of log1pmx
 # Accuracy within ~2ulps for -0.227 < x < 0.315

From c984a99fd35860ac48d5aa744a23bfd90a0c16f1 Mon Sep 17 00:00:00 2001
From: Simone Surace <simone.surace@artificialy.com>
Date: Mon, 2 May 2022 11:47:00 +0200
Subject: [PATCH 3/4] Add some tests

---
 test/basicfuns.jl | 16 ++++++++++++++++
 1 file changed, 16 insertions(+)

diff --git a/test/basicfuns.jl b/test/basicfuns.jl
index 4479080d..cccaf876 100644
--- a/test/basicfuns.jl
+++ b/test/basicfuns.jl
@@ -177,12 +177,28 @@ end
     @test iszero(log1pmx(0.0))
     @test log1pmx(1.0) ≈ log(2.0) - 1.0
     @test log1pmx(2.0) ≈ log(3.0) - 2.0
+
+    @test iszero(log1pmx(0f0))
+    @test log1pmx(1f0) ≈ log(2f0) - 1f0
+    @test log1pmx(2f0) ≈ log(3f0) - 2f0
+
+    for x in -0.5:0.1:10
+        @test log1pmx(Float32(x)) ≈ Float32(log1pmx(x))
+    end
 end
 
 @testset "logmxp1" begin
     @test iszero(logmxp1(1.0))
     @test logmxp1(2.0) ≈ log(2.0) - 1.0
     @test logmxp1(3.0) ≈ log(3.0) - 2.0
+
+    @test iszero(logmxp1(1f0))
+    @test logmxp1(2f0) ≈ log(2f0) - 1f0
+    @test logmxp1(3f0) ≈ log(3f0) - 2f0
+
+    for x in 0.1:0.1:10
+        @test logmxp1(Float32(x)) ≈ Float32(logmxp1(x))
+    end
 end
 
 @testset "logsumexp" begin

From 145d976bfded6c0d85769287595ed75d19998bc1 Mon Sep 17 00:00:00 2001
From: Simone Surace <simone.surace@artificialy.com>
Date: Mon, 2 May 2022 11:47:15 +0200
Subject: [PATCH 4/4] Bump version

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index c4b03a24..f1febd46 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "LogExpFunctions"
 uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
 authors = ["StatsFun.jl contributors, Tamas K. Papp <tkpapp@gmail.com>"]
-version = "0.3.12"
+version = "0.3.13"
 
 [deps]
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"