3

Sieve of Eratosthenes on the segment: Sometimes you need to find all the primes that are in the range [L...R] and not in [1...N], where R is a large number.

Conditions: You are allowed to create an array of integers with size (R−L+1).

Implementation:

bool isPrime[r - l + 1]; //filled by true
for (long long i = 2; i * i <= r; ++i) {
    for (long long j = max(i * i, (l + (i - 1)) / i  * i); j <= r; j += i) {
        isPrime[j - l] = false;
    }
}
for (long long i = max(l, 2); i <= r; ++i) {
    if (isPrime[i - l]) {
        //then i is prime
    }
}

What is the logic behind setting the lower limit of 'j' in second for loop??

Thanks in advance!!

Will Ness
  • 70,110
  • 9
  • 98
  • 181
Anand Singh
  • 87
  • 10
  • `((n+p-1) // p ) * p` is the smallest natural number that is 1. a multiple of `p` and 2. not less than `n`. You can subtract the `l` from it just once, at `j`'s initialization. – Will Ness Jul 20 '17 at 13:26

1 Answers1

2

Think about what we want to find. Ignore the i*i part. We have only (L + (i - 1)) / i * i) to consider. (I wrote the L capital since l and 1 look quite similar)

What should it be? Obviously it should be the smallest number within L..R that is divisible by i. That's when we want to start to sieve out.

The last part of the formula, / i * i finds the next lower number that is divisible by i by using the properties of integer division.

Example: 35 div 4 * 4 = 8 * 4 = 32, 32 is the highest number that is (equal or) lower than 35 which is divisible by 4.

The L is where we want to start, obviously, and the + (i-1) makes sure that we don't find the highest number equal or lower than but the smallest number equal or bigger than L that is divisible by i.

Example: (459 + (4-1)) div 4 * 4 = 462 div 4 * 4 = 115 * 4 = 460.

460 >= 459, 460 | 4, smallest number with that property

(the max( i*i, ...) is only so that i is not sieved out itself if it is within L..R, I think, although I wonder why it's not 2 * i)

For reasons of readability, I'd made this an inline function next_divisible(number, divisor) or the like. And I'd make it clear that integer division is used. If not, somebody clever might change it to regular division, with which it wouldn't work.

Also, I strongly recommend to wrap the array. It is not obvious to the outside that the property for a number X is stored at position X - L. Something like a class RangedArray that does that shift for you, allowing you a direct input of X instead of X - L, could easily take the responsibility. If you don't do that, at least make it a vector, outside of a innermost class, you shouldn't use raw arrays in C++.

Aziuth
  • 3,652
  • 3
  • 18
  • 36
  • I understood the logic, thank you......but in your last paragraph you said "It is not obvious to the outside that the property for a number X is stored at position X - L"....but I think it always works. – Anand Singh Jul 20 '17 at 10:06
  • @AnandSingh The array could be given to some other place which does only know that it represents prime number but not that it has the shift. Granted, this is a rare case, but in C++, things should be responsible for themselves. Technically, it works at is is, this is just about some standards that one should take on. – Aziuth Jul 20 '17 at 10:09
  • 2
    it's `max( ... i*i ...)` and not `... 2*i ...`, because such multiples of `i` (i.e. `k*i` where `k < i`) will be *already eliminated* as multiples of those `k`'s, `k = 2,3..i-1`. because we process the `i`'s in increasing order, so `k < i` is already processed, before `i` because a smaller `i` was equal to that `k` at some point, previously, already, in the loop. – Will Ness Jul 20 '17 at 13:28
  • 1 more doubt....Why do we have max( ... i*i ....). Can't we just start with j=(l+(i-1))/i *i . – Anand Singh Jul 20 '17 at 16:12
  • As explained, i could be within L..R and also be a prime number. In this case, it would falsely mark itself as not being a prime. Read the comment by Will Ness why it's i^2 in particular. – Aziuth Jul 20 '17 at 16:20
  • @Aziuth yes; and even if that `j` is above `i`, it would be inefficient to start there. we only need to start from `i*i`, or *higher* (inside the segment). – Will Ness Jul 23 '17 at 17:11