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:
scale_resultrebuilds a result below the smallest normal from integer bits. Witht = |x|\,2^{\,e - e_{\min} + p}(e_min = min_exponent,p = digits), the result isk\,2^{\,e_{\min} - p}fork = \operatorname{rint}(t) \le 2^{p-1}, and the integerkis the bit pattern of its magnitude: one rounding, aspldexprounds. Only those lanes are scaled, so the other lanes of the packet raise neitherFE_OVERFLOWnorFE_INVALID.- 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.
- So is a real
x^-1orx^2whose result falls below the smallest normal: the scalar path rounds those once, as1/xandx * xdo, where the double word would roundhia second time. - 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). - The
pdivshortcut forx^-1is skipped, as NEON'spdivisvrecpeplus two Newton steps rather than a correctly rounded division. Thepmulshortcut forx^2stays.
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/xandx.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)^nand(1+i, 1+2i, \dots)^2, raises neitherFE_INVALIDnorFE_OVERFLOW. This is checked where packets flush; elsewhere the double word can raiseFE_INVALIDon 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 3a subnormal result is rounded fromhialone, on every platform and on the scalar path, so it can be rounded twice at a tie. The comment aboveint_pow_double_wordsays correctly rounded, andunary_powallows 3 ulp there. Breaking the tie withloinscale_resultwould fix it everywhere; that is a separate change. - Non-integer exponents, and integer exponents beyond
max_squaring_exponent, go throughgeneric_pow, which still flushes on ARMv7. - On ARMv7,
pldexpitself returns zero for some normal results whose intermediate product flushes, such aspldexp(2^-124f, -1), for every caller.scale_resultleaves such lanes aspldexpreturns them, since it rebuilds onlyt < 2^{p-1}; nopowinput 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 of1.18e-38, inint_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 forn = -3(37 cases). There the scalar tail takes the unscaled loop, and the two differ by one ulp in a component about2^{-36}below the magnitude. Along doublereference found the lane correctly rounded in 2 of 3 sampled cases and the tail in 1.
- real: 20000 random bases for each
- Mutations: removing any one of items 1–5 fails
unary_pow_1orunary_pow_3. - Other targets' code:
bench_unary_powcompiles 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) measuredfloatx^2andx^35.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.pyandci/scripts/run-clang-tidy.sh(clang-tidy 18) are clean.