批改娘 20021. Dynamic Range Sum

題目描述

在二維平面上有 $N$ 個點,每一個點 $p_i(x_i, y_i)$各自帶權重 $w_i$,這些點不時會移動和改變權重。

現在 Morris 希望你幫忙撰寫線性搜索的函數,好讓他專注數據結構上的調整。函數詢問包含

  • $N$ 個點的資訊 (以 SoA (Structure of Array) 的方式儲存,以達到最好的快取使用率)
  • 詢問的矩形 $\text{Rect}$ (正交於兩軸)

函數必須回傳在矩形內部的點權重和。

1
int32_t search_range(Rect rect, int32_t x[], int32_t y[], int32_t w[], int32_t n);

main.c (測試用)

  • 一開始,在二維空間 $[0, R) \times [0, R)$ 之間產生 $N$ 個點的資訊
  • 接著,模擬 $M$ 次點的變化,並且呼叫 search_range
  • 最後,將所有答案 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
#include <stdio.h>
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <inttypes.h>
#include "DRS.h"
static uint32_t seed = 0;
static void p_srand(uint32_t x) { seed = x;}
static uint32_t p_rand() {return seed = (seed*9301 + 49297);}
static void swap(int *x, int *y) {
int tmp = *x;
*x = *y;
*y = tmp;
}
static inline Rect rand_rect(int R) {
Rect r;
r.lx = p_rand()%R;
r.ly = p_rand()%R;
r.rx = p_rand()%R;
r.ry = p_rand()%R;
if (r.lx > r.rx) swap(&r.lx, &r.rx);
if (r.ly > r.ry) swap(&r.ly, &r.ry);
return r;
}
static void init(int N, int R, int32_t x[], int32_t y[], int32_t w[]) {
for (int i = 0; i < N; i++) {
x[i] = p_rand()%R;
y[i] = p_rand()%R;
w[i] = p_rand()%R;
}
}
static void tick(int N, int R, int32_t x[], int32_t y[], int32_t w[]) {
for (int i = 0; i < 5; i++) {
int idx = p_rand()%N;
x[idx] = p_rand()%R;
y[idx] = p_rand()%R;
}
for (int i = 0; i < 5; i++) {
int idx = p_rand()%N;
w[idx] = p_rand()%R;
}
}
#define MAXN 1048576
int main() {
p_srand(0);
static int32_t x[MAXN], y[MAXN], w[MAXN];
int N = 1000, M = 10000, R = 100;
init(N, R, x, y, w);
int32_t hash = 0;
for (int it = 0; it < M; it++) {
Rect rect = rand_rect(R);
int32_t ret = search_range(rect, x, y, w, N);
hash ^= ret;
tick(N, R, x, y, w);
}
printf("%" PRIi32 "\n", hash);
return 0;
}

DRS.h

1
2
3
4
5
6
7
8
9
10
11
12
13
#ifndef __DRS_H
#define __DRS_H
#include <stdint.h>
typedef struct Rect {
int32_t lx, ly, rx, ry;
} Rect;
int32_t search_range(Rect rect, int32_t x[], int32_t y[],
int32_t w[], int32_t n);
#endif

DRS.c

你的目標是要加速下述函數的計算,通過最低要求加速 2 倍以上。

1
2
3
4
5
6
7
8
9
10
11
12
13
#include "DRS.h"
int32_t search_range(Rect rect, int32_t x[], int32_t y[],
int32_t w[], int32_t n) {
int32_t ret = 0;
for (int i = 0; i < n; i++) {
if (rect.lx <= x[i] && x[i] <= rect.rx &&
rect.ly <= y[i] && y[i] <= rect.ry) {
ret += w[i];
}
}
return ret;
}

測資限制

  • $N \le 131072$
  • $R \le 32768$

範例輸入

no input

範例輸出

1
8967

編譯參數

1
2
3
4
5
all: main
main: main.c
gcc -std=c99 -Ofast -march=native DRS.c -c -o DRS.o
gcc -std=c99 -Ofast -march=native main.c DRS.o -o main

Solution

一般的線性作法,透過 rect.lx <= x[i] && x[i] <= rect.rx && rect.ly <= y[i] && y[i] <= rect.ry 操作,翻譯成組合語言時,四次的 branch & jump,在管線處理上的效能易受到影響。為了使這種 branch 減少而不使管線處理的效能下降,通常會使用查表法來完成,藉由并行的數值計算,在最後一步才進行 load/store 完成指定區塊的指令。

在這一題中,唯有條件式皆成立才執行,由於分枝操作上只有一種情況,故查表法就不適用於此。我們仍可以并行數個矩形判斷,並且存在至少一個矩形成立再運行 load 加總所需的資料,將會給予效能上的大幅提升。

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
#include "DRS.h"
#include <x86intrin.h>
/*
int32_t search_range(Rect rect, int32_t x[], int32_t y[],
int32_t w[], int32_t n) {
int32_t ret = 0;
for (int i = 0; i < n; i++) {
if (rect.lx <= x[i] && x[i] <= rect.rx &&
rect.ly <= y[i] && y[i] <= rect.ry) {
ret += w[i];
}
}
return ret;
}
*/
int32_t search_range(Rect rect, int32_t x[], int32_t y[], int32_t w[], int32_t n) {
__m128i ret = _mm_set_epi32(0, 0, 0, 0);
rect.lx--, rect.ly--;
rect.rx++, rect.ry++;
__m128i lx = _mm_broadcastd_epi32(*((__m128i *) &rect.lx));
__m128i ly = _mm_broadcastd_epi32(*((__m128i *) &rect.ly));
__m128i rx = _mm_broadcastd_epi32(*((__m128i *) &rect.rx));
__m128i ry = _mm_broadcastd_epi32(*((__m128i *) &rect.ry));
__m128i zo = _mm_set_epi32(0, 0, 0, 0);
__m128i ic = _mm_set_epi32(3, 2, 1, 0);
for (int i = 0; i+4 <= n; i += 4) {
__m128i sx = _mm_load_si128((__m128i *) (x+i));
__m128i sy = _mm_load_si128((__m128i *) (y+i));
__m128i c1 = _mm_and_si128(_mm_cmplt_epi32(lx, sx), _mm_cmplt_epi32(sx, rx));
__m128i c2 = _mm_and_si128(_mm_cmplt_epi32(ly, sy), _mm_cmplt_epi32(sy, ry));
if (_mm_testz_si128(c1, c2) == 0) {
__m128i cc = _mm_and_si128(c1, c2);
__m128i vi = _mm_add_epi32(ic, _mm_set_epi32(i, i, i, i));
__m128i rs = _mm_mask_i32gather_epi32(zo, w+i, ic, cc, 4);
ret = _mm_add_epi32(ret, rs);
}
}
int32_t sum = 0;
for (int i = (n>>2)<<2; i < n; i++) {
if (rect.lx <= x[i] && x[i] <= rect.rx &&
rect.ly <= y[i] && y[i] <= rect.ry) {
sum += w[i];
}
}
static int32_t tmp[4] __attribute__ ((aligned (16)));
_mm_store_si128((__m128i*) &tmp[0], ret);
sum += tmp[0] + tmp[1] + tmp[2] + tmp[3];
return sum;
}
Read More +

