[devmotion]

This PR fixes several numerical problems in the ChainRules definitions, makes some of them faster, and adds rules for the four functions that had none.

Each new formula was compared against 256-bit BigFloat references in Float64 (and in Float32 for logistic). Errors are given in ulps and timings in ns per call (Apple Silicon, Julia 1.13). The benchmarks used standalone copies of the formulas, not the extension code itself.

Rule	Old → new	Max error (ulp)	Time (ns)
logistic	Ω * (1 - Ω) → e / (1 + e)^2, e = exp(-abs(x))	∞ (returns 0 at x = 40) → 3.7; Float32: 1.7e7 → 3.3	4.4 → 3.8
logabstanh	inv(cosh(x) * sinh(x)) → 2 / sinh(2x)	3.2 → 2.2	14.0 → 7.2
log1psq	2x / (1 + x^2) → same for abs(x) < 1, else 2 / (x + inv(x))	overflows for |x| ≳ 1e154 → 1.4	–
log1mexp	-exp(x - Ω) → -inv(expm1(-x))	31 → 1.5	24.5 → 23.9
logexpm1	exp(x - Ω) → -inv(expm1(-x))	25 → 1.8; NaN at Inf → 1	21.1 → 19.5
logaddexp	exp(x - Ω), exp(y - Ω) → a single logistic call	511 → 2.0; NaN at (Inf, 1) → (1, 0)	27.6 → 24.9
logsubexp	branched ±exp(· - Ω) → a single expm1 call	386 → 1.5; NaN at (Inf, 1) → (1, -0)	29.0 → 23.2

logaddexp and logsubexp. Both satisfy f(x + c, y + c) = f(x, y) + c, so ∂x + ∂y = 1. The rules compute the partial t for the smaller argument directly (logistic(-Δ) or -inv(expm1(Δ)), where Δ = |x - y|) and set the other partial to 1 - t. Since t ≤ 1/2, 1 - t does not suffer from cancellation. As in the primal implementations, Δ is set to zero if x - y is NaN (x == y == ±Inf for logaddexp, x == y == -Inf for logsubexp). Unlike the primal implementations, this replacement is not applied for finite x == y: there x - y is already zero, and replacing it would drop higher-order derivatives when the rules are evaluated on dual or tracked numbers whose == only compares primal values.

Behaviour changes at x == y, where the derivative isn't defined:

• logaddexp(±Inf, ±Inf) now returns (1/2, 1/2) instead of NaN (the ChainRules counterpart of logaddexp introduces NaNs for non-finite values #61).
• logsubexp(x, x) now returns (Inf, -Inf), the limit as x approaches y from above, instead of (-Inf, Inf).

New rules. These functions had no rules: loglogistic (logistic(-x)), log1mlogistic (-logistic(x)), logitexp (-inv(expm1(x))) and logit1mexp (inv(expm1(x))).

Unchanged. The log2mexp rule stays as it is: the function is ill-conditioned near log(2), and the expm1-based alternative is no more accurate. The cloglog and cexpexp rules are now restricted to Real arguments; their formulas are unchanged.

Tests. I added the following tests:

• test_scalar for the four new rules.
• Comparisons against BigFloat references with rtol = 4eps. The old rules fail each of these:
  • logistic at large |x|, in Float64 and Float32
  • log1mexp(-1e-12)
  • logexpm1(6.04e-8)
  • logaddexp(300, 301)
  • logsubexp with arguments 1e-6 apart, where finite differences are unreliable
  • log1psq(1e200), which overflowed before
• Cases at Inf, checked by hand: logexpm1, and logaddexp/logsubexp with one infinite argument.
• Equal and NaN arguments of logaddexp and logsubexp, in Float64 and Float32: logaddexp(x, x) for x = -Inf, 1, Inf; logsubexp(x, x) for x = -Inf, 1 (and NaN for x = Inf); NaN propagation.

The full test suite passes locally.

The same rules are ported to DiffRules (used by ForwardDiff and ReverseDiff) in JuliaDiff/DiffRules.jl#121, which fixes #61.

🤖 Generated with Claude Code

[tpapp]

I would keep the original comments in the code, more or less, and if the concern is NaNs, just code it directly, instead of figuring out what will be NaN from the inputs.

[tpapp]

why not simply

Suggested change

	Δ = x == y && isinf(x) ? zero(x - y) : abs(x - y)
	Δ = abs(x - y)
	if isnan(Δ)
	Δ = zero(Δ)
	end

which conveys the intention?

[devmotion]

If x or y are NaN, I think one could argue that it should return NaN, not 0. Hence I made the zero branch a bit more restrictive.

[tpapp]

similarly as above

[tpapp]

lgtm