미러에서 풀었어야 했는데 괜히 능지 더 쓰려다가 못 풀었다. 3시간 째 실제 대회에서 푼 사람이 안 나오길래 엄청 어려운 줄... 

그래도 문자열 문제 어떻게 푸는지 기억이 아예 사라진 것은 아니어서 조금 신기했다. 가뭄에 콩 나듯 하는 PS....

 

Suffix Array + LCP를 만들자. 여기서 만약에 $$\max(lcp_{i-1}, lcp_{j+1}) < \min(lcp_i, lcp_{i+1}, \cdots, lcp_j)$$가 성립한다면, $i-1$번째 substring 부터 $j$번째 substring까지 총 $j-i+2$개의 substring이 $$\min(lcp_i, lcp_{i+1}, \cdots, lcp_j)$$ 길이의 common prefix를 가지며, 특히 $$[\max(lcp_{i-1}, lcp_{j+1}) + 1, \min(lcp_i, lcp_{i+1}, \cdots, lcp_j)]$$ 구간에 속하는 길이의 prefix 들은 정확히 $j-i+2$번 등장하게 된다. 즉, 이들은 $f(j-i+2)$에 들어가는 후보군이라고 볼 수 있다. 

 

그래서 일단 저런 경우들을 모두 관리해야 하는데, 이건 stack을 써서 할 수 있다. 대충 $st$번째 문자열부터 $en$번째 문자열까지 common prefix가 $l$이라는 것을 하나의 struct로, 즉 $(st, en, l)$ 형태로 관리한다고 치자. 그러면 $lcp_i$를 보면서 얻는 것은 $(i-1, i, lcp_i)$라는 struct다.

 

이제 stack을 머리속으로 그리면서 생각을 해보자. 스택의 top을 $S$라고 생각을 하고, 지금 $(i-1, i, lcp_i)$를 본다고 하자.

  • $S_l > lcp_i$인 경우를 생각해보면, 현재 [$S$ 이전의 원소인 $T$], [$S$], [$(i-1, i, lcp_i)$]가 순서대로 stack에 대기중이다.
    • stack 상 invariant를 하나 유지하자. 이제부터 stack에 쌓인 순서대로 $l$ 값이 증가하도록 한다. 자연스러운 생각이다.
    • 그러면 $S_l > \max(T_l, lcp_i)$이므로, 여기서 $f$ 값을 update 할 케이스가 발생한다. 
    • 즉, $[\max(T_l, lcp_i) + 1, S_l]$ 구간에 속하는 prefix들은 정확히 $S_{en} - S_{st} + 1$번 등장하게 된다. 
    • 이제 $f$를 계산하려면 disjoint occurrence 횟수를 계산해야 한다. 
      • 이는 $S$에 대응되는 suffix array 값들을 set으로 관리하고, naive 하게 lower bound를 사용하면서 계산하면 된다.
      • 결국 "최대 occurrence를 유지하는 최대 길이"를 계산해야 하는데, 이는 $[\max(T_l, lcp_i) + 1, S_l]$에서 이분탐색을 하면 된다.
    • 이제 다시 $T_l \ge lcp_i$인 경우와 $T_l < lcp_i$인 경우로 나뉘게 된다. 
    • $T_l \ge lcp_i$인 경우를 생각해보면, 이제 $S_l = T_l$인 것으로 생각해도 무방하다는 것을 알 수 있다.
      • 그러므로, $S_l \leftarrow T_l$로 $S$를 바꾸고, $S$와 $T$를 하나로 합친다.
      • 여기서 합친다는 것은, $(T_{st}, S_{en}, S_l = T_l)$을 만든다는 것이다. 
      • 만약 여전히 $T_l > lcp_i$이면 다시 맨 처음부터 반복하면 된다.
      • 앞서 확인했듯이 $(st, en, l)$에 대응되는 suffix array 값들의 집합을 set으로 관리해야 한다. 
        • 이는 $S$에 대응되는 집합과 $T$에 대응되는 집합을 합치는 것과 같다. 그러니 small-to-large 하면 ok.
    • $T_l < lcp_i$인 경우, 이제 $S_l < lcp_i$인 경우와 같은 방식으로 대응하면 된다. 
  • $S_l = lcp_i$인 경우, 단순히 $S \leftarrow (S_{st}, i = S_{en} + 1, S_l)$으로 바꾸고 set을 관리하면 된다.
  • $S_l < lcp_i$인 경우, stack에 $(i-1, i, lcp_i)$를 push 하면 된다. 

모든 과정이 끝나면, stack에 원소들이 남을 것이며 $l$이 증가하는 순서일 것이다.

위 과정과 마찬가지로 pop하고 prefix들 handle 하고 스택의 다음 원소와 합치는 과정을 반복하면 답을 모두 구할 수 있다.

 

disjoint occurrence 계산을 naive 하게 할 수 있을 거라고 koosaga가 말해줬는데 안 믿은 게 문제다 ㅠㅠ 

엄청나게 어려운 자료구조를 관리해야 할 거라고 생각해서 망했다. 시간복잡도 증명이 매우 궁금한 문제.

 

꾸베릴과 알대인

 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
#include <stdio.h>
#include <cstring>
#include <iostream>
#include <vector>
#include <algorithm>
#include <string>
#include <map>
#include <stack>
#include <set>
using namespace std;
typedef long long int ll;
 
bool debug = false;
int ord[51111], nord[51111], cnt[51111], aux[51111];
int sfx[51111], rev[51111], lcp[51111];
string s;
int ans[51111]; 
int app[51111];
 
struct L {
    int len_com;
    int st; int en; int idx;
    L(int len_com, int st, int en, int idx): len_com(len_com), st(st), en(en), idx(idx) {}
};
 
int iidx = 0;
set<ll> VALUES[51111];
stack<L> stk;
 