批改娘 20020. Dot Product

題目描述

請嘗試使用 SIMD 技術 AVX/SSE/MMX 來加速以下的純數值計算。

main.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
#include <stdio.h>
#include <assert.h>
#include <inttypes.h>
#include <stdint.h>
static inline uint32_t rotate_left(uint32_t x, uint32_t n) {
return (x << n) | (x >> (32-n));
}
static inline uint32_t encrypt(uint32_t m, uint32_t key) {
return (rotate_left(m, key&31) + key)^key;
}
static uint32_t f(int N, int off, uint32_t key1, uint32_t key2) {
uint32_t sum = 0;
for (int i = 0, j = off; i < N; i++, j++)
sum += encrypt(j, key1) * encrypt(j, key2), i++, j++;
return sum;
}
int main() {
int N;
uint32_t key1, key2;
while (scanf("%d %" PRIu32 " %" PRIu32, &N, &key1, &key2) == 3) {
uint32_t sum = f(N, 0, key1, key2);
printf("%" PRIu32 "\n", sum);
}
return 0;
}

輸入格式

有多組測資,每組一行包含三個整數 $N, \; \text{key1}, \; \text{key2}$,表示向量長度 $N$、向量 $\vec{A}$ 由亂數種子 $\text{key1}$ 產生、向量 $\vec{B}$ 由亂數種子 $\text{key2}$ 產生。

  • $1 \le N \le 16777216$

輸出格式

對於每組測資輸出一行整數,為 $\vec{A} \cdot \vec{B}$ 的 unsigned 32-bit integer 結果。

範例輸入

1
2
16777216 1 2
16777216 3 5

範例輸出

1
2
2885681152
2147483648

編譯參數

1
gcc -std=c99 -O3 -march=native main.c -lm

參考資料

Solution

對於數值計算時,SIMD 能充分地加速程序,每一個元素皆經過一連串的數學函數計算,那麼把一連串的數學式拆分,化成最簡的邏輯計算,並且找出常數向量存放到暫存器中。如在影像處理的程序,即使沒有 GPU 幫忙,搭配 SIMD 也是個不錯的選擇,可以加速 2~4 倍之多。

為了凸顯效能差異,題目設計時必須在計算單一元素結果上複雜些,防止大部分的加速效果是來自於減少 branch 操作 (loop unrolling)。

AVX

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
#include <stdio.h>
#include <assert.h>
#include <inttypes.h>
#include <stdint.h>
#include <x86intrin.h>
static inline uint32_t rotate_left(uint32_t x, uint32_t n) {
return (x << n) | (x >> (32-n));
}
static inline uint32_t encrypt(uint32_t m, uint32_t key) {
return (rotate_left(m, key&31) + key)^key;
}
static uint32_t SSE(int N, int off, uint32_t key1, uint32_t key2) {
uint32_t sum = 0;
for (int i = (N>>3)<<3; i < N; i++)
sum += encrypt(i+off, key1) * encrypt(i+off, key2);
__m256i s_i = _mm256_set_epi32(off, off+1, off+2, off+3, off+4, off+5, off+6, off+7);
__m256i s_4 = _mm256_set_epi32(8, 8, 8, 8, 8, 8, 8, 8);
__m256i s_k1 = _mm256_set_epi32(key1, key1, key1, key1, key1, key1, key1, key1);
__m256i s_k2 = _mm256_set_epi32(key2, key2, key2, key2, key2, key2, key2, key2);
uint32_t modk1 = key1&31;
uint32_t modk2 = key2&31;
uint32_t cmodk1 = 32 - modk1;
uint32_t cmodk2 = 32 - modk2;
__m256i s_ret = _mm256_set_epi32(0, 0, 0, 0, 0, 0, 0, 0);
N >>= 3;
for (int it = 0; it < N; it++) {
__m256i r_1 = _mm256_or_si256(_mm256_slli_epi32(s_i, modk1), _mm256_srli_epi32(s_i, cmodk1));
r_1 = _mm256_xor_si256(_mm256_add_epi32(r_1, s_k1), s_k1);
__m256i r_2 = _mm256_or_si256(_mm256_slli_epi32(s_i, modk2), _mm256_srli_epi32(s_i, cmodk2));
r_2 = _mm256_xor_si256(_mm256_add_epi32(r_2, s_k2), s_k2);
__m256i r_m = _mm256_mullo_epi32(r_1, r_2);
s_ret = _mm256_add_epi32(s_ret, r_m);
s_i = _mm256_add_epi32(s_i, s_4);
}
{
static int32_t tmp[8] __attribute__ ((aligned (32)));
_mm256_store_si256((__m256i*) &tmp[0], s_ret);
sum += tmp[0] + tmp[1] + tmp[2] + tmp[3];
sum += tmp[4] + tmp[5] + tmp[6] + tmp[7];
}
return sum;
}
int main() {
int N;
uint32_t key1, key2;
while (scanf("%d %" PRIu32 " %" PRIu32, &N, &key1, &key2) == 3) {
uint32_t sum = SSE(N, 0, key1, key2);
printf("%" PRIu32 "\n", sum);
}
return 0;
}

