Algorithm is incorrect
The algorithm you are using is incorrect because you do not find the proper range of primes to use for your segmented sieve. You are currently finding and using primes in the range 0..(n-m+1)
. But you actually need to be finding primes in the range 0..sqrt(n)
. Here is a trivial test case that demonstrates the problem:
1
9 9
For this test case, you find primes in the range 0..1
, which finds no primes, which leads to an incorrect result where you print 9
as a prime. The correct range of primes should be 0..3
for this test case.
Find primes once
Given the maximum \$n\$ of \$10^9\$, you should do a single sieve to find all the primes under \$\sqrt {10^9}\$. There should be 3401 primes under 31623. In fact, you could even hardcode this list of primes into your program to avoid performing the sieve at runtime. It should only take 14KB to hardcode that list (7KB if you store them as short
).
Optimizing the segmented sieve
Let's take a look at your main loop:
for (const auto &prime : primes) {
for (int i = (m / prime) * prime; i <= n; i += prime) {
if (i >= m && i != prime) {
is_prime[i - m] = false;
}
}
}
There are several things I notice that are suboptimal:
The loop increments by prime
. You could achieve 2x speeds if you only handled odd numbers and incremented by 2*prime
instead. The only even prime is 2 and you can check for that as a special case outside the sieve.
You have a if (i >= m)
check which could be removed if you just started the loop from the right starting value.
You have a if (i != prime)
check which could be removed if you just started the loop from the right starting value.
You do a subtraction i-m
on every iteration.
So what is the correct starting value? You were close in your (m / prime) * prime)
expression. The problem is that your expression could result in a value smaller than m
. You should just round up like this:
int start = ((m + prime - 1) / prime) * prime;
This gets rid of the i >= m
check. But also, start
should be at least prime * prime
, because any multiple of prime
less than that will already have been handled by a lower prime.
start = max(start, prime * prime);
This gets rid of the case where you might reach i == prime
. Thirdly, start
should be odd, because if you are doing #1 above, you should start at an odd index and increment by an even increment.
if ((start & 1) == 0)
start += prime;
Lastly, by subtracting m
from start and making your loop go until n-m
, you can get rid of the subtraction inside the loop.
start -= m;
Rewrite
Here is a rewrite of your whole program using the above ideas:
#include <cmath>
#include <cstdio>
#include <vector>
using namespace std;
#define MAX 1000000000
typedef vector<bool> VB;
typedef vector<int> VI;
// This holds the primes from 3..sqrt(MAX)
VI primes;
// Finds the primes from 3..n and stores them in the global "primes".
static void findPrimes(int n)
{
int i;
int sqrtn = sqrt(n);
VB is_prime(n+1, true);
for (i = 3; i <= sqrtn; i += 2) {
if (is_prime[i]) {
int increment = i+i;
primes.push_back(i);
for (int j = i * i; j <= n; j += increment) {
is_prime[j] = false;
}
}
}
for (; i <= n; i += 2) {
if (is_prime[i])
primes.push_back(i);
}
}
// Finds the primes in the range [m..n] and prints them.
void segmented_sieve(int m, int n)
{
int range = n - m;
VB is_prime(range + 1, true);
for (int prime : primes) {
int increment = prime + prime;
// Start should be the first multiple of prime >= m, or prime^2,
// whichever is larger. Also make sure start is odd.
int start = max(((m + prime - 1) / prime) * prime, prime * prime);
if ((start & 1) == 0)
start += prime;
// Adjust to range so we don't have to subtract m in the loop.
start -= m;
for (int i = start; i <= range; i += increment) {
is_prime[i] = false;
}
}
if (m <= 2 && n >= 2)
printf("2\n");
for (int k = (m & 1) ? 0 : 1; k <= range; k += 2) {
if (is_prime[k]) {
printf("%d\n", m + k);
}
}
}
int main()
{
int t;
scanf("%d", &t);
findPrimes(sqrt(MAX));
while (t--) {
int m, n;
scanf("%d%d", &m, &n);
if (m == 1) m++;
segmented_sieve(m, n);
}
return 0;
}
using namespace std
. \$\endgroup\$ – Cody Gray Jan 24 '17 at 15:58