Core: fix pldexp intermediate overflow for |e| >= 2^ExponentBits

Rasmus Munk Larsenrequested to merge
rmlarsen1/eigen-new:fix-pldexp-split-overflow into master

What

pldexp(0, 256) returns NaN instead of 0 in single precision, and denormals scale to inf instead of their exact finite result, for exponent magnitudes >= 2^ExponentBits (256 for float/bfloat16, 2048 for double). Regression from !2527 (merged).

Why

The depth-3 multiply tree introduced in !2527 (merged) forms c1*c1 = 2^(2b) as a standalone factor of the 4-way scale split. That factor is independent of a, so for e >= 2^ExponentBits it overflows to inf on its own even when the exact result is finite or zero:

  • pldexp(0.f, 256) → 0 * inf = NaN (should be 0)
  • pldexp(2^-149f, 256) → 2^-85 * inf = inf (should be 2^107)

The pre-!2527 multiply order (((a*c1)*c1)*c1)*c2 is safe because every partial product contains a: the intermediates grow (or shrink) monotonically towards the final result and cannot spuriously over- or underflow.

The packetmath test missed this because its hardest case (-min_exponent + max_exponent) probes e = 2045 for double and e = 253 for float — just below the broken region at 2048/256.

How

  • Restore the sequential multiply order in pldexp_generic and in the SSE/AVX/AVX512 double specializations. The other !2527 (merged) improvements (single-instruction vpmovsxdq factor construction, pnegate removal, b_remainder computation) are kept; the cost of the fix is one extra dependent multiply in the full-range pldexp only — pldexp_fast, which pexp/pexp2 use on their hot path, is untouched.
  • Document in pldexp_generic why the multiplies must not be re-associated.
  • Extend the packetmath pldexp coverage to exponent magnitudes >= 2 * max_exponent with zero, denormal, and unit lanes (denormal inputs are replaced by min() on 32-bit ARM, matching the existing FTZ workarounds).

Found while writing tests for a follow-up MR that exposes pldexp as ArrayBase::ldexp() (as promised in the !2694 (merged) review discussion); that MR will be rebased on this fix.

Not released: the regression is on master only, the 5.0 branch predates !2527 (merged).

🤖 Generated with Claude Code

Merge request reports