淺談旅行推銷員問題 (Travelling Salesman Problem) 搜索加速

學弟們修課〈演算法設計與分析〉內容中的一環,看他們的講義進行加速比較,也許有人會認為我太閒,但這只是追求真理。看看維基百科上,就有那上萬點的極佳近似解算法,而接下來要說的相當於玩具等級。

Task Description

給一組城市和每一個城市之間的單向距離,請找出從起點城市出發,經過所有城市且不重複回到起點城市的最短距離和。

例如,走訪順序為 $1-2-4-3-1$. 距離總和為 $10+25+30+15 = 80$.

Subtask

$N \le 100$

Input Format

只有一組測資,每組第一行有一個整數 $N$ ,表示有 $N$ 座城市,接下來會有 $N$ 行,每行上有 $N$ 個非負整數 $D_{i, j}$ 表示節點 $i$ 到節點 $j$ 的距離,如果 $D_{i, j} = 0$ ,視同沒有邊。

你可以假設每一條邊可能是單向的,保證至少一條 TSP 路徑。

  • $3 < N \leq 100$
  • $0 < \textit{distance between any pair of cities} \le 10000$

Output Format

輸出一行最小距離和。

Sample Input 1

1
2
3
4
5
4
0 10 15 20
10 0 35 25
15 35 0 30
20 25 30 0

Sample Output 1

1
80

Solution

初學第一步-窮舉

從最直觀的想法中,我們可以使用窮舉 $O(N!)$ 解決,再搭配簡單的剪枝 (加上當前點 $u$ 返回起點的最短路徑長是否大於已知解)。當 $N \ge 12$ 時,理論運算次數達 4 億之多,運算時間相當於 1 秒左右 (粗估一般 CPU 可達 4 GHz),接下來窮舉速度就不堪負荷。

時間複雜度 $O(N!)$ , 空間複雜度 $O(N)$

踏出第一步-動態規劃

當學會動態規劃後,解決這一 NP-Complete 問題時,我們使用記憶體空間去掉重複換取時間,定義狀態 dp[i][j] 表示當前走過的點集合 $i$ ,停留最後一個點位置 $j$ ,點集合可以用 bitmask 的技巧壓縮,集合狀態個數 $2^N$ ,窮舉所有停留點轉移需要 $O(N)$ 次。當 $N \ge 25$ 時,理論次數達 8 億之多,再加上記憶體存取效能,導致運行約數十秒。

時間複雜度 $O(N^2 \cdot 2^N)$ , 空間複雜度 $O(N \cdot 2^N)$

更進一步

我們發現到基礎窮舉中可以剪枝卻又不夠局觀,在動態規劃中擁有移除重複卻又不能剪枝,動態規劃又佔據大量的記憶體空間,那麼使用更厲害的 branch-and-reduce 算法吧。簡單介紹一下 branch-and-reduce 算法的特性:

  • 屬於窮舉的一環
  • 在每一步中,花費較多的時間估算整體基礎花費,用來剪枝
  • 窮舉一步之後,將盤面選擇縮小,甚至做到同構移除來減少窮舉分支,這個步驟稱為 reduce
  • 使用的記憶體空間少 (與一般窮舉相當)

運用在 TSP 問題時,我們首先將輸入轉換成 $N \times N$ 矩陣,當移動從 $u$$v$ 時,矩陣就會變成 $(N-1) \times (N-1)$ 大小。

計算基礎花費

對於一個 TSP 而言,每一個點必然有入邊和出邊,因此基礎花費必然每一個節點權重最小的入邊和出邊總和。範例操作如下:

$$\begin{align*} \text{cost matrix } M_5 = \begin{bmatrix} \infty & 20 & 30 & 10 & 11 \\ 15 & \infty & 16 & 4 & 2 \\ 3 & 5 & \infty & 2 & 4 \\ 19 & 6 & 18 & \infty & 3 \\ 16 & 4 & 7 & 16 & \infty \end{bmatrix} \end{align*}$$

計算出邊最小值,將每一行的最小值算出,例如第一行 $[\infty \; 20 \; 30 \; 10 \; 11]$ 為 10,最後分別扣掉 10, 2, 2, 3, 4,得到出邊基礎花費 $10+2+2+3+4 = 21$ ,接著我們扣除掉基礎花費得到

$$\begin{align*} \text{reduce row, cost matrix } M_5 = \begin{bmatrix} \infty & 10 & 20 & 0 & 1 \\ 13 & \infty & 14 & 2 & 0 \\ 1 & 3 & \infty & 0 & 2 \\ 16 & 3 & 15 & \infty & 0 \\ 12 & 0 & 3 & 12 & \infty \end{bmatrix} \end{align*}$$

此時,我們再針對入邊找到基礎花費,將每一個列的最小值算出,分別為 1, 0, 3, 0, 0,因此得到基礎花費 $21+1+3 = 25$ ,如此一來當下限高於當前的最佳解就可以退出。下一階段的矩陣就會被改造成如下所示

$$\begin{align*} \text{reduce row/column, cost matrix }M_5 = \begin{bmatrix} \infty & 10 & 17 & 0 & 1 \\ 12 & \infty & 11 & 2 & 0 \\ 0 & 3 & \infty & 0 & 2 \\ 15 & 3 & 12 & \infty & 0 \\ 11 & 0 & 0 & 12 & \infty \end{bmatrix} \end{align*}$$

當我們從某一點 $u$ 移動到 $v$ 時,移除掉矩陣中的某一行 $u$ 和某一列 $v$ 即可。特別注意到,上述範例為完全圖,對於稀疏圖會發生無解情況,需要特別判斷。

優化一夜

用陣列實作不算困難,當抽出一行一列時,可能會需要複製整個矩陣,又或者只用標記來防止計算到經過的點,不管哪一種寫法都不算漂亮。這裡提供一個解決方法:在處理基礎花費時,對整個矩陣每行每列扣一個定值,實作時不真正修改數值,而是把差值記錄下來,用 $O(N)$ 空間 $Dx[], \; Dy[]$ 記錄差值,如此一來就不用耗費 $O(N^2)$ 空間和時間進行複製。

相對地,這會造成計算真正值 $M'_{x, y} = M_{x, y} - Dx[x] - Dy[y]$ 需要數次的記憶體存取。無妨地,我們已經將了每一層的空間和搬運次數,又因為陣列小到可以全塞進快取中。

優化二夜

對於稀疏圖而言,用矩陣實作過於浪費空間,而且很容易拖累速度,因此可以使用雙十字鏈表 Dancing Links 來模擬上述的矩陣,高效地降低運行常數。實作複雜度指數上升。

優化三夜

遞迴窮舉必然會遇到復原操作,除非我們完整地複製每一層狀態。通常我們只會修改差異,這裡需要對整個矩陣修改,這導致差異數量非常多。為了加速復原操作,遞迴時維護自己 stack 來復原修改差異,而且僅當 $Dx[], Dy[]$ 不為零的時候才進行儲存,因為一部分的 reduce 結果會是零,這些影響的效能非常多。

除了上述外,計算基礎花費仍佔據了 $O(N^2)$ 的時間並佔大部分的運行時間。為了減少計算量,可以做到以下幾點:

  • 在每一行/列計算最小值時,走訪 Dancing Links 每一行每一列,如果當前最小值已經被更新到 0 時,可以直接退出。
  • 傳遞當前已知最佳解,當累積計算花費總和高於已知解就退出。

優化四夜

從啟發的觀點來看,每一階段有 $n$ 種選擇,各別算出預期最少花費後,排序後優先從花費小的開始拓展。

下述程式只支援頂點個數 30 以內,在頂點個數 30 時,只需要耗費 161 ms 即可完成。
當我將頂點個數放大 50 時,小夥伴們的程序花了 2 個小時還沒結果,而我寫的版本跑了 1X 秒就結束。測資的極限就只到這了,暫時先擱著吧。題目鏈結 20007. Fast Travelling Salesman Problem

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
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
#include <bits/stdc++.h>
using namespace std;

#define MAXNODE 5000
#define MAXCOL 32
#define MAXN 32
#define INF 0x3f3f3f
struct DancingNode {
int left, right, up, down;
int ch, rh, w;
} DL[MAXNODE];
struct HelperNode {
int head, size, next, prev;
} HN[MAXNODE];
int help_head;
int s[MAXCOL];
int head, size, Ans;
void Remove(int c) {
DL[DL[c].right].left = DL[c].left;
DL[DL[c].left].right = DL[c].right;
for (int i = DL[c].down; i != c; i = DL[i].down) {
if (HN[DL[i].rh].head == i)
HN[DL[i].rh].head = DL[i].left;
HN[DL[i].rh].size--;
DL[DL[i].right].left = DL[i].left;
DL[DL[i].left].right = DL[i].right;
s[DL[i].ch]--;
}
}
void Resume(int c) {
for (int i = DL[c].down; i != c; i = DL[i].down) {
if (HN[DL[i].rh].head == i)
HN[DL[i].rh].head = DL[i].right;
HN[DL[i].rh].size++;
DL[DL[i].right].left = i;
DL[DL[i].left].right = i;
s[DL[i].ch]++;
}
DL[DL[c].right].left = c;
DL[DL[c].left].right = c;
}
void Reduce(int i) {
int t = DL[i].rh;
HN[HN[t].next].prev = HN[t].prev;
HN[HN[t].prev].next = HN[t].next;
s[DL[i].ch]--;

DL[DL[i].down].up = DL[i].up;
DL[DL[i].up].down = DL[i].down;
for (int k = DL[i].right; k != i; k = DL[k].right) {
DL[DL[k].down].up = DL[k].up;
DL[DL[k].up].down = DL[k].down;
s[DL[k].ch]--;
}
}
void Recover(int i) {
int t = DL[i].rh;
HN[HN[t].next].prev = t;
HN[HN[t].prev].next = t;
s[DL[i].ch]++;

DL[DL[i].down].up = i;
DL[DL[i].up].down = i;
for (int k = DL[i].right; k != i; k = DL[k].right) {
DL[DL[k].down].up = k;
DL[DL[k].up].down = k;
s[DL[k].ch]++;
}
}

