Precision loss for BigFloat in logexpm1, logcosh and log1pmx

Open
  • Is something wrong?
#135 0 comments 0 reactions 0 assignees View on GitHub

@devmotion is already working on this.

Since Nov 9, 2026.

  • #136 by @devmotion — open

Assessment

Difficulty
4/5
Estimated time
3-5 days
Newbie friendliness
30/100
Issue type
Bug
Clarity
Clearly specified
Activity status
Stale
Tech stack
julia
Domain
backend

Research direction

Start by reading the implementations of logexpm1, logcosh, and log1pmx, along with the wrappers named in the issue. Run the provided BigFloat comparison script to reproduce the precision losses and check the existing tests for these functions. Done means the affected functions and wrappers retain precision across the reported input ranges; pull request #136 is already open, so check its work before starting.

Written by the indexing model from the issue text.

Description

All exported single-argument functions accept BigFloat, but a few of them lose a large part of the precision because the implementations contain thresholds or branch choices that are tuned for Float64.

Method

Each function was evaluated at the default 256-bit BigFloat precision on a grid of inputs spanning its domain (tiny, moderate and huge magnitudes). The reference is the naive mathematical definition evaluated at 4096 bits. "Digits lost" is log10(relerr / eps(BigFloat)), so 0 means correctly rounded to within a few ulps and 77 means no correct digits. The condition number of each function was computed as well and subtracted out, so the numbers below are algorithmic loss only, not inherent ill-conditioning (e.g. near roots). Measured on master (5346ca8).

Results

Function Digits lost (algorithmic) Where
xlogx 0
xexpx 0
logistic 0
logit 0
log1psq 0
log1pexp 0
log1mexp 0
log2mexp 0
logexpm1 up to 61 18 < x < 170
softplus 0
invsoftplus up to 61 18 < x < 170
log1pmx 0.7 at 0.1, 2 at 1e-3, 10 at 1e-10, 30 at 1e-30, all 77 at 1e-100 small |x|
logmxp1 0
logcosh 1.4 at 0.1, 6 at 1e-3, 19 at 1e-10, all 77 at 1e-100 small |x|
logabssinh 0
logabstanh 0
cloglog 0
cexpexp 0
loglogistic 0
logitexp up to 61 -170 < x < -18
log1mlogistic 0
logit1mexp up to 61 -170 < x < -18

Causes

