diff --git a/src/LogExpFunctions.jl b/src/LogExpFunctions.jl index 739e005..bc5fefd 100644 --- a/src/LogExpFunctions.jl +++ b/src/LogExpFunctions.jl @@ -8,7 +8,7 @@ import LinearAlgebra export xlogx, xlogy, xlog1py, xexpx, xexpy, logistic, logit, log1psq, log1pexp, log1mexp, log2mexp, logexpm1, softplus, invsoftplus, log1pmx, logmxp1, logaddexp, logsubexp, logsumexp, logsumexp!, softmax, - softmax!, logcosh, logabssinh, cloglog, cexpexp, + softmax!, logcosh, logabssinh, logabstanh, cloglog, cexpexp, loglogistic, logitexp, log1mlogistic, logit1mexp # expm1(::Float16) is not defined in older Julia versions, diff --git a/src/basicfuns.jl b/src/basicfuns.jl index e013adf..e2828f1 100644 --- a/src/basicfuns.jl +++ b/src/basicfuns.jl @@ -149,6 +149,31 @@ end """ $(SIGNATURES) +Return `log(abs(tanh(x)))`, carefully evaluated without intermediate calculation of `tanh(x)`. + +The implementation ensures `logabstanh(-x) = logabstanh(x)`. +""" +function logabstanh(x::Real) + return log1p(-2/(exp(2*abs(x))+1)) +end +function logabstanh(x::Float32) + abs_x = abs(x) + if abs_x < 0.0625f0 + return log(abs_x) - x^2/3 + end + return log1p(-2/(exp(2*abs_x)+1)) +end +function logabstanh(x::Float64) + abs_x = abs(x) + if abs_x < 0x1p-5 + return log(abs_x) + evalpoly(x*x, (0, -1/3, 7/90, -62/2835)) + end + return log1p(-2/(exp(2*abs_x)+1)) +end + +""" +$(SIGNATURES) + Return `log(1+x^2)` evaluated carefully for `abs(x)` very small or very large. """ log1psq(x::Real) = log1p(abs2(x))