Compute integer powers by double-word repeated squaring, real and complex

Plain repeated squaring doubles its accumulated rounding error at every squaring, so x.pow(n) with an integer exponent came out about 0.7 n ulp off: 5 ulp at n = 8, 74 ulp at n = 100, 700 ulp at n = 1000. The floating-point exponent path limited it to [-3, 7] and used generic_pow (log/exp) beyond, at 4.5 ns (float) to 19 ns (double) per element on AVX2; the integer exponent path had no limit at all. Complex bases always went through std::pow, which for an integer-valued floating exponent is exp(y log z): (1+i)^2.0 = (1.2e-16, 2), hundreds of ulp for larger n.

This MR computes integer powers of float and double bases, real and complex, by repeated squaring in double-word arithmetic, and vectorizes complex pow for the first time.

Method

• The running power is a double word {hi, lo}; each step is a fast_twoprod (Joldes, Muller & Popescu 2017, DWTimesDW1, 7 u² per step), so the result is correctly rounded unless the exact power lies within about 7 n u² of a rounding boundary. A negative exponent reciprocates the base as a double word first.
• The power is scaled by powers of two throughout: the base is normalized to [2, 4) once and the running power renormalized every fourth step through its exponent bits, with the binary exponent tracked separately and applied by one pldexp at the end. Without this, a result below 2^-100 (float) makes the residuals subnormal and the microcode assists made ordinary inputs like (1e-4f)^8 a hundred times slower than master. Nothing under- or overflows before the final pldexp, so complex bases no longer produce NaN or a spurious infinity when a lane product overflows while the result is finite.
• Complex bases run the same code on separate real and imaginary double words sharing one exponent (z² = (a² − b²) + 2ab·i is three products); a complex packet's interleaved real view is split once into a packet per component and joined at the end.
• A packet whose bases all lie within 2^(±budget/|n|), budget = −(min_exponent + digits), cannot leave the normal range at any step, so it runs the loop with no normalization, exponent, renormalization or final pldexp at all (complex: both nonzero components are tested, one less in the exponent for √2ⁿ and at most half the range for the reciprocal's |w|²). Zero, inf and NaN fail that test. For n ≤ 100 this covers |x| ∈ [1/512, 512] for double, [0.5, 2] for float.
• Integer-valued floating exponents use squaring up to 2^12 (float) or 2^20 (double), the measured crossover with generic_pow, and generic_pow beyond. Integer-typed exponents beyond the base's significand always square, since generic_pow would round them and lose the parity of pow(-0, INT_MIN). x^1, x^-1 and a real x^2 are single correctly rounded operations and short-circuit. half, bfloat16, long double and custom scalars keep plain squaring with the previous cutoffs.
• The double-word path is used only where the packet provides what it needs (exponent-bit scaling, pldexp, ordered compares, predux_any), marked by has_exponent_bit_ops: an integer_packet typedef, or for AVX Packet4d without AVX2 an SSE2 emulation on its 128-bit halves in arch/AVX. The GPU, MSA and HVX packets have none and keep plain squaring.
• A nonzero complex component s far below the larger one L has z^n = L^n(1 + n\delta) to working precision, \delta = i\,b/a or -i\,a/b: the loop runs on L alone and adds L^n n\delta with a double-word coefficient, n included, since it is inexact in float beyond 2^{24}. The cutoff is |s/L| < 2^{e_{\min} + p + 63}, not the normal range: before the first renormalization a reciprocal power falls to 2^{-62} of its scale, which below that cutoff underflows the smaller component's residuals. First order drops (nr)^2/2 \le u/4 for |n| \le 2^{25} (float) or any n (double); larger float n renormalize at every step instead. The unscaled loop also requires the smaller nonzero component in range, so an element's result does not depend on its packet neighbours.
• Real zero, infinity and NaN: NaN propagates; ±0 and ±inf become NaN in the input scaling and are substituted at the end by the base (or its reciprocal) with the sign dropped for an even exponent. Every combination of the 21 values {±0, ±inf, ±NaN, ±denorm_min, ±min, ±max, ±1, ±2, ±0.5, ±3, ε} with 277 integer exponents (all of [−64, 64], ±(2ᵏ − 1, 2ᵏ, 2ᵏ + 1) for k = 7 … 30, and INT_MAX, INT_MIN and their neighbours), as int and as float/double, matches std::pow (in long double for int exponents) including the sign of zero, finite results within 1 ulp: 11 634 cases per type on AVX2+FMA, SSE2 without FMA, and without vectorization.
• Complex special values follow the rules and table under Special values below, without calling std::pow. Zero, infinite and NaN bases are replaced by 1 + i before the loop and their results set at the end, so they neither force their packet into the scaled loop nor come back NaN. None of them is in range, so only a packet already bound for the scaled loop is checked for them. A base with one zero component gives the zero component of its power the sign of the first-order limit: z = L(1 + i t) gives z^n = L^n(1 + i n t), so the sign of i\,n\,t\,L^n.
• Includes the functor_traits fix from !3144 (merged) (no packet access for a complex exponent), since the same line changes. The double-word twosum, which DPR1EigenSolver also uses, becomes the branch-free SloppyDWPlusDW (11 operations instead of 16, no less accurate than the previous version) and gains a twodiff; fast_twoprod is added next to twoprod. Special values in generic_pow and handle_nonint_nonint_errors are selected with por/pand on the all-ones comparison mask (1492 → 1458 instructions for the two float pow overloads, AVX2).

Special values

Neither standard pins these down. C Annex G gives cpow no special-value table ("may be computed as cexp(y*clog(x))"), and C++ defines complex pow through exp and log. The value is mathematically determined in some cases; in the rest only its representation (signs of zeros, which component is infinite or NaN) is left to the implementation, and libstdc++'s follows its evaluation order.

Eigen treats zero and infinite bases as the points 0 and ∞ of the extended complex plane. As std::proj does, every value with an infinite component, a NaN one included, is the one infinity +inf + 0i. So a zero base gives +0 + 0i for n > 1 and +inf + 0i for n < 0, and an infinite base the reverse, whatever the signs of their zeros. A base with a NaN component and no infinite one gives NaN + NaN i. n = 1 returns the base unchanged. None of this calls std::pow, so the results are the same with every standard library.

complex<float>, AVX2+FMA; the packet and the scalar tail give the same results. Each cell is Eigen / libstdc++ std::pow, and bold marks where they differ.

base	n = 2	n = 3	n = −1	n = −2	defined?
(1e30, 1e30), finite, power overflows	(0, inf) / (−inf, inf)	(−inf, inf) / (−inf, nan)	(5e-31, −5e-31) / same	(0, −0) / (−0, −0)	yes: the correctly rounded power. libstdc++ overflows in a^2 - b^2
(2, −0), one zero component	(4, −0) / (4, 0)	(8, −0) / same	(0.5, 0) / same	(0.25, 0) / same	the nonzero component yes; the zero's sign no. Eigen: first-order limit
(−0, 2), one zero component	(−4, −0) / same	(0, −8) / same	(−0, −0.5) / (0, −0.5)	(−0.25, 0) / same	as above
(0, 0)	(0, 0) / same	(0, 0) / same	(inf, 0) / (inf, nan)	(inf, 0) / (inf, nan)	0 for n > 0, ∞ for n < 0; not their representation
(−0, −0)	(0, 0) / same	(0, 0) / (0, −0)	(inf, 0) / (−inf, nan)	(inf, 0) / (inf, nan)	as above
(inf, 0)	(inf, 0) / (inf, nan)	(inf, 0) / (inf, nan)	(0, 0) / same	(0, 0) / same	∞ for n > 0, 0 for n < 0; not their representation
(−inf, 2)	(inf, 0) / (inf, −inf)	(inf, 0) / (nan, inf)	(0, 0) / (−0, −0)	(0, 0) / same	as above
(0, inf)	(inf, 0) / (−inf, nan)	(inf, 0) / (nan, −inf)	(0, 0) / (0, −0)	(0, 0) / (−0, 0)	as above
(inf, inf)	(inf, 0) / (nan, inf)	(inf, 0) / (−inf, inf)	(0, 0) / (0, −0)	(0, 0) / (0, −0)	as above
(inf, nan)	(inf, 0) / (inf, nan)	(inf, 0) / (inf, nan)	(0, 0) / same	(0, 0) / same	an infinity (Annex G counts it as one); representation no
(nan, 1), (1, nan), (nan, nan)	(nan, nan) / same	(nan, nan) / same	(nan, nan) / same	(nan, nan) / same	yes: NaN

The libstdc++ column comes from its repeated multiplication, and other standard libraries are free to differ.

Accuracy and speed relative to generic_pow

AVX2+FMA (i7-13700KF, GCC 13.3, -O3), dynamic arrays of 4097 elements with bases log-uniform in |xⁿ| ≤ 2⁶⁰, best of 25 runs. ulp is the worst over 20 000 log-uniform bases of either sign whose powers are normal and at least 2^digits above the smallest normal, against __float128 powq, identical for int and floating exponents. generic_pow is 4.4 ns/elem for float and 18.3 for double.

n	float ns	float speedup	float ulp	double ns	double speedup	double ulp
2	0.12	36.8×	0.50	0.15	121.8×	0.50
3	0.36	12.3×	0.50	0.62	29.5×	0.50
4	0.39	11.3×	0.50	0.68	26.9×	0.50
5	0.47	9.4×	0.50	0.85	21.5×	0.50
7	0.57	7.8×	0.50	1.06	17.2×	0.50
8	0.50	8.8×	0.50	0.91	20.1×	0.50
12	0.62	7.1×	0.50	1.18	15.5×	0.50
16	0.66	6.7×	0.50	1.20	15.2×	0.50
32	0.80	5.5×	0.50	1.55	11.8×	0.50
64	0.98	4.5×	0.50	1.86	9.8×	0.50
100	1.37	3.2×	0.50	2.64	6.9×	0.50
256	3.28	1.3×	0.50	2.71	6.7×	0.50
1000	5.50	0.8×	0.50	11.11	1.6×	0.50
4096	4.83	0.9×	0.50	9.83	1.9×	0.50
65536	generic		1.9	14.06	1.3×	0.50
2^20	generic		1.9	19.06	1.0×	0.50
−2	0.41	10.8×	0.50	0.86	21.2×	0.50
−3	0.51	8.7×	0.50	0.98	18.6×	0.50
−4	0.55	8.0×	0.50	1.02	17.9×	0.50
−8	0.71	6.2×	0.50	1.38	13.2×	0.50
−16	0.89	5.0×	0.50	1.79	10.2×	0.50
−100	1.88	2.4×	0.50	3.69	5.0×	0.50
−1000	6.62	0.7×	0.50	13.28	1.4×	0.50

Cost tracks the bit length plus the popcount of n. Results within 2^digits above the smallest normal stay at 0.50 ulp; results that are themselves subnormal are rounded a second time by the final pldexp and reach 0.75 ulp. Float exponents beyond 2^12 go to generic_pow.

Complex, integer exponent, vectorized: at most 0.70 ulp of |zⁿ| for n = ±1 … ±1000 (4000 bases per n, random angle, log-uniform magnitude keeping zⁿ in range, against a __float128 product), against 0.7–870 ulp for libstdc++'s std::pow(complex, int) and 58–4800 for the exp/log path. The bound is normwise, relative to |zⁿ|, as a^2 - b^2 and ac - bd can cancel, and its n u^2 term grows with |n|: (0.6 + 0.8i)^n in complex is 0.17 u off at n = 4096, 0.53 u at 2^{24} and 4.2 u at 2^{31} - 1. Speed at |z| = 1.05: complex n = 2 / 3 / 8 / 100 at 1.6 / 3.4 / 4.2 / 17 ns per element against 1.2 / 1.2 / 2.4 / 5.4 for libstdc++'s repeated squaring, which master uses, and 54–57 for the exp/log path; complex 3.4 / 7.0 / 8.6 / 28 against 1.3 / 1.4 / 2.7 / 6.7 and 93–102. A complex base with a floating exponent type stays on the scalar path (no vectorized complex pow for the non-integer case): 9–74 ns instead of 54–105.

What gets slower

Everything that previously used plain squaring with a small integer exponent: x.pow(3) on a dynamic float array goes from 0.08 to 0.38 ns per element (double 0.16 to 0.64), and on a fixed-size array with a literal exponent, where master folds to two multiplies, from 0.02 to 0.16 (float[8]) and 0.05 to 0.32 (double[4]). That is the price of 0.5 ulp against 1.3 (n = 3) to 4 (n = 7) ulp; with use_repeated_squaring it is a one-line change to keep plain squaring for the smallest exponents if that is preferred. Complex integer powers at |z| ≈ 1 are 1.3–5.2× slower than master's scalar libstdc++ squaring for positive n (above), for 0.7 ulp against up to 870. A complex base whose smaller nonzero component lies below 2^(−budget/|n|) takes the scaled loop: +28% for random unit-square complex at n = 17, about 5% at small n. Real-valued complex data, (a, ±0), runs the zero-sign fixup on every packet: +16–38% at small n, +5% at n = 100.

Validation

test/unary_pow.cpp (new): float and double against pow in the next wider format, complex against a double-word reference, exceptional values, signed zeros, subnormal bases and results, exponents on both sides of every cutoff up to INT_MAX, half and bfloat16 against float, complex at 2^±13107, widely separated complex components per component (subnormal and 2^(min_exponent+4) smaller components up to n = ±1003, float n = 2³⁰ + 1 against the binomial series), independence from packet neighbours, complex special values against std::pow exactly (scalar tail included) and zero signs against the perturbed-base limit, and exact (1+i)^(8k+7) at the top of the range. tensor_gpu covers integer exponents on aligned GPU maps (run on an RTX 4070 Ti, CUDA 13.3). Passed with array_cwise, packetmath, mixingtypes, array_for_matrix and numext on AVX2+FMA (with and without F16C), AVX, SSE2 without FMA (Veltkamp two-product), and without vectorization; clang 18 builds and passes AVX2 and scalar; AVX-512 compiles without warnings but could not be run here. Where long double is binary64 (MSVC, Apple arm64) the complex reference cannot form negative powers, so those are skipped; checked with -mlong-double-64 on SSE2 and AVX2+FMA. benchmarks/Core/bench_unary_pow.cpp compares against generic_pow and std::pow. MSA's complex packets gain the as_real typedef the complex path needs; that architecture was not built.

Update (d9d55f8f): zero complex bases no longer go through std::pow, reversing part of the fix for @rmlarsen1's P2 on special values. Their result is implementation-defined there, and routing them through it was expensive: a zero forced its whole packet into the scaled loop, and each zero lane paid a store, a scalar std::pow and a reload. They now give +0 + 0i (n ≥ 2) and +inf + 0i (n < 0), documented on ArrayBase::pow and tabulated under Special values. Infinite and NaN bases still match std::pow. AVX2+FMA, 4096 elements, ns per element:

	zeros	before	after
complex, n = 3	0% / 10% / 50% / 100%	3.11 / 4.30 / 8.10 / 6.49	3.08 / 3.37 / 3.84 / 3.69
complex, n = −2	0% / 10% / 50% / 100%	3.50 / 5.12 / 11.53 / 10.54	3.50 / 4.23 / 4.12 / 4.02
complex, n = 3	0% / 10% / 50% / 100%	6.16 / 7.61 / 12.32 / 12.22	6.25 / 6.89 / 8.03 / 7.62
complex, n = −2	0% / 10% / 50% / 100%	7.23 / 9.22 / 16.70 / 17.83	7.46 / 8.20 / 8.95 / 8.31

Update (a6048978): infinite and NaN complex bases no longer take std::pow either. Their results depended on the standard library, and libstdc++'s aren't even directionally meaningful: (-\infty + 2i)^3 = nan + inf i. They now follow the extended-plane rules under Special values. Nothing in the double-word path calls std::pow any more, and the per-lane store, scalar call and reload are gone. Data without special values and always-scaled data (n = 1000) are unchanged, within ±2%. Data with 50% zero bases is 10–13% slower than with the zero-only check of d9d55f8f, but still 2–2.7× faster than before it.

🤖 Generated with Claude Code