All right, so I get the principle behind the Sieve of Eratosthenes. I tried, and largely failed, to write one myself (I wrote a functional prime number generator, it just wasn't efficient or a sieve). I think I'm more having an issue understanding the math, but the programming has me mixed 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's generating a list of truth values that will progressively be marked false as they are determined to be composite. I think it's saying no more than 1/3 of the values less than n will be prime, but what does adding a modolus operation achieve?
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 is setting the range to check. No number has a product greater than its square root. Whats up with the dividing 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
Ok so, if sieve[i]
is true (so prime / not yet marked as composite?) then 3*i + 1 or 1
? How does this work exactly? It seems to be generating primes that it will be multiplied later to remove all the subsequent products, but I don't know how.
sieve[ ((k*k)/3) ::2*k] = False
sieve[(k*k+4*k-2*k*(i&1))/3::2*k] = False
ok so both of these are taking the primes k
and marking all further multiples of them? if n=5
wouldn't that make it sieve[8.33::10] = False
? And the other one like sieve[41.3::10]
? I just don't get it.
print(np.r_[2,3,((3*np.nonzero(sieve)[0]+1)|1)])
All right, and finally actually generating the list of primes. Again, what's up with multiplying it by three? Obviously there is something fundamental about 3 in this code that I just don't get. Also I'm failing to grasp the |1
concept again.
Oh and just for fun here's my effective and horribly inefficient 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 57s to calculate the primes up to 500,000. I was doing the Euler sum of primes up to 2,000,000 problem.