Sieve of Eratosthenes in C++ with wheel factorization

Question

I have written a completely working C++ program for the first time, and I have managed to compile it. I have decided to learn C++ and I have written a C++ program to familiarize myself with C++, this is my first ever C++ program and I wrote it all by myself, having only started learning a few hours prior. (And I haven't written a "Hello, World." program.)

The purpose is to practice my C++ skills and put what I have learned into use, to test if I really know the concepts I used in the script.

My idea is simple, the program accepts two command-line arguments (no std::cin and std::cout, I dislike interactivity), the first argument is the limit up to which primality of numbers should be checked, the second is the file path to which the prime numbers should be saved to.

The program should find all prime numbers up to the given limit using Sieve of Eratosthenes with Wheel Factorization optimization and write the prime numbers to the file line by line.

Here is the code:

#include <cstdlib>
#include <fstream>
#include <iostream>
#include <stdexcept>
#include <vector>

const int wheel[8] = { 4, 2, 4, 2, 4, 6, 2, 6 };
const int triple[3][2] = { {4, 2}, {9, 6}, {25, 10} };

std::vector<uint64_t> wheel_sieve(uint64_t limit)
{
    limit++;
    std::vector<bool>is_prime(limit, true);
    is_prime[0] = is_prime[1] = false;
    for (auto& pair : triple)
    {
        int square = pair[0];
        int multiple = pair[1];
        for (uint64_t i = square; i < limit; i += multiple)
        {
            is_prime[i] = false;
        }
    }
    uint64_t k = 7;
    int i = 0;
    while (true)
    {
        uint64_t square = k * k;
        if (square > limit)
        {
            break;
        }
        if (is_prime[k])
        {
            uint64_t twice = 2 * k;
            for (uint64_t i = square; i < limit; i += twice)
            {
                is_prime[i] = false;
            }
        }
        k += wheel[i];
        if (++i == 8)
        {
            i = 0;
        }
    }
    std::vector<uint64_t> primes;
    for (uint64_t i = 0; i < limit; i++)
    {
        if (is_prime[i])
        {
            primes.push_back(i);
        }
    }
    return primes;
}

int main(int argc, char* argv[])
{
    if (argc != 3)
    {
        throw std::invalid_argument("number of arguments must be 2");
    }
    uint64_t limit = atoi(argv[1]);
    std::vector<uint64_t> primes = wheel_sieve(limit);
    std::ofstream file;
    file.open(argv[2]);
    for (uint64_t& prime : primes)
    {
        file << prime << std::endl;
    }
    file.close();
    return 0;
}

I have managed to compile it and run it, and I have confirmed its correctness. I compiled it using Visual Studio 2022 in x64 Release mode.

I guess it can be faster?:

PS C:\Users\Xeni> measure-command {C:\Users\Xeni\source\repos\Wheel_Sieve\x64\Release\Wheel_Sieve.exe 1048576 D:\primes-under-1048576.txt}

Days              : 0
Hours             : 0
Minutes           : 0
Seconds           : 0
Milliseconds      : 295
Ticks             : 2950071
TotalDays         : 3.41443402777778E-06
TotalHours        : 8.19464166666667E-05
TotalMinutes      : 0.004916785
TotalSeconds      : 0.2950071
TotalMilliseconds : 295.0071

How can I make it more efficient?

Answer 1

Stop using std::endl

Use '\n' instead of std::endl; the latter is equivalent to the former, but also forces the output to be flushed. This is often unnecessary, and hurts performance.

Reserve space for vectors

If you use std::vector's reserve() member function, it can allocate all the memory it will need in one go. If you don't do it, it will have to reallocate several times while you push_back() items to it.

Of course, you don't know exactly how many primes there are up to limit without actually calculating the primes, but there are prime counting functions that can give you an estimate. It doesn't hurt to reserve a bit more memory than you'll need.

You don't need the vector primes

You actually don't need to create a std::vector<uint64_t> for holding the prime numbers, all the information is already in is_prime. Instead of creating the vector primes, just loop over is_prime and immediately print the numbers to file.

Missing error checking

Writing to a file can fail for various reasons; either it cannot be opened (maybe you don't have write permissions), or something happens while writing to an already opened file (maybe the disk got full, or it's a USB stick and you removed it before your program completed, and so on). You want to avoid your program pretending that everything is fine in that case.

While you could check every single I/O operation, you can just do the error checking at the end. After calling file.close(), check if file.good() == true; if not then something happened previously. Make sure you then write an error message to std::cerr, and return EXIT_FAILURE. Or throw an exception like you did for argc != 3; this will have a similar effect.

Answer 2

• There is no need to #include <cstdlib>. You don't use anything it provides.

• Throwing an exception on invalid command line argument seems drastic. Printing an error message and exiting is much more usual.

• file.open may fail. Check file.is_open() before proceeding.

• No naked loops (aka more functions please). Every loop implements an important algorithm, and therefore deserves a name. For sure, break wheel_sieve into ramp_up (not sure what the right name is), sieve and filter_primes.

  BTW, filter_primes should be a direct application od std::copy_if.

• A sieve loop (while (true) and break) looks dubious. Consider

    while ((square = k*k) < limit)
  

• if is, generally speaking, a performance killer. I trust the compiler to rewrite

    if (++i == 8) {
        i = 0;
    }
  

  into an equivalent if-less form

    i = (i + 1) & 0x07;
  

  but I'd rather write it explicitly.

• Benchmarking is definitely biased. The execution time is dominated by IO (that is, formatting, and writing to the file). Consider timing just the sieving.

Answer 3

uint64_t

The 64-bit unsigned type std::uint64_t is declared in <cstdint>, but this program is missing the necessary #include line.

I guess that this is because your platform incidentally declares this type as a side-effect of some other standard header that you include. While implementations are permitted to transitively include more of the standard library than you ask for, you can't portably depend on which declarations you'll get unless you include the correct (documented) header for each identifier you use.

Another platform-specific aspect of this code is that it uses this type without its namespace qualification. When you include a standard library header, your implementation is allowed (but not required) to additionally declare versions of its identifiers in the global namespace. As you might guess, it's not portable to depend on this behaviour, and our programs should refer to the std namespace type - either in full each time or by using std::uint64_t;, or (probably best) by declaring our own alias such as using integer_type = std::uint64_t; so that there's a single place to change if we want to use a different type for these integers.

(As an aside: there is a header <stdint.h> which is specified to define the fixed-width integer types in the global namespace, and may - optionally - also declare the std-namespace versions. This is intended for use in headers that are also to be included in C programs, and is not recommended for writing C++ code.)

Finally, nothing in this code requires our integer type to be exactly 64 bits (remember that this type is optional in C++, and provided only if the target has a suitable type with no padding). We could still usefully use it on platforms that provide a larger type instead. I would recommend using std::uint_fast64_t, which is defined to be the fastest available type of 64 bits or more, and is available on all targets.

Whatever type is used, think carefully about overflow. For example, square = k * k is risky, and effectively limits the allowable range of k to that of std::uint32_t.