void Select(int i, int j) {
int s = DL[i].rh;
Reduce(i);
Remove(j);
HN[HN[s].next].prev = HN[s].prev;
HN[HN[s].prev].next = HN[s].next;
}

void Cancel(int i, int j) {
int s = DL[i].rh;
Resume(j);
Recover(i);
HN[HN[s].next].prev = s;
HN[HN[s].prev].next = s;
}

int new_node(int up, int down, int left, int right) {
DL[size].up = up, DL[size].down = down;
DL[size].left = left, DL[size].right = right;
DL[up].down = DL[down].up = DL[left].right = DL[right].left = size;
return size++;
}

void new_row(int n, int R[][2], int rh) {
int r, row = -1, k;
int h = size;
for (int i = 0; i < n; i++) {
r = R[i][0];
DL[size].ch = r, s[r]++;
DL[size].rh = rh;
DL[size].w = R[i][1];
if (row == -1) {
row = new_node(DL[DL[r].ch].up, DL[r].ch, size, size);
} else {
k = new_node(DL[DL[r].ch].up, DL[r].ch, DL[row].left, row);
}
}
HN[rh].size = n;
HN[rh].head = h;
HN[rh].next = HN[help_head].next;
HN[rh].prev = help_head;
HN[HN[help_head].next].prev = rh;
HN[help_head].next = rh;
}
void init(int m) {
size = 0;
help_head = 0, HN[help_head].next = HN[help_head].prev = help_head;
head = new_node(0, 0, 0, 0);
for (int i = 1; i <= m; i++) {
new_node(i, i, DL[head].left, head);
DL[i].ch = i, s[i] = 0;
DL[i].rh = 0; // important pointer (fail pointer)
}
}

int Dx[MAXN], Dy[MAXN];
int markRStk[MAXNODE], markRidx = -1;
int markCStk[MAXNODE], markCidx = -1;

void printDL(int n) {
for (int i = HN[help_head].prev; i != help_head; i = HN[i].prev) {
int j = HN[i].head, prev = 1;
printf("%2d ", DL[j].rh);
do {
while (prev < DL[j].ch)
prev++, printf(" -- ");
prev = DL[j].ch+1;
int rh = DL[j].rh, ch = DL[j].ch;
printf("[%2d]", DL[j].w - Dx[rh] - Dy[ch]);
j = DL[j].right;
} while (j != HN[i].head);
while (prev <= n)
prev++, printf(" -- ");
puts("");
}
puts("----------------------");
}
int cH(int limit = INF) {
int ins = 0;
for (int i = HN[help_head].prev; i != help_head; i = HN[i].prev) {
if (HN[i].size == 0)
return INF;
int j = HN[i].head, v = INF;
int rh = DL[j].rh, ch;
do {
ch = DL[j].ch;
v = min(v, DL[j].w - Dx[rh] - Dy[ch]);
j = DL[j].right;
} while (j != HN[i].head && v);
if (v == INF || ins+v >= limit)
return INF;
if (v) {
ins += v;
Dx[rh] += v;
markRStk[++markRidx] = v;
markRStk[++markRidx] = rh;
}
}
/*
printf("cH first\n");
printDL(4);
*/

for (int i = DL[head].right; i != head; i = DL[i].right) {
if (DL[i].down == i)
return INF;
int j = DL[i].down, v = INF;
int ch = DL[i].ch, rh;
do {
rh = DL[j].rh;
v = min(v, DL[j].w - Dx[rh] - Dy[ch]);
j = DL[j].down;
} while (j != i && v);
if (v == INF || ins+v >= limit)
return INF;
if (v) {
ins += v;
Dy[ch] += v;
markCStk[++markCidx] = v;
markCStk[++markCidx] = ch;
}
}
return ins;
}
void rH(int backRidx, int backCidx) {
while (markRidx > backRidx) {
int a = markRStk[markRidx--];
int b = markRStk[markRidx--];
Dx[a] -= b;
}
while (markCidx > backCidx) {
int a = markCStk[markCidx--];
int b = markCStk[markCidx--];
Dy[a] -= b;
}
}

#ifdef SAME_CUT
int isSame(int vx, int vy) {
int i = HN[vx].head, j = HN[vy].head;
do {
if (DL[i].ch != DL[j].ch && (DL[i].ch != vy && DL[j].ch != vx))
return 0;
int e1 = DL[i].w - Dx[vx] - Dy[DL[i].ch];
int e2 = DL[j].w - Dx[vy] - Dy[DL[j].ch];
if (e1 != e2)
return 0;
i = DL[i].right;
j = DL[j].right;
} while (i != HN[vx].head && j != HN[vy].head);
return i == HN[vx].head && j == HN[vy].head;
}
#endif

#ifdef MST_CUT
int parent[MAXN];
int findp(int x) {
return x == parent[x] ? x : (parent[x]=findp(parent[x]));
}
int joint(int x, int y) {
x = findp(x), y = findp(y);
if (x == y) return 0;
parent[y] = x;
return 1;
}
int MST(int st, int ed, int limit = INF) {
static int used[MAXN] = {}, cases = 0;
if (HN[HN[help_head].prev].prev == help_head)
return 0;
int i = HN[help_head].prev;
int n = 0;
int ret = 0, in = INF, out = INF;
cases++;
do {
int j = HN[i].head;
do {
int x = DL[j].rh, y = DL[j].ch;
int w = DL[j].w - Dx[x] - Dy[y];
if (x == st) {
in = min(in, w);
} else if (y == ed) {
out = min(out, w);
} else {
if (used[x] != cases)
used[x] = cases, parent[x] = x, n++;
if (used[y] != cases)
used[y] = cases, parent[y] = y, n++;
n -= joint(x, y);
}
j = DL[j].right;
} while (j != HN[i].head);
i = HN[i].prev;
} while (i != help_head);
ret = in + out;
if (ret >= limit)
return limit;
if (n > 1) {
fprintf(stderr, "cc cut\n");
return INF;
}
return ret;
}
#endif
void DFS(int u, int lower_bound) {
#ifdef MST_CUT
if (lower_bound + MST(u, 1, Ans - lower_bound) >= Ans) {
return ;
}
#else
if (lower_bound >= Ans) {
return ;
}
#endif

if (HN[HN[help_head].prev].prev == help_head) {
Ans = lower_bound;
return ;
}

int i = HN[u].head;
vector< pair<int, int> > pick;
do {
int v = DL[i].ch;
int runnable = 1;
if (v == 1)
runnable = 0;
#ifdef SAME_CUT
if (runnable) {
int e1 = DL[i].w - Dx[u] - Dy[v];
for (int j = HN[u].head; j != i && runnable; j = DL[j].right) {
int e2 = DL[j].w - Dx[u] - Dy[DL[j].ch];
if (e1 != e2)
continue;
if (isSame(v, DL[j].ch))
runnable = 0;
}
}
#endif

if (runnable) {
int backRidx = markRidx, backCidx = markCidx;
int tu = HN[u].head;
Select(tu, v);
int c = cH(Ans - lower_bound);
if (c != INF) {
pick.push_back(make_pair(lower_bound + c + DL[i].w-Dx[u]-Dy[v], v));
// DFS(v, lower_bound + c + DL[i].w-Dx[u]-Dy[v]);
}
rH(backRidx, backCidx);
Cancel(tu, v);
}
i = DL[i].right;
} while (i != HN[u].head);

sort(pick.begin(), pick.end());

for (int i = 0; i < pick.size(); i++) {
int v = pick[i].second;
/*
printf("%d -> %d, %d\n", u, v, pick[i].first);
printDL(4);
*/

int backRidx = markRidx, backCidx = markCidx;
int tu = HN[u].head;
Select(tu, v);
/*
printf("Select\n");
printDL(4);
*/

int c = cH(Ans - lower_bound);
/*
printf("Reduce\n");
printDL(4);
*/

if (c != INF)
DFS(v, pick[i].first);
rH(backRidx, backCidx);
Cancel(tu, v);
}
}

int main() {
int G[MAXN][MAXN];
int n;
while (scanf("%d", &n) == 1) {
init(n);
for (int i = 1; i <= n; i++) {
int R[MAXN][2], Rcnt = 0;
for (int j = 1; j <= n; j++) {
assert(scanf("%d", &G[i][j]) == 1);
if (G[i][j])
R[Rcnt][0] = j, R[Rcnt][1] = G[i][j], Rcnt++;
}
new_row(Rcnt, R, i);
}

int backRidx = markRidx, backCidx = markCidx;
Ans = INF;
DFS(1, cH());
rH(backRidx, backCidx);
printf("%d\n", Ans);

}
return 0;
}
Read More +

淺談多重背包問題 (Multiple Knapsack Problem) 優化那些事

收錄於 批改娘 20008. Fast Multiple Knapsack Problem

之所以出了這一題,源自於實驗室另一名同學跑實驗太久,進而撰寫優化程序。聽說原本跑了十分鐘的實驗,改善後提升了到一分鐘跑完。

輸入格式

