Counting the number of primes in C and C++11

Question

I found this old piece of code that I wrote while I was learning C a while back:

#include <stdio.h>
#include <stdlib.h>
#include <math.h>


int main(int argc, char *argv[]) {
    if (argc < 2) {
        fprintf(stderr, "Usage: %s <N>\n", argv[0]);
        return 1;
    }

    size_t limit = strtoul(argv[1], (char **) NULL, 10);
    size_t i, j, res = 0;
    size_t isqrt = (size_t) sqrt(limit);
    char* numbers = calloc(limit, 1);

    for (i = 2; i <= isqrt; i++) {
        if (numbers[i]) continue;
        for (j = i*i; j < limit; j += i) numbers[j] = 1;
    }

    for (i = 2; i < limit; i++) {
        if (!numbers[i]) ++res;
    }
    printf("%lu\n", res);
    free(numbers);
    return 0;
}

It's a sieve of Erathosthenes. Since now I started learning C++, I decided to rewrite it in C++11:

#include <iostream>
#include <memory>
#include <vector>
#include <string>
#include <algorithm>
#include <math.h>


int main(int argc, char *argv[]) {
    if (argc < 2) {
        std::cerr << "Usage: " << argv[0] << " <N>" << std::endl;
        return 1;
    }

    size_t limit = std::stoul(argv[1]);
    size_t i, j, res = 0;
    size_t isqrt = size_t (sqrt(limit));
    std::unique_ptr<std::vector<bool> > sieve(new std::vector<bool>(limit));

    for (i = 2; i <= isqrt; ++i) {
        if (sieve->at(i)) continue;
        for (j = i*i; j < limit; j += i) sieve->at(j) = 1;
    }

    for_each(sieve->cbegin() + 2, sieve->cend(), [&res](bool i) { res += !i; });
    std::cout << res << std::endl;
    return 0;
}

Is this the most idiomatic way to code in C++? I also tried to replace the first for loop with a for_each and a functor that accepts sieve, but in the end I decided it wasn't worth it and left it like that.

In both versions I am allocating on the heap because the sieve array/vector gets too big for the stack.

Interestingly, the C++ version is faster! I thought they were pretty much identical, but instead the C++ version wins by a large margin:

$ gcc --std=c99 -O3 erat.c -lm -o erat
$ time ./erat 1000000000
50847534
./erat 1000000000  11.50s user 0.41s system 99% cpu 11.927 total
$ g++ --std=c++14 -O2 erat.cpp -o erat2
$ time ./erat2 1000000000
50847534
./erat2 1000000000  8.52s user 0.02s system 99% cpu 8.554 total

Feedback on the C version is also welcomed!

Answer 1

Idiomatic C++ Review

#include <math.h>

C++ has its own version of math header.

#include <cmath>

Difference is all the functions are placed in the std namespace.

Prefer not to use std::endl.

        std::cerr << "Usage: " << argv[0] << " <N>" << std::endl;

Forcing a flush is not normally what you want to do. The streams will automatically flush when required and doing it manually will only make the code less efficient.

Prefer to put one variable per line. Also for loop counters can be declared inside the for statement (see below).

    size_t i, j, res = 0;

Don't dynamical allocate when a local variable will do:

    std::unique_ptr<std::vector<bool> > sieve(new std::vector<bool>(limit));

    // Should be
    std::vector<bool> sieve(limit);

Also std::vector<bool> is special. It is optimized for space not speed. It also has some other quirks that make it undesirable, so prefer std::vector<char> unless you really really need to save space. BUT because of the space saving you can get better caching which may improve time. In this case presumably because of the size of the vector keeping it as bool makes it much more efficient.

see: Comparing std::vector to std::vector

Declare loop variables inside the for:

    for (i = 2; i <= isqrt; ++i) {

    // Usually written like this:
    for(int i = 2; i <= isqrt; ++i) {

The at() method is a checked access to the member. It validates your index is in range (you don't do this in the C code).

        if (sieve->at(i)) continue;

Unless you are using unvalidated user input there is little need to use the at() method, prefer to use operator[].

In C++14 the functions std::begin() and std::end() were introduced to get the iterators for containers. You should prefer to use these as they work with arrays as well as the standard containers.

    for_each(sieve->cbegin() + 2, sieve->cend(), [&res](bool i) { res += !i; });


    // Can be written as:
    auto begin = std::begin(sieve);
    std::advance(begin, 2);
    for_each(begin, std::end(sieve), [&res](bool i) { res += !i; });

Don't bother with a return at the end of main.

    return 0;

Answer 2

I would say that no, it's not idiomatic C++. There's no reason to use dynamic allocation for your vector<bool>, for one obvious example.

As for its being faster, yes, that's expected and fits well with previous tests.

If you care about speed, you could at least consider using operator[] instead of at. I rarely find at very useful, personally.

I'd advise against using std::endl as a general rules. If you just want a new-line, just write a new-line. If you want a new-line followed by a flush (which is what endl actually does) do that explicitly, such as std::cout << res << "\n" << std::flush;.

While res isn't exactly a terrible name, I think something more descriptive (e.g., result) would be better. If (for whatever reason) you're concerned about the length of the name, I'd consider using sum instead.

Answer 3

1. Non-portable format specifier with size_t. Use z

   size_t res = 0;
   ...
   // printf("%lu\n", res);
   printf("%zu\n", res);
   // if z not available
   printf("%llu\n", (unsigned long long) res);
   

2. size_t may be wider than unsigned long

   // size_t limit = strtoul(argv[1], (char **) NULL, 10);
   size_t limit = strtoumax(argv[1], (char **) NULL, 10);
   // or 
   size_t limit = strtoull(argv[1], (char **) NULL, 10);
   

3. size_t isqrt = (size_t) sqrt(limit); is weak code. FP sqrt() can easily have less precision that size_t and sqrt() itself need not return integer values exactly as expected. Better to use an integer sqrt(). Sample follows

   unsigned isqrt(unsigned num) {
     unsigned res = 0;
     // The second-to-top bit is set: 1 << 30 for 32 bits
     unsigned bit = 1u << (sizeof(num) * CHAR_BIT - 2); 
   
     // "bit" starts at the highest power of four <= the argument.
     while (bit > num)
       bit >>= 2;
   
     while (bit > 0) {
       if (num >= res + bit) {
         num -= res + bit;
         res = (res >> 1) + bit;
       } else res >>= 1;
       bit >>= 2;
     }
     return res;
   }
   

4. Avoid magic numbers. Why 1? Using sizeof *pointer_type is correct regardless the type of the pointer numbers.

   // numbers = calloc(limit, 1);
   numbers = calloc(limit, sizeof *numbers);