UVa 13154 - Extreme XOR Sum

Problem

給一長度為 $N$ 的整數序列,詢問任意區間的極端異或和 Extreme XOR Sum。

定義 Extreme XOR Sum 為一系列操作的最後一個值

  • 當長度 $n>1$ 時,將陣列縮小為 $n-1$
  • $[a_0, a_1, a_2, \cdots, a_{n-1}]$,每一個元素與後一個元素運行互斥或,將會被轉換成 $[a_0 \oplus a_1, a_1\oplus a_2, a_2 \oplus a_3, \cdots, a_{n-2}\oplus a_{n-1}]$
  • 直到只剩下一個元素,即為 Extreme XOR Sum

Sample Input

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

Sample Output

1
2
3
4
Case 1:
1
5
14

Solution

這題詢問次數非常多,一般運行將對每一個詢問達到 $O(N)$ 的複雜度,這很容易得到 TLE。從大多的數據結構,如線段樹、塊狀表 … 等,他們提供高效率的查找效能,但也必須符合某些條件才能使用。因此,在此題若要符合結合律將變得相當困難。

觀察

假設要計算陣列 $[1, 4, 6, 7, 8]$ 的值時

  • 第一步,$[1 \oplus 4, 4 \oplus 6, 6 \oplus 7, 7 \oplus 8]$
  • 第二步,$[1 \oplus 4 \oplus 4 \oplus 6, \cdots]$
  • 如此類推下去,XOR 有結合律,我們發現到各別使用了 1 次 $a_0$、4 次 $a_1$、6 次 $a_2$、4 次 $a_3$ 和 1 次 $a_4$

對於不同的長度,我們發現到是二項係數的配對情況。由於偶數次的 XOR 會互消,只需要計算出現奇數次的即可,因此我們列出二項次數模二的情況,進而得到 Sierpinski triangle/Sierpinski sieve。即使知道 Sierpinski sieve 是二項係數模二的結果,我們仍不知道要怎麼套用結合律達到剖分加速的條件。

二項係數的公式如下

$$\begin{align*} \binom{n}{m} &= \binom{n-1}{m-1} + \binom{n-1}{m} \\ &= \frac{n!}{m!(n-m)!} \end{align*}$$

階層運算在數學運算上的性質並不多,逼得我們只好往碎形上觀察,以下列出前幾項的結果

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
1
1 1
1 0 1
1 1 1 1
1 0 0 0 1
1 1 0 0 1 1
1 0 1 0 1 0 1
1 1 1 1 1 1 1 1
1 0 0 0 0 0 0 0 1
1 1 0 0 0 0 0 0 1 1
1 0 1 0 0 0 0 0 1 0 1
1 1 1 1 0 0 0 0 1 1 1 1
1 0 0 0 1 0 0 0 1 0 0 0 1
1 1 0 0 1 1 0 0 1 1 0 0 1 1
1 0 1 0 1 0 1 0 1 0 1 0 1 0 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

發現它是一個很有趣的碎形,每個三角形大小都是以二的冪次的。我們按照 $2^3 = 8$ 切割一下上圖,並且把右邊斜的補上 0 得到下圖。

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
1
1 1
1 0 1
1 1 1 1
1 0 0 0 1
1 1 0 0 1 1
1 0 1 0 1 0 1
1 1 1 1 1 1 1 1
^---------------
1 0 0 0 0 0 0 0 |1 0 0 0 0 0 0 0
1 1 0 0 0 0 0 0 |1 1 0 0 0 0 0 0
1 0 1 0 0 0 0 0 |1 0 1 0 0 0 0 0
1 1 1 1 0 0 0 0 |1 1 1 1 0 0 0 0
1 0 0 0 1 0 0 0 |1 0 0 0 1 0 0 0
1 1 0 0 1 1 0 0 |1 1 0 0 1 1 0 0
1 0 1 0 1 0 1 0 |1 0 1 0 1 0 1 0
1 1 1 1 1 1 1 1 |1 1 1 1 1 1 1 1
^----------------^--------------
1 0 0 0 0 0 0 0 |0 0 0 0 0 0 0 0 |1 0 0 0 0 0 0 0
1 1 0 0 0 0 0 0 |0 0 0 0 0 0 0 0 |1 1 0 0 0 0 0 0
1 0 1 0 0 0 0 0 |0 0 0 0 0 0 0 0 |1 0 1 0 0 0 0 0
1 1 1 1 0 0 0 0 |0 0 0 0 0 0 0 0 |1 1 1 1 0 0 0 0
1 0 0 0 1 0 0 0 |0 0 0 0 0 0 0 0 |1 0 0 0 1 0 0 0
1 1 0 0 1 1 0 0 |0 0 0 0 0 0 0 0 |1 1 0 0 1 1 0 0
1 0 1 0 1 0 1 0 |0 0 0 0 0 0 0 0 |1 0 1 0 1 0 1 0
1 1 1 1 1 1 1 1 |0 0 0 0 0 0 0 0 |1 1 1 1 1 1 1 1
^ ^ ^
箭頭表示本身也是 Sierpinski sieve

區塊縮影得到 miniature pattern 也是 Sierpinski sieve
1
1 1
1 0 1

得到數個一模一樣的子圖,上述全零和非零的區塊,又恰好構成 Sierpinski sieve。這告訴我們任何操作全都要以二的冪次為基準,且合併區段時須以二項係數為係數。設定 pattern 大小為 $M=2^{\lceil \log_2 N\rceil}$,最後得到 miniature pattern。在同一層中,非零構成的條紋都是相同的模式,例如上述得圖中,最後一層的箭號組合必然是 101 或者是 000,最後得到下列公式計算條紋。

$A_{i, j} = A_{i-1}{j} \oplus A_{i-1,j+M}$

接下來,我們將需要確定每一個條紋 (stripe) 是否使用全零或者非零,只需要查找 miniature pattern 相應的係數即可。

如何在 Sierpinski sieve 找到非零係數的位置

$\binom{n}{i} \mod 2 = 1$,必滿足 $n\&i = i$。其證明從數學歸納法來,由二冪次的長度碎形著手,移除最高位的 1 得到 $i'$,從 $i'$ 舊有位置集合,保留此集合,並對每一個元素增加二的冪次得到碎形的另一邊。

故可利用下述算法,準確地找到每一個非零的係數位置

1
2
for (int pos = n; pos; pos = (pos-1)&n)
C[n][pos] mod 2 = 1

最後附上優化後得到 Rank 1 的程序 0.040 s

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

static const int M = (1<<7);
static const int MAXN = 10005;
static int A[M+1][MAXN];
void miniature(int n) {
for (int i = 1; i*M < n; i++) {
for (int j = 0; j+i*M < n; j++)
A[i][j] = A[i-1][j] ^ A[i-1][j+M];
}
}
int extract(int l, int r) {
const int n = r-l;
const int m = n/M;
const int o = n%M;
int ret = A[m][l];
for (int i = o; i; i = (i-1)&o)
ret ^= A[m][l+i];
return ret;
}
namespace MM {
inline int readchar() {
const int N = 1048576;
static char buf[N];
static char *p = buf, *end = buf;
if(p == end) {
if((end = buf + fread(buf, 1, N, stdin)) == buf) return EOF;
p = buf;
}
return *p++;
}
inline int ReadInt(int *x) {
static char c, neg;
while((c = readchar()) < '-') {if(c == EOF) return 0;}
neg = (c == '-') ? -1 : 1;
*x = (neg == 1) ? c-'0' : 0;
while((c = readchar()) >= '0')
*x = (*x << 3) + (*x << 1) + c-'0';
*x *= neg;
return 1;
}
class Print {
public:
static const int N = 1048576;
char buf[N], *p, *end;
Print() {
p = buf, end = buf + N - 1;
}
void printInt(int x, char padding) {
static char stk[16];
int idx = 0;
stk[idx++] = padding;
if (!x)
stk[idx++] = '0';
while (x)
stk[idx++] = x%10 + '0', x /= 10;
while (idx) {
if (p == end) {
*p = '\0';
printf("%s", buf), p = buf;
}
*p = stk[--idx], p++;
}
}
void flush() {
*p = '\0', p = buf;
printf("%s", buf);
}
static inline void online_printInt(int x) {
static char ch[16];
static int idx;
idx = 0;
if (x == 0) ch[++idx] = 0;
while (x > 0) ch[++idx] = x % 10, x /= 10;
while (idx)
putchar(ch[idx--]+48);
}
~Print() {
*p = '\0';
printf("%s", buf);
}
} bprint;
}
int main() {
int testcase, cases = 0;
int n, m;
// scanf("%d", &testcase);
MM::ReadInt(&testcase);
while (testcase--) {
// scanf("%d", &n);
MM::ReadInt(&n);
for (int i = 0; i < n; i++)
// scanf("%d", &A[0][i]);
MM::ReadInt(&A[0][i]);
miniature(n);
// scanf("%d", &m);
MM::ReadInt(&m);
printf("Case %d:\n", ++cases);
for (int i = 0; i < m; i++) {
int l, r;
// scanf("%d %d", &l, &r);
MM::ReadInt(&l), MM::ReadInt(&r);
// printf("%d\n", extract(l, r));
MM::bprint.printInt(extract(l, r), '\n');
}
MM::bprint.flush();
}
return 0;
}
/*
1
5
1 4 6 7 8
3
0 0
0 1
2 4
*/