每組測資第一行包含兩個正整數,分別代表背包大小 $M$ ($\leq 10^6$) 和物品個數 $N$ ($\leq 1000$),下一行開始每行包含兩個正整數,分別代表物品價值 $P_i$ ($\leq 10^3$)、物品重量 $W_i$ ($ \leq 10^5$) 以及物品最多可以挑 $C_i$ 個 ($\le 100$)。

輸出格式

對於每組測資,請輸出最大收益。

範例輸入 1

1
2
3
4
5
6
7
8
50 7
66 31 1
232 10 4
49 20 1
54 19 1
426 4 3
589 3 10
10 6 4

範例輸出 1

1
7178

Solution

不管是 0/1 背包或者多重背包,兩者都屬於 bounded knapsack problem 問題。即便如此,優化上仍有些許的不同,請讓我緩緩道來。

在此之前,您必須先理解上一篇《淺談背包問題 (0/1 Knapsack Problem) 優化那些事》的部分,不然會造成閱讀上的困難。

多重背包有一個二進制優化,也就是當物品限制最多拿 $C$ 個時,我們可以利用二進制組合的方式,轉換到 0/1 背包問題,因此我們會得到新的 $N \log C$ 個物品跑一次 0/1 背包,因此複雜度落在 $O(N \log C \; W)$

然而,從公式定義上,在好幾年前的論文中,使用斜率優化降到 $O(N \; W)$,推倒過程如下,

$j = k \; w_i + r, \; 0 \le r \le w_i - 1$ $$\begin{align*} dp[i][j] &= \max\left\{dp[i-1][j], dp[i-1][j-w_i]+p_i, \cdots, dp[i-1][j-c_i \; w_i] + c_i \; p_i\right\} \\ &= \max\left\{dp[i-1][k \; w_i + r], dp[i-1][(k-1) \; w_i + r] + p_i, \cdots , dp[i-1][(k-c_i) \; w_i + r] + c_i p_i\right\} \\ &= \max\left\{dp[i-1][k \; w_i + r] - k \; p_i, dp[i-1][(k-1) \; w_i + r] + (k-1) \; p_i, \cdots , dp[i-1][(k-c_i) \; w_i + r] - (k-c_i) p_i\right\} + k \; p_i\\ \end{align*}$$

隨著式子的轉移,我們發現每一個取值將不依賴相對位置,只跟自身的位置有關,那麼可以使用單調堆 (monotone queue) 運行 $O(1)$ 的 sliding windows 查找極值。最後,將相同餘數分堆處理,單調堆中最多存在 $O(c)$ 個元素。

優化初夜

如果只使用二進制優化,套上我們的 0/1 優化方案,將有大幅度地提升。

加入 0/1 背包的優化策略,再套上最簡單的斜率優化算法,得到下面的程式。這裡很懶惰地,由於單調堆最多入隊 $W$ 次,不外乎地直接只用大小為 $W$ 的方式實作。

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
#include <bits/stdc++.h>
using namespace std;

namespace {

static const int MAXW = 1000005;
static const int MAXN = 1005;

struct BB {
int w, v, c;
BB(int w = 0, int v = 0, int c = 0):
w(w), v(v), c(c) {}
bool operator<(const BB &x) const {
return w * c < x.w * x.c;
}
};
static int run(BB A[], int dp[], int W, int N) {
static int MQ[MAXW][2];
for (int i = 0, sum = 0; i < N; i++) {
int w = A[i].w, v = A[i].v, c = A[i].c;
sum = min(sum + w*c, W);
for (int j = 0; j < w; j++) {
int l = 0, r = 0;
MQ[l][0] = 0, MQ[l][1] = dp[j];
for (int k = 1; k*w+j <= sum; k++) {
if (k - MQ[l][0] > c)
l++;
int dpv = dp[k*w+j] - k*v;
while (l <= r && MQ[r][1] <= dpv)
r--;
r++;
MQ[r][0] = k, MQ[r][1] = dpv;
dp[k*w+j] = max(dp[k*w+j], MQ[l][1] + k*v);
}
}
}
}

static int knapsack(int C[][3], int N, int W) {
vector<BB> A;
for (int i = 0; i < N; i++) {
int w = C[i][0], v = C[i][1], c = C[i][2];
A.push_back(BB(w, v, c));
}
assert(N < MAXN);
static int dp1[MAXW+1], dp2[MAXW+1];
BB Ar[2][MAXN];
int ArN[2] = {};
memset(dp1, 0, sizeof(dp1[0])*(W+1));
memset(dp2, 0, sizeof(dp2[0])*(W+1));
sort(A.begin(), A.end());
int sum[2] = {};
for (int i = 0; i < N; i++) {
int ch = sum[1] < sum[0];
Ar[ch][ArN[ch]] = A[i];
ArN[ch]++;
sum[ch] = min(sum[ch] + A[i].w*A[i].c, W);
}

run(Ar[0], dp1, W, ArN[0]);
run(Ar[1], dp2, W, ArN[1]);

int ret = 0;
for (int i = 0, j = W, mx = 0; i <= W; i++, j--) {
mx = max(mx, dp2[i]);
ret = max(ret, dp1[j] + mx);
}
return ret;
}
}

int main() {
int W, N;
assert(scanf("%d %d", &W, &N) == 2);
int C[MAXN][3];
for (int i = 0; i < N; i++)
assert(scanf("%d %d %d", &C[i][1], &C[i][0], &C[i][2]) == 3);
printf("%d\n", knapsack(C, N, W));
return 0;
}

不幸地,相較於一般的斜率優化寫法,並沒有太大的改善。

優化二夜

運行 sliding windows 操作時,前 $c$ 次,是不會進行 pop_front() 操作的,因此把迴圈分兩堆處理,增加 branch predict。以及在乘數運算上,使用強度減少 (strength reduction) 的技術,將乘法換成加法。

只能些許地改善 5% 的效能。

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
static int run(BB A[], int dp[], int W, int N) {
static int MQ[MAXW][2];
for (int i = 0, sum = 0; i < N; i++) {
int w = A[i].w, v = A[i].v, c = A[i].c;
sum = min(sum + w*c, W);
for (int j = 0; j < w; j++) {
int l = 0, r = 0;
MQ[l][0] = 0, MQ[l][1] = dp[j];
for (int k = 1, tw = w+j, tv = v; tw <= sum && k <= c; k++, tw += w, tv += v) {
int dpv = dp[tw] - tv;
while (l <= r && MQ[r][1] <= dpv)
r--;
r++;
MQ[r][0] = k, MQ[r][1] = dpv;
dp[tw] = max(dp[tw], MQ[l][1] + tv);
}
for (int k = c+1, tw = (c+1)*w+j, tv = (c+1)*v; tw <= sum; k++, tw += w, tv += v) {
if (k - MQ[l][0] > c)
l++;
int dpv = dp[tw] - tv;
while (l <= r && MQ[r][1] <= dpv)
r--;
r++;
MQ[r][0] = k, MQ[r][1] = dpv;
dp[tw] = max(dp[tw], MQ[l][1] + tv);
}
}
}
}

優化三夜

後來發現,sliding windows 滑動時,我們常常看前看後,因此常常會發生 cache miss,因為他要跳躍一大段記憶體空間查找數值,所以可以考慮花點操作將極值放在 stack 上,視為一種 software cache 來加速,來減少 cache miss 的懲罰。

接著,在迴圈邊界比較時,我們可以算得更精準些,回到一般的 i++ 的 format pattern,讓編譯器幫我們做常見的迴圈優化。

改善了 10% 效能

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
 static int run(BB A[], int dp[], int W, int N) {
static int MQ[MAXW][2];
for (int i = 0, sum = 0; i < N; i++) {
int w = A[i].w, v = A[i].v, c = A[i].c;
sum = min(sum + w*c, W);
for (int j = 0; j < w; j++) {
int l = 0, r = 0;
MQ[l][0] = 0, MQ[l][1] = dp[j];
int cache_max = MQ[l][1], cache_idx = MQ[l][0];
int k_bound;
k_bound = min((sum-j)/w, c);
for (int k = 1, tw = w+j, tv = v; k <= k_bound; k++, tw += w, tv += v) {
// tw = k*w+j, tv = k*v;
int dpv = dp[tw] - tv;
while (l <= r && MQ[r][1] <= dpv)
r--;
r++;
MQ[r][0] = k, MQ[r][1] = dpv;
if (r == l) cache_max = dpv, cache_idx = k;
dp[tw] = max(dp[tw], cache_max + tv);
}
k_bound = (sum-j)/w;
for (int k = c+1, tw = (c+1)*w+j, tv = (c+1)*v; k <= k_bound; k++, tw += w, tv += v) {
int dpv = dp[tw] - tv;
while (l <= r && MQ[r][1] <= dpv)
r--;
r++;
MQ[r][0] = k, MQ[r][1] = dpv;
if (r == l)
cache_max = dpv, cache_idx = k;
else if (k - cache_idx > c)
l++, cache_idx = MQ[l][0], cache_max = MQ[l][1];
dp[tw] = max(dp[tw], cache_max + tv);
}
}
}
}

優化四夜

儘管上面使用的 software cache 的方式減少 cache miss,但 DP table 仍與數據結構的記憶體位置相當遙遠,為了使他們貼近,應使用環狀隊列的實作,空間從 $O(W)$ 將到 $O(N)$,實作時,將大小限制在 $2^k$,方便運行時使用 AND 運算取代耗時的模數運算。

由於限制個數分佈上,很容易造成貪心算法有解,因此先跑一次貪心,如果貪心沒辦法達到剛好大小,那麼再跑 DP 找解。DP 找解時,可以將物品嘗試進行二進制轉換,將等價物品合併,來觸發計算邊界的優化。完成的程序如下:

