## DEV Community

Antidisestablishmentarianism

Posted on • Updated on

# 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 = 2;
primes = 3;

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

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

// Start pointer at prime number 3 (x) 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 != ULONG_MAX && x <= 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;
}

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;
}

``````