Read More +

關於高效大數除法的那些事

前情提要

從國小開始學起加減乘除,腦海裡計算加法比減法簡單,減法可以換成加法,乘法則可以換成數個加法。即使到了中國最強的利器——珠心算,我們學習這古老的快速運算時,也都採取這類方法,而除法最為複雜,需要去估算商,嘗試去扣完,再做細微調整。然而,當整數位數為 $N$ 時,加減法效能為 $N$ 基礎運算,而乘除法為 $N^2$ 次基礎運算,一次基礎運算通常定義在查表,例如背誦的九九乘法表,使得我們在借位和乘積運算達到非常快速。

這些方法在計算機發展後,硬體實作大致使用這些算法,直到了快速傅立葉 (FFT) 算法能在 $O(N \log N)$ 時間內完成旋積 (卷積) 計算,順利地解決了多項式乘法的問題,使得大數乘法能在 $O(N \log N)$ 時間內完成。

一開始我們知道乘法和除法都可以在 $O(N^2)$ 時間內完成,有了 FFT 之後,除法是不是也能跟乘法一樣在 $O(N \log N)$ 內完成?

大數除法

就目前研究來看,除法可以轉換成乘法運算,在數論的整數域中,若存在乘法反元素,除一個數相當於成一個數,而我們發展的計算機,有效運算為 32/64-bit,其超過的部分視為溢位 (overflow),溢位則可視為模數。因此,早期 CPU 設計中,除法需要的運行次數比乘法多一些,編譯器會將除法轉換乘法運算 (企圖找到反元素),來加速運算。現在,由於 intel 的黑魔法,導致幾乎一模一樣快而沒人注意到這些瑣事。

回過頭來,我們介紹一個藉由 FFT 達到的除法運算 $O(N \log N \log N)$ 的算法。而那多的 $O(\log N)$ 來自於牛頓迭代法。目標算出 $C = \lfloor A/B \rfloor$,轉換成 $C = \lfloor A \ast (1/B) \rfloor$,快速算出 $1/B$ 的值,採用牛頓迭代法。

額外要求整數長度 $n = \text{length}(A)$$m = \text{length}(B)$,當計算 $1/B$ 時,要求精準度到小數點後 $n$ 位才行。

牛頓迭代法

目標求 $f(x) = 0$ 的一組解。

$$\begin{align} x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)} \end{align}$$

不斷地迭代逼近找解。若有多組解,則依賴初始值 $x_0$ 決定優先找到哪一組解。

這部分也應用在快速開根號——反平方根快速演算法 Fast inverse square root

倒數計算

藉由下述定義,找到 $x = 1/B$

$$\begin{align} f(x) = \frac{1}{x} - B, \; f'(x) = -\frac{1}{x^2} \end{align}$$

接著推導牛頓迭代法的公式

$$\begin{align} x_{n+1} &= x_n - \frac{f(x_n)}{f'(x_n)} \\ &= x_n - \frac{\frac{1}{x} - B}{-\frac{1}{x^2}} = (2 - x_n B) x_n \end{align}$$

運行範例

$A = 7123456, \; B = 123$,計算倒數 $1/B = 1/123$,由於 $A$$n=7$ 位數,小數點後精準 7 位,意即迭代比較小數點後 7 位不再變動時停止。

以下描述皆以 十進制 (在程式運行時我們可能會偏向萬進制,甚至億進制來充分利用 32-bit 暫存器)

  • 決定 $x_0$ 是很重要的一步,這嚴重影響著迭代次數。
  • $x_0$$B$ 的最高位數 1 進行大數除小數,得到 $x_0 = 0.01$
  • $x_1 = (2 - x_0 \ast B) \ast x_0 = 0.0077$
  • $x_2 = (2 - x_1 \ast B) \ast x_1 = 0.0081073$
  • $x_3 = (2 - x_2 \ast B) \ast x_2 = 0.00813001$
  • $x_4 = (2 - x_3 \ast B) \ast x_3 = 0.00813008$

接下來進行移位到整數部分,進行乘法後再移回去即可。

實作細節

精準度

使用浮點數實作的 FFT,需要小心精準度,網路上有些代碼利用合角公式取代建表。這類型的優化,在精準度不需要很高的時候提供較好的效能,卻無法提供精準的值,誤差修正成了一道難題。

牛頓迭代

由於有 FFT 造成的誤差,檢查收斂必須更加地保守,我們可能需要保留小數點下 $n+10$ 位,在收斂時檢查 $n+5$ 位確保收斂。通常收斂可以在 20 次左右完成,若發現迭代次數過多,有可能是第一個精準度造成的影響,盡可能使用內建函數得到的三角函數值。

加速優化

一般使用 IEEE-754 的浮點數格式,根據 FFT 的長度,要避開超過的震級,因此,在 Zerojudge b960 中,最多使用十萬進制進行加速。

誤差修正

儘管算出的商 $C=A \ast (1/B)$,得到的 $C$ 可能會差個 1,需要在後續乘法中檢驗。

參考題目/資料

b960 的亂碼

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

template<typename T> class TOOL_FFT {
public:
struct complex {
T x, y;
complex(T x = 0, T y = 0):
x(x), y(y) {}
complex operator+(const complex &A) {
return complex(x+A.x,y+A.y);
}
complex operator-(const complex &A) {
return complex(x-A.x,y-A.y);
}
complex operator*(const complex &A) {
return complex(x*A.x-y*A.y,x*A.y+y*A.x);
}
};
T PI;
static const int MAXN = 1<<17;
complex p[2][MAXN];
int reversePos[MAXN];
TOOL_FFT() {
PI = acos(-1);
preprocessing();
}
int NumberOfBitsNeeded(int PowerOfTwo) {
for (int i = 0;; ++i) {
if (PowerOfTwo & (1 << i)) {
return i;
}
}
}
inline uint32_t FastReverseBits(uint32_t a, int NumBits) {
a = ( ( a & 0x55555555U ) << 1 ) | ( ( a & 0xAAAAAAAAU ) >> 1 ) ;
a = ( ( a & 0x33333333U ) << 2 ) | ( ( a & 0xCCCCCCCCU ) >> 2 ) ;
a = ( ( a & 0x0F0F0F0FU ) << 4 ) | ( ( a & 0xF0F0F0F0U ) >> 4 ) ;
a = ( ( a & 0x00FF00FFU ) << 8 ) | ( ( a & 0xFF00FF00U ) >> 8 ) ;
a = ( ( a & 0x0000FFFFU ) << 16 ) | ( ( a & 0xFFFF0000U ) >> 16 ) ;
return a >> (32 - NumBits);
}

void FFT(bool InverseTransform, complex In[], complex Out[], int n) {
// simultaneous data copy and bit-reversal ordering into outputs
int NumSamples = n;
int NumBits = NumberOfBitsNeeded(NumSamples);
for (int i = 0; i < NumSamples; ++i)
Out[reversePos[i]] = In[i];

// the FFT process
for (int i = 1; i <= NumBits; i++) {
int BlockSize = 1<<i, BlockEnd = BlockSize>>1, BlockCnt = NumSamples/BlockSize;
for (int j = 0; j < NumSamples; j += BlockSize) {
complex *t = p[InverseTransform];
int k = 0;
#define UNLOOP_SIZE 8
for (; k+UNLOOP_SIZE < BlockEnd; ) {
#define UNLOOP { \
complex a = (*t) * Out[k+j+BlockEnd]; \
Out[k+j+BlockEnd] = Out[k+j] - a; \
Out[k+j] = Out[k+j] + a;\
k++, t += BlockCnt;\
}

#define UNLOOP4 {UNLOOP UNLOOP UNLOOP UNLOOP;}
#define UNLOOP8 {UNLOOP4 UNLOOP4;}
UNLOOP8;
}
for (; k < BlockEnd;)
UNLOOP;
}
}
// normalize if inverse transform
if (InverseTransform) {
for (int i = 0; i < NumSamples; ++i) {
Out[i] = Out[i].x / NumSamples;
}
}
}
void convolution(T *a, T *b, int n, T *c) {
static complex s[MAXN], d1[MAXN], d2[MAXN], y[MAXN];
n = MAXN;
for (int i = 0; i < n; ++i)
s[i] = complex(a[i], 0);
FFT(false, s, d1, n);
s[0] = complex(b[0], 0);
for (int i = 1; i < n; ++i)
s[i] = complex(b[n - i], 0);
FFT(false, s, d2, n);
for (int i = 0; i < n; ++i)
y[i] = d1[i] * d2[i];
FFT(true, y, s, n);
for (int i = 0; i < n; ++i)
c[i] = s[i].x;
}
void preprocessing() {
int n = MAXN;
for (int i = 0; i < n; i ++) {
p[0][i] = complex(cos(2*i*PI / n), sin(2*i*PI / n));
p[1][i] = complex(p[0][i].x, -p[0][i].y);
}

int NumBits = NumberOfBitsNeeded(n);
for (int i = 0; i < n; i++)
reversePos[i] = FastReverseBits(i, NumBits);
}
};