SSE

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
#include <stdio.h>
#include <assert.h>
#include <inttypes.h>
#include <stdint.h>
#include <x86intrin.h>
static inline uint32_t rotate_left(uint32_t x, uint32_t n) {
return (x << n) | (x >> (32-n));
}
static inline uint32_t encrypt(uint32_t m, uint32_t key) {
return (rotate_left(m, key&31) + key)^key;
}
static uint32_t SSE(int N, int off, uint32_t key1, uint32_t key2) {
uint32_t sum = 0;
for (int i = N/4*4; i < N; i++)
sum += encrypt(i+off, key1) * encrypt(i+off, key2);
__m128i s_i = _mm_set_epi32(off, off+1, off+2, off+3);
__m128i s_4 = _mm_set_epi32(4, 4, 4, 4);
__m128i s_k1 = _mm_set_epi32(key1, key1, key1, key1);
__m128i s_k2 = _mm_set_epi32(key2, key2, key2, key2);
uint32_t modk1 = key1&31;
uint32_t modk2 = key2&31;
uint32_t cmodk1 = 32 - modk1;
uint32_t cmodk2 = 32 - modk2;
__m128i s_ret = _mm_set_epi32(0, 0, 0, 0);
N >>= 2;
for (int it = 0; it < N; it++) {
__m128i r_1 = _mm_or_si128(_mm_slli_epi32(s_i, modk1), _mm_srli_epi32(s_i, cmodk1));
r_1 = _mm_xor_si128(_mm_add_epi32(r_1, s_k1), s_k1);
__m128i r_2 = _mm_or_si128(_mm_slli_epi32(s_i, modk2), _mm_srli_epi32(s_i, cmodk2));
r_2 = _mm_xor_si128(_mm_add_epi32(r_2, s_k2), s_k2);
__m128i r_m = _mm_mullo_epi32(r_1, r_2);
s_ret = _mm_add_epi32(s_ret, r_m);
s_i = _mm_add_epi32(s_i, s_4);
}
{
static int32_t tmp[4] __attribute__ ((aligned (16)));
_mm_store_si128((__m128i*) &tmp[0], s_ret);
sum += tmp[0] + tmp[1] + tmp[2] + tmp[3];
}
return sum;
}
int main() {
int N;
uint32_t key1, key2;
while (scanf("%d %" PRIu32 " %" PRIu32, &N, &key1, &key2) == 3) {
uint32_t sum = SSE(N, 0, key1, key2);
printf("%" PRIu32 "\n", sum);
}
return 0;
}
Read More +

批改娘 20018. Square Root of Vector Elements

Background

Intel 公司在 X86 指令集上提供了 AVX (Advanced Vector Extensions)、SSE (Streaming SIMD Extensions) 的特殊指令,這些在 SIMD 指令上提供強力的加速。

Problem

現在給你長度為 $N$ 的 32-bit floating point,分別將其開根號回傳。

main.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
#include <stdio.h>
#include <memory.h>
#include <math.h>
#include <time.h>
#include "VSQRT.h"
static float a[1048576], b[1048576];
int main() {
int N = 1048576;
for (int i = 0; i < N; ++i)
a[i] = i;
for (int it = 0; it < 20; it++)
{
memcpy(b, a, sizeof(a[0])*N);
clock_t t = clock();
for (int i = 0; i < 10; i++)
sqrt2(b, b+N);
t = clock() - t;
fprintf(stderr, "It took me %f seconds.\n", ((float) t)/CLOCKS_PER_SEC);
float sum = 0;
for (int i = 0; i < N; i++)
sum += b[i];
printf("%f\n", sum);
}
return 0;
}

VSQRT.h

1
2
3
4
5
6
#ifndef __VSQRT_H
#define __VSQRT_H
void sqrt2(float *begin, float *end);
#endif

VSQRT.c

請嘗試加速以下程式碼。

1
2
3
4
5
6
7
8
#include "VSQRT.h"
#include <math.h>
void sqrt2(float *begin, float *end) {
for (; begin != end; begin++)
*begin = sqrt(*begin);
}

Sample Input (stdin)

no input

Sample Output (stdout)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
1051776.125000
1051776.125000
1051776.125000
1051776.125000
1051776.125000
1051776.125000
1051776.125000
1051776.125000
1051776.125000
1051776.125000
1051776.125000
1051776.125000
1051776.125000
1051776.125000
1051776.125000
1051776.125000
1051776.125000
1051776.125000
1051776.125000
1051776.125000

Solution

這裡我們使用 Intel CPU 撰寫 SIMD 指令時,透過編譯參數 -msse -mavx ... 指定相關的指令集進行調適。這些單一指令多資料處理的操作,使得在指定層級得到較高的平行度。

為了使用它們,我們需要較多的前置處理,例如要把資料載入和儲存時,必然要符合對齊的標準,若記憶體位址不是 32 的倍數 (AVX 的 256 bits) 或者 16 的倍數 (SSE 的 128 bits),則系統會產生 trap,進而導致程式運行中斷而發生錯誤。這方面必須特別小心,有時編譯器恰好幫你對齊而正常運行,若增加其他變數時,對齊就有可能跑掉,這時候再去找錯誤會變得相當瑣碎。

