Beating std::pow: Up to 7x Faster Exponentiation in C++
This started, as these things often do, with a profiler output I didn't expect.
I was working on a Method of Moments solver at CMAP (Ecole Polytechnique) — the kind of code that assembles dense interaction matrices from integral equations. If you've worked with MoM before, you know the pattern: the kernel evaluations involve Green's functions, and those are full of power computations. Millions of base^exp calls per matrix assembly, most of them with integer exponents.
When I profiled the hot loop, std::pow was right there at the top. Not the matrix algebra, not the integration quadrature — the power function. That felt wrong, so I started digging.
What std::pow Actually Does
std::pow(x, y) handles the fully general case: arbitrary floating-point base and exponent. Under the hood, it computes exp(y * log(x)) — two transcendental function evaluations. That's the mathematically correct approach, and it gives you full precision for any real-valued exponent.
But when your exponent is 3, or 7, or 12... you're paying for log and exp just to do a few multiplications. On my test machine (Intel 14-core, 5.2 GHz), that costs about 213-288 ns per call depending on the type. For a single pow. In an inner loop that runs millions of times.
I wanted to see how far I could push this.
Going Recursive: The Hierarchical Approach
The idea is simple: if n is even, x^n = (x^2)^(n/2). If odd, x^n = x * (x^2)^((n-1)/2). Recurse until n is 0 or 1. That's O(log n) multiplications instead of whatever std::pow is doing internally.
template <typename BaseType, typename ExpType> inline BaseType pow_hierarchical(BaseType base, ExpType exp) { if (exp == 0) return static_cast<BaseType>(1); if (exp == 1) return base; BaseType half = pow_hierarchical( static_cast<BaseType>(base * base), static_cast<ExpType>(exp >> 1) ); return (exp & 1u) ? base * half : half; }
Nothing fancy algorithmically — this is textbook. What's interesting is what happens at the CPU level. Let's look at the actual assembly GCC generates with -Ofast -march=native:
; test_hierarchical (uint64, uint64) — GCC 13, -Ofast ; First levels are fully inlined — no call instruction imul r8, rdi ; base * base (squaring) shr r9 ; exp >> 1 cmp r9, 1 je .L79 ; base case reached imul r10, r8 ; (base*base) * (base*base) — next level shr r11, 2 ; exp >> 2 cmp r11, 1 je .L80 ; ... continues unrolling up to 9 levels deep
Looking at the generated assembly (g++ -S -Ofast), GCC unrolls the recursion up to 9 levels before falling back to an actual call — I counted 9 copies of the imul+shr+cmp pattern in the output. And the conditional "multiply by base if odd" at each level:
and r9d, 1 ; test if exp bit is set je .L82 ; skip if even imul r8, r10 ; multiply result by base
No cmov here actually — I initially expected conditional moves, but GCC prefers branch-and-jump. On a modern CPU with good branch prediction, this is fine because the branch pattern is highly predictable (exponent bits don't change between calls). The key is that the entire thing is straight-line integer code: imul, shr, and, je. No call pow@PLT, no FPU.
Compare with what std::pow generates:
test_std_pow: jmp pow@PLT ; that's it — tail call into libm
One single jmp into the math library, which then does the full exp(y*log(x)) dance with x87/SSE transcendentals. The 200+ ns difference is the cost of that generality.
Result: ~38 ns. That's a 5.6 to 7.7x improvement depending on the integer type, and the code is 6 lines of C++ that the compiler turns into about 30 instructions of pure integer arithmetic.
What Else I Tried
I didn't stop at one implementation. Part of the fun was seeing which approaches the compiler likes and which ones it doesn't. Spoiler: none of them beat hierarchical.
Binary Exponentiation (Iterative)
The classic square-and-multiply loop — same algorithm, iterative form:
template <typename BaseType, typename ExpType> inline BaseType pow_binary(BaseType base, ExpType exp) { BaseType result = 1; while (exp > 0) { if (exp & 1) result *= base; base *= base; exp >>= 1; } return result; }
~39-55 ns. The loop overhead costs about 15 ns compared to the recursive version. I expected these to be equivalent, but the assembly tells the story:
; test_binary — the inner loop .L90: test sil, 1 ; if (exp & 1) je .L89 imul rax, rdi ; result *= base .L89: imul rdi, rdi ; base *= base shr rsi ; exp >>= 1 jne .L90 ; while (exp > 0)
Clean and tight — 5 instructions per iteration. But it's a loop, so the CPU pays for the branch back to .L90 on every iteration. The hierarchical version avoids this entirely: GCC unrolls the recursion into straight-line code, so the instruction pointer only moves forward. On a modern out-of-order CPU, that's a meaningful difference.
Hardcoded Switch (ultra_fast)
For common small exponents, just spell it out:
switch (exp) { case 2: return base * base; case 3: return base * base * base; case 4: { auto sq = base * base; return sq * sq; } case 8: { auto sq = base*base; auto q = sq*sq; return q*q; } default: // fall back to binary }
~43-53 ns. I thought this would win for small exponents. The assembly is optimal for known cases — case 4 is literally two imul instructions:
.L104: imul rdi, rdi ; sq = base * base .L116: mov rax, rdi ; rax = sq imul rax, rdi ; return sq * sq ret
But the switch dispatch adds overhead — a jump table lookup via movsx + indirect jmp:
lea rdx, .L103[rip] ; jump table base movsx rax, DWORD PTR [rdx+rsi*4] ; table[exp] add rax, rdx notrack jmp rax ; indirect jump
That indirect jump is expensive (branch predictor has to guess the target) and it only fires for exp <= 8. For larger exponents, it falls through to the same binary loop anyway. The hierarchical version just does straight imul+shr all the way — no jump table, no indirect branches. Not worth the complexity.
The Numbers
All compiled with -Ofast -march=native, measured with Google Benchmark on an Intel 14-core (5.2 GHz):
| Base^Exponent | std::pow | hierarchical | binary | ultra_fast |
|---|---|---|---|---|
| UInt16^UInt16 | 214 ns | 38 ns | 39 ns | 43 ns |
| UInt32^UInt32 | 213 ns | 38 ns | 55 ns | 52 ns |
| UInt64^UInt64 | 288 ns | 38 ns | 46 ns | 53 ns |
One thing I found reassuring: GCC and Clang produce nearly identical timings on the custom kernels. The performance doesn't depend on compiler-specific magic.
Does It Break Anything?
The whole point of the MoM solver is to compute physically meaningful results. A fast pow that introduces numerical drift would be worse than useless — the matrix conditioning would suffer silently.
So I measured error systematically, using std::pow as reference:
| Method | MaxAbsErr | MaxRelErr |
|---|---|---|
std::pow | 0 | 0 (reference) |
pow_hierarchical (int base) | 4.18e-12 | 1.0e-15 (~1 ULP) |
pow_hierarchical (float base) | 1.83e-7 | 7.89e-10 |
For integer bases with integer exponents, the hierarchical version is exact — it's just multiplications, no rounding involved. For float bases, errors stay within 1-2 ULP. The solver didn't notice the difference.
Going Further: SIMD Batch Exponentiation
The MoM solver doesn't compute one pow at a time — it computes thousands with the same exponent but different bases. That's a perfect SIMD scenario: every lane does the same square-and-multiply sequence, only the data differs.
#include <immintrin.h> inline __m256d pow_avx2_d(__m256d bases, unsigned int n) { if (n == 0) return _mm256_set1_pd(1.0); if (n == 1) return bases; __m256d result = _mm256_set1_pd(1.0); __m256d current = bases; while (n > 0) { if (n & 1) result = _mm256_mul_pd(result, current); current = _mm256_mul_pd(current, current); n >>= 1; } return result; }
The if (n & 1) branch is a scalar decision on the exponent bits — the same for all 4 lanes. No lane divergence, no masking. One _mm256_mul_pd does the work of four scalar multiplies.
I benchmarked across three base types (all with uniform integer exponent, batch of 1024 elements):
double bases (AVX2, 4 lanes):
| Exponent | std::pow | hierarchical | SIMD AVX2 | vs hierarchical | vs std::pow |
|---|---|---|---|---|---|
| 3 | 3.27 ns/elem | 0.65 ns/elem | 0.44 ns/elem | 1.5x | 7.4x |
| 7 | 6.79 ns/elem | 1.20 ns/elem | 0.58 ns/elem | 2.1x | 11.7x |
| 13 | 6.11 ns/elem | 1.68 ns/elem | 0.78 ns/elem | 2.2x | 7.8x |
float bases (AVX2, 8 lanes):
| Exponent | std::pow | hierarchical | SIMD AVX2 | vs hierarchical | vs std::pow |
|---|---|---|---|---|---|
| 3 | 2.76 ns/elem | 0.61 ns/elem | 0.28 ns/elem | 2.2x | 10.0x |
| 7 | 2.89 ns/elem | 1.08 ns/elem | 0.54 ns/elem | 2.0x | 5.3x |
| 13 | 3.09 ns/elem | 1.44 ns/elem | 0.48 ns/elem | 3.0x | 6.4x |
uint32 bases (AVX2, 8 lanes, _mm256_mullo_epi32):
| Exponent | hierarchical | SIMD AVX2 | Speedup |
|---|---|---|---|
| 3 | 0.71 ns/elem | 0.27 ns/elem | 2.6x |
| 7 | 1.42 ns/elem | 0.40 ns/elem | 3.6x |
| 13 | 1.78 ns/elem | 0.55 ns/elem | 3.2x |
The pattern is clear: more lanes = more speedup. Float gets 8 lanes in AVX2 vs 4 for double, and the speedup reflects that (2-3x vs 1.5-2.2x). Integer 32-bit also gets 8 lanes via _mm256_mullo_epi32. For uint64, there is no native AVX2 multiply instruction (_mm256_mullo_epi64 doesn't exist) — you'd need AVX-512DQ or an emulated 64-bit multiply that negates the SIMD benefit.
Accuracy is identical to scalar in all cases — same multiply sequence, just multiple lanes.
Can We Beat the Compiler? Hand-Written Assembly
Looking at the assembly GCC generates for pow_hierarchical, I noticed something: GCC never uses cmov. The odd-exponent check is always a branch:
; GCC's binary loop — branches on every iteration test sil, 1 je .skip ; branch: skip if even imul rax, rdi ; result *= base .skip: imul rdi, rdi ; base *= base shr rsi ; exp >>= 1 jne .loop
For random exponents, this branch mispredicts ~50% of the time. Each mispredict costs ~15 cycles on modern Intel. I tried multiple C++ formulations to coax GCC into emitting cmov — ternary operators, branchless arithmetic — GCC refused every time.
So I wrote the loop by hand:
; Hand-crafted: branchless with cmov .loop: mov rcx, rax ; temp = result imul rcx, rdi ; temp *= base test sil, 1 ; exp & 1? cmovnz rax, rcx ; if odd: result = temp (no branch) imul rdi, rdi ; base *= base shr rsi ; exp >>= 1 jnz .loop
One more instruction (the mov to create a temp), but zero mispredicts. The cmov is a data dependency, not a control dependency — the CPU never has to guess.
To test this properly, I benchmarked with random exponents drawn from different ranges (not a single fixed exponent like the earlier benchmarks). This is where branch misprediction actually hurts — with a fixed exponent, the branch predictor learns the pattern quickly. The times here are higher than the earlier ~38 ns because of that randomness:
| Implementation | Small exp (0-10) | Medium (10-30) | Large (30-63) |
|---|---|---|---|
| Hierarchical (GCC) | 71.9 ns | 82.0 ns | 68.5 ns |
| Hand asm (cmov) | 64.2 ns | 73.5 ns | 57.4 ns |
| Hand asm (unrolled) | 65.6 ns | 66.5 ns | 56.8 ns |
The hand-written assembly is 10-23% faster than what GCC produces. The gap widens with larger exponents because there are more iterations and more branch mispredictions avoided.
The other thing GCC does poorly: register spills. The hierarchical version pushes 6 callee-saved registers and allocates 56 bytes of stack. The algorithm only needs 4 values live at a time. Our asm version uses 4 registers and zero stack.
This is where the compiler optimization boundary lies for this pattern: GCC's heuristic assumes the branch is cheaper (because when not taken, you skip the imul). That's true when branches predict well. It's wrong for varied exponents.
But What About Double Bases?
I expected the same trick to work for doubles. It doesn't. The x86 ISA has cmov for general-purpose registers (integers), but there is no cmov for floating-point registers. The cheapest branchless FP selection requires andpd/andnpd/orpd with a constructed mask — about 5 extra instructions per iteration. That's far more expensive than the branch it replaces.
I benchmarked it — the hand-written asm for doubles is 2.2-2.8x slower than GCC's branchy code. The branch predictor handles the exp & 1 pattern well enough for FP, and the cost of avoiding it is too high.
| Base type | Hand asm vs GCC |
|---|---|
| uint64 | Asm wins 10-23% (cmov = 1 instruction) |
| double | Asm loses 2.2-2.8x (FP blend = 5 instructions) |
The lesson: know your ISA. cmov is an integer instruction. For FP branchless code, you need either SIMD masking (which is what the batch version does naturally) or you accept the branch.
What Didn't Work
Not every idea panned out. Two approaches I tried that went nowhere:
Fractional exponents (x^(2/3)). The MoM kernel also needs these for the Green's function singularity extraction. I tried exp(2/3 * log(x)), cbrt(x*x), and a 10-term binomial series. On this machine, std::pow at 77.5 ns beat all of them — cbrt at 91 ns, exp_log at 99 ns, series at 166 ns. The standard library is well-optimized for the general case here. "Roll your own" doesn't always win.
Memoization. I tested std::map, unordered_map, vector<optional>, and a static array. Cache lookup overhead ranged from 80 to 134 ns — slower than the hierarchical kernel itself at 38 ns. The function is too cheap to cache.
What to Use: A Decision Guide
After all this exploration, here is what I'd recommend depending on your scenario:
Single scalar call — pick the right algorithm
| Your scenario | Best choice | ns/pow | std::pow | Speedup |
|---|---|---|---|---|
uint64^uint64 (fixed exponent) | pow_hierarchical | 38 ns | 288 ns | 7.6x |
uint64^uint64 (random exponents) | Hand asm (cmov) | 57-64 ns | 288 ns | 4.5-5.0x |
uint32^uint32 (fixed exponent) | pow_hierarchical | 38 ns | 213 ns | 5.6x |
double^int (fixed exponent) | pow_hierarchical | ~38 ns | 327 ns | 8.6x |
double^(2/3) or fractional | std::pow | 77 ns | 77 ns | 1x (can't beat it) |
Note: with fixed/predictable exponents, the branch predictor learns the pattern and hierarchical wins. With random exponents, branches mispredict and the hand-written cmov asm gains 10-23%.
Batch with uniform exponent — use SIMD
All benchmarked over 1024 elements with exponent = 7:
| Your scenario | Best choice | ns/pow | Scalar hierarchical | Speedup |
|---|---|---|---|---|
double batch | AVX2 (4 lanes) | 0.58 ns | 1.20 ns | 2.1x |
float batch | AVX2 (8 lanes) | 0.54 ns | 1.08 ns | 2.0x |
uint32 batch | AVX2 mullo_epi32 (8 lanes) | 0.40 ns | 1.42 ns | 3.6x |
uint64 batch | Scalar (no AVX2 mullo_epi64) | 1.20 ns | — | needs AVX-512DQ |
What NOT to do
- Don't write branchless asm for floats. No
cmovfor FP registers — the blend workaround costs 5 instructions and is 2-3x slower than a branch. - Don't memoize. Cache lookup (80-134 ns) is slower than the kernel itself (38 ns). The function is too cheap to cache.
- Don't use a switch table for small exponents. The indirect jump costs more than the hierarchical approach which the compiler already unrolls.
- Don't use
std::powfor integer exponents. Ever. It callsexp(y*log(x))— two transcendentals for what should be a few multiplies.
The full picture, for one scenario: 1024 doubles, exponent = 7
| Approach | ns/pow | vs naive |
|---|---|---|
std::pow in a loop (naive) | 6.79 ns | 1x |
pow_hierarchical in a loop | 1.20 ns | 5.7x |
| AVX2 SIMD batch | 0.58 ns | 11.7x |
For the MoM solver that started this whole thing, switching from std::pow in a scalar loop to a SIMD batch path turned a coffee break into an interactive simulation.
The code is on GitHub: powerix.
A Note on Hardware Dependence
I ran these benchmarks on two machines: an Intel 22-thread at 4.7 GHz (CMAP lab workstation) and a 14-core Intel at 5.2 GHz. The integer exponent story is consistent — hierarchical wins by 5-8x everywhere. But the fractional exponent results flipped: on the first machine, exp_log beat std::pow by 1.6x. On the second, std::pow was actually faster.
The takeaway: for integer exponents, the win is architectural — you're replacing two transcendentals with a few multiplies, and that wins on any x86 chip. For fractional exponents, the answer depends on how well your specific CPU pipelines the math library functions. Don't trust anyone's benchmarks (including mine) — run them on your target hardware.