TOOL_FFT<double> tool;
struct BigInt {
long long *v;
int size;
static const int DIGITS = 5;
static const int MAXN = 1<<17;
static int compare(const BigInt a, const BigInt b) {
for (int i = MAXN-1; i >= 0; i--) {
if (a.v[i] < b.v[i])
return -1;
if (a.v[i] > b.v[i])
return 1;
}
return 0;
}
void str2int(char *s, long long buf[]) {
int n = strlen(s);
size = (n+DIGITS-1) / DIGITS;
int cnt = n%DIGITS == 0 ? DIGITS : n%DIGITS;
int x = 0, pos = size-1;
v = buf;
for (int i = 0; i < n; i++) {
x = x*10 + s[i] - '0';
if (--cnt == 0) {
cnt = DIGITS;
v[pos--] = x, x = 0;
}
}
}
void println() {
printf("%lld", v[size-1]);
for (int pos = size-2; pos >= 0; pos--)
printf("%05lld", v[pos]);
puts("");
}
BigInt multiply(const BigInt &other, long long buf[]) const {
int m = MAXN;
static double na[MAXN], nb[MAXN];
static double tmp[MAXN];
memset(na+size, 0, sizeof(v[0])*(m-size));
for (int i = 0; i < size; i++)
na[i] = v[i];
memset(nb+1, 0, sizeof(v[0])*(m-other.size));
for (int i = 1, j = m-1; i < other.size; i++, j--)
nb[j] = other.v[i];
nb[0] = other.v[0];
tool.convolution(na, nb, m, tmp);
BigInt ret;
ret.size = m;
ret.v = buf;
for (int i = 0; i < m; i++)
buf[i] = (long long) (tmp[i] + 1.5e-1);
for (int i = 0; i < m; i++) {
if (buf[i] >= 100000)
buf[i+1] += buf[i]/100000, buf[i] %= 100000;
}
ret.reduce();
return ret;
}
void reduce() {
while (size > 1 && fabs(v[size-1]) < 5e-1)
size--;
}
BigInt divide(const BigInt &other, long long buf[]) const {
{
int cmp = compare(*this, other);
BigInt ret;
ret.size = MAXN-1, ret.v = buf;
if (cmp == 0) {
memset(buf, 0, sizeof(buf[0])*MAXN);
buf[0] = 1;
ret.reduce();
return ret;
} else if (cmp < 0) {
memset(buf, 0, sizeof(buf[0])*MAXN);
buf[0] = 0;
ret.reduce();
return ret;
}
}
// A / B = A * (1/B)
// x' = (2 - x * B) * x
int m = MAXN;
static double na[MAXN], nb[MAXN];
static double invB[MAXN], netB[MAXN], tmpB[MAXN];
static long long bufB[MAXN];
int PRECISION = size+10;
memset(nb+1, 0, sizeof(v[0])*(m-other.size));
for (int i = 1, j = m-1; i < other.size; i++, j--)
nb[j] = other.v[i];
nb[0] = other.v[0];

memset(invB, 0, sizeof(invB[0])*m);
{
long long t = 100000, a = other.v[other.size-1];
if (other.size-2 >= 0)
t = t*100000, a = a*100000+other.v[other.size-2];
for (int i = PRECISION-other.size; i >= 0; i--) {
invB[i] = t/a;
t = (t%a)*100000;
}
}
for (int it = 0; it < 100; it++) {
// netB = xi * B
tool.convolution(invB, nb, m, netB);
long long carry = 0;
for (int i = 0; i <= PRECISION*2; i++) {
bufB[i] = carry + (long long) (netB[i] + 1.5e-1);
if (bufB[i] >= 100000)
carry = bufB[i]/100000, bufB[i] %= 100000;
else
carry = 0;
bufB[i] = -bufB[i];
}
// tmpB = 2 - xi * B
bufB[PRECISION] += 2;
memset(tmpB, 0, sizeof(tmpB[0])*m);
for (int i = 0; i <= PRECISION*2; i++) {
if (bufB[i] < 0)
bufB[i] += 100000, bufB[i+1]--;
if (i != 0)
tmpB[m-i] = bufB[i];
else
tmpB[i] = bufB[i];
}
// netB = tmpB * invB
tool.convolution(invB, tmpB, m, netB);
{
long long carry = 0;
memset(bufB, 0, sizeof(bufB[0])*m);
for (int i = 0; i <= PRECISION*2; i++) {
bufB[i] = carry + (long long) (netB[i] + 1.5e-1);
if (bufB[i] >= 100000)
carry = bufB[i]/100000, bufB[i] %= 100000;
else
carry = 0;
}
}

{
int same = 1;
for (int i = PRECISION; same && i >= 5; i--)
same &= ((long long) (invB[i]) == bufB[i+PRECISION]);
if (same)
break;
}
memset(invB, 0, sizeof(invB[0])*m);
for (int i = 0; i+PRECISION <= PRECISION*2; i++)
invB[i] = bufB[i+PRECISION];
}

memset(na+1, 0, sizeof(v[0])*(m-size));
for (int i = 1, j = m-1; i < size; i++, j--)
na[j] = v[i];
na[0] = v[0];
tool.convolution(invB, na, m, netB);

BigInt ret;
ret.size = m-1;
ret.v = buf;
long long carry = 0;
for (int i = 0; i < m; i++) {
buf[i] = carry + (long long) (netB[i] + 1.5e-1);
if (buf[i] >= 100000)
carry = buf[i]/100000, buf[i] %= 100000;
else
carry = 0;
}
for (int i = 0; i+PRECISION < m; i++)
buf[i] = buf[i+PRECISION];
memset(buf+PRECISION, 0, sizeof(buf[0])*(m-PRECISION));
{
memset(na, 0, sizeof(na[0])*m);
for (int i = 1, j = m-1; i < m-1; i++, j--)
na[j] = buf[i];
na[0] = buf[0];
for (int i = 0; i < m; i++)
nb[i] = other.v[i];
tool.convolution(nb, na, m, netB);
long long carry = 0;
for (int i = 0; i < m; i++) {
bufB[i] = (carry + (long long) (netB[i] + 1e-1));
if (bufB[i] >= 100000)
carry = bufB[i]/100000, bufB[i] %= 100000;
else
carry = 0;
}
carry = 0;
for (int i = 0; i < m; i++) {
bufB[i] = v[i] - bufB[i] + carry;
if (bufB[i] < 0)
bufB[i] += 100000, carry = -1;
else
carry = 0;
}
BigInt R;
R.size = m-1, R.v = bufB;
if (compare(other, R) <= 0) {
buf[0]++;
for (int i = 0; buf[i] >= 100000; i++)
buf[i+1]++, buf[i] -= 100000;
}
}
ret.reduce();
return ret;
}
};

int main() {
static char sa[1<<20], sb[1<<20];
while (scanf("%s %s", sa, sb) == 2) {
static long long da[1<<19], db[1<<19], dc[1<<19];

BigInt a, b, c;
a.str2int(sa, da);
b.str2int(sb, db);
c = a.divide(b, dc);
c.println();
}
return 0;
}
Read More +

2017 Google Code Jam Round QR

「在 24 小時內想不到解法,就註定這一生為處男」-《解題暴君》

感謝小夥伴們、茵可、妮可的提醒和翻譯,提醒要比賽,結果題目意思猜了一個下午 …

A. Oversized Pancake Flipper

$N$ 個煎餅排成一列,每個煎餅有正反兩面,分別以 +- 表示,現在有一個超大的煎餅翻轉氣,一次可以翻轉恰好連續 $K$ 個煎餅的狀態,若一次不滿 $K$ 視為無效操作。給予起始盤面,請問至少要操作幾次才能使其全部都成為正面。

由邊界開始觀察,最左的正面一定不再被翻,如果翻了必然為多一次操作,不斷地推移下去,只有最左邊的反面要依序翻轉,貪心模擬即可。

B. Tidy Numbers

要找到一種特別的數字,這些數字滿足從高位到低位的每一位呈現非嚴格遞減,請問從數字 $1$ 判斷到 $N$ 最後一個滿足的數字為何。

最直觀的方案就 $N$ 開始倒退過去找第一個符合的解。為解決大測資,$N$ 大到不太可能窮舉,即便很容易出現 ???999999 的解,我們可以用更簡單的貪心算法找到這個小於等於 $N$ 的解。

  • 例如輸入為 423
  • 從最高位開始窮舉,第一位不能放 4,其原因在於貪心 444 已經大於 423,退而求其次 333 是可行的一組解
  • 接著再考慮下一位時,我們不必要求小於等於 2,因為前一次窮舉保證小於 423,直接從 9 開始嘗試 399 是否符合解 … 類推。

這樣的解法效能為 $O(m^2)$,由於位數很小,就可以偷懶地去寫它。

C. Bathroom Stalls

