SoE

//    https://codereview.stackexchange.com/q/189083

//    https://pastebin.com/JMdTxbeJ
//    Segmented SoE (primes < 2^32 in 4.6 seconds).

/*                    old-ns-new               old-s-new
      primes <   25     47    47    p < 2^31   7.4   7.0
             <   50     62    63      < 2^32  15.8  15.0
             <  100    109    78
             <  200    187   141
             <  400    312   218
             <  800    578   375
             < 1600   1092   718
             < 3200   2152  1435
             < 6400   4274  2855                           */

#include "stdafx.h"
#include <vector>
#include <iostream>
#include <ctime>
#include <Windows.h>
using namespace std;
typedef unsigned int u32;
vector<u32> x; vector<u32> p;

void buildX(u32 m)  // mark odd composites
{
    m -= m / ~0u; m += m & 1; m >>= 1;
    u32 a = 1, b = 2, c = 3, d = m >> 5;
    x.resize(d + 1); x[0] = 0x9b4b3491;
    for (; a <= d; a += 3) x[a] = 0x24924924;
    for (; b <= d; b += 3) x[b] = 0x49249249;
    for (; c <= d; c += 3) x[c] = 0x92492492;
    for (a = 1; a <= d; a += 5) x[a] |= 0x42108421;
    for (a = 2; a <= d; a += 5) x[a] |= 0x10842108;
    for (a = 3; a <= d; a += 5) x[a] |= 0x84210842;
    for (a = 4; a <= d; a += 5) x[a] |= 0x21084210;
    for (a = 5; a <= d; a += 5) x[a] |= 0x08421084;
    for (a = 7, b = 24, c = 12; b < m; a += 2, b += c += 4)
        if ((x[a >> 6] & 1 << (a >> 1)) == 0)
            for (d = b; d < m; d += a) x[d >> 5] |= 1 << d;
    x[m >> 5] |= ~0u << m;
}

void countPrimes()
{
    u32 c = 1; int i = x.size() - 1;
    while (i >= 0) c += __popcnt(~x[i--]);
    p.resize(c);
}

void primes(u32 m)  // primes <= m
{
    if (m > 1)
    {
        buildX(m); countPrimes(); p[0] = 2;
        u32 u = 1, v = 1, xi; DWORD tz;
        for (int i = 0, n = 1, c = p.size(); n < c; v = u += 64)
            for (xi = ~x[i++]; xi;)
            {
                _BitScanForward(&tz, xi); xi >>= tz; xi >>= 1;
                p[n++] = v += tz << 1; v += 2;
            }
    }
}

int main()
{
    for (u32 m = 25; m <= 6400; m <<= 1)  // for (u32 m = 0; m < 9; m++)
    {
        primes(m);
        if (m > 1)
        {
            u32 maxP = p[p.size() - 1];
            cout << "largest prime <= " << m << " : " << maxP << "  ";
        }
        clock_t clock0 = clock();
        for (int i = 1000000; i; i--)
        {
            x.resize(0); p.resize(0); primes(m);
        }
        cout << clock() - clock0 << " ns: ";
        cout << p.size() << " primes" << endl << endl;
    }
    x.resize(0); p.resize(0); u32 m = ~0u >> 1;
    clock_t clock0 = clock(); primes(m);
    cout << clock() - clock0 << " ms: ";
    cout << p.size() << " primes <= " << m << endl;
    getchar();
}