Sieve of Eratosthenes for small limits

Question

I coded a sieve of Eratosthenes, I think it's quite fast for small limits. Any suggestions how to improve it? Especially for small limits, say limits < 50,000,000. (Also coded a parallel sieve (in C#) which finds the primes < 2^32 in 1.6 seconds, but it's relatively slow for small limits).

How it works

Let's find the primes <= 75. Odd composites are marked in an unsigned int array "x", each bit represents an odd number. The first (hexadecimal) number in x:

x[0] = 0x9B4B3491

• Binary: 1001 1011 0100 1011 0011 0100 1001 0001
• The lowest bit is set: 1 is not prime. The next bit isn't set: 3 is prime, etc.

x[1] = 0xFFFFFFE5

• Binary: 1111 1111 1111 1111 1111 1111 1110 0101
• Primes: 67, 71, 73

Count the primes, the 0-bits in x, with __popcnt(unsigned int ...), (population count, number of set bits). __popcnt(x[0]) gives the number of set bits, the number of odd composites. __popcnt(~x[0]) gives the number of zero bits, the number of odd primes.

~x[0] = 0x64B4CB6E

• Binary: 0110 0100 1011 0100 1100 1011 0110 1110

~x[1] = 0x0000001A

• Binary: 0000 0000 0000 0000 0000 0000 0001 1010

~x[0] has 17 set bits, ~x[1] has 3, 20 odd primes. Resize the prime array "p" to 21. That's 20 odd primes and 1 even prime (2).

// Put the primes in "p", p[0] = 2
// Initialize: v = 1
// Count trailing zeros (tz), shift xi right by tz + 1:
//       xi = ~x[0] = ....101101101110  tz = 1, p[1] = v += tz << 1 =  3
//                                                 xi >>= tz + 1, v += 2
//                          1011011011  tz = 0, p[2] = v += tz << 1 =  5
//                                                 xi >>= tz + 1, v += 2
//                           101101101  tz = 0, p[3] = v += tz << 1 =  7
//                                                 xi >>= tz + 1, v += 2
//                            10110110  tz = 1, p[4] = v += tz << 1 = 11
//                                                 xi >>= tz + 1, v += 2
//                              101101  tz = 0, p[5] = v += tz << 1 = 13
//                                                 xi >>= tz + 1, v += 2
//                               10110  tz = 1, p[6] = v += tz << 1 = 17
//                                                 xi >>= tz + 1, v += 2
// 
// Results (Intel i7-4790 CPU @3.6 GHz from 2015):
// primes <= 100   109 ns
//        <= 200   187 ns
//        <= 400   312 ns
//        <= 800   578 ns
//        < 2^31   7.4 seconds
//        < 2^32  15.8 seconds

#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
{
    // make it safe:  1 << (int)(a >> 1)  =>    1 << (int)(a >> 1 & 31) ?
    //                1 << (int)d         =>    1 << (int)(d & 31)      ?
    //              ~0u << (int)m         =>  ~0u << (int)(m & 31)      ?

    m -= m / ~0u; m += m & 1; m >>= 1;
    x.resize((m >> 5) + 1); x[0] = 1;
    for (u32 a = 3, b = 4, c = 4, d; b < m; a += 2, b += c += 4)
        if ((x[a >> 6] & 1 << (int)(a >> 1)) == 0)
            for (d = b; d < m; d += a) x[d >> 5] |= 1 << (int)d;
    x[m >> 5] |= ~0u << (int)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, j = p.size(), n = 1;;)
        {
            xi = ~x[i++];
            while (xi)
            {
                _BitScanForward(&tz, xi); xi >>= tz; xi >>= 1;
                p[n++] = v += tz << 1; v += 2;
            }
            if (n >= j) break;
            v = u += 64;
        }
    }
}

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;
    clock_t clock0 = clock(); primes(m);
    cout << (clock() - clock0) * 1e-3 << " s: ";      //  15.709 s
    cout << p.size() << " primes <= " << m << endl;

    x.resize(0); p.resize(0);
    clock0 = clock(); primes(m);
    cout << (clock() - clock0) * 1e-3 << " s: ";      //  15.522 s

    // 15.709 - 15.522 = 0.187 => memory allocation makes it ~ 1.2 % slower?

    cout << p.size() << " primes <= " << m << endl;
    for (int j = p.size(), i = j - 5; i < j;) cout << p[i++] << endl;
    getchar();
}

Answer 1

Overall the code is a bit hard to read with the multiple statements per source line, multiple variable declarations (of different types) on one source line, and non-descriptive variable names.

You use a lot of hard coded constants, assuming that you're working with 32 bit integers. Using named constants for these would make it clearer where some of your numbers come from, as well as allowing for an easier expansion to a 64 bit compilation.

The expression 1 << (int)d used in buildX results in Undefined Behavior when d is larger than the number of bits used in the left operand (32 in this case). Many hardware implementations of shift left reduce the shift count modulo the number of bits in the operand being shifted, so this is likely doing what you expect, but that is not guaranteed.

Since your vectors already have all the memory they need allocated when you enter your 1 million iteration timing loop, your times are solely for your algorithm and do not consider the time required for the memory allocations as the vectors expand. The timing for your maximum case, since it only runs once, includes time spent growing the vectors. (It's just 2 allocations, but is a difference in the timings.)

Answer 2

A faster "buildX(u32 m)". Odd multiples of 3 are written with three 32-bits words in "x", starting at x[1], a simple solution for the first word: x[0] = 0x9b4b3491; In a similar way the next prime(s) 5,7.. could be marked, but with deminishing returns.

/*                   old-ns-new               old-s-new
     primes <   25     47    47    p < 2^31   7.4   7.1
            <   50     62    62      < 2^32  15.8  15.2
            <  100    109    94
            <  200    187   156
            <  400    312   250
            <  800    578   436
            < 1600   1092   843
            < 3200   2152  1700
            < 6400   4274  3385                          */

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 = 5, b = 12, c = 8; 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;
}

Latest version: https://pastebin.com/JMdTxbeJ
Primes < 2^32 in 4.2 seconds.