在公共澡堂中,有 $N$ 個位子排成一列,最左和最右有守衛。現在 $M$ 名使用者依照順序進入澡堂,每一個人會遠則離其他人盡可能遠的位置,意即往左和往右到第一個人的距離最小值最大化,往左距離為 $L_s$ 往右為 $R_s$,這樣的位置也許會有好幾個,這時候優先選擇 $\max(L_s, R_s)$ 最大的那一個 (也就是盡可能靠左側)。請問最後一個人的 $\max(L_s, R_s)$$\min(L_s, R_s)$ 分別為多少?

若使用純模擬,時間複雜度可能高達 $O(NM)$,大測資是無法在時間內通過的。仔細觀察到題目並未要求找到最後一個人的所在位置,在數學推敲上並不要求太高,那麼我們轉向繼續連續為空的區段分別有多少個。

  • 例如一開始長度為 $N=10$
  • 有一段長度為 10,第一個人必然會將這一個分成兩段 4 和 5
  • 下一個人則會把 5 再分成 2 和 2
  • 下一個人則會把 4 再分成 1 和 2

計算過程中,我們按照長度高到低依序挑出,便可以知道最終的結果。為了加速運算,需要紀錄某個長度有多少個,一次統計會使用相同挑選方案的人,來逼近最終的 $M$ 個人。時間複雜度 $O(\log N)$

D. Fashion Show

在一個方型時尚展場中,我們可以放置 3 種不同的服裝,分別為 +, x, o 三種,其中 +, x 可以為展場增添 1 分,o 可以增添 2 分,我們希望按照下述的規則使得展場總分最高:

  • 對於同行或者同列的放置服裝,任兩個之中,至少要有一個為 +
  • 對於斜對角 45 度的放置服裝,任兩個之中,至少要有一個為 x

現在給予一個空的展場和一些已經放置服裝的位置,你只能升級 +, x 服裝為 o 和增加新的服裝在空位置上,使得分數最高,輸出其中一組解即可。

如果只看行列,通常我們可以轉換成二分匹配,其中 (行 i) -- (列 j) 表明了一種放置方案,而中間的匹配邊成為最終的放置選擇。現在多了斜對角,匹配成為最難處理的部份,

題目的規則可以轉換成,行列最多放置一個非 +,同理對角最多放置一個非 x。按照原先的匹配思路,如果放置 x,則相當於匹配某行某列。同理,放置 +,則相當於匹配某斜和某反斜。這裡我們發現如果同時匹配,則視為 o,其匹配數恰好符合題目要求的得分。

對於已經擺設上去的點,預先將其匹配,並且那些行、列、斜線及反斜的代表點則不再與其他節點匹配,最後,我們構造節點個數為 $6N$ 的二分匹配圖,運行匈牙利或者是最大流算法都可解決。

Solution

A

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

int main() {
int testcase, cases = 0, K;
char s[1024];
scanf("%d", &testcase);
while (testcase--) {
scanf("%s %d", s, &K);
int n = strlen(s);
int ret = 0, valid = 0;
for (int i = 0; i < n; i++) {
if (s[i] == '+')
continue;
if (i+K <= n) {
ret++;
for (int j = 0; j < K; j++) {
if (s[i+j] == '-')
s[i+j] = '+';
else
s[i+j] = '-';
}
}
}
valid = 1;
for (int i = 0; i < n; i++)
valid &= s[i] == '+';
if (valid)
printf("Case #%d: %d\n", ++cases, ret);
else
printf("Case #%d: IMPOSSIBLE\n", ++cases);
}
return 0;
}

B

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

int main() {
int testcase, cases = 0;
uint64_t N;
scanf("%d", &testcase);
while (testcase--) {
char s[1024];
scanf("%s", s);
sscanf(s, "%llu", &N);
int n = strlen(s);
uint64_t ret = 0;
int begin_small = 0;
for (int i = 0; i < n; i++) {
int prev = ret%10;
int st = begin_small ? 9 : (s[i]-'0');
for (int j = st; j >= prev; j--) {
uint64_t test = ret;
for (int k = i; k < n; k++)
test = test*10 + j;
if (test <= N) {
if (!begin_small && j != s[i]-'0')
begin_small = 1;
ret = ret*10 + j;
break;
}
}
}
printf("Case #%d: %llu\n", ++cases, ret);
}
return 0;
}

C

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

int main() {
int testcase, cases = 0;
scanf("%d", &testcase);
while (testcase--) {
long long n, m;
scanf("%lld %lld", &n, &m);
map<long long, long long> S;
S[n] = 1;
long long ret1 = -1, ret2 = -1;
while (m) {
pair<long long, long long> e = *S.rbegin();
S.erase(e.first);
if (e.second >= m) {
ret1 = e.first/2, ret2 = (e.first-1)/2;
break;
}
m -= e.second;
if (e.first%2) {
long long t = e.first/2;
S[t] += 2LL*e.second;
} else {
long long t = e.first/2;
S[t] += e.second;
S[t-1] += e.second;
}
}
printf("Case #%d: %lld %lld\n", ++cases, ret1, ret2);
}
return 0;
}

D

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

const int MAXV = 40010;
const int MAXE = MAXV * 200 * 2;
const int INF = 1<<29;
typedef struct Edge {
int v, cap, flow;
Edge *next, *re;
} Edge;
class MaxFlow {
public:
Edge edge[MAXE], *adj[MAXV], *pre[MAXV], *arc[MAXV];
int e, n, level[MAXV], lvCnt[MAXV], Q[MAXV];
void Init(int x) {
n = x, e = 0;
for (int i = 0; i < n; ++i) adj[i] = NULL;
}
void Addedge(int x, int y, int flow){
edge[e].v = y, edge[e].cap = flow, edge[e].next = adj[x];
edge[e].re = &edge[e+1], adj[x] = &edge[e++];
edge[e].v = x, edge[e].cap = 0, edge[e].next = adj[y];
edge[e].re = &edge[e-1], adj[y] = &edge[e++];
}
void Bfs(int v){
int front = 0, rear = 0, r = 0, dis = 0;
for (int i = 0; i < n; ++i) level[i] = n, lvCnt[i] = 0;
level[v] = 0, ++lvCnt[0];
Q[rear++] = v;
while (front != rear){
if (front == r) ++dis, r = rear;
v = Q[front++];
for (Edge *i = adj[v]; i != NULL; i = i->next) {
int t = i->v;
if (level[t] == n) level[t] = dis, Q[rear++] = t, ++lvCnt[dis];
}
}
}
int Maxflow(int s, int t){
int ret = 0, i, j;
Bfs(t);
for (i = 0; i < n; ++i) pre[i] = NULL, arc[i] = adj[i];
for (i = 0; i < e; ++i) edge[i].flow = edge[i].cap;
i = s;
while (level[s] < n){
while (arc[i] && (level[i] != level[arc[i]->v]+1 || !arc[i]->flow))
arc[i] = arc[i]->next;
if (arc[i]){
j = arc[i]->v;
pre[j] = arc[i];
i = j;
if (i == t){
int update = INF;
for (Edge *p = pre[t]; p != NULL; p = pre[p->re->v])
if (update > p->flow) update = p->flow;
ret += update;
for (Edge *p = pre[t]; p != NULL; p = pre[p->re->v])
p->flow -= update, p->re->flow += update;
i = s;
}
}
else{
int depth = n-1;
for (Edge *p = adj[i]; p != NULL; p = p->next)
if (p->flow && depth > level[p->v]) depth = level[p->v];
if (--lvCnt[level[i]] == 0) return ret;
level[i] = depth+1;
++lvCnt[level[i]];
arc[i] = adj[i];
if (i != s) i = pre[i]->re->v;
}
}
return ret;
}
} g;

int main() {
int testcase, cases = 0;
scanf("%d", &testcase);
while (testcase--) {
int n, m;
scanf("%d %d", &n, &m);

g.Init(6*n+5);
static int mm[800][800];
static int mmPos[800][800][2];
int ret = 0;
memset(mm, 0, sizeof(mm));
int source = 6*n+1, sink = 6*n+2;
int board[128][128] = {};
int board_t[128][128] = {};
char board_s[128][128] = {};
int ban[800] = {};
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
mm[i][j+n] = 1;
mm[(i+j)+2*n][(i-j+n-1)+4*n] = 1;
mmPos[i][j+n][0] = i;
mmPos[i][j+n][1] = j;
mmPos[(i+j)+2*n][(i-j+n-1)+4*n][0] = i;
mmPos[(i+j)+2*n][(i-j+n-1)+4*n][1] = j;
}
}
for (int i = 0; i < m; i++) {
char s[16];
int x, y;
scanf("%s %d %d", s, &x, &y);
x--, y--;
board_s[x][y] = s[0];
if (s[0] == 'o') {
ret += 2, board[x][y] = 2;
ban[x] = 1, ban[y+n] = 1;
ban[(x+y)+(2*n)] = 1;
ban[(x-y+n-1)+4*n] = 1;
} else if (s[0] == '+') {
ret += 1, board[x][y] = 1;
ban[(x+y)+(2*n)] = 1;
ban[(x-y+n-1)+4*n] = 1;
// printf("%d %d\n", x+y+2*n, x-y+n-1+4*n);
} else if (s[0] == 'x') {
ret += 1, board[x][y] = 1;
ban[x] = 1, ban[y+n] = 1;
}
}
memcpy(board_t, board, sizeof(board));

