[DimitriChrysafis]

logaddexp computes the second input's derivative as 1 - sigmoid(a - b).
When the sigmoid rounds to one, this drops a derivative that is still
representable. Swapping the inputs changes gradient accuracy even though the
forward operation is symmetric.

import mlx.core as mx

x = mx.array(-20.0)
mx.grad(lambda x: mx.logaddexp(0.0, x))(x)
# main: 0.0; expected: approximately 2.06115e-9

The error also affects JVPs, second derivatives, and binary cross-entropy with
negative logits and zero targets. It appears at smaller input gaps in float16
and bfloat16.

For a single differentiated input, compute its sigmoid weight directly. When
both inputs need derivatives, an internal two-output primitive computes both
weights using one exponential. With e = exp(-abs(a - b)), the smaller weight is
e / (1 + e), rather than a subtraction from one. CPU, Metal, and CUDA backends
share this rule. The CPU implementation vectorizes contiguous inputs; the
primitive also supports higher derivatives, vmap, and export/import.

Tests cover both argument positions, all supported floating types on CPU,
broadcasting and strided layouts, empty arrays, non-unit cotangents, nonfinite
inputs, higher derivatives, batching, compiled gradients, export/import, and
the loss-function regression. Small-gradient checks use zero absolute tolerance.

Validation on an Apple M3 Max, macOS 27.2.1, Python 3.12.9:

• Full Metal Python suite: 1,040 passed, 9 skipped.
• CPU-only Python suite: 911 passed, 120 skipped; 18 Metal-only DLPack tests deselected and covered by the Metal run.
• CPU C++ suite: 276 test cases and 4,434 assertions passed.
• All pre-commit hooks pass. CUDA code is included but was not built or run locally.

The benchmark below uses 262,144 elements, 10 warmups and 100 timed evaluations,
with the median result across three fresh processes. The comparison restores
main's VJP/JVP rules in the same release build, then restores the revised rules;
the added primitive is unused in the baseline.

Device / dtype	Mode	Operation	Original rule (μs)	Revised (μs)
CPU / float32	compiled	both_vjp	331	286
CPU / bfloat16	compiled	both_vjp	986	725
CPU / float32	compiled	both_jvp	752	281
GPU / float32	compiled	both_vjp	176	176
CPU / float16	eager	both_vjp	188	214
GPU / float32	compiled	both_jvp	158	175

Some cases remain slower: eager float16 CPU VJP/JVP by about 14%, and compiled
Metal JVP by about 8–11% in this benchmark. The script also measures 4,096 elements,
single-input gradients, and float16/bfloat16/float32 on both devices.

Reproduce with python benchmarks/python/logaddexp_grad_bench.py --device cpu
and repeat with --device gpu or --dtype float16 / --dtype bfloat16.