Three distinct problems account for every nonzero row:

  • logexpm1 has Float64 thresholds baked in. Above 18 it switches to x - exp(-x), which drops a term of size exp(-2x)/2, and above 33.3 it returns x unchanged, dropping exp(-x). Both cutoffs are right for 53 bits and wrong for 256. The loss peaks around 61 digits near x = 34 and only fades to zero past x ≈ 170, where exp(-x) drops below the 256-bit epsilon. invsoftplus, logitexp and logit1mexp are thin wrappers around logexpm1 and inherit the problem.
  • logcosh has no small-argument kernel for non-hardware floats. The polynomial kernel is gated on Float16/Float32/Float64; everything else evaluates |x| + log1pexp(-2|x|) - log(2), where the two large terms cancel. The loss grows as about 2*log10(1/|x|) digits.
  • log1pmx falls back to log1p(x) - x, which loses about log10(1/|x|) digits for small x. (logmxp1's generic fallback is fine because of its branch at x = 1/2.) The analogous question for Float32 and Float16 is tracked in #46.

Possible fixes: precision-aware thresholds in logexpm1 based on eps(x) instead of constants, and a series or direct log(cosh(x)) branch in logcosh and a series for log1pmx for types without a dedicated kernel.

Script 
using LogExpFunctions, Printf, ForwardDiff

# (function, reference definition, test points)
L(a) = BigFloat.(a)
small = [1e-100, 1e-30, 1e-10, 1e-3]
cases = [
 (xlogx,        x -> iszero(x) ? zero(x) : x*log(x),     L([0; small; 0.5; 1; 2; 1e10; 1e100])),
 (xexpx,        x -> x*exp(x),                           L([-1e100; -1e3; -10; -1; -small; 0; small; 1; 10; 500])),
 (logistic,     x -> 1/(1+exp(-x)),                      L([-1e5; -800; -40; -1; -small; 0; small; 1; 40; 800; 1e5])),
 (logit,        x -> log(x/(1-x)),                       L([1e-100; 1e-10; 0.1; 0.5; 0.5+1e-20; 0.9; 1-1e-10; 1-1e-15])),
 (log1psq,      x -> log1p(x^2),                         L([0; small; 1; 10; 1e10; 1e100])),
 (log1pexp,     x -> log1p(exp(x)),                      L([-1e5; -800; -40; -1; 0; 1; 40; 800; 1e5])),
 (log1mexp,     x -> log1p(-exp(x)),                     L([-1e5; -800; -40; -1; -0.5; -small])),
 (log2mexp,     x -> log(2-exp(x)),                      L([-1e5; -40; -1; -small; 0; small; 0.5; log(big(2))-1e-20])),
 (logexpm1,     x -> log(expm1(x)),                      L([small; 0.5; 1; 40; 800; 1e5])),
 (softplus,     x -> log1p(exp(x)),                      L([-1e5; -800; -40; -1; 0; 1; 40; 800; 1e5])),
 (invsoftplus,  x -> log(expm1(x)),                      L([small; 0.5; 1; 40; 800; 1e5])),
 (log1pmx,      x -> log1p(x)-x,                         L([-0.9; -0.5; -small; 0; small; 0.5; 1; 10; 1e10])),
 (logmxp1,      x -> log(x)-x+1,                         L([1e-10; 0.1; 0.5; 1-1e-20; 1; 1+1e-20; 2; 10; 1e10])),
 (logcosh,      x -> log(cosh(x)),                       L([-1e5; -40; -1; -small; 0; small; 1; 40; 1e5])),
 (logabssinh,   x -> log(abs(sinh(x))),                  L([-1e5; -40; -1; -small; small; 1; 40; 1e5])),
 (logabstanh,   x -> log(abs(tanh(x))),                  L([-1e5; -40; -1; -small; small; 1; 40; 1e5])),
 (cloglog,      x -> log(-log(1-x)),                     L([1e-100; 1e-10; 0.1; 0.5; 0.9; 1-1e-10; 1-1e-20])),
 (cexpexp,      x -> 1-exp(-exp(x)),                     L([-1e5; -800; -40; -1; 0; 1; 5; 40])),
 (loglogistic,  x -> -log1p(exp(-x)),                    L([-1e5; -800; -40; -1; 0; 1; 40; 800; 1e5])),
 (logitexp,     x -> log(exp(x)/(1-exp(x))),             L([-1e5; -800; -40; -1; -log(big(2)); -0.5; -small])),
 (log1mlogistic,x -> -log1p(exp(x)),                     L([-1e5; -800; -40; -1; 0; 1; 40; 800; 1e5])),
 (logit1mexp,   x -> log((1-exp(x))/exp(x)),             L([-1e5; -800; -40; -1; -log(big(2)); -0.5; -small])),
]

const REFPREC = 4096
println("| function | works | return type | max digits lost | of which inherent (cond) | excess (algorithmic) | worst input |")
println("|---|---|---|---|---|---|---|")
for (f, ref, xs) in cases
    worst = -Inf; worstx = nothing; rt = nothing; err = nothing; worstcond = 0.0; worstlost = 0.0
    for x in xs
        y = try f(x) catch e; err = e; break end
        rt = typeof(y)
        r = setprecision(REFPREC) do; ref(BigFloat(x)) end
        isfinite(r) && !iszero(r) || continue
        relerr = Float64(abs(big(y) - r) / abs(r))
        lost = relerr == 0 ? 0.0 : log10(relerr / Float64(eps(BigFloat)))
        cond = setprecision(REFPREC) do
            xb = BigFloat(x); d = ForwardDiff.derivative(ref, xb)
            iszero(xb) ? 0.0 : max(0.0, log10(Float64(abs(xb*d/r))))
        end
        excess = max(0.0, lost - cond)
        if excess > worst; worst = excess; worstx = x; worstcond = cond; worstlost = lost; end
    end
    if err !== nothing
        @printf("| `%s` | no (%s) | | | |\n", f, nameof(typeof(err)))
    else
        wx = worstx === nothing ? "" : @sprintf("%.3g", Float64(worstx))
        @printf("| `%s` | yes | %s | %.1f | %.1f | %.1f | %s |\n", f, rt, max(worstlost, 0.0), worstcond, max(worst, 0.0), wx)
    end
end

This analysis and the text of this issue were written by Claude Code at my request. I have reviewed it.

Dominant language
Julia
Stars
86
Forks
25
Avg merge
4d 7h
Merged PRs (30d)
3

Getting set up

This project ships no dev container, Dockerfile or contributing guide, so setting up is up to you: start from its README, and see our first-contribution guide for the general steps.

First steps

  1. Read the whole issue, then the project's contributing guide.
  2. Comment on the issue to say you are picking it up — it saves two people doing the same work.
  3. Fork the repository and make your change on a branch.
  4. Open a pull request that references the issue number.

More from JuliaStats/LogExpFunctions.jl

All issues in JuliaStats/LogExpFunctions.jl

Similar issues

More Julia issues

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.