最終加速,改善 10%,期待你我的分享增進。根據鴿籠原理,cp 直種類不多時,可以高達 10 倍以上的加速。

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
#include <bits/stdc++.h>
using namespace std;

namespace {

static const int MAXW = 1000005;
static const int MAXN = 1005;
static const int MAXC = 1<<12;

struct BB {
int w, v, c;
BB(int w = 0, int v = 0, int c = 0):
w(w), v(v), c(c) {}
bool operator<(const BB &x) const {
return w * c < x.w * x.c;
}
};
static bool cmpByWeight(BB x, BB y) {
return x.w < y.w;
}
static int run(BB A[], int dp[], int W, int N) {
static int MQ[MAXC][2];

for (int i = 0, sum = 0; i < N; i++) {
int w = A[i].w, v = A[i].v, c = A[i].c;
assert(c < MAXC);
sum = min(sum + w*c, W);
if (c != 1) {
for (int j = 0; j < w; j++) {
int l = 0, r = 0;
MQ[r][0] = 0, MQ[r][1] = dp[j];
int cache_max = MQ[r][1], cache_idx = MQ[r][0];
int k_bound;
r = (r+1)&(MAXC-1);
k_bound = min((sum-j)/w, c);
for (int k = 1, tw = w+j, tv = v; k <= k_bound; k++, tw += w, tv += v) {
// tw = k*w+j, tv = k*v;
int dpv = dp[tw] - tv;
while (l != r && MQ[(r-1+MAXC)&(MAXC-1)][1] <= dpv)
r = (r-1+MAXC)&(MAXC-1);
MQ[r][0] = k, MQ[r][1] = dpv;
if (l == r) cache_max = dpv, cache_idx = k;
r = (r+1)&(MAXC-1);
dp[tw] = max(dp[tw], cache_max + tv);
}
k_bound = (sum-j)/w;
for (int k = c+1, tw = (c+1)*w+j, tv = (c+1)*v; k <= k_bound; k++, tw += w, tv += v) {
int dpv = dp[tw] - tv;
while (l != r && MQ[(r-1+MAXC)&(MAXC-1)][1] <= dpv)
r--;
MQ[r][0] = k, MQ[r][1] = dpv;
if (l == r) cache_max = dpv, cache_idx = k;
else if (k - cache_idx > c)
l = (l+1)&(MAXC-1), cache_idx = MQ[l][0], cache_max = MQ[l][1];
r = (r+1)&(MAXC-1);
dp[tw] = max(dp[tw], cache_max + tv);
}
}
} else if (c == 1) {
for (int j = sum; j >= w; j--)
dp[j] = max(dp[j], dp[j-w]+v);
}
}
}
static int greedy(int C[][3], int N, int W) {
struct GB {
int w, v, c;
GB(int w = 0, int v = 0, int c = 0):
w(w), v(v), c(c) {}
bool operator<(const GB &x) const {
if (v * x.w != x.v * w)
return v * x.w > x.v * w;
return c > x.c;
}
};
vector<GB> A;
for (int i = 0; i < N; i++) {
int w = C[i][0], v = C[i][1], c = C[i][2];
A.push_back(GB(w, v, c));
}
sort(A.begin(), A.end());
int ret = 0;
for (int i = 0; i < N; i++) {
int t = min(A[i].c, W/A[i].w);
if (t == 0)
return -1;
W -= t*A[i].w;
ret += t*A[i].v;
if (W == 0)
return ret;
}
return ret;
}
static int knapsack(int C[][3], int N, int W) {
// filter
{
int filter = greedy(C, N, W);
if (filter != -1)
return filter;
}
vector<BB> A;
for (int i = 0; i < N; i++) {
int w = C[i][0], v = C[i][1], c = C[i][2];
A.push_back(BB(w, v, c));
}
// reduce
{
sort(A.begin(), A.end(), cmpByWeight);
map<pair<int, int>, int> R;
for (int i = 0; i < N; i++)
R[make_pair(A[i].w, A[i].v)] = i;
for (int i = 0; i < N; i++) {
int c = A[i].c;
map<pair<int, int>, int>::iterator it;
for (int k = 1; k <= c; k <<= 1) {
int w = A[i].w * k, v = A[i].v * k;
it = R.find(make_pair(w, v));
if (it != R.end() && i != it->second) {
int j = it->second;
A[j].c ++;
A[i].c -= k;
}
c -= k;
}
if (c > 0) {
int w = A[i].w * c, v = A[i].v * c;
it = R.find(make_pair(w, v));
if (it != R.end() && i != it->second) {
int j = it->second;
A[j].c ++;
A[i].c -= c;
}
}
}
}
static int dp1[MAXW+1], dp2[MAXW+1];
BB Ar[2][MAXN];
int ArN[2] = {};
memset(dp1, 0, sizeof(dp1[0])*(W+1));
memset(dp2, 0, sizeof(dp2[0])*(W+1));
sort(A.begin(), A.end());
int sum[2] = {};
for (int i = 0; i < N; i++) {
int ch = sum[1] < sum[0];
Ar[ch][ArN[ch]] = A[i];
ArN[ch]++;
sum[ch] = min(sum[ch] + A[i].w*A[i].c, W);
}

run(Ar[0], dp1, W, ArN[0]);
run(Ar[1], dp2, W, ArN[1]);

int ret = 0;
for (int i = 0, j = W, mx = 0; i <= W; i++, j--) {
mx = max(mx, dp2[i]);
ret = max(ret, dp1[j] + mx);
}
return ret;
}
}

int main() {
int W, N;
assert(scanf("%d %d", &W, &N) == 2);
int C[MAXN][3];
for (int i = 0; i < N; i++)
assert(scanf("%d %d %d", &C[i][1], &C[i][0], &C[i][2]) == 3);
printf("%d\n", knapsack(C, N, W));
return 0;
}
Read More +

淺談背包問題 (0/1 Knapsack Problem) 優化那些事

收錄於 批改娘 20005. 0/1 Knapsack Problem。之所以有機會談到這個問題,其原因於早期的背包問題,大多都是用 branch-and-bound 算法來完成,也因此學弟課程出了這一份作業,大部分的測資,使用 branch-and-bound 能跑得比一般記憶體化 DP 快上非常多。當然,作為一個 Morris(?) 怎能允許這樣的事情發生。

現在回顧一下背包問題的模型吧!

輸入格式

每組測資第一行包含兩個正整數,分別代表背包大小 $M$ ($\leq 5×10^6$) 和物品個數 $N$ ($\leq 1000$),下一行開始每行包含兩個正整數,分別代表物品價值 $P_i$ ($\leq 10^5$)和物品重量 $W_i$ ($ \leq 10^5$)。

輸出格式

對於每組測資,請輸出最大收益。

範例輸入 1

1
2
3
4
5
6
7
8
50 7
70 31
20 10
39 20
37 19
7 4
5 3
10 6

範例輸出 1

1
107

範例輸入 2

1
2
3
4
5
6
7
8
170 7
442 41
525 50
511 49
593 59
546 55
564 57
617 60

範例輸出 2

1
1735

Solution

Branch-and-bound

bound knapscak problem 古耕竹同學提供

如果物品可以被切割,那麼可以利用物品的 CP 值 ($\textit{cp}_i = p_i/w_i$)排序,使用貪心算法在 $O(N \log N)$ 找到最佳解。然而,背包問題在於物品只能挑或不挑,一旦無法切割物品,那麼貪心算法無法將剩餘的部分填滿,進而可能產生更好的一組解填滿剩餘部分。

想要更快嗎?多看論文且實作它吧!

branch-and-bound 基本核心操作為

  • 按照 CP 值由大到小排序
  • 貪心法最佳解 $\textit{bound}$
  • 進行深度優先搜索,優先挑 CP 值大的入選
    • 當前挑選方案最佳解 $g$ + 剩餘物品使用貪心法最佳解 $g$ 小於等於當前最佳解 $\textit{bound}$,則退出搜索。
    • 更新 $\textit{bound} = \max(\textit{bound}, g)$
    • 選擇加入 或 不加入 (註:學弟說我的某些測資,通通優先不選可以快個數十倍)

結論

branch-and-bound 空間複雜度只跟 $N$ 有關,使用記憶體空間小,相較於記憶化搜索有較少的 cache miss,速度取決於搜索順序。對於同一 CP 的等價處理薄弱,一遇到這種情況,搜尋時間瞬間指數次方上去,可以等個昏天暗地。

Dynamic Programming

前言

請不要忘記背包問題屬於 NP-Complete,我們能做的事情只能優化計算,最慘的情況仍要面對,優化常數是可以努力的方向,讓我們嘗試變得更快吧。

基礎寫法解說請參考 DJWS - Bounded Knapsack Problem 的說明。

從定義上,我們通常會宣告 dp[i][j] 表示放入前 $i$ 個物品時,總共最多為 $j$ 的最大價值為何。這樣空間宣告使用 $\mathcal{O}(NW)$。在實作上,我們可以藉由運算順序將空間降為 $\mathcal{O}(W)$。因此,寫出以下代碼並不難

1
2
3
4
5
int dp[MAXW];
memset(dp, 0, sizeof(dp));
for (int i = 0; i < N; i++)
for (int j = W; j >= w[i]; j--)
dp[j] = max(dp[j], dp[j-w[i]]+v[i]);

優化初夜

從實際運行上,我們可以發現每次跑 $\mathcal{\theta}(W)$ 非常浪費,只需要跑 $\min(W, \sum w_i)$ 即可。因此,第一份計算量優化如下

