소수 구하는 방법 (Sieve of Eratosthenes)
그동안 내가 써왔던 소수 구하는 방법을 한번 정리해보았다..
우선 어떤 수 n이 소수인지 아닌지 알아보려면 n보다 작고 1보다 큰 모든 수로 나눠보면 된다..
그러나.. 사실 잘 생각해보면 n보다 작은 소수로만 나누어보면 된다..
그러나.. 더 곰곰히 생각해보면 sqrt(n) 보다 작은 소수로만 나눠보면 된다.. 왜 그럴까.. ㅎㅎ
다음은 50000보다 작은 소수를 모두 구해서 저장하는 함수..
어떤 수가 소수인지 아닌지 알아보기 위해 현재까지 구해진 소수로 다 나누어보고..
한번도 나누어 떨어지지 않으면 소수이므로 소수 리스트에 저장한다..
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의 배수 다 지우고.. 등등..
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에서 배낀거같은데.. 코드의 원리는 모름.. -_-
#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 |