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.
...
...
BAD Bad_Speculation: 63.49 % Slots [100.00%] BN
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
add rax, r14
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];
uint64_t bitmap = read64(BYTE_BITMAPS[i], offset);
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.
Larger Contiguous Reads
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.
lea [r9 + r9 ...]line), but to handle 2 and 3, you need to do something different for theP mod 6 == 1andP mod 6 == 5cases. Unroll by 2 or table lookup could both work. \$\endgroup\$