void SFSA(string str) {
    int n = str.length();
    int p = 1;
    memset(ord, 0sizeof(ord));
    for(int i=0; i<n; i++){
        sfx[i] = i;
        ord[i] = str[i];
    }
    int pnt = 1;
    while(1) {
        memset(cnt, 0sizeof(cnt));
        for(int i=0; i<n; i++) cnt[ord[min(i+p, n)]]++;
        for(int i=1; i<=|| i<=255; i++) cnt[i] += cnt[i-1];
        for(int i=n-1; i>=0; i--)
        aux[--cnt[ord[min(i+p, n)]]] = i;
        memset(cnt, 0sizeof(cnt));
        for(int i=0; i<n; i++) cnt[ord[i]]++;
        for(int i=1; i<=|| i<=255; i++) cnt[i] += cnt[i-1];
        for(int i=n-1; i>=0; i--)
        sfx[--cnt[ord[aux[i]]]] = aux[i];
        if(pnt == n) break;
        pnt = 1;
        nord[sfx[0]] = 1;
        for(int i=1; i<n; i++){
            if(ord[sfx[i-1]] != ord[sfx[i]] || ord[sfx[i-1+ p] != ord[sfx[i] + p]){
            pnt++;
        }
        nord[sfx[i]] = pnt;
        }
        memcpy(ord, nord, sizeof(int* n);
        p *= 2;
    }
    for(int i=0; i<n; i++) rev[sfx[i]] = i;
    int h = 0;
    for(int i=0; i<n; i++){
        if(rev[i]){
            int prv = sfx[rev[i] - 1];
            while(str[prv + h] == str[i + h]) h++;
            lcp[rev[i]] = h;
        }
        h = max(h-10);
    }
}
 
int get_count(int idx, int jump) {
    int ret = 0;
    auto it = VALUES[idx].begin();
    while(it != VALUES[idx].end()) {
        ret += 1;
        it = VALUES[idx].lower_bound((*it) + jump);
    }
    return ret;
}
 
void printL(L S) {
    cout << "len_com: " << S.len_com << endl;
    cout << "st: " << S.st << endl;
    cout << "en: " << S.en << endl;
    cout << "idx: " << S.idx << endl;
}
 
void finish_handle(int upper_cut, L S) {
    if(S.len_com == 0return;
    if(debug) {
        cout << "finish_handle " << upper_cut << endl;
        printL(S);
    }
    int num_times = S.en - S.st + 1;
    int num_app = get_count(S.idx, upper_cut);
    if(num_app < app[num_times]) return;
    int l = upper_cut;
    int r = S.len_com;
    int best = upper_cut, mid;
    while(l <= r) {
        mid = (l + r) / 2;
        if(get_count(S.idx, mid) == num_app) {
            best = mid;
            l = mid + 1;
        }
        else r = mid - 1;
    }
    if(debug) {
        cout << "num_times: " << num_times << endl;
        cout << "num_app: " << num_app << endl;
        cout << "best: " << best << endl;
    }
    if(num_app > app[num_times]) {
        app[num_times] = num_app;
        ans[num_times] = best;
    }
    else if(num_app == app[num_times]) {
        ans[num_times] = max(ans[num_times], best);
    }
}
 
L merge(L S1, L S2) {
    int from, to;
    if(VALUES[S1.idx].size() < VALUES[S2.idx].size()) {
        from = S1.idx; to = S2.idx;
    }
    else {
        from = S2.idx; to = S1.idx;
    }
    int new_st = min(S1.st, S2.st);
    int new_en = max(S1.en, S2.en);
    int new_len = min(S1.len_com, S2.len_com);
    int new_idx = to;
    for(auto it = VALUES[from].begin() ; it != VALUES[from].end() ; it++) {
        VALUES[to].insert((*it));
    }
    return L(new_len, new_st, new_en, new_idx);
}
 
int main(void) {
    cin >> s; SFSA(s);
    ans[1= s.length();
    for(int i = 1 ; i < s.length() ; i++) {
        if(debug) cout << "on " << i << endl;
        while(!stk.empty() && lcp[i] < stk.top().len_com) {
            L S = stk.top(); stk.pop();
            int upper_cut = 0;
            if(stk.empty() || stk.top().len_com < lcp[i]) {
                upper_cut = lcp[i] + 1;
                finish_handle(upper_cut, S);
                stk.push(L(lcp[i], S.st, S.en, S.idx));
            }
            else {
                upper_cut = stk.top().len_com + 1;
                finish_handle(upper_cut, S);
                L S2 = stk.top(); stk.pop();
                L TT = merge(S2, S);
                stk.push(TT);
 
            }
        }
            
        if(stk.empty() || lcp[i] > stk.top().len_com) {
            iidx++;
            VALUES[iidx].insert(sfx[i-1]); VALUES[iidx].insert(sfx[i]);
            L Lv = L(lcp[i], i - 1, i, iidx);
            stk.push(Lv); continue;
        }
 
        if(!stk.empty() && stk.top().len_com == lcp[i]) {
            VALUES[stk.top().idx].insert(sfx[i]);
            L stk_top = stk.top(); stk.pop();
            stk.push(L(stk_top.len_com, stk_top.st, i, stk_top.idx));
        } 
    }
 
    while(!stk.empty()) {
        L S = stk.top(); stk.pop();
        int upper_cut = 1;
        if(!stk.empty()) { upper_cut = max(upper_cut, stk.top().len_com + 1); }
        finish_handle(upper_cut, S);
        if(!stk.empty()) {
            L S2 = stk.top(); stk.pop();
            L TT = merge(S2, S);
            stk.push(TT);
        }
    }
 
    for(int i = 1 ; i <= s.length() ; i++cout << ans[i] << " ";
    cout << endlreturn 0;    
}
 
cs

 

 

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

10, 11월의 PS 일지 - Part 1  (0) 2020.11.08
8월의 PS 일지 - Part 6  (0) 2020.08.31
8월의 PS 일지 - Part 4  (0) 2020.08.12
8월의 PS 일지 - Part 2  (0) 2020.08.09
7월의 PS 일지 - Part 2  (0) 2020.08.04

main.pdf
3.41MB

쓰기 시작한지는 꽤 되었는데

  • 회사 + 연구 + CTF하면서 시간이 잘 안나서 + 게을러져서 (...) 진도가 안빠지고
  • 이대로 두었다가는 지금까지 쓴 부분도 공개가 너무 늦어질 것 같고
  • 사실 지금까지 쓴 부분만 잘 이해해도 문제 해결에 지장이 없어서

미리 올립니다. 솔직히 머리 속에는 어떻게 글을 써야하는지 아는데 실제로 글을 쓰려면 막막해서 슬프네요 :(

 

그래도 이거 올리고 사람들이 좀 읽어주면 이 글을 다시 쓸 motivation을 얻을 것 같기도 하고 잘 모르겠습니다

 

LaTeX 형식은 Evan Chen의 Napkin에서 가져왔습니다. 참고하세요.

4번 문제를 5시가 되기 직전에 해결했다. 사진은 7시 쯤 찍었다.

5번은 긁기는 쉬운데 만점을 받는 사람이 없길래 그냥 건드리지 않았다. 4번도 풀었고 ^__^

 

1번 : 원 안의 점 

naive 하게 $-R \le x, y \le R$인 점을 모두 시도하면 $\mathcal{O}(R^2)$ 풀이가 나온다. 

답을 $\mathcal{O}(R)$에 구하기 위해서는, $x$ 하나를 고정하고 가능한 $y$의 개수를 $\mathcal{O}(1)$에 구하면 된다.

$x^2+y^2 \le R^2-1 \iff -\sqrt{R^2-1-x^2} \le y \le \sqrt{R^2 -1-x^2}$임을 계산하면, 가능한 $y$의 개수가 $$2 \lfloor \sqrt{R^2 - 1 - x^2} \rfloor + 1$$임을 알 수 있다. 이는 sqrt 함수를 사용하거나 이분탐색을 해서 빠르게 계산할 수 있다. 

 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
#include <bits/stdc++.h>
#define fio ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
using namespace std;
typedef long long int ll;
typedef unsigned long long int ull;
typedef long double ldb;
const ll mod = 1e9 + 7;
 
ll calc(ll x)
{
    ll lef=0, rig=1e9, best=0, mid;
    while(lef<=rig)
    {
        mid = (lef + rig) / 2;
        if(mid * mid <= x) 
        {
            best = mid;
            lef = mid + 1;
        }
        else rig = mid - 1;
    }
    return best;
}
 
void solve(void)
{
    ll R, ans=0cin >> R;
    for(ll i=-R+1 ; i<=R-1 ; i++)
        ans += 2 * calc(R * R - 1 - i * i) + 1;
    cout << ans << endl;
}
 
 
int main(void)
{
    fio; ll i, tc; cin >> tc;
    for(i=1 ; i<=tc ; i++
    {
        cout << "Case #" << i << endl;
        solve(); // endl in solve
    }
    return 0;
}
cs

 

2번 : 직8각형

점 $(x_1, y_1), \cdots (x_8, y_8)$이 있을 때, 적당한 $X, Y$를 골라서 이 점들을 $$(X, Y), (X, Y+K), (X+K, Y+2K), (X+2K, Y+2K)$$ $$(X+3K, Y+K), (X+3K, Y), (X+2K, Y-K), (X+K, Y-K)$$와 같도록 해야 한다. 각 8개의 점을 직8각형의 8개의 점에 대응시키는 방법에는 총 $8!$가지가 있다. 이러한 방법을 하나 고정시키고 생각하자.

 

예를 들어, $(x_1, y_1), \cdots , (x_8, y_8)$이 위 8개의 점과 순서대로 대응된다고 가정하자. 이 경우, 움직여야 하는 총 거리는 $$|x_1 - X| + |y_1 - Y| + |x_2 - X| + |y_2 - (Y+K)| + |x_3 - (X+K)| + |y_3 - (Y+2K)|$$ $$+|x_4 - (X+2K)| + |y_4 - (Y+2K)| + |x_5 - (X+3K)| + |y_5 - (Y+K)| + |x_6 - (X+3K)| + |y_6 - Y|$$ $$+|x_7 - (X+2K)| + |y_7 - (Y-K)| + |x_8 - (X+K)| + |y_8 - (Y-K)|$$가 되며, 우리의 목표는 이 식의 최솟값을 구하는 것이다.

첫 번째 관찰은 위 식에서 $X, Y$가 각각 독립적으로 나온다는 것이다. 위 식에서 $X$가 등장하는 부분만 살펴보면, $$|X-x_1| + |X-x_2| + |X - (x_3-K)| + |X - (x_4-2K)|$$ $$+ |X - (x_5-3K)| + |X - (x_6-3K)| + |X - (x_7 - 2K)| + |X - (x_8 -K)|$$ 두 번째 관찰은, $x_i$들과 $K$가 이미 알고 있는 값이므로, 이 식은 이미 알고 있는 값 $c_1, c_2, \cdots, c_8$에 대해 $$\sum_{i=1}^8 |X - c_i|$$ 형태로 쓸 수 있다는 것이다. 이 식은 $X$가 $c_i$들의 "중간"에 있을 때 최솟값을 가진다.

확인하고 싶다면, 위 식을 그래프로 그려보자. 또한, 이 때 최솟값은 $c_i$들 중 최대 4개에서 최소 4개를 뺀 값이 됨을 확인할 수 있다.

$Y$에 대한 부분도 이렇게 최솟값을 구할 수 있으며, 두 결과를 합치면 식의 최솟값을 구할 수 있다. 이를 $8!$번 반복하면 해결.

 

시간복잡도는 $n$이 점의 개수라고 하면, $\mathcal{O}(n! \cdot n \log n)$ 정도로 볼 수 있을 것이다.

 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
#include <bits/stdc++.h>
#define fio ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
using namespace std;
typedef long long int ll;
typedef unsigned long long int ull;
typedef long double ldb;
const ll mod = 1e9 + 7;
 
pair<ll, ll> pts[8]; ll K;
ll x[8], y[8];
ll difx[8= {00123321};
ll dify[8= {012210-1-1};
 
void solve(void)
{
    ll i, ans=1e18cin >> K; 
    for(i=0 ; i<8 ; i++cin >> pts[i].first >> pts[i].second;
    sort(pts, pts+8);
    do
    {   
        for(i=0 ; i<8 ; i++) x[i] = pts[i].first - difx[i] * K;
        for(i=0 ; i<8 ; i++) y[i] = pts[i].second - dify[i] * K;
        sort(x, x+8); sort(y, y+8);
        ll cur = 0;
        for(i=4 ; i<8 ; i++) cur += (x[i] + y[i]);
        for(i=0 ; i<4 ; i++) cur -= (x[i] + y[i]);
        ans = min(ans, cur);
    } while(next_permutation(pts, pts+8));
    cout << ans << endl;
}
 
int main(void)
{
    fio; ll i, tc; cin >> tc;
    for(i=1 ; i<=tc ; i++
    {
        cout << "Case #" << i << endl;
        solve(); // endl in solve
    }
    return 0;
}
cs

 

3번 : 산탄총

정말 피곤한 문제였다. 진짜 부분합만 알면 풀 수 있는데, 깔끔하게 안 풀어서 그런지 계산이 너무 귀찮았다.

이제부터 모든 것을 수식으로 설명한다. 하지만 실제로 문제를 풀거나 이 풀이를 읽을 때는 그림을 그려가면서 생각하는 것이 좋을 것이다.

 

배열은 1-index를 사용하여, $1 \le i, j \le N$에 배열 $a$가 채워졌다고 생각하도록 하겠다.

 

우선, 목표로 하는 함수인 $$ \text{score}(Y, X) = \sum_{|y-Y| + |x-X| \le K-1} (K-|y-Y|-|x-X|) \cdot a[y][x]$$를 정의하면, 우리의 목표는 $\text{score}$의 최댓값을 찾는 것이다. 우선 $Y$ 또는 $X$가 $-K$ 미만이거나 $N+K$ 초과면, $\text{score}$ 함수가 $0$이 됨을 확인할 수 있다.

그러니 $-K \le X, Y \le N+K$인 경우만 계산하면 된다. 우리의 궁극적인 목표는 $\text{score}$ 함수를 이 범위 전체에 대해서 $\mathcal{O}((N+K)^2) = \mathcal{O}(N^2)$ 시간에 계산하는 것이다. 이를 위해서는 $\text{score}$ 함수가 $(Y, X)$가 한 칸 움직였을 때 어떻게 변화하는지 알아볼 필요가 있다. 

$$\text{score}(Y, X) - \text{score}(Y, X-1)$$ $$= \sum_{|y-Y| + |x-X| \le K-1} (K-|y-Y|-|x-X|) \cdot a[y][x] - \sum_{|y-Y| + |x-X+1| \le K-1} (K - |y-Y| - |x-X+1|) \cdot a[y][x]$$

이 값을 분석하는 가장 좋은 방법은 각 $(y, x)$에 대해 $a[y][x]$의 계수를 찾는 것이다. 그림을 그려 풀 때도 비슷하다.

 

우선 생각해보면 $|y-Y|+|x-X|$와 $|y-Y|+|x-(X-1)|$이 모두 $K$ 이하인 경우, $a[y][x]$의 계수는 $$|x-X+1| - |x-X|$$가 된다. 이 값은 $x \ge X$에서 $1$이고 $x \le X-1$에서 $-1$임을 알 수 있다. 

$|y-Y|+|x-X|$나 $|y-Y|+|x-(X-1)|$ 중 하나라도 $K+1$ 이상인 경우, 둘 다 $K$ 이상이 되어 계수는 $0$이다. 즉,

$$\text{score}(Y, X) - \text{score}(Y, X-1) = \sum_{|y-Y| + |x-X| \le K-1, x \ge X} a[y][x] - \sum_{|y-Y| + |x-X+1| \le K-1, x \le X-1} a[y][x] $$가 되고, 비슷한 원리로 계산하면 $$\text{score}(Y, X) - \text{score}(Y-1, X) = \sum_{|y-Y|+|x-X| \le K-1, y \ge Y} a[y][x] - \sum_{|y-Y+1|+|x-X| \le K-1,  y \le Y-1} a[y][x]$$를 얻는다. 그러니까 이제 필요한 것은 네 종류의 "삼각형 부분합"이며, 이들 역시 같은 방법으로 계산이 가능하다.

 

예를 들어, "1사분면에 대한 삼각형 부분합"을 $$Q1(Y, X) = \sum_{|y-Y| + |x-X| \le K-1, y \le Y, x \ge X} a[y][x]$$라 하자. 이 값을 계산하기 위하여, 다시 인접한 $Q1$ 값의 차이를 계산해보면, $$Q1(Y, X) - Q1(Y, X-1) = \sum_{|y-Y| + |x-X| \le K-1, y \le Y, x \ge X} a[y][x] - \sum_{|y-Y| + |x-X+1| \le K-1, y \le Y, x \ge X-1} a[y][x]$$인데, $a[y][x]$의 계수를 생각해보면 다음과 같은 결과를 얻을 수 있다. 그림을 그려서 생각하는 게 편하다.

  • $x \ge X, y \le Y$이고 $(Y-y) + (x-X) = |y-Y|+|x-X| = K-1$이면 계수가 $+1$
  • $x = X-1$이고 $Y-K+1 \le y \le Y$면 계수가 $-1$
  • 나머지 경우에 대해서는 전부 계수가 $0$

첫 번째 경우에 대한 합은 "대각선 부분합"이고, 두 번째 경우에 대한 합은 "일직선 부분합"이니, 전부 부분합 테크닉으로 빠르게 구할 수 있다.

그러니 $Q1(Y, X-1)$이 있으면 $Q1(Y, X)$를 구할 수 있고, 비슷한 계산으로 $Q1(Y-1, X)$이 있으면 $Q1(Y, X)$를 빠르게 구하는 식을 얻을 수 있다.

 

$Q1(-K, -K)$를 naive 하게 직접 구한 후 위 방법을 적용하면 모든 $Y, X$에 대해 $Q1(Y, X)$를 $\mathcal{O}(N^2)$에 구할 수 있다.

이를 각 4개의 사분면에 대해서 적용할 수 있고, 이제 $\text{score}$ 역시 전부 $\mathcal{O}(N^2)$에 구할 수 있다. 

 

당연하지만 이 문제에서는 index가 범위를 벗어나는 것에 대한 처리가 매우 귀찮고 중요하다.

 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
#include <bits/stdc++.h>
#define fio ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
using namespace std;
typedef long long int ll;
typedef unsigned long long int ull;
typedef long double ldb;
const ll mod = 1e9 + 7;
 
ll K, N;
ll val[1811][1811];
ll DIAG1[1811][1811]; // NW SE
ll DIAG2[1811][1811]; // NE SW
ll WE[1811][1811];
ll NS[1811][1811];
ll CALC[1811][1811];
ll QUAD[4][1811][1811];
 
void rset(void)
{
    memset(val, 0sizeof(val));
    memset(DIAG1, 0sizeof(DIAG1));
    memset(DIAG2, 0sizeof(DIAG2));
    memset(WE, 0sizeof(WE));
    memset(NS, 0sizeof(NS));
    memset(CALC, 0sizeof(CALC));
    memset(QUAD, 0sizeof(QUAD));
}
 
ll NS_p(ll Y1, ll Y2, ll X)
{
    Y1 = min(Y1, N + 2 * K);
    Y2 = max(Y2, 0LL);
    return NS[Y1][X] - NS[Y2][X];
}
 
ll WE_p(ll Y, ll X1, ll X2)
{
    X1 = min(X1, N + 2 * K);
    X2 = max(X2, 0LL);
    return WE[Y][X1] - WE[Y][X2];
}
 
ll getdiag1(ll Y1, ll X1)
{
    if(Y1 < 0 || X1 < 0return 0;
    else if(Y1 <= N+2*&& X1 <= N+2*K) return DIAG1[Y1][X1];
    else
    {
        ll offset = max(Y1 - (N + 2 * K), X1 - (N + 2 * K));
        Y1 -= offset; X1 -= offset;
        if(Y1 < 0 || X1 < 0return 0;
        else return DIAG1[Y1][X1];
    }
}
 
ll DIAG1_p(ll Y1, ll X1, ll Y2, ll X2)
{
    return getdiag1(Y1, X1) - getdiag1(Y2, X2);
}
 
ll getdiag2(ll Y1, ll X1)
{
    if(Y1 < 0 || X1 > N + 2 * K) return 0;
    else if(Y1 <= N+2*&& X1 >= 0return DIAG2[Y1][X1];
    else 
    {
        ll offset = max(Y1 - (N + 2 * K), -X1);
        Y1 -= offset; X1 += offset;
        if(Y1 < 0 || X1 < 0 || X1 > N + 2 * K) return 0;
        else return DIAG2[Y1][X1];
    }
}
 
ll DIAG2_p(ll Y1, ll X1, ll Y2, ll X2)
{
    return getdiag2(Y1, X1) - getdiag2(Y2, X2);
}
 
ll getval(ll Y, ll X)
{
    if(X<=0 || X>N+2*|| Y<=0 || Y>N+2*K) return 0;
    return val[Y][X];
}
 
void finish_QUAD0(void)
{
    ll i, j;
    for(i=0 ; i<=K-1 && i<=1 ; i++)
        for(j=0 ; i+j<=K-1 ; j++)
            QUAD[0][1][1+= getval(1-i, 1+j);
    for(i=1 ; i<=N+2*K ; i++)
    {
        for(j=1 ; j<=N+2*K ; j++)
        {
            if(i == 1 && j == 1continue;
            if(j == 1// get from top
            {
                QUAD[0][i][j] = QUAD[0][i-1][j] 
                               + WE_p(i, j+K-1, j-1)
                               - DIAG1_p(i-1, j+K-1, i-K-1, j-1);
 
            }
            else // get from left
            {
                QUAD[0][i][j] = QUAD[0][i][j-1]
                               + DIAG1_p(i, j+K-1, i-K, j-1
                               - NS_p(i, i-K, j-1);
            }
        }
    }
}
 
void finish_QUAD1(void)
{
    ll i, j;
    for(i=0 ; i<=K-1 && i<=1 ; i++)
        for(j=0 ; i+j<=K-1 && j<=1; j++)
            QUAD[1][1][1+= getval(1-i, 1-j);
    for(i=1 ; i<=N+2*K ; i++)
    {
        for(j=1 ; j<=N+2*K ; j++)
        {
            if(i == 1 && j == 1continue;
            if(j == 1// get from top
            {
                QUAD[1][i][j] = QUAD[1][i-1][j]
                              + WE_p(i, j, j-K)
                              - DIAG2_p(i-1, j-K+1, i-K-1, j+1);
            }
            else // get from left
            {
                QUAD[1][i][j] = QUAD[1][i][j-1]
                              + NS_p(i, i-K, j)
                              - DIAG2_p(i, j-K, i-K, j);
            }
        }
    }
}
 
void finish_QUAD2(void)
{
    ll i, j;
    for(i=0 ; i<=K-1 ; i++)
        for(j=0 ; i+j<=K-1 && j<=1 ; j++)
            QUAD[2][1][1+= getval(1+i, 1-j);
    for(i=1 ; i<=N+2*K ; i++)
    {
        for(j=1 ; j<=N+2*K ; j++)
        {
            if(i == 1 && j == 1continue;
            if(j == 1// get from top
            {
                QUAD[2][i][j] = QUAD[2][i-1][j]
                              + DIAG1_p(i+K-1, j, i-1, j-K)
                              - WE_p(i-1, j, j-K);
            }
            else // get from left
            {
                QUAD[2][i][j] = QUAD[2][i][j-1]
                              + NS_p(i+K-1, i-1, j) 
                              - DIAG1_p(i+K-1, j-1, i-1, j-K-1);
            }
        }
    }
}
 
void finish_QUAD3(void)
{
    ll i, j;
    for(i=0 ; i<=K-1 ; i++)
        for(j=0 ; i+j<=K-1 ; j++)
            QUAD[3][1][1+= val[1+i][1+j];
    for(i=1 ; i<=N+2*K ; i++)
    {
        for(j=1 ; j<=N+2*K ; j++)
        {
            if(i == 1 && j == 1continue;
            if(j == 1// get from top
            {
                QUAD[3][i][j] = QUAD[3][i-1][j]
                              + DIAG2_p(i+K-1, j, i-1, j+K)
                              - WE_p(i-1, j+K-1, j-1);
            }
            else // get from left
            {
                QUAD[3][i][j] = QUAD[3][i][j-1]
                              + DIAG2_p(i+K-1, j, i-1, j+K)
                              - NS_p(i+K-1, i-1, j-1);
            }
        }
    }
}
 
void solve(void)
{
    ll i, j; rset(); cin >> N >> K;
    for(i=1 ; i<=N ; i++)
        for(j=1 ; j<=N ; j++)
             cin >> val[K+i][K+j];
    for(i=0 ; i<=N+2*K ; i++)
    {
        for(j=0 ; j<=N+2*K ; j++)
        {
            if(i!=0 && j!=0) DIAG1[i][j] = DIAG1[i-1][j-1+ val[i][j];
            if(i!=0) DIAG2[i][j] = DIAG2[i-1][j+1+ val[i][j];
            if(j!=0) WE[i][j] = WE[i][j-1+ val[i][j];
            if(i!=0) NS[i][j] = NS[i-1][j] + val[i][j];
        }
    }
    finish_QUAD0(); finish_QUAD1(); finish_QUAD2(); finish_QUAD3();
    for(i=1 ; i<=K ; i++)
    {
        for(j=1 ; j<=K ; j++)
        {
            ll dist = abs(i-1+ abs(j-1);
            if(dist <= K-1) CALC[1][1+= (K - dist) * val[i][j];
        }
    }
    for(i=1 ; i<=N+2*K ; i++)
    {
        for(j=1 ; j<=N+2*K ; j++)
        {
            if(i == 1 && j == 1continue;
            if(j == 1// get from top
            {
                CALC[i][j] = CALC[i-1][j];
                CALC[i][j] += (QUAD[2][i][j] + QUAD[3][i][j] - NS_p(i+K-1, i-1, j));
                CALC[i][j] -= (QUAD[0][i-1][j] + QUAD[1][i-1][j] - NS_p(i-1, i-K-1, j));
            }
            else // get from left 
            {
                CALC[i][j] = CALC[i][j-1];
                CALC[i][j] += (QUAD[0][i][j] + QUAD[3][i][j] - WE_p(i, j+K-1, j-1));
                CALC[i][j] -= (QUAD[1][i][j-1+ QUAD[2][i][j-1- WE_p(i, j-1, j-K-1));
            }
        }
    }
    ll ans = 0;
    for(i=1 ; i<=N+2*K ; i++)
        for(j=1 ; j<=N+2*K ; j++)
            ans = max(ans, CALC[i][j]);
    cout << ans << endl;
}
 
int main(void)
{
    fio; ll i, tc; cin >> tc;
    for(i=1 ; i<=tc ; i++
    {
        cout << "Case #" << i << endl;
        solve(); // endl in solve
    }
    return 0;
}
cs

 

4번 : 패턴 매칭

부분문제 2까지는 최근 IOI 선발고사에서 나온 문제와 같고, 이 아이디어를 확장하면 문제를 해결할 수 있다.

 

저 동치 조건을 깔끔하게 나타낼 방법을 찾는 것이 이 문제를 해결하는 핵심이다.

결론부터 말하자면, 각 문자열의 각 문자에 대해서 "마지막으로 그 문자가 나온 위치까지의 거리"를 생각하면 동치 조건이 깔끔해진다. 예를 들어, 

 

superguesser의 경우 순서대로 마지막으로 그 문자가 나온 위치까지 거리가 -1, -1, -1, -1, -1, -1, 5, 4, 8, 1, 3, 7가 된다.

abcdefbdaade도 역시 순서대로 마지막으로 그 문자가 나온 위치까지 거리가 -1, -1, -1, -1, -1, -1, 5, 4, 8, 1, 3, 7가 된다.

 

단, -1은 이전에 그 문자가 나오지 않았음을 의미한다.  

또한, 이 두 문자열은 동치임을 직접 확인할 수 있으며, 일반적으로도 이렇게 두 문자열의 동치 여부를 확인할 수 있음을 증명할 수 있다.

 

이제 가장 "자연스러운" 접근은, 문자열을 위와 같이 숫자 형태로 전환시킨 다음, KMP를 쓰는 것이다. 

이 접근은 다 좋은데, -1을 처리하는 것에 약간의 신경을 써야 한다. 예를 들어, KMP를 쓰면 기본적으로 하는 접근이 $[i-fail[i]+1, i]$가 전체 문자열의 prefix와 같다는 것이다. 여기서 다음 문자를 확인하여 이 prefix를 연장시킬 수 있는지 확인해야 한다. 만약 실제 prefix에서 $fail[i]$번째 문자에 대응되는 값이 -1이라면, 이는 이 문자가 prefix에서 지금까지 등장하지 않은 문자임을 의미한다. 이를 확인하기 위해서는 단순히 $i+1$번째 문자에 대응되는 값이 -1임을 확인하면 안되고, 그 값이 -1이거나 $fail[i]$를 초과하는지 확인해야 한다. 즉, KMP에서 지금 확인하고 있는 suffix의 범위를 넘어가는 경우를 역시 고려해주어야 한다.

 

여기까지 생각하면 KMP 코드를 거의 그대로 가져와서 부분문제 2까지 해결할 수 있다.

 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
#include <bits/stdc++.h>
#define fio ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
using namespace std;
typedef long long int ll;
typedef unsigned long long int ull;
typedef long double ldb;
const ll mod = 1e9 + 7;
 
string s; ll N, K;
string pat[31111];
ll pogs[2222222];
ll fail[511], pog[511];
ll rec_s[26], rec_p[26];
 
void precompute(void)
{
    ll i;
    for(i=0 ; i<26 ; i++) rec_s[i] = -1;
    for(i=0 ; i<s.length() ; i++)
    {
        int cur = s[i] - 'a';
        if(rec_s[cur] == -1) pogs[i] = -1;
        else pogs[i] = i - rec_s[cur];
        rec_s[cur] = i;
    }
}
 
ll calc(ll idx)
{
    ll i, j=0
    for(i=0 ; i<26 ; i++) rec_p[i] = -1;
    memset(fail, 0sizeof(fail));
    for(i=0 ; i<pat[idx].length() ; i++)
    {   
        int cur = pat[idx][i] - 'a';
        if(rec_p[cur] == -1) pog[i] = -1;
        else pog[i] = i - rec_p[cur];
        rec_p[cur] = i;
    }
    for(i=1 ; i<pat[idx].length() ; i++)
    {
        while(1)
        {
            if(j == 0break;
            if(pog[i] == pog[j]) break;
            if(pog[j] == -1 && pog[i] > j) break;
            j = fail[j-1];
        }
        if((pog[i] == pog[j]) || (pog[j] == -1 && pog[i] > j)) fail[i]=++j;
    }
    ll ret=0; j=0;
    for(i=0 ; i<s.length() ; i++)
    {
        while(j>0 && !(pogs[i] == pog[j] || (pog[j] == -1 && pogs[i] > j))) j=fail[j-1];
        if((pogs[i] == pog[j]) || (pog[j] == -1 && pogs[i] > j))
        {
            if(j==pat[idx].length()-1) { ret++; j=fail[j]; }
            else j++;
        }
    }
    return ret;
}
 
void solve(void)
{
    cin >> N >> K; ll i, ans = 0cin >> s; precompute();
    for(i=1 ; i<=K ; i++cin >> pat[i];
    for(i=1 ; i<=K ; i++) ans += i * calc(i);
    cout << ans << endl;
}
 
int main(void)
{
    fio; ll i, tc; cin >> tc;
    for(i=1 ; i<=tc ; i++
    {
        cout << "Case #" << i << endl;
        solve(); // endl in solve
    }
    return 0;
}
cs

 

이제 패턴이 여러 개 존재하니, KMP를 Aho-Corasick으로 바꾸어주면 된다. Trie의 각 노드에는

  • 기본적인 정보인 failure link, output link, count, 끝 정점인지 여부 
  • 현재 보고 있는 문자열의 길이 (-1을 처리하기 위함이다)
  • 그리고 각 경우에 대응되는 다음 노드의 포인터

를 저장한다. 이때, 다음 노드를 확인하기 위해 사용하는 것은 알파벳 자체가 아니라 해당 알파벳의 마지막 등장 위치까지의 거리다. 

각 노드에 저장되어 있는 문자열의 길이 정보에 따라서, 내가 -1을 사용해야 하는지 실제 등장 위치까지의 거리를 사용해야 하는지가 달라진다.  

 

위 KMP 풀이와 Aho-Corasick 알고리즘을 잘 이해하고 있다면, 풀이를 변형하는 것은 어렵지 않다. 코드는 https://blog.myungwoo.kr/101를 참고했다. 

 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
#include <bits/stdc++.h>
#define fio ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
using namespace std;
typedef long long int ll;
typedef unsigned long long int ull;
typedef long double ldb;
const ll mod = 1e9 + 7;
 
string s; ll N, K;
string pat[31111];
struct Trie
{
    map<ll, Trie*> S;
    Trie* fail;
    Trie* output;
    ll cur_len;
    ll count; bool is_end;
    Trie() { S.clear(); fail=nullptr; output=nullptr; count=0; cur_len=0; is_end=false; }
    ~Trie() { S.clear(); }
};
Trie* end_node[31111];
queue<Trie*> Q;
vector<Trie*> ord;
 
ll rec_p[26], rec_s[26];
ll pog[31111], pogs[2222222];
map<ll, Trie*>::iterator it;
 
void solve(void)
{
    ord.clear(); cin >> N >> K; ll i, j, ans = 0cin >> s;
    memset(pog, 0sizeof(pog));
    memset(pogs, 0sizeof(pogs));
    for(i=1 ; i<=K ; i++cin >> pat[i];
    Trie *root = new Trie;
    // Step 1 : Build the Trie
    for(i=1 ; i<=K ; i++)
    {
        Trie *now = root;
        for(j=0 ; j<26 ; j++) rec_p[j] = -1;
        for(j=0 ; j<pat[i].length() ; j++)
        {   
            int cur = pat[i][j] - 'a';
            if(rec_p[cur] == -1) pog[j] = -1;
            else pog[j] = j - rec_p[cur];
            rec_p[cur] = j;
        }
        for(j=0 ; j<pat[i].length() ; j++)
        {
            if(now->S.find(pog[j]) == now->S.end()) {
                now->S[pog[j]] = new Trie;
            }
            now->S[pog[j]]->cur_len = now->cur_len + 1;
            now = now->S[pog[j]];
        }
        now->is_end = true;
        end_node[i] = now;
    }
    // Queue Setup
    for(it = root->S.begin() ; it != root->S.end() ; it++) {
        it->second->fail = root;
        Q.push(it->second);
    }
    // Fail Setup
    while(!Q.empty())
    {
        Trie *cur = Q.front(); Q.pop();
        if(cur->is_end) cur->output = cur;
        else cur->output = cur->fail->output;
        for(it = cur->S.begin() ; it != cur->S.end() ; it++) {
            ll cur_step = it->first;
            Trie *nxt = it->second;
            nxt->fail = cur->fail;
            while(1)
            {
                if(nxt->fail == root) break;
                ll len = nxt->fail->cur_len;
                ll true_step = cur_step;
                if(true_step > len) true_step = -1;
                if(nxt->fail->S.find(true_step) != nxt->fail->S.end()) break;
                nxt->fail = nxt->fail->fail;
            }
            ll true_step = cur_step;
            ll len = nxt->fail->cur_len;
            if(true_step > len) true_step = -1;
            if(nxt->fail->S.find(true_step) != nxt->fail->S.end()) 
                nxt->fail = nxt->fail->S[true_step];
            Q.push(nxt);
        }
    }
    // start finding
    Trie *now = root;
    for(i=0 ; i<26 ; i++) rec_s[i] = -1;
    for(i=0 ; i<s.length() ; i++)
    {
        int cur = s[i] - 'a';
        if(rec_s[cur] == -1) pogs[i] = -1;
        else pogs[i] = i - rec_s[cur];
        rec_s[cur] = i;
    }
    for(i=0 ; i<s.length() ; i++)
    {
        ll cur_step = pogs[i];
        while(1)
        {
            if(now == root) break;
            ll true_step = cur_step;
            if(true_step > now->cur_len) true_step = -1;
            if(now->S.find(true_step) != now->S.end()) break;
            now = now->fail;
        }
        ll true_step = cur_step;
        if(true_step > now->cur_len) true_step = -1;
        if(now->S.find(true_step) != now->S.end()) now=now->S[true_step];
        if(now->output) now->output->count++;
    }
    // finish
    Q.push(root);
    while(!Q.empty())
    {
        Trie *cur=Q.front(); Q.pop();
        ord.push_back(cur);
        for(it = cur->S.begin() ; it != cur->S.end() ; it++) Q.push(it->second);
    }
    reverse(ord.begin(), ord.end());
    for(i=0 ; i<ord.size() ; i++)
        if(ord[i]->is_end && ord[i]->fail->output)
            ord[i]->fail->output->count+=ord[i]->count;
    for(i=1 ; i<=K ; i++) ans += i * end_node[i]->count;
    delete root;
    cout << ans << endl;
}
 
int main(void)
{
    fio; ll i, tc; cin >> tc;
    for(i=1 ; i<=tc ; i++
    {
        cout << "Case #" << i << endl;
        solve(); // endl in solve
    }
    return 0;
}
cs

 

'PS > 대회 후기' 카테고리의 다른 글

SCPC 2021 1차 예선 풀이  (0) 2021.07.16
ACM-ICPC Seoul Regional 2020  (0) 2020.11.20
SCPC 2020 본선 후기  (5) 2020.11.09
ACM-ICPC Seoul Regional Preliminary 2020  (0) 2020.10.19
SCPC 2020 2차 예선 풀이  (1) 2020.09.05

대충 1시간 30분 정도 쓴 것 같다.

 

1번. 친구들

친구 관계가 일부 주어졌고, 친구의 친구 역시 친구로 간주할 때 "친구 관계인 maximal 그룹"의 개수를 구하는 문제다.

이는 사람을 정점으로, 친구 관계를 간선으로 볼 때 연결 컴포넌트의 개수를 구하는 문제와 같다. disjoint set 자료구조로 쉽게 해결할 수 있다.

 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
#include <bits/stdc++.h>
#define fio ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
using namespace std;
typedef long long int ll;
typedef unsigned long long int ull;
typedef long double ldb;
const ll mod = 1e9 + 7;
 
int n;
int d[111111];
int par[111111];
 
int find(int x)
{
    if(par[x] == x) return x;
    return par[x] = find(par[x]);
}
 
void unite(int u, int v)
{
    u = find(u);
    v = find(v);
    if(u == v) return;
    par[u] = v;
}
 
void solve(void)
{
    int i, ans = 0cin >> n;
    for(i=1 ; i<=n ; i++cin >> d[i];
    for(i=1 ; i<=n ; i++) par[i] = i;
    for(i=1 ; i<=n ; i++)
        if(i+d[i]<=n) unite(i, i+d[i]);
    for(i=1 ; i<=n ; i++if(i == find(i)) ans++;
    cout << ans;
}
 
int main(void)
{
    fio; int tc; cin >> tc;
    for(int i=1 ; i<=tc ; i++
    {
        cout << "Case #" << i << endl;
        solve(); cout << endl;
    }
    return 0;
}
cs

 

2번. 이진수

$b_i = a_{i-t} | a_{i+t}$로 정의된 bit string $b$가 주어졌을 때, 사전순으로 가장 앞선 $a$를 구하는 문제다.

단, index가 가능한 범위 밖인 경우 값은 $0$이라고 생각하도록 한다. 또한, 조건을 만족하는 $a$가 존재함은 보장된다. 

 

먼저 $b_i = 0$이면 $a_{i-t} = a_{i+t} = 0$임이 강제되며, 이렇게 둘 경우 $b_i = 0$ 역시 만족함을 알 수 있다. 

이에 해당되는 $a$의 bit에 0을 넣는 것을 확정하고 나면, 이제 $b_i = 1$ 형태의 조건들만 남았다. 

이들은 결국 $a$의 두 bit 중 적어도 하나가 1이라는 의미고, 이는 직관적으로 봤을 때 1을 더 넣으면 넣을수록 만족하기 쉬운 조건이다.

그러니, 앞서 0으로 확정된 bit가 아닌 임의의 bit에 대해서는 그냥 1로 두면 조건을 만족하는 $a$가 될 것이다.

 

이제 문제는 사전순으로 가장 앞선 $a$를 구하는 건데, 이런 문제가 그렇듯 앞부터 보면서 0을 넣을 수 있는지 판별하면 된다.

$a_i = 0$이 가능한지 고려해보자. 먼저 $a_i = 0$이 앞서 확정되었다면 굳이 고려를 할 필요도 없다.

만약 $i-t$가 범위에 있고 $b_{i-t} = 1$이라면, $a_{i-2t}$ 또는 $a_i$가 $1$이 되어야 한다. $a_{i-2t}=0$이라면, 여기서 $a_i = 1$이 강제된다. 

만약 $i+t$가 범위에 있고 $b_{i+t}=1$이라면, $a_i$ 또는 $a_{i+2t}$가 $1$이 되어야 한다. $a_{i+2t}$가 이미 $0$으로 확정되었다면, 여기서 $a_i = 1$이 강제된다.

이를 제외한 경우에서는, $a_i = 0$을 해도 문제가 되지 않는다. $a_i$가 영향을 주는 $b$의 값은 $b_{i-t}$, $b_{i+t}$이며, 이들의 값을 $a_i = 0$으로 두면서도 정확하게 맞춰줄 수 있기 때문이다. 특히, 마음에 걸릴 수 있는 유일한 부분이 $a_i = 0$으로 두면서 후에 $a_{i+2t}$를 $1$로 둔 것인데, 애초에

  • $a_i = 1$로 두면 무조건 $a_i = 0$으로 두는 것보다 사전순으로 뒤고, 그러니 $a_i = 0$인게 이득이고
  • $a_{i+2t}$는 확정되지 않은 bit이므로 $a_{i+2t} = 1$으로 두었다고 해서 조건을 만족하는 게 "어려워지지 않는다"

애초에 이 경우에는 $a_i = 0$으로 두고 확정된 bit를 제외한 뒤의 모든 bit를 1로 두면 조건을 만족하게 되어있다.

어쨌든 이런 식으로 bit를 계산하면 $\mathcal{O}(n)$에 모든 계산을 할 수 있다. 2번 문제인 것 치고는 약간 까다로웠다. 

 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
#include <bits/stdc++.h>
#define fio ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
using namespace std;
typedef long long int ll;
typedef unsigned long long int ull;
typedef long double ldb;
const ll mod = 1e9 + 7;
 
// b[i] = a[i-t] | a[i+t]
 
int n, t, a[55555], b[55555];
string s;
 
void solve(void)
{
    int i; cin >> n >> t >> s;
    for(i=1 ; i<=n ; i++) b[i] = s[i-1- '0';
    for(i=1 ; i<=n ; i++) a[i] = -1;
    for(i=1 ; i<=n ; i++)
    {
        if(b[i] == 0)
        {
            if(i+t<=n) a[i+t] = 0;
            if(i-t>=1) a[i-t] = 0;
        }
    }
    for(i=1 ; i<=n ; i++)
    {
        if(a[i] != -1continue;
        if(i-t>=1 && (i<=2*|| a[i-2*t] == 0)) a[i] = 1;
        else if(i+t<=&& (i+2*t>|| a[i+2*t]==0)) a[i] = 1;
        else a[i] = 0;
    }
    for(i=1 ; i<=n ; i++cout << a[i];
}
 
int main(void)
{
    fio; int tc; cin >> tc;
    for(int i=1 ; i<=tc ; i++
    {
        cout << "Case #" << i << endl;
        solve(); cout << endl;
    }
    return 0;
}
cs

 

3번. No Cycle

사이클이 생기지 않도록 간선의 방향을 결정하는데, 사전순 조건이 있다. 즉, 앞쪽에서 입력받은 간선들은 최대한 입력받은 순서를 유지하도록 해야한다. 

이 문제 역시 사전순 문제에서 항상 등장하는 방법인, 앞부터 보면서 답을 맞춰나가는 방법을 사용한다. 결국 우리가 풀어야하는 문제는 

  • 방향 그래프 $G$와 지금 당장 방향을 정해야 하는 간선 $e$, 그리고 뒤에 추가할 무방향 간선의 배열 $E$가 주어져있다.
  • $e$를 입력받은 순서로 방향을 줄 때, $E$의 간선들의 방향을 잘 설정하여 $e, E$를 전부 추가한 뒤에도 그래프에 사이클이 없도록 할 수 있는가? 

이 문제를 해결할 수 있다면, 본 문제 역시 위 문제를 각 간선에 대하여 해결하는 것을 반복하여 해결할 수 있다.

즉, 입력으로 방향 그래프 $G$와 무방향 간선 $e_1, e_2, \cdots, e_k$가 들어온다면,

  • $G, e_1, [e_2, \cdots, e_k]$에 대해 문제를 해결. 그 후 결정된 방향에 맞게 $G_1 \gets G + \vec{e_1}$.
  • $G_1, e_2, [e_3, \cdots , e_k]$에 대해 문제를 해결. 그 후 결정된 방향에 맞게 $G_2 \gets G + \vec{e_2}$.
  • 반복하면서 최종적으로 $G_{k-1}, e_k, []$에 대하여 문제를 해결하는 것으로 마무리.

물론, 위 과정에서 각 간선은 그리디하게 가능하면 입력 순서대로 방향을 잡아야 한다. 

 

Claim : $G$가 DAG고, $E$가 추가할 무방향 간선들이라고 하자. $E$의 간선들의 방향을 잘 정하여 $E$를 추가한 뒤에도 $G$가 DAG이도록 할 수 있다.

Proof of Claim : $G$를 위상정렬하고, 위상정렬 순서로 방향을 주면 된다. 증명 끝.

 

결국 위에서 논의한 $G, e, E$에 대한 문제는 사실 $G + e$가 DAG인지, 즉 사이클이 없는지 판별하는 문제와 같다.

$G$ 역시 DAG이므로, $G + e$에 사이클이 있는지 판별하는 문제는 $G$에 특정 경로가 존재하는지 판별하는 문제와 같다.

즉, $e$가 $u$에서 $v$로 가는 간선이라면, $G + e$의 사이클 존재 유무는 $G$에서 $v$에서 $u$로 가는 경로의 존재 여부와 같다.

 

결론적으로, 간선 $u \rightarrow v$가 주어졌다면 현재 그래프에서 $v \rightarrow u$ 경로가 존재하는지 보고,

  • 존재한다면 순서를 뒤집어서 $v \rightarrow u$로 두고 이 간선을 그래프에 추가
  • 그렇지 않다면 순서를 그대로 두고 간선을 그래프에 추가

하는 것을 반복하면 된다. 대충 $\mathcal{O}(K(N+M+K))$ 정도에 돌 것이다. 아마 자료구조 빡세게 쓰면 subquadratic 될듯? 아니면 말고

 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
#include <bits/stdc++.h>
#define fio ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
using namespace std;
typedef long long int ll;
typedef unsigned long long int ull;
typedef long double ldb;
const ll mod = 1e9 + 7;
 
vector<int> edge[511];
pair<intint> E[4333];
queue<int> Q;
int vis[555];
int n, m, k;
 
int work(int u, int v)
{
    int i; for(i=1 ; i<=n ; i++) vis[i] = 0;
    vis[v] = 1; Q.push(v); 
    while(!Q.empty())
    {
        int c = Q.front(); Q.pop();
        for(i=0 ; i<edge[c].size() ; i++)
        {
            if(vis[edge[c][i]] == 0)
            {
                vis[edge[c][i]] = 1;
                Q.push(edge[c][i]);
            }
        }
    }
    if(vis[u]) return 1;
    return 0;
}
 
void solve(void)
{
    int i, u, v;
    cin >> n >> m >> k;
    for(i=1 ; i<=n ; i++) edge[i].clear();
    for(i=1 ; i<=m ; i++)
    {
        cin >> u >> v;
        edge[u].push_back(v);
    }
    for(i=1 ; i<=k ; i++)
    {
        cin >> u >> v;
        E[i] = make_pair(u, v);
    }
    for(i=1 ; i<=k ; i++)
    {
        u = E[i].first;
        v = E[i].second;
        int res = work(u, v);
        cout << res;
        if(res == 0) edge[u].push_back(v);
        else edge[v].push_back(u);
    }
}
 
 
int main(void)
{
    fio; int tc; cin >> tc;
    for(int i=1 ; i<=tc ; i++
    {
        cout << "Case #" << i << endl;
        solve(); cout << endl;
    }
    return 0;
}
cs

 

4번. 예약 시스템

$2 \times m$ 직사각형에 사람들을 배치한다. 사람들은 총 $n$개의 그룹으로 나뉘어 있으며, 각 그룹에는 5명 이상이 있다. 

각 사람들은 스트레스 지수 $w$를 갖는데, 자신의 이웃에 자신과 다른 그룹에 속하는 사람이 있다면 그런 사람의 수와 $w$의 곱만큼 스트레스가 오른다.

우리의 목표는 스트레스 지수의 총합을 최소화하는 것이다. 제한은 전체적으로 큰 편이다. 

 

동적계획법을 쓰기 어렵게 생긴 제한인만큼, 한방에 문제를 푸는 방법이 필요해보인다. 

 

그룹 $S$ 하나를 고정해보자. $S$의 사람들이 스트레스를 느낄 때는, $S$와 $S^C$가 한 변을 공유할 때다.

그러니까 $S$의 둘레의 (단, 직사각형의 둘레는 제외한다) 길이가 길면 길수록 $S$ 전체의 스트레스가 커진다.

이는 $S$의 사람들이 뭉쳐있어야 스트레스가 적을 것 같다는 직관과도 거의 같은 결과다. 

 

"$S$의 둘레"를 생각해보면, 전체 직사각형의 길이 2의 세로변을 잠깐 무시하고 생각해보면

  • $S$의 크기가 짝수인 경우, 직사각형 형태를 이루도록 하면 둘레가 4가 된다.
  • 특히, 이 경우 $S$에서 서로 다른 4명의 사람들이 정확히 $S^C$의 한 사람과 만나게 된다. 이 사람들은 스트레스 지수가 작을수록 좋다.
  • $S$의 크기가 홀수인 경우, 직사각형에 하나 튀어나온 형태를 만들면 둘레가 5가 된다.
  • 특히, 이 경우 $S$에서 서로 다른 4명의 사람들이 $S^C$의 사람과 만나게 된다. 4명 중 정확히 한 명은 $S^C$를 두 명 만나고, 나머지는 한 명 만난다.
  • 물론, 두 명 만나는 사람의 스트레스 지수가 가장 작아야하며, 나머지의 스트레스 지수 역시 작을수록 좋다.

그리고 조금 열심히 생각해보면 이게 진짜 최적임을 납득할 수 있다. 또한, 홀수 크기의 그룹이 짝수개 존재하므로, 홀수 크기의 그룹끼리 서로 맞닿아 직사각형을 이루도록 하면 배치 역시 어렵지 않다. 이제 문제는 길이 2의 세로변을 처리하는 것이다. 길이 2의 세로변에 걸쳐서 $S$를 배치한다면, $S$에서 스트레스를 받는 4명 중 2명은 (이들은 4명 중에서 스트레스가 많은 두 명일 것이다) 전체 직사각형의 세로변의 도움을 받아 스트레스를 받지 않아도 된다. 

 

단순하게 생각하면, 여기서도 "가장 이득을 보는 그룹"을 양쪽 끝에 배치하여 세로변의 도움을 받을 수 있도록 하면 된다.

그런데 이게 이렇게 단순하지는 않은 게, 배치를 양쪽 끝에 하면 앞서 언급한 최적의 배치를 만드는 것이 어려울 수 있기 때문이다.

홀수 크기의 그룹이 정확히 2개 있고, 짝수 크기의 그룹이 여러 개 있다고 하자. 만약 홀수 크기의 그룹 2개가 양쪽 끝에 놓인다면, 짝수 크기 그룹을 직사각형 형태로 사용할 수 없고, 직사각형에서 양쪽에 한 개씩 튀어나온 형태로 사용해야 한다. 이때 짝수 크기 그룹에서 희생을 더 해야한다. 나머지 경우에서는 어떤 크기의 그룹을 양쪽에 배치해도, 우리가 생각하는 "최적의 배치"가 (직사각형 최대한 사용, 홀수 크기 2개 직사각형으로 만들기) 가능하다. 이제 열심히 케이스 분류.

 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
#include <bits/stdc++.h>
#define fio ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
using namespace std;
typedef long long int ll;
typedef unsigned long long int ull;
typedef long double ldb;
const ll mod = 1e9 + 7;
 
 
ll n, m;
vector<ll> even_save;
vector<ll> odd_save;
 
void solve(void)
{
    int i, j; 
    ll sz, val, ans=0, subt=0, odd=0, even=0, saved, ex = 0
    cin >> n >> m;
    vector<ll> cc;
    for(i=1 ; i<=n ; i++)
    {  
        cin >> sz; cc.clear(); cc.resize(sz);
        for(j=0 ; j<sz ; j++cin >> cc[j];
        sort(cc.begin(), cc.end());
        if(!(sz & 1))
        {
            val = cc[0+ cc[1+ cc[2+ cc[3];
            saved = cc[2+ cc[3];
            ex += cc[0+ cc[1];
            ans += val; even++;
            even_save.push_back(saved);
        }
        else
        {
             val = 2 * cc[0+ cc[1+ cc[2+ cc[3];
            saved = cc[2+ cc[3]; 
            ans += val; odd++
            odd_save.push_back(saved);   
        }
    }
    sort(even_save.begin(), even_save.end());
    reverse(even_save.begin(), even_save.end());
    sort(odd_save.begin(), odd_save.end());
    reverse(odd_save.begin(), odd_save.end());
    if(even >= 2) subt = max(subt, even_save[0+ even_save[1]);
    if(even >= 1 && odd >= 1) subt = max(subt, even_save[0+ odd_save[0]);
    if(odd >= 2 && (odd >= 4 || n == 2)) subt = max(subt, odd_save[0+ odd_save[1]);
    if(odd == 2 && n >= 3)
        subt = max(subt, odd_save[0+ odd_save[1- ex);
    even_save.clear(); odd_save.clear();
    cout << ans - subt;
}
 
int main(void)
{
    fio; int tc; cin >> tc;
    for(int i=1 ; i<=tc ; i++
    {
        cout << "Case #" << i << endl;
        solve(); cout << endl;
    }
    return 0;
}
cs

 

5번. 차이

주어진 조건을 두 변수 사이의 간선이라고 보면, 연결 컴포넌트에 안에서는 모든 차이 계산이 가능하다. 

문제가 되는 것은 사이클, 특히 한바퀴 돌았는데 합이 0이 되지 않는 사이클이다. 이 경우, 사이클에 도달할 수 있는 모든 정점들은 모순된 상태가 된다.

 

그러니, disjoint set 자료구조를 가져왔을 때 대강

  • 두 정점이 같은 disjoint set에 속하지 않으면 Not Connected
  • 두 정점이 같은 disjoint set에 속하는데 여기에 모순된 cycle이 있으면 Conflict
  • 두 정점이 같은 disjoint set에 속하고 모순도 없으면 계산이 가능

함을 알 수 있다. 연결관계 자체만 봤을 때 disjoint set은 자명하고, 모순을 찾거나 계산을 하는 것이 문제다.

 

계산을 하기 위해서, 각 disjoint set의 대표 노드를 하나 두고 각 원소가 그 대표 노드의 값보다 얼마나 큰지 저장한다. 

원소들 사이의 차이는 알지만 원소 자체의 값은 모르니, 대표 노드의 값을 하나의 변수로 보고 차이만 저장하자는 의미다. 

 

이제 간선, 즉 하나의 등식이 추가되었다고 가정하자. 

  • 만약 두 정점이 같은 disjoint set에 있다면, 모순이 있는지 바로 확인할 수 있다. 저장된 값의 차이가 입력된 값과 다르면 모순이다.
  • 두 정점이 다른 disjoint set에 있다면, 이번에 추가된 등식은 두 disjoint set의 대표 노드 값 사이의 차이를 제공하는 정보가 된다. 이 경우, 두 disjoint set 중 하나의 대표를 전체의 새로운 대표로 정한다. 그 후, 흡수되는 집합의 각 원소들에 대해서 새로운 대표의 값과의 차이를 계산해야 한다. 이는 새로 얻어진 등식을 활용하여 할 수 있을 것이다. 물론, 두 disjoint set 중 이미 모순된 집합이 있다면 새로운 집합 역시 모순이 있다. 

두 disjoint set을 합치는 과정에서 걸리는 시간이 흡수되는 집합의 크기에 비례함을 알 수 있다. 그러므로, 작은 집합을 큰 집합에 흡수시키는 small-to-large 알고리즘을 사용하면 빠른 시간에 문제를 해결할 수 있다. 또한, 각 집합을 std::set 자료구조로 저장해야 집합의 순회가 쉽다. 

 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
#include <bits/stdc++.h>
#define fio ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
using namespace std;
typedef long long int ll;
typedef unsigned long long int ull;
typedef long double ldb;
const ll mod = 1e9 + 7;
 
ll ans[111111];
int par[111111];
int sz[111111];
int doom[111111];
set<int> idx[111111];
set<int>::iterator it;
ll n, k;
 
int find(int x)
{
    if(par[x] == x) return x;
    return par[x] = find(par[x]);
}
 
void unite(int uu, int vv, int u, int v, ll val)
{
    if(sz[uu] < sz[vv])
    {
        swap(uu, vv);
        swap(u, v);
        val = -val;
    }
    sz[uu] += sz[vv];
    doom[uu] |= doom[vv];
    ll dif = -val - ans[v] + ans[u];
    for(it = idx[vv].begin() ; it != idx[vv].end() ; it++)
    {
        int loc = (*it);
        ans[loc] += dif;
        par[loc] = uu;
        idx[uu].insert(loc);
    }
    idx[vv].clear();
}
 
void solve(void)
{
    int i, whi, u, v, uu, vv; ll val;
    cin >> n >> k;
    for(i=1 ; i<=n ; i++) idx[i].clear();
    for(i=1 ; i<=n ; i++) ans[i] = 0, par[i] = i, sz[i] = 1, doom[i] = 0;
    for(i=1 ; i<=n ; i++) idx[i].insert(i);
    for(i=1 ; i<=k ; i++)
    {
        cin >> whi >> u >> v;
        if(whi == 1)
        {
            cin >> val;
            uu = find(u);
            vv = find(v);
            if(uu == vv && ans[u] - ans[v] != val) doom[uu] = 1;
            if(uu != vv) unite(uu, vv, u, v, val);
        }
        else if(whi == 2)
        {
            uu = find(u);
            vv = find(v);
            if(uu != vv) cout << "NC\n";
            if(uu == vv && doom[uu] == 1cout << "CF\n";
            if(uu == vv && doom[uu] == 0cout << ans[u] - ans[v] << "\n";
        }
    }
}
 
int main(void)
{
    fio; int tc; cin >> tc;
    for(int i=1 ; i<=tc ; i++
    {
        cout << "Case #" << i << endl;
        solve(); cout << endl;
    }
    return 0;
}
cs

'PS > 대회 후기' 카테고리의 다른 글

SCPC 2021 2차 예선 풀이  (1) 2021.08.07
ACM-ICPC Seoul Regional 2020  (0) 2020.11.20
SCPC 2020 본선 후기  (5) 2020.11.09
ACM-ICPC Seoul Regional Preliminary 2020  (0) 2020.10.19
SCPC 2020 2차 예선 풀이  (1) 2020.09.05

600문제는 아마 영원히 못 풀지 않을까?

'PS > 수학 계열 PS' 카테고리의 다른 글

8월의 PS 일지 - Part 5  (0) 2020.08.12
8월의 PS 일지 - Part 3  (0) 2020.08.10
8월의 PS 일지 - Part 1  (0) 2020.08.04
Brief Introduction to Hackenbush Games  (1) 2020.05.17
5월의 PS 일지 - Part 1  (0) 2020.05.16

11단원의 내용으로 내가 알고 있는 PS/CP 정수론 분야 지식은 (특히, 문제풀이에 직접적인 내용은) 다 정리한 것 같다.

처음에는 그래도 친절하게 설명하려고 노력했는데, 후반부에 가니까 슬슬 다른 일의 우선순위가 높아져서 급하게 작성한 것 같아 아쉽다.

  • 초반 부분 : 이미 자료가 있으니까 간략하게 써야지 ㅋㅋㅋ
  • 후반 부분 : 내가 다른 일이 있으니까 간략하게 써야지 ㅋㅋㅋ

가 된 것 같아서 조금 그렇긴 한데 ㅋㅋ; 그래도 시간투자를 충분히 하면 이해할 수 있을 정도로는 쓴 것 같다. 나중에 코드, 연습문제도 올릴거고..


일단 이 내용을 다 이해했다면 (또는 관심있는 부분까지 다 이해했다면) 할 수 있는 것은 대강

  • 백준 문제풀기. 문제집 제작을 내가 하면 더 제대로 할 수 있을 것이고, 아니어도 solved.ac나 이 블로그 수학 PS일지를 보면서 문제를 고를 수 있다.
  • 진짜 고수들이 쓴 자료 읽기. 여러 의미를 담고 있는 말인데, 정리하자면
  • 1. 정수론 교과서 읽기. 이유를 설명하지 않고 넘긴게 많으니, 이를 이해하고 싶다면 책을 사서 읽는 것이 제일 빠르다.
  • 2. PS/CP 정수론 진짜 고수가 쓴 자료 읽기. 예로, kipa00님의 NTA, HYEA님의 소멤 글, min_25의 블로그, adamant의 블로그
  • adamant의 블로그는 PS 수학 전반에 대한 다양한 내용이 있다. 특히 FFT 및 다항식 처리 관련 자료가 어마어마하다.
  • HYEA님의 글과 min_25의 블로그에는 더 난이도 높은 알고리즘들이 많이 나와있다. 또한, HYEA님의 글 중 min_25의 글 번역본도 있다.
  • kipa00님의 NTA는 수학적 깊이가 상당하고, computational number theory의 여러 부분을 다룬다. 수학하는 사람이면 읽는 게 좋을듯.
  • 혹시 NTA를 다 읽고 컨텐츠가 부족하다고 느낀다면 (ㅋㅋ;) Victor Shoup의 computational number theory 책을 읽자.
  • 3. Project Euler 풀기. 여기에는 고수들이 가득하고, thread 등에서 풀이나 아이디어 공유가 많다. 
  • 특히, 10단원, 11단원의 내용을 확장/강화하는 토론들이 많이 이루어졌다. 대표적인 예시로 baihacker의 블로그를 참고하자. 
  • 개인적으로 중국/일본에 비해 한국에서 Project Euler를 푸는 사람이 적다는 것에 많은 아쉬움을 느끼고 있다.
  • 4. PS/CP scope 벗어나기. 이제 $2^{64}$ 미만의 수만 다루고, 짧은 시간 안에 답을 내야한다는 가정도 버리자.
  • 아예 Computational Number Theory를 공부할 수도 있고, (NTA 읽는 것도 방법) 그 대표적인 활용처인 암호학을 공부하는 것도 좋다.
  • 개인적으로는 PS/CP 정수론 및 PS/CP에서 겪은 문제풀이 경험이 CTF 암호학으로 넘어가는 것에 정말정말 많은 도움을 주었다.


시간이 없어서 글을 제대로 못 쓰기는 했는데 그래도 질문 받을 시간은 있으니 필요할 때 댓글 달아주세요 :)

관련 코드는 github에서 찾아볼 수 있다.


이 글에서는 다음 내용을 다룬다.

  • multiplicative function의 정의, 성질
  • Dirichlet Convolution의 정의와 성질
  • xudyh's sieve, min_25 sieve를 이용한 multiplicative function의 prefix sum

이 내용을 다루는 튜토리얼 중 가장 접근 가능한 것은 역시 이 Codeforces 블로그. 나도 여기서 배웠다.


multiplicative function

  • 기존 정의를 생각하면, multiplicative function이란 $\text{gcd}(a, b)=1$일 때 $f(ab)=f(a)f(b)$인 함수 $f$를 말한다.
  • 특히, 우리가 다룰 함수들은 $f(1)=1$을 만족한다. ($f=0$ 같은 특수한 경우를 고려하지 않겠다는 것이다)
  • 이제 $n$을 소인수분해하고, $n=p_1^{e_1} \cdots p_k^{e_k}$라고 쓰자. 그러면 $f(n) = f(p_1^{e_1}) \cdots f(p_k^{e_k})$가 성립한다.
  • 이는 prime power에 대한 $f$ 값들만 결정하면 모든 자연수에 대한 $f$ 값이 자동 결정이라는 뜻이다.
  • 이 방법으로 생각하는 것이 아마 multiplicative function을 이해하는 가장 쉬운 방법일 것이다. 예시를 보자.
  • 오일러 파이 함수 : $f(p^k) = p^{k-1}(p-1)$으로 정의된 multiplicative function
  • mobius 함수 : $f(p)=-1$, $f(p^k)=  0$ for $k \ge 2$로 정의된 multiplicative function
  • 약수의 개수 함수 : $f(p^k) = k+1$로 정의된 multiplicative function
  • 약수의 합 함수 : $f(p^k) = 1 + p + \cdots + p^k$로 정의된 multiplicative function
  • identity 함수 : $f(p^k) = p^k$로 정의된 multiplicative function
  • $[n=1]$ 함수 : $f(p^k) = 0$으로 정의된 multiplicative function
  • $f=1$ 함수 : $f(p^k) = 1$으로 정의된 multiplicative function
등등, prime power만 보고 모든 값을 얻을 수 있는 함수들을 나열할 수 있을 것이다.

이제 multiplicative function의 여러 성질을 어렵지 않게 증명할 수 있다. 예시를 하나 보인다.

사실. $f(n)$이 multiplicative라면, $g(n) = \sum_{d|n} f(d)$ 역시 multiplicative.
증명. $g(p^k) = 1 + f(p) + \cdots + f(p^k)$인 multiplicative function이 됨을 보일 수 있다. 예를 들어, $$ g(p_1^{e_1} \cdots p_k^{e_k}) = \sum_{d_i \le e_i} f(p_1^{d_1} \cdots p_k^{d_k}) = \sum_{d_i \le e_i} f(p_1^{d_1}) \cdots f(p_k^{d_k}) = \sum_{d_1 \le e_1} f(p_1^{d_1}) \cdot \sum_{d_2 \le e_2} f(p_2^{d_2}) \cdots \sum_{d_k \le e_k} f(p_k^{d_k})$$이고 이는 정의상 $g(p_1^{e_1}) \cdots g(p_k^{e_k})$와 동일하다. 증명 끝. 이제 본격적인 주제로 들어가자.


Dirichlet Convolution

  • Dirichlet Convolution -> Dirichlet Series -> Analytic Number Theory의 연결 관계로 중요한 주제인 것으로 안다.
  • 하지만 PS/CP에서는 Dirichlet Convolution 그 자체가 중요한 경우는 거의 못 봤고, 아래에서 다룰 xudyh's sieve에서 등장한다.
  • 하지만, 실제로 Dirichlet Series와 같은 방식으로 PS/CP 문제풀이를 하는 경우를 Project Euler에서 봤다. 이건 나도 잘 모른다 ㅋㅋ
  • 정의. 두 함수 $f, g$에 대하여, 그 Dirichlet Convolution $(f * g)$를 $(f * g)(n) = \sum_{d|n} f(d)g(n/d)$로 정의한다.
  • 사실. $f, g$가 multiplicative 하다면, 그 Dirichlet Convolution $(f * g)$ 역시 multiplicative 하다.
  • 증명. 또 $(f * g)(p^k) = f(1)g(p^k) + f(p)g(p^{k-1}) + \cdots + f(p^k)g(1)$로 두고, 위에서 보인 예시처럼 하면 된다.
예시를 몇 개 들어보자.
  • $f(n) = n$, $g(n)=1$이면 $(f * g)(n)= \sigma(n) = \sum_{d|n} d$ - 약수의 합
  • $f(n)=1$, $g(n)=1$이면 $(f*g)(n) = \tau(n) = \sum_{d|n} 1$ - 약수의 개수
  • $f(n) = \mu(n)$, $g(n)=1$이면 $(f*g)(n) = \sum_{d|n} \mu(d) = [n=1]$ - mobius function
  • $f(n) = \phi(n)$, $g(n)=1$이면 $(f * g)(n) = \sum_{d|n} \phi(d) = n$
  • $f(n) = \mu(n)$, $g(n) = n$이면 $(f * g)(n) = \sum_{d|n} \mu(d) (n/d) = \phi(n)$
  • 특히 - mobius inversion을 제대로 표현하는 방법은 역시 $h = (f * 1)$이면 $f = (h * \mu)$라는 것이다.
  • $f(n) = n^k \phi(n)$, $g(n) = n^k$면 $(f * g)(n) = n^{2k+1}$ 등.

이제 메인 주제인, multiplicative function의 prefix sum을 계산하는 것으로 넘어가자.

xudyh's sieve
  • multiplicative function $f$가 있다. 우리의 목표는 $\sum_{i=1}^n f(i)$를 $\mathcal{O}(n^{2/3})$에 계산하는 것이다.
  • 적당한 multiplicative function $g$가 있어, $g$의 prefix sum과 $(f * g)$의 prefix sum이 계산하기 쉽다고 가정하자.
  • 함수 $f$에 대하여, $s_f(n) = \sum_{i=1}^n f(i)$라고 정의한다. 이제 준비 끝. 열심히 식을 조작해보자.

우선 Dirichlet Convolution의 정의를 사용하여 $$s_{f * g}(n) = \sum_{i=1}^n (f * g)(i) = \sum_{i=1}^n \sum_{d|i} f(i/d)g(d)$$를 얻는다. 이제 $d$에 대하여 식을 정리하면, 가능한 $i/d$의 값이 $1$부터 $\lfloor n/d \rfloor$이므로 $$ s_{f * g}(n) = \sum_{d=1}^n g(d) \sum_{k=1}^{\lfloor n/d \rfloor} f(k) = \sum_{d=1}^n g(d) s_f(\lfloor n/d \rfloor)$$이다. $g(1)=1$이니, 이제 이 식을 $s_f(n)$에 대하여 정리하면 결과적으로 $$s_f(n) =  s_{f * g}(n) - \sum_{d=2}^n g(d) s_f(\lfloor n/d \rfloor) $$

  • 계속해서 등장하지만, $\lfloor n/d \rfloor$ 형태를 갖는 서로 다른 자연수의 개수는 $\mathcal{O}(\sqrt{n})$개이다.
  • 또한, $\lfloor n/d \rfloor$의 특정 값에 대응되는 $d$의 값은 하나의 구간을 이룬다.
  • 특히, $g$의 prefix sum이 쉽게 계산될 수 있다고 가정했으므로, 그 구간에 대한 $g$ 값의 합 역시 쉽게 구할 수 있다.
  • 여기서는 당연히 $\sum_{i=l}^r g(i) = s_g(r) - s_g(l-1)$이라는 식을 사용한다. 
  • 그러니, 우리가 실제로 계산해야 하는 것은 $s_f(\lfloor n/d \rfloor)$ 값이다.

첫 방식은 $\lfloor n/d \rfloor$ 형태의 값들에 대해 $s_f$를 정직하게 DP로 계산하는 것이다.

$s_f(x)$를 계산하기 위해 드는 시간복잡도가 $\mathcal{O}(\sqrt{x})$이니, 우리에게 필요한 시간은 총 $$ \sum_{i=1}^{\sqrt{n}} \sqrt{i} + \sum_{i=1}^{\sqrt{n}} \sqrt{n/i} = \mathcal{O}(n^{3/4})$$임을 알 수 있다. 이제 지수를 또 $3/4$에서 $2/3$으로 줄여보자. 이를 위해서는 $f$가 multiplicative임을 이용한다.

  • 즉, 에라토스테네스의 체나 Linear Sieve를 사용하여, 작은 $i$들에 대한 $f(i)$를 전부 전처리하자.
  • 그러면 자연스럽게, 부분합을 계산하여 작은 $i$들에 대한 $s_f(i)$ 역시 전부 전처리할 수 있다.
  • 이번에는 $n^{2/3}$ 이하의 $i$에 대한 $s_f(i)$를 전처리하고, $s_f(n)$을 구하기 위한 점화식을 top-down으로 적용하자.
  • 이미 전처리된 값은 바로 반환하고, 메모이제이션을 이용하여 중복 계산을 막는다. 그러면 시간복잡도는 $$ \mathcal{O}(n^{2/3}) + \sum_{i=1}^{n^{1/3}} \sqrt{n/i} = \mathcal{O}(n^{2/3})$$이 된다. 이 방법이 xudyh's sieve이다. 이제 문제는 $s_f$를 구하고 싶을 때, $g$를 어떻게 잡냐는 것이다.

일단 $g$를 잡는 가장 좋은 방법은 위에서 선보인 예시들을 다 참고하는 것이다. 구글 검색하면 예시는 더 나온다.

특히, $1$, $n$, $n^k$ 등 $n$에 대한 다항식으로 나타나는 함수들은 그 prefix sum을 구하기가 상대적으로 쉽다는 점을 이용한다.

어려운 문제들을 풀기 시작하면, $g$를 잡기 위해서 전략적으로 각 prime power들에 대해 $g$의 값을 부여해야 할 수도 있다.

일단 위 문제의 가장 자주 보이는 예시는, $\phi(n)$의 합을 구하기 위하여 $g(n) = 1$, $(f * g) (n) =  n$을 이용하는 것이다.


구현 노트.

  • 구현 자체에 크게 어려운 점은 없으나, 역시 $\lfloor n/d \rfloor$가 매우 클 수 있으니 C++의 std::map을 이용해야 한다.
  • 하지만 std::map을 사용하면 로그가 하나 더 붙는다. 로그를 붙이지 않고 계산을 할 수는 없을까?
  • 애초에 $\lfloor n/d \rfloor \le n^{2/3}$이면 바로 값을 리턴하므로, $\lfloor n/d \rfloor > n^{2/3}$인 경우만 생각하자. 이 경우, $\lfloor n/d \rfloor$ 대신 그냥 $d$를 key 값으로 사용할 수 있다. 
  • 그러면 key 값이 $n^{1/3}$ 이하로 매우 작으니, std::map 대신에 그냥 배열을 사용해도 무방하다. 물론, 이 경우 초기화 문제에 주의해야 할 것이다.


min_25 sieve

  • 우선 필요한 정의부터 나열하고 시작하도록 하자. 모든 설명은 이 글을 따른다.
  • 이제부터 모든 나눗셈은 C++ 방식의 정수 나눗셈이다. 자연수만 다루니 부호는 신경쓰지 말자. (10장에서도 이럴 걸 ㅋㅋ)
  • $p_k$는 $k$번째 소수. 1-index를 사용한다. 즉, $p_1, p_2, p_3, \cdots = 2, 3, 5, \cdots$.
  • $lpf(n)$은 $n$이 갖는 최소의 소인수. 단, $n=1$인 경우 $lpf(1)=1$로 정의한다.
  • $F_{prime}(n)$은 $\sum_{2 \le p \le n} f(p)$이며, 이 합은 소수 $p$에 대한 것이다.
  • $F_k(n)$은 $\sum_{p_k \le lpf(i), 2 \le i \le n} f(i)$이다. 특히, $\sum_{i=1}^n f(i) = F_1(n) + f(1) = F_1(n)+1$이다.
  • 즉, $F_k(n)$은 최소 소인수가 $p_k$ 이상인 자연수에 대한 $f$ 값의 합이다. 이제 본격적으로 시작.

식 조작을 하자. 각 자연수의 최소 소인수가 무엇인가, 그리고 그 소인수를 몇 개 가지고 있는가로 분류하면, $$F_k(n) = \sum_{p_k \le lpf(i), 2 \le i \le n} f(i) = \sum_{k \le i, p_i \le n} \sum_{c \ge 1, p_i^c \le n} f(p_i^c) \left( 1 + F_{i+1}(n/p_i^c) \right) $$과 같다. 즉, $p_i^c$와 $1$의 곱, 또는 최소 소인수가 $p_{i+1}$이면서 $n/p_i^c$ 이하인 자연수의 곱으로 분류한다. 


한편, $n < p_k$면 $F_k(n) = 0$임을 사용하고, 적당히 식을 정리하면 $$ F_k(n) = \left( \sum_{k \le i, p_i^2 \le n} \sum_{c \ge 1, p_i^{c+1} \le n} \left( f(p_i^{c+1}) + f(p_i^c) F_{i+1}(n/p_i^c) \right) \right) + \sum_{k \le i, p_i \le n} f(p_i) $$를 얻는다. 여기서는 기존 식에서 $p_i^2 > n$이면 나오는 항이 $f(p_i)$ 하나임을 이용했다. $F_{prime}$ 표현을 이용하면, $$ F_k(n) = \left( \sum_{k \le i, p_i^2 \le n} \sum_{c \ge 1, p_i^{c+1} \le n} \left( f(p_i^{c+1}) + f(p_i^c) F_{i+1}(n/p_i^c) \right) \right) + F_{prime}(n) - F_{prime}(p_{k-1})$$을 얻는다. 이제 이 점화식을 써먹는 방법을 고민해보자.

  • 위 점화식을 top-down으로 계산한다고 가정하자.
  • 또 $\lfloor n/i \rfloor$ 형태의 값들만 등장한다. 이제 $\mathcal{O}(\sqrt{n})$개의 값들만 등장함은 익숙하다.
  • 재귀 $k$ 값들을 보면, $p_i^2 \le n$을 아는 상태에서 $F_{i+1}(*)$ 값들을 호출함을 볼 수 있다. 
  • 그러니 $F_k(n)$을 구해야 한다는 재귀호출이 나왔다면, $p_{k-1}$이 큰 값은 아니라는 것이다. 
  • 왜 이게 중요하냐면, $F_{prime}(p_{k-1})$의 값이 전처리의 영역에 있다는 것을 설명하기 때문이다. 
  • 즉, $n^{1/2}$ 정도의 값까지 $F_{prime}$의 값을 전처리하면, $F_{prime}(p_{k-1})$ 값은 $\mathcal{O}(1)$에 정리된다.
  • 이제 $\lfloor n/i \rfloor$ 형태의 값들에 대해 $F_{prime}$ 값을 전부 구해야 한다는 문제가 있다.
  • 만약 $f(p)$가 $p$에 대한 다항식이라면, 앞에서 배운 Lucy-Hedgehog 알고리즘으로 이를 전부 전처리할 수 있다.

시간복잡도 분석은 자료에 따르면 $\mathcal{O}(n^{1-\epsilon})$이다. 하지만 시간복잡도에 비해 매우 좋은 성능을 보인다.

이 알고리즘의 강점은 xudyh's sieve처럼 $g$를 고민해서 잡을 필요가 없다는 점, 그리고 $f$에 대해 가정하는 점이 적다는 점이다.


구현 노트.

  • 중요한 사실인데, $p_i < \sqrt{n} < p_{i+1}$이라고 해서 $n/p_i < p_{i+1}$인 것은 아니다.
  • 그러니 전처리를 하는 소수/체의 범위를 조금 조심해서 고를 필요가 있다. 마음 편하게 $\sqrt{n} + 500$ 정도로 두자.


관련 코드는 github에서 찾아볼 수 있다.


이 글에서는 다음 알고리즘을 다룬다.

  • Lucy-Hedgehog의 Algorithm 소개 및 최적화.
  • Meissel-Lehmer Algorithm 소개 및 최적화.


Lucy-Hedgehog Algorithm

  • 글을 읽기 전에, 경건한 마음을 가지고 Project Euler에 가입하자. 그 다음, (쉬운) 10번 문제를 풀자.
  • 그러면 여러 사람이 자신의 접근과 풀이, 코드를 공유하는 thread에 접근할 수 있다.
  • 이제 https://projecteuler.net/thread=10;page=5#111677로 가면, 이 알고리즘의 기원을 볼 수 있다.
  • Note. 동적계획법만 알고 있어도 이 알고리즘을 충분히 이해할 수 있다.


$n$ 이하의 소수의 개수를 $\mathcal{O}(n^{3/4})$에 계산하는 방법을 소개한다. 먼저, $n^{1/2}$ 이하의 소수들을 에라토스테네스의 체를 이용하여 전처리한다. 

이제 $dp[v][m]$을 에라토스테네스 체에서 $m$까지 보았을 때, $1$ 이상 $v$ 이하 자연수 중 지워지지 않은 것의 개수라고 하자.

  • 즉, $dp[v][m]$은 소수거나, 모든 소인수가 $m$ 초과인 $1$ 이상 $v$ 이하 자연수의 개수이다.
  • 특히, $(m+1)^2 > v$라면 모든 소인수가 $m$ 초과인 $v$ 이하 자연수는 무조건 소수이므로, $dp[v][m]$은 $v$ 이하 소수의 개수다.


이제 $dp$ 식을 생각해보자.

  • 만약 $m$이 소수가 아니라면, 애초에 에라토스테네스 체에서 하는 일이 없다. 이때는 $dp[v][m]=dp[v][m-1]$.
  • 만약 $m^2 > v$라면, $m$이 소수여도 이번에 새로 지워지는 값은 없을 것이다. 여전히 $dp[v][m]=dp[v][m-1]$.
  • $m$이 소수라고 가정하자. 이번에 제거되는 값들은 최소 소인수가 $m$인 값들이다.
  • 즉, $m$과 "최소 소인수가 $m-1$ 초과인 자연수"을 곱한 값들이 지워진다. 
  • 이 값의 개수를 구하려면, $\lfloor v/m \rfloor$ 이하 자연수 중 최소 소인수가 $m-1$ 초과인 자연수의 개수를 구하면 된다.
  • 그런데 $dp[\lfloor v/m \rfloor][m-1]$에 소수이거나, 최소 소인수가 $m-1$ 초과인 값의 개수가 저장되어 있다.
  • 이 값에 $m-1$ 이하 소수의 개수를 빼면 원하는 값이 나올 것이다. 그 값은 $dp[m-1][m-1]$에 저장되어 있다.
  • 그러므로, $dp[v][m] = dp[v][m-1] - (dp[\lfloor v/m \rfloor][m-1] - dp[m-1][m-1])$이라는 결과가 나온다.


이제 필요한 $dp$ 상태를 생각해보자. 우리는 $dp[n][\lfloor \sqrt{n} \rfloor]$가 필요하다.

  • 먼저 두 번째 인자 $m$을 보자. 우리의 관찰에 의해, $m^2 > v$인 경우는 추가적으로 고려할 필요가 없다. 
  • 특히, 우리가 다룰 모든 두 번째 인자 $m$ 값은 $\sqrt{n}$ 이하이다.
  • 이제 첫 인자인 $n$을 생각하자. 필요한 첫 번째 인자의 값은 $\lfloor n/i \rfloor$ 형태의 값을 가짐을 알 수 있다.
  • 즉, 필요한 첫 번째 인자의 값은 $1, 2, \cdots, \lfloor \sqrt{n} \rfloor$과 $\lfloor n/1 \rfloor, \lfloor n/2 \rfloor, \cdots , \lfloor n/\lfloor \sqrt{n} \rfloor \rfloor$다.
  • 중요한 점은, 다루어야 할 첫 번째 인자의 값의 종류가 $\mathcal{O}(\sqrt{n})$개라는 점이다. 이제 준비 완료.


이제 본격적으로 알고리즘을 구체화한다.

  • DP에 익숙한 사람들은 바로 눈치를 챘겠지만, 위 DP 식은 정말 토글링을 하기 좋은 식이다. 
  • 그러니 $dp$를 1차원 배열로 정의하고, 순서대로 $m$ 값을 늘려나가보자. $dp$ 배열은 $dp[v][m]$을 갖고 있을 것이다.
  • 그러기 위해서는 사실 $dp[v][1]$을 정의해야 하는데, 이는 시작하자마자 체에서 1을 제거했다고 치고 $dp[v][1]=v-1$로 둔다.
  • 이제 큰 $v$값부터 순서대로 $dp[v][m]=dp[v][m-1]-(dp[\lfloor v/m \rfloor][m-1] - dp[m-1][m-1])$을 적용한다.
  • 그러다가, $v \le m^2$인 시점이 오면 업데이트를 중단하고 다음 $m$으로 넘어가도록 한다. 
  • 물론, 우리가 실제로 업데이트를 해야 하는 $v$의 값들은 앞서 언급한 $\mathcal{O}(\sqrt{n})$ 개의 값들이다.
  • 마지막으로 얻은 $dp[n]$의 값이 우리가 원하는 $n$ 이하의 소수 개수다. 이게 알고리즘 끝이다. 예상보다 간단했을 것.
  • 특히, 실제로 관심이 있던 각 값 $v$에 대하여, $dp[v]$는 $v$ 이하의 소수 개수다. 즉, 예를 들어, $n/2$, $n/3$ 이하의 소수 개수도 동시에 구해진다.


시간복잡도 분석을 하자.

  • 이 알고리즘은 "실제로 계산해야 하는" $dp[v][m]$의 개수만큼 계산을 한다.
  • 특히, 우리는 $m^2 \le v$인 경우에만 추가적인 계산을 하니, 우리가 다루는 각 $v$값에 대해 $\sqrt{v}$를 더한만큼 계산을 한다. 이는 $$ \sum_{i=1}^{\sqrt{n}} \sqrt{i} + \sum_{i=1}^{\sqrt{n}} \sqrt{n/i} = \mathcal{O}(\sqrt{n} \cdot \sqrt{n}^{1/2}) + \mathcal{O}(\sqrt{n} \cdot \sqrt{n}^{1/2}) = \mathcal{O}(n^{3/4})$$다. 여기서 $1+1/\sqrt{2} + 1/\sqrt{3} + \cdots + 1/\sqrt{n} = \mathcal{O}(\sqrt{n})$을 이용하였다.


구현 노트.

  • $dp$ 배열을 관리하기 위하여 가장 쉽게 생각할 수 있는 방법은 역시 C++의 std::map일 것이다. 하지만 로그가 붙는다.
  • 로그를 붙이지 않고 구현을 하려면, $dp$를 말 그대로 배열에 저장해야 한다. 이를 위해서는, $v = 1, 2, \cdots, \lfloor \sqrt{n} \rfloor$에 대응되는 길이 $\lfloor \sqrt{n} \rfloor$의 배열과, $v = \lfloor n/1 \rfloor, \lfloor n/2 \rfloor, \cdots ,  \lfloor n/\lfloor \sqrt{n} \rfloor \rfloor$에 대응되는 길이 $\lfloor \sqrt{n} \rfloor$의 배열을 따로 두면 된다. 자세한 디테일은 github 참고.


추가 코멘트.

  • 필자는 이 알고리즘을 많이 써먹었다. 이 알고리즘의 장점은, 확장이 쉽게 가능하다는 것이다.
  • 예를 들어, 소수의 개수가 아니라 소수의 합, 소수의 제곱의 합 등을 원하는 경우도 $dp$ 식만 살짝 바꿔주면 해결된다.
  • 심지어, $n$ 이하 $4k+1$ 꼴 소수의 개수/합, $4k+3$ 꼴 소수의 개수/합 역시 약간의 아이디어를 추가하면 할 수 있다.
  • 확장이 쉽게 가능한 이유는, 역시 기본 원리인 dp 자체가 간단한 편이기 때문이다. 이해하기도 쉬운 편인 알고리즘이라 생각. 
  • 사실 시간복잡도 분석을 빡빡하게 하지 않았다. 실제로는 $m$이 소수인 경우만 중요하다는 사실도 역시 써먹을 수 있다.
  • 빡빡한 시간복잡도 분석을 위해서는, $n$ 이하 소수의 개수가 약 $n/\ln n$개 존재한다는 Prime Number Theorem을 참고하라.


이 알고리즘을 $\mathcal{O}(n^{2/3})$ 수준으로 최적화시킬 수 있다. 여기서 "로그가 시간복잡도에 얼마나 붙는지"는 따지지 않겠다.

  • 지수를 3/4에서 2/3으로 내리는 것은 11단원에서도 등장한다. 보통은, "에라토스테네스의 체를 더 써먹자"의 방식으로 최적화가 된다.
  • 사실 우리는 기존 알고리즘에서 에라토스테네스의 체로 $1$부터 $n^{1/2}$까지의 소수 판별만 했다. 체를 더 써먹자.

이제부터, 기존 점화식 $dp[v][m] = dp[v][m-1] - (dp[\lfloor v/m \rfloor][m-1] - dp[m-1][m-1])$을 top-down으로 돌려보자.

목표는 실제로 계산을 전부 하는 것이 아니라, 어떤 $dp$ 값들이 필요한지를 탐색하는 것이다. 시작은 물론 $dp[n][\lfloor \sqrt{n} \rfloor]$을 호출하는 것이다.

  • $dp[m-1][m-1]$은 $m$ 미만의 소수 개수와 같다. 이는 에라토스테네스의 체로 처리할 수 있으니, $dp$ 값으로 간주하지 않는다.
  • 만약 $m=1$이라면, $dp[v][m]=v-1$임을 이미 알고 있으니, 이를 필요한 $dp$ 값으로 간주하지 않는다.
  • 만약 $v \le n^{2/3}$이라면, "$dp[v][m]$을 계산해야 한다"라는 정보를 저장하고 리턴한다.
  • 재귀 과정에서 이미 방문한 위치 $v, m$에 다시 도달했다면, 이미 $dp[v][m]$을 계산하기 위해 필요한 $dp$ 위치를 알고 있는 상태이므로 바로 리턴.
  • 즉, $(v,m)$의 자손이 $(v,m-1)$과 $(\lfloor v/m \rfloor ,m-1)$인 이진트리를 생각하고, $m=1$이거나 $v \le n^{2/3}$인 경우를 리프 노드로 둔다.


시간복잡도를 분석하자.

  • 사실 1. 위 재귀 루틴의 시간복잡도는 $\mathcal{O}(n^{2/3})$이다. 
  • 증명. $dp$ 상태 중, $v \ge n^{2/3}$인 것의 개수를 구하면 충분하다. 실제로 관심 있는 $v$만 고려하면, 이는 $$ \sum_{i=1}^{n^{1/3}} \sqrt{n/i} = \mathcal{O}(\sqrt{n} \cdot \sqrt{n^{1/3}}) = \mathcal{O}(n^{2/3})$$
  • 사실 2. 결과적으로 우리가 "계산해야 하는 $dp$ 값"들의 위치는 $\mathcal{O}(n^{2/3})$개이다.
  • 증명. 물론, 이는 사실 1에서 바로 따라온다. 애초에 시간복잡도가 저러니 저장한 위치의 개수도 마찬가지.

이제 우리는 $\mathcal{O}(n^{2/3})$개의 $dp$ 값만 구해주면 된다. 이 값들을 다 구하면, 다시 재귀를 돌려 답을 얻을 수 있다.

즉, 이전에 "$dp[v][m]$을 계산해야 한다"고 하고 리턴한 경우, 이번에는 실제로 계산된 $dp[v][m]$의 값을 리턴하면 된다!

이제부터 우리는 알고리즘의 영역으로 온다. 처음으로 이 가이드에서 C++ STL에 없는 자료구조를 사용한다.

  • 오프라인으로 $dp$ 값들을 계산한다. $m$의 값이 증가하는 순서로 필요한 값들을 계산할 것이다.
  • $dp[v][m]$을 계산하는 것은, 결국 정의상 $m$까지 체를 돌렸을 때 $1$부터 $v$까지 살아남은 자연수의 개수를 구하는 것이다.
  • 현재 자연수가 남아있는지 여부를 0/1로 표현하면, 이는 결국 $[1, v]$에서 구간합을 구하는 것과 마찬가지다.
  • 에라토스테네스의 체에서 값을 지우는 것은 결국 해당 자연수에 대응되는 값을 1에서 0으로 바꾸는 것과 같다.
  • 그러니 이 문제는 사실 Point Update + Segment Query. 세그먼트 트리로 해결할 수 있다. (속도를 위해 펜윅을 쓰는 것을 추천한다.)

이제 전체 문제가 $\mathcal{O}(n^{2/3})$ 시간에 해결된다. 필요한 에라토스테네스의 체 크기가 $n^{2/3}$이었음에 주목하자.

  • 이 최적화 방법은 앞서 언급한 Lucy-Hedgehog 알고리즘의 여러 확장에도 바로 적용된다.
  • 예를 들어, 소수의 합을 구하고 싶다면 자연수 $k$가 남아있는지 여부를 0/$k$로 나타내면 된다.


Meissel-Lehmer Algorithm

  • 사실 Lucy-Hedgehog와 큰 차이는 없다. 개인적으로 많이 사용하지 않아, 얼마나 확장이 쉽게 되는지는 잘 모르겠다.
  • 하지만 알고리즘의 원리가 근본적으로 비슷하니, 아마 크게 다르지는 않을 것이라고 생각한다. 독자에게 맡긴다.


이번에는 notation을 살짝 바꾼다. 그래도 Lucy-Hedgehog 알고리즘과 비슷한 점이 많으니, 비교해가며 읽자.

  • $p_k$를 $k$번째 소수로 정의한다. 단, $1$-index를 이용한다. 즉, $p_1, p_2, p_3 = 2, 3, 5$.
  • $\pi(x)$를 $x$ 이하의 소수의 개수로 정의한다. 이는 매우 standard 한 notation이다. 예를 들어, $\pi(x) \approx x / \ln x$.
  • $\phi(x, a)$를 모든 소인수가 $p_{a+1}$ 이상인 $x$ 이하 자연수 개수라 하자. 중요한 점은, $p_1, p_2, \cdots ,p_a$ 역시 제거했다는 것이다.
  • 또한, $\phi(x, a)$를 생각할 때 $1$도 포함한다. 그러니, 정확하게 말하자면 $p_1, p_2, \cdots , p_a$의 배수를 제거한 결과를 보는 것이다.
  • 이어서, $P_k(x, a)$를 모든 소인수가 $p_{a+1}$ 이상인 $x$ 이하 자연수 중, 소인수 개수가 정확히 $k$개인 것의 개수라 하자.
  • 이때, "소인수 개수"란, $n=p_1^{e_1} \cdots p_k^{e_k}$라 할 때 $e_1 + e_2 + \cdots + e_k$로 정의된다. 
  • 특히, $p_{a+1}^k > x$라면 $P_k(x, a) = 0$이 성립한다. 이제 필요한 용어의 정의는 끝났다.


$\pi(x)$를 계산하려면, 필요한 값은

  • 우선 억울하게 제거당한 소수들 $p_1, p_2, \cdots , p_a$가 있다. 총 $a$개.
  • 이제 모든 소인수가 $p_{a+1}$ 이상인 $x$ 이하 자연수가 필요하다. 총 $\phi(x, a)$개.
  • 그런데 $1$은 소수가 아니니, $1$을 빼주어야 할 것이다.
  • 또한, 소인수 개수가 $2$ 이상인 것에는 관심이 없으니, $P_2(x, a) + P_3(x,a) + \cdots $를 빼야한다. 즉, $$ \pi(x) = a + \phi(x, a) - 1 - \sum_{k=2}^\infty P_k(x, a)$$
  • 물론, $p_{a+1}^k > x$면 $P_k(x, a)=0$이므로, 실제로는 무한합을 계산할 필요가 없다. 


특히, 우리는 $a = \pi(x^{1/3})$을 선택할 것이다. 그러면 $p_{a+1}^3 > x$이므로, $P_3$부터 쓸모가 없다.

그러므로, $\pi(x) = a + \phi(x, a) - 1 - P_2(x, a)$를 계산하면 된다. 이제 각 항을 계산할 전략을 설계해보자.

  • $a$의 계산. 단순히 $x^{1/3}$ 이하의 소수의 개수이므로, 에라토스테네스의 체로 정리된다.
  • $\phi(x, a)$의 계산. 이 부분은 Lucy-Hedgehog 알고리즘의 dp와 크게 다르지 않다.
  • $\phi(x, a)$를 계산해보자. 일단 $\phi(x, a-1)$에서 시작하면, $a-1$에서 $a$로 넘어갈 때 새로 제거되는 수의 개수를 세면 된다.
  • 이 값을 계산하려면, "최소 소인수가 정확히 $p_a$인 $x$ 이하의 자연수"의 개수를 세면 된다.
  • 이는 결국 "모든 소인수가 $p_{a}$ 이상인 $\lfloor x/p_a \rfloor$ 이하의 자연수"의 개수와 같다. 이는 $\phi(\lfloor x/p_a \rfloor, a-1)$.
  • 그러므로, dp 식 $\phi(x, a) = \phi(x, a-1) - \phi(\lfloor x/p_a \rfloor, a-1)$이 성립함을 알 수 있다.
  • 우선 Lucy-Hedgehog에서 본 DP와 여러 성질이 겹친다. $\lfloor x/i \rfloor$ 형태의 값만 필요하다는 점 역시 똑같다.
  • 토글링이 가능한 DP 형태라는 점도 똑같다. 그러니 이미 살펴본 접근 방식과 거의 동일한 방법이 여기서도 먹힌다.
  • 하지만 이 방법의 중요한 점은, $p_a^2 > x$여도 $\phi(x, a)$와 $\phi(x, a-1)$이 다르다는 것이다. 
  • 에라토스테네스의 체와 다르게 여기서는 $p_a$도 제거하므로, $\phi(x, a) = \phi(x, a-1) - 1$이 성립하겠다.
  • 물론, $x < p_a$면 $\phi(x, a) = \phi(x, a-1)$이 성립함은 당연하다.
  • 그래서 기존 $\mathcal{O}(n^{3/4})$ 알고리즘을 적용하기가 조금 귀찮아진다. $1$을 진짜 정직하게 빼주면 당연히 시간복잡도가 유지되지 않는다.
  • 그러니 $p_a^2 > x$면 $\phi(x, a)$를 업데이트 하지 않고, 대신 "실제로 $1$이 몇 번 더 빠졌어야 하는가"를 필요할 때 계산해주자.
  • 재밌는 점은 $\mathcal{O}(n^{2/3})$ 최적화 알고리즘은 그대로 적용된다는 것이다.
  • 언급한 사실 1, 사실 2 역시 그대로 유지되며, 오프라인 쿼리 + 세그트리/펜윅트리로 $\phi$를 계산할 수 있음도 동일하다.
  • $P_2(x, a)$의 계산. 이 값은 다행히도, 정직하게 계산할 수 있다.
  • $P_2(x, a)$는 결국 $a<b \le c$를 만족하면서 $p_bp_c \le x$인 $(b, c)$의 순서쌍의 개수를 구하는 것이다.
  • $b$를 고정하고, $c$의 개수를 세자. 그러면 $p_c \le x/p_b$와 $b \le c$가 동시에 성립해야 하니, $\pi(\lfloor x/p_b \rfloor) - (b-1)$이다.
  • 가능한 $b$의 범위는 물론 $a+1$부터 $\pi(x^{1/2})$까지다. 즉, 이 결과를 합치면 우리는 $$P_2(x, a) = \sum_{b=a+1}^{\pi(x^{1/2})} \left( \pi(\lfloor x/p_b \rfloor) - (b-1) \right)$$를 얻는다. 이제 이 값을 빠르게 계산하면 된다. 여기서 핵심은 주어진 범위의 $b$에 대하여, $x/p_b < x^{2/3}$이 성립한다는 것이다. 
  • 그러니 실제로 필요한 것은 $x^{2/3}$까지의 에라토스테네스의 체 결과다. 그러니, 모든 계산을 $\mathcal{O}(n^{2/3})$ 시간에 할 수 있다.


이를 종합하면 우리가 원하는 알고리즘이 나온다. 최적화로 $\mathcal{O}(n^{2/3})$까지 가는 방법도 설명하였다. 그런데,

  • 소수 카운팅 알고리즘의 속도를 측정할 수 있는 가장 쉬운 방법은, Codeforces의 Four Divisors 문제를 해결하는 것이다.
  • 당장 Editorial에서도 이 글에서 다룬 $\mathcal{O}(n^{2/3})$ 테크닉을 언급한다. 자세한 설명이 있으니 읽어보는 것도 좋다.
  • 그런데, 정작 submission을 실행 속도로 정렬한 다음 몇 개의 제출을 읽어보면 오프라인이고 세그트리고 다 없다.
  • 즉, 시간복잡도로 표현하는 속도와 실제 실행 속도와의 괴리감이 약간 있는 부분이라고 볼 수 있겠다. 이제 실제로 빠른 알고리즘을 간략하게 설명한다.
  • (물론, 효율적인 코드를 매우 잘 작성하는 분들이라면 $\mathcal{O}(n^{2/3})$ 알고리즘으로 저 분들을 이길 수도 있을 것이다.)
  • 먼저, 에라토스테네스의 체를 많이 돌린다. $n=10^{11}$ 기준으로, 500만 정도 크기로 돌린 것 같다.
  • 앞서 우리는 $\pi(x^{1/3})$에서 식을 자르는 방식으로 $P_3(x, a)$ 항을 날렸다. 이번에는 $\pi(x^{1/4})$에서 식을 잘라서 $P_3(x, a)$ 항을 살린다.
  • 즉, $a = \pi(x^{1/4})$를 잡고, $\pi(x) = a + \phi(x, a) - 1 - P_2(x, a) - P_3(x, a)$를 계산한다.
  • 이제 추가적으로 $P_3$에 대한 전개식을 계산해야 한다. $p_bp_cp_d$ 형태를 보면 되는데, 진짜 정직하게 $b, c$에 대한 이중 for 문을 돌린다.
  • 이제 $\phi(x, a)$를 계산해야 한다. 이 값을 계산하기 위해서, 상당한 최적화를 한다.
  • $\phi(x, a)$의 context에서, 특정 값이 살아있는지 여부는 그 값을 $p_1p_2 \cdots p_a$로 나눈 나머지에 의해 결정된다.
  • 그러니, $a$가 작은 경우 ($a \le 7$을 이용한다) 아예 어떤 나머지가 살아남는지를 통째로 전처리를 한다. 크기 $p_1 \cdots p_a$짜리 전처리다.
  • 이러면 $a \le 7$인 입력이 들어온 경우, 아예 $\mathcal{O}(1)$에 $\phi(x, a)$를 반환할 수 있다. 이것이 첫 최적화.
  • 만약 $x < p_a^2$이라면, $\phi(x, a)$는 사실 $p_{a+1}$ 이상 $x$ 이하 소수의 개수이다. 
  • 그러니, 만약 $x$ 이하의 소수 개수가 체로 전처리된 상황이라면, 바로 답을 반환할 수 있다. 이것이 두 번째 최적화.
  • 만약 $x < p_a^3$이라면, $\phi(x, a) = \pi(x) - a + 1 + P_2(x, a)$다. 이때, $\pi(x)$의 값이 이미 체로 전처리된 상황이라 하자.
  • 이 경우에는, $P_2(x, a)$를 기존 방식과 동일한 방법으로 구해 $\phi(x, a)$를 얻는다. 이것이 세 번째 최적화.
  • 만약 이 모든 경우에 해당하지 않는다면, 그제서야 점화식 $\phi(x, a) = \phi(x, a-1) - \phi(\lfloor x/p_a \rfloor, a-1)$을 호출한다.
  • 사실 기존 Meissel-Lehmer 알고리즘에서 $a = \pi (x^{1/2})$를 선택하는 것도 방법이다. 방법은 정말 많다..
실제로 어떻게 구현되었는지 보려면, 이 링크를 타고 여러 코드를 탐색해보자. 상당히 빠르다 :) 

$\mathcal{O}(n^{2/3})$ 알고리즘들이 메모리를 많이 잡아먹는데 (물론, 적게 먹는 variant가 있는 것으로 안다) 주의해야 한다.

필자는 가능하면 Lucy-Hedgehog를 쓰고 (특히 Project Euler와 같은 시간제한이 없는 경우에는 더더욱) 만약 시간제한이 빡빡하다면 Meissel-Lehmer의 최적화된 구현체, 즉 위 Four Divisors 문제의 링크에 있는 알고리즘을 사용하여 문제를 푼다. 모두 흥미로운 알고리즘들이다. 이 정도로 이번 내용은 끝.