add logerf(a, b) (#229)

Project.toml

@@ -1,6 +1,6 @@
 name = "SpecialFunctions"
 uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
-version = "0.10.2"
+version = "0.10.3"
 
 [deps]
 OpenSpecFun_jll = "efe28fd5-8261-553b-a9e1-b2916fc3738e"

src/SpecialFunctions.jl

@@ -36,6 +36,7 @@ export
     erfcx,
     erfi,
     erfinv,
+    logerf,
     logerfc,
     logerfcx,
     eta,

src/erf.jl

@@ -506,3 +506,27 @@ end
 
 logerfcx(x::Real) = logerfcx(float(x))
 logerfcx(x::AbstractFloat) = throw(MethodError(logerfcx, x))
+
+@doc raw"""
+    logerf(x, y)
+
+Compute the natural logarithm of two-argument error function. This is an accurate version of
+ `log(erf(x, y))`, which works for large `x, y`.
+
+External links: [Wikipedia](https://en.wikipedia.org/wiki/Error_function).
+
+See also: [`erf(x,y)`](@ref SpecialFunctions.erf).
+"""
+function logerf(a::Real, b::Real)
+    if abs(a) ≤ 1/√2 && abs(b) ≤ 1/√2
+        return log(erf(a, b))
+    elseif b > a > 0
+        return logerfc(a) + log1mexp(logerfc(b) - logerfc(a))
+    elseif a < b < 0
+        return logerfc(-b) + log1mexp(logerfc(-a) - logerfc(-b))
+    else
+        return log(erf(a, b))
+    end
+end
+
+log1mexp(x::Real) = x < -log(oftype(x, 2)) ? log1p(-exp(x)) : log(-expm1(x))

test/erf.jl

@@ -141,5 +141,22 @@
         @test isnan(erf(NaN, 0))
         @test isnan(erf(0, NaN))
         @test isnan(erf(NaN, NaN))
+
+        @test logerf(-5, 35) ≈ 0.693147180559176579520017808560
+        @test logerf(30, 35) ≈ -903.974117110643878079600243618
+        @test logerf(-35, -30) ≈ -903.974117110643878079600243618
+        @test logerf(10, Inf) ≈ -102.879889024844888574804787140
+        @test logerf(-Inf, Inf) ≈ 0.693147180559945309417232121458
+        @test logerf(Inf, Inf) == -Inf
+        @test logerf(-Inf, -Inf) == -Inf
+        @test logerf(-Inf, Inf) ≈ log(2)
+        @test logerf(-1e-6, 1e-6) ≈ -13.0015811397694169056785314498
+        @test isnan(logerf(NaN, 0))
+        @test isnan(logerf(0, NaN))
+        @test isnan(logerf(NaN, NaN))
+        @test logerf(-1e-30, 1e-30) ≈ -68.2636233716261799887769930733
+        @test logerf(1e-30, 2e-30) ≈ -68.9567705521861252981942251947   
+        @test logerf(-2e-30, -1e-30) ≈ -68.9567705521861252981942251947
+        @test_throws DomainError logerf(2, 1)
     end
 end