通常些 SSE/AVX 加速指令,將會在 gcc -Ofast 中使用,而在 -O3 下仍有可能尚未開啟。而在 LLVM 中的 clang -O2 編譯下就會自動開啟 SSE。無論如何,要明白開啟這些的并行指令的代價,數據小造成運行效能低落,隨著資料量遽增才會明顯加速,若牽涉到浮點數計算,務必注意精準度的需求。

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
#include "VSQRT.h"
#include <stdint.h>
#include <math.h>
#include <x86intrin.h>
//#undef __AVX__
//#undef __SSE__
#if defined(__AVX__)
void sqrt2(float *a, float *end) {
uint32_t offset = ((void *) a - (void *) 0);
for (; a < end && (offset&31); a++, offset += sizeof(float))
*a = sqrt(*a);
__m256* ptr = (__m256*) a;
for (; a + 8 <= end; a += 8, ptr++)
_mm256_store_ps(a, _mm256_sqrt_ps(*ptr));
for (; a < end; a++)
*a = sqrt(*a);
}
#elif defined(__SSE__)
void sqrt2(float *a, float *end) {
uint32_t offset = ((void *) a - (void *) 0);
for (; a < end && (offset&15); a++, offset += sizeof(float))
*a = sqrt(*a);
__m128* ptr = (__m128*) a;
for (; a + 4 <= end; a += 4, ptr++)
_mm_store_ps(a, _mm_sqrt_ps(*ptr));
for (; a < end; a++)
*a = sqrt(*a);
}
#else
void sqrt2(float *begin, float *end) {
for (; begin != end; begin++)
*begin = sqrtf(*begin);
}
#endif
Read More +

助教生涯-平行總結

從碩班入學開始的這兩年內,助教一職圍繞在身邊許久,不擅長教別人的我卻非得接下這份工作-「教練,我可以不當助教嗎?」一轉眼兩年過去,當了三個學期半的助教,兩年間也不過四個學期,我想這些經歷肯定對於其他人是特別的,甚至會覺得有博班學長跟你一起當助教。不幸地,沒人可以給你經驗學習當助教-「從今天開始,你就是課程助教了!」

關於作業

大部分的課程都會使用的相同作業,助教只需要發下去給同學們寫就好。若加點變化,這就讓助教苦了,那是為什麼呢?一個小小的更動將更新參考答案、增加檢測方法、準備應付學生難以理解題目的諮詢,單兵作業的情況下,要如何做到全面防守便是一場長期消耗戰。

在學校課程中,也只有幾堂課是週週出作業、交作業,那週週改作業是怎麼一回事!那麼就是一個完全的管線處理,一下子從期初忙到期末,作業格式要找老師檢查、課堂中督導學生撰寫、課後收作業批閱、收到補交作業批改、接獲隔週作業需求設計、進行初步實驗和作業規劃、郵件通知學生 … 等,每週這樣子循環,有時候在想要是突然生病或消失,會出什麼事情呢?好想知道,希望是個可以找個人接替的日子,那麼就能好好地做自己的事情。

永無止盡地思考——這樣的設計會不會出事?誰能來提供更好的解法,能不能再設計得更好了一些?你是否也有相同的煩惱呢?

關於考試

考試不像作業有一陣子的緩衝期,出題面向與準備方向需要跟老師確定,一旦沒有喬好,就會在學生和老師間被踢來踢去,要是讓學生準備方向錯誤,那麼一定會件很可怕的事情。檢查再檢查,還是會有一些描述在當下造成學生無法理解的部分。更擔心的是,知道考題內容也沒有附上參考答案和配分方法。當問答題時更為頭痛,總會有那些出乎意料的答題方法,學生考試當下詢問時,這要怎麼辦呢?例外處理要怎麼給分?

考卷批改又是另一回事,由老師改還是助教改呢?改完就要準備讓學生檢討「助教?為什麼這會錯?為什麼他寫那樣就會對,我寫就不對?」有時候,文字描述總是很可怕的,當我們充分理解這位仁兄時,看考卷的觀點會變成「假設你懂」相反地就是「假設你不懂」,通常都要以後者為主,但如果偏向實作課程時,對於那些實作能力很好,卻不想背誦專有名詞定義的高手而感到惋惜。

最好的解決方法,就是助教都認得出來每個學生和表現,這樣也許會好上一些,一個學期能不能認識一個班級的人呢?有些很少出現和表現的同學,到底要怎麼認識啊!

關於環境架設

在程式撰寫方面,替學生進行環境設定相當苦惱,因為有分成標準與不標準,像我這種非正規教育出來的,只知道一堆網路搜索的資訊結果,哪一種解決方法才是最好的,能帶到下一個運行版本下持續運作,長遠性的規劃!只發現以前課程所遺留下來的資料,大多都有一些毛病在,還有一些尚未描述清楚的部分。「咱很笨,能教教我嗎?上上上屆學長們,你們消失到哪了?」

「原來我之前什麼都不明白啊」—《正解的卡多 KADO》

雜言

「無言」-《從零開始的魔法書》

當收到同學問程式,時常忘記給予任何嘗試的足跡,只有一份有錯的程序,接著就要助教幫忙通靈找錯誤,你先可以問問一起修課的同學嗎?每次都這樣子的話,咱已經受不了啊,一年、兩年過去,也疲於應付這些,我也想要屬於我自己的救贖啊!

「」-《為美好的世界獻上祝福》

每週跟課三個多小時,又要確定作業題目、生測資和測試,我自己也想鑽研作業,所追求的作業是沒有底的標準,能更好就更好,多給一點意見,讓我也能從中學習吧,為什麼都只有作業寫好就完事?有點傷心,不甘如此地結束,你能驚艷我的,對吧?

從系統環境架設、考試準備、通知學生和課程時替學生除錯,製作課程簡報,修正幾代學長遺留下來的錯誤,撰寫參考解答,追求更好的效能與算法,行政事項的準備,還有外來的委託架設。該經歷的都經歷了,只是一個人忙不太過來,也一直撐不住這麼多龐大的壓力。

然而,要準備離開學校的當下,信箱不時有一些待申請的大一新生信件,其中還會有外系學生想要使用我們課程中使用的系統來學習。系上的轉系考要我弄環境、考題,看來我的人生很不平靜呢 … 一些邏輯說不通的事情發生在身上時,過得是多麼無力呀。

Read More +

純粹 Part 2

接續上一篇所講,繼續說下去吧。