1
2
3
4
5
6
7
int dp[MAXW];
memset(dp, 0, sizeof(dp));
for (int i = 0, sum = 0; i < N; i++) {
sum += w[i];
for (int j = min(W, sum); j >= w[i]; j--)
dp[j] = max(dp[j], dp[j-w[i]]+v[i]);
}

計算邊界優化通常可以達到 2x 加速

如此一來,在數量多權重小時,剛啟動的效能時可以賺到非常多。然而,不乏第一次就給權重的大的,目標最小化 $\sum \text{sum}_i$,從數學觀念很明顯地瞭解,只要一開始將權重 $w_i$ 由小到大排序即可,這樣能保證最小化計算量!

1
2
3
4
5
6
7
8
int dp[MAXW];
memset(dp, 0, sizeof(dp));
sort (w, v) by w
for (int i = 0, sum = 0; i < N; i++) {

sum += w[i];
for (int j = min(W, sum); j >= w[i]; j--)
dp[j] = max(dp[j], dp[j-w[i]]+v[i]);
}

數學使得我們更進一步,達到 1.5x 加速

優化二夜

經由平行的訓練,也許我們可以更往上一層優化。

接下來,打算把物品拆成兩堆,再利用優化初夜學到的技巧,就能引爆更多計算邊界優化。如果拆成三堆以上,合併操作變得相當複雜,當只有兩堆時,保證合併效能一定在 $\mathcal{\theta}(W)$ 完成。

如何合併兩堆的計算結果,假設 dp1[i] 表示其中一堆重量小於等於 $i$ 的最佳解,同理 dp2[j] 的計算結果。

當要湊出重量為 $W$ 的最佳解時,窮舉其中一堆的重量 $i$ 維護其中一堆的前綴最大值 $j$,相當於使用掃描線算法在線性時間內合併。合併操作如下:

1
2
3
4
5
int ret = 0;
for (int i = 0, j = W, mx = 0; i <= W; i++, j--) {
mx = max(mx, dp2[i]);
ret = max(ret, dp1[j] + mx);
}

Divide-and-Conquer,使得我們更快再更快!達到 1.5 加速

然而,優化問題將轉移到最佳分堆策略,好的分堆策略將使得計算量下降更多。目標分兩堆,使得 $\sum \text{sum1}_i + \sum \text{sum2}_i$ 最小化。明顯地,由小到大排序物品重量,依序將物品放到總和最小的那一堆即可。最後,我們整合每一夜的結果如下:

一個好的分堆,達到 1.2x 加速

相信在不久之後,還有更好的優化策略,也許不是延伸,而是全新的面貌。

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
#include <bits/stdc++.h>
using namespace std;

namespace {

static const int MAXW = 5000005;
static const int MAXN = 1005;

static void run(int A[][2], int dp[], int W, int N) {
for (int i = 0, sum = 0; i < N; i++) {
int w = A[i][0], v = A[i][1];
for (int j = min(W, sum + w); j >= w; j--)
dp[j] = max(dp[j], dp[j-w]+v);
sum += w;
}
}

static int knapsack(int C[][2], int N, int W) {
vector< pair<int, int> > A;
for (int i = 0; i < N; i++)
A.push_back(make_pair(C[i][0], C[i][1]));
N = A.size();
assert(N < MAXN);
static int dp1[MAXW+1], dp2[MAXW+1];
int Ar[2][MAXN][2], ArN[2] = {};
memset(dp1, 0, sizeof(dp1[0])*(W+1));
memset(dp2, 0, sizeof(dp2[0])*(W+1));
sort(A.begin(), A.end());
int sum[2] = {};
for (int i = 0; i < N; i++) {
int ch = sum[1] < sum[0];
Ar[ch][ArN[ch]][0] = A[i].first;
Ar[ch][ArN[ch]][1] = A[i].second;
ArN[ch]++;
sum[ch] += A[i].first;
}

run(Ar[0], dp1, W, ArN[0]);
run(Ar[1], dp2, W, ArN[1]);

int ret = 0;
for (int i = 0, j = W, mx = 0; i <= W; i++, j--) {
mx = max(mx, dp2[i]);
ret = max(ret, dp1[j] + mx);
}
return ret;
}
}

int main() {
int W, N, C[MAXN][2];
while (scanf("%d %d", &W, &N) == 2) {
for (int i = 0; i < N; i++)
assert(scanf("%d %d", &C[i][1], &C[i][0]) == 2);
printf("%d\n", knapsack(C, N, W));
}
return 0;
}
Read More +

[讀檔操作] 有限記憶體排序數組

背景

現在的系統架構中,內存 (Memory,台灣翻譯:記憶體,大陸翻譯:內存) 的大小仍無法完全像硬碟機一樣。如檔案大小 64GB,內存只有 4GB,處理檔案無法全部 In Memory,當然最近的硬體技術也逐漸朝著 All In Memory 的方式進行加速。回過頭來,在內存不足的情況下,來練習如何處理檔案大小遠大於內存的情況吧。

題目描述

給予一個 binary file 的檔案名稱,裡面用 4 個位元組為一個 signed 32-bit 整數,有數個直到 EOF,請排序後以 binary mode 導入 stdout 輸出。

輸入格式

標準輸入串流只有一行一個字串,表示檔案名稱 $F$。檔案大小最多 16 MB,而你的程序被限制最多使用 4 MB 的內存。

每個整數 $x$,限制條件 $0 \le x \le 10^9$

輸出格式

輸出在標準串流以 binary mode 下,請避開使用 console 輸出,會因為特殊字元 (如 「嗶」一聲) 導致系統嚴重當機。

範例輸入

1
1.dat

範例輸出

1
... binary mode 無法顯示

提示

download p10068-in.dat

1
2
3
4
5
6
7
8
9
morris821028@morris821028-PC MINGW64 ~/Desktop/10068
$ ./sort2.exe >test.out
1.dat

morris821028@morris821028-PC MINGW64 ~/Desktop/10068
$ xxd test.out
00000000: 0000 0000 0100 0000 0100 0000 0100 0000 ................
00000010: 0300 0000 0400 0000 0500 0000 0500 0000 ................
00000020: 0800 0000 0900 0000 ........

1.dat 以 binary file 儲存 10 個 4 bytes 整數,依序為 5, 9, 3, 5, 0, 1, 1, 8, 4, 1,排序後即為 0, 1, 1, 1, 3, 4, 5, 5, 8, 9。

由於限制內存使用量,無法寫入暫存檔案,可利用多次讀檔完成。

Solution

乍看之下,這一題是最常見的 external sort (外部排序),外部排序需要額外寫檔案,由於記憶體用量計算,一寫檔案會計算到使用的記憶體中,實際上這一題不寫檔也是能解決的。

  • 給定要排序的數據範圍
  • 二分搜尋可容納的排序範圍,利用平衡樹或者 hash 來完成計算 <某整數, 某整數出現個數>
  • 將可容納到 hash 的所有數字排序,再將其直接寫到輸出檔案。
  • 不斷地遞增二分搜尋的左邊邊界。
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
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <limits.h>

#define swap(x, y) {int t; t = x, x = y, y = t;}
int fsize(FILE *fp) {
int prev = ftell(fp);
fseek(fp, 0L, SEEK_END);
int size = ftell(fp);
fseek(fp, prev, SEEK_SET);
return size;
}
#define USAGEMEM (2<<20)
#define MAXELE ((2<<20)/8)
int A[MAXELE];
#define HSIZE 100007
#define HNODE 50000
typedef struct HashNode {
int f, cnt;
struct HashNode *next;
} HashNode;
HashNode *hhead[HSIZE], hnodes[HNODE];
int nodesize = 0;
int cmp(const void *x, const void *y) {
return ((HashNode*) x)->f - ((HashNode*) y)->f;
}
unsigned int hash(int f) {
unsigned int a = 63689, b = 378551;
unsigned int value = 0;
value = value * a + f, a = a * b;
return value;
}
void initHash() {
memset(hhead, 0, sizeof(hhead));
nodesize = 0;
}
int insertHash(int f) {
unsigned int hidx = hash(f)%HSIZE;
for (HashNode *p = hhead[hidx]; p != NULL; p = p->next) {
if (f == p->f) {
p->cnt++;
return 1;
}
}
if (nodesize >= HNODE) return 0;
HashNode s = {.f = f, .cnt = 1, .next = NULL};
hnodes[nodesize] = s;
hnodes[nodesize].next = hhead[hidx];
hhead[hidx] = &hnodes[nodesize];
nodesize++;
return 1;
}
int merge_int(FILE *fp, int l, int r) {
initHash();
int ret = 0, x, n = 0;
fseek(fp, 0, SEEK_SET);
while ((n = fread(A, sizeof(int), MAXELE, fp))) {
for (int i = 0; i < n; i++) {
if (A[i] >= l && A[i] <= r) {
if (!insertHash(A[i]))
return 0;
}
}
}
return 1;
}
#ifdef _WIN32
#include <fcntl.h>
#endif
int block_process(FILE *fp) {
int base = 0;
int l, r, mid, ret;
l = base, r = 1000000000, ret = 0;
#ifdef _WIN32
if (setmode(fileno(stdout), O_BINARY) == -1) {
perror("Cannot set stdout to binary mode");
return 2;
}
#endif
#ifdef __linux__
FILE *const out = fdopen(dup(fileno(stdout)), "wb");
#endif
while (l <= r) {
mid = (l + r)/2;
int status = merge_int(fp, base, mid);
if (status == 1) {
l = mid + 1, r = 1000000000;
base = mid + 1;
qsort(hnodes, nodesize, sizeof(HashNode), cmp);
for (int i = 0; i < nodesize; i++) {
int x = hnodes[i].f;
for (int j = hnodes[i].cnt-1; j >= 0; j--) {
#ifdef _WIN32
fwrite(&x, sizeof(int), 1, stdout);
#endif
#ifdef __linux__
fwrite(&x, sizeof(int), 1, out);
#endif
}
ret += hnodes[i].cnt;
}
} else {
r = mid - 1;
}
}
return ret;
}
int main() {
char fName[128];
scanf("%s", fName);
FILE *fin = fopen(fName, "rb");

int n = fsize(fin) / sizeof(int);
block_process(fin);
return 0;
}
Read More +