for (int i = 0; i < n; i++)
g.Addedge(source, i, 1);
for (int i = 2*n; i < 4*n; i++)
g.Addedge(source, i, 1);
for (int i = n; i < 2*n; i++)
g.Addedge(i, sink, 1);
for (int i = 4*n; i < 6*n; i++)
g.Addedge(i, sink, 1);

for (int i = 0; i < n; i++) {
for (int j = n; j < 2*n; j++) {
if (ban[i] || ban[j] || !mm[i][j])
continue;
g.Addedge(i, j, 1);
}
}
for (int i = 2*n; i < 4*n-1; i++) {
for (int j = 4*n; j < 6*n-1; j++) {
if (ban[i] || ban[j] || !mm[i][j])
continue;
g.Addedge(i, j, 1);
}
}

int flow = g.Maxflow(source, sink);


for (int i = 0; i < n; i++) {
for (Edge *p = g.adj[i]; p != NULL; p = p->next) {
if (p->v >= n && p->v < 2*n && p->flow == 0) {
board_t[mmPos[i][p->v][0]][mmPos[i][p->v][1]]++;
// printf("---- %d %d\n", mmPos[i][p->v][0], mmPos[i][p->v][1]);
board_s[mmPos[i][p->v][0]][mmPos[i][p->v][1]] = 'x';
}
}
}
for (int i = 2*n; i < 4*n; i++) {
for (Edge *p = g.adj[i]; p != NULL; p = p->next) {
if (p->v >= 4*n && p->v < 6*n && p->flow == 0) {
board_t[mmPos[i][p->v][0]][mmPos[i][p->v][1]]++;
// printf("---- %d %d\n", mmPos[i][p->v][0], mmPos[i][p->v][1]);
board_s[mmPos[i][p->v][0]][mmPos[i][p->v][1]] = '+';
}
}
}

struct DD {
int x, y;
char val;
DD(int x, int y, char val): x(x), y(y), val(val) {}
};
vector<DD> V;
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
if (board_t[i][j] != board[i][j]) {
if (board_t[i][j] == 2)
V.push_back(DD(i+1, j+1, 'o'));
// printf("o %d %d\n", i+1, j+1);
else if (board_t[i][j] == 1)
V.push_back(DD(i+1, j+1, board_s[i][j]));
// printf("%c %d %d\n", board_s[i][j], i+1, j+1);
}
}
}
printf("Case #%d: %d %d\n", ++cases, ret+flow, V.size());
for (int i = 0; i < V.size(); i++)
printf("%c %d %d\n", V[i].val, V[i].x, V[i].y);
}
return 0;
}
Read More +

