Can someone help me understand this Sieve of Eratosthenes script? The last couple of lines have puzzled me

Question

Well, that's why I understand the principle of the sieve of Eratosthenes. I tried and basically couldn't write it myself (I wrote a functional prime generator, it's just inefficient or a sieve). I think I have a bigger problem with math, but the programming has mixed me up too.

import numpy as np def primesfrom2to(n): sieve = np.ones(n/3 + (n%6==2), dtype=np.bool) # print(sieve) for n = 10 returns [True, True, True] 

Ok, write off the bat, I'm a little confused. It generates a list of truth values ​​that will be gradually marked as false because they are defined as compound. I think he says that no more than 1/3 of the values ​​less than n will be simple, but what does adding the modolus operation do?

  sieve[0] = False 

marks 1 as false?

  for i in range(int(n**0.5)//3+1): # print(i) for n = 10 returns 0 and 1 on separate lines 

This sets the range to check. No product has a product larger than its square root. What happens to division by 3+1 ?

  if sieve[i]: k=3*i+1|1 #print(k) for n = 10 returns '5' doesn't change till it adds 7 at n = 50 

So, if sieve[i] true (why is simple / not yet marked as compound?), Then 3*i + 1 or 1 ? How it works? It seems to be generating prime numbers that will subsequently be multiplied to remove all subsequent products, but I don't know how to do this.

  sieve[ ((k*k)/3) ::2*k] = False sieve[(k*k+4*k-2*k*(i&1))/3::2*k] = False 

ok, so they both take primes k and mark all further multiples of them? if n=5 does not do it sieve[8.33::10] = False ? And another one like sieve[41.3::10] ? I just do not understand.

 print(np.r_[2,3,((3*np.nonzero(sieve)[0]+1)|1)]) 

Everything is in order and, finally, it actually generates a list of primes. Again, what with multiplication by three? Obviously, there is something basic in this code that I just don't understand. Also, I again do not understand the concept of |1 .

Oh, and just for fun, here is my effective and terribly ineffective attempt.

 import numpy def sieve(num): master = numpy.array([i for i in range(2, num+1)]) run = [] y=2 while y <= ((num+1)**(1/2)): thing = [x for x in range(y, num+1, y) if x > 5 or x == 4] run = run + thing y = y+1 alist = numpy.in1d(master, run, invert = True) blist = (master[alist]) print(blist) 

It took 57 seconds to calculate prime numbers up to 500,000. I did the Euler sum of primes up to 2,000,000 problems.

Answer 1

Here is a simple sump algorithm in numpy:

 import numpy as np def sieve(Nmax): """Calculate prime numbers between 0 and Nmax-1.""" is_prime = np.ones(Nmax, dtype=bool) Pmax = int(np.sqrt(Nmax)+1) is_prime[0] = is_prime[1] = False p = 2 dp = 1 while p <= Pmax: if is_prime[p]: is_prime[p*2::p] = False p += dp dp = 2 primes = np.argwhere(is_prime)[:,0] print(primes) sieve(100) 

If I delete the print statement, it starts for Nmax = 10 ^ 8 in 2.5 seconds on my modest laptop.

This algorithm is a little inefficient, since it should keep the simplicity of even values ​​and multiple values ​​of 3. You can save disk space by filtering them to track the simplicity of n * 6 + 1 and n * 6 + 5, for any integer n> = 1 This saves you two-thirds of the storage space, due to a slightly larger amount of accounting. It seems you were trying to start with a complicated version. Moreover, all accounting activities tend to become expensive if they include a Python interpreter that evaluates if and performs memory management on your lists. Slicing a Python / numpy array is a much more natural way to do this, and it is probably fast enough for your purposes.

Regarding your questions:

• (n% 6 == 2) is logical, will be interpreted as 0 or 1 and will add one more element to prevent the error "in one".
• int(n**0.5)//3+1 should be read as int(int(n**0.5)/3) + 1 . The department goes before adding. The division by three is that you allocate space only for the values ​​6n + 1 and 6n + 5.
• k=3*i+1|1 means k=3*i+1 , add one if it is even (is it bit wise or). The array element 'i' corresponds to a potential prime number (3*i+1)|1 . Therefore, if the array indices are [0, 1, 2, 3, 4, 5, ...] , they represent the numbers [1, 5, 7, 11, 13, 17, ...] .
• Setting sieve elements to False is done separately for elements representing 6n + 1 and 6n + 5, so there are two purposes.
• If you use this in Python 2.7, integer division is always truncated, so 9/2 will result in 4. In Python 3, it would be better to use the // operator for integer divisions, i.e. (k*k)/3 must be (k*k)//3 .

Edit: you can specify the algorithm you asked for: The quickest way to list all primes below N.