Prime Number Generation Part 1: Optimized Trial Division With AVX512 SIMD Intrinsics.

#include <iostream>
#include <fstream>
#include <immintrin.h>

// Returns # of primes up to limit and optionally saves them to a file. Limit must be an integer below 4294967295 (ULONG_MAX) and greater than 3;
static int Primes(uint32_t limit, bool save = false)
{
    uint32_t total = 0;

    // Allocate memory for expected # of primes. Formula from Pierre Dusart
    int end = (limit / log(limit)) * (1 + 1.2762 / log(limit));
    uint32_t* primes = new uint32_t[end];

    // Fill expected array maximum integer value. Needed for one of the conditions for stopping the trial division loop.
    std::fill(primes, primes + end - 1, ULONG_MAX);

    // Seed prime array with first two primes
    primes[0] = 2;
    primes[1] = 3;

    //Create vector of all zeros to compare remainders with.
    __m512i ZeroVector = _mm512_setzero_si512();

    for (uint32_t index = 3; index <= limit; index += 2)
    {
        __mmask16 p = 0;
        uint32_t sq = sqrt(index);

        // Start pointer at prime number 3 (x[1]) since for loop only iterates over odd numbers.
        uint32_t* x = primes + 1;

        // Create vector filled with current index.
        __m512i VectorOfIndex = _mm512_set1_epi32(index);

        // Check for remainders (mod) from 16 primes in parallel for each loop iteration.
        while (x[0] != ULONG_MAX && x[0] <= sq && !(p = _mm512_cmpeq_epi32_mask(_mm512_rem_epu32(VectorOfIndex, _mm512_loadu_epi32((__m512i*)x)), ZeroVector)))
            x += 16;

        // If number is prime add it our prime array to compute further primes.
        if (!p)
        {
            total++;
            primes[total] = index;
        }
    }

    if (save)
    {
        std::ofstream myfile("primes.txt");

        for (int count = 0; count <= total; count++)
            myfile << primes[count] << " ";

        myfile.close();
    }

    //for (int count = 0; count <= total; count++) std::cout << primes[count] << std::endl;

    delete[] primes;
    return total + 1;
}


int main()
{
    uint32_t limit = 2000000000;
    time_t t1 = time(0);

    int total = Primes(limit, false);

    std::cout << "There are " << total << " primes from 0 to " << limit << "." << std::endl;
    std::cout << "Computation completed in " << time(0) - t1 << " seconds." << std::endl;
}