Find primes using “naive” algorithm

I know I can use Sieves to make this faster, but I wondered if there's anything I'm missing with this implementation of the "naive" algorithm.

On my PC running the command (first million primes):

g++ -std=c++11 -Wall -Wextra -Ofast primes.cpp -o primes && time ./primes 15485863

this takes approximately 9.3 secs.

#include <climits>
#include <cstdio>
#include <string>
#include <vector>

typedef unsigned long long ull;

int main(int argc, char *argv[])
{
ull limit = ULLONG_MAX;
if (argc > 1) limit = std::stoi(argv[1]);

std::vector<ull> primes;
primes.push_back(2);
printf("2 is prime\n");
primes.push_back(3);
printf("3 is prime\n");
primes.push_back(5);
printf("5 is prime\n");
int inc = 2;
for (ull i = 7; i <= limit; i += inc) {
bool isprime = true;
for (ull j = 0; primes[j] * primes[j] <= i && j < primes.size(); j++) {
if (i % primes[j] == 0) {
isprime = false;
break;
}
}

if (isprime) {
printf("%llu is prime\n", i);
primes.push_back(i);
}
inc == 2 ? inc = 4 : inc = 2;
}
}

• Don't be sure the sieve is faster. This loop runs in the L1 cache, any sieve algorithm runs off main memory. A while back I ran a test and for the number in question the "naive" approach blew the sieve out of the water. – Loren Pechtel Oct 14 '14 at 19:47
• What was the number? 42? – Andrew Lazarus Oct 15 '14 at 1:50

A quick detail

You've written inc == 2 ? inc = 4 : inc = 2; : this is not the ""right"" way to use the "ternary" operator. Don't get me wrong, it will do what you expect it to do. However, the whole point of the operator is to return a value. In your case, it is no different that : if (inc == 2) inc = 4; else inc = 2;.

This calls for a ternary operator but to write it in a simple way : inc = (inc == 2) ? 4 : 2;

• Or just: inc = -inc + 6. Probably slightly faster too. – agtoever Oct 14 '14 at 18:55
• @agtoever faster, but very unclear as to see what it's supposed to do for future readers (personally, I'd wonder where the '6' constant comes from and why it's hardcoded...) – Sanchises Oct 14 '14 at 19:11
• @sanchises OP's using optimization flags for speed. When it comes to performance, I'd rather write a paragraph of comments than have code that's easier for a human to read than a machine. If this were anywhere else where milliseconds don't matter, I'd let the code do the talking in a more legible way. – Corey Ogburn Oct 14 '14 at 21:54
• @sanchises 6 = 2*3. and next comes 30=2*3*5. – Will Ness Oct 20 '14 at 20:43

Code style and good practices:

On a C++ program you should give preference to the standard C++ library. I know the syntax of printf and friends might be tempting some times, but std::cout is the proper C++ way of printing to the standard output.

You should also prefer to use std::numeric_limits over the macros of <climits> (ULLONG_MAX).

This block could be rewritten to use both numeric_limits and const:

ull limit = ULLONG_MAX;
if (argc > 1) limit = std::stoi(argv[1]);


Using the ternary operator:

const ull limit = (argc > 1) ? std::stoi(argv[1]) : std::numeric_limits<ull>::max();


Try to take advantage of const whenever you can.

Possible otimization:

Some time is probably being spent reallocating the vector primes when it has to increase its capacity. If you could pre-allocate some estimated space with vector::reserve() before entering the loops, it could potentially give a performance boost to your program.

Miscellaneous:

Judging by your use of std::stoi, your compiler is probably C++11 compliant. You can then try the new using alias to replace the older style typedef:

using ull = unsigned long long;

