에라토스테네스의 체
일정 범위 내의 소수를 찾는 가장 보편적인 방법은 에라토스테네스의 체이다.
에라토스테네스의 체에 대해 간략히 설명하자면, 기본적으로 2는 소수로 처리해놓고 [2..Range)를 순회하며 소수인 것들의 배수들에 대해 not prime으로 만들며 반복하는 알고리즘이다.
const int RANGE = 1e7;
int np[RANGE]; // np[x] = x가 not prime인지?
void getEratosthenes() {
for(int i = 2; i*i < RANGE; ++i) {
if(!np[i]) {
for(int j = i*i; j < RANGE; j += i) {
np[j] = 1;
}
}
}
}
소수들의 배수들에 대해 not prime으로 처리 할 때 i*i에서부터 시작한다는 최적화를 했지만, 결국 중복되는 수에 대해 not prime으로 처리하는 경우가 생긴다.
그 예로 $30=2\cdot3\cdot5$에서 $i=2$일 때 $2\cdot15$에서 $30$을 not prime 처리하고, $i=3$일 때 $30$을 중복으로 not prime 처리하게 된다.
이로 인해 총 시간복잡도는 $O(n\log{\log{n}})$이 된다. 증명
오일러의 체의 경우는?
에라토스테네스의 체도 충분히 빠르고 간편한데 이보다 더 최적화를 할 수 있을까?
놀랍게도 오일러의 체를 사용하면 $O(n)$만에 소수들을 구할 수 있다.
또한 이 알고리즘은 구현을 위해 각 수의 minimum prime factor(최소 소인수 = spf)를 구하기 때문에 $O(\log{n})$만에 범위 내의 수를 소인수분해 할 수 있다.
어떻게 줄여야 할까?
에라토스테네스의 체에서 중복이 나오는 경우를 잘 살펴보면 소수의 배수들에 대해 not prime 처리를 하는 부분에서 소수 * 합성수형태인 경우 중복된다는 것을 알 수 있다.
이것을 뒤집어서 생각해 “그냥 수 * 소수에 대해 not prime 처리를 하면 안 될까?” 라는 발상을 해볼 수 있다.
그냥 소수를 곱해버리면 소인수가 3개 이상인 경우에는 또 중복이 나올 게 뻔하니까 모든 수는 $x = spf * number$와 같은 구조로 되어 있음을 이용해 수 * spf인 경우에만 not prime을 처리하면 중복을 피할 수 있다는 아이디어가 나온다.
spf는 어떻게 구해야 할까?
$x = number * spf(x)$인 x들에 대해 not prime을 처리하는 것은 에라토스테네스의 체와 비슷한 로직으로 생각하면 편하다.
대략적인 알고리즘은 다음과 같다.
일단 수들을
에라토스테네스의 체와 비슷하게 [2..RANGE)를 순회하며 그들의 spf배들에 대해 not prime을 칠한다.순회 중 현재 not prime이 칠해지지 않은 수들은 소수임이 자명하다.
not prime을 칠할 때는 boolean 배열을 통해 칠하지 말고 해당 수의
spf를 저장하는 식으로 구현해 각 수의spf를 저장해준다.
이 때 문제가 현재 순회중인 수가 $x$라고 했을 때 $y = x*spf(y)$를 not prime으로 칠할 때 $spf(y)$가 어떻게 $y$의 $spf$인지 알고 칠하냐는 것이다.
상식적으로 생각하면 쉬운 문제인데, 현재 $x$에 대해 $y$가 될 수 있는 $spf(y)$를 고르는 조건은 아래 두 가지를 만족시키면 된다.
- 소수여야한다.
- $spf(y) \leq spf(x)$
즉, 현재까지 구해진 소수들을 $x$에 곱해가며 $y$를 칠해가는데, 소수가 $x$보다 커지면 루프를 끊는 식으로 구현하면 된다.
코드
#define FOR(i, f, n) for(int (i) = (f); (i) < (int)(n); ++(i))
const int RANGE = 1e7;
int pn, spf[RANGE], pr[RANGE];
void eulerSieve() {
FOR(x, 2, RANGE) {
if(!spf[x]) spf[x] = pr[pn++] = x; // 현재 x에 대해 spf가 정해지지 않았다는 것은 spf[x] = x라는 의미이고, x가 소수라는 의미이다.
for(int j = 0; x*pr[j] < RANGE; ++j) {
spf[x*pr[j]] = pr[j]; // 수 * spf들에 대해 spf[x*spf] = spf로 칠해준다.
if(x % pr[j] == 0) break; // spf[x] == pr[j]이면 그만 둔다. (이 이상 pr이 늘어나면 x * pr에 대해 더 이상 spf가 pr이 아닌 x의 spf이니까)
}
}
}
시간복잡도 분석
분석이라고 할 것도 없이 어떤 합성수 $y = x * spf(y)$를 not prime이라고 칠할 때 가능한 x는 딱 한 가지이므로 무조건 한 번 밖에 칠해지지 않는다.
실제로 $Range = 10^8$에 대해서 not prime이 칠해지는 경우를 세어보면 다음과 같다.
Euler’s Sieve: 94,238,543회
Eratosthenes’s Sieve: 242,570,202회
spf를 이용한 소인수분해
spf[x]자체가 x의 소인수를 담고 있기 때문에 단순히 구현만 하면 된다.
x의 소인수개수 만큼 순회하기 때문에 시간복잡도는 $O(n)$이다.
const int RANGE = 1e7;
int spf[RANGE];
void f(int x) {
assert(x > 0);
while(x > 1) {
cout << spf[x] << ' ';
x /= spf[x];
}
}
관련 문제
코드포스 Div2-C 문제이다.
이 문제를 통해서 오일러의 체에 대해서 알게 되었는데, 일반적인 에라토스테네스의 체나 소인수 분해를 사용하면 시간초과가 나는 문제이다.
$n$개의 숫자들이 주어지는데, 이 중 몇 개를 지워야 원래 숫자들의 $gcd$보다 남은 숫자들의 $gcd$가 커질까?
$(2 \leq n \leq 3 \cdot 10^5)$, $(1 \leq {a_i} \leq 1.5 \cdot 10^7)$
$gcd$란 모든 숫자들의 공통 인수 중 가장 큰 것이다.
따라서 모든 수들을 $gcd$로 나누고 난 나머지 수들의 공통 소인수는 없을 것이고, 몇 개의 수를 지워 공통 소인수를 만드는 순간 $gcd$가 증가하게 될 것이다.
즉, 모든 수들을 소인수분해해서 cnt[x]배열에 각 소인수가 몇 번 등장했는지를 세고, $n - max(cnt[x])$를 구해주면 된다.
이 때 소인수분해를 일반적인 방법으로 해서는 $n = 3 \cdot 10^5$이고, 수의 범위 ${a_i} = 1.5 \cdot 10^7$이기 때문에 시간 초과가 난다.
실제로 $x$를 소인수분해 할 때 $\sqrt{x}$까지의 소수들로 나눠보는 방법을 사용해봤지만 시간초과가 났다.
하지만 위에서 설명한 오일러의 체를 사용하면 $O(A + n \cdot \log{A})$ 만에 문제를 해결 할 수 있다.
#ifndef ONLINE_JUDGE
#define _CRT_SECURE_NO_WARNINGS
#endif
#include <bits/stdc++.h>
using namespace std;
// ........................macro.......................... //
// [i, n)
#define FOR(i, f, n) for(int (i) = (f); (i) < (int)(n); ++(i))
// [i, n]
#define RFOR(i, f, n) for(int (i) = (f); (i) >= (int)(n); --(i))
#define pb push_back
#define emb emplace_back
#define fi first
#define se second
#define ENDL '\n'
#define sz(A) (int)(A).size()
#define ALL(A) A.begin(), A.end()
#define UNIQUE(c) (c).resize(unique(ALL(c)) - (c).begin())
#define next next9876
#define prev prev1234
typedef vector<int> vi;
typedef vector<vector<int> > vvi;
typedef vector<pair<int, int> > vii;
typedef vector<vector<pair<int, int> > > vvii;
typedef pair<int, int> ii;
typedef pair<int, pair<int, int> > iii;
typedef long long i64;
typedef unsigned long long ui64;
// inline i64 GCD(i64 a, i64 b) { if(b == 0) return a; return GCD(b, a % b); }
inline int getidx(const vi& ar, int x) { return lower_bound(ALL(ar), x) - ar.begin(); } // 좌표 압축에 사용: 정렬된 ar에서 x의 idx를 찾음
inline i64 GCD(i64 a, i64 b){ i64 n; if(a<b) swap(a, b); while(b!=0) { n = a%b; a = b; b = n; } return a; }
inline i64 LCM(i64 a, i64 b) { if (a == 0 || b == 0) return GCD(a, b); return a / GCD(a, b) * b; }
inline i64 CEIL(i64 n, i64 d) { return n / d + (i64)(n % d != 0); } // 음수일 때 이상하게 작동할 수 있음.
inline i64 ROUND(i64 n, i64 d) { return n / d + (i64)((n % d)*2 >= d); }
inline i64 POW(i64 a, i64 n) {
if (n < 0) return 0;
i64 ret = 1;
while (n) { if (n % 2) ret *= a; a = a * a; n /= 2; }
return ret;
}
template <class T>
ostream& operator<<(ostream& os, vector<T> v) {
os << "[";
int cnt = 0;
for (auto vv : v) { os << vv; if(++cnt < v.size()) os << ","; }
return os << "]";
}
template <class T>
ostream& operator<<(ostream& os, set<T> v) {
os << "[";
int cnt = 0;
for (auto vv : v) { os << vv; if(++cnt < v.size()) os << ","; }
return os << "]";
}
template <class L, class R>
ostream& operator<<(ostream& os, pair<L, R> p) {
return os << "(" << p.fi << "," << p.se << ")";
}
void debug_out() { cerr << endl; }
template <typename Head, typename... Tail>
void debug_out(Head H, Tail... T) {
cerr << " " << H;
debug_out(T...);
}
#ifndef ONLINE_JUDGE
#define debug(...) cerr << "[" << #__VA_ARGS__ << "]:", debug_out(__VA_ARGS__)
#else
#define debug(...) 42
#endif
// ....................................................... //
const int MAXN = 3e5, RANGE = 1.5e7+10;
int n, g, pn, ar[MAXN], spf[RANGE], pr[RANGE], cnt[RANGE];
void input() {
cin >> n;
FOR(i, 0, n) cin >> ar[i];
}
int solve() {
g = ar[0];
FOR(i, 1, n) g = __gcd(ar[i], g);
FOR(i, 0, n) ar[i] /= g;
// euler's sieve
FOR(x, 2, RANGE) {
if(!spf[x]) spf[x] = pr[pn++] = x;
for(int j = 0; x*pr[j] < RANGE; ++j) {
spf[x*pr[j]] = pr[j];
if(spf[x] == pr[j]) break;
}
}
// factional decomposition
FOR(i, 0, n) {
while(ar[i] > 1) {
int x = spf[ar[i]];
++cnt[x];
while(ar[i] % x == 0) ar[i] /= x;
}
}
int ans = *max_element(cnt, cnt + RANGE);
if(ans == 0) cout << -1 << ENDL;
else cout << n-ans << ENDL;
return 0;
}
// ................. main .................. //
void execute() {
input(), solve();
}
int main(void) {
#ifndef ONLINE_JUDGE
freopen("in.txt", "r", stdin);
// freopen("out.txt", "w", stdout);
#endif
cin.tie(0); ios_base::sync_with_stdio(false);
execute();
return 0;
}
// ......................................... //
링크
PS
사실 이 알고리즘의 이름이 Euler's Sieve인지는 확실하지 않다.
1978년 David Gries, Jayadev Mirsa의 A Linear Sieve Algorithm for Finding Prime Numbers에서 나온 알고리즘이라는데, 코드포스 튜토리얼에서는 Euler's Sieve라고 하길래 그냥 그렇게 표현했다.