淺談旅行推銷員問題 (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 +

2017 Facebook Hacker Cup Round 1

感謝小夥伴妮可、茵可熱情支援

Facebook 2017 Hacker Cup Round 1

A. Pie Progress

單身狗的 $N$ 天日子中 (娛樂性質翻譯),每天晚餐想要一道點心派搭配。每天早晨決定到當地的餅舖採購,每天一定會生產 $M$ 派,每一種派的價格也有所不同,不用考慮派會壞掉的情況,預先庫存保留著吃。為防止不當商人購買數量過多,當天若購買 $K$ 個派,需要額外支付 $K^2$ 的交易手續費,請問採購花費最少為何。

明顯地,每一天的狀態就是採購了多少個派,得到狀態 $\text{dp}[i][j]$ 表示前 $i$ 天總共採購 $j$ 個餅,轉移過程中要保證數量足夠支付每一天,意即只對 $i \le j$ 進行轉移。每一天窮舉購買的數量,窮舉採購的花費時,勢必要先排序每塊派的價格,每次只挑選前幾個小的,時間複雜度 $O(N^2 M)$

$$\begin{align*} dp[i][j] = \left\{\begin{matrix} 0 && i = 0\;, j = 0\\ \min(dp[i-1][j-k-1] + \text{SumC}[k] + (k+1)^2) && i \le j \\ \infty && \text{otherwise} \end{matrix}\right. \end{align*}$$

B. Fighting the Zombies

在 D&D 遊戲,身為一個魔法師要消滅地圖上的殭屍們。一次操作有兩個步驟,第一步驟圈選任意半徑內的所有殭屍,不改變其相對位置將他們進行轉移,第二步驟選擇長寬為 $R$ 的方形內的所有殭屍,請問一次操作最多可以消滅多少殭屍。

從第二步驟中觀察到消滅的大小是固定的,因此圈選半徑會被約束在 $R$ 內,實際上也不用考慮圓,因為方形被包含在圓裏。最後,我們直接求第一步驟的所有方形情況,將內部的殭屍全部移除後,再窮舉一次方形範圍內部的其他殭屍,所有可能取最大值即可。時間複雜度 $O(N^6)$ 。由於 $N \le 50$ ,六分鐘內是可以容忍的。

C. Manic Moving

搬家公司在 $N$ 個城鎮之間服務,貨車司機打算用最小的油量花費依序完成公司給定 $K$ 個訂單。第 $i$ 名客戶要求從 $S_i$ 地搬到 $D_i$ ,貨車一次可以載運兩名客戶的量。根據訂單順序,先來的就要載貨,同理也要先卸貨。

從題目中發現對於順序要求非常嚴苛,定出每一階段的狀態 $dp[i][j][2]$ 表示完成前 $i$ 個訂單、最後停留位置在 $j$ 地,最後的 [2] 表示前一個階段是否已經卸貨。分成兩種方式討論,時間複雜度 $O(KN)$

D. Beach Umbrellas

$N$ 個人各自帶著半徑 $R_i$ 的降落傘,在海岸進行降落,岸上有 $M$ 個降落點,每個降落點間隔一公尺,請問有多少種降落方式使得他們不會碰撞。

從題目給的說明中,我們發現到左右兩側的降落點比較特別,因為他們的傘的一部份可能會落在 $M$ 點之外,因此考慮窮舉降落在左右側的所有方法數 $N^2$ ,若要計算固定左右兩側的方法數,可以使用重複組合 H 得到 (滿足 $x_1 + x_2 + \cdots + x_n = Y$ 且每個數皆為非負整數的方法數)。然而,這樣計算方法缺少順序,最後補上排列計數 $(N-2)!$

來講講窮舉左右兩側之後怎麼算出方法數

  • 左右兩側分別為 $R_i$$R_j$ 的情況
  • 海岸左寬度增加 $R_i$,同理右寬度增加 $R_j$
  • 如此一來,左右變數的情況就能套用重複組合分配 $N$ 個變數,總和為 $M + R_i + R_j$ ,每個變數至少大於等於 $R_i$

特別地,變數 $M$ 過大。在窮舉所有情況中,組合類型最多 $2R$ 種,而非 $N^2$ 種。計算量多到必須預先建表,每一個組合數 $C^{M+?}_{N}$ ,由於底數是固定的,利用區間滑動在 $O(1)$ 轉換 (需要乘法反元素支援)。預先建表的時間 $O(R)$,窮舉部分 $O(N^2 \log R)$

Solution A

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

int main() {
int testcase, cases = 0;
scanf("%d", &testcase);
while (testcase--) {
int N, M;
scanf("%d %d", &N, &M);
int dp[305][305] = {};
const int INF = 0x3f3f3f3f;
for (int i = 0; i <= N; i++) {
for (int j = 0; j <= N; j++)
dp[i][j] = INF;
}

dp[0][0] = 0;
for (int i = 0; i < N; i++) {
int A[305];
for (int j = 0; j < M; j++)
scanf("%d", &A[j]);
sort(A, A+M);
for (int j = 0, sum = 0; j < M; j++) {
sum += A[j];
A[j] = sum;
}
for (int j = i; j < N; j++) {
if (dp[i][j] == INF)
continue;
for (int k = 0; k < M && k+j <= N; k++) {
dp[i+1][j+k+1] = min(dp[i+1][j+k+1], dp[i][j] + A[k] + (k+1)*(k+1));
}
}
for (int j = i+1; j <= N; j++)
dp[i+1][j] = min(dp[i+1][j], dp[i][j]);
}
printf("Case #%d: %d\n", ++cases, dp[N][N]);
}
return 0;
}

Solution B

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

int main() {
int testcase, cases = 0;
scanf("%d", &testcase);
while (testcase--) {
int N, R;
scanf("%d %d", &N, &R);
vector< pair<int, int> > A;
set<int> SX, SY;
for (int i = 0; i < N; i++) {
int x, y;
scanf("%d %d", &x, &y);
A.push_back(make_pair(x, y));
SX.insert(x), SY.insert(y);
}
sort(A.begin(), A.end());

int ret = 0;
for (auto LX : SX) {
for (auto LY : SY) {
int cnt = 0;
vector<int> used(N, 0);
for (int i = 0; i < N; i++) {
if (A[i].first >= LX && A[i].first <= LX+R
&& A[i].second >= LY && A[i].second <= LY+R)
cnt++, used[i] = 1;
}
for (auto TX : SX) {
for (auto TY: SY) {
int dd = 0;
for (int i = 0; i < N; i++) {
if (used[i])
continue;
if (A[i].first >= TX && A[i].first <= TX+R
&& A[i].second >= TY && A[i].second <= TY+R)
dd++;
}
ret = max(ret, dd+cnt);
}
}
}
}
printf("Case #%d: %d\n", ++cases, ret);
}
return 0;
}

Solution C

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;

int main() {
int testcase, cases = 0;
scanf("%d", &testcase);
while (testcase--) {
int N, M, K;
scanf("%d %d %d", &N, &M, &K);

long long g[105][105] = {};
const long long INF = 1LL<<60;
for (int i = 1; i <= N; i++) {
for (int j = 1; j <= N; j++)
g[i][j] = INF;
g[i][i] = 0;
}
for (int i = 0; i < M; i++) {
int A, B;
long long G;
scanf("%d %d %lld", &A, &B, &G);
g[A][B] = min(g[A][B], G);
g[B][A] = min(g[B][A], G);
}

for (int k = 1; k <= N; k++) {
for (int i = 1; i <= N; i++) {
for (int j = 1; j <= N; j++)
g[i][j] = min(g[i][j], g[i][k]+g[k][j]);
}
}

// for (int i = 1; i <= N; i++) {
// for (int j = 1; j <= N; j++)
// printf("%lld ", g[i][j]);
// puts("");
// }

int S[5005], D[5005];
for (int i = 0; i < K; i++)
scanf("%d %d", &S[i], &D[i]);
static long long dp[5005][105][2] = {};
for (int i = 0; i <= K; i++) {
for (int j = 0; j <= N; j++)
dp[i][j][0] = INF, dp[i][j][1] = INF;;
}
dp[0][1][0] = 0;

for (int i = 0; i < K; i++) {
int s1 = S[i], d1 = D[i];
for (int j = 1; j <= N; j++) {
long long cc;
cc = g[j][s1]+g[s1][d1];
dp[i+1][d1][0] = min(dp[i+1][d1][0], dp[i][j][0] + cc);
cc = g[j][s1];
dp[i+1][s1][1] = min(dp[i+1][s1][1], dp[i][j][0] + cc);

if (dp[i][j][1] != INF && i > 0) {
int sP = S[i-1], dP = D[i-1];
cc = g[j][s1]+g[s1][dP];
dp[i+1][dP][1] = min(dp[i+1][dP][1], dp[i][j][1] + cc);
cc = g[j][s1]+g[s1][dP]+g[dP][d1];
dp[i+1][d1][0] = min(dp[i+1][d1][0], dp[i][j][1] + cc);
}
}

// for (int j = 1; j <= N; j++)
// printf("%lld ", dp[i+1][j][0]);
// puts("");
}

long long ret = -1;
for (int i = 1; i <= N; i++) {
if (dp[K][i][0] != INF)
ret = max(ret, dp[K][i][0]);
}
printf("Case #%d: %lld\n", ++cases, ret);
}
return 0;
}

Solution D

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;

const long long MOD = 1000000007LL;
void exgcd(long long x, long long y, long long &g,
long long &a, long long &b) {
if (y == 0)
g = x, a = 1, b = 0;
else
exgcd(y, x%y, g, b, a), b -= (x/y) * a;
}
long long inverse(long long x, long long p) {
long long g, b, r;
exgcd(x, p, g, r, b);
if (g < 0) r = -r;
return (r%p + p)%p;
}

int main() {
int testcase, cases = 0;
scanf("%d", &testcase);
while (testcase--) {
int N, M, R[2048], S = 0, mxR = 0;
scanf("%d %d", &N, &M);
for (int i = 0; i < N; i++)
scanf("%d", &R[i]), S += R[i], mxR = max(mxR, R[i]);
if (N == 1) {
printf("Case #%d: %lld\n", ++cases, M);
continue;
}
fprintf(stderr, "%d %d %d\n", N, S, M);
long long invNplus = 1;
map<long long, long long> C;
{
long long f = 1;
for (int i = 1; i <= N; i++)
f = (f * i)%MOD;
invNplus = inverse(f, MOD);

int l = 1, r = 1;
f = 1;
for (int i = M - 2*S; i <= M - 2*S + 2*mxR; i++) {
int tM = i+N-1;
if (tM < N)
continue;
int L = tM-N+1, R = tM;
if (r < L)
l = r = f = L;
while (l < L)
f = (f*inverse(l, MOD))%MOD, l++;
while (r < R)
r++, f = (f * r)%MOD;
C[tM] = (f * invNplus)%MOD;
// printf("C(%lld %d) = %lld, %lld\n", tM, N, C[tM], f);
}
}
long long ret = 0;
for (int i = 0; i < N; i++) {
for (int j = 0; j < N; j++) {
if (i == j)
continue;
int tM = M + R[i] + R[j] - 2*S;
if (tM+N-1 < N)
continue;
// printf("add C(%d %d)\n", tM+N-1, N);
ret += C[tM+N-1];
ret %= MOD;
}
}

long long f = 1;
for (int i = 1; i <= N-2; i++)
f = (f * i)%MOD;
ret = ret * f;
ret %= MOD;
assert(ret >= 0);
printf("Case #%d: %lld\n", ++cases, ret);
}
return 0;
}
Read More +

2017 Facebook Hacker Cup 資格賽

原本只是想推碩一學弟去寫,學弟邀著邀著,我這個老骨頭只好跟著寫

Facebook Hacker Cup 2017 Qualification Round

A. Progress Pie

給一個落在 $(0, 0)\; , (100, 100)$ 矩形內部的圓餅圖,並且從垂直十二點方向開始,順時針繞一圈 $P \%$ ,又額外給一座標,請問該點是黑色還是白色,並且保證任何一組詢問點,鄰近 $10^{-6}$ 都屬於相同顏色。

從最後一個條件來看,我們處理邊界條件的誤差是可以容忍的。由於輸入都是整數,完全在整數上操作的部分尚未想到,但我們可以透過內積外積得到詢問點是順時針轉了 $R, \; 0 \le R < 2 \pi$ ,只需要判斷 $R$ 是否小於等於 $P$ 即可。

十二點鐘的方向向量為 $\vec{v} = (0, 50)$ ,詢問點與圓心的向量為 $\vec{u} = (X-50, Y-50)$ ,計算這兩個向量的夾角 $\theta = \cos^{-1}(\frac{u \cdot v}{|u| |v|})$ ,這樣算出來的角度只會落在 $[0, \pi)$ ,接著透過外積決定順時針還是逆時針,補回來即可。

B. Lazy Loading

搬家公司的工人要搬運 $N$ 個重量不同的傢俱,主管要求每次搬運至少 50 磅,工人為了偷懶,每次只搬運一部份的傢俱,然而主管不會準確計算工人搬運的總重,只會問一次搬運的最大重量和個數,工人想藉由分配方法來增加工作天數,請問要怎麼符合需求達到最多搬運天數。

可想而知,我們只需要貪心計算即可,每次挑選最重的那一個,接著搭配當前最輕的來湊數,一超過 50 磅就當作一天的搬運方案,直到沒有物品。一開始排序好 $O(N \log N)$ ,接著只需要掃描一次 $O(N)$ 即可完成。

C. Fighting the Zombie

在 D&D 遊戲中,我們需要施放技能攻擊血量為 $H$ 的殭屍,施放採用擲骰子的方式,骰一個 $Y$ 面骰 $X$ 次得到的點數總和加上固定值 $Z$ ,請問一擊必殺的機率最高為何,由於盤面上有許多骰子可以挑選,請輸出機率最高的那個骰子的機率為何。

首先,我們必須先瞭解最基礎的六面骰,投擲 $X$ 總和的方法數怎麼計算,定義 $\text{dp}[i][j]$ 表示投擲 $i$ 次,點數總和為 $j$ 的方法數。我們得到

$$\begin{align*} dp[i][j] = \left\{\begin{matrix} 1 && i = 0\;, j = 0\\ dp[i-1][j-1] + dp[i-1][j-2] + \cdots + dp[i-1][j-6] && i \le j \\ 0 && \text{otherwise} \end{matrix}\right. \end{align*}$$

上述的遞迴考慮 $i-1$ 個骰子的總和方法數,再決定第 $i$ 個骰子擲出哪一種點數。然而,這種方法不適用此題計算機率,很容易發生 overflow,方法數的總和為 $Y^X$ ,所以一開始我們就採用機率的方式統計。

$$\begin{align*} dp[i][j] = \left\{\begin{matrix} 1 && i = 0\;, j = 0\\ (dp[i-1][j-1] + dp[i-1][j-2] + \cdots + dp[i-1][j-6])/6 && i \le j \\ 0 && \text{otherwise} \end{matrix}\right. \end{align*}$$

這樣的 DP 計算消耗時間 $O(X^2 Y^2)$ ,加上滑動窗口統計總和則可以落在 $O(X^2 Y)$

比賽當下寫的,思路不是很清楚,變數命名和邏輯判斷會有點醜,沒有好好整理。

Solution A

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

int main() {
const double eps = 1e-8;
const double pi = acos(-1);
int testcase, cases = 0;
scanf("%d", &testcase);
while (testcase--) {
int P, X, Y;
scanf("%d %d %d", &P, &X, &Y);
int ret = 0;

if (X == 50 && Y == 50) {
ret = 1;
} else if (hypot(X-50, Y-50) > 50) {
} else if (P == 100) {
ret = 1;
} else if (P == 0) {

} else {
int vx = X-50, vy = Y-50;
int tx = 0, ty = 50;
double theta = acos((vx*tx+vy*ty)/hypot(vx, vy)/hypot(tx, ty));
if (tx*vy - ty*vx > 0)
theta = 2*pi-theta;
double t = (double) P/100.0*2*pi;
if (theta <= t)
ret = 1;
}
printf("Case #%d: %s\n", ++cases, ret ? "black" : "white");
}
return 0;
}

Solution B

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

int main() {
int testcase, cases = 0;
scanf("%d", &testcase);
while (testcase--) {
int n;
vector<int> A;
scanf("%d", &n);
for (int i = 0; i < n; i++) {
int x;
scanf("%d", &x);
A.push_back(x);
}
sort(A.begin(), A.end());
int ret = 0;
int r = n-1, l = 0;
while (r >= l) {
int need = ((50 + A[r]-1) / A[r]);
if (l + need-1 > r)
break;
l += need-1, r--;
ret++;
}
printf("Case #%d: %d\n", ++cases, ret);
}
return 0;
}

Solution C

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

int main() {
int testcase, cases = 0;
scanf("%d", &testcase);
while (testcase--) {
int H, S;
scanf("%d %d", &H, &S);
vector< pair<int, int> > A;
double ret = 0;
for (int i = 0; i < S; i++) {
char s[128];
scanf("%s", s);
int X = 0, Y = 0, Z = 0;
for (int j = 0, x = 0, sign = 1, idx = 0; s[j]; j++) {
if (isdigit(s[j]))
x = x * 10 + s[j] - '0';
if (s[j+1] == '\0' || !isdigit(s[j])) {
x = x * sign;
if (idx == 0)
X = x;
else if (idx == 1)
Y = x;
else
Z = x;
if (s[j] == '-')
sign = -1;
else
sign = 1;
x = 0, idx++;
}
}
int l = X+Z, r = X*Y+Z;
static const int OFF = 1024;
static double dp[2][OFF];
for (int i = 0; i < OFF; i++)
dp[0][i] = dp[1][i] = 0;
dp[0][0] = 1;
for (int i = 0; i < X; i++) {
int p = i&1, q = i&1^1;
for (int j = 0; j < OFF; j++)
dp[q][j] = 0;
for (int j = 0; j < OFF; j++) {
if (dp[p][j] <= 0) continue;
for (int k = 1; k <= Y; k++)
dp[q][j+k] += dp[p][j]*((double) 1.f/Y);
}
}
double sum = 0;
for (int j = 0; j <= X*Y; j++) {
if (j+Z >= H)
sum += dp[(X-1)&1^1][j];
}
ret = max(ret, sum);
}
printf("Case #%d: %.6lf\n",++cases, ret);
}
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 +

UVa 13100 - Painting the Wall

Problem

給一張 $N \times M$ 的圖,用最少筆畫勾勒出指定圖像,每一筆只能平行於兩軸,並且輸出任意一組解。

Sample Input

1
2
3
4
5
6
5 7
.*...*.
.*...*.
.*****.
.*...*.
.*...*.

Sample Output

1
2
3
4
3
vline 2 1 5
vline 6 1 5
hline 3 2 6

Solution

明顯地,要繪製的座標 $(x, y)$,可以視為一條邊 $x \Rightarrow y$,只要挑選某些點就可以覆蓋相鄰邊即可,在二分圖上找最少點集覆蓋是 P 多項式時間內可解的,並且最少點集大小為最大匹配數。但這一題要求要連續的筆畫,因此需要針對給定的圖形進行重新編號,將不連續的 x 和 y 座標重新定義。

這裡採用 Dinic’s algorithm 解決二分匹配,最後利用貪心回溯的方式可以找到一組解。

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
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <algorithm>
#include <queue>
#include <stack>
#include <assert.h>
#include <set>
using namespace std;
const int MAXV = 40010;
const int MAXE = MAXV * 200 * 2;
const int INF = 1<<29;
typedef struct Edge {
int v, cap, flow;
Edge *next, *re;
} Edge;
class MaxFlow {
public:
Edge edge[MAXE], *adj[MAXV], *pre[MAXV], *arc[MAXV];
int e, n, level[MAXV], lvCnt[MAXV], Q[MAXV];
void Init(int x) {
n = x, e = 0;
assert(n < MAXV);
for (int i = 0; i < n; ++i) adj[i] = NULL;
}
void Addedge(int x, int y, int flow){
edge[e].v = y, edge[e].cap = flow, edge[e].next = adj[x];
edge[e].re = &edge[e+1], adj[x] = &edge[e++];
edge[e].v = x, edge[e].cap = 0, edge[e].next = adj[y];
edge[e].re = &edge[e-1], adj[y] = &edge[e++];
assert(e < MAXE);
}
void Bfs(int v){
int front = 0, rear = 0, r = 0, dis = 0;
for (int i = 0; i < n; ++i) level[i] = n, lvCnt[i] = 0;
level[v] = 0, ++lvCnt[0];
Q[rear++] = v;
while (front != rear){
if (front == r) ++dis, r = rear;
v = Q[front++];
for (Edge *i = adj[v]; i != NULL; i = i->next) {
int t = i->v;
if (level[t] == n) level[t] = dis, Q[rear++] = t, ++lvCnt[dis];
}
}
}
int Maxflow(int s, int t){
int ret = 0, i, j;
Bfs(t);
for (i = 0; i < n; ++i) pre[i] = NULL, arc[i] = adj[i];
for (i = 0; i < e; ++i) edge[i].flow = edge[i].cap;
i = s;
while (level[s] < n){
while (arc[i] && (level[i] != level[arc[i]->v]+1 || !arc[i]->flow))
arc[i] = arc[i]->next;
if (arc[i]){
j = arc[i]->v;
pre[j] = arc[i];
i = j;
if (i == t){
int update = INF;
for (Edge *p = pre[t]; p != NULL; p = pre[p->re->v])
if (update > p->flow) update = p->flow;
ret += update;
for (Edge *p = pre[t]; p != NULL; p = pre[p->re->v])
p->flow -= update, p->re->flow += update;
i = s;
}
}
else{
int depth = n-1;
for (Edge *p = adj[i]; p != NULL; p = p->next)
if (p->flow && depth > level[p->v]) depth = level[p->v];
if (--lvCnt[level[i]] == 0) return ret;
level[i] = depth+1;
++lvCnt[level[i]];
arc[i] = adj[i];
if (i != s) i = pre[i]->re->v;
}
}
return ret;
}
} g;
const int MAXN = 1024*1024;
int mx[MAXN], my[MAXN];
int X[MAXN], Y[MAXN];
int n, m, VX, VY;
int C[MAXN][3], R[MAXN][3];
void dfs(int u) {
if (X[u]) return ;
X[u] = 1;
for (Edge *p = g.adj[u]; p != NULL; p = p->next) {
if (p->v - VX >= 1 && p->v - VX <= VY) {
if (my[p->v - VX] != -1 && Y[p->v - VX] == 0) {
Y[p->v - VX] = 1;
dfs(my[p->v - VX]);
}
}
}
}
int main() {
static char s[1024][1024];
while (scanf("%d %d", &n, &m) == 2) {
for (int i = 1; i <= n; i++)
scanf("%s", s[i]+1);

static int xg[1024][1024] = {};
static int yg[1024][1024] = {};
memset(xg, 0, sizeof(xg));
memset(yg, 0, sizeof(yg));
VX = 0, VY = 0;
for (int i = 1; i <= n; i++) {
for (int j = 1; j <= m; j++) {
if (s[i][j] == '*') {
if (xg[i][j] == 0) {
xg[i][j] = ++VX;
int k = j;
while (k <= m && s[i][k] == '*')
xg[i][k] = VX, k++;
C[VX][0] = i, C[VX][1] = j, C[VX][2] = k-1;
}
if (yg[i][j] == 0) {
yg[i][j] = ++VY;
int k = i;
while (k <= n && s[k][j] == '*')
yg[k][j] = VY, k++;
R[VY][0] = j, R[VY][1] = i, R[VY][2] = k-1;
}
}
}
}

int source = 0, sink = VX+VY+1;
g.Init(VX+VY+2);
for (int i = 1; i <= VX; i++)
g.Addedge(source, i, 1);
for (int i = 1; i <= VY; i++)
g.Addedge(VX+i, sink, 1);
for (int i = 1; i <= n; i++) {
for (int j = 1; j <= m; j++) {
if (s[i][j] == '*') {
// printf("%d - %d\n", xg[i][j], yg[i][j]);
g.Addedge(xg[i][j], yg[i][j] + VX, 1);
}
}
}
int mxflow = g.Maxflow(source, sink);
memset(mx, -1, sizeof(mx));
memset(my, -1, sizeof(my));
memset(X, 0, sizeof(X));
memset(Y, 0, sizeof(Y));
int match = 0;
for (int i = 1; i <= VX; i++) {
for (Edge *p = g.adj[i]; p != NULL; p = p->next) {
int x = i, y = p->v, flow = p->flow;
if (flow == 0 && y - VX >= 1 && y - VX <= VY) {
// match x - (y - VX) in bipartite graph
int r = x, c = y - VX;
mx[r] = c;
my[c] = r;
match++;
}
}
}
assert(match == mxflow);

for (int i = 1; i <= VX; i++) {
if (mx[i] == -1)
dfs(i);
}

printf("%d\n", mxflow);
for (int i = 1; i <= VX; i++) {
if (!X[i] && mx[i] != -1) {
printf("hline %d %d %d\n", C[i][0], C[i][1], C[i][2]);
mxflow--;
}
}
for (int i = 1; i <= VY; i++) {
if (Y[i]) {
printf("vline %d %d %d\n", R[i][0], R[i][1], R[i][2]);
mxflow--;
}
}
assert(mxflow == 0);
}
return 0;
}
Read More +

UVa 13102 - Tobby Stones

Problem

$N$ 個石頭排成一列,一開始都是白色石頭,經過 $M$ 詢問,分別統計區間內的黑色、白色石頭個數。

操作分別有以下幾種:

1. 統計區間 [l, r] 的黑色、白色石頭個數。
2. 將區間 [l, r] 的石頭的排列反轉。
3. 將區間 [l, r] 的石頭的顏色反轉。
4. 將區間 [l, r] 的石頭的顏色全部改成 $c$

## Sample Input ##

1
2
3
4
5
6
7
8
9
10
10  7
0 0 9
3 0 4 0
0 0 4
1 0 9
0 5 9
2 5 9
0 3 9
100 1
0 0 50

Sample Output

1
2
3
4
5
0 10
5 0
5 0
0 7
0 51

Solution

用 Splay Tree 完成區間反轉,打個懶標記就可以完成,每個操作都可以在 $\mathcal{O}(\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
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
#include <bits/stdc++.h> 
using namespace std;
const int MAXN = 1048576;

class SPLAY_TREE { // Splay Tree
public:
static const int MAXN = 1048576;
struct Node {
static Node *EMPTY;
Node *ch[2], *fa;
char rev, set, inv, val;
int size;
int sumw, sumb;
Node() {
ch[0] = ch[1] = fa = NULL;
rev = set = inv = 0;
val = sumw = sumb = 0, size = 1;
}
bool is_root() {
return fa->ch[0] != this && fa->ch[1] != this;
}
void pushdown() {
if (rev) {
if (ch[0] != EMPTY) ch[0]->rev ^= 1;
if (ch[1] != EMPTY) ch[1]->rev ^= 1;
swap(ch[0], ch[1]);
rev ^= 1;
}
if (set) {
if (ch[0] != EMPTY) ch[0]->set = set, ch[0]->inv = 0;
if (ch[1] != EMPTY) ch[1]->set = set, ch[1]->inv = 0;
if (set == 1) // white
val = 0, sumw = sumw + sumb, sumb = 0;
else
val = 1, sumb = sumw + sumb, sumw = 0;
set = 0, inv = 0;
}
if (inv) {
if (ch[0] != EMPTY) {
if (ch[0]->set)
ch[0]->set = 3 - ch[0]->set;
else
ch[0]->inv ^= 1;
}
if (ch[1] != EMPTY) {
if (ch[1]->set)
ch[1]->set = 3 - ch[1]->set;
else
ch[1]->inv ^= 1;
}
swap(sumw, sumb), val = !val;
inv = 0;
}
}
void pushup() {
if (ch[0] != EMPTY) ch[0]->pushdown();
if (ch[1] != EMPTY) ch[1]->pushdown();
sumw = sumb = 0;
if (val == 0) sumw++;
else sumb++;
sumw += ch[0]->sumw + ch[1]->sumw;
sumb += ch[0]->sumb + ch[1]->sumb;
size = 1 + ch[0]->size + ch[1]->size;
}
} _mem[MAXN];

int bufIdx;
SPLAY_TREE::Node *root;
SPLAY_TREE() {
Node::EMPTY = &_mem[0];
Node::EMPTY->fa = Node::EMPTY->ch[0] = Node::EMPTY->ch[1] = Node::EMPTY;
Node::EMPTY->size = Node::EMPTY->val = 0;
bufIdx = 1;
}
void init() {
bufIdx = 1;
}
Node* newNode() {
Node *u = &_mem[bufIdx++];
*u = Node();
u->fa = u->ch[0] = u->ch[1] = Node::EMPTY;
return u;
}
void rotate(Node *x) {
Node *y;
int d;
y = x->fa, d = y->ch[1] == x ? 1 : 0;
x->ch[d^1]->fa = y, y->ch[d] = x->ch[d^1];
x->ch[d^1] = y;
if (!y->is_root())
y->fa->ch[y->fa->ch[1] == y] = x;
x->fa = y->fa, y->fa = x;
y->pushup();
}
void deal(Node *x) {
if (!x->is_root()) deal(x->fa);
x->pushdown();
}
Node* find_rt(Node *x) {
for (; x->fa != Node::EMPTY; x = x->fa);
return x;
}
void splay(Node *x, Node *below) {
Node *y, *z;
deal(x);
while (!x->is_root() && x->fa != below) {
y = x->fa, z = y->fa;
if (!y->is_root() && y->fa != below) {
if (y->ch[0] == x ^ z->ch[0] == y)
rotate(x);
else
rotate(y);
}
rotate(x);
}
x->pushup();
if (x->fa == Node::EMPTY) root = x;
}
Node* build(int l, int r, Node *p, char s[]) {
if (l > r) return Node::EMPTY;
int mid = (l + r)/2;
Node *t = newNode();
t->val = s[mid], t->fa = p;
if (t->val == 0) t->sumw ++;
else t->sumb ++;
t->ch[0] = build(l, mid-1, t, s);
t->ch[1] = build(mid+1, r, t, s);
t->pushup();
if (p == Node::EMPTY) root = t;
return t;
}
void debug(Node *u){
if (u == Node::EMPTY) return;
u->pushdown();
debug(u->ch[0]);
printf("%d", u->val);
debug(u->ch[1]);
}
Node* kthNode(int pos) {
Node *u = root;
for (int t; u != Node::EMPTY;) {
u->pushdown();
t = u->ch[0]->size;
if (t+1 == pos) return u;
if (t >= pos)
u = u->ch[0];
else
pos -= t+1, u = u->ch[1];
}
return Node::EMPTY;
}
void insert(int pos, int val) {
Node *p, *q, *r;
p = kthNode(pos), q = kthNode(pos+1);
r = newNode();
splay(p, Node::EMPTY);
splay(q, root);
r->val = val, q->ch[0] = r, r->fa = q;
splay(r, Node::EMPTY);
}
void erase(int pos) {
Node *p, *q;
p = kthNode(pos-1), q = kthNode(pos+1);
splay(p, Node::EMPTY), splay(q, root);
q->ch[0] = Node::EMPTY;
splay(q, Node::EMPTY);
}
void reverse(int l, int r) {
Node *p, *q;
p = kthNode(l), q = kthNode(r+2);
splay(p, Node::EMPTY), splay(q, root);
q->ch[0]->rev ^= 1;
splay(q->ch[0], Node::EMPTY);
}
void inverse(int l, int r) {
Node *p, *q;
p = kthNode(l), q = kthNode(r+2);
splay(p, Node::EMPTY), splay(q, root);
q->ch[0]->inv ^= 1;
splay(q->ch[0], Node::EMPTY);
}
void reset(int l, int r, int c) {
Node *p, *q;
p = kthNode(l), q = kthNode(r+2);
splay(p, Node::EMPTY), splay(q, root);
if (c == 1) {
q->ch[0]->set = 1;
} else {
q->ch[0]->set = 2;
}
splay(q->ch[0], Node::EMPTY);
}
pair<int, int> stat(int l, int r) {
Node *p, *q;
p = kthNode(l), q = kthNode(r+2);
splay(p, Node::EMPTY), splay(q, root);
return make_pair(q->ch[0]->sumb, q->ch[0]->sumw);
}

} tree;
SPLAY_TREE::Node *SPLAY_TREE::Node::EMPTY;

int main() {
static char s[1048576] = {}; // white
int n, m, cmd, l, r, c;
while (scanf("%d %d", &n, &m) == 2) {
for (int i = 0; i <= n+1; i++)
s[i] = 0;
tree.init();
tree.build(0, n+1, SPLAY_TREE::Node::EMPTY, s);
pair<int, int> ret;
for (int i = 0; i < m; i++) {
scanf("%d", &cmd);
if (cmd == 0) {
scanf("%d %d", &l, &r);
ret = tree.stat(l+1, r+1);
printf("%d %d\n", ret.first, ret.second);
} else if (cmd == 1) {
scanf("%d %d", &l, &r);
tree.reverse(l+1, r+1);
} else if (cmd == 2) {
scanf("%d %d", &l, &r);
tree.inverse(l+1, r+1);
} else if (cmd == 3) {
scanf("%d %d %d", &l, &r, &c);
tree.reset(l+1, r+1, c);
}
}
}
return 0;
}
Read More +