[讀檔操作] 有限記憶體找中位數

背景

現在的系統架構中,內存 (Memory,台灣翻譯:記憶體,大陸翻譯:內存) 的大小仍無法完全像硬碟機一樣。如檔案大小 64GB,內存只有 4GB,處理檔案無法全部 In Memory,當然最近的硬體技術也逐漸朝著 All In Memory 的方式進行加速。回過頭來,在內存不足的情況下,來練習如何處理檔案大小遠大於內存的情況吧。

題目描述

給予一個 binary file 的檔案名稱,裡面用 4 個位元組為一個 signed 32-bits 整數,有數個直到 EOF,保證恰好有奇數個,請輸出中位數為何?

輸入格式

標準輸入串流只有一行一個字串,表示檔案名稱 $F$。檔案大小最多 16 MB,而你的程序被限制最多使用 4 MB 的內存。

每個整數 $x$,限制條件 $0 \le x \le 10^9$

輸出格式

輸出一行整數 $Median$ 為 binary file 的中位數。

範例輸入

1
p10067-in.dat

範例輸出

1
97537111

提示

二分搜尋,分桶

download p10067-in.dat

p10067-in.dat

1
2
3
10825098
97537111
208681850

Solution

這一個問題很久以前見過,對岸通常叫做「海量資料找中位數」約束在記憶體不夠的情況下,如何高效率找到中位數,操作時允許重複讀檔。

遲遲沒有 Online Judge 進行記憶體用量檢測以及開放讀寫檔案,如今有機會擔任台灣大學計算機程式設計的課程助教,架了一個允許亂來的 Online Judge - Judge Girl (中文翻譯:批改娘系統),讀寫檔案也不會被封鎖!洽詢 網站,目前只針對課堂學生開放。

回到問題,題目給訂 16MB 的整數序列,最多 $n = 4 \times 10^6$,由於記憶體不足沒辦法直接宣告陣列大小為 $n$ 的整數陣列,又由於數字可能會重複,就沒辦法利用 bitmask 進行壓縮,只能倚靠不斷地讀檔案進行計算。若以數值範圍 $[0, V]$,大致放分成兩種策略

  • 二分搜尋 (AC, 250ms),效率 $\mathcal{O}(V \log V)$,開檔次數 $\mathcal{O}(\log V)$
  • 分桶算法 (AC, 70ms),效率 $\mathcal{O}(N)$,開檔次數 $\mathcal{O}(1)$

讀取檔案時,在有限記憶體空間,配置一塊作為緩衝區,一次讀入一坨子序列,一個一個讀取效率非常差,別輕忽從磁碟讀取資料的速度。剩餘空間再進行紀錄用途。從效能比較 二分搜尋 慢於 分桶算法

二分搜尋

二分答案 $x$,接著線性掃描序列,計算有多少個整數小於等於 $x$,藉此找到中位數。

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 <stdio.h>

int fsize(FILE *fp) {
int prev = ftell(fp);
fseek(fp, 0L, SEEK_END);
int size = ftell(fp);
fseek(fp, prev, SEEK_SET);
return size;
}

#define MAXELE (3<<20)/4
int countLess(FILE *fp, int ans) {
static int A[MAXELE];
int ret = 0, x, n = 0;
fseek(fp, 0, SEEK_SET);
while ((n = fread(A, sizeof(int), MAXELE, fp))) {
for (int i = 0; i < n; i++)
ret += A[i] <= ans;
}
return ret;
}

int main() {
char fName[128];
scanf("%s", fName);
FILE *fin = fopen(fName, "rb");

int n = fsize(fin) / sizeof(int);

int m = n/2 + 1;
int l = 0, r = 1000000000, mid, ret = 0;

while (l <= r) {
mid = (l + r)/2;
int tn = countLess(fin, mid);
if (tn < m)
l = mid + 1, ret = mid+1;
else
r = mid - 1;
}
printf("%d\n", ret);
return 0;
}

分桶算法

分桶算法分成兩次掃瞄,假設內存分配最多 $I$ 個 32bits 整數 ($I \ll N$)。

  • 由於記憶體限制大小關係,將數值範圍分成 $I$ 堆,如 [0, V/I-1][V/I, 2V/I-1]、… 等,計算區間範圍內的個數,讀取一次檔案,可以得到中位數介於 [iV/I, (i+1)V/I-1]
  • 接著再劃分一次,直到區間長度為 1,中位數便可得到。
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
#include <stdio.h>
#include <assert.h>
#include <string.h>
int fsize(FILE *fp) {
int prev = ftell(fp);
fseek(fp, 0L, SEEK_END);
int size = ftell(fp);
fseek(fp, prev, SEEK_SET);
return size;
}

#define MAXELE ((1<<20)/4)
static int A[MAXELE];
static int B[MAXELE];
int countLess(FILE *fp, int ans) {
int ret = 0, x, n = 0;
fseek(fp, 0, SEEK_SET);
while ((n = fread(A, sizeof(int), MAXELE, fp))) {
for (int i = 0; i < n; i++)
ret += A[i] <= ans;
}
return ret;
}

int main() {
char fName[128];
scanf("%s", fName);
FILE *fin = fopen(fName, "rb");

int n = fsize(fin) / sizeof(int);

int m = n/2 + 1;
int l = 0, r = 1000000000, mid, ret = 0;

int WIDTH = r, SHIFT;

for (SHIFT = 1; (1<<SHIFT) < MAXELE; SHIFT++);
WIDTH = r / (1<<SHIFT);

memset(B, 0, sizeof(B));
fseek(fin, 0, SEEK_SET);
while ((n = fread(A, sizeof(int), MAXELE, fin))) {
for (int i = 0; i < n; i++) {
assert((A[i]>>SHIFT) < MAXELE);
B[A[i]>>SHIFT]++;
}
}

int BASE = 0, SELECT;
for (int i = 0, sum = m; i < MAXELE; i++) {
sum -= B[i];
if (sum <= 0) {
BASE = i * (1<<SHIFT);
SELECT = sum + B[i];
break;
}
}

memset(B, 0, sizeof(B));
fseek(fin, 0, SEEK_SET);
while ((n = fread(A, sizeof(int), MAXELE, fin))) {
for (int i = 0; i < n; i++) {
if (A[i] - BASE < MAXELE && A[i] >= BASE)
B[A[i] - BASE]++;
}
}
for (int i = 0, sum = SELECT; i < MAXELE; i++) {
sum -= B[i];
if (sum <= 0) {
printf("%d\n", i + BASE);
break;
}
}
return 0;
}
Read More +

[ACM 題目] 障礙物轉換

Problem

在計算幾何領域中,機器人規劃路線是一個實用的議題, 由於機器人是一個有體積的物體,要避開其他的障礙物,障礙物碰撞計算會變得相當複雜。由於某 M 沒有學好這一部分,知道一種方法可以將地圖轉換,把機器人轉換成一個點,計算得到新的障礙物,如此一來就不用處理障礙物碰撞。

從上圖中,一個箭頭型機器人要從左下角移動到右上角,假設不考慮物體可以旋轉,請問是否能移動過去。若把機器人當作一點看待,則必須將物體邊緣增加,變成中間那張圖 configuration space,只剩下一個點和數個多邊形障礙物。接著可以轉換成 visibility graph,把剩下的白色空間進行梯形剖分或者三角剖分,將區域表示成一個點,相鄰區域拉一條邊,就成 dual graph,就能套用一般的最短路徑算法,來協助機器人行走。

處理這問題對於某 M 過於困難,於是將問題更簡化些,只需要把中間那一張圖其中一個障礙物變化求出即可。

現在給定一個凸多邊形障礙物以及移動凸多邊形,假設定位點在深黑色點方形部分,藉由貼著障礙物可以構出新的障礙物,請問新的障礙物為何?

Input

第一行會有一個整數 $T$,表示有多少組測資。

每一組測資會包含兩個凸多邊形。第一行有一個整數 $N$ 表示移動凸多邊形的頂點數量,接下來會有 $N$ 行,每一行會有兩個整數 $x, y$ 表示頂點座標。在 $N+1$ 行之後,會有一個整數 $M$ 表示障礙物凸多邊形的頂點數量,接下來會有 $M$ 行,每一行會有兩個整數 $x, y$ 表示頂點座標。

  • $2 < N, M < 512$
  • $-32767 < x, y < 32767$
  • 假設定位點 (深黑色點方形部分) 都設在原點 (0, 0)。
  • 輸入順序會按照順時針或逆時針給定。

Sample Input

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
2 

5
-1 -1
-1 1
0 3
1 1
1 -1

4
1 1
8 2
12 -1
6 -4

3
-1 -4
-1 1
1 1

5
1 0
6 5
10 5
10 -3
1 -3

Sample Output

1
2
Case #1: 89.0
Case #2: 115.5

Hint

Solution