「神明曾向某個青年下了詛咒,那是看到悲傷之人就會全身感到痛苦的詛咒。青年為了避免自己痛苦,於是向所有悲傷的人伸出了援手 —— 神明著著他的樣子,仿製了一個外貌、動作與他一模一樣的假貨,他也像真正的青年一樣對所有悲傷的人伸出了援手。神明分別為他們起了名字 —— 一方為善,另一方為偽善。你覺得哪個是善,哪個是偽善?」-《重啟咲良田》

純粹究竟是什麼樣子?不帶有任何的自我意識,完全地無意識到周遭,更不會去評斷做了什麼,就只是去做,最終是好是壞那又是什麼?不是我們給予一個行為或結果的解釋嗎,都只是解釋的話,真相有必要探究嗎?若是純粹,那還可以被稱作人?還是被稱為機器人合適?

純粹地,想知道。

日常篇

原本只是想做一個簡單的平行算法,從趕 Conference 論文換成 Transaction 論文,原本 Conference 論文頁數不得超過十頁,導致很多提及的數據結構都不太想仔細說明,咱的指導老師只好先改其他人的論文投 Conference,而我被安排到投 Transaction,沒有任何時間壓力。

三月中是個分水嶺,同屆的只剩下我要寫論文,其他人已經不用到校 (只要不接助教工作,事實上不必總是到實驗室),而我還接下了平行助教的工作,到學校等著論文修改的 on call,平常實驗室只剩下我和學弟們,不時跟他們討些作業來做,亦或者討論計畫所需要的算法。也許有人會想問,那實驗室的博班學長呢?只有 meeting 那一天有機會在實驗室見著,平常是完全的零交流,是非常平行的實驗室呢!

當助教的日子裡,就像去年一樣「我不知道該怎麼辦,但覺得這樣下去不太行。」,若跟其他人一樣顧好自己畢業,其他事都不管的話,這裡的許許多多事情終究會重複上演。設計題目出測資,看相關論文後,做個簡單的簡報,在每週跟課的時候跟同學說明,為了設計新的題目而調整系統,按照這個行程走,每天都有一些工作需要處理。

「和其他人一起努力的話,很多問題可以迎刃而解」-《3月的獅子》

然而,做的一切還會被同學嫌。要說這裡是台大嗎?讓我這個半吊子當助教本身就是錯的,我是來學習的,不是來工作的!真正值得學習的對象卻很少出現在我眼前給予我指導,曾經的我也好想在這裡永遠追隨某人、想和其他人一起努力,但每個人都有自己的規則和約束,能一起努力的機會就更少了,有著各自的工作與理念,即使在第一學府眾多的人群中,也未必能找到能共事的同事。

助教工作時看著自己專研過的算法和代碼,不時地仍會看著其他人的程序學習,儘管在運行上可能不是很有效率、邏輯不漂亮。別因為追不到我的程序就放棄努力,我仍然可以向妳學習「就差那麼一步,你的想法將會推翻整個局面,本該如此的美麗,卻少了最後那份信心。」《3月的獅子》,別過於相信我,我並沒你們想像中地強大。

「本應如此地美麗」-《3月的獅子》

論文從十頁放寬到二十五頁的限制,能放置的圖表更加地豐富,老師要我從祖宗十八代的算法開始描述,必且要講得可以不查閱其他論文和資料就能閱讀,大部分的內容都被要求重寫,提及約十個算法、快二十張圖片說明,製圖耗費的經歷更加地多 (原因都來自於版本控制的需求,畫圖都像寫程式一般),調整描寫順序,用最簡單的方法講清楚一件事情,於是一下子就到了快二十五頁。改我論文的痛苦已經充分展現在老師的臉書塗鴉牆,那大概是老師有生之年第一次這麼辛苦。

「為了創作出傑作,我正在養精蓄銳」-《我的妹妹是黃漫老師》

別太期待會有什麼特別的,自己覺得挺簡單的,可能還沒有什麼用呢。

計畫篇

因為要回學校當助教又要改論文,暫時辭去了實習工作 (留職停薪?為什麼不打算讓我離職,怎麼想都很奇怪,這麼照顧我,真的辛苦你們了),畢業之後又會到哪裡,做些什麼事情,不經讓我想起《秒速五釐米》的那段話,

「總之這幾年裡,我很想向前邁進,想觸及那無法觸及的事物,儘管不知道那具體指什麼。我不知道這份勉強的感情是從何處如何孕育而生,只能一味地工作,等回過神來,日漸喪失彈性的心靈是如此傷痛。接著某天早上,我意識到過去那份銘刻於心的感情已消失得無影無踪,我知道自己到了極限,於是辭去了工作。」─《秒速五釐米》

那份心情到底要如何醞釀?為了什麼而活的心情?想等待著什麼樣的日子?太多太多的問題想要知道,還是不需要去想意義,若是每一件事情都去找意義,那就不能做事了,那先做看看吧,誰也不知道答案!當然,接受包養的情況也是可以接受的!

「請包養我!我真的不想去工作」-《不正經的魔術講師雨禁忌教典》

理念篇

在台大讀研究所,有時會覺得佔了某人的位置,畢竟一開始的我並不在這裡,而應該出現在這裡的人卻又不怎麼出現,我出現在這裡合適嗎?如果合適的話,為什麼要這樣對我?如果不合適的話,又為什麼要這樣要求我?不得不去思考這些問題,否則就很難邁出下一步。

「有時我會覺得佔了某人的位置,這讓我感到難過。比方說,活在世上就是一種浪費,霸佔了其他人的位置,他們也許能活出更有意義的人生」-《愛在來世》

很想追祢,也很想追妳,在同一個次元也罷,不同次元也行,那是成為一般人且獲得認同的證明,對於我是如此地渴望。

「就算是距離無法縮短,也不能成為我不前進的理由」-《3月的獅子》

人們說那種東西只是個雙保險,也許,只是害怕一個人待著,那晚的事情我還記得很清楚,晚上十點多寫完程式闔上筆電,跟家人說聲晚安,走上樓睡覺時,覺得今晚是如此地奇怪,總覺得從另一個視角看著自己,但什麼事情也沒發生,無疑地躺上床後,馬上被寫程式的疲累帶進夢裡。

