소수 구하기.(에라토스테네스의 체, 앳킨의 체, 등등)

소수를 구할때 프로그래밍을 처음 하는 사람이 아니고서는 에라토스테네스의 체를 쓴다.

쉬운 개념이므로 링크만 걸고 코드만 올려놓는다.

https://en.wikipedia.org/wiki/Sieve_of_Eratosthenes

# -*- encoding: cp949 -*-
def prime_number(max_n):
    p = [0 for i in xrange(max_n+1)]
    for i in xrange(2,max_n+1):
        for j in xrange(i*2,max_n+1,i):
            p[j] = 1
    prime = [i for i in xrange(2,max_n+1) if p[i]==0]
    return prime

print prime_number(100) #2~100까지 소수 목록을 리스트로 반환. 100 포함한다.

그런데 문제를 풀다 찾아보니 앳킨의 체라는 더 빠른게 있다. 얼마나 차이가 있겠어.. 했는데

엄청나게 빠르다. 구글에서 코드를 가져왔다. (원리는 찾아봐도 이해를 못하겠다;;)

https://programmingpraxis.com/2010/02/19/sieve-of-atkin-improved/

에라토스테네스의 체와 앳킨의 체 둘다 2000000까지 소수를 찾기로 하고 시간 비교를 해봤을 때

압도적으로 시간 차이가 나는걸 볼 수 있다.

# -*- encoding: cp949 -*-
import math, time
from math import *
def atkins13(limit=2000000): #앳킨의 체
    '''use sieve of atkins to find primes <= limit.'''
    primes = [0] * limit
 
    # n = 3x^2 + y^2 section
    xx3 = 3
    for dxx in xrange( 0, 12*int( sqrt( ( limit - 1 ) / 3 ) ), 24 ):
        xx3 += dxx
        y_limit = int(12*sqrt(limit - xx3) - 36)
        n = xx3 + 16
        for dn in xrange( -12, y_limit + 1, 72 ):
            n += dn
            primes[n] = not primes[n]
 
        n = xx3 + 4
        for dn in xrange( 12, y_limit + 1, 72 ):
            n += dn
            primes[n] = not primes[n]
 
    # n = 4x^2 + y^2 section
    xx4 = 0
    for dxx4 in xrange( 4, 8*int(sqrt((limit - 1 )/4)) + 4, 8 ):
        xx4 += dxx4
        n = xx4+1
 
        if xx4%3:
            for dn in xrange( 0, 4*int( sqrt( limit - xx4 ) ) - 3, 8 ):
                n += dn
                primes[n] = not primes[n]
 
        else:
            y_limit = 12 * int( sqrt( limit - xx4 ) ) - 36
 
            n = xx4 + 25
            for dn in xrange( -24, y_limit + 1, 72 ):
                n += dn
                primes[n] = not primes[n]
 
            n = xx4+1
            for dn in xrange( 24, y_limit + 1, 72 ):
                n += dn
                primes[n] = not primes[n]
 
    # n = 3x^2 - y^2 section
    xx = 1
    for x in xrange( 3, int( sqrt( limit / 2 ) ) + 1, 2 ):
        xx += 4*x - 4
        n = 3*xx
 
        if n > limit:
            min_y = ((int(sqrt(n - limit))>>2)<<2)
            yy = min_y*min_y
            n -= yy
            s = 4*min_y + 4
 
        else:
            s = 4
 
        for dn in xrange( s, 4*x, 8 ):
            n -= dn
            if n <= limit and n%12 == 11:
                    primes[n] = not primes[n]
 
    xx = 0
    for x in xrange( 2, int( sqrt( limit / 2 ) ) + 1, 2 ):
        xx += 4*x - 4
        n = 3*xx
 
        if n > limit:
            min_y = ((int(sqrt(n - limit))>>2)<<2)-1
            yy = min_y*min_y
            n -= yy
            s = 4*min_y + 4
 
        else:
            n -= 1
            s = 0
 
        for dn in xrange( s, 4*x, 8 ):
            n -= dn
            if n <= limit and n%12 == 11:
                    primes[n] = not primes[n]
 
    # eliminate squares        
    for n in xrange(5,int(sqrt(limit))+1,2):
        if primes[n]:
            for k in range(n*n, limit, n*n):
                primes[k] = False
 
    return [2,3] + filter(primes.__getitem__, xrange(5,limit,2))


def prime_number(max_n): #에라토스테네스의 체
    p = [0 for i in xrange(max_n+1)]
    for i in xrange(2,max_n+1):
        for j in xrange(i*2,max_n+1,i):
            p[j] = 1
    prime = [i for i in xrange(2,max_n+1) if p[i]==0]
    return prime

s = time.time()
prime_number(2000000)
print '에라토스테네스의 체 :',time.time() - s
s = time.time()
atkins13(2000000)
print '앳킨의 체 :',time.time() - s

---

2016.10.24 추가.

위의 엣킨의 체보다 더 빠른걸 찾았다. 1억까지의 소수를 구할때 위의 엣킨의 체로는 16초가 걸리는데

아래 코드로는 4초가 걸린다. 코드도 더 짧고 좋은듯

# -*- encoding: cp949 -*-
def primes(n):
    correction = (n%6>1)
    n = {0:n,1:n-1,2:n+4,3:n+3,4:n+2,5:n+1}[n%6]
    sieve = [True] * (n/3)
    sieve[0] = False
    for i in xrange(int(n**0.5)/3+1):
      if sieve[i]:
        k=3*i+1|1
        sieve[      ((k*k)/3)      ::2*k]=[False]*((n/6-(k*k)/6-1)/k+1)
        sieve[(k*k+4*k-2*k*(i&1))/3::2*k]=[False]*((n/6-(k*k+4*k-2*k*(i&1))/6-1)/k+1)
    return [2,3] + [3*i+1|1 for i in xrange(1,n/3-correction) if sieve[i]]

또 여기 링크에는 어떤분이 속도비교를 올려놓으셨다.

http://hashcode.co.kr/questions/2353/n-%EC%9D%B4%ED%95%98%EC%9D%98-%EC%86%8C%EC%88%98%EB%A5%BC-%EB%82%98%EC%97%B4%ED%95%98%EB%8A%94-%EA%B0%80%EC%9E%A5-%EB%B9%A0%EB%A5%B8-%EB%B0%A9%EB%B2%95