Non-Sieve Prime Number generator

#include <iostream>
#include <vector>
#include <cmath>
using namespace std;
int main()
{
int l=0;
bool isprime;
vector<int>primes;
primes.push_back(2);
for (int nur = 3;true;nur+=2)
{

isprime=true;
double ilimit= sqrt(nur);
for (int primecount=0 ;primes[primecount]<=ilimit;primecount++)
{
if (nur % primes[primecount] == 0)
{
isprime = false;
break;
}

}

if (isprime)
{
cout<<nur<<endl;
primes.push_back(nur);
// l++;
}

if (nur<0)
break;

}
cout<<"Finished"<<endl;
//cout<<primes[l];

return 0;
}

Is there anyway i can optimize this further without using a sieve?

• You could save on sqrt() calls by computing the next number when that actually has to be computed: if (nur >= next_limit) { ilimit = sqrt(nur); next_limit = (ilimit+1)*(ilimit+1); } (or something similar; been a while since I used C++). Jun 13 '16 at 21:23
• I dont really get what your saying. Could you please Elaborate? Jun 13 '16 at 21:31
• sqrt(15) = 3.87, sqrt(16) = 4.00, sqrt(17) = 4.12, sqrt(18) = 4.24, sqrt(19) = 4.36, sqrt(20) = 4.47, sqrt(21) = 4.58, sqrt(22) = 4.69, sqrt(23) = 4.80, sqrt(24) = 4.90, sqrt(25) = 5.00, sqrt(26) = 5.10 ---> Instead of setting ilimit = sqrt(nur) for every nur, you could set it only when needed (in this case, when nur = 16 and when nur = 25, but not for the rest of the numbers that I've listed). Jun 13 '16 at 21:39
• Ok anything else? Jun 13 '16 at 21:40
• You might find this interesting: codereview.stackexchange.com/a/90503/36120 Jun 15 '16 at 12:40

Your code includes the optimization of evaluating only odd-valued candidates (by iterating with a step size of 2) because the even-valued candidates are known to contain a factor of 2, and are therefore always composite.

This candidate-skipping optimization can be generalized to include prime factors greater than 2, by varying the step size in a certain way.

I made a version of your code which skips candidates (values of nur) containing factors of 2 and/or 3. This is done fairly cheaply by alternating step size between 2 and 4.

Extending this optimization to even greater primes gives diminishing returns for two reasons: The step size adjustment calculation becomes more complex (and slower), eventually costing more time than what is saved by skipping candidates; and higher prime factors exist more sparsely than low ones in the candidate population. Skipping 2 alone catches 1/2 of the composite candidates, and skipping 3 alone would catch 1/3. But if used together, half of the composites 3 would skip were already skipped by the 2 skip, so skipping 3 just catches an additional 1/6 of the composite candidates. Including a 5 skip would only eliminate an additional 1/30 of the composite candidates, and a 7 skip 1/210.

I experimented with replacing the floating-point sqrt call with an integer approximation, but it turns out that sqrt is very fast (at least on my machine). And then I gave up entirely after seeing Vedran Šego's optimization, which is brilliant, especially because as nur gets bigger, the gaps between calls to sqrt become larger, to the point where it doesn't really matter much how slow your sqrt calculation is. So, I had to add this to my example, thank you Vedran Šego.

One other thing, the nur loop must terminate. If nur is signed (like int), it will overflow and wrap around to a negative value. This is a completely pathological situation -- calling sqrt with a negative number, plus other things, depedning on optimizations employed.

Even for unsigned nur, duplicate primes will be added to the primes vector until the process crashes when it runs out of memory.

So my example terminates before nur overflows. I also defined a prime_t type, to make it easy to switch the integer type used for calculating and storing primes. I use unsigned int in my example, because it gives results a larger maximum value at no extra cost compared to int. Try unsigned long long if you really want to see the program continue for much longer than you will be willing to wait.

#include <iostream>
#include <vector>
#include <cmath>

typedef unsigned int prime_t;
// max_prime_t is the maximum value which may be held in a variable of type prime_t
// this assumes prime_t is an integer type (signed or unsigned)
// It dupicates information available from std::numeric_limits
const prime_t max_prime_t =
prime_t(0)-prime_t(1) > prime_t(0) ?
~prime_t(0) :
~prime_t(0) ^ (prime_t(1) << (sizeof(prime_t)*8-1));

int main() {
std::vector<prime_t> primes = { 2, 3 };
prime_t ilimit = 0;
prime_t ilimit_limit = 0;

static const prime_t max_nur = max_prime_t - 4;

// "step" alternates between 2 and 4,
// to skip values containing factors of 2 and/or 3
for(prime_t nur=5, step=4; nur<=max_nur; nur += (step ^= 6)) {
if(nur >= ilimit_limit) {
// Vedran Sego's optimization
ilimit = static_cast<prime_t>(std::sqrt(nur));
ilimit_limit = (ilimit + 1) * (ilimit + 1);
}
for(auto one_prime : primes) {
if(one_prime > ilimit) {
std::cout << nur << '\n';
primes.push_back(nur);
break;
}
if(nur % one_prime == 0) {
break;
}
}
}
}
• Thanks, I'll remove 5 from the skipping claim and try to figure out what made me think I was skipping 5 when I wrote the old Javascript code this was taken from -- maybe I missed something. Jun 15 '16 at 10:51
• Forget integer approximations to sqrt: it's almost completely unnecessary, and if you initialise ilimit and ilimit_limit correctly it becomes completely unnecessary. With the exception of the first time round the loop (when it sets ilimit = 2; ilimit_limit = 9;) the sqrt always evaluates to ilimit + 1. Jun 15 '16 at 13:45

Optimization

Probably the single biggest optimization on this code has nothing to do with the prime computation itself at all.

For the sake of testing, I changed it to generate primes up to 10000000, by modifying the outer loop to become:

for (int nur = 3; nur<10000000; nur += 2)

With the rest of the code as originally written, this ran in ~2300 ms. I then changed:

cout << nur << endl;

to:

cout << nur << "\n";

...and this reduced the time to ~1200 ms (approximately doubled the speed).

Another fairly minor change was to change from computing the square root, to computing the square of the trial divisor, and seeing if that was less than the number being tested. This doesn't make a huge difference in speed, but does help a bit.

The obvious step after that (at least to me) was to rewrite the code to encapsulate the logic into a class. I didn't expect this to give a huge increase in speed by itself, but did expect to get cleaner code that might be a little easier to optimize further. That gave this:

#include <iostream>
#include <vector>
#include <cmath>
#include <algorithm>
#include <iterator>

class isprime {
std::vector<int> primes;
int n;

bool check(int n) {
auto div = primes.begin();
for ( ; *div * *div <= n; ++div)
if (n % *div == 0)
return false;
primes.push_back(n);
return true;
}

public:
isprime() : primes{ 2 }, n{ 3 } {}

int operator()() {
while (!check(n))
n += 2;
auto ret = n;
n += 2;
return ret;
}
};

int main() {
isprime p;

int prime;
while ((prime = p()) < 10000000)
std::cout << prime << '\n';
}

Somewhat surprisingly, this actually gave another fairly healthy improvement in speed, so it ran in about 900 ms. I decided that was enough, at least for now.

As far as reviewing the code otherwise goes:

indentation

Right now your indentation renders the code considerably less readable than it could be.

vertical whitespace

Much like your indentation, you have a fair (unfair?) amount of seemingly random vertical white space (i.e., blank lines) that don't seem to serve an purpose, and reduce visibility. I expect a blank line to be a little like a separation between paragraphs in text--a separation that signals a small change of subject, or something on that order.

Implement functions as functions

You don't need to stuff everything into main. Some people do this in a mistaken attempt at optimization.

I wouldn't be afraid to implement the code to test whether a number is prime as a function instead of inline code. Compilers are usually pretty decent at figuring out where you'll get significant benefit from having the code generated inline, and doing so where it'll help (though MS does support __forceline for those cases where it is possible, and for whatever reason the compiler chooses not to).

avoid std::endl

As shown above, std::end can and does lead to significant slow-downs in some cases, especially when substantial amounts of output are being produced. If you just want a new-line, just write a new-line.

Don't be afraid of classes either

As shown here, encapsulating the code into a class didn't slow it down at all--in fact, it improved speed by 25% or so.

The obvious next step would be to make the code multi-threaded. One obvious way to do this would be to start by using one thread to compute the primes up to the square root of the maximum you care about. From there, split the job into N threads, each of which gets a copy of the primes found so far, and uses it to check every Nth candidate for primacy.

Perhaps less obviously (but only slightly so) the same process can be used recursively--that is, to compute the primes up to the square root of the limit, compute the primes up to the square root of that limit in a single thread, then split the job out into multiple threads from there (continue the same basic idea until you get such a small group of primes to find that splitting it among threads is no longer worthwhile).

Of course, after this you're left with a number of separate collections of primes. You'll need to merge those together to get a single collection that's all in order.

This potentially gives nearly perfect scaling, so with N cores it can run about N times faster than the previous (but even a really simple sieve is still probably going to be faster).

This is not really a review but rather a benchmark of the original code, Christopher's, Jerry's code and my take on the code for a comparison.

My approach is to compute and store the square of each prime as I go as opposed to calculating the square every time as Jerry's code does. I also only check numbers of the form $6k\pm 1$ which saves some time. Even so it is only 4% faster and I'm leaning towards calling it statistically insignificant.

I don't think we're going to get much faster with simple trial division without applying more advanced primality tests.

#include <iostream>
#include <vector>
#include <cmath>
#include <iterator>
using namespace std;

const auto MAX_PRIMES = 50000000;

#if 0
// Original: 10.727s
int main(){
bool isprime;
vector<int>primes;
primes.push_back(2);
for (int nur = 3;nur<MAX_PRIMES;nur+=2){

isprime=true;
double ilimit= sqrt(nur);
for (int primecount=0 ;primes[primecount]<=ilimit;primecount++){
if (nur % primes[primecount] == 0){
isprime = false;
break;
}
}

if (isprime){
cout<<nur<<endl;
primes.push_back(nur);
}

if (nur<0)
break;
}
cout<<"Finished"<<endl;
return 0;
}
#endif

#if 1
// Mine: 6.575
// value is 6k +-1. hence 2 and 3 are not factors of value,
// skipp them with p=2
for (size_t p = 2; primesAndSquares[p].second <= value; ++p){
if (value % primesAndSquares[p].first == 0){
return false;
}
}
primesAndSquares.emplace_back(value, value*value);
cout<<value<<'\n';
return true;
}