半夜突然驚醒,正要拿起手機看時間,坐在床上都會倒下,在半夜裡無助地喊著,一邊思考著「是不是要結束人生了,我的腦已經要停止運作了嗎?」還好那時在家裡,半夜還能被送去急診,無力的我還得被推著輪椅進去,醫師只跟我說明是「腸胃炎」,那時的我不怎麼相信,畢竟腸胃炎都得過多少次了,都還沒遇過這種情況。

躺在上休息的我,吃不下也喝不下,就只是發呆著,等著可以分清楚方向,那時朋友仍用著我的帳號在臉書發圖文,家人還以為我已經好到可以用手機發文,誰都不知道、也沒發現早已陣亡的我,一如往常地回覆我。這才發現不需要我也可以正常運作,這世界的容錯能力很好,一瞬間就可以修正我這個錯誤。之後的我更沒想過會失去平衡一個多月,車也不敢騎,頭也不敢回,快速轉個頭可能就要跌倒,那究竟是什麼情況,至今還不清楚。

「我十七歲的時候還在妄想著『自己會不會成為很厲害的人啊』」-《3月的獅子》

曾經的幻想「自己會不會成為很厲害的人啊」,到了考大學的結果才發現一點也不行,有些東西就是學不來,真的學不來呀,怎麼樣追逐的速度都跟不上,只能在一旁一味地看著自己離大家越來越遠,最後吊了車尾。

「想到要為了賺錢,埋首研究如何讓日常生活更便利,就受不了,我要為了思考而思考」-《世紀天才》

我想要為了思考而思考,找了工作後,身旁的人就很在意薪水多高,不斷地吹捧那薪水價值,但我想那都不是重點,回顧一個常見的小故事——

〈說話的藝術〉的小故事「如果你說一個女大學生,晚上去夜總會陪酒,聽起來就不太好,可如果你說一個夜總會小姐,白天堅持去大學聽課,就滿滿的正能量了。如果你說你是一個學者,開了個公司,會被鄙視,認為你俗,真是斯文敗類。可是如果你說你是一個商人,經商之餘還專研學術,別人會肅然起敬,尊稱你為儒商。所以說話的時候,順序特別重要。」

我想說得差不多對照類似的意思,強調的順序反了!理念上感到非常奇怪的點。

貼了圖不久,就收到了訊息
「你可以為了思考而思考」
「然後錢給我,幫你花!」

此時,我的心情寫照

「你這麼說的話,不就是在表達想跟我永遠在一起嗎」-《從零開始的魔法書》

儘管知道這只是純粹的話,仍必須好好地建造高牆,否則容易受到動搖。

Read More +

純粹 Part 1

從二月到五月這段日子,發生了許多事情,也有些不變的日常。其中,有一段疲於應付的日子過得水生火熱,又在其中穿插了幾段小趣事,獻給想到聽八卦的朋友們,接下來將分成感情、日常、計畫與理念四篇。

感情篇——在一月時已提及了大部分感情篇的要素,那些沒有後續的故事,通常不太想去提及、不想去談論,那是純粹的還是虛偽的感情?這些問題隨著這段日子不斷反覆思考著。

日常篇——講著大部分的實驗室的情況,我的日常不是一般的日常,回家就是睡,醒來就會在實驗室,說是一般的娛樂,也都是在學校度過。其原因可以歸納在「反正我沒有朋友、沒有女朋友,有著大把大把的時間可以去思考。」

計畫篇——計畫著下一階段要怎麼運行,一旦失去學生身份,就得開始履行當兵義務,而當完兵又是怎麼樣的日子,我應該要以什麼的理念活下去?

理念篇——從國小曾思考著「世界是否從我眼界下才開始構築」,到現在不斷地加深自己的理念,與其拓展到個方位上的思考轉換,那些思考有什麼意義?

感情篇

明白和揣測對方的想法,就像寫程序就能想到執行的預期結果是怎麼樣子,只要構築環境與每一步的指令都相當明確,那麼結果必然呼之欲出。這個能力用於寫程式很方便,一旦用到了感情上就不被接受。

當朋友身邊有著不錯的對象,也有著像漫畫劇情般那樣發展,在對方尚未開口的處境下,喜歡與不喜歡就決定於下一步,有可能錯過一個表明心意的機會,人生就差了十萬八千里遠。基於上述,學弟和朋友們開始喧鬧「不試著找個機會告白嗎?」「她應該在等你了吧」… 等。

「其實我不善與人交談」-《風夏》

「別鬧了,你們不夠理解她,她本來就是這樣子的人 …」那些無法用語言表述的推測,明確地告訴我,對方本質加上輸入後的結果-不會成功、也不會失敗,沒有機會告白,也不會有機會告白,我試著將這些跟朋友說到,但是誰也不信我。

然後到了某天早晨,在意識不清的情況下,傳了一張圖片過去,上頭這麼寫著「我有著一種能力,喜歡的對象總會有男朋友」想著考慮要交男朋友的她,我的能力也許能幫上她,「妳交上男朋友了嗎?」如此地問道,收拾著包包去學校,在等十字路口紅綠燈時,看到走在一旁男女國中生們,才想起稍早些到底做了什麼蠢事 …

看著中午時,手機傳來了訊息通知,在開啟之前,內心感覺就像是考卷在發的過程中,在台下能明確地幻想到下五分鐘的未來,那種心境和感覺在此刻重現,「求你了,我所想的預測不要成真」忍了很久,開啟的第一句話大致上說明故事已經落幕,心裡早就知道的事情,真實呈現時的痛苦啊!

「雖然知道會這樣,不過還是很痛苦啊」-《清戀》

沒事兒,至少在情人節之前挑戰過了,失敗就可以給各位朋友一個交代,心中也鬆了一口氣

「我已經不要緊了」「連我的心靈支柱 … 我以為的心靈支柱都這麼說 …」

