# SSE loop to walk likely primes

This is a continuation of a discussion that was started here. While there are some interesting points there about instruction timing and latency, it is not necessary to read that Question to understand this one.

### What I'm looking for in a review

My project is heavily bottlenecked on my (6 asm instruction) inner loop. I'm hoping someone can suggest ways to improve this code.

I have trimmed out as much of the surrounding code as I could, to focus on the part I'm hoping to get help optimizing.

### Specs

To run this code you must have a processor that supports AVX2 instructions. My target is Kaby Lake, but I'm told that Skylake has identical instruction timing, so that would work too. It is an x64 program written using VS 2017, running on Windows 7 (although there isn't anything OS-specific here, so it should run on linux with minor tweaks).

Since the code of interest is written in asm, it doesn't much matter what build options you use (including debug).

### What this code does

The purpose of this code is to walk 'likely' primes. When walking consecutive numbers to see if they are prime(-ish), there are tricks you can use.

Starting simply, you know that every 3rd number is going to be evenly divisible by 3 (and thus not prime). Likewise for 5, 7, 11, etc.

This algorithm identifies likely primes by eliminating all numbers which are a multiple of the first 30 primes (from 3 through 127, as 2 is handled implicitly). The basic idea is to maintain periodic 30 counters, one for each prime, each of which has a period equal to the associated prime. For each candidate value, we decrement each counter by 1, and check if any counter is zero, indicating the candidate is a multiple of that prime. If no counter is zero, we have a likely prime.

Using SSE, you can put all 30 counter values in a single ymm register and do the decrement mod P on all counters simultaneous, wrapping using vpminub.

ymm2: 2 4 6 10 12...


And then perform math on them all:

ymm0: 1 1 1 1 1

vpsubb ymm2, ymm2, ymm0        ; Decrement all the bytes by 1


Since we are counting by ones, an underflow will set the byte to 0xff, which can be detected and tested:

vpmovmskb r10d, ymm2

test r10d, r10d
jnz nottop                     ; No underflows means prime-ish, otherwise keep looking


When the underflow occurs, you need to reset the values. Since the underflow results in a really big number, the SSE vpminub instruction can be used:

ymm1: 2 4 6 10 12...

vpminub ymm2, ymm2, ymm1       ; Adjust any underflows


Given the size of a ymm register, you can check ~30 factors with just a couple of instructions and 3 ymm registers.

A few items of note:

• My code doesn't start checking at 1, so it needs to calculate the initial value of ymm2. This is only done once (per thread), so it is not a bottleneck. That's why I have omitted it and provided a few hard-coded values (see the C code).

• While a ymm register could hold 32 primes, I'm only using 30 (plus the number 2 which is handled implicitly). That's because the next two prime numbers (131 and 137) require 8 bits, which looks like an underflow and mucks up the algorithm.

• My target platform is Kaby Lake, which has 4 hyperthreaded cores. Running 8 threads doesn't give me 8 times the performance of a single thread, but it does give better total throughput than using 4,5,6,7 or 9. Since CPU execution ports are shared between the threads running on a core, using 2 threads per core helps keep the ports fully occupied (see Peter's detailed comments on this here). While the C code below doesn't use threading, it's structured to support it.

• I have omitted the actual work I do for each likely prime from this post, since it would make the code way more difficult to understand (as well as requiring libraries, more hw, etc). Maybe next time.

• (For obvious reasons) the asm code doesn't print out any status messages while running. So it's just going to sit there blinking a cursor at you for a minute or so (YMMV):

    Count for thread 0: 7945632366
Elapsed: 66
Elapsed: 65
Elapsed: 0

• One of the possible optimizations mentioned on the original thread was based on how many times the innermost loop is expected to run before exiting to the outer loop. While the loop can run up to ~100 times, it is going to be weighted heavily toward 1 rather than being evenly distributed. I wasn't able to make that suggestion work, but I did produce some typical stats if someone else has some ideas:
    loop    cumulative
count   percentages
1          15%
2          30%
3          53%
4          61%
5          71%
6          81%
7          86%
8          89%
9          93%
10         95%


### The C code

As mentioned, this is just a (drastically simplified) test harness to drive the asm (which is the code I care about).

#include <stdio.h>
#include <time.h>

typedef unsigned __int64 ULONGLONG;
typedef unsigned char BYTE;

extern "C" ULONGLONG ProcessA(ULONGLONG start, ULONGLONG end, const BYTE *state);

extern "C"
{
ULONGLONG b;

{
case 0:
{
__declspec(align(32))
BYTE state[32] = {
0x00, 0x03, 0x00, 0x00, 0x07, 0x08, 0x0c, 0x0f,
0x0d, 0x0d, 0x1c, 0x00, 0x10, 0x06, 0x20, 0x24,
0x29, 0x18, 0x3d, 0x1c, 0x10, 0x04, 0x37, 0x53,
0x3d, 0x40, 0x45, 0x39, 0x55, 0x1e, 0x62, 0x40 };

b = ProcessA(69780348563, 131, state);
break;
}
case 1:
{
__declspec(align(32))
BYTE state[32] = {
0x00, 0x02, 0x05, 0x0a, 0x0b, 0x0e, 0x0a, 0x03,
0x08, 0x0b, 0x0d, 0x00, 0x23, 0x15, 0x1d, 0x28,
0x33, 0x14, 0x19, 0x1f, 0x0a, 0x47, 0x04, 0x38,
0x0a, 0x0f, 0x17, 0x6b, 0x34, 0x3e, 0x45, 0x87 };

b = ProcessA(139560697001, 69780348563, state);
break;
}
case 2:
{
__declspec(align(32))
BYTE state[32] = {
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x40,
0x20, 0x10, 0x08, 0x04, 0x02, 0x01, 0x00, 0x00 };

// Note: this isn't the correct start/end for this state, so
// the 'likely' primes produced in ProcessA won't actually
// be primes.  However this is still useful for testing since
// it produces the largest loop count I know how to create
// using 30 primes.
b = ProcessA(501, 305, state);
break;
}
default:
{
b = 0;
break;
}
}

return 0;
}

int main(void)
{
time_t before, after;

for (int x = 0; x < 3; x++)
{
time(&before);
RunCpu(x);
time(&after);

printf("Elapsed: %d\n", (int)(after - before));
}
}


### The asm code

This is the heart of the code. The loop with the bottleneck should be easy to spot...

CONST SEGMENT READONLY ALIGN(32) 'CONST'

; We are only using 30 primes, so the last 2 are left as zero
allones:
db 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1
db 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0

; The first 30 prime numbers (minus 1)
asmprimes:
db   2,   4,   6,  10,  12,  16,  18,  22
db  28,  30,  36,  40,  42,  46,  52,  58
db  60,  66,  70,  72,  78,  82,  88,  96
db 100, 102, 106, 108, 112, 126,   0,   0

CONST ENDS

.code

; Call with:
; extern "C" ULONGLONG ProcessA(ULONGLONG start, ULONGLONG end, const BYTE *state);

ProcessA proc

; Windows calling convention:

; rcx - start
; rdx - end
; r8 - current state
; r9, r10, r11 - scratch
; rax - return value

vmovups ymm0, YMMWORD PTR [allones]     ; Set all bytes to 1
vmovdqa ymm1, YMMWORD PTR [asmprimes]   ; Load primes
vmovdqa ymm2, YMMWORD PTR [r8]          ; Load the current state

xor eax, eax                            ; Returns # of prime-ish found

xor r9d, r9d                            ; How many steps to next prime-ish?
vpsubb ymm2, ymm2, ymm0                 ; Decrement all the bytes by 1
jmp top

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; This loop is the bottleneck.
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
align 16

nottop:

vpminub ymm2, ymm2, ymm1                ; Adjust any underflows
vpsubb ymm2, ymm2, ymm0                 ; Decrement all the bytes by 1

align 8

top:
vpmovmskb r10d, ymm2                    ; Check the upper bits for underflows