可以用 $O(nm \log nm)$ 來完成此題,方法很簡單,先將機器人的圖形按照原點旋轉 180 度,接著將障礙物每一個頂點座標加上機器人的座標,總共會得到 $nm$ 的頂點,套用凸包算法求一次即可。

這一題應用 Minkowski Sum 來完成,對於兩個凸多邊形的 Minkowski Sum,可以在 $O(n+m)$ 時間內完成,將兩個凸包採用逆時針儲存,挑選最左,等價時抓下方 y 座標的點,可以知道新凸包的頂點一定是由兩個凸包的頂點$A_i, B_j$ 疊加而成,一開始抓到的左下方頂點疊加可以得到第一個頂點,接著要按照凸包順序求出,則比較$(A_i, A_{i+1})$$(B_j, B_{j+1})$ 哪一個偏的角度小,就挑選哪一個往前推動 i = i+1 或者是 j = j+1。最後就能在 $O(n+m)$ 時間內完成。

假設不是凸多邊形,兩個簡單多邊形最慘會到 $O(n^2 m^2)$

1

Read More +

[ACM 題目] 優勢商品

Problem

一份產品有很多屬性,把每一個屬性當作一個維度,當成 $D$ 個維度的向量,當維度越多時絕對支配 (各方面都是最好的) 的機率就越低。若把 $D = 2$,就是在二維空間中的支配點 (Zerojudge d555 平面上的極大點)。由於絕對優勢的產品隨著維度增加而變少,定義出一種相對優勢產品,給定一個 $K$,若 $p_i$ 支配 (k-dominate) $p_j$,必須滿足下列三條:

  • $\exists S' \subseteq S, |S'| = K$
  • $\forall s_r \in S', p_i.s_r \geq p_j.s_r$
  • $\exists s_t \in S', p_i.s_t > p_j.s_t$

用中文解釋這幾條數學公式,假設有 6 個維度的產品,目標 $K = 5$,則表示要在 6 的維度中任意挑選 5 個維度比較,若存在一種挑法可以支配,則可以說 $p_i$ k-dominate $p_j$

例如現在有 8 台筆記型電腦,共有 6 種屬性可以比較,數值越大越好。

Notebook CPU RAM HDD HDD S. Q. VRAM
$N_1$ 3 3 5 6 6 8
$N_2$ 9 4 9 7 7 9
$N_3$ 8 4 7 9 2 7
$N_4$ 5 6 8 9 5 9
$N_5$ 9 7 9 6 2 4
$N_6$ 6 6 6 5 3 5
$N_7$ 5 7 3 8 4 6
$N_8$ 4 4 8 6 6 4

若要判斷 $N_2$ 是否支配 $N_3$,從中挑取 (CPU, RAM, HDD, HDD S., Q.) 這 5 個屬性比較,得到 $N_2$ 勝過於 $N_3$。最後可以發現到 $N_2$ 這產品可以支配 $N_1, N_3, N_5, N_6, N_8$

現在問題來了,要找到不被任何產品支配的產品 (k-dominant skyline)。如上表 8 項產品的情況,只會有 $N_2, N_4$

Input

第一行會有一個整數 $T$ 表示有多少組測資。

每一組測資第一行會有三個整數 $N, D, K$,分別表示產品數量、產品屬性種類、以及支配比較的屬性種類。接下來會有 $N$ 行數據,每一行上會有 $D$ 個整數 $s_j$ 表示一個產品的屬性,產品按照輸入順序編號。

  • $N < 1024$
  • $0 < D, K < 8$
  • $0 < s_j < 1000$

Output

對於每組測資輸出一行,輸出該組測資所有的 k-dominant skyline 的產品編號,由小到大輸出。若不存在則輸出 NONE

Sample Input

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
4

2 3 1
20 20 20
20 20 20

2 5 1
20 5 20 20 20
5 20 20 20 20

4 6 5
9 4 9 7 7 9
9 7 9 6 5 4
9 6 9 8 3 3
1 2 3 9 2 2

8 6 5
3 3 5 6 6 8
9 4 9 7 7 9
8 4 7 9 2 7
5 6 8 9 5 9
9 7 9 6 2 4
6 6 6 5 3 5
5 7 3 8 4 6
4 4 8 6 6 4

Sample Output

1
2
3
4
Case #1: 1 2
Case #2: NONE
Case #3: 1
Case #4: 2 4

Solution

這是一篇論文題目 k-dominant skyline set,網路上提出不少算法,看起來都是有容錯率。

這一題只是介紹題,沒有強求高效算法計算,直接 $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
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
#include <stdio.h>
#include <stdlib.h>
#include <set>
#include <queue>
#include <string.h>
#include <algorithm>
using namespace std;
const int MAXN = 1024;
int N, D, K;
vector< vector<int> > C;
struct Product {
int v[8], id;
bool operator<(const Product &a) const {
return sum() > a.sum();
}
int sum() const {
int s = 0;
for (int i = 0; i < D; i++)
s += v[i];
return s;
}
int domain(Product &u) {
for (auto &x : C) {
int h1 = 1, h2 = 0;
for (auto &e : x) {
h1 &= v[e] >= u.v[e];
h2 |= v[e] > u.v[e];
}
if (h1 && h2) return 1;
}
return 0;
}
} item[MAXN];
int used[MAXN];
int main() {
int testcase, cases = 0;
scanf("%d", &testcase);
while (testcase--) {
scanf("%d %d %d", &N, &D, &K);
for (int i = 0; i < N; i++) {
for (int j = 0; j < D; j++)
scanf("%d", &item[i].v[j]);
item[i].id = i;
}
sort(item, item+N);

C.clear();
for (int i = 0; i < (1<<D); i++) {
if (__builtin_popcount(i) == K) {
vector<int> A;
for (int j = 0; j < D; j++) {
if ((i>>j)&1)
A.push_back(j);
}
C.push_back(A);
}
}

vector<Product> filter;
vector<int> ret;
for (int i = 0; i < N; i++)
used[i] = 0;
for (int i = 0; i < N; i++) {
if (used[i] == 0) {
filter.push_back(item[i]);
for (int j = i+1; j < N; j++) {
if (used[j] == 0 && item[i].domain(item[j]))
used[j] = 1;
}
}
}
for (int i = 0; i < filter.size(); i++) {
int ok = 1;
for (int j = 0; j < N && ok; j++) {
if (item[j].domain(filter[i]))
ok = 0;
}
if (ok == 1)
ret.push_back(filter[i].id);
}

sort(ret.begin(), ret.end());
printf("Case #%d:", ++cases);
for (int i = 0; i < ret.size(); i++)
printf(" %d", ret[i]+1);
puts(ret.size() ? "" : " NONE");
}
return 0;
}
Read More +

[ACM 題目] 竹馬不敵天降

Problem

「話說,剛剛明明說的是合乎聖誕節的商品。這不是 H-GAME 嗎?哪裡有聖誕元素啊!」
「嗯,真虧你發現這點啊,但還是有點可惜,這是全齡版,沒有工口哦!」
「不管怎麼樣,一點也沒合聖誕節啊。」
「雖然你有好好學習許多東西,但對本質上完全沒有理解。」
「最近開始覺得有點明白了 …」
「那麼就來說明下這個作品吧」
「是一年前發售的大人氣美少女遊戲的續篇吧?」
「是的,那為什麼你沒有察覺到-『一年前』這句話所包含的意義」
「對於期盼著的人那就意味著 … 與戀人時隔一年的再會!」

在 GALGAME 故事發展過程中,竹馬幾乎沒勝過天降,走哪一條線進行發展正是玩 GALGAME 的樂趣。

然而有一款極為糟糕的 GALGAME 的初始方式則是在一開始選定座標下,系統會根據座標找到附近最鄰近少女,並且在隨後的故事情節都會圍繞這名少女。

遊戲的地圖設計很簡單,用一個矩形表示,現在已知 N 名少女的所在地 (都在矩形內部)。一開始選定的座標也必須落在矩形內,請問玩哪每個線路機率分布為何?

Input

輸入有多組測資。

每組第一行會有三個整數 N, A, B,表示在$[0:A] \times [0:B]$ 地圖中有 N 個已知座標。
接下來會有 N 行,每行上會有兩個整數 x, y,表示每名少女所在的座標$(x, y)$ 保證不會有重疊的情況。

$(0 < N \le 100, 0 < A, B \le 1000, 0 \le x \le A, 0 \le y \le B)$

Output

對於每組測資輸出一行,根據輸入順序輸出每一個故事線的機率。

Sample Input

1
2
3
4
5
6
7
8
9
10
11
12
2 5 3
1 2
4 2

2 5 3
1 1
2 2

3 5 3
1 1
2 2
4 1

Sample Output

1
2
3
Case 1: 0.500 0.500
Case 2: 0.300 0.700
Case 3: 0.292 0.313 0.396

Solution

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
#include <stdio.h> 
#include <math.h>
#include <algorithm>
#include <set>
#include <vector>
using namespace std;

