[tpapp]

Fixes #102.

Co-coded with Claude.

(Incidental trailing whitespace cleanup)

[devmotion]

I also did a literature search. The hi/lo reduction here is the standard Cody–Waite argument reduction (Cody & Waite 1980; Muller et al., Elementary Functions, range-reduction chapter). Some form of it is unavoidable, because near log 2 the result depends on bits of log 2 beyond the working precision. Other software handles irrational singular points the same way, with hi/lo constants: Julia Base exp (LogBo256U/L) and rem_pio2 (pio2_1/pio2_1t), and Boost.Math's digamma near its positive root (root1 + root2 + root3).

However, I'm not convinced by the implementation:

1. Not generic. Only Float16/32/64 are fixed. BigFloat, Float128, DoubleFloats etc. keep the old behaviour, e.g. BigFloat at 64 bits is off by 7e16 ulp just below log 2, and the error grows to Inf at higher precision. Also, unlike log1pexp/logexpm1/logcosh, there's no log2mexp(x::Real) = _log2mexp(float(x)) funnel.
2. @assume_effects :foldable on BigFloat/MPFR/setprecision code promises the compiler something that isn't true. It does fold today, but whenever it doesn't (interpreter, --compile=min, GPU/overlay method tables) we either run MPFR with allocations on every call or fail to compile. This file already has a robust pattern for this: hard-coded per-type constants plus a generic fallback (_log1pexp_thresholds_*).
3. NaN instead of DomainError for x > log 2. That's a behaviour change, and it's inconsistent with log, log1mexp and all other number types. The check isn't needed: log of a negative argument already throws.
4. _log2_hilo reads like base-2 log2.

Suggestion (generic, simple, no compiler heuristics):

log2mexp(x::Real) = _log2mexp(float(x))

function _log2mexp(x::Real)
    if x <= oftype(x, 0.5)
        return log1p(-expm1(x))
    else
        # log(2 - exp(x)) = log(-2 expm1(x - log(2))), without cancellation if x - log(2) is accurate
        return log(-2 * expm1(_minus_logtwo(x)))
    end
end

# x - log(2) (Cody & Waite): x - hi is exact (Sterbenz), only lo = log(2) - hi needs extra precision
_minus_logtwo(x::Real) = (x - oftype(x, IrrationalConstants.logtwo)) - _logtwo_lo(x)
# MPFR rounds x - l correctly, so a single subtraction of a more precise log(2) suffices
function _minus_logtwo(x::BigFloat)
    prec = precision(x) + precision(BigFloat) + 32
    return x - BigFloat(IrrationalConstants.logtwo; precision = prec)
end

_logtwo_lo(::Float64) = 2.3190468138462996e-17
_logtwo_lo(::Float32) = -1.9046542f-9
_logtwo_lo(::Float16) = Float16(-0.0002122)
function _logtwo_lo(x::Real)
    l = BigFloat(IrrationalConstants.logtwo; precision = 2 * precision(typeof(x)) + 32)
    return oftype(x, l - oftype(x, IrrationalConstants.logtwo))
end

This needs IrrationalConstants = "0.2.5" for the precision keyword.

Quick prototype numbers (max ulp error against a 2048-bit reference, last 2000 floats below log 2 / x ∈ [0.5, log 2)):

	master	this PR	suggestion
Float16	31.5 / 11.2	1.64 / 1.27	1.64 / 1.27
Float32	3.4e4 / 1.7e4	0.58 / 1.02	0.58 / 1.02
Float64	Inf / 1.4e13	0.52 / 1.11	0.52 / 1.11
BigFloat (64 … 1000 bits)	7e16 … Inf	unchanged	≤ 0.5

Float32/Float64 timings are unchanged (~9.5/13 ns). The generic _logtwo_lo reproduces the three literals exactly, so a test can compare them. The remaining 1.64 ulp for Float16 comes from the unchanged x ≤ 0.5 branch (double rounding). Float16(log2mexp(Float32(x))) would give 0.5 ulp, but that's a separate issue.

Side note, not really part of this PR: for Float16/Float32/Float64 there's no performance issue. oftype(x, logtwo) is a literal cached by the @irrational macro (Float16 is converted from the Float32 value), and _logtwo_lo is hard-coded. Only BigFloat (and other generic types) are affected: IrrationalConstants defines logtwo via log(big(2)), so the high-precision value is recomputed on every call (~2× slower for x > 0.5 at 256 bits). Upstream deliberately moved away from MPFR's cached mpfr_const_log2 because of MPFR breakage on Julia master. Maybe that change should be reverted at some point.

It might also be worth a docstring note: near log(2) the function is ill-conditioned. Users who have d = log(2) − x directly should use logtwo + log1mexp(-d).

[devmotion]

:foldable asserts :consistent, :effect_free and :terminates, but this allocates BigFloats with finalizers and calls MPFR under setprecision. If const-eval doesn't happen, this runs on every call. If the assumption is wrong, it's UB. Hard-coded constants (like _log1pexp_thresholds_*) avoid all of this. Also: log2 in the name suggests the base-2 logarithm; maybe _logtwo_*? And 2p + a few guard bits is enough, rather than 4p + 64.

[devmotion]

Restricting this to Union{Float16,Float32,Float64} leaves BigFloat & co. (and the float(x) funnel for Int/Rational) on the inaccurate path. The algorithm itself is type-agnostic; only the lo constant is type-specific.

[devmotion]

This changes the behaviour for x > log 2 from DomainError to NaN, only for these three types. I'd drop the check: log(-2 * expm1(-d)) already throws a DomainError for d < 0.

[devmotion]

• ≈ uses rtol = √eps (0.03 for Float16), which is weak.
• The reference reuses the code under test; log(2 - exp(big(x))) would be independent.
• Coverage is missing for m = 0 (T(log(2)) lies above log 2 for Float16/Float32 and below for Float64), x within 1 ulp, the DomainError for x > log 2, and BigFloat.
• Could use the existing ulp_error/ulp_error_maximum helpers, e.g. exhaustively for Float16.
• The BigFloat accuracy testset scales tolerances by the condition number, which is exactly why it never caught log2mexp(x) accuracy vanishes around the right edge of the domain, at log(2) #102. Maybe add an unscaled check near log 2.

[tpapp]

Thanks, please close this and submit your code as a PR, you can do a much better job of this.

[tpapp]

Closing in favor of #140.

@devmotion: Incidentally, I sent you a PM on the Julia discourse forum, asking about Claude, wonder if you have seen it.