FordPrefect
FordPrefect

Reputation: 406

Problem implementing wheels to sieve of Eratosthenes

I am a bit struggling with further optimizing my prime calculating function.

So far I ended up with the sieve of Eratosthenes.

I found on https://primesieve.org/ a hint to further optimize this with the implementation of wheels and a link to this article: ftp://ftp.cs.wisc.edu/pub/techreports/1991/TR1028.pdf

I tried to translate this pseudocode into python, but it's not properly working. I have the feeling that the iterations in step B are not correct. When calculating prime_sieve_fast(100, 3), 91 is not removed. This is logical since the running variables never reach 7*13 or 13*7. What did I get wrong?

def prime_sieve(n):
    prime_list=[0,0]
    for i in range(2,n+1):
        prime_list.append(1)
    for p in range(2,int(n**(1/2))+1):
        for j in range(p**2,n+1,p):
            prime_list[j]=0
    primenumbers=[]
    for i in range(n):
        if prime_list[i+1]==1:
            primenumbers.append(i+1)
    return prime_list,primenumbers    

def prime_sieve_faster(n,n_wheel):
    primes=prime_sieve(100)[1][:n_wheel+1]
    w=wheel(n_wheel,primes[:-1])
    M=len(w)
    prime_list=[1]*(n+1)
    for r in range(M):
        if w[r%M]==0:
            b=0
        else:
            b=1
        for i in range(r,n+1,M):
            prime_list[i]=b
    for i in range(n_wheel):
        prime_list[primes[i]]=1
    prime_list[1]=0
    for p in range(primes[n_wheel],int(n**(1/2))+1):
        print(p)
        step=w[p%M]
        if step==0:
            prime_list[p]=0
        else:
            for f in range(p,p+M+1,step):
                for x in range(p*f,n+1,M*p):
                    prime_list[x]=0
                    print(p,f,x,M*p)
    primenumbers=[]
    for i in range(n):
        if prime_list[i+1]==1:
            primenumbers.append(i+1)
    return prime_list,primenumbers


def wheel(k,primes):
    M=1
    for prime in primes:
        M*=prime
    W=[1]*M
    for prime in primes:
        for x in range(0,M,prime):
            W[x]=0
    W[M-1]=2
    prev=M-1
    for x in range (M-2,0,-1):
        if W[x]!=0:
            W[x]=prev-x
            prev=x
    return W

Upvotes: 1

Views: 232

Answers (0)

Related Questions