#define eps 1e-6
#define MAXN 32767
struct Pt {
double x, y;
Pt(double a = 0, double b = 0):
x(a), y(b) {}
Pt operator-(const Pt &a) const {
return Pt(x - a.x, y - a.y);
}
bool operator<(const Pt &a) const {
if (fabs(x - a.x) > eps)
return x < a.x;
if (fabs(y - a.y) > eps)
return y < a.y;
return false;
}
};
double dot(Pt a, Pt b) {
return a.x * b.x + a.y * b.y;
}
double cross(Pt o, Pt a, Pt b) {
return (a.x-o.x)*(b.y-o.y)-(a.y-o.y)*(b.x-o.x);
}
int between(Pt a, Pt b, Pt c) {
return dot(c - a, b - a) >= -eps && dot(c - b, a - b) >= -eps;
}
int onSeg(Pt a, Pt b, Pt c) {
return between(a, b, c) && fabs(cross(a, b, c)) < eps;
}
struct Seg {
Pt s, e;
double angle;
int label;
Seg(Pt a = Pt(), Pt b = Pt(), int l=0):s(a), e(b), label(l) {
angle = atan2(e.y - s.y, e.x - s.x);
}
bool operator<(const Seg &other) const {
if (fabs(angle - other.angle) > eps)
return angle > other.angle;
if (cross(other.s, other.e, s) > -eps)
return true;
return false;
}
};
Pt getIntersect(Seg a, Seg b) {
double a1, b1, c1, a2, b2, c2;
double dx, dy, d;
a1 = a.s.y - a.e.y, b1 = a.e.x - a.s.x;
c1 = a1 * a.s.x + b1 * a.s.y;
a2 = b.s.y - b.e.y, b2 = b.e.x - b.s.x;
c2 = a2 * b.s.x + b2 * b.s.y;
dx = c1 * b2 - c2 * b1, dy = a1 * c2 - a2 * c1;
d = a1 * b2 - a2 * b1;
return Pt(dx/d, dy/d);
}
Seg deq[MAXN];
vector<Pt> halfPlaneIntersect(vector<Seg> segs) {
sort(segs.begin(), segs.end());
int n = segs.size(), m = 1;
int front = 0, rear = -1;
for (int i = 1; i < n; i++) {
if (fabs(segs[i].angle - segs[m-1].angle) > eps)
segs[m++] = segs[i];
}
n = m;
deq[++rear] = segs[0];
deq[++rear] = segs[1];
for (int i = 2; i < n; i++) {
while (front < rear && cross(segs[i].s, segs[i].e, getIntersect(deq[rear], deq[rear-1])) < -eps)
rear--;
while (front < rear && cross(segs[i].s, segs[i].e, getIntersect(deq[front], deq[front+1])) < -eps)
front++;
deq[++rear] = segs[i];
}
while (front < rear && cross(deq[front].s, deq[front].e, getIntersect(deq[rear], deq[rear-1])) < -eps)
rear--;
while (front < rear && cross(deq[rear].s, deq[rear].e, getIntersect(deq[front], deq[front+1])) < -eps)
front++;

vector<Pt> ret;
for (int i = front; i < rear; i++) {
Pt p = getIntersect(deq[i], deq[i+1]);
ret.push_back(p);
}
if (rear > front + 1) {
Pt p = getIntersect(deq[front], deq[rear]);
ret.push_back(p);
}
return ret;
}
double calc_convexarea(const vector<Pt> &p) {
double ret = 0;
int n = p.size();
for(int i = 0, j = n - 1; i < n; j = i++)
ret += p[j].x * p[i].y - p[j].y * p[i].x;
return fabs(ret /2.0);
}
int main() {
int n, A, B, cases = 0;
Pt p[128], mid, vij;
while (scanf("%d %d %d", &n, &A, &B) == 3) {
for (int i = 0; i < n; i++)
scanf("%lf %lf", &p[i].x, &p[i].y);
printf("Case %d:", ++cases);
for (int i = 0; i < n; i++) {
int m = 0;
vector<Seg> segs;
segs.resize(n - 1 + 4);
for (int j = 0; j < n; j++) {
if (i == j) continue;
mid.x = (p[i].x + p[j].x) /2.0;
mid.y = (p[i].y + p[j].y) /2.0;
vij.x = mid.x - (p[j].y - p[i].y);
vij.y = mid.y + (p[j].x - p[i].x);
segs[m] = Seg(mid, vij), m++;
}
segs[m++] = Seg(Pt(0, 0), Pt(A, 0));
segs[m++] = Seg(Pt(A, 0), Pt(A, B));
segs[m++] = Seg(Pt(A, B), Pt(0, B));
segs[m++] = Seg(Pt(0, B), Pt(0, 0));
vector<Pt> convex = halfPlaneIntersect(segs);
// for (int j = 0; j < convex.size(); j++)
// printf("%lf %lf\n", convex[j].x, convex[j].y);
// puts("--");
printf(" %.3lf", calc_convexarea(convex) / (A * B));
}
puts("");
}
return 0;
}
/*
2 5 3
1 2
4 2

2 5 3
1 1
2 2

3 5 3
1 1
2 2
4 1
*/
Read More +

[ACM 題目] 順來逆受

Problem

「其實雙親不准我在家裡畫漫畫。」
「是嗎?你都畫些什麼啊」
「乃是在下 x 大哥的原創本是也。」
「還有,薰是普通人」
「薰姐,功的反義詞是什麼?」
「什麼?是守嗎?」
「薰姐,剛才我說的話請你全忘記吧。」

在一個 N 個城市的國家,每個城市之間用 M 條邊相連。隨著時間,它們開始分裂,彼此會開始斷訊,無法聯絡到相鄰的城市。定時回報某個城市可以藉由間接或直接連絡的城市個數。

操作:

  • D i:移除輸入順序 i 個連接邊。
  • Q i:輸出城市 i 能聯絡的城市個數。

Sample Input

1
2
3
4
5
6
7
8
9
10
11
5 5
1 2
1 3
2 3
3 4
4 5
4
Q 1
D 4
Q 5
Q 2

Sample Output

1
2
3
5
2
3

Solution

1

Read More +

[ACM 題目] 動態前綴

Problem

某 M 「給一個字串 S,接著詢問操作 [l, r] 告訴區間內是否剛好為一個迴文?」
學弟 「用最長迴文的方式去想嗎?」
某 M 「我沒有標準答案哦 wwww,吾等還在想能不能離線 O(1)。」

暗自敲著 UVa 另外一題類似題說著,就算能 O(1) 檢查迴文,區間還是要窮舉 O(n),在此宣告某 M 陣亡- Time Limit

某 M 心想,假解也是不錯選擇,單純詢問迴文肯定會被 manacher’s algorithm 秒掉 (O(1) - O(n))。如果詢問 [l, r] 的子字串是否符合 pattern AA (例如 abcabc,其中 A = abc),也會被 suffix array、suffix tree 解決 (O(log n) - O(n log n)/O(n))。

百般無奈之下,也許這題就這樣定好了:

C p x:將字串位置 p 修改成字符 x。
Q p q:求子字串 S.substr(p) 和 S.substr(q) 的最長共同前綴長度。

1
2
3
4
5
6
7
8
void change(char s[], int p, int x) {
s[p] = x;
}
int query(char s[], int p, int q) {
int ret = 0;
for (; s[p] == s[q] && s[p] && s[q]; p++, q++, ret++);
return ret;
}

Sample Input

1
2
3
4
5
abcdabcc
3
Q 0 4
C 3 c
Q 0 4

Sample Output

1
2
3
4

Solution

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
#include <stdio.h>
#include <string.h>
#include <algorithm>
#include <assert.h>
using namespace std;
const long long mod = 100000007LL;
#define MAXN (1<<17)
char S[MAXN];
long long base[MAXN];
long long tree[MAXN<<1];
long long build(int k, int l, int r) {
assert(l <= r);
if (l == r)
return tree[k] = (S[l] % mod);
int m = (l + r)/2;
build(k<<1, l, m);
build(k<<1|1, m+1, r);
tree[k] = (tree[k<<1] * base[r - m] + tree[k<<1|1]) %mod;
}
void modify(int k, int l, int r, int x, int v) {
if (l == r) {
tree[k] = v;
return;
}
int m = (l + r)/2;
if (x <= m)
modify(k<<1, l, m, x, v);
else
modify(k<<1|1, m+1, r, x, v);
tree[k] = (tree[k<<1] * base[r - m]%mod + tree[k<<1|1]) %mod;
}
long long query(int k, int l, int r, int x, int y) {
assert(l >= 0 && l <= r);
if (x <= l && r <= y) {
return tree[k];
}
int m = (l + r)/2;
if (y <= m)
return query(k<<1, l, m, x, y);
else if(x > m)
return query(k<<1|1, m+1, r, x, y);
else {
long long p, q;
p = query(k<<1, l, m, x, y);
q = query(k<<1|1, m+1, r, x, y);
return (p * base[min(y, r) - m]%mod + q) %mod;
}
}
int main() {
freopen("in.txt", "r+t", stdin);
freopen("out.txt", "w+t", stdout);
base[0] = 1;
for (int i = 1; i < MAXN; i++)
base[i] = (base[i-1] * 2)%mod;
char cmd[8], s[8];
int Q, p, q, n;
while (scanf("%s", S) == 1) {
n = strlen(S);
build(1, 0, n - 1);
scanf("%d", &Q);
for (int i = 0; i < Q; i++) {
scanf("%s", cmd);
if (cmd[0] == 'Q') {
scanf("%d %d", &p, &q);
if (p == q) {
printf("%d\n", n - p);
continue;
} else if (S[p] != S[q]) {
puts("0");
continue;
}
int l = 0, r = min(n - p, n - q) - 1, m, ret = 0;
long long hp, hq;
while (l <= r) {
m = (l + r)/2;
hp = query(1, 0, n-1, p, p + m);
hq = query(1, 0, n-1, q, q + m);
if (hp == hq)
l = m + 1, ret = m;
else
r = m - 1;
}
printf("%d\n", ret + 1);
} else {
scanf("%d %s", &p, s);
S[p] = s[0];
modify(1, 0, n - 1, p, s[0]);
}
}
}
return 0;
}
Read More +