test r10d, r10d
jnz nottop                              ; No underflows means prime-ish
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; End of bottleneck loop
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

; r9 contains the count of iterations (0-97 when using 30 primes)

lea r9, [r9 + r9 + 2]                   ; Account for even numbers and starting at 0

; Decrement the current value to the next 'likely' prime.
; This may require sbb for values greater than 64 bits.
sub rcx, r9

; At this point, rcx contains the likely prime.

cmp rcx, rdx                            ; if (start < end) goto done;
jb done

; The vpsubb is done here in the (mostly vain) hope that it can be scheduled
; to run concurrently with other instructions.

xor r9d, r9d                            ; Reset for next loop
vpsubb ymm2, ymm2, ymm0                 ; Decrement all the bytes by 1

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; This is a placeholder for a more involved routine.
;
; There's no simple way to provide a runable version
; of it, since it makes use of a lookup table that
; requires libraries to create and is too big (8gig)
;
; OTOH, if the loop above can be optimized to the
; point where *this* code becomes the bottleneck, I'll
; see what I can do.
;
; FYI it makes no use of any xmm or ymm register, nor
; does it 'call' any other routines (it's just math).
; It does make use of most GP registers.

inc rax ; Count how many we find

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

jmp top

done:

; rax already contains the return value.
; All the registers used are 'scratch', so didn't get saved and
; don't need to be restored.
; All that leaves is clearing the ymm registers.
vzeroall

ret

ProcessA endp

end

• In practice, I think anyone with an AVX2 processor can take a decent shot at it. Yes, Haswell as some different timings, and there are edge cases where one algorithm would be faster on Haswell but slower on Skylake, but what works on one generally works on another - and you can confirm this by running on your Kaby Lake. Jul 21, 2017 at 16:29
• @BeeOnRope: posted an answer. I haven't tried to look for algorithmic improvements, just things that might be friendlier for branch-prediction. I suspect this spends a lot of time on branch-mispredicts given how many different iteration counts are common. Jul 21, 2017 at 22:30
• @BeeOnRope I'm pretty sure. See comments after peter's answer. Jul 22, 2017 at 0:20
• "Not too shabby" - Thanks, I do try. I've been experimenting with an alternative which is just starting to work. It trades precision for speed. Consider: Skipping 2's cuts the numbers to be checked in half. Skipping 3's would cut it by a third, except that the evens were already omitted, so it's only a sixth. 5's would be a fifth, except... By the time I get to prime #30, I'm probably not trimming much. Since the outer loop turns out to be relatively cheap, I can use a looser definition of 'likely.' And I've still got Peter's improvement for when outer (inevitably) becomes more expensive. Jul 22, 2017 at 3:07
• Indeed, you can encode more of the smaller primes in the algorithm itself, which kind of makes them free, at a code size cost in the inner loop or a table lookup, as you move beyond 2 however, the pattern is no longer "even" (formally, having a period of P with exactly 1 likely prime in P). With 2 it was just a matter of using the normal loop (and adjust the initial state) and then doubling the value (your lea [r9 + r9 ...] line), but to handle 2 and 3, you need to do something different for the P mod 6 == 1 and P mod 6 == 5 cases. Unroll by 2 or table lookup could both work. Jul 22, 2017 at 16:03

Let's review this code purely from a performance angle, without a focus on style or anything else (in addition to optimization suggestions, Peter already mentioned several things in areas other than performance).

First, you can play with all the algorithm discussed here in this github repo. I compiled it on Linux but it should approximately work on Windows if you have nasm or yasm - if you add a thunk to adjust the calling convention. If someone really wants it I'll do it.

## Profiling

If you've ever asked for performance help, no doubt someone has told you to profile, profile, then profile some more. So sure, let's start there. I'm going to use Linux's perf since it is awesome and free and available on Linux, but you can get the same information presented nicely in Windows using VTune or maybe with this stuff.

Anyway, let's run perf stat on the original ProcessA algorithm to get a feel for any high-level issues:

$perf stat ./likely-primes time ProcessA 69780348563 Finding all likely primes between 69780348563 and 131... Count 7945632366, density=0.114 elapsed 110.689s, cycles/candidate 4.12, cycles/prime 36.22 Performance counter stats for './likely-primes time ProcessA 69780348563': 110689.294210 task-clock (msec) # 1.000 CPUs utilized 227 context-switches # 0.002 K/sec 0 cpu-migrations # 0.000 K/sec 99 page-faults # 0.001 K/sec 286,208,237,077 cycles # 2.586 GHz 249,386,067,831 instructions # 0.87 insn per cycle 50,840,454,194 branches # 459.308 M/sec 9,257,714,721 branch-misses # 18.21% of all branches 110.691877920 seconds time elapsed  ### Branch Mispredicts This code is being crushed by branch mispredictions: 20% of branches are mispredicted, and there are a lot of branches about 1 of every 5 instructions1. Usually you are looking for a value less than 1% to avoid a big impact. Off the top of your head this is going to have an impact of about 9.3 billion misses * 15 cycles = ~140 billion cycles, or about half the total running time. This also explains the poor IPC of 0.87. ### Toplev.py We can try toplev.py to confirm: $ ../pmu-tools/toplev.py ./likely-primes time ProcessA 697803485
Using level 1.
...
...
This category represents slots fraction wasted due to
incorrect speculations...


So toplev.py is telling us that we're wasting more than 63% of our CPU's potential due to "bad speculation" - which in this case is 100% branch mispredicts.

Now this behavior isn't unexpected given the frequent and random nature of the exits from the inner loop. In fact, we expect about 1 misprediction per found likely-prime, since the jnz nottop is usually jumping (because only about 23% of examined values are likely prime2), so it will be predicted taken and mispredict on loop exit. We can check by comparing the mispredictions to the found primes: 9,257,714,721 misses / 7,945,632,366 primes = ~1.16 misses / prime. So even slightly worse than we expected, probably because the predictor will constantly being finding "patterns" that simply aren't due to the random nature of the likely-prime sequence.

All that said, note that we still are still processing candidate values in 4.12 cycles (which works out to ~36 cycles per found likely prime), which doesn't sound too shabby...

# Making It Fast

So the obvious question after profiling is: can we do better?

## Optimizing the Existing Algorithm

Let's start by seeing what we can do with modest optimizations to the existing algorithm.

Evidently, a first line of attack would be to get rid the branch mispredictions. This means we need to do something inside the inner loop which records without branching the likely primes. For example, since we are just counting the primes, let's just conditionally increment the sum, in our next version of the routine, ProcessA2:

.top:

vpminub ymm2, ymm2, ymm1                ; Adjust any underflows
vpsubb ymm2, ymm2, ymm0                 ; Decrement all the bytes by 1

vpmovmskb r10d, ymm2                    ; Check the upper bits for underflows

test r10d, r10d

setz r14l

dec rcx
jnz .top


The inner loop has now grown to 8 instructions from 6, but we have banished the branch mispredictions:

perfc stat ./likely-primes time ProcessA2 69780348563 Finding all likely primes between 69780348563 and 131... Count 7945632366, density=0.114 elapsed 36.285s, cycles/candidate 1.35, cycles/prime 11.87 Performance counter stats for './likely-primes time ProcessA2 69780348563': 36285.661492 task-clock (msec) # 0.999 CPUs utilized 1,338 context-switches # 0.037 K/sec 0 cpu-migrations # 0.000 K/sec 100 page-faults # 0.003 K/sec 94,117,387,177 cycles # 2.594 GHz 279,182,098,024 instructions # 2.97 insn per cycle 34,901,500,964 branches # 961.854 M/sec 72,383 branch-misses # 0.00% of all branches 36.305622051 seconds time elapsed  What a difference a removed branch makes! It's about three times faster, at 1.35 cycles per candidate. That's despite the fact that we are executing more instructions: about 280 billion versus 240 million. The upside is all due to removing branch-misses, which are now reported at 0.00%, and the IPC has increased to ~3 instructions per cycle. ## Loop Splitting Of course in the real world, you don't want to just count the primes, you want to do something with them. That's fine: it's a slight modification to the above to generate a bitmap indicating which values are likely primes, rather than simply counting them, without slowing down much. So to avoid the mis-predictions, you process some fixed number of candidates with the loop above, generating a bitmap, and then you iterate in a branch-prediction aware way over the bitmap (e.g., using tzcnt and popcnt) to generate the likely prime values for your "secondary processing". I won't actually flesh this out fully for ProcessA2 since we are about to move into the fast lane with a different approach entirely (which still ultimately uses the same "bitmap" output format. ## Bitmaps FTW Let's step back a moment and understand what the core of the algorithm is doing. Basically it implements 30 periodic counters incremented in sync and tries to determine if at least one counter has "wrapped" on every iteration. To do this, it uses 30 byte counters in a ymm. Since AVX2 lets us do 32 byte operations per operation, it means we can do 30 operations on this counter per instruction (and perhaps up to 3*30 = 90 if we use all 3 vector ports fully). We get lucky that the vpminub instruction works well to do a 30-way mod or "wrap" operation! What if instead of using byte counters, we use a series of bitmaps with one bit per candidate, which encode the same periodic behavior as the counters? That is, to replace the prime = 3 counter which goes 2 -> 1 -> 0 -> 2 -> ... we use the bitmap 0b100100100... with every 3rd bit set? Well now a 256-bit ymm register holds 256 counter values, not 30. Of course, the correspondence isn't exact, since everything is transposed: one register contains a lot of state, but for one prime. You'll need 30 such registers for all the primes. Still, we can ballpark this: first note that combining registers is simply a matter of vporing them together - the remaining zeros are the likely primes. So it will take 29 vpor instructions to combine 30 registers, and the result will cover 256 candidate values, so that's ~8.5 candidates per instruction, versus ~1 for the counter approach. Now this is a very rough analysis and leaves out a lot of details like how do you get the bitmaps all aligned properly, but it seems like we may get about an order of magnitude improvement with this approach even over the branchless counter version. ### A C++ Prototype Let's try to prototype this in C++. Here's the core loop: #define NUM_PRIME 30 uint64_t accum = 0; for (int i = 0; i < NUM_PRIME; i++) { unsigned s = shift[i]; uint64_t bitmap = s <= 63 ? BITMAPS[i] << s : 0; accum |= bitmap; uint8_t prime = primes[i] + 1; shift[i] = (s + (prime - 64 % prime)) % prime; }  The BITMAP array is just an array of prime-specific bitmaps as discussed above: static const uint64_t BITMAPS[] = { /* 3 */ 0x9249249249249249 , /* 5 */ 0x1084210842108421 , /* 7 */ 0x8102040810204081 , /* 11 */ 0x0080100200400801 , /* 13 */ 0x0010008004002001 , /* 17 */ 0x0008000400020001 , ... /* 113 */ 0x0000000000000001 , /* 127 */ 0x0000000000000001 , };  Each bitmap is "normalized" such that the LSB is always 1. Then the loop is very simple: it loops over all 30 primes, and shifts the bitmap by the right amount to "align" it, and ORs all the results together. In this way, the loop handles 64 candidates. The shift amount is simply the amount need to correct stitch the bitmap from the last iteration so that it is periodic. Using a 16-bit example for sanity, the bitmap for prime == 3 in binary is 10010010 01001001. In the next iteration, you can use the same one since the effective 32-bit bitmap would be 10010010 01001001 10010010 01001001. Oops, two adjancent 1s! You just need to shift it over by 1: 10010010 01001001 << 1 == 00100100 10010011 and now it stitches fine. In general the stitch amounts for any prime p have period p and take all the values between 0 and p-1 inclusive. The last line in the loop calculates them. Let's try this guy:  perfc stat ./likely-primes time Bitmap1 69780348563
Finding all likely primes between 69780348563 and 131...
Count 7945632366, density=0.114 elapsed 139.631s, cycles/candidate 5.20, cycles/prime 45.69

Performance counter stats for './likely-primes time Bitmap1 69780348563':

139640.604802      task-clock (msec)         #    0.999 CPUs utilized
8,191      context-switches          #    0.059 K/sec
0      cpu-migrations            #    0.000 K/sec
99      page-faults               #    0.001 K/sec
362,181,261,309      cycles                    #    2.594 GHz
353,558,955,926      instructions              #    0.98  insn per cycle
34,926,378,434      branches                  #  250.116 M/sec
3,011,758,815      branch-misses             #    8.62% of all branches

139.741954558 seconds time elapsed


Huh, well that kind of sucked. At 5.2 cycles per candidate it's a bit slower than the algorithm in the OP. For one thing it still has ~9% branch mispredictions. It turns out the main culprit is this line:

uint64_t bitmap = s <= 63 ? BITMAPS[i] << s : 0;


This is just a "saturating" shift, which returns zero if the shift amount is 64 or more3 and which compiles to a branch4. It also turns out that the loop has no less than two very slow division instructions every time around, coming from the two % operators in this line:

shift[i] = (s + (prime - 64 % prime)) % prime;


The first term 64 % prime is really a constant (per prime) that could just be looked up in an array, but the second is more fundamental.

### Lookup Tables

They say every problem in programming can be solved by another layer of indirection, and every performance problem can be solved with a suitable LUT, so let's add at two layers of indirection and at least one LUT to solve all our problems.

Rather than calculating the shift amount each time and actually doing the shift, let's just do a byte-aligned load which already embeds the right shift amount. As it turns out, such a byte will occur in the bitmap within the first p bytes of the start. That is, for the prime 11, all possible alignments of the pattern appear starting in the first 11 bytes5. We keep track of the periodic pattern of the offset to load at using a small lookup table, and a "wrapping" counter approach. Here's the core loop:

    for (unsigned i = 0; i < NUM_PRIME; i++) {
uint8_t oidx = offset_indexes[i];
uint8_t offset = BYTE_OFFSETS[i][oidx];
accum |= bitmap;
oidx++;
offset_indexes[i] = (oidx == OFFSET_PERIODS[i] ? 0 : oidx);
}


The oidx is the "magic location" to load from the bitmap array which will have the correct alignment without needing shifting. We precalculate the offset values (they have a period p like almost everything else) and store them in BYTE_OFFSETS, and take the next offset on each iteration (the last line wraps when we hit the period).

How does this one do?

$perfc stat ./likely-primes time Bitmap2 69780348563 Finding all likely primes between 69780348563 and 131... Count 7945632366, density=0.114 elapsed 26.981s, cycles/candidate 1.01, cycles/prime 8.83 Performance counter stats for './likely-primes time Bitmap2 69780348563': 26982.332312 task-clock (msec) # 1.000 CPUs utilized 63 context-switches # 0.002 K/sec 0 cpu-migrations # 0.000 K/sec 100 page-faults # 0.004 K/sec 69,989,143,836 cycles # 2.594 GHz 242,084,078,467 instructions # 3.46 insn per cycle 18,541,869,190 branches # 687.186 M/sec 524,344,827 branch-misses # 2.83% of all branches 26.984215432 seconds time elapsed  Just 1 cycle per candidate, the fastest we've seen yet and more than 4 times faster than the original algorithm. All this in an non-SIMD algorithm that gcc doesn't even do an awesome job of optimizing (although not terrible either) - here's the above loop:  401860: movzx r10d,BYTE PTR [rsp+rdx*1+0x10] 401866: mov esi,edx 401868: shl rsi,0x7 40186c: movzx esi,BYTE PTR [r10+rsi*1+0x40ace0] 401875: mov rax,r10 401878: add eax,0x1 40187b: or rbp,QWORD PTR [rdi+rsi*1] 40187f: cmp al,BYTE PTR [rdx+0x407c20] 401885: cmove eax,r14d 401889: add rdi,0x86 401890: mov BYTE PTR [rsp+rdx*1+0x10],al 401894: add rdx,0x1 401898: cmp rdx,0x1e 40189c: jne 401860 <Bitmap2()+0x80>  ### Bring back the Vectors The next step is to vectorize this. This is getting long so we'll skip the first version (asm256 which clocks in at 0.27 cycles per candidate) and just go to my final version, asm512: .top: mov r15d, BYTE_BITMAPS512 xor ecx,ecx xor esi,esi vpxor xmm0, xmm0 vpxor xmm1, xmm1 .inner: movzx r14d,BYTE [rdi + rcx] mov eax,ecx shl rax,0x7 movzx edx,BYTE [r14 + rax + BYTE_OFFSETS512] vpor ymm0, ymm0, [r15 + rdx] vpor ymm1, ymm1, [r15 + rdx + 32] lea r14d,[r14+0x1] cmp r14l,BYTE [OFFSET_PERIODS512 + rcx] cmove r14d,ebx add r15,190 mov BYTE [rdi+rcx*1],r14l add rcx,1 cmp rcx,30 jne .inner vmovups [r12], ymm0 vmovups [r12+32], ymm1 add r12,64 cmp r13,r12 jne .top  The inner loop here runs once per prime and processes 64 bytes per iteration (512 odd candidates), in 14 instructions. The the two vpor instructions are doing the heavy lifting of combining the bitmmaks into the two ymm accumulators, and the rest is mostly just managing the indexes. The outer loop runs when the bitmap for all 30 primes have have been accumulated and stores the bitmap into a temporary buffer provided by the caller. We periodically break out of the asm code to examine the generated primes (in this case, simply counting) - see WrapBitmap for some details. In a real implementation, you would still want to do the handling periodically, but you might inline it right into the function. Let's time this guy: $ perfc stat ./likely-primes time asm512 69780348563
Finding all likely primes between 69780348563 and 131...
Count 7945632366, density=0.114 elapsed 4.387s, cycles/candidate 0.16, cycles/prime 1.44

Performance counter stats for './likely-primes time asm512 69780348563':

4388.416039      task-clock (msec)         #    1.000 CPUs utilized
18      context-switches          #    0.004 K/sec
1      cpu-migrations            #    0.000 K/sec
99      page-faults               #    0.023 K/sec
11,383,019,158      cycles                    #    2.594 GHz
34,299,568,928      instructions              #    3.01  insn per cycle
2,667,914,347      branches                  #  607.945 M/sec
72,436,051      branch-misses             #    2.72% of all branches

4.388903849 seconds time elapsed


We are down to 0.16 cycles per candidate! That's fully 25 times faster than the original algorithm, and if you measure it by cycles per prime, we are finding a prime every 1.44 cycles. Unless you are doing almost "zero work" per found prime, it's very likely that the other work will start to dominate here.

## Further Optimizations

If you are so inclined, this can still be made much faster, probably by a factor of 5 at least. Of course, before you pursue that, you would need to benchmark your full application, since it is highly likely that the unspecified work you do per prime is what is slowing this down now.

### Minor Optimizations

The loop above directly admits some minor optimizations. For example, rcx which counts off the 30 primes could be inverted so that it counts down to zero (or from -29 up to 0) allowing use to remove the cmp rcx,30 check at the end (we use the flag from the prior add instead). The shl rax,0x7 could be changed to a 3-argument shlx, avoiding the prior mov, or this whole calculation could be removed by using the induction counter r15 instead by making the row size of the two involved tables BYTE_OFFSETS512 and OFFSET_PERIODS512 consistent (right now one has an inner dimension of 128 and the other 190).

These may shave another small fraction of a cycle off of the existing time, but the ones below are much bigger.

The above asm512 algorithm reads uses vpor on two consecutive ymm registers worth of data (64 bytes) from the calculated index. It is in fact the slightly bigger brother of the not-shown asm256 variant, which only reads one 32B value in the inner loop. That guy ran at 0.27 cycles/candidate, so just doubling the read size in the loop nearly doubled the speed.

It's easy to see why: it only took one extra vpor instruction to do that, while the other 12 instructions in the loop are pure index calculation overhead which are now doing double work. So by increasing the loop by one instruction it does double the work.

You can just carry this idea to its logical conclusion, reading 4, 8 or however many values per loop. There is no particular reason it has to be a power of two, either. These will give very fast and easy speedups: I guess it is easy to get below 0.1 cycles/candidate using this approach.

The larger reads come at a size cost for the BYTE_BITMAPS - larger reads mean a larger table6. This optimization is probably the best and easiest one if you want performance. The code is already kind of half-generic.

I call this "unrolling horizontally" based on my mental model of each prime being a long horizontal bitmap, with primes stacked vertically one above another. So the vpor is accumulating in vertical slices (column-wise) and this unrolling moves in the horizontal direction.

### Unroll the Inner Loop

This is the "usual" unrolling and the counterpart of the horizontal unrolling discussed above. Currently the inner loop iterates over all 30 primes. This loop has a fixed trip-count, and it could be completely unrolled.

Several instructions in the inner loop would just disappear, such as all of the loop control, the instructions dealing with r15 and the shl.

This should give a reasonable one-time gain and the loop should still easily fit in the uop cache. It's less appealing than the horizontal unrolling since you can only do it once!

### Unroll the Outer Loop

Once you've unrolled the inner loop, you may want to unroll the outer loop as well. Unrolling this by N would result N copies of the unrolled inner loop so, the code would get big, fast, but I think you could probably unroll it by 3 or 4 and still fit it in the uop cache.

This allows some very interesting optimization since by unrolling the inner loop you now have unique sections of code handling each prime. When you unroll the outer loop, you may now be handling several reads for the same prime, in explicit unrolled code. The big win here is that you can directly hardcode the "offset sequence" that normally has to be painstakingly calculated by the generic code.

For example here's the start of BYTE_OFFSETS512 table for consecutive 64-byte reads:

extern "C" const uint8_t BYTE_OFFSETS512[][128] = {
/*   3 */  {  0,  1,  2,  0,},
/*   5 */  {  0,  4,  3,  2,  1,  0,},
/*   7 */  {  0,  1,  2,  3,  4,  5,  6,  0,},
/*  11 */  {  0,  9,  7,  5,  3,  1, 10,  8,  6,  4,  2,  0,},


The first row (for prime 3) means that after doing a 64 byte read at position 0, the next read should start at position 1 to properly stitch the bitmaps together, the next at 2, and so on. The row for 5 says jump from 0 to 4, then 3. Notice though that these are simple increments mod the prime. So the series for 5 is just "plus 4, mod 5": (0 mod 5), (4 mod 5), (8 mod 5), (12 mod 5). All the rows are similar.

So if you are reading 4 consecutive 64-byte bitmaps for prime 3, you can do it directly like this, without any index calculations at all (assuming rsi has the base of the LUT for the prime):

vpor ymm0, ymm0, [rsi]; jump 0
vpor ymm1, ymm1, [rsi + 32]; jump 0
vpor ymm0, ymm0, [rsi +  1]; jump 1
vpor ymm1, ymm1, [rsi + 33]; jump 1
vpor ymm0, ymm0, [rsi +  2]; jump 2
vpor ymm1, ymm1, [rsi + 34]; jump 2


You increased from 2 reads to 6 without adding any index calculation overhead! It's quite similar to the "horizontal unrolling" discussed above, where you just read more consecutive bytes (2 x 32B reads in this example), but that approach nearly doubles the table size (since you need to accommodate a 64B read at all positions), while this approach increases it very little (by 2 bytes, probably, since you just need to accomodate the maximum offset of +2 at every position).

The code is different per-prime: the code for 5 adds 4, or else subtracts 1 (they are equivalent, but you need to design your index handling to account for which direction you are going).

This approach is promising because it lets you get more reads without greatly increasing the table size. It works best for the smaller primes since the jump amounts (and hence required table padding) are smaller, while the horizontal unrolling described above works best for larger primes since the table increas is relatively less (since it is fixed to the read size and the large primes already have larger lookup tables).

Unrolling both the inner and outer loops allows you to even pick and choose different strategies for different primes.

### Optimize the Lookup Tables

I just made the tables 2D arrays, for simplicity, but this wastes a lot of space, since the rows for the smaller primes are often very short (I padded out some of them with 0xFF to help catch bugs). To optimize this, you'd probably want to first pack the tables more tightly, either as a jagged matrix (i.e., an array of pointers to rows), or just as one large packed 1D array. The latter approach is great if you've done a lot of unrolling as described above since the various vpor instructions can just directly embed the offset into the array in their memory operand "for free".

I didn't make any effort to align any of the lookup tables at all (and they won't have any natural alignment since they are all bytes. No doubt about 50% of the 32B loads will be "split loads" that cross a cache line. With a bit of care you can reduce that to 0% for the smaller tables with very little size increase. For the larger primes I think you can reduce it to 0% but at a 100% size increase (just thinking about it, I haven't checked), which may not be worth it.

Furthermore, depending on the number of consecutive bytes you are reading (see PB3_READ_BYTES) there are opportunities to dramatically reduce the larger tables. For example, if you are doing 64-bit reads like the non-SIMD bitmap algorithms, for any prime >= 67 every 64-bit read returns either an all-zeros value or a value with exactly one bit set. Yet such each such prime is using their own large lookup tables which are mostly zero. To support any possible 64-bit read for all primes >= 67 you need only 8 zero all bytes, and 8 other 15-byte regions with the all bytes zero other than the middle byte which 1 of the possible 8 bits set. You can overlap this all nicely so it takes about 72 bytes. So you can replace all so you can replace the 13 * 134 byte lookup tables for the primes from 67 to 127 by 72 bytes: a reduction of about 25 times!

Even better, this scales as you add larger primes: even if you want to add 100 more primes, you don't need any additional lookup tables for the bitmap.

For the fully generic algorithm which uses the BYTE_OFFSET table on every calculation, this transformation is free. For the unrolled versions that encode prime-specific knowledge into the reads it doesn't work as well. It also doesn't work as well for the larger reads: the version that reads 64B (512 bits) never gets close the "zero or 1 bit" set case for the first 30 primes, so you can't use it there. It would be useful if you wanted to use more primes, however and since this algorithm is so fast, it makes sense to do so.

### Combine Small Primes

Currently every prime is handled separately: although there could be some bitmap sharing for larger primes as described above, each prime still implies at least one vpor to incorporate it. There is nothing particular special about one-prime-per-bitmap, however: why not simply combine several primes together into one pre-calculated bitmap? Instead of having two bitmaps for 3 and 5 like:

p == 3  : 10010010010...
p == 5  : 10000100001...


Just use a combined bitmap with is the or of the two original bitmaps:

p == 3,5: 10010110011...


Now 3 and 5 are handled with half of the work. There is no need to stop at 2 primes either, you could include any number of primes in the pre-calculated bitmap.

So if it's somehow "free" to combined together primes, why am I mentioning this last? Can't we basically make the sieve as fast as we want by combining more and more primes? Not really. The main problem is that combining together primes, the period of bitmap increases to the product of all the primes. For example, by combining 3 and 5 into one bitmap, the new period is 3 * 5 == 15.

For larger primes or combing more than a couple primes, the period quickly becomes very large, requiring a large lookup table. To combine the first 4 primes, you'd have a period of 3 * 5 * 7 * 11 == 1155, much larger than the largest prime (127) in the original set. Furthermore, unlike single large primes, such bitmaps aren't very sparse (mostly zeros) so you can't optimize the tables in the same way as described above. For larger primes, like 113 and 127 the period for only those two primes is 14351, so it essentially can't effectively be used for primes of that size.

Still, it might be worth combining several of the small primes for a small boost if you've exhausted the avenues above. This technique would work very well if you want to use less than 30 primes, since the relative boost across a few small primes could be very big.

1 This makes total sense when you eyeball the code: the inner loop has 6 instructions, and the trips to the outer loop increase the branch density a bit.

2 Note that the output indicates the prime density is half that: 11.5% - and that's the true prime density - but the algorithm only examines half the numbers since it skips all even values, so from the point of view of the looping structure the prime density is 23%.

3 Shifts of a uint64_t by amounts larger than 63 are famously undefined behavior in C++, so this is needed for correctness, but even at the x86 assembly level, we'd need something because x86 shifts are "mod 64", so a shift by 64 is the same as a shift by zero, not what we want.

4 This would be much better as a conditional, but gcc doesn't do it that way, perhaps because there is a read of the BITMAP array on one branch and gcc doesn't want to do that in the case the value isn't used (even though it can probably prove that BITMAPS[i] is always in bounds).

5 This is just a consequence of 11 and 8 being relatively prime, and indeed since we are only dealing with odd prime numbers on the one hand and powers of two (for the various bitmap arrays), this useful properly will occur repeatedly.

6 At the limit of very large reads, the size of the table goes up proportionally to the read size, but for smaller values it is sub-linear. For example, when I moved from 1 ymm read in the avx256 algorithm to 2 reads in avx512, the size increased only from 158 bytes per prime to 190 bytes per prime. You can see the behavior by adjusting the PB3_READ_BYTES constant and running the likely-primes table command.

• re: table size reductions by taking any 64-bit run of zeros: If you're unrolling the loops enough, you can omit an OR with all-zero (although code-size is probably a problem at that point). An OR with a single set bit can be done with bt rax, 12 instead of actually loading anything. (Or movzx edx, byte [bitpos + rcx] / bt rax, rdx to store it as a table of bit-indices instead of a table of masks, but that costs two fused-domain uops instead of one for or rax, [masks+rcx*8]. Except on SnB/IvB where an indexed addressing mode will un-laminate even for or) Jul 28, 2017 at 0:14
• @PeterCordes The problem is that you mostly don't know when values are going to be zero or not based on simple unrolling because the behavior is periodic based on the prime. So if you are dealing with all primes in one loop, there is no unroll amount that is really useful. Sure, if you runrolled by 127, you'd always get zero in the same place for that prime, but not for any of the others. If you want to handle the next prime, 113 the zero behavior now has period 14351! That's why the indirection provided by the table lookups is so valuable. Jul 28, 2017 at 0:20
• Ah I see. I started skimming half way through, so I lost track of how the unrolling was working. Got it now. Jul 28, 2017 at 0:22
• They allow you to handle all the different periods in one loop, as long as you periodically handle the wrapping back to the start of the tables. Unrolling helps with hardcoding offsets into those tables (e.g., in the +1, +2, +3 example under "Unroll the Outer Loop") but it's harder to encode direct info about the primes. You can do it for the small ones (e.g, 3 * 5 = period 15), but the period quickly spirals out of control and it only shaves off a few primes: better to use your code size on other types of unrolling. Jul 28, 2017 at 0:23
• ... something I didn't add yet though, is that things change as you get to higher primes. The bitmaps are really great for the smaller primes - and you can combine several into one, which I didn't mention yet, e.g., 3 * 5 * 7 = 105 into one row whose bitmap reflects all the primes. They are probably still "great" for the primes like 127 which are the largest here, since at least you typically have several bits set in a 512-bit read, but overall this technique scales (runtime-wise) as O(N*M) where N are the number of primes you are checking and M is the number of candidates. Jul 28, 2017 at 0:25

I think you're underestimating the importance of the code outside the inner loop, and the mis-predicts of exiting the loop with different numbers of counts. With that iteration-count profile, it's very unlikely that regular throughput bottlenecks like dep chains or ALU uops are the major factor in your total run-time.

You're almost certainly getting a lot of branch mispredicts, and should use a profiler with perf counters to see exactly how bad it is. That will give us a better idea how much it's worth doing extra work to avoid branching. (See the branchless unrolling strategy at the end of this answer.)

Keep in mind that HW performance counters will account cycles to the instruction waiting for a slow result, rather than the instruction that was slow producing it. With an 8GB LUT, maybe your inner loop is waiting on cache-miss loads.

Asm style: indent your code one level farther than labels. I also like to indent operands to a mostly-consistent column, like you're already doing for comments. But if you end up using a really long mnemonic, I wouldn't always go and re-format all surrounding lines.

;; ugly, hard to see branch-targets
nottop:

vpminub ymm2, ymm2, ymm1                ; Adjust any underflows
vpsubb ymm2, ymm2, ymm0                 ; Decrement all the bytes by 1


;; tidier
nottop:
vpminub ymm2, ymm2, ymm1            ; Adjust any underflows
vpsubb  ymm2, ymm2, ymm0            ; -= 1


Some people even like to indent instructions inside loops, but most people don't.

Don't put alignment padding inside your critical loop! Especially on modern Intel CPUs with a uop cache, the align on the branch target probably gains you nothing. But the NOP hurts your front-end throughput every time through the loop.

The main reason for aligning branch targets is to avoid the branch target being near the end of a fetch-block, so fewer than the usual amount can be fetched and/or decoded in the first cycle. For CPUs without a uop cache, this is may mean a 16B naturally-aligned block of x86 code. e.g. 1 or 0 instructions might be decoded in the first cycle if the branch target was within a couple bytes of a 16B boundary. But even Pentium M / Core2 isn't restricted to fetching in naturally-aligned chunks, so the rules are very complicated. (And yes, front-end bottlenecks can often be an issue there. SnB's uop-cache solved an important problem for a lot of code.)

For CPUs with a uop-cache, things are very different. Usually it's not worth putting a nop anywhere that will actually execute. It's still a good idea to align separate functions, but usually aligning a short loop doesn't gain you anything (because the loop buffer handles it). Aligning branch targets / functions to not be right at the end of a 32B boundary can help, because the CPU can't read more than one cache line per clock, and uop cache lines end at 32B boundaries of the x86 code. An unconditional jmp / call always ends a uop cache-line.

Each uop cache line holds up to 6 uops, and a single 32B block can have at most 3 uop cache lines holding its decoded instructions. (So e.g. align 32 / times 19 nop will force that block to be re-decoded every time.)

Since you want to jump into the loop anyway, there might possibly be something to gain from aligning this way:

jmp top

align 32     ; increase this from 16 to 32
nottop:

vpminub ymm2, ymm2, ymm1                ; Adjust any underflows
vpsubb ymm2, ymm2, ymm0                 ; Decrement all the bytes by 1
;; align 8    ; REMOVE THIS
top:
vpmovmskb r10d, ymm2                    ; Check the upper bits for underflows

test r10d, r10d
jnz nottop                              ; No underflows means prime-ish


The whole inner loop will probably be in one line of the uop cache. The first fetch cycle after jumping to top will probably only fetch the last two uops (vpmovmskb and macro-fused test/jnz), but at least you can be pretty sure they won't be in separate cache lines. IDK if that was likely with align 8, but it certainly wasn't helping. Even with a long-NOP in the middle, the loop was only 6 uops, so it still would have fit in a single cache line. Aligning to 32 is overkill, but since code-size probably isn't an issue for your whole program (including omitted code in the outer loop), I would just go for align 32 there since you're always jumping over it. (Fun fact: the code up to / including the first jmp top is 7 uops, so saving one by folding a load (see below) is likely beneficial to uop-cache footprint).

If x86 code size is an issue (pushing you over a 32B boundary and ending a uop cache line at an inconvenient place), look at using esi or edi instead of r9d/r10d, to avoid a REX prefix. And replacing add 1 with inc would also save another byte. There is no performance downside to inc for HSW/SKL/KBL. Or even for KNL / Silvermont, because you're not reading those flags.

HSW can have uop-cache throughput problems more easily than SKL/KBL. SKL increased uop-cache read throughput to 6 uops per clock (read from uop-cache aka DSB, added to IDQ). I'm not sure I'm understanding this right, but it sounds to me like pre-Skylake could end up reading 4 uops, then 2 uops, then 4, then 2, when reading from uop cache lines that are each filled with 6 uops. (This doesn't usually happen; 32B boundaries and multi-uop instructions end uop cache lines early.)

All-ones can be a broadcast-load, or cheaply generated on the fly with ALU ops (actually just one is enough in your case). Hard to guess which is more efficient for the surrounding code. If you branch-mispredict a lot, OOO execution might have a hard time hiding the ~5 cycle load-use latency before the loop needs it. But if you're really just bottlenecked on ALU throughput, then vpcmpeqb still needs a vector ALU uop (but it is dependency-breaking, so can run as soon as OOO execution sees it and has a spare cycle).

;; load version: only needs a db 1,1,1,1 constant

;; ALU version:
vpcmpeqb     ymm0, ymm0, ymm0    ; set1(-1)


In fact, even if you broadcast-load, you should still load set1(-1), so you can fold this:

vmovdqa ymm2, [r8]
...
vpsubb ymm2, ymm2, ymm0


into this:

vpaddb ymm2, ymm0, [r8]     # 1 micro-fused uop for the front-end


lea r9, [r9 + r9 + 2] is unnecessary. You can just initialize with mov r9d, 2, and use add r9d, 2 inside the loop.

Also, the sub rcx, r9 (+ potentially sbb) / cmp rcx, rdx / jb done check could be simplified. If you compute end - start once at the top of the function, a single 64b register can hold the difference, for any reasonable problem-size. (2^63 / 97 is longer than you'll ever want to let the program run for, to say the least.)

Then you can just check cmp r9, rcx / jae done. Or just decrement end-start register until it wraps (which you could detect with test rcx,rcx / jl. (Not js because that doesn't macro-fuse). Or look at flags set by the last decrement. But probably don't put a dec&branch into the hot loop in place of the add r9d,1, because then you might bottleneck on port6.

vzeroall is (sometimes much) slower than vzeroupper. On SKL/KBL: 34 uops (in 64-bit mode) vs. 4. (vzeroupper kinda sucks on AMD, but it's still faster than vzeroall. They're equally slow on KNL.)

You can and should just leave garbage in the low 128 when you're done, and just use vzeroupper before calling or returning to something that might execute legacy SSE instructions.

It would be nicer if you could leave this for the compiler to do or not, by using _mm256_zeroupper() at the end of your calling C++ function. That way, it might possibly optimize away if called in a loop or something. Or if compiled for KNL. I'm not sure whether the intrinsic guarantees that you'll get a vzeroupper or if it's just a hint to the compiler that AVX-upper state is dirty at this point. It's not usually a problem when using intrinsics instead of calling an asm function, since then the compiler knows how the ymm regs were used.

Alternate strategy to handle the counter: do it in one of the vector elements you aren't using. Getting the counter out for use in the done check would suck though. If your loop typically ran more than 20 iterations, it would probably be worth the extra vmovd + movzx to extract and vpinsrw or vpand to zero between invocations of the inner loop, to avoid an extra integer add uop in the loop. Or have a separate vpminub with zeros instead of 255 in those two elements, for use in a peeled loop-entry iteration.

In your case, I think the extra add will be better than trying to mess with the last 2 vector elements. Especially since out-of-order execution can have that value ready well ahead of the vector being ready, if the loop-exit branch-predicts correctly.

## Unrolling

With the overhead of an integer add in your loop, you should look at unrolling a bit. With multiple branch-exits, you can work out the odd increments by having the jumps out of the loop go to different places than the fall-through.

Branch prediction may be better or worse overall. On the plus side, the pattern of how many iterations the loop runs is spread over more branches. e.g. when unrolling by 2, a pattern of 1, 2, 1, 2 would mean one branch is always taken, the other branch is never taken. But OTOH, it may make the worst case worse: if the loop had recently run mostly low-iteration-count loops, multiple branches might all mispredict while getting into a 97-iteration loop. Of course, the same branch could predict wrong multiple times in a row in a fully-rolled loop. Also, modern branch predictors can look at patterns of which way other preceding branches went.

So overall I would hope that unrolling doesn't make branch-prediction any worse, and might help. But check with perf counters.

;; untested
;; comments mostly stripped for brevity.  Keep them in for the real version

align 16        ; probably already aligned, but this shouldn't hurt
ProcessA proc

vpcmpeqd   ymm0, ymm0,ymm0     ; set1(-1)
vmovdqa    ymm1, YMMWORD PTR [asmprimes]

; This may require sbb for values greater than 64 bits, but the result should be 64b
sub   rdx, rcx                          ; start-end

xor   eax, eax                          ; Returns # of prime-ish found
mov  r9d, 2                             ; How many steps to next prime-ish?  We already do one while setting up
jmp loop_entry

align 32   ; try with / without this
innerloop:

add       r9d, 2*2                       ; increment by 2 times the unroll factor
vpminub   ymm2, ymm2, ymm1               ; Adjust any underflows
vpaddb    ymm2, ymm2, ymm0               ; += -1
loop_entry:
vpmovmskb r10d, ymm2                     ; Check the upper bits for underflows
test r10d, r10d
jz loop_exit_odd                         ; No underflows means prime-ish

vpminub   ymm2, ymm2, ymm1
vpmovmskb r10d, ymm2
test r10d, r10d
jnz innerloop
;; fall out into loop_exit_even:
add  r9d, 2*1               ; Add another 2, since we start the loop with r9d=2
loop_exit_odd:
; r9d = 2 + 4*n because we skip the final +=2 when exiting the loop from the middle

; Decrement how far above end we still are
sub   rdx, r9
jc done               ; if (start < end) goto done;

;; peeling this out of the loop makes some sense.
mov     r9d, 2                           ; Reset for next loop
vpaddb  ymm2, ymm2, ymm0                 ; += -1

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; This is a placeholder for a more involved block of code.
inc rax ; Count how many we find
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

jmp loop_entry

done:

vzeroupper
ret


The loop-entry point is before a branch that's not-taken unless the loop exits right away (when the peeled single iteration outside the loop was enough). This is probably better for the front-end. Peeling the vpmovmskb / test/jcc for that case might be a good idea, e.g. to replace the jmp loop_entry right before done with a jcc to the sub rdx,r9. Or not, because it would only be taken 15% of the time, which is enough to mispredict a lot. But anyway, that would let you jump to inner_loop instead of right to a loop-condition, which is maybe better for the front-end.

For larger iteration counts, the add r9d, 4 is obviously amortized nicely. Arranging things to minimize the add instructions for 0, 1 and 2 iterations was the interesting trick. Initializing to 2 helps, but then figuring out where to place the loop entry point and the add r9d, 4 was tricky. Also arranging the loop-exit code to not run too many extra instructions. If the loop always ran many iterations, I wouldn't worry about having one way of leaving the loop end up running an add 4 and a sub 2, but we could avoid it here with no downsides.

Unrolling by more than 2 with this technique would result in a chain of add r9d, 2 instructions after the loop. It might be better to have the non-fall-through loop exits jump to blocks at the end of the function like add r9d, 6 / jmp loop_exit. You can do some code-duplication here to make each of those blocks do a cmp rcx,r9d/jcc big_block/vzeroupper/ret, so any given path doesn't have extra unconditional branches. You're replacing the jmp loop_exit with code equivalent to the jcc done that you do in that path anyway. After that, jmp done or actually duplicating the vzeroupper/ret doesn't matter.

## Make some of the unrolled inner loop branchless

For example, do two iterations before checking for exit, then sort out afterwards when you should have exited (and use that counter value). Since we need to optimize for very low iteration counts, we can't have too much sorting-out overhead. But if we do it branchlessly, we might significantly reduce branch mispredicts by having 0 and 1 iterations taking the same path through the outer loop. So instead of 11 different paths covering 95% of the cases, we maybe only have 6.

This technique works great for something like strlen, where you can OR together a whole cache-line of pcmpeqb results before branching on there being any zeros, so you bottleneck on loads instead of ALU uops even with data hot in L1D cache.

Give this a try; if it helps branch-prediction a lot, then look at finding ways to optimize this more. Maybe branching outside the loop (on odd or even) wouldn't be bad, although that's probably no better than just having two loop exits. Maybe using a branch to fixup ymm2 after loading from the LUT would be fine. A branch mispredict during load-use latency of a cache miss should be much better than a mispredict before the load instruction. Or not, because the line will probably still end up in L1 cache and be there when the code runs with the correct values in vector regs (since that code doesn't use them). Branchless handling of just r9d is cheap enough.

    mov  r9d, 4                          ; 1 peeled iteration + 1 in the bottom half of the loop.
jmp loop_entry

align 32   ; try with / without this
innerloop:

add       r9d, 2*2                   ; increment by 2 times the unroll factor
vpminub   ymm2, ymm2, ymm1           ; Adjust any underflows
vpaddb    ymm2, ymm2, ymm0           ; += -1
loop_entry:
vpmovmskb esi, ymm2                 ; Check the upper bits for underflows
;; unconditionally do another iteration in different regs, and resolve it later

vpminub   ymm3, ymm2, ymm1     ; use ymm3 for 2nd half, leaving the ymm2 result from the previous iteration
vpmovmskb r10d, ymm3           ; 32 bytes -> 32 bits
imul      r10, rsi             ; r10==0 only if either vector was all-zero.
test      r10, r10             ; needed because imul leaves ZF undefined.  Basically free because it macro-fuses, though
jnz  innerloop

;; r9d = 4 + 4*n
;; If esi == 0, we should have stopped before doing the last step
;; and need to fix ymm2 and do r9d -= 2

sub   esi, 1        ; sets CF only if esi==0
sbb   esi,esi       ; -1 or 0.  False dep on esi, but it has to be ready if flags are.  2 uops on Intel pre-Broadwell, but this should be good on SKL
;; Note that esi is not sign-extended into rsi, but that's fine as long as we're doing 32-bit math only.
lea   r9d, [r9 + rsi*2]    ; r9d += 2 * (-went_too_far)
;; The final r9d can't be negative, so we can still use it as a 64-bit unsigned value even though that LEA couldn't wrap.  This is why it's safe to skip the REX on sbb rsi,rsi

; fix ymm2 using the all-ones/zeros from SBB as a vector blend mask
vmovd ymm4, esi
vpbroadcastd ymm4, xmm4     ; AVX512 has integer->vector broadcast in a single instruction, but AVX2 only has xmm/m32 source

vpblendvb  ymm2,  ymm3,ymm2,  ymm4   ; take src2 when the bit in ymm4 is 1


Detecting a mispredict early is important, and imul is 3c latency, so that's not ideal. It may be better than more branches, though. imul only runs on p1, so it competes with vector ops. It's a pretty good way to check if either of two arbitrary 32b values are zero. AND/TEST don't work (e.g. 0xf0 & 0x0f is 0). Using setcc to get two boolean registers and then testing them against each other would work, but take even more uops.

The code after the loop is 3 uops (sub/sbb/lea) to fixup r9d + 4 uops (vmovd/vpbroadcastd/vpblendvb) to fixup ymm2 on Skylake/Kaby Lake. I used sbb instead of setcc because it was useful to have an all-ones to get into vectors. And also useful to have a negative value that I could multiply by 2 with an LEA. I could have done it the other way, and used setcc + lea to add an extra 2 if esi was non-zero. But that wouldn't help with the fixup for ymm2.

Doing the check in vector regs would be nice, instead of needing to go int->ymm to get a mask for vpblendvb. But we effectively need a horizontal-OR before we combine two vectors with anything like AND. Not even ptest with its ZF and CF outputs can check two arbitrary registers for either being zero.

• Whoa. When you create an answer, you really create an answer. Hard to respond to a post this long using just comments, but I'll do my best. "underestimating outer loop" - While I'm fully prepared to believe that VS's profiler is lying to me, in this case I'm pretty sure. Excising the actual work from the loop still leaves the majority of runtime. "Asm style" Fair enough. I just hate fighting VS's handling of tabs, but you're right, it's time. Jul 22, 2017 at 0:18
• "Align" - Aligning jump targets is supposed to be a "good thing," but I never saw (either of) the alignments make any measurable difference (plus or minus). I left it there to remind myself that I'd thought about it. "broadcast/vzeroall" - Useful tips for SSE newbies such as myself, but probably not significant here. "add r9d, 2" the lea is a leftover from in iteration that used xor & inc. The idea being that xor is 'free' and a single lea satisfies the need for both mov r9,1 and the shl r9,1. None of which made any difference. Jul 22, 2017 at 0:18
• "overhead of an integer add" - In theory, I can use one of the 2 unused bytes in ymm2 to do the add. In practice, the vpextrb I used to get the value out combined the the vpand to reset it turned out to be more expensive. Possibly the vpsubb/vpminub leave enough latency that the add is free? "end - start" While I don't use any xmm/ymm registers in the placeholder, I do use 'start' (rcx). Jul 22, 2017 at 0:19
• Ok, when I fixed a problem related to calculating the start/end point, I broke something else. Fixing the fix gives the right answer. It now calculates the right answer, and does so ~10% faster than the OP using the test framework. Upvote, but I'm holding off on accepting til I see if someone else has an alternative. Jul 22, 2017 at 1:35
• @PeterCordes - I was able to get perf numbers for the code. It is crushed by branch misprediction: 20% mis-predict rate, and probably 50% of the runtime due directly to mispredicts (both as reported by toplev.py and ballpark given it's about 40 cycles per prime and ~1 mispredict per prime at lets say 20 cycles). I'll add more in an answer. Jul 22, 2017 at 22:59

When I was composing this question, I kept thinking about what questions people were going to ask me: Did you think about A? Did you try B? Why didn't you go with the obvious thing and do C?

Normally, this is a good thing, since I can often answer my own questions without ever having to hit the "Send" button.

In this case, there were 2 questions buzzing in my head as I was writing the post. Should I wait until I answered them before I post? Heck, they might not even work, right? I decided to ignore them and post what I had. But they were still bugging me.

So while I was nervously awaiting the "Hey dummy, why don't you just..." response I was sure was about to come (I HATE looking foolish in public), I decided to pursue them:

1. If the "outer loop" is so fast that the SSE enumeration of primes is the limiting factor, why not go with something drastically simpler? At its simplest, I could just walk all odd numbers. Somewhere between those two extremes there's got to be a balance.

2. Using this SSE count-by-ones thing, there's a fixed number of permutations, right? For example, if I were only using the first two primes (3, 5), then there are only 15 possible values for those two, and I'm going to keep cycling thru those 15 permutations over and over again. Heck, I could just save those values and avoid constantly computing them.

Now, for 30 primes, there are a few more than 15 permutations (closer to 2,007,238,469,666,518,094,547,220,599,513,022,568,322,942,623,865).

But what if I pick something in between? At 9 primes, there are 3,234,846,615 permutations. And I don't really need to store all 3gig. I only need the likely primes and that's just 1,021,870,080. So if I were to store the counts (aka r9 from ProcessA) for each permutation, I would always have the "next" count instantly available from a (sequential) memory read. No computation required.

And the great thing is, the contents of the table are not dependent upon what prime number I want to start computing from. Ok, there's a bit of a trick figuring out where in the table to start for a given prime, but once I start, I just keep going, reading the counts one byte at a time, and looping back to the beginning of the table as necessary. So the table can be computed once, and stored on disk.

The down side is that I'm only checking 9 primes. That's a disappointing letdown from 30.

Curious, I checked the values for 10:

Permutations: 100,280,245,065 (100 billion)
Liklies: 30,656,102,400 (30 billion)


Big, but still small enough to fit in my 64gig of RAM. Hmm. What's the next one? Maybe I can use a memory mapped file:

Permutations: 3,710,369,067,405 (3.7 trillion)
Liklies: 1,108,988,592,192 (1.1 trillion)


Ok, maybe not. 10. 10 isn't so bad. Yeah, 10 is good.

I could post the code for creating the file. But it's still pretty ugly. And there's not much point in cleaning it up, since I only ever need to run it once.

The code for "figuring out where in the table to start" is a bit trickier. I've got some thoughts for speeding it up, but it only runs once per thread at startup. Relatively speaking, performance is not a big consideration here.

And the code for walking the table? It's almost too stupid to post, and certainly not worth writing in assembly anymore:

while (start >= end)
{
start -= steps[pos];
pos++;

if (pos >= iBuffSize)
pos = 0;

res++; // This is a placeholder for a more involved routine.
}


Yes, I'm planning to reverse the buffer so I can count down instead of up.

This has chopped the time to walk the primes down to < 10 seconds. Damn, my "involved routine" has probably just become my bottleneck!

Does that mean the effort by Peter and others to tighten up the SSE code was a waste of time? Not a bit of it.

First of all, the fact that I kept thinking about it is what led me to this new approach.

Secondly, while the "involved routine" is currently lightweight, I feel quite confident that this is the one part of the code that will grow more complex over time. As the cost of those computations goes up, the benefit from filtering out more primes will become more important.

And lastly, I can already hear someone asking me (see @BeeOnRope's comment about p mod 6): Is there some way to combine these two? Hmm.

I'm not going to 'accept' this answer, since it doesn't really use "SSE loop to walk likely primes." But it is (currently) the direction I intend to pursue.

I put it here in case it helps the next guy.

• If you can vectorize the "more involved" code at all, you can easily update a vector of start values. (Especially if your steps table entries store the sum of the 4 preceding values, so you don't have to do any prefix-sum horizontal stuff and can just do startvec = _mm256_sub_epi32(startvec, pmovzxbd(load(steps+pos))) instead of doing horizontal ops after loading.) Depending on what the more involved routine is like, vectorizing with 128b vectors might make it easier (cross-lane shuffling costs more).) Jul 22, 2017 at 23:37
• To be fair, it is hard to evaluate the performance of this answer because it depends heavily on the details of the "other work" done for each likely prime which is "too complex to reveal". Then closing the door by saying "this is the direction I'll pursue" is certainly discouraging to other posters! Jul 23, 2017 at 4:38
• "hard to evaluate" That's fair. A proposal to "check every odd" might be faster given this dinky placeholder. OTOH, if the SSE is the bottleneck, a solution that eases this code by putting more load elsewhere seems (in retrospect) an obvious path. "closing the door" I did say "currently," which implies (to me) "subject to future events." My Q was about SSE, so if someone's got a better solution there, I'm totally prepared to both upvote and accept it. I'd also be prepared to upvote (and use) an answer like mine that uses more primes or is (somehow) even faster. Jul 23, 2017 at 7:29
• "too complex to reveal" - (Mostly) this isn't due to my desire to be mysterious. ;) Programs tend to be a collection of interconnected pieces. The code that follows depends on a table, which requires some other code, which requires CUDA, etc. As I'm not looking for someone to re-write my entire app for me, I selected the one piece that seemed to most important and created a framework to demonstrate it. Yes, this makes the relative expense of computing 'N' more likelies impossible to define. Not sure what I can do about that, esp since the expense depends upon the value being processed. Jul 23, 2017 at 7:30
• Fair enough. I would note, however, that you could perhaps replace the code with more realistic dummy code, e.g., - the table lookup could still be there but full of zeros (depending on exactly what is done with that value, etc). Because I think I can make the "likely prime" finding much faster, but a lot the interesting interaction happens at the boundary between the inner (likely prime checking) and outer (processing found likely-primes) loops. Jul 23, 2017 at 20:51