• I usually do stick to the stl, in fact i was using cout until just before i submitted this. But I've read before that cout has issues with speed? Or is that just endl vs \n ? – LordAro Oct 14 '14 at 14:46
• I've not been convinced by the new using vs typedef - it just seems unecessary to me, and in some ways less clear – LordAro Oct 14 '14 at 14:48
• @LordAro, funny, the using alias always seemed more natural to me, as it is written like dest = source vs set source dest in the typedef version. – glampert Oct 14 '14 at 14:52
• I didn't test extensively, but i seemed to get about 0.5secs increase in speed from switching to printf (including endls) but i'll try couts without endl. When you put it like that, the using seems so much nicer! Maybe I'm just not used to it – LordAro Oct 14 '14 at 14:56
• If you're serious about optimizing I/O, you should do std::ios_base::sync_with_stdio(false); to turn off the synchronization between cout and STDOUT (and the other streams). By default, synchronization is turned on, which can seriously impact the performance of writing to cout. (but with the benefit of there being no surprises if you mix cout and printf) – user14393 Oct 14 '14 at 17:39

In terms of performance, the first thing to note is that writing to a terminal is slow. time ./primes gives 7.995 seconds on my machine, but time ./primes > primes.out gives 5.107 seconds, already 36% faster.

The next thing is that 64-bit integer division/modulo operations are horribly slow, something like 50 cycles each compared to around 5 for a 32-bit division. If you swap unsigned long long for unsigned int or std::uint32_t (from <cstdint>), the runtime drops to 2.383 seconds.

By the way, this line:

for (ull j = 0; primes[j] * primes[j] <= i && j < primes.size(); j++) {


is particularly scary because the primes[j] * primes[j] <= i check happens before the j < primes.size() check, potentially accessing outside the vector. But it turns out not to cause undefined behaviour because the first inequality will always become false before the second.

• I'm curious as to why 64 bit division is so slow, especially on 64 bit CPUs. Could you explain? – LordAro Nov 28 '14 at 12:34
• I'm not really a hardware/electronics guy but basically, circuits for computing division are much larger and deeper than other arithmetic operations and thus have many more gate delays. – Tavian Barnes Nov 28 '14 at 14:52

A simple way to make this faster: You check primes[j] * primes[j] <= i all the time, over and over again. In reality, let's say i = 1001, then you will check the primes up to 31. And you will check the primes up to 31 until i = 37*37 = 1369. So have a local variable n that indicates how many primes you might be checking at most, and if i == primes[n] * primes[n] then increase n by 1.

Also, start the loop with j = 2 since i is never divisible by 2 or 3.

• primes [i] * primes [i] <= i er, no? unless you mean primes[j], but j is changing in each value of the loop – LordAro Oct 14 '14 at 15:53
• The point is that the possible number of iterations is usually the same, and can be found with a single test before the iteration starts, so we get a loop "for (j = 2; j < count; ++j) ... – gnasher729 Oct 14 '14 at 17:46
• Given that sqrt() is much slower than integer multiplication, have you identified the crossover point where it becomes faster (on your system)? – Toby Speight Sep 11 '18 at 13:54

Though many hinted towards it but haven't explicitly mentioned it.
Pre-calculate the square root of i, say sqrt_i and run this loop:

for (ull j = 0; primes[j] <= sqrt_i && j < primes.size(); j++)


You really should add more comments to your code. I'd also write out increment rather than abbreviating inc. After reading the other answer, I was wondering what inc did. It took me a little bit to find where it was actually used. In general, any time you do something clever, you should comment.

// 2 and 3 are the only primes that do not fit the pattern 6 * n +/- 1
// so we don't have to check every number, we can start with 7 (6*1+1)
// increment by 4 to 11 (6*2-1); increment by 2 to 13 (6*2+1);
// so on and so forth.  We alternate incrementing by 2 and 4 every time.
// we update the increment during the loop, so setting to 2 now
// means that we will set to 4 for the first iteration
int increment = 2;
for (ull i = 7; i <= limit; i += increment) {
// 6 - 2 = 4 and 6 - 4 = 2, so subtracting increment from 6
// alternates between 2 and 4, as we want
increment = 6 - increment;


Your default limit won't work. The variable will never be greater than its maximum. This is especially bad because you are using <= and an increment greater than 1. It might barely work when incrementing by 1 with <. As is, you need a maximum limit of something like ULLONG_MAX - 7.

• Your concerns about the limit are valid but it's not like ULLONG_MAX could ever be reached in this way anyway. – Tavian Barnes Oct 15 '14 at 15:51