Comment by gjm11

sieve[i] is nonzero precisely when 2i+1 is an odd prime or i=0 (because it's initialized to all-1 and nothing ever zeros the first element).

So np.nonzero(sieve) is all the indices for which i=0 or 2i+1 is an odd prime.

Except that the behaviour of np.nonzero is a bit weird: it returns a tuple of arrays, one for each axis. (So, e.g., if you feed it a 2d array it gives you (a,b) such that the nonzero elements are at positions (a[0],b[0]), (a[1],b[1]), etc.) That's kinda reasonable behaviour, but it means that for the simple case of a 1-d array you get a 1-element tuple.

So np.nonzero(sieve)[0] is actually all the indices i for which i=0 or 2i+1 is an odd prime, and now you know why I didn't write "all the indices i for which ..." above.

The case i=0 would correspond to 1 being prime, which it isn't considered to be. So we throw that one out.

So np.nonzero(sieve)[0][1::] is all the indices i for which 2i+1 is an odd prime.

And so 2*np.nonzero(sieve)[0][1::]+1 is all the 2i+1 for which 2i+1 is an odd prime; that is, all the odd primes.

I don't think I've encountered np.r_ before, but it's pretty obvious what it's doing: prepending 2 to that array, so now it's "2 followed by all the odd primes" or, more simply, all the primes.

(Obviously, in the above phrases like "all the primes" mean "all the primes up to the limit defined by the size of the array".)