[devmotion]

Fixes #102. Supersedes #139: this is the generic version suggested in the review there, with tests.

Changes

• For x ≤ log(3/2), log2mexp(x) still computes log1p(-expm1(x)).
• For x > log(3/2), it computes log(-2 * expm1(x - log(2))). x - log(2) uses Cody–Waite argument reduction:
  • x - hi with hi = T(log(2)) is exact (Sterbenz lemma, since x ≥ hi/2), and lo = log(2) - hi is subtracted separately.
  • lo is hard-coded for Float16/32/64. Other types use a generic fallback, and a test checks that the hard-coded values match it.
  • BigFloat subtracts a higher-precision log(2) instead, because MPFR rounds that subtraction correctly.
• Threshold log(3/2): at u = -expm1(x) = -1/2 both formulas are equally well conditioned. This cutoff gives lower maximum errors than 0.5 for all three float types.
• No change to the domain (x ≤ T(log(2))): T(log(2)) rounds above log(2) for Float16/Float32. It is what log(T(2)) and log2mexp(T(-Inf)) return, so it still maps to -Inf, as before. This also keeps the InverseFunctions round trip working.
  • For Float64, T(log(2)) is below log(2), so log2mexp(log(2.0)) now returns the correct ≈ -37.6 instead of -Inf.
• Inputs outside the domain throw a DomainError with a clearer message. Before, the error came from log. NaN still returns NaN.

Accuracy

Maximum error in ulps, measured against a 2048-bit reference:

	new
Float16, all finite x ≤ log(2)	1.29
Float32, last 2000 floats ≤ log(2) / [log(1.5), log(2)] / [-20, 0.5]	0.58 / 1.25 / 1.33
Float64, same ranges as Float32	0.52 / 1.72 / 1.53
BigFloat (64, 256, 1000 bits) near log(2)	≤ 1 (master: 7e16 to 2e297)

Performance

• Float32/Float64: about the same as master, within roughly 0.7 ns per call.
• BigFloat with x > log(3/2): about 2× slower, because log(2) is recomputed at higher precision. On master the results in this range were wrong.

Version bumped to 1.0.6. The IrrationalConstants compat (0.1, 0.2) doesn't need to change: its lowest supported version already provides BigFloat(logtwo; precision=...).

🤖 Generated with Claude Code

[tpapp]

Thanks, this looks good to me.

That said, I am not sure that improving the precision of this function is the right way to tackle the problem. Apparently log2mexp is used in two places: the log(c)cdf of the Laplace function, and in tdistinvlogcdf in StatsFuns.

I am pretty sure that the first application at least should be handled as log1mexp(-z-loghalf), before rearranging the formula to use log2mexp. I am still looking at the second one.

[devmotion]

Since it's part of the public API of the package, I guess it would be good to improve accuracy. But indeed, it might not be the most relevant part of the package and of course downstream packages shouldn't use it if there's numerically better behaved alternatives.