Core: Compute subnormal integer powers exactly on ARMv7 NEON

This makes integer powers of float and complex<float> return IEEE subnormal results on ARMv7 NEON, which fixes the unary_pow_1 and unary_pow_3 failures on the 32-bit ARM jobs (for example job 16765449075, in !3181 (merged)'s pipeline). ARMv7 NEON flushes subnormal operands and results to zero, so the double-word repeated squaring from !3149 (merged) returned pow(1.10466f, -1000) = 0 instead of 5.9e-44, returned a subnormal base unchanged, and dropped subnormal components of complex powers. Other targets compile to the same code as before: the new dispatch is compiled only under #if EIGEN_ARCH_ARM.

What changes

Everything is gated on flushes_subnormals<Packet> (32-bit ARM, a non-scalar packet, real scalar float) through tag dispatch, and both dispatch points sit under #if EIGEN_ARCH_ARM, so the new bodies are instantiated only on 32-bit ARM:

  1. scale_result rebuilds a result below the smallest normal from integer bits. With t = |x|\,2^{\,e - e_{\min} + p} (e_min = min_exponent, p = digits), the result is k\,2^{\,e_{\min} - p} for k = \operatorname{rint}(t) \le 2^{p-1}, and the integer k is the bit pattern of its magnitude: one rounding, as pldexp rounds. Only those lanes are scaled, so the other lanes of the packet raise neither FE_OVERFLOW nor FE_INVALID.
  2. A packet holding a subnormal base, found by testing its bits because a NEON comparison reads it as zero, is computed lane by lane on scalars, whose VFP arithmetic does not flush.
  3. So is a real x^-1 or x^2 whose result falls below the smallest normal: the scalar path rounds those once, as 1/x and x * x do, where the double word would round hi a second time.
  4. The unscaled loop is skipped. Its budget keeps the power normal, but not all of its residuals, nor the norm and components of a complex reciprocal: \big(1.5\cdot 2^{62}(1+i)\big)^{-1} came back as (0, -0).
  5. The pdiv shortcut for x^-1 is skipped, as NEON's pdiv is vrecpe plus two Newton steps rather than a correctly rounded division. The pmul shortcut for x^2 stays.

The cost falls on ARMv7 only: the always-scaled loop, a bit test per packet, and a scalar recomputation of the rare packets in items 2 and 3. It is not benchmarked, since QEMU timings say nothing about hardware.

Tests

The existing unary_pow test already fails without the fix. The new checks are sized from numeric_limits, so double covers the same regimes:

  • complex powers with both components subnormal, exactly: 2^{5s}(1+i)^{\pm 5};
  • complex reciprocals whose norm exceeds 1/min, or whose smaller component is subnormal;
  • x.pow(-1) == 1/x and x.pow(2) == x * x, bitwise, where the results are subnormal;
  • a packet mixing a subnormal or zero result with normal ones, (\tfrac12, 1, \dots, 1)^n and (1+i, 1+2i, \dots)^2, raises neither FE_INVALID nor FE_OVERFLOW. This is checked where packets flush; elsewhere the double word can raise FE_INVALID on the zero padding of a partial packet, as on AVX2.

They compare coefficient by coefficient: a vectorized == on NEON flushes the expected subnormal as well, and an (a == b).all() form of the first check passed with the bug in place.

Left undone

  • For |n| \ge 3 a subnormal result is rounded from hi alone, on every platform and on the scalar path, so it can be rounded twice at a tie. The comment above int_pow_double_word says correctly rounded, and unary_pow allows 3 ulp there. Breaking the tie with lo in scale_result would fix it everywhere; that is a separate change.
  • Non-integer exponents, and integer exponents beyond max_squaring_exponent, go through generic_pow, which still flushes on ARMv7.
  • On ARMv7, pldexp itself returns zero for some normal results whose intermediate product flushes, such as pldexp(2^-124f, -1), for every caller. scale_result leaves such lanes as pldexp returns them, since it rebuilds only t < 2^{p-1}; no pow input reaching them has been found.
  • unary_pow_6 (bfloat16) fails on ARMv7 for seed 2904837800 at master as well: pow(-55040bf16, -8) returns 0 instead of 1.18e-38, in int_pow_plain, which this MR does not touch.
Appendix: validation

Host: 13th Gen Intel(R) Core(TM) i7-13700HX under WSL2. ARM32 runs under qemu-arm-static 8.2.2 with the CI flags -march=armv7-a -mfpu=neon-vfpv4.

Configuration unary_pow_1..4 Random seeds (r2)
arm-linux-gnueabihf-g++-10 10.5.0 pass 100 runs, 0 failures
clang++-14 14.0.6, --target=arm-linux-gnueabihf pass 100 runs, 0 failures
x86-64 g++ 13.3.0, default (SSE2) pass —
x86-64 g++ 13.3.0, -mavx -mno-avx2 pass —
x86-64 g++ 13.3.0, -mavx2 -mfma pass —
  • Packet lanes against the scalar tail (ARM32, g++-10, arrays of one repeated base, compared bitwise):
    • real: 20000 random bases for each n \in \{-1, 2, \pm 3, 5, \pm 101, \pm 1000\}, no lane differs;
    • complex: 20000 random bases for each n \in \{-1, 2, \pm 3, 5, \pm 8, \pm 101, \pm 1000\}, no lane differs except for n = -3 (37 cases). There the scalar tail takes the unscaled loop, and the two differ by one ulp in a component about 2^{-36} below the magnitude. A long double reference found the lane correctly rounded in 2 of 3 sampled cases and the tail in 1.
  • Mutations: removing any one of items 1–5 fails unary_pow_1 or unary_pow_3.
  • Other targets' code: bench_unary_pow compiles to the same code at the merge base and at this MR, compared function by function with addresses normalized: GCC 13 and clang 18 on x86-64 (SSE2, SSE4.2, AVX, AVX2, AVX-512; C++14 and C++17) and aarch64, and GCC 14 on ppc64le (VSX), loongarch64 (LSX), and s390x and riscv64 without their vector units, whose vector backends fail to compile this benchmark without this MR as well. With tag dispatch alone, GCC's code differed in 3 to 35 functions per target (clang's did not). In one clean round of two (a CI job disturbed the other), a Zen 5 (AMD Ryzen AI MAX+ 395, WSL2, GCC 13, -march=native) measured float x^2 and x^3 5.7–7.6% slower; on the i7-13700HX with AVX2 no benchmark family was slower beyond run-to-run noise. Hence the #if.
  • Lint: clang-format-17, scripts/check_style.py and ci/scripts/run-clang-tidy.sh (clang-tidy 18) are clean.
Edited by Rasmus Munk Larsen

Merge request reports

Loading
Loading