You've picked one of the slowest possible ways to check
c*c == a*a + b*b // assuming c is non-negative
That compiles to three integer multiplications (one of which can be hoisted out of the loop). Even without pow()
, you're still converting to double
and taking a square root, which is terrible for throughput. (And also latency, but branch prediction + speculative execution on modern CPUs means that latency isn't a factor here).
Intel Haswell's SQRTSD instruction has a throughput of one per 8-14 cycles (source: Agner Fog's instruction tables), so even if your sqrt()
version keeps the FP sqrt execution unit saturated, it's still about 4 times slower than what I got gcc to emit (below).
You can also optimize the loop condition to break out of the loop when the b < c
part of the condition becomes false, so the compiler only has to do one version of that check.
void foo_optimized()
{
for (int a = 1; a <= SUMTOTAL; a++) {
for (int b = a+1; b < SUMTOTAL-a-b; b++) {
// int c = SUMTOTAL-(a+b); // gcc won't always transform signed-integer math, so this prevents hoisting (SUMTOTAL-a) :(
int c = (SUMTOTAL-a) - b;
// if (b >= c) break; // just changed the loop condition instead
// the compiler can hoist a*a out of the loop for us
if (/* b < c && */ c*c == a*a + b*b) {
// Just print a newline. std::endl also flushes, which bloats the asm
std::cout << "a: " << a << " b: " << b << " c: "<< c << '\n';
std::cout << a * b * c << '\n';
}
}
}
}
This compiles (with gcc6.2 -O3 -mtune=haswell
) to code with this inner loop. See the full code on the Godbolt compiler explorer.
# a*a is hoisted out of the loop. It's in r15d
.L6:
add ebp, 1 # b++
sub ebx, 1 # c--
add r12d, r14d # ivtmp.36, ivtmp.43 # not sure what this is or why it's in the loop, would have to look again at the asm outside
cmp ebp, ebx # b, _39
jg .L13 ## This is the loop-exit branch, not-taken until the end
## .L13 is the rest of the outer loop.
## It sets up for the next entry to this inner loop.
.L8:
mov eax, ebp # multiply a copy of the counters
mov edx, ebx
imul eax, ebp # b*b
imul edx, ebx # c*c
add eax, r15d # a*a + b*b
cmp edx, eax # tmp137, tmp139
jne .L6
## Fall-through into the cout print code when we find a match
## extremely rare, so should predict near-perfectly
On Intel Haswell, all these instructions are 1 uop each. (And the cmp/jcc pairs macro-fuse into compare-and-branch uops.) So that's 10 fused-domain uops, which can issue at one iteration per 2.5 cycles.
Haswell runs imul r32, r32
with a throughput of one iteration per clock, so the two multiplies inside the inner loop aren't saturating port 1 at two multiplies per 2.5c. This leaves room to soak up the inevitable resource conflicts from ADD and SUB stealing port 1.
We're not even close to any other execution-port bottlenecks, so the front-end bottleneck is the only issue, and this should run at one iteration per 2.5 cycles on Intel Haswell and later.
Loop-unrolling could help here to reduce the number of uops per check. e.g. use lea ecx, [rbx+1]
to compute b+1 for the next iteration, so we can imul ebx, ebx
without using a MOV to make it non-destructive.
A strength-reduction is also possible: Given b*b
we could try to compute (b-1) * (b-1)
without an IMUL. (b-1) * (b-1) = b*b - 2*b + 1
, so maybe we can do an lea ecx, [rbx*2 - 1]
and then subtract that from b*b
. (There are no addressing-modes that subtract instead of add. Hmm, maybe we could keep -b
in a register, and count up towards zero, so we could use lea ecx, [rcx + rbx*2 - 1]
to update b*b
in ECX, given -b
in EBX).
Unless you actually bottleneck on IMUL throughput, this might end up taking more uops and not be a win. It might be fun to see how well a compiler would do with this strength-reduction in the C++ source.
You could probably also vectorize this with SSE or AVX, checking 4 or 8 consecutive b
values in parallel. Since hits are really rare, you just check if any of the 8 had a hit and then sort out which one it was in the rare case that there was a match.
See also the x86 tag wiki for more optimization stuff.
a*a
is probably 1 instruction.pow
is at least a function call, plus whatever work the function does.pow
is a floating point function. Second, posting how much time a function takes only makes sense if you're running an optimized build. If you're running an unoptimized of "debug" build, the time is meaningless. And last, but not least, don't use pow if the exponent is an integerc*c = a*a + b*b
: multiplication, especially integer multiplication, is faster than square root. But it's only correct ifc*c
doesn't overflow.c*c
overflows, thena*a + b*b
would also overflow (assuming that they are in fact equal), so it probably should not matter much.