Irmen's blog

Priemgetallen

Beetje aan het prutsen geweest met een algoritme dat in Python op een efficiente manier priemgetallen berekent. De beste tot nu toe is deze:

import math

def primes(largest):
    bound=(largest+1)/2 # only use space for the odd numbers
    sieve=[True]*bound
    # sieve index i = odd integer 2i+1
    sieve[0]=False
    sq=int(round(largest**0.5))
    for i in xrange(1,sq/2+1):
        if sieve[i]:
            #r=len(range(i*2*(i+1), bound, 2*i+1))
            r=int(math.ceil((bound-i*2*(i+1))/(2.0*i+1)))
            sieve[i*2*(i+1)::2*i+1]=[False]*r
    result=[2]
    for i in xrange(bound):
        if sieve[i]: result.append(i+i+1)
    return result

Zeef van Eratosthenes, gebruikt alleen geheugen voor de oneven getallen (omdat alle even getallen behalve 2 zeker geen priemgetal zijn). Toch zijn er nog 2 problemen mee:

• neemt bij grote getallen nog behoorlijk veel geheugen in. Zo veel dat hij crasht met een memoryerror bij Spoj.
• hij begint altijd bij 2. Ben er nog niet achter of er wel een handige manier is om vanaf een ander getal te beginnen (niet met dit algoritme in elk geval).

edit: Ik heb wat extra algoritmes toegevoegd: een factorisatie functie, een functie om te checken of een gegeven getal priem is, en een nieuwe functie die de priemgetallen tussen min en max berekent. Deze laatste maakt gebruik van de voorberekende sieve als de range klein genoeg is, en anders gaat hij met factorisatie aan de slag. Het is nog niet voldoende om de Spoj te halen (hij is simpelweg te traag) maar voor niet al te grote getallen volstaat het wel denk ik. Voor meer info zie bijvoorbeeld Sieve en voor een snelle generator in C: Primegen.

Code volgt hieronder...

import math,sys,bisect

def sieveprimes(largest):
    bound=(largest+1)/2 # only use space for the odd numbers
    sieve=[True]*bound
    # sieve index i = odd integer 2i+1
    sieve[0]=False
    sq=int(round(largest**0.5))
    for i in xrange(1,sq/2+1):
        if sieve[i]:
            r=int(math.ceil((bound-i*2*(i+1))/(2.0*i+1)))
            sieve[i*2*(i+1)::2*i+1]=[False]*r
    result=[2]
    for i in xrange(bound):
        if sieve[i]: result.append(i+i+1)
    return result

# prepare the base array of primes
primes=sieveprimes(1000000)   # ~ 78k primes

def factorize(n):
    result=[]
    sq=int(round(n**0.5))
    if sq>primes[-1]:
        raise ArithmeticError("sieve too small")
    for f in primes:
        if f>sq or n==1:
            break
        while n%f==0:
            yield f
            n=n/f
    if n!=1:
        yield n

def is_prime(n):
    # first try some obvious division tests
    if n%2==0 or n%3==0 or n%5==0 or n%7==0 or n%11==0:
        return n in (2,3,5,7,11)
    if n<=primes[-1]:
        # try to find it in the sieve
        i=bisect.bisect_left(primes,n)
        return primes[i]==n
    else:
        # check if it has any prime factors other than itself
        return factorize(n).next()==n

def calcprimes(min,max):
    if max>primes[-1]:
        # use factorization (of only the odd numbers)
        if min==1 or min==2:
            yield 2
            min=3
        if min%2==0:
            min+=1
        if max%2==0:
            max-=1
        for n in xrange(min,max+1,2):
            if is_prime(n):
                yield n
    else:
        # just return part of the sieve
        i=bisect.bisect_left(primes, min)
        pl=len(primes)
        while primes[i]<=max:
            yield primes[i]
            i+=1
            if i>=pl:
                break

if __name__=="__main__":
    count=int(sys.stdin.readline())
    for i in range(count):
        min,max=sys.stdin.readline().split()
        min,max=int(min),int(max)
        for p in calcprimes(min,max):
            print p
        print