int main(){
vector<pair<int,int>> primesAndSquares; // Store the prime and its square
primesAndSquares.reserve(1000);
primesAndSquares.emplace_back(2, 4);
primesAndSquares.emplace_back(3, 9);
primesAndSquares.emplace_back(5, 25);
primesAndSquares.emplace_back(7, 49); // Needed so that we can start with p=2 above

for (int nur = 12; nur < MAX_PRIMES; nur += 6){
}

cout<<"Finished\n";
return 0;
}

#endif

#if 0
// Jerry 6.803
class isprime {
std::vector<int> primes;
int n;

bool check(int n) {
auto div = primes.begin();
for ( ; *div * *div <= n; ++div)
if (n % *div == 0)
return false;
primes.push_back(n);
return true;
}

public:
isprime() : primes{ 2 }, n{ 3 } {}

int operator()() {
while (!check(n))
n += 2;
auto ret = n;
n += 2;
return ret;
}
};

int main() {
isprime p;

int prime;
while ((prime = p()) < MAX_PRIMES)
std::cout << prime << '\n';

cout<<"Finished\n";
}

#endif

#if 0
// Christopher: 9.376 s
typedef unsigned int prime_t;
// max_prime_t is the maximum value which may be held in a variable of type prime_t
// this assumes prime_t is an integer type (signed or unsigned)
// It dupicates information available from std::numeric_limits
const prime_t max_prime_t =
prime_t(0)-prime_t(1) > prime_t(0) ?
~prime_t(0) :
~prime_t(0) ^ (prime_t(1) << (sizeof(prime_t)*8-1));

