diff --git a/Project.toml b/Project.toml
index c9b1c1a0a..84090560e 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "Distributions"
 uuid = "31c24e10-a181-5473-b8eb-7969acd0382f"
 authors = ["JuliaStats"]
-version = "0.25.94"
+version = "0.25.95"
 
 [deps]
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
diff --git a/src/univariate/continuous/normal.jl b/src/univariate/continuous/normal.jl
index f652a27d7..1dcadadac 100644
--- a/src/univariate/continuous/normal.jl
+++ b/src/univariate/continuous/normal.jl
@@ -89,6 +89,15 @@ end
 # Use Julia implementations in StatsFuns
 @_delegate_statsfuns Normal norm μ σ
 
+# `logerf(...)` is more accurate for arguments in the tails than `logsubexp(logcdf(...), logcdf(...))`
+function logdiffcdf(d::Normal, x::Real, y::Real)
+    x < y && throw(ArgumentError("requires x >= y."))
+    μ, σ = params(d)
+    _x, _y, _μ, _σ = promote(x, y, μ, σ)
+    s = sqrt2 * _σ
+    return logerf((_y - _μ) / s, (_x - _μ) / s) - logtwo
+end
+
 gradlogpdf(d::Normal, x::Real) = (d.μ - x) / d.σ^2
 
 mgf(d::Normal, t::Real) = exp(t * d.μ + d.σ^2 / 2 * t^2)
diff --git a/test/univariate/continuous/normal.jl b/test/univariate/continuous/normal.jl
index 6f494d938..ccc6266db 100644
--- a/test/univariate/continuous/normal.jl
+++ b/test/univariate/continuous/normal.jl
@@ -1,4 +1,4 @@
-using Test, Distributions, ForwardDiff
+using Test, Distributions, StatsFuns, ForwardDiff
 
 isnan_type(::Type{T}, v) where {T} = isnan(v) && v isa T
 
@@ -17,9 +17,19 @@ isnan_type(::Type{T}, v) where {T} = isnan(v) && v isa T
     @test -Inf === logpdf(Normal(), Inf)
     @test iszero(logcdf(Normal(0, 0), 0))
     @test iszero(logcdf(Normal(), Inf))
-    @test logdiffcdf(Normal(), Float32(5), Float32(3)) ≈ -6.607938594596893 rtol=1e-12
-    @test logdiffcdf(Normal(), Float32(5), Float64(3)) ≈ -6.607938594596893 rtol=1e-12
-    @test logdiffcdf(Normal(), Float64(5), Float64(3)) ≈ -6.607938594596893 rtol=1e-12
+    @test @inferred(logdiffcdf(Normal(), 5f0, 3f0)) ≈ -6.607938594596893 rtol=1e-12
+    @test @inferred(logdiffcdf(Normal(), 5f0, 3.0)) ≈ -6.607938594596893 rtol=1e-12
+    @test @inferred(logdiffcdf(Normal(), 5.0, 3.0)) ≈ -6.607938594596893 rtol=1e-12
+    @test_throws ArgumentError logdiffcdf(Normal(), 3, 5)
+
+    # Arguments in the tails
+    logdiffcdf_big(d::Normal, x::Real, y::Real) = logsubexp(logcdf(d, big(y)), logcdf(d, big(x)))
+    for d in (Normal(), Normal(2.1, 0.1)), (a, b) in ((15, 10), (115, 100), (1015, 1000))
+        for (x, y) in ((a, b), (-b, -a))
+            @test isfinite(@inferred(logdiffcdf(d, x, y)))
+            @test logdiffcdf(d, x, y) ≈ logdiffcdf_big(d, x, y)
+        end
+    end
     let d = Normal(Float64(0), Float64(1)), x = Float64(-60), y = Float64(-60.001)
         float_res = logdiffcdf(d, x, y)
         big_x = BigFloat(x; precision=100)