Add logerfc function (#194)

* Add prototype logerfc function.
* Add unit tests (and change src/erf.jl argument types).
* Add docstring.

Author: c123w

docs/src/functions_list.md

@@ -8,6 +8,7 @@ CurrentModule = SpecialFunctions
 SpecialFunctions.erf
 SpecialFunctions.erfc
 SpecialFunctions.erfcx
+SpecialFunctions.logerfc
 SpecialFunctions.erfi
 SpecialFunctions.dawson
 SpecialFunctions.erfinv

docs/src/functions_overview.md

@@ -37,6 +37,7 @@ Here the *Special Functions* are listed according to the structure of [NIST Digi
 | [`erfc(x)`](@ref SpecialFunctions.erfc)       | complementary error function, i.e. the accurate version of ``1-\operatorname{erf}(x)`` for large ``x`` |
 | [`erfcinv(x)`](@ref SpecialFunctions.erfcinv) | inverse function to [`erfc()`](@ref SpecialFunctions.erfc) |
 | [`erfcx(x)`](@ref SpecialFunctions.erfcx)     | scaled complementary error function, i.e. accurate ``e^{x^2} \operatorname{erfc}(x)`` for large ``x`` |
+| [`logerfc(x)`](@ref SpecialFunctions.logerfc) | log of the complementary error function, i.e. accurate ``\operatorname{ln}(\operatorname{erfc}(x))`` for large ``x`` |
 | [`erfi(x)`](@ref SpecialFunctions.erfi)       | imaginary error function defined as ``-i \operatorname{erf}(ix)`` |
 | [`erfinv(x)`](@ref SpecialFunctions.erfinv)   | inverse function to [`erf()`](@ref SpecialFunctions.erf) |
 | [`dawson(x)`](@ref SpecialFunctions.dawson)   | scaled imaginary error function, a.k.a. Dawson function, i.e. accurate ``\frac{\sqrt{\pi}}{2} e^{-x^2} \operatorname{erfi}(x)`` for large ``x`` |

src/SpecialFunctions.jl

@@ -34,6 +34,7 @@ export
     erfcx,
     erfi,
     erfinv,
+    logerfc,
     eta,
     digamma,
     invdigamma,
@@ -64,7 +65,7 @@ include("betanc.jl")
 include("beta_inc.jl")
 include("deprecated.jl")
 
-for f in (:digamma, :erf, :erfc, :erfcinv, :erfcx, :erfi, :erfinv,
+for f in (:digamma, :erf, :erfc, :erfcinv, :erfcx, :erfi, :erfinv, :logerfc,
           :eta, :gamma, :invdigamma, :logfactorial, :lgamma, :trigamma, :ellipk, :ellipe)
     @eval $(f)(::Missing) = missing
 end

src/erf.jl

@@ -426,3 +426,36 @@ function erfcx(x::BigFloat)
         return (1+s)/(x*sqrt(oftype(x,pi)))
     end
 end
+
+@doc raw"""
+    logerfc(x)
+
+Compute the natural logarithm of the complementary error function of ``x``, that is
+
+```math
+\operatorname{logerfc}(x) = \operatorname{ln}(\operatorname{erfc}(x))
+\quad \text{for} \quad x \in \mathbb{R} \, .
+```
+
+This is the accurate version of ``\operatorname{ln}(\operatorname{erfc}(x))`` for large ``x``.
+
+External links:
+[Wikipedia](https://en.wikipedia.org/wiki/Error_function).
+
+See also: [`erfcx(x)`](@ref SpecialFunctions.erfcx).
+
+# Implementation
+Based on the [`erfc(x)`](@ref SpecialFunctions.erfc) and [`erfcx(x)`](@ref SpecialFunctions.erfcx) functions.
+Currently only implemented for `Float32`, `Float64`, and `BigFloat`.
+"""
+function logerfc(x::Union{Float32, Float64, BigFloat})
+    # Don't include Float16 in the Union, otherwise logerfc would currently work for x <= 0.0, but not x > 0.0
+    if x > 0.0
+        return log(erfcx(x)) - x^2
+    else
+        return log(erfc(x))
+    end
+end
+
+logerfc(x::Real) = logerfc(float(x))
+logerfc(x::AbstractFloat) = throw(MethodError(logerfc, x))

test/erf.jl

@@ -9,8 +9,18 @@
         @test erfc(Float64(1))  ≈ 0.15729920705028513066 rtol=2*eps(Float64)
 
         @test_throws MethodError erfcx(Float16(1))
-        @test erfcx(Float32(1)) ≈ 0.42758357615580700442    rtol=2*eps(Float32)
-        @test erfcx(Float64(1)) ≈ 0.42758357615580700442    rtol=2*eps(Float64)
+        @test erfcx(Float32(1)) ≈ 0.42758357615580700442 rtol=2*eps(Float32)
+        @test erfcx(Float64(1)) ≈ 0.42758357615580700442 rtol=2*eps(Float64)
+
+        @test_throws MethodError logerfc(Float16(1))
+        @test_throws MethodError logerfc(Float16(-1))
+        @test logerfc(Float32(-100)) ≈ 0.6931471805599453 rtol=2*eps(Float32)
+        @test logerfc(Float64(-100)) ≈ 0.6931471805599453 rtol=2*eps(Float64)
+        @test logerfc(Float32(1000)) ≈ -1.0000074801207219e6 rtol=2*eps(Float32)
+        @test logerfc(Float64(1000)) ≈ -1.0000074801207219e6 rtol=2*eps(Float64)
+        @test logerfc(1000) ≈ -1.0000074801207219e6 rtol=2*eps(Float32)
+        @test logerfc(Float32(10000)) ≈ log(erfc(BigFloat(10000, 100))) rtol=2*eps(Float32)
+        @test logerfc(Float64(10000)) ≈ log(erfc(BigFloat(10000, 100))) rtol=2*eps(Float64)
 
         @test_throws MethodError erfi(Float16(1))
         @test erfi(Float32(1)) ≈ 1.6504257587975428760 rtol=2*eps(Float32)
@@ -69,6 +79,9 @@
         @test erfcx(BigFloat(1.8e88))   ≈ erfcx(1.8e88)             rtol=4*eps()
         @test isnan(erfcx(BigFloat(NaN)))
 
+        @test logerfc(BigFloat(1000, 100)) ≈ -1.0000074801207219e6 rtol=2*eps(Float64)
+        @test isnan(logerfc(BigFloat(NaN)))
+
         @test_throws MethodError erfi(big(1.0))
 
         @test_throws MethodError erfinv(BigFloat(1))