int main() {
std::vector<prime_t> primes = { 2, 3 };
prime_t ilimit = 0;
prime_t ilimit_limit = 0;

static const prime_t max_nur = MAX_PRIMES;//max_prime_t - 4;

// "step" alternates between 2 and 4,
// to skip values containing factors of 2 and/or 3
for(prime_t nur=5, step=4; nur<=max_nur; nur += (step ^= 6)) {
if(nur >= ilimit_limit) {
// Vedran Sego's optimization
ilimit = static_cast<prime_t>(std::sqrt(nur));
ilimit_limit = (ilimit + 1) * (ilimit + 1);
}
for(auto one_prime : primes) {
if(one_prime > ilimit) {
std::cout << nur << '\n';
primes.push_back(nur);
break;
}
if(nur % one_prime == 0) {
break;
}
}
}
}

#endif

Benchmarks ran with g++ --std=c++11 -O3 primes.cc -o primes.exe && time ./primes.exe > foo under mingw64 on Windows 7, i5-3570K (3.4GHz). Compiler was gcc 5.4.0.

• I'm interested in the benchmark results without including the output each time a prime is found. I/O is a big bottleneck and is subject to variability in its blocking duration and just shouldn't be included for benchmarking something unrelated to I/O. If each test performs the same output operation and many iterations were performed, at best it just adds the same big chunk of time to each result. This should tend to make the results closer to each other than a more proportional test of the performance of each algorithm would produce. Jun 15 '16 at 21:34
• Oh, one more thing, if you remove the output from inside the loop, you will need to make sure the calculations still have a side-effect, because the slim possibility exists that the compiler will apply an optimization which skips something which should be reflected in the results. I think printing out the primes in the vector (or some kind of hash derived from them), after the benchmark timer has stopped, should suffice. Jun 15 '16 at 21:48
• @Christopher the amount of output is linear in the number of values checked while the algorithms are O(n^2). The running time of the algorithm dominates as n is large. I also verified this by printing the sum of the primes only and the time didn't change in any significant way. Jun 16 '16 at 5:18