문제 풀다보면 소수를 구해야하는 경우가 많이 생긴다..
그동안 내가 써왔던 소수 구하는 방법을 한번 정리해보았다..

우선 어떤 수 n이 소수인지 아닌지 알아보려면 n보다 작고 1보다 큰 모든 수로 나눠보면 된다..
그러나.. 사실 잘 생각해보면 n보다 작은 소수로만 나누어보면 된다..
그러나.. 더 곰곰히 생각해보면 sqrt(n) 보다 작은 소수로만 나눠보면 된다.. 왜 그럴까.. ㅎㅎ


다음은 50000보다 작은 소수를 모두 구해서 저장하는 함수..
어떤 수가 소수인지 아닌지 알아보기 위해 현재까지 구해진 소수로 다 나누어보고..
한번도 나누어 떨어지지 않으면 소수이므로 소수 리스트에 저장한다..

int prime_cnt;
int prime[50000];
void init_prime()
{
    int i, j, n;
    n = 50000;
    prime[0] = 2;
    prime_cnt = 1;
    for (i = 3; i < n; i += 2) {
        for (j = 1; j < prime_cnt; j++) {
            if (i % prime[j] == 0)
                break;
        }
        if (j == prime_cnt)
            prime[prime_cnt++] = i;
    }
}

예전에는 위의 코드를 자주 사용했는데.. 이제는 비효율적이어서 잘 안쓰고.. 에라토스테네스의 체를 사용한다..
2의 배수 다 지우고.. 3의 배수 다 지우고.. 5의 배수 다 지우고.. 등등..

/* Sieve of Eratosthenes */
char is_prime[1000001];
int prime[100000];
int prime_cnt;
void sieve()
{
    int i, j;
    memset(is_prime, -1, sizeof(is_prime));
    prime_cnt = 0;
    is_prime[0] = is_prime[1] = 0;
    for (i = 2; i <= 1000000; i++) {
        if (is_prime[i]) {
            for (j = 2; i * j <= 1000000; j++) {
                is_prime[i*j] = 0;
            }
            prime[prime_cnt++] = i;
        }
    }
}

사실 is_prime[] 배열만 구하려면 i 루프는 sqrt(1000000) 까지만 돌면 된다..
굳이 100만까지 도는 이유는 100만 보다 작은 모든 소수의 리스트를 구하기 위해서이다..
그리고 memset에서 -1로 초기화하는 이유는 나의 단순한 습관때문이다..
배열이 int형이거나 long long 형이라도 항상 원하는 값으로 초기화하려고 -1로 초기화 한다..
속도를 조금이나마 더 빠르게하려고.. 2의 배수부터 다 지우고.. i를 2씩 증가시키는 사람도 있다..
나는 그렇게 안한다.. 귀찮아서.. -_-

어떤 수가 소수인지 아닌지를 저장하기 위해 한바이트도 아까워서..
한바이트에 여러 정보를 넣는 bitwise sieve도 있다..
나는 거의 쓰지는 않지만.. 속도도 빠르고 메모리도 적게쓰는 방법..
예전에 어디 UVa Board에서 배낀거같은데.. 코드의 원리는 모름.. -_-

// Yarin's sieve
#define MAXSIEVE 1000000 // All prime numbers up to this
#define MAXSIEVEHALF (MAXSIEVE/2)
#define MAXSQRT sqrt((double)MAXSIEVE)/2
char a[MAXSIEVE/16+2];
#define isprime(n) (a[(n)>>4]&(1<<(((n)>>1)&7))) // Works when n is odd
void initprime(){
    int i,j;
    memset(a,255,sizeof(a));
    a[0]=0xFE;
    for(i=1;i<MAXSQRT;i++)
        if (a[i>>3]&(1<<(i&7)))
            for(j=2*i*(i+1);j<MAXSIEVEHALF;j+=2*i+1)
                a[j>>3]&=~(1<<(j&7));
}


Sieve 중에는 특정 구간만 체로 치는 방법도 있다..
예를들어 십억부터 십억백만까지의 소수를 다 구하고자 할때.. 배열을 십억개를 잡을 수 없으므로 난감하다..
이때는 segmented sieve를 사용한다..

/* Segmented Sieve */
/* if 'i' is prime, is_prime[i-L] = 1 */
int is_prime[1000001];
void segmented_sieve(int l, int u)
{
    int i, j, d, sq;
    d = u - l + 1;
    for (i = 0; i < d; i++)
        is_prime[i] = 1;
    for (i = (l%2 != 0); i < d; i+= 2)
        is_prime[i] = 0;

    /* sieve by prime factors starting from 3 till sqrt(u) */
    sq = (int)sqrt(u);
    for (i = 3; i <= sq; i+= 2) {
        if (i > l && !is_prime[i-l])
            continue;
        /* choose the first number to be sieved..
           >= l && divisible by i, and not i itself */
        j = (l / i) * i;
        if (j < l)
            j += i;
        if (j == i)
            j += i;
        j -= l;
        for (; j < d; j+=i)
            is_prime[j] = 0;
    }
    /* mark 1 as false, 2 as true */
    if (l <= 1)
        is_prime[1-l] = 0;
    if (l <= 2)
        is_prime[2-l] = 1;
}

뭐.. 내가 짠거는 아니고.. 다 어디서 배껴온 코드들.. -_-;;

사실 prime number에 관한 이슈로는 primality testing 같은게 더 있는데.. 아직 거기까지 공부를 안해서..;;
나중에 또 포스팅할 기회가 있기를..

'Problem Solving > Algorithm notes' 카테고리의 다른 글

Finding Minimum Path Cover in DAG  (0) 2009.06.15
Plane Equation (평면의 방정식)  (0) 2009.04.16
Ellipse (타원)  (0) 2009.04.15
Combination 개수 구하기 (Pascal's Triangle)  (2) 2009.02.07
Number of Swap Operations  (0) 2008.07.24
Horner's Rule  (0) 2008.05.04
GCD SUM  (0) 2008.03.18
Erdos & Gallai  (0) 2008.03.04
Misère Nim  (2) 2007.12.16
BSP Tree  (0) 2007.08.28

to Top