「失敗就意味著挑戰過了」-《3月的獅子》


那件事情也過了一陣子,開學前就要忙於論文報告和助教工作的準備,這段忙碌的日子過了好一陣子

「Hello, 我來了」在晚上九點多的實驗室門口傳來熟悉的聲音

「Hello, 妳怎麼突然來了?」還來不及反應過來,只能隨手抓個詞應答道

「剛 meeting 結束,過來看看你們」
「你還記得我的論文做些什麼嗎?」

「還記得些,我記得是運行 XXX YYYY 達到 XXXX」對於記憶有點信心,至少認為特別的或重要的事情不會忘

「沒錯!今天老師跟我說 …」

憑藉著演算法的那些小伎倆,在不擅長的領域到指出自己的想法,希望我還可以派上用場!在這一晚,我想在一旁的學弟又會開始幻想些什麼,儘管在失敗後告誡他們別再拿我的蠢事開玩笑。想著沒有信心地講完後,又要無奈地應對他們,這一晚的關卡好難過。

又過了一陣子,某個晚上

「Hello, 我來了」這種開場白,不知道為什麼總是這麼有趣
「你們實驗室有 XXX 可以借用一下嗎?急用」
「我們不是做那一領域的,不會有的,你有問過其他實驗室嗎?」
「沒耶,第一個想到的就是問你」

此時,我心中的寫照 《廢天使加百列》
『妳這樣子說,等等我要怎麼跟學弟解釋啊』這個誤會嚴重影響到擅長腦補的朋友們啊!

四處問了下相關領域的學長們是否能租借,最後只能安慰著「明天老實跟老師說吧,沒關係的」

又過了一陣子,某個下午

「Hi」又聽到熟悉的聲音,但下午寫程序到了一半,暫時不想停止手邊的工作。
「你們實驗室有 XXX 可以借用一下嗎?」

學弟起身,翻著好久沒有整理的紙箱仍沒有結果。工作告一段落的我,加入搜尋的行列
「看起來是沒有,就算有,那可以用嗎?」

「應該可以吧?」

繼續翻著箱子說道「看起來是沒有耶」

「ㄟㄟ,你在寫些什麼?」走到我的位子上,稍微看著螢幕上的程序說道

「那個 … 我在寫學弟的作業」身為學長寫著學弟的作業,完全感受不到威嚴

「我的論文重要?還是學弟作業重要!」

突然間,實驗室的人都在竊笑,想必此生的我已經結束了,這跳進黃河也洗不清。
《廢天使加百列》


四月初仍忙著論文撰寫,在這修改的煎熬日子裡,不是等待回應,就是將描述手法改來改去。然後,之前在遊戲《楓之谷》認識的交大某數已經提早把論文寫完等畢業,週末突然問我要不要去拜月老,一個剛跟女友分手的正常人和一個從未正常交過女友的肥宅就這麼約著網友聚會去龍山寺拜月老。

不知道那天是怎了,想把手上的事情放掉,鐵了心想做了不同的事情,但是關於「出門拜月老?」這個出遊議題,在內心掙扎著「出門就輸了,若去的話就驗證我二十多年來的失敗 …」最終,由於沒有任何不去的理由,隔天一早就出發了。

當日頂著大太陽,什麼事前調查都沒準備,那麼拜也沒什麼效用,想著當作觀光之旅來看待。至於月老求姻緣的部分,秉持著平常心看待囉!好久沒有拿著香在寺廟裡晃來晃去,有拜學業、拜健康、拜姻緣還有拜商業 … 等。想到自己要來拜個姻緣?嗯,這種無法明確努力的項目,就像天生註定的,我想行為只是基因的載體,意識到自己從未有過,再想著有過與沒有過的差別,我的確不是那一型的個性,那麼沒有也是很正常的吧。

然而,都來了一趟,仍硬是要我問問關於姻緣,心想求姻緣要連續三個聖筊,這個的機率這麼低,如果沒有對於對象任何要求的話,我有機會嗎?先來個二分搜尋法,

「五年內有機會遇的到嗎?」
「不會」神如此回應
「那 … 十年內遇的到嗎?」大膽突破點地詢問
「沒辦法」神明搖著頭說道
「二十年內呢?」作為人而言,我的想法太奢求了
「沒有」
「這一生有可能嗎?」人生要果斷點,也許才有機會
「沒有」

內心憔悴,先來個中場休息吧。若是再擲下去,探究真理的我也許會想問「還需要活下去嗎?神」看著一旁有人頂著大太陽折騰了好幾十分鐘「兄弟呀,我理解機率是多麼地困難」

「她跟我不同次元對吧?」換個說法看看吧,希望是找出來的
「對」
「這樣的我可以活得開心吧」淚水不斷地在內心流淌,原來神告訴我的真理是這個樣子,至今我都沒好好對待。
「可以」

Read More +

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 +

探討平行與優化技術於熱輻射法 (下)

接續上一篇

平行設計 Parallel Design

回顧原始算法,為了加快收斂速度,每一次會挑選單位面積能量最多的三角形 $f$ ,之後便拿所有三角形 $t$ 進行傳導,直到傳導能量小於閥值,迭代將停止。算法如下所述:

1
2
3
4
5
6
while (not converge) {
f = PickMaxRadiosityTriangle()
foreach triangle t in model
shader(f, t);
clearRadiosity(f)
}

在嘗試傳遞的過程中,若三角形 $f$ 的三頂點的能量差異大,則選擇自適應切割三角形,直到三頂點能量差異小,這麼計算 Form-Factor 才會正確。在自適應部份,切割方法如下圖所示:

當偵測到綠色三角形 $A$ 頂點之間的 Form-Factor 差異過大時,使用最長邊的中點切割,這麼做是盡可能產生銳角三角形,為了圖形完整,必然也要對鄰居切割。為減少計算量,只算新增中心點 $P$ 的 Form-Factor。對於下方的三角形 $B$ 而言,分成兩種情況,已在這一輪完成計算,則重新計算 Form-Factor;相反地,不做任何事。

Single-Source Parallel Algorithm

我們發現計算 Form-Factor 相當獨立的,但自適應處理需要遞迴切割,因此選用多核心平台,而非通用圖形處理器平台,因為目前的 GPU 實作遞迴所需的 stack 使用 global memory 作為存取位址,所以一般多核心平台效果會更好。在這一次報告中,我們選用 OpenMP 這套跨平台多執行緒 API 進行實驗。

著手將多個三角形平行處理,意即每一個執行緒負責多個三角形的 Form-Factor 計算。

1
2
3
4
5
6
// Single-Source Parallel Algorithm
while (not converage) {
f = PickMaxRadiosityTriangle()
parallel foreach triangle t in model
shader(s, t);
clearRadiosity(f)

Multi-Source Parallel Algorithm

從原始算法中,我們發現到每一次迭代將只有一個熱源輻射到場景中,當場景有多個高能量熱源時,場景必須經過好幾次迭代才能近似最終結果。藉由平行處理的效能,我們可以一次迭代多熱源,便可將低執行緒之間分配工作不均的情況,不僅僅前幾次迭代就能近似最終結果,同時也能加速運算。

1
2
3
4
5
6
7
8
9
10
// Multi-Source Parallel Algorithm
while (not converge) {
set<Triangle> f = RadiosityTriangleCandiateCandidate();
parallel foreach triangle t in model
if (f.find(t))
continue;
foreach s in f
shader(s, t);
clearRadiosity(f);
}

我們所用的平行方無法搭配上述自適應的切割方案,其原因在於分裂過程中,同時也要對鄰居三角形分裂,整個圖形產生的節點與邊的關係才會正確,無法保證鄰居在同一執行緒內處理,若沒有做好空間切割,我們便無法處理這部份。若模型格式會是數個獨立的物體,而非單純的三角形資訊,可分配每一個執行緒處理多個獨立物體,我們預期可以達到更好的效果。

平行過程中每一個執行緒共享和衝突的區段越少越好,這意味著我們必須在運行輻射前就必須將模型切得相當細緻。特別注意到,切得細緻與否對於光線投射 (Ray Casting) 複雜度不變,因為邏輯上他們處理同一平面。

一旦切得細緻,傳遞的效果就不是這麼好,在邊界的陰影更加顯著。根據理論和實作層面推測,其一原因是能傳遞的總能量隨著迭代減少,那麼從分裂過程中傳遞能量採用較多的加法完成,相較於多個 32-bit floating point 誤差就少了許多。

我們也試著使用獨立的切割方案─重心切割,切割的結果不依賴鄰居,只需要在加入三角形清單部份使用 critical section 即可,效能影響並不大。

根據重心切割,下述實驗中,從 156 個三角形,自動分裂到 30000 個三角形後進行輻射的結果如下圖所示:

明顯地,根據重心的切割方法容易產生鈍角三角形,看起來就會像很多紡錘體。在眾多數學性質中,只使用重心也許不是好的解決方案,這是値得探討的一部份,由於製作上的時間限制,我們並沒有去探討各個不同切割方案,所對應的自適應的效果如何。

Longest Edge Center

展示結果 Demo

只使用優化技術渲染結果

blocks room
hall church

平行效能比較

在 Intel Xeon E5-2620 v3 上,我們測試不同平行度帶來的影響,由於只有兩個實體 CPU,每一個 CPU 有 6 個核心,每個核心皆有 Hyper-threading 技術,故可產生 24 個執行緒。

Single-Source Parallel Algorithm - Scalability

我們對模型 room.tri 以預先切割 14977 個三角形後,根據先前提到的平行算法 Single-Source Parallel Algorithm,即是迭代一次只取一個熱源,平行計算所有三角形到此熱源的 Form-Factor 値,針對不同的執行緒個數和運行時間記錄,結果如上圖。在由於過多的執行緒可能會帶來更多的 false sharing,造成資料在不同的 CPU 之間運行 data transmission,所以效果就逐漸不明顯。

Multi-Source Parallel Algorithm - Scalability

接著,我們測試 Multi-Source Parallel Algorithm,以 room.tri 預先切割 50017 個三角形後,每次迭代皆取數個熱源,平行計算所有三角形到所有被選取熱源的 Form-Factor 的總和。同樣地,因為查找的資料重複存取的模式造成不好的影響,類似上述所提到的 false sharing 影響,故會呈現一種陡坡。

參考論文 Reference

後記 Note

當我們進行平行效能比較時,要特別小心編譯器行為的差異,意即平行處理 $P$ 控制時,當 $P=1$ 時,使用平行版本比較合適。因為有時候,平行部分的函數被編譯器偵測到,由於分享記憶體的關係,部分代碼無法優化,導致效能瞬間慢個四到五倍都是有可能的,再加上 false sharing 的關係,更有可能在發生密集計算時,效能更加低落。

在不同平台上的情況也有所不同,例如在一般 server 上運行時,CPU 頻率通常都會比一般 PC 慢上許多,又因為很多個 CPU 導致總共的快取大小遠比一般 PC 多,所以平行效能將會受限於運行的應用行為,是否需要時常存取大量的資料。

致謝 Thanks

一開始在挑選主題相當困惑,畢竟使用 Unity 可以做出更生動的作品,但作品的元素和創意相當重要,相當迷惘於要選哪一種類型才好。如果要兩人以上一起做,那麼主題又不能太過狹隘。百般思慮下,還是由我拉選了這個主題下來做。

「我可能會扯你後腿,還是我們各別做?」組員擔心地說道

不用擔心,我自己也沒信心將所有程序都看完且修改更好更快,每天都焦頭爛額地煩惱整份程序的運作,深怕來不及在時限內完成足夠的報告份量。

「快來幫幫我啊,在身邊一起 trace 程序也好,咱的記憶體不足啊。」內心如此吶喊

一起修課的博班學長給了我們一些意見與鼓勵,而學弟們只會在一旁扯後腿問「學姊今天會來嗎?」最後,我們完成了整份程序的理解與討論。

Read More +