logerfcx (#203)

Author: cossio

docs/src/functions_list.md

@@ -9,6 +9,7 @@ SpecialFunctions.erf
 SpecialFunctions.erfc
 SpecialFunctions.erfcx
 SpecialFunctions.logerfc
+SpecialFunctions.logerfcx
 SpecialFunctions.erfi
 SpecialFunctions.dawson
 SpecialFunctions.erfinv

docs/src/functions_overview.md

@@ -38,6 +38,7 @@ Here the *Special Functions* are listed according to the structure of [NIST Digi
 | [`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`` |
+| [`logerfcx(x)`](@ref SpecialFunctions.logerfcx) | log of the scaled complementary error function, i.e. accurate ``\operatorname{ln}(\operatorname{erfcx}(x))`` for large negative ``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

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

src/erf.jl

@@ -458,3 +458,35 @@ end
 
 logerfc(x::Real) = logerfc(float(x))
 logerfc(x::AbstractFloat) = throw(MethodError(logerfc, x))
+
+@doc raw"""
+    logerfcx(x)
+
+Compute the natural logarithm of the scaled complementary error function of ``x``, that is
+
+```math
+\operatorname{logerfcx}(x) = \operatorname{ln}(\operatorname{erfcx}(x))
+\quad \text{for} \quad x \in \mathbb{R} \, .
+```
+
+This is the accurate version of ``\operatorname{ln}(\operatorname{erfcx}(x))`` for large and negative ``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 logerfcx(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(erfc(x)) + x^2
+    else
+        return log(erfcx(x))
+    end
+end
+
+logerfcx(x::Real) = logerfcx(float(x))
+logerfcx(x::AbstractFloat) = throw(MethodError(logerfcx, x))

test/erf.jl

@@ -22,6 +22,22 @@
         @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 logerfcx(Float16(1))
+        @test_throws MethodError logerfcx(Float16(-1))
+        @test iszero(logerfcx(0))
+        @test logerfcx(Float32(1)) ≈ -0.849605509933248248576017509499 rtol=2eps(Float32)
+        @test logerfcx(Float64(1)) ≈ -0.849605509933248248576017509499 rtol=2eps(Float32)
+        @test logerfcx(Float32(-1)) ≈ 1.61123231767807049464268192445 rtol=2eps(Float32)
+        @test logerfcx(Float64(-1)) ≈ 1.61123231767807049464268192445 rtol=2eps(Float32)
+        @test logerfcx(Float32(-100)) ≈ 10000.6931471805599453094172321 rtol=2eps(Float32)
+        @test logerfcx(Float64(-100)) ≈ 10000.6931471805599453094172321 rtol=2eps(Float64)
+        @test logerfcx(Float32(100)) ≈ -5.17758512266433257046678208395 rtol=2eps(Float32)
+        @test logerfcx(Float64(100)) ≈ -5.17758512266433257046678208395 rtol=2eps(Float64)
+        @test logerfcx(Float32(-1000)) ≈ 1.00000069314718055994530941723e6 rtol=2eps(Float32)
+        @test logerfcx(Float64(-1000)) ≈ 1.00000069314718055994530941723e6 rtol=2eps(Float64)
+        @test logerfcx(Float32(1000)) ≈ -7.48012072190621214066734919080 rtol=2eps(Float32)
+        @test logerfcx(Float64(1000)) ≈ -7.48012072190621214066734919080 rtol=2eps(Float64)
+
         @test_throws MethodError erfi(Float16(1))
         @test erfi(Float32(1)) ≈ 1.6504257587975428760 rtol=2*eps(Float32)
         @test erfi(Float64(1)) ≈ 1.6504257587975428760 rtol=2*eps(Float64)