UVa 12254 - Electricity Connection

Problem

給一個 $8 \times 8$ 平面圖,上面至多有八個住家,我們目標要從發電廠出發,牽電到所有的住家,拉線跨過水路的花費 $pw$、一般陸路為 $pl$,求最少花費。

經典的斯坦納樹問題,但有別於一般的平面圖,使用歐基里德距離或者曼哈頓距離作為花費函數。接著讓我們細談如何進行常數優化。

Sample Input

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
2
0 10
H.W.WH..
..W.W...
..WGW...
........
........
........
........
........
0 0
H.W.WH..
..W.W...
..WGW...
........
........
........
........
........

Sample Output

1
2
Case 1: 12
Case 2: 7

Solution

斯坦納樹為 NP-hard 問題,因此沒有多項式解法,而要得到確切的最小化解,則必須透過動態規劃來完成。由於題目給定的範圍很小,首先將目標要連通完成的點集,壓縮成 $N$ 個位元,接著紀錄這個聯通元件的其中一個節點視為根。最後,得到狀態數為 $M \cdot 2^N$,可以參考 《「Steiner tree problem in graphs」斯坦納樹》 -日月卦長 的說明。

公式可以拆成兩種情況,第一種為從子集合併中著手,另一種為拓展連通元件 (替換根節點,但不改變目前已經連到的目標集合)。定義 dp[s][i] 為根 i,連通集合 s 的最小花費

  • $dp[S][i] = \min(dp[T][j]+dp[S−T][j]+\text{dist}(i, j):j \in V,T \subset S)$
  • $dp[S][i] = \min(dp[S][k] + \text{dist}(S, k))$

由上述的公式,我們便可知道複雜度為 $O(M \cdot 3^N)$

實作細節

  • 對於內存布局,有兩個選擇 dp[2^N][M] 或者是 dp[M][2^N],其中以 dp[2^N][M] 最為適合,在撰寫迴圈的時候,最內層的迴圈為替換根,這麼一來 cache miss 的機會就非常低。更容易透過 unroll loop 和向量化來運作。

    • 如果內存布局使用 dp[M][2^N],在撰寫向量化時,需要使用 gather 相關的指令,這部分只有在 AVX2 有,並不是每一個 online judge 都支援,而 latency 也算挺高的,等到哪天 CPU 架構換了,這解法可能才會快得起來。
  • 當我們窮舉子集合時,發現到公式有對稱性,便可只窮舉上半部。這樣可以加速 20% 的效能。

1
2
3
4
5
6
7
8
9
10
for (int i = 1, h; i < (1<<n); i++) {
if ((i&(-i)) == i)
h = i-1;
for (int k = (i-1)&i; k > h; k = (k-1)&i) {
for (int j = 0; j < 64; j++)
dp[i][j] = min(dp[i][j], dp[k][j]+dp[i^k][j]-w[j]);
}
relax(i);
}

基礎篇

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
#pragma GCC target("avx")
#pragma GCC optimize ("O3")
#include <bits/stdc++.h>
using namespace std;
char g[8][16], w[64];
static const int dx[] = {0, 0, 1, -1};
static const int dy[] = {1, -1, 0, 0};
static int32_t dp[1<<8][64];
const int INF = 0x3f3f3f3f;
void relax(int s) {
static int8_t Q[1024];
uint64_t inq = 0;
int Qn = 0;
for (int i = 0; i < 64; i++) {
if (dp[s][i] != INF)
Q[Qn++] = i, inq |= 1ULL<<i;
}
for (int i = 0; i < Qn; i++) {
int u = Q[i];
int x = u>>3, y = u&7;
inq ^= 1ULL<<u;
for (int k = 0; k < 4; k++) {
int tx = x+dx[k], ty = y+dy[k];
if (tx < 0 || ty < 0 || tx >= 8 || ty >= 8)
continue;
int v = tx<<3|ty;
if (dp[s][v] > dp[s][u]+1+w[v]) {
dp[s][v] = dp[s][u]+1+w[v];
if (((inq>>v)&1) == 0) {
inq |= 1ULL<<v;
Q[Qn++] = v;
}
}
}
}
}
int main() {
int testcase, cases = 0;
int pl, pw;
scanf("%d", &testcase);
while (testcase--) {
scanf("%d %d", &pl, &pw);
for (int i = 0; i < 8; i++)
scanf("%s", &g[i]);
int n = 0, root = 0;
int A[8];
for (int i = 0; i < 8; i++) {
for (int j = 0; j < 8; j++) {
int u = i<<3|j;
if (g[i][j] == 'H')
A[n++] = u, w[u] = 0;
else if (g[i][j] == 'G')
w[u] = 0, root = u;
else if (g[i][j] == 'W')
w[u] = pw;
else if (g[i][j] == '.')
w[u] = pl;
}
}
memset(dp, 0x3f, sizeof(dp));
for (int i = 0; i < n; i++)
dp[1<<i][A[i]] = 0;
for (int i = 1, h; i < (1<<n); i++) {
if ((i&(-i)) == i)
h = i-1;
for (int k = (i-1)&i; k > h; k = (k-1)&i) {
for (int j = 0; j < 64; j++)
dp[i][j] = min(dp[i][j], dp[k][j]+dp[i^k][j]-w[j]);
}
relax(i);
}
int ret = dp[(1<<n)-1][root];
printf("Case %d: %d\n", ++cases, ret);
}
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
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
#pragma GCC target("avx")
#pragma GCC optimize ("O3")
#include <bits/stdc++.h>
#include <x86intrin.h>
using namespace std;
char g[8][16];
static const int dx[] = {0, 0, 1, -1};
static const int dy[] = {1, -1, 0, 0};
static int32_t w[64] __attribute__ ((aligned(16)));
static int32_t dp[1<<8][64] __attribute__ ((aligned(16)));
const int INF = 0x3f3f3f3f;
void relax(int s) {
static int8_t Q[1024];
uint64_t inq = 0;
int Qn = 0;
for (int i = 0; i < 64; i++) {
if (dp[s][i] != INF)
Q[Qn++] = i, inq |= 1ULL<<i;
}
for (int i = 0; i < Qn; i++) {
int u = Q[i];
int x = u>>3, y = u&7;
inq ^= 1ULL<<u;
for (int k = 0; k < 4; k++) {
int tx = x+dx[k], ty = y+dy[k];
if (tx < 0 || ty < 0 || tx >= 8 || ty >= 8)
continue;
int v = tx<<3|ty;
if (dp[s][v] > dp[s][u]+1+w[v]) {
dp[s][v] = dp[s][u]+1+w[v];
if (((inq>>v)&1) == 0) {
inq |= 1ULL<<v;
Q[Qn++] = v;
}
}
}
}
}
int main() {
int testcase, cases = 0;
int pl, pw;
scanf("%d", &testcase);
while (testcase--) {
scanf("%d %d", &pl, &pw);
for (int i = 0; i < 8; i++)
scanf("%s", &g[i]);
int n = 0, root = 0;
int A[8];
for (int i = 0; i < 8; i++) {
for (int j = 0; j < 8; j++) {
int u = i<<3|j;
if (g[i][j] == 'H')
A[n++] = u, w[u] = 0;
else if (g[i][j] == 'G')
w[u] = 0, root = u;
else if (g[i][j] == 'W')
w[u] = pw;
else if (g[i][j] == '.')
w[u] = pl;
}
}
memset(dp, 0x3f, sizeof(dp));
for (int i = 0; i < n; i++)
dp[1<<i][A[i]] = 0;
for (int i = 1, h; i < (1<<n); i++) {
if ((i&(-i)) == i)
h = i-1;
for (int k = (i-1)&i; k > h; k = (k-1)&i) {
for (int j = 0; j+4 <= 64; j += 4) {
__m128i mv = _mm_load_si128((__m128i*) (dp[i]+j));
__m128i a = _mm_load_si128((__m128i*) (dp[k]+j));
__m128i b = _mm_load_si128((__m128i*) (dp[i^k]+j));
__m128i tm = _mm_add_epi32(a, b);
__m128i c = _mm_load_si128((__m128i*) (w+j));
__m128i tn = _mm_sub_epi32(tm, c);
__m128i mn = _mm_min_epi32(mv, tn);
_mm_store_si128((__m128i*) (dp[i]+j), mn);
}
// dp[i][j] = min(dp[i][j], dp[k][j]+dp[i^k][j]-w[j]);
}
relax(i);
}
int ret = dp[(1<<n)-1][root];
printf("Case %d: %d\n", ++cases, ret);
}
return 0;
}

AVX2 Gather

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
/*
It doesn't work at UVA 2018/08/13 because server CPU does not support
AVX2 instruction set. Although we could pass the compiler, you still
get runtime error during executing an illegal instruction.
*/
#pragma GCC target("avx")
#pragma GCC target("avx2")
#pragma GCC optimize ("O3")
#include <bits/stdc++.h>
#include <x86intrin.h>
#include <avx2intrin.h>
using namespace std;
char g[8][16], w[64];
static const int dx[] = {0, 0, 1, -1};
static const int dy[] = {1, -1, 0, 0};
static int32_t dp[64][1<<8] __attribute__ ((aligned(16)));
const int INF = 0x3f3f3f3f;
void relax(int s) {
static int8_t Q[1024];
uint64_t inq = 0;
int Qn = 0;
for (int i = 0; i < 64; i++) {
if (dp[i][s] != INF)
Q[Qn++] = i, inq |= 1ULL<<i;
}
for (int i = 0; i < Qn; i++) {
int u = Q[i];
int x = u>>3, y = u&7;
inq ^= 1ULL<<u;
for (int k = 0; k < 4; k++) {
int tx = x+dx[k], ty = y+dy[k];
if (tx < 0 || ty < 0 || tx >= 8 || ty >= 8)
continue;
int v = tx<<3|ty;
if (dp[v][s] > dp[u][s]+1+w[v]) {
dp[v][s] = dp[u][s]+1+w[v];
if (((inq>>v)&1) == 0) {
inq |= 1ULL<<v;
Q[Qn++] = v;
}
}
}
}
}
int main() {
int testcase, cases = 0;
int pl, pw;
scanf("%d", &testcase);
while (testcase--) {
scanf("%d %d", &pl, &pw);
for (int i = 0; i < 8; i++)
scanf("%s", &g[i]);
int n = 0, root = 0;
int A[8];
for (int i = 0; i < 8; i++) {
for (int j = 0; j < 8; j++) {
int u = i<<3|j;
if (g[i][j] == 'H')
A[n++] = u, w[u] = 0;
else if (g[i][j] == 'G')
w[u] = 0, root = u;
else if (g[i][j] == 'W')
w[u] = pw;
else if (g[i][j] == '.')
w[u] = pl;
}
}
memset(dp, 0x3f, sizeof(dp));
for (int i = 0; i < n; i++)
dp[A[i]][1<<i] = 0;
for (int i = 1, h; i < (1<<n); i++) {
if ((i&(-i)) == i)
h = i-1;
__attribute__ ((aligned(16))) static int subset1[1<<8] = {};
__attribute__ ((aligned(16))) static int subset2[1<<8] = {};
int sn = 0;
for (int k = (i-1)&i; k > h; k = (k-1)&i)
subset1[sn] = k, subset2[sn] = i^k, sn++;
while (sn&4)
subset1[sn] = 0, subset2[sn] = i, sn++;
for (int j = 0; j < 64; j++) {
int32_t mn = dp[j][i]+w[j];
__m128i mv = _mm_setr_epi32(mn, mn, mn, mn);
for (int t = 0; t <= sn; t += 4) {
int k;
__m128i a = _mm_load_si128((__m128i*) subset1+t);
__m128i b = _mm_load_si128((__m128i*) subset2+t);
__m128i t1 = _mm_i32gather_epi32(dp[j], a, 4);
__m128i t2 = _mm_i32gather_epi32(dp[j], b, 4);
__m128i tm = _mm_add_epi32(t1, t2);
mv = _mm_min_epi32(mv, tm);
}
__attribute__ ((aligned(16))) int32_t mr[4];
_mm_store_si128((__m128i*) mr, mv);
mn = min(min(mr[0], mr[1]), min(mr[2], mr[3]));
dp[j][i] = mn-w[j];
// mn = min(mn, dp[j][k]+dp[j][i^k]-ww);
}
relax(i);
}
int ret = dp[root][(1<<n)-1];
printf("Case %d: %d\n", ++cases, ret);
}
return 0;
}
Read More +

優化技巧 - 排序/優先隊列

主要問題

在算法問題中,時常結合到排序,而排序的常數也大大影響到我們整體的效能。

基數排序降低常數

若是一般原始型別的排序,可以透過基數排序 (radix sort)。從排序的範圍來決定是否要劃分 8-bit 一組一組。若範圍介於可容忍的 $[0, v]$,當 $v < n$ 時,直接開 $n$ 個 bucket 的方式更好。因為大多數的複雜度都大於線性 $O(n)$

非負整數基數排序

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
void _radix_sort(Pt *A, int n) {
static Pt _tmp[MAXN];
const int CHUNK = 256;
static int C[CHUNK];
Pt *B = _tmp, *T;
for (int x = 0; x < 8; x++) {
const int d = x*8;
memset(C, 0, sizeof(C));
for (int i = 0; i < n; i++)
C[(A[i].x>>d)&(CHUNK-1)]++;
for (int i = 1; i < CHUNK; i++)
C[i] += C[i-1];
for (int i = n-1; i >= 0; i--)
B[--C[(A[i].x>>d)&(CHUNK-1)]] = A[i];
T = A, A = B, B = T;
}
}

浮點數基數排序

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
void radix_sort(Pt *A, int n) {
for (int i = 0; i < n; i++) {
int32_t &v = *((int32_t *) &(A[i].s));
if ((v>>31)&1)
v = ~v;
else
v = v | 0x80000000;
}
static Pt _tmp[MAXN*2];
static const int CHUNK = 256;
static int C[1<<8];
Pt *B = _tmp, *T;
for (int x = 0; x < 4; x++) {
const int d = x*8;
memset(C, 0, sizeof(C));
for (int i = 0; i < n; i++)
C[((*((int32_t *) &(A[i].s)))>>d)&(CHUNK-1)]++;
for (int i = 1; i < CHUNK; i++)
C[i] += C[i-1];
for (int i = n-1; i >= 0; i--)
B[--C[((*((int32_t *) &(A[i].s)))>>d)&(CHUNK-1)]] = A[i];
T = A, A = B, B = T;
}
for (int i = 0; i < n; i++) {
int32_t &v = *((int32_t *) &(A[i].s));
if ((v>>31)&1)
v = v & 0x7fffffff;
else
v = ~v;
}
}

慎選內建排序

內建排序常見的有 qsort, sort, stable_sort,我不推薦使用 qsort,因為它在很久以前找到了一個退化情況 A Killer Adversary for Quicksort - 1995,有些題目的測資會出這種特別案例,導致當年的萌新內心受創。除非只有 C 的情況,否則請避開 qsort。如果算法屬於調整某些元素後,再對整體進行排序,這時候 sortstable_sort 慢很多。當情況非常接近已經排序的時候,就使用 stable_sort

優先隊列

使用 priority queue 的方法通常有三種 set, priority_queue, heap

  • 功能最多的 set、次少的 priority_queue,最少的 heap
  • 效能方面則是 heap 最快、次著 priority_queue、最後 set
  • 代碼維護最容易的 set、最差的 heap

之前很排斥使用 priority_queue,原因在於撰寫 operator< 與我的邏輯相反,因此都偏愛使用 set,但是效能被拉出來的時候,又必須退回類似 C 的 make_heappop_heap … 的接口。

Read More +

UVa 12310 - Point Location

Problem

給定一平面上 $n$ 個點,接著拉 $m$ 個線段拼湊成不相交的區域,並將標記數個點位置所在的區域。接著有數個詢問「點在哪一個先前標記的同一區域內」。

這一問題常在幾何處理中出現,詳細可查閱 Wikipedia - Point location,在此不額外說明。

Sample Input

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
14 16 5 5
0 0
30 0
40 0
60 0
60 50
0 50
20 40
10 10
30 10
30 20
40 20
50 30
50 40
30 40
1 2
2 9
9 8
8 7
7 9
9 10
10 11
11 3
3 4
4 5
5 6
6 1
12 13
13 14
14 12
2 3
20 20
10 20
35 10
45 39
1 60
28 11
29 14
34 7
40 38
70 1
0 0 0 0

Sample Output

1
2
3
4
5
1
2
3
4
5

Solution

Previously

在大學修計算幾何的時候,知道這一類的題目要透過梯形剖分的方式,搭配的資料結構類似四分樹,分隔的基準是一個線段,將線段兩個端點分別拉出無限延長的垂直線,但這方法對於線段本身是垂直時,處理情況比較複雜,那時構思了好幾天沒個底,只好放棄。現在,在 Cadence 公司上班快滿一年,解決幾何操作不在少數,成長的我應該能解決這一題了吧。

Detail

題目已經保證給予的線段除了端點外,彼此之間不會相交任何一點。便可以對所有線段建立空間分割相關的資料結構,好讓我們對於解決,從詢問點出發往 $y$ 無窮遠的射線,最後會打到哪一個線段。

處理的流程如下:

Step 1. 輸入線段

Step 2. 套上邊界,需要能包含所有詢問點

Step 3. 對每一個相連的線段集合,找出 $y$ 軸最高的點,放出輔助射線

Step 4. 對每一個頂點進行極角排序,進行繞行分組

首先,對於幾何剖分的問題,為處理方便,都會套上一個無限大的外框,把所有頂點包在裡面。你可以透過 GeoGebra 這套數學軟體,將處理過程中的資料丟進去,方便你除錯。如果要批次輸入一大筆,請透過上方的建立按鈕 (Button),將 command 的語法以行為單位放入,語法類似 javascript。

接著,我們需要對每一個線段集合的最高點放出輔助射線 (這部分已經縮減了不少,原則上對每一個頂點放出射線也可以),這些輔助射線用來解決內部的孤立形,要把整個圖串成一個連通圖。否則,當詢問屬於內部的孤立形狀的外部時,我們會缺少足夠的資訊連到包含這個孤立形的外框。

最後,對每一個頂點相連的邊進行極角排序,接著才能決定相鄰的邊,而任兩個相鄰的邊所夾的區域屬於同一個集合,因此我們需要對每一個有方向的邊進行編號,將邊與邊合併到同一個集合。對於每一個詢問,只需要對詢問點放出射線找到接觸的線段,便可以知道所在的集合。

  • 若使用 KD tree,即使交替選擇 $x$-$y$ 軸分割,也沒辦法保證樹高。因為處理的是線段,而不是點,每一次挑選的分隔軸,將會產生三個子樹,中間的子樹為與分隔軸相交的線段。在大部分的情況下,整體時間複雜度落在 $O(n \log n)$,空間複雜度 $O(n)$
  • 若使用線段樹,將在 $x$ 軸上劃分,每一個線段至多被拆成 $O(\log n)$ 個,詢問射線第一個碰觸的線段時,單一詢問的時間複雜度 $O(\log n)$,常數相較於 KD tree 多。整體時間複雜度落在 $O(n \log n)$,空間複雜度 $O(n \log n)$

還有許多 spatial data structure 可以考慮,這裡只挑選兩個出來實驗,上述皆為在線算法。也可以將所有操作離線,這時候將套用掃描線算法,但對於垂直線段如何在算法中維護二元樹,目前遇到一些實作上的問題,這將在未來才能給各位參考。

KD tree

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
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
#include <bits/stdc++.h>
using namespace std;
const float eps = 1e-6;
const float BOX_MX = 2000000;
struct Pt {
float x, y;
Pt() {}
Pt(float a, float b): x(a), y(b) {}
Pt operator-(const Pt &a) const { return Pt(x - a.x, y - a.y); }
Pt operator+(const Pt &a) const { return Pt(x + a.x, y + a.y); }
Pt operator*(const double a) const { return Pt(x * a, y * a); }
bool operator<(const Pt &a) const {
if (fabs(x - a.x) > eps) return x < a.x;
if (fabs(y - a.y) > eps) return y < a.y;
return false;
}
bool operator==(const Pt &a) const {
return fabs(x - a.x) < eps && fabs(y - a.y) < eps;
}
void println() const {
printf("(%lf, %lf)\n", x, y);
}
};
struct Seg {
Pt s, e;
int i;
Seg() {}
Seg(Pt a, Pt b, int i):s(a), e(b), i(i) {}
void println() {
printf("Segment((%lf, %lf), (%lf, %lf))\n", s.x, s.y, e.x, e.y);
}
};
double dot(Pt a, Pt b) {
return a.x * b.x + a.y * b.y;
}
double cross(Pt o, Pt a, Pt b) {
return (a.x-o.x)*(b.y-o.y)-(a.y-o.y)*(b.x-o.x);
}
double cross2(Pt a, Pt b) {
return a.x * b.y - a.y * b.x;
}
int between(Pt a, Pt b, Pt c) {
return dot(c - a, b - a) >= -eps && dot(c - b, a - b) >= -eps;
}
int onSeg(Pt a, Pt b, Pt c) {
return between(a, b, c) && fabs(cross(a, b, c)) < eps;
}
int cmpZero(float v) {
if (fabs(v) > eps) return v > 0 ? 1 : -1;
return 0;
}
Pt getIntersect(Seg a, Seg b) {
Pt u = a.s - b.s;
double t = cross2(b.e - b.s, u)/cross2(a.e - a.s, b.e - b.s);
return a.s + (a.e - a.s) * t;
}
struct AngleCmp {
Pt o;
AngleCmp(Pt o = Pt()):o(o) {}
bool operator() (const pair<Pt, int>& ppa, const pair<Pt, int>& ppb) {
Pt pa = ppa.first, pb = ppb.first;
Pt p1 = pa - o, p2 = pb - o;
if (p1.y == 0 && p2.y == 0 && p1.x * p2.x <= 0) return p1.x > p2.x;
if (p1.y == 0 && p1.x >= 0 && p2.y != 0) return true;
if (p2.y == 0 && p2.x >= 0 && p1.y != 0) return false;
if (p1.y * p2.y < 0) return p1.y > p2.y;
double c = cross2(p1, p2);
return c > 0 || (c == 0 && fabs(p1.x) < fabs(p2.x));
}
};
static Pt pts[10005];
static Seg segs[30005];
struct BSP {
static const int MAXN = 60005;
static const int MAXNODE = 60005;
struct Node {
float lx, ly, rx, ry;
Node *ls, *ms, *rs;
void extend(Node *u) {
if (u == NULL)
return ;
lx = min(lx, u->lx), ly = min(ly, u->ly);
rx = max(rx, u->rx), ry = max(ry, u->ry);
}
} nodes[MAXNODE];
int sn[MAXNODE];
Seg *seg[MAXNODE];
float axis[MAXN];
Node *root;
int size;
Node* newNode() {
Node *p = &nodes[size++];
assert(size < MAXNODE);
p->ls = p->ms = p->rs = NULL;
sn[p-nodes] = 0;
return p;
}
Node* _build(int k, Seg segs[], int n) {
if (n == 0)
return NULL;
if (k == 2)
k = 0;
int Lsize = 0, Msize = 0, Rsize = 0;
Seg *L = NULL, *M = NULL, *R = NULL;
if (k == 0) {
for (int i = 0; i < n; i++)
axis[i<<1] = segs[i].s.x, axis[i<<1|1] = segs[i].e.x;
nth_element(axis, axis+n, axis+2*n);
const float mval = axis[n];
L = segs;
R = std::partition(segs, segs+n, [mval](const Seg &s) {
return max(s.s.x, s.e.x) <= mval;
});
M = std::partition(R, segs+n, [mval](const Seg &s) {
return min(s.s.x, s.e.x) <= mval;
});
Msize = segs+n - M;
Rsize = M - R;
Lsize = R - segs;
} else {
for (int i = 0; i < n; i++)
axis[i<<1] = segs[i].s.y, axis[i<<1|1] = segs[i].e.y;
nth_element(axis, axis+n, axis+2*n);
const float mval = axis[n];
L = segs;
R = std::partition(segs, segs+n, [mval](const Seg &s) {
return max(s.s.y, s.e.y) <= mval;
});
M = std::partition(R, segs+n, [mval](const Seg &s) {
return min(s.s.y, s.e.y) <= mval;
});
Msize = segs+n - M;
Rsize = M - R;
Lsize = R - segs;
}
Node *u = newNode();
u->lx = BOX_MX, u->ly = BOX_MX;
u->rx = -BOX_MX, u->ry = -BOX_MX;
if (Lsize == n || Rsize == n || Msize == n) {
sn[u - nodes] = n, seg[u - nodes] = segs;
for (int i = 0; i < n; i++) {
u->lx = min(u->lx, min(segs[i].s.x, segs[i].e.x));
u->rx = max(u->rx, max(segs[i].s.x, segs[i].e.x));
u->ly = min(u->ly, min(segs[i].s.y, segs[i].e.y));
u->ry = max(u->ry, max(segs[i].s.y, segs[i].e.y));
}
} else {
u->ls = _build(k+1, L, Lsize), u->extend(u->ls);
u->ms = _build(k+1, M, Msize), u->extend(u->ms);
u->rs = _build(k+1, R, Rsize), u->extend(u->rs);
}
return u;
}
void build_tree(Seg s[], int m) {
size = 0;
root = _build(0, s, m);
}
Pt q_st, q_ed;
int q_si;
void rayhit(Seg &seg) {
if (seg.s.x == seg.e.x) {
if (cmpZero(seg.s.x - q_st.x) == 0) {
double low = min(seg.s.y, seg.e.y);
if (low > q_st.y && low < q_ed.y) {
q_ed.y = low;
q_si = seg.i;
}
}
return ;
}
if (max(seg.s.x, seg.e.x) < q_st.x || min(seg.s.x, seg.e.x) > q_st.x)
return ;
float y = seg.s.y + (float) (seg.e.y - seg.s.y) * (q_st.x - seg.s.x) / (seg.e.x - seg.s.x);
if (y > q_st.y && y < q_ed.y) {
q_ed.y = y;
q_si = seg.i;
}
}
void search(Node *u) {
if (u == NULL)
return ;
if (u->lx > q_st.x || u->rx < q_st.x || u->ry <= q_st.y || u->ly >= q_ed.y)
return ;
for (int i = 0; i < sn[u - nodes]; i++)
rayhit(seg[u - nodes][i]);
search(u->ls);
search(u->ms);
search(u->rs);
}
pair<int, Pt> raycast(Pt st) {
q_st = st;
q_ed = Pt(st.x, BOX_MX+1);
q_si = -1;
search(root);
return {q_si, q_ed};
}
} tree;
struct Disjoint {
static const int MAXN = 65536;
uint16_t parent[MAXN], weight[MAXN];
void init(int n) {
if (n >= MAXN)
exit(0);
for (int i = 0; i <= n; i++)
parent[i] = i, weight[i] = 1;
}
int findp(int x) {
return parent[x] == x ? x : parent[x] = findp(parent[x]);
}
int joint(int x, int y) {
x = findp(x), y = findp(y);
if (weight[x] >= weight[y])
parent[y] = x, weight[x] += weight[y];
else
parent[x] = y, weight[y] += weight[x];
}
} egroup, sgroup;
int main() {
int n, m, p, q;
while (scanf("%d %d %d %d", &n, &m, &p, &q) == 4 && n) {
for (int i = 0; i < n; i++) {
int x, y;
scanf("%d %d", &x, &y);
pts[i] = Pt(x, y);
}
sgroup.init(n);
for (int i = 0; i < m; i++) {
int st_i, ed_i;
scanf("%d %d", &st_i, &ed_i);
st_i--, ed_i--;
segs[i] = Seg(pts[st_i], pts[ed_i], i);
sgroup.joint(st_i, ed_i);
}
segs[m] = Seg(Pt(BOX_MX, BOX_MX), Pt(BOX_MX, -BOX_MX), m), m++;
segs[m] = Seg(Pt(BOX_MX, -BOX_MX), Pt(-BOX_MX, -BOX_MX), m), m++;
segs[m] = Seg(Pt(-BOX_MX, -BOX_MX), Pt(-BOX_MX, BOX_MX), m), m++;
segs[m] = Seg(Pt(-BOX_MX, BOX_MX), Pt(BOX_MX, BOX_MX), m), m++;
static map<Pt, vector<pair<Pt, uint8_t>>> g; g.clear();
static vector<vector<Pt>> on_seg; on_seg.clear(); on_seg.resize(m);
for (int i = 0; i < m; i++) {
on_seg[segs[i].i].reserve(4);
on_seg[segs[i].i].push_back(segs[i].s);
on_seg[segs[i].i].push_back(segs[i].e);
}
// for (int i = 0; i < n; i++)
// pts[i].println();
// for (int i = 0; i < m; i++)
// segs[i].println();
tree.build_tree(segs, m);
{
Pt top[10005];
for (int i = 0; i < n; i++)
top[i] = pts[i];
for (int i = 0; i < n; i++) {
int gid = sgroup.findp(i);
if (top[gid].y < pts[i].y)
top[gid] = pts[i];
}
for (int i = 0; i < n; i++) {
if (sgroup.findp(i) != i)
continue;
auto p = top[i];
auto hit = tree.raycast(p);
if (hit.first >= 0) {
on_seg[hit.first].emplace_back(hit.second);
g[p].emplace_back(hit.second, 1);
g[hit.second].emplace_back(p, 1);
}
}
}
for (int i = 0; i < m; i++) {
vector<Pt> &a = on_seg[i];
sort(a.begin(), a.end());
a.resize(unique(a.begin(), a.end()) - a.begin());
auto *prev = &g[a[0]];
for (int j = 1; j < a.size(); j++) {
prev->emplace_back(a[j], 0);
prev = &g[a[j]];
prev->emplace_back(a[j-1], 0);
}
}
for (auto &e : g)
sort(e.second.begin(), e.second.end(), AngleCmp(e.first));
static map<Pt, map<Pt, int>> R; R.clear();
int Rsize = 0;
for (auto &e : g) {
int sz = e.second.size();
map<Pt, int> &Rg = R[e.first];
for (auto &f : e.second) {
int &eid = Rg[f.first];
if (eid == 0)
eid = ++Rsize;
}
}
egroup.init(Rsize);
for (auto &e : g) {
int sz = e.second.size();
map<Pt, int> &Rg = R[e.first];
for (int i = sz-1, j = 0; j < sz; i = j++) {
int l = R[e.second[i].first][e.first];
int r = Rg[e.second[j].first];
egroup.joint(l, r);
if (e.second[i].second != 0) {
r = Rg[e.second[i].first];
assert(l > 0 && r > 0);
egroup.joint(l, r);
}
}
}
for (auto &e : g) {
int sz = e.second.size();
int n = 0;
for (int i = 0; i < sz; i++) {
if (e.second[i].second == 0)
e.second[n++] = e.second[i];
}
e.second.resize(n);
}
static int region[65536]; memset(region, 0, sizeof(0)*Rsize);
for (int i = 0; i < p; i++) {
float x, y;
scanf("%f %f", &x, &y);
pair<int, Pt> hit = tree.raycast(Pt(x, y));
if (hit.first < 0)
continue;
if (g.find(hit.second) != g.end()) {
auto &adj = g[hit.second];
AngleCmp cmp(hit.second);
int pos = 0;
pair<Pt, int> q(Pt(x, y), i+1);
pos = lower_bound(adj.begin(), adj.end(), q, cmp) - adj.begin() - 1;
assert(pos >= 0);
int l = R[adj[pos].first][hit.second];
assert(l > 0);
l = egroup.findp(l);
region[l] = i+1;
} else {
auto &adj = on_seg[hit.first];
Pt lpt = adj[0], rpt = adj[1];
int l;
if (cmpZero(cross(lpt, rpt, Pt(x, y))) <= 0)
l = R[lpt][rpt];
else
l = R[rpt][lpt];
assert(l > 0);
l = egroup.findp(l);
region[l] = i+1;
}
}
for (int i = 0; i < q; i++) {
float x, y;
scanf("%f %f", &x, &y);
pair<int, Pt> hit = tree.raycast(Pt(x, y));
if (hit.first < 0) {
printf("0\n");
continue;
}
if (g.find(hit.second) != g.end()) {
auto &adj = g[hit.second];
AngleCmp cmp(hit.second);
int pos = 0;
pair<Pt, int> q(Pt(x, y), i+1);
pos = lower_bound(adj.begin(), adj.end(), q, cmp) - adj.begin() - 1;
assert(pos >= 0);
int l = R[adj[pos].first][hit.second];
assert(l > 0);
l = egroup.findp(l);
printf("%d\n", region[l]);
} else {
auto &adj = on_seg[hit.first];
Pt lpt = adj[0], rpt = adj[1];
int l;
if (cmpZero(cross(lpt, rpt, Pt(x, y))) <= 0)
l = R[lpt][rpt];
else
l = R[rpt][lpt];
assert(l > 0);
l = egroup.findp(l);
printf("%d\n", region[l]);
}
}
}
return 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
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
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
#include <bits/stdc++.h>
using namespace std;
const double eps = 1e-8;
const double BOX_MX = (3000000);
const int AUXID = 200000;
struct Pt {
double x, y;
Pt(double a = 0, double b = 0): x(a), y(b) {}
Pt operator-(const Pt &a) const { return Pt(x - a.x, y - a.y); }
Pt operator+(const Pt &a) const { return Pt(x + a.x, y + a.y); }
Pt operator*(const double a) const { return Pt(x * a, y * a); }
bool operator<(const Pt &a) const {
if (fabs(x - a.x) > eps) return x < a.x;
if (fabs(y - a.y) > eps) return y < a.y;
return false;
}
bool operator==(const Pt &a) const {
return fabs(x - a.x) < eps && fabs(y - a.y) < eps;
}
void println() const {
printf("(%lf, %lf)\n", x, y);
}
};
struct Seg {
Pt s, e;
int i;
Seg() {}
Seg(Pt a, Pt b, int i):s(a), e(b), i(i) {}
void println() const {
printf("Segment((%lf, %lf), (%lf, %lf))\n", s.x, s.y, e.x, e.y);
}
};
double dot(Pt a, Pt b) {
return a.x * b.x + a.y * b.y;
}
double cross(Pt o, Pt a, Pt b) {
return (a.x-o.x)*(b.y-o.y)-(a.y-o.y)*(b.x-o.x);
}
double cross2(Pt a, Pt b) {
return a.x * b.y - a.y * b.x;
}
int between(Pt a, Pt b, Pt c) {
return dot(c - a, b - a) >= -eps && dot(c - b, a - b) >= -eps;
}
int onSeg(Pt a, Pt b, Pt c) {
return between(a, b, c) && fabs(cross(a, b, c)) < eps;
}
int cmpZero(double v) {
if (fabs(v) > eps) return v > 0 ? 1 : -1;
return 0;
}
Pt getIntersect(Seg a, Seg b) {
Pt u = a.s - b.s;
double t = cross2(b.e - b.s, u)/cross2(a.e - a.s, b.e - b.s);
return a.s + (a.e - a.s) * t;
}
struct AngleCmp {
Pt o;
AngleCmp(Pt o = Pt()):o(o) {}
bool operator() (const pair<Pt, int>& ppa, const pair<Pt, int>& ppb) {
Pt pa = ppa.first, pb = ppb.first;
Pt p1 = pa - o, p2 = pb - o;
if (p1.y == 0 && p2.y == 0 && p1.x * p2.x <= 0) return p1.x > p2.x;
if (p1.y == 0 && p1.x >= 0 && p2.y != 0) return true;
if (p2.y == 0 && p2.x >= 0 && p1.y != 0) return false;
if (p1.y * p2.y < 0) return p1.y > p2.y;
double c = cross2(p1, p2);
return c > 0 || (c == 0 && fabs(p1.x) < fabs(p2.x));
}
};
static Pt pts[10005];
static Seg segs[30005];
static double interpolate(const Pt& p1, const Pt& p2, double& x) {
if (cmpZero(p1.x - p2.x) == 0) return min(p1.y, p2.y);
return p1.y + (double)(p2.y - p1.y) / (p2.x - p1.x) * (x - p1.x);
}
struct CMP {
static double x;
bool operator() (const Seg &i, const Seg &j) {
double l = interpolate(i.s, i.e, x);
double r = interpolate(j.s, j.e, x);
return l < r;
}
};
double CMP::x;
struct BSP {
struct Node {
double ly, ry;
void init() {
ly = BOX_MX, ry = -BOX_MX;
}
void extend(Node *u) {
if (u == NULL)
return ;
ly = min(ly, u->ly);
ry = max(ry, u->ry);
}
} bbox[262144];
vector<double> x;
set<Seg, CMP> tree[524288];
Seg segs[262144];
void clear(int k, int l, int r) {
tree[k].clear();
bbox[k].init();
if (l+1 >= r)
return ;
int m = (l+r)/2;
clear(k<<1, l, m);
clear(k<<1|1, m, r);
}
void insert(Seg &s, int k, int l, int r) {
if (s.s.x <= x[l] && x[r] <= s.e.x) {
CMP::x = (x[l] + x[r])/2;
tree[k].insert(s);
double ly = interpolate(s.s, s.e, x[l]);
double ry = interpolate(s.s, s.e, x[r]);
bbox[k].ly = min(bbox[k].ly, min(ly, ry));
bbox[k].ry = max(bbox[k].ry, max(ly, ry));
return ;
}
if (l+1 >= r)
return ;
int m = (l+r)/2;
if (s.s.x <= x[m]) {
insert(s, k<<1, l, m);
bbox[k].extend(&bbox[k<<1]);
}
if (s.e.x > x[m]) {
insert(s, k<<1|1, m, r);
bbox[k].extend(&bbox[k<<1|1]);
}
}
void build_tree(Seg s[], int m) {
memcpy(segs, s, sizeof(segs[0])*m);
x.clear();
for (int i = 0; i < m; i++)
x.push_back(segs[i].s.x), x.push_back(segs[i].e.x);
sort(x.begin(), x.end());
x.resize(unique(x.begin(), x.end()) - x.begin());
clear(1, 0, x.size()-1);
for (int i = 0; i < m; i++) {
if (segs[i].s.x > segs[i].e.x)
swap(segs[i].s, segs[i].e);
if (segs[i].s.x < segs[i].e.x)
insert(segs[i], 1, 0, x.size()-1);
}
}
Pt q_st, q_ed;
int q_si;
void search(int k, int l, int r) {
if (bbox[k].ly >= q_ed.y || bbox[k].ry <= q_st.y)
return ;
if (x[l] <= q_st.x && q_st.x <= x[r]) {
CMP::x = q_st.x;
auto it = tree[k].upper_bound(Seg(q_st, q_st, 0));
while (it != tree[k].end()) {
double y = interpolate(it->s, it->e, CMP::x);
if (y > q_st.y) {
if (y < q_ed.y) {
q_ed.y = y;
q_si = it->i;
}
break;
}
it++;
}
}
if (l+1 >= r)
return ;
int m = (l+r)/2;
if (q_st.x <= x[m])
search(k<<1, l, m);
if (q_st.x >= x[m])
search(k<<1|1, m, r);
}
pair<int, Pt> raycast(Pt st) {
q_st = st;
q_ed = Pt(st.x, BOX_MX+1);
q_si = -1;
search(1, 0, x.size()-1);
return {q_si, q_ed};
}
} tree;
struct Disjoint {
static const int MAXN = 65536;
uint16_t parent[MAXN], weight[MAXN];
void init(int n) {
if (n >= MAXN)
exit(0);
for (int i = 0; i <= n; i++)
parent[i] = i, weight[i] = 1;
}
int findp(int x) {
return parent[x] == x ? x : parent[x] = findp(parent[x]);
}
int joint(int x, int y) {
x = findp(x), y = findp(y);
if (weight[x] >= weight[y])
parent[y] = x, weight[x] += weight[y];
else
parent[x] = y, weight[y] += weight[x];
}
} egroup, sgroup;
int main() {
int n, m, p, q;
while (scanf("%d %d %d %d", &n, &m, &p, &q) == 4 && n) {
for (int i = 0; i < n; i++) {
int x, y;
scanf("%d %d", &x, &y);
pts[i] = Pt(x, y);
}
sgroup.init(n);
for (int i = 0; i < m; i++) {
int st_i, ed_i;
scanf("%d %d", &st_i, &ed_i);
st_i--, ed_i--;
segs[i] = Seg(pts[st_i], pts[ed_i], i);
sgroup.joint(st_i, ed_i);
}
segs[m] = Seg(Pt(BOX_MX, BOX_MX), Pt(BOX_MX, -BOX_MX), m), m++;
segs[m] = Seg(Pt(BOX_MX, -BOX_MX), Pt(-BOX_MX, -BOX_MX), m), m++;
segs[m] = Seg(Pt(-BOX_MX, -BOX_MX), Pt(-BOX_MX, BOX_MX), m), m++;
segs[m] = Seg(Pt(-BOX_MX, BOX_MX), Pt(BOX_MX, BOX_MX), m), m++;
static map<Pt, vector<pair<Pt, uint8_t>>> g; g.clear();
static vector<vector<Pt>> on_seg; on_seg.clear(); on_seg.resize(m);
for (int i = 0; i < m; i++) {
on_seg[segs[i].i].reserve(4);
on_seg[segs[i].i].push_back(segs[i].s);
on_seg[segs[i].i].push_back(segs[i].e);
}
// for (int i = 0; i < n; i++)
// pts[i].println();
// for (int i = 0; i < m; i++)
// segs[i].println();
tree.build_tree(segs, m);
{
Pt top[10005];
for (int i = 0; i < n; i++)
top[i] = pts[i];
for (int i = 0; i < n; i++) {
int gid = sgroup.findp(i);
if (top[gid].y < pts[i].y)
top[gid] = pts[i];
}
for (int i = 0; i < n; i++) {
if (sgroup.findp(i) != i)
continue;
auto p = top[i];
auto hit = tree.raycast(p);
if (hit.first >= 0) {
on_seg[hit.first].emplace_back(hit.second);
g[p].emplace_back(hit.second, 1);
g[hit.second].emplace_back(p, 1);
}
}
}
for (int i = 0; i < m; i++) {
vector<Pt> &a = on_seg[i];
sort(a.begin(), a.end());
a.resize(unique(a.begin(), a.end()) - a.begin());
auto *prev = &g[a[0]];
for (int j = 1; j < a.size(); j++) {
prev->emplace_back(a[j], 0);
prev = &g[a[j]];
prev->emplace_back(a[j-1], 0);
}
}
for (auto &e : g)
sort(e.second.begin(), e.second.end(), AngleCmp(e.first));
static map<Pt, map<Pt, int>> R; R.clear();
int Rsize = 0;
for (auto &e : g) {
int sz = e.second.size();
map<Pt, int> &Rg = R[e.first];
for (auto &f : e.second) {
int &eid = Rg[f.first];
if (eid == 0)
eid = ++Rsize;
}
}
egroup.init(Rsize);
for (auto &e : g) {
int sz = e.second.size();
map<Pt, int> &Rg = R[e.first];
for (int i = sz-1, j = 0; j < sz; i = j++) {
int l = R[e.second[i].first][e.first];
int r = Rg[e.second[j].first];
egroup.joint(l, r);
if (e.second[i].second != 0) {
r = Rg[e.second[i].first];
assert(l > 0 && r > 0);
egroup.joint(l, r);
}
}
}
for (auto &e : g) {
int sz = e.second.size();
int n = 0;
for (int i = 0; i < sz; i++) {
if (e.second[i].second == 0)
e.second[n++] = e.second[i];
}
e.second.resize(n);
}
static int region[65536]; memset(region, 0, sizeof(0)*Rsize);
for (int i = 0; i < p; i++) {
float x, y;
scanf("%f %f", &x, &y);
pair<int, Pt> hit = tree.raycast(Pt(x, y));
if (hit.first < 0)
continue;
if (g.find(hit.second) != g.end()) {
auto &adj = g[hit.second];
AngleCmp cmp(hit.second);
int pos = 0;
pair<Pt, int> q(Pt(x, y), i+1);
pos = lower_bound(adj.begin(), adj.end(), q, cmp) - adj.begin() - 1;
assert(pos >= 0);
int l = R[adj[pos].first][hit.second];
assert(l > 0);
l = egroup.findp(l);
region[l] = i+1;
} else {
auto &adj = on_seg[hit.first];
Pt lpt = adj[0], rpt = adj[1];
int l;
if (cmpZero(cross(lpt, rpt, Pt(x, y))) <= 0)
l = R[lpt][rpt];
else
l = R[rpt][lpt];
assert(l > 0);
l = egroup.findp(l);
region[l] = i+1;
}
}
for (int i = 0; i < q; i++) {
float x, y;
scanf("%f %f", &x, &y);
pair<int, Pt> hit = tree.raycast(Pt(x, y));
if (hit.first < 0) {
printf("0\n");
continue;
}
if (g.find(hit.second) != g.end()) {
auto &adj = g[hit.second];
AngleCmp cmp(hit.second);
int pos = 0;
pair<Pt, int> q(Pt(x, y), i+1);
pos = lower_bound(adj.begin(), adj.end(), q, cmp) - adj.begin() - 1;
assert(pos >= 0);
int l = R[adj[pos].first][hit.second];
assert(l > 0);
l = egroup.findp(l);
printf("%d\n", region[l]);
} else {
auto &adj = on_seg[hit.first];
Pt lpt = adj[0], rpt = adj[1];
int l;
if (cmpZero(cross(lpt, rpt, Pt(x, y))) <= 0)
l = R[lpt][rpt];
else
l = R[rpt][lpt];
assert(l > 0);
l = egroup.findp(l);
printf("%d\n", region[l]);
}
}
}
return 0;
}
Read More +

十二天軍旅生活 (前)

畢業至當兵前

結束研究所離校手續後,便回到家裡等兵單。某方面來說不算等兵單,研發替代役可以上網填入營日期,選擇適合自己的梯次入伍即可。彌補好幾個寒暑假都沒回家的罪過,於是跑去報名汽車駕訓班,沒想到現在學開車外加考照要一個多月,那麼先前填的入伍時間來不及,在確定報名的日期前修改即可,往後延了一個梯次,直到九月中才進去服十二天的新訓。

駕訓班

在這一個多月間,駕訓班平均兩天上一次課,每天排課差不多都是早上八、九點,其中突然有一陣子沒辦法練車,因為有某一梯次的學員要進行考試,好幾天就空著練筆試。上班時期才去學開車也許是個艱難的選項,而我們花蓮的駕訓班排課好像只到晚上六點,想必下班也會來不及。與公司、役政署那邊確定好時間,也要把這一人生階段任務完成。

大家都說以前駕訓班就學手排車,於是一開始學手排車,後來發現對離合器的協調性不足,換自排車繼續學。開了手排車幾堂課,轉換到自排車就覺得相當簡單,自排車比較不會打錯檔位,也許是因為檔位設計的問題吧!哪天教練車高檔一點說不定就好辦了。咱這點自尊先捨棄吧。

不曉得其他駕訓班是怎麼教開車,背口訣打方向盤的圈數、看後照鏡與標線之間的關係,一開始總是憑感覺打方向盤,結果常常被教練說嘴。心想那些口訣換了車就不一樣,而且每個人坐在駕駛坐看到的視角關係也不致相同,我到底該信什麼才好?內心好煎熬啊。

有時開車兩個人輪流開一台車,兩個學員相互看開車的情況,也不知道為什麼幾乎都是女孩子,而且開車技術都還算不錯,像頭文字 D 那般的快速打檔,感覺平常就有開車,坐在一旁都好佩服。後來才發現教練幾乎都帶女生,怪不得安排上都是跟女生居多。開 S 灣的時候壓到線時,當一旁人的一言不發,腦補病患可承受不住的。

直到考試當日,當日早上八點集合帶去監理站考試,筆試背罰金額度的總是相當惱人,那一天死記活拚也弄了個 97.5 分,明明不該是這麼困難的項目,帶著沒拿到滿分的遺憾迎來下午的場內考試。場內考試依舊是各教練分開安排,排隊上車進行測驗,由於監理站人手不足,其他教練評分時仍有全程攝影,待考生坐在後座看前一個考生跑完全程,六十多個人坐在安全島上等待考試。各自帶開後,不出所料這島都是女生,而另一頭都是男生,教練們的合格率賭盤就開始了!不得不說,大部分人都是考自排車,非常少數的人考手排車,大部分上坡起步退場的也都是手排車。

一般而言,大家都以為這樣就結束。「不是的」,像我這種總是可以走到特殊路線的人而言,恰好多了道路駕駛考試,也就是在今年六月開始強制上路的考試,必須在外頭繞一圈合格後才能拿到駕照。好在是在花蓮考試,如果在台北考道路駕駛也許不堪設想。開車出去前要口頭報出車況,我永遠不能明白當初定說要檢查車底的高就在想什麼,要蹲下去看車底高喊「車底無異物」,並且是車子四個方位都要走過去喊,想必在外頭沒有一個人可以做到這一點。為了安全考量看一次,但看到四次只有強迫症患者做得到吧。

說起如何練習道路考試,那是一件很有趣的故事,一般人會想說教練會帶著進行第一次的練習,結果我們教練要求讓家長帶我們去開過一次考試路線,不是家長開給你開,而是你開給家長看。多嚇過幾輪後,教練才會帶你出去練車。然而,礙於其他因素就沒拜託老爸帶我去開。於是-第一次開車上路,沒有經驗的我開車載妹子出去啦,教練還說「把這當載老婆小孩,一句話都不會說,你得一個人開車。」,幸好平安落幕。隨後幾次練車,曾跟著奇葩一起,坐在後頭都會嚇死,教練說道「輕踩壓車減速」,心中吶喊「等等,咱們突然停在路上,停了,真的停了」,原來他沒有騎過機車、腳踏車上路的經驗,當初還以為大家都會騎機車才考汽車,這世界還是很奇妙的。

拉著老爸練了一周的車,考試當天等了三個多小時才輪到最後一個我,坐在前一個考生後頭時,從頭到尾一直在喊「饒了咱吧,快把你的雙手放開,疊在一起放在方向盤上搖來搖去,咱的心臟要受不了」一個急轉彎壓了慢車道過去,全程超速駕駛,考官全程也不能說話,不知道那樣子還算過嗎。等到我開完路線,監考官叫我多開一段路回駕訓班,那時教練們聚著才開始抱怨「現在的年輕人啊…」

老爸拉著我去開台9、台11線,中間繞光豐公路,練習在山路上轉方向盤,開車還要閃遊覽車之類的驚悚橋段,一旁就是山溝,還有道路坍方的單向道,想起來還是得再多多加強。要練到開山路不踩剎車,再開始載人?

練身體

雖然這年頭研發替代役只有十二天的新訓,人生這一階段還是練練身體。每天下午跑去附近公園跑個六、七圈,順便拉著老爸去運動。

覺得慢跑的最重要因素還是腿部和腰部肌肉。高中時間肥到 76kg 時,跑個 1600m 都覺得辛苦,剛上大學覺得自己是個一無是處廢物,大一的生活就有固定跑到宿舍旁邊的操場運動,「因為太廢了,這八圈是懲罰」
。然而,轉學後就一陣子沒運動,直到大三開始跟轉學生學長一起住,瘋了才去跑校園兩圈練半半馬 10km,後來還真的跑去參加半半馬,一個人在人群中擠來擠去。但上了研究所兩年完全沒運動過,每天走三十分鐘往返學校而已。突然跑起 3000m 還是很刺激的,調個一周才讓大腿適應。

後來被老爸抓去花蓮太平洋公園旁的道路跑步,第一次在下午三點太陽正熱,即使有海風吹著跑起來還是相當折磨,單程約有 2km,來回跑就有四公里多,在海邊跑步如此現充?不對,觀光客騎腳踏車、機車超輕鬆的從你旁邊呼嘯而過,而你只有看著狹長的海岸線估算要以何種步伐才得以抵達盡頭的份。

Read More +

吶,還記得我們的「約定」嗎

長途跋涉

碩二下的剩餘幾個月裡,一般情況在口試前一個月開始把實驗補齊,前一週才把論文寫完,隨後口試完,再花個兩三週把論文丟到圖書館後離校,以上是一般碩士生的的常態生活,而我可能就有一點特別。碩二上學期結束時把論文概要想完,過年把中文論文寫個草稿,下學期便耗費數幾個月與老師來來回回地改了好幾輪,一下子就來到了學期結束。

「長篇化就拜託你了」-《情色漫畫老師》

不得不說,內容改了好幾個月,把一些不太想說明、為什麼要說的部分釐清變得相當困難,大部分不想說明的內容、晦暗不明的語句依依剝離出來,反反覆覆地增加了四十幾張圖片,使用 Tikz 這等套件畫圖,可說是相當痛苦得一件事情,一張圖片可能就要花一整個早上來完成。

細節實作的「道理」與其實驗結果,逐一補上「為什麼要這麼做」、「為什麼會快」、「為什麼選用這種方法」諸如此類的問題。在撰寫時只覺得這麼做是對的,感覺上沒有錯,實驗上可以驗證是對的,一被問起來仍是一臉茫然地不知如何說明,還得回去思考到底要怎麼說才能容易明白。

「這不是超過 100 頁了嗎?」-《情色漫畫老師》

一路寫到五月底,而六月初口試的當下,口試現場來了四個教授來聽,相較於學長那時的三人排場的確更盛大點。壓縮簡報內容至一百頁內的投影片,才可能在四十分鐘內講完呢。那些不想細講的內容,在論文裡不說又會被老師要求,口試的時候能不能不說?百感交集下,硬著頭皮一路殺下去講純資料結構和算法的改進,對於一般做計畫而生的論文,我的內容肯定很無趣吧。要是我能強一點就好了。

「要是我能更強就好了」-《夏目友人帳》

講完近一個小時的成果發表,講累的我還一度沒聽懂要開離教室,等待委員商議是否能通過的會議,口試的評分表不知道是以怎麼方式運行,坐在教室外的我焦慮地等待,心想「是不是做的內容有些瑕疵,還是已經有人做過了?還是內容不足以作為一個碩士研究?」

教室那鎖不好的木門聲響起「嘣」一聲,老師從門口向著我說「下去把 XX 一起叫上來」,走下樓把一起口試的同學都叫來教室,才聽到「恭喜你們通過口試」,讀了兩年就等待這一句話,那究竟是怎麼樣的體驗,感覺有點不切實際,就這樣子了?

「這是奇蹟」-《正解的卡多》

儘管已經提早口試,內心再怎麼雀躍,也無法順心地馬上離開這裡,每天依舊來來學校做點雜事,「鈴~鈴~」實驗室電話一響,電話那一頭「有其他人在嗎?那個誰在嗎?麻煩上來一下」,早上十點半的我只能無奈地「誰也不在,現在只有我」,收到「好吧,那你上來」的日常似乎沒有太大的差異。

口試通過後,就該加緊腳步離校了對吧?離校的兩大關卡——口試審定書和離校同意書,過了一關後,還有一關無法通過,那要怎麼通過老師給的離校同意書呢?方法可說是千奇百種啊,男生的最終武器應該是兵單吧,不然總是問最晚什麼時候你才會離開而避開要簽離校,也許還有另一種可能,讓老師覺得「留你也沒用」,幻想著跟老師說「老師啊,就算沒有我,地球也照樣繼續轉,留我做什麼呢?」

「就算沒有我,地球也照樣繼續轉」—《櫻花任務》

六月初那時,別老是問我「為什麼不走?」,心想「如果能走的話,早就離開了」,這麼能輕易明白的道理,在我身上就這麼難以理解,萬事皆有因,而螢幕前的你是否能理解我呢?看著行事曆上的工作吧,「六月十九日,校定期末考週」這下明白了吧,原本以為可以託付給學弟撐著弄分散式部分,沒想到碩一期末時期這麼忙碌,弄個作業就兩三天不見人影,實驗室好幾週進入三班工作制。

撐過了幾週的工作準備和實驗環境修整,期末考週終於可以好好休息,考題總要給老師弄好,助教輕鬆地去監考就行了。期末考當天監考完,由於沒有確切的參考答案,當下也不知道怎麼向同學描述怎麼寫才是對的,要提示到哪種地步,當下感慨萬分,「對不起,這裡你要自己想,我不清楚」。

原以為老師很快就可以改好考卷,沒想到突如其來的論文 deadline,使得改考卷一事被拖延,然而送成績又是一週後的事情,先放個三天讓老師去改。催一下進度,發現一題也沒改完,先改了兩題,剩下的題目再給老師改。在放個兩天,還是沒有太多的進展,再跟老師要了兩題來改,到最後還是由我改了大部分的題目,為 … 為什麼會變成那樣?

「為 ... 為什麼會變成那樣」-《情色漫畫老師》

七月初到,仍然沒有辦法離開學校,事情接踵而來,弄完轉系考才能離開,又有實驗室經費要花,交代事情的當下只有我在學校,「嗯?怎麼?」對於研究生要在一早十點九點到實驗室,想必是個相當難的議題,學長們可是曾經全員不在而被叫上去罵呢,而我嘗試支持了一年多,電話從我座位旁響起到接起已成為最佳化的結果。在空無一人的早晨接電話,對於指定任務交代給尚未清醒的學弟等。買伺服器時,得再三催廠商寄送發票;買耗材時,得催學弟別忘記。

「太奇怪了」-《愛麗絲與藏六》

終於到辦公室跟老師提了兵單,才從老師口中說出「下週會簽離校,順便回去實驗室跟他們說」一回到實驗室說,便聽到「為什麼老師不願意簽我離校啦」,說出「下週可以簽離校」,感覺同學得知這消息倍感高興。而在那週開會的結束時,跟老師報告做的事情時,僅剩下我一人了,順帶把離校同意書一起簽。

回到實驗室,滿是「可以離校了」「什麼時候清東西」的議題,這些話題平時離我太遠,使得我無法參與討論。學弟反過來問「那你什麼時候打算簽?」說起這個傷心事,原來我總是比較特別的部分對吧!我不知道要提些什麼,悄敲地回應「為什麼覺得我沒簽呢?」然而,在簽之前,比我晚口試的古學長早在前幾天拿到畢業證書,而我才剛拿到手。不經開始懷疑自己的能力與價值。要從口中說出「拿到了」,內心可說是五味雜陳地說不出啊。

「我拿到了」-《情色漫畫老師》

突襲再訪

敬請期待

再續故事

敬請期待

可能走向

敬請期待

誌謝

以下收錄於論文中

能順利地完成這篇論文,首先,感謝劉邦鋒和吳真貞老師的指導,在論文架構和描
述手法的教導,才使得篇幅雜亂的初稿便得更加地易於理解,討論過程中更加地精
練專業知識,學生在此衷心感謝老師。

特別感謝中央大學的郭人維學長,拉拔在演算法及資料結構領域上的研究,其給予
的助力使得論文開花結果。更感謝在網路上許許多多來自各方的朋友分享研究心得
,讓彼此切磋茁壯。

在進入臺灣大學研究所的這兩年中,感謝實驗室學長們的鼓勵與支持,本對於學術
研究文化和自身能力的迷茫,接受學長們的啟示後,最終得以撐過第一學年。在第
二年中,感謝共同奮戰論文的古耕竹、古君葳、鄭以琳、吳軒衡等同學,彼此加油
打氣,使得在撰寫論文的路上並不孤單,督促進展、協助撰寫與驗證想法更加地順
利,願你們也能順利畢業、研究出滿意的成果。

在研究實驗上,感謝實驗室學弟張逸寧、林明璟、蔡慶源、朱清福的參與,被迫實
作出放置於批改娘系統上題目,這些題目原本為論文的一小部分,可透過不同資料
結構與算法解決,在本篇追求效能極致的路上貢獻了一份心力,為本實驗結果給予
更有信心的立論基礎。願你們在接下來的一年裡,經過老師指導與同學們相互引領
下順利畢業。

此外,特別感謝高中時期帶入門的溫健順老師,在選擇領域分組時,相信我在資訊
領域上的發展,拉近資訊組培養,經歷三年的教導後,才能順利走向這一條道路。

最後,感謝家人們一路相伴,進入資訊工程領域後,經歷大學轉學、延畢到研究所
的路上,對我的選擇給予支持。

感謝上述的各位與師長們一路上的支持與資助。

Read More +

批改娘 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 +