I’m writing a program to count prime numbers from 2 to N using OpenMP. I wrote this code and while testing it, for example, with N = 10 (OMP_NUM_THREADS=2 ./main 10) I get mostly 4-s but sometimes I get 5-s in my output. It seems like there is a race condition that I cannot see and I don’t now where it is. Can you help me to find this condition and get rid of it?

#include <iostream>
#include <omp.h>
#include <vector>

int main(int argc, char *argv[])
{
    if (argc != 2)
        return -1;

    const unsigned long N = std::atoi(argv[1]);

    unsigned counter = 0;
    std::vector<bool> numbers(N + 1, true);
    numbers[0] = numbers[1] = false;

#pragma omp parallel shared(numbers, counter)
    {
#pragma omp single
        {
            for (size_t i = 2; i * i <= N; ++i)
            {
                if (numbers[i] == true)
#pragma omp task
                {
                    for (size_t j = i * i; j <= N; j += i)
                        numbers[j] = false;
                }
            }
#pragma omp taskwait
        }
#pragma omp for reduction(+ 
                          : counter)
        for (size_t i = 1; i <= N; ++i)
            if (numbers[i] == true)
                counter++;
    }

    std::cout << counter << std::endl;

    return 0;
}