手札日記

純粹 Part 1

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

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

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

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

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

感情篇

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

又過了一陣子，某個晚上

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

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

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

又過了一陣子，某個下午

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

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

「應該可以吧？」

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

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

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

「我的論文重要？還是學弟作業重要！」

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

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

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

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

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

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

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

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

Read More +
解題區/解題區 - UVa

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 +
學校課程/計算機圖學

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

此為計算機圖學課程中的一環，自由挑選主題。而這個 Radiosity 有好幾年前的原始碼，程式碼部分由老師提供。看到幾處寫得不是很好的地方，於是就拿來加速。在 tracing code 耗費了相當大的力氣，雖然才幾千行的 C 程序，完全不熟的領域，慢慢做起來也別有一番風趣。morris821028/hw-radiosity

Parallel Computing and Optimization Skills for Radiosity

  • R04922067 楊翔雲、R04922133 古君葳

介紹 Introduction

熱輻射法 (Radiosity) 是一種渲染的技術。相較於光線追蹤法 (Ray Tracing)，熱輻射法可以產生更接近於現實場景中光亮的變化。當場景使用光線追蹤法時，物體的陰影的邊緣相對銳利，但在現實情況下，物體陰影漸層呈現，因此使用熱輻射法可以更貼近我們想要的。

$$\begin{align*} B_i dA_i = E_i dA_i + R_i \int_j B_j F_{ji} dA_j \end{align*}$$
  • $A_i$ : Area of element i (computable)
  • $B_i$ : Radiosity of element i (unknown)
  • $E_i$ : Radient emitted flux density of element i (given)
  • $R_i$ : Refletance of element i (given)
  • $F_{ji}$ : Form Factor from j to i (computable)

假設整個場景中有 $N$ 個三角形，每一次迭代選擇一個最亮的三角形當作光源，由這個光源計算與場景中其他三角形的 Radiosity 之值。其中，判斷光源是否可以輻射到某個三角形之複雜度介於 $O(\log N)$ 到 $O(N)$ (視Data structure而定)，而計算 Form-Factor的花費可以視為常數 $O(1)$ ，因此每次迭代的複雜度介於 $O(N \log N)$ 到 $O(N^2)$。

其中佔據效能的因素是 Form-Factor 估算，因此有像 Hemicube之類的近似逼近，大幅度減少計算量，但投影回到原本物件上會失真。

$$\begin{align*} F_{ij} = \frac{1}{A_i} \int_{A_i}\int_{A_j} \frac{\cos \phi_2 \cos \phi_1}{\pi r^2} dA_i dA_j \end{align*}$$

進入優化主題吧

優化技術 Code Review & Optimization

首先，我們先對助教提供的程式碼加速，分成以下幾個部分討論

  • 減少光線投射計算量 Strength Reduction for Ray Casting
  • 減少 Form-Factor計算量 Strength Reduction for Form-Factor
  • 改善資料局部性 Improve Data Locality
  • 其他優化 Other Optimization:
    • Improve I-cache Miss
    • Short Circuit Design
    • Clipping Algorithm
    • Strength Reduction for Float-Point
    • Shrink the Scope of Variables
    • Reduce the Number of Arguments
    • Remove Implications of Pointer Aliasing
    • Copy Optimization

減少光線投射計算量 Strength Reduction for Ray Casting

判斷射線 (Ray) 是否能打到三角形 $A$ 上，先用 bounding box 包住 $A$ ，計算 $p$ 到 bounding box 的時間 $t$ ，若 $t$ 大於目前的最小 $t_{\min}$ ，則退出。相反地，再計算更精準的 $t$ 。加入利用已知結果 $v = p + t_{\min} \cdot d, t_{\min} > 0$

1
2
3
4
5
6
7
8
9
10
int TriangleHitted(Vector p, Vector d, TrianglePtr tp, float *t_min) {
    float t = /* time t from p to bounding box of Triangle tp */;
    if (t < eps)
        return false;
    if (t >= *t_min)    /* important !! */
        return false;
    /* ... */
    *t_min = t;
    return true;
}

減少 FF 計算量 Strength Reduction for Form-Factor

根據公式 $F_{ij} = \frac{1}{A_i} \int_{A_i}\int_{A_j} \frac{\cos \phi_2 \cos \phi_1}{\pi r^2} dA_i dA_j$ ，一般我們的判斷順序會得如下：

1
2
3
4
5
6
7
8
9
float computeFormFactor(TrianglePtr srcTri, int logSrc, TrianglePtr desTri, int logDes, Vector p) {
    Vector dir = srcTri->c - p;
    float ff = 0;
    if (RayHitted(p, dir, logDes) == logDes) {
        float theta1 = CosTheta(dir, srcTri->normal), theta2 = CosTheta(dir, desTri->normal);
        ff = theta1 * theta2 * srcTri->area / (norm2(dir) * PI);
    }
    return max(ff, 0.f);
}

效能考量的因素：

  • RayHitted() 需要大量的計算
  • Form-Factor 在 float 儲存格式下可能無法得到貢獻，改採優先計算 Form-Factor 的值，再運行 RayHitted 判斷。調整加速了 2 倍多。
1
2
3
4
5
6
7
8
9
10
11
float computeFormFactor(TrianglePtr srcTri, int logSrc, TrianglePtr desTri, int logDes, Vector p)
    Vector dir = srcTri->c - p;
    float theta1 = CosTheta(dir, srcTri->normal), 
theta2 = CosTheta(dir, desTri->normal);
    float ff = theta1 * theta2;
    if (ff <= 0)    return 0.f;
    ff *= srcTri->area / (norm2(dir) * PI);
    if (ff <= 0)    return 0.f;
    if (RayHitted(p, dir, logDes) == logDes)
        return ff;
    return 0.f;

改善資料局部性 Improve Data Locality

程式碼中使用 3D-DDA Line Algorithm/Bresenham’s Line Algorithm 搭配 Octree，在找尋某個射線與哪個最近三角形相交。

  • 只需要儲存葉節點代表的立方體中，所有可能相交的三角形編號
  • 移除掉中間產生的編號，讓每一次 access 的 cache-miss 下降

在 3D-DDA 中，我們需要反查找空間中一點 p 在哪一個葉節點中，藉由固定的長寬高切割長度，可以在 $O(1)$ 時間內得知 [i][j][k] 各別的值。若限制大小，則建立陣列 [i][j][k] 查找。若自適應大小，則建立 hash 表查找，但根據實驗結果，效能並沒有改善，因為三角形個數過多導致命中機率過低。

對於 Static Tree 的 Memory Layout，大致上分成四種DFS Layout、Inorder Layout、BFS Layout、和 van Emde Boas Layout，目前程式使用的是 DFS Layout，這方面會影響到存取的效能。若有更多的時間，我們也可以測試平行處理的細粒度與這些 Memory Layout 的影響。

其他優化 Other Optimization

Improve I-cache Miss

壓縮程式碼長度以改善 I-cache miss，因為大部分的初始化只會運行一次，不應該交錯在時常被呼叫的函數之間，指令載入效能才會提高，同時也要做好 Code Layout，就能改善執行效能。

原始版本如下：

1
2
3
4
5
6
7
8
9
x, y = compute(0)
buildTree(x, y)
x, y = compute(1)
buildTree()
x, y = compute(2)
buildTree()
x, y = compute(3)
buildTree()
...

壓縮和整理後

1
2
3
for i from 0 to n
    x, y = compute(i)
    buildTree()

Short Circuit Design

判斷三角形與一個正交立方體是否相交，使用投影到二維空間中與一線段的交點是否皆存在。投影方法有 3 種 x-y, y-z, z-x，線段投影共有 2 種，共計 6 種情況。原先程式沒有做好短路設計，只要其中一種不符就應退出。

1
2
3
4
5
6
7
8
9
10
11
12
int CrossOver(TrianglePtr tri, 
            Vector g0, Vector g1) {
    for (xyz = 0; xyz < 3; xyz++) {
        // front face project
        if (!test())
            return false;
        // back face project
        if (!test())
            return false;
    }
    return true;
}

Clipping Algorithm

我們實作課程中所描述的 Cohen–Sutherland Algorithm 降低 branch 次數，使用 bitwise 操作引出 SSE (Streaming SIMD Extensions)。儘管 compiler -O2 替我們優化，為減少 stack push/pop 的次數，實作時請不要使用的 procedure call，否則會慢上許多。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
void computeRadiosity(TrianglePtr srcTri, TrianglePtr desTri, 
float ff[3]) {
    char mask1 = (p0[x] < g0[x])<<0 |
        (p0[x] > g1[x])<<1 |
        (p0[y] < g0[y])<<2 |
        (p0[y] > g1[y])<<3 ;
    char mask2 = (p1[x] < g0[x])<<0 |
        (p1[x] > g1[x])<<1 |
        (p1[y] < g0[y])<<2 |
        (p1[y] > g1[y])<<3 ;
    if (mask1&mask2)    
        return false;
    if (!(mask1|mask2))
    return true;
    // ... test

Strength Reduction for Float-Point

兩個外積結果相乘小於零，減少 instruction cycle 量，盡量用整數作為運算型態。

在現在的 Intel CPU 中，32-bit 浮點數運算基本上跟整數一樣快

1
2
3
4
5
6
7
8
float a = cross(/* */);
float b = cross(/* */);
if (a * b < 0)
    return false;
b = cross(/* */);
if (a * b < 0)
    return false;
...

事先判斷正負號，同時也防止溢位。

1
2
3
4
5
6
7
8
int a = cross(/* */) < 0;
int b = cross(/* */) < 0;
if (a != b)
    return false;
b = cross(/* */) < 0;
if (a != b)
    return false;
...

Shrink the Scope of Variables

減少變數生命週期的長度以增加放入暫存器的機會，而非 stack 上。

1
2
3
4
5
6
7
8
9
10
float rgb[3];
for (int i = 0; i < 3; i++)
    rgb[f(i)] = g(i);
/* ... */
if (maybe) {
    for (int i = 0; i < 3; i++) {
        rgb[h(i)] = g(i);
    }
    /* ... */
}

當邏輯很複雜時，編譯器不太能幫忙做分析，所以自己手動優化的效果會比較好，在 C/C++ 語言中，可以利用大括弧進行區域變數的設定。

1
2
3
4
5
6
7
8
9
10
11
12
13
{
    float rgb[3];
    for (int i = 0; i < 3; i++)
        rgb[f(i)] = g(i);
    /* ... */
}
if (maybe) {
    float rgb[3];
    for (int i = 0; i < 3; i++) {
        rgb[h(i)] = g(i);
    }
    /* ... */
}

Reduce the Number of Arguments

減少 stack push/pop 次數

1
2
3
4
5
6
struct Arg {
    int a0, a1;
};    // p1.a1 = p2.a0
int f(Arg p1, Arg p2) {
    /* ... */
}

如何修改合適的參數個數，必須看使用的機率和次數，才能達到最大效益。

1
2
3
4
5
6
struct Arg {
    int a0, a1, a2;
};
int f(Arg p1p2) {
    /* ... */
}

Remove Implications of Pointer Aliasing

移除指標 Aliasing，意指可能會指向相同記憶體位址，導致每次計算都要重新載入，不能放進暫存器中。如下述的寫法，編譯器無法判定 srcTri 是否與 desTri 相同，在累加時則重新載入 srcTri->deltaB[] 的數值，計算上可能會產生數次的 cache miss，隨著迴圈次數不斷突顯效能差異。

1
2
3
4
5
6
7
8
9
10
11
void computeRadiosity(TrianglePtr srcTri, TrianglePtr desTri,
 float ff[3]) {
    for (int v = 0; v < 3; v++) {    // vertex
        for (int c = 0; c < 3; c++) {    // color RGB
            flaot deltaB = desTri->Frgb[c]/255.0*RefRatio*srcTri->deltaB[c]*ff[v]/3;
            desTri->deltaB[c] += deltaB; 
            desTri->accB[v][c] += deltaB; 
            desTri->deltaAccB[v][c] += deltaB; 
        }
    }
}
  • 方法 1: 加入 if (srcTri != desTri) 判斷，讓編譯器在 Function Pass 階段著手的 Dependency Analysis 更好
  • 方法 2: 使用 Copy Optimization，同時把重複計算搬到 stack 上，或者使用 Polyhedal 表示法進行 Reordering Accesses for Loop Nest。這裡我們選擇前者，更容易引出 SSE
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
void computeRadiosity(TrianglePtr srcTri, TrianglePtr desTri, 
float ff[3]) {
    const float k = RefRatio / 255.0;
    float lo[3] = { desTri->Frgb[0]*k*(srctri->deltaB[0]), 
                    desTri->Frgb[1]*k*(srctri->deltaB[1]), 
                    desTri->Frgb[2]*k*(srctri->deltaB[2])};
    for (int v = 0; v < 3; v++)    {    // vertex
        for (int c = 0; c < 3; c++) {    // color RGB
            /* calculate the reflectiveness */
            float deltaB = lo[c] * ff[v] / 3;
            desTri->deltaB[c] += deltaB;
            desTri->accB[v][c] += deltaB;
            desTri->deltaaccB[v][c] = deltaB;
        }
    }
}

小結

  • 使用洽當的編譯器參數可加速 2 倍
  • 減少 Form-Factor 計算加速 2 倍
  • 剩餘優化部份改善 10% ~ 20% 效能。

至今，我們加速了 4 倍，在下一篇文章中，我們將繼續探討平行處理。

Experiment

Model Origin (sec.) Our v0.1 (sec.) Speedup
room.tri 10.27 4.70 2.18
hall.tri 176.92 38.50 4.59
church.tri 72.32 42.64 1.69
Read More +
手札日記

那段過往的蠢事

我一定誤會了什麼，才冒冒失失地進入這個世界，大家所期待的那個人明明不是我－《三月的獅子》

實驗室調侃

「什麼時候要簽博？」

學校年底寄了封向碩士一學年修畢、在學成績前 30% 且被教授認可有研究能力的學生，提供 逕行修讀博士學位 ，從那陣子之後，這句話不斷地從旁人口說出，可想而知是在玩弄我。大概是從某次開會中，老師三番兩次地說誰誰看起來適合簽博，又朝著我說能力怎樣怎樣地很安心，必然是學弟們的光芒尚未展現，目光才莫名其妙地投射在我身上。

有趣的是隔了幾週後，系上放寬到一學期修畢即可，心想這下子總算能反攻了，但論談話辯論技巧，一時還是贏不了。

「學弟們來簽博吧！」充滿嘲諷的反擊道
「可是老師只認同你」學弟緩緩地說道
「…」到底是怎麼想到那一塊去的，心中滿滿的困惑

順水推舟

「什麼 什麼」－ 《吹響吧！上低音號》

年底的事情仍尚未結束，在旁人的催促下，感覺需要紀錄點什麼－那段有趣的誤會。

跨年那一夜

在台北的第二次跨年，想去年在生測資完後，回家的路上看著一群攝影愛好者在圖書館後頭拍攝，而自己只能拖著疲累的身軀回家休息。2017 年稍微特別一點，在靠北工程師看到一句話「人家有閃光才叫跨年，我們這群單身狗不跨年，叫熬夜。」尾牙聚餐因學弟一句話衝動發出邀請後，果然還是一個人熬夜吧，這回是學弟們一起在實驗室熬夜，別有樂趣。

這一夜原本想約著一起寫作業，自己一個人做做還可以，嘗試各種修改，絞盡腦力讓自己無法去想那原本的計劃，差點就在 TLE 下跨年了。

「的確 成就感真是活下去所需要的能量啊」－《三月的獅子》

「楓之谷倒數計時 3, 2, 1, …」
「限時販售商店在哪？」

一群肥宅就這麼在楓之谷上跨年。這麼說也許突然，學弟們看我有玩楓之谷，一個個都來跟我一起玩，一次要養這麼多人相當不容易，也因此那陣子完全不能自由活動，一上線就是當工具人。

那晚走在回租屋處的路上，傳了封訊息說聲「新年快樂」這讓我不禁想起，曾經也做過類似的事情，對著四年沒聯絡的手機號碼，到了那天日子傳了封「祝你生日快樂」這麼說起來，就像動畫《秒速五釐米》心境，平凡地過著日子，卻時時刻刻想著那過往的蠢事。

跨年那一夜的清晨走回家，平常一起租屋的高中同學們都會在家，而今晚一個都沒有，房裡、客廳裡只有從路燈照來的光，摸著冰冷的樓梯扶手走回房裡，收到封「抱歉，整天沒都沒看訊息，現在在喝酒」

然而，不只是現實生活中，連遊戲中都感受到一種奇怪的氛圍，難道沒長時間在線上掛網，連網友都和這一夥人都想到一塊去了！素未謀面的網友，傳了些訊息給我。

「吾 你今年會脫魯」遞上作業的截圖說道
「吾 你今年X月X日前會脫魯」

心中滿滿地無奈，處處都有人在督促我。

新年那一日

三天連假必然要為了平行加速的期末專題奮戰，所以約了小夥伴一起努力。一早八點坐在床頭滑著手機，看看昨夜的廢文有什麼樣的反應，看著未讀訊息中一封也沒有，收拾著包包走到學校，在實驗室等著沒有約定的討論。看著什麼都沒進展的程序，不知不覺就一個人待到中午，心想那傢伙肯定宿醉了。

「抱歉，剛剛才醒。」
「宿醉了對吧？」
「沒宿！」

比起對於算法分析和優化技術，也許應該從更基本的算法開始講起。腦海中大量的資訊嘗試從平常不怎麼說話的口中說出，一時什麼都講不出來

「沒關係，你慢慢說」
「是這個意思吧？」描繪著投影片上的圖示說道
「抱歉，這講起來有點枯燥」
「沒關係，聽你講著講著，我開始對這個有點興趣。」

小組討論就這麼持續到晚上，學弟們還不跟我吃飯，不知道這群傢伙在想些什麼，你們這樣玩弄學長對嗎？就是要我們倆個一起吃。

心情複雜

一想到在自己喜歡的事物面前，思緒總沒辦法順暢且正常地表達出來！覺得自己遜斃了。

「遜斃了」－《風夏》

趕緊逃進期末專題程式編寫的修改，不再去理會旁人的總總，心想事情根本不是這樣子的，這種事情怎麼會落在我身上呢。

優化計算機圖學浮點數輸出轉文字檔部份，轉換 json 就可以更快，加速 4 倍之多。參考論文說道 16 倍加速，提供的程序 Warning 編譯訊息停不下來，修一修 Bug 送上 pull request，沒想到 Lemire 教授也在貢獻者裡，略感興奮。清晨就收到信，修改部份被 Merge！

一個小小的 CPU 中，突然多了一條程序需要執行，帶著這樣的心情去工作。沒想到這一天是老闆來，一整天全是全英開會，而且還有一場是咱的報告，事前都不知道啊！時間不斷消逝，終於迎來那一刻 … 原來已經沒救了，好丟臉啊。

「太丟臉了」－《風夏》

在當下心中不冒出某晚一起討論的閒聊

指著螢幕上一行一行的程序，試圖解釋自己的想法
「我改了這裡，好讓整個三角形更加地完整，不會破圖，不過也因此不太能平行 … 這有點麻煩」
在腦海中盤算著時程和算法
「如果要平行的話，要改成這個樣子，不過要花點時間，再弄個幾天吧」

「每個人都有自己不擅長的，不會每件事情都做得好，沒・關・係・的」
「我們已經做很多了呀，先把這些弄好吧」

我啊，對於自己不擅長的項目太過在意，對於做不好一件可以完成的事覺得相當愧疚，總是投入好幾倍的時間來完成，想想自己這麼疲憊到底是為了什麼？在追求什麼？聽到那段話，有如放下了大石，旁人要不什麼都不說，要不鼓勵盡量往前衝，距離上一次對我傳達「你可以休息一下」又是多麼久遠的事了。

不斷地敲打鍵盤，嘗試將自己的想法加入並進行測試，突然冒出了一句話
「欸欸，你爸有缺老婆嗎？」
「？」頓時落下手上的工作，腦海中喀喀聲瞬間即逝，『等等，這是什麼話題，現在我們在弄期末專題吧？』
「沒啦，我媽離婚一陣子，最近又交了另一半。」
「嗯？」對於這話題到底從哪裡來，心中滿滿的困惑
「她居然嗆我說『我比妳還有市場耶』，我是不是應該隨便交個男友？」
『咦？咦咦咦咦？』現在應該要沉住氣吧，擔心一不小心就要進警局
「沒事，你繼續寫吧。」

繼續敲著鍵盤的同時，心想好像錯過了什麼選項。

角色轉換

「只是一味地爬升，爬升，最終抵達的地方沒有回去的路」－《三月的獅子》

「今天學姊會來嗎？」學弟如此問道
「哪個學姊？你想說什麼」不耐煩地詢問

終於來到這一天，跟學弟談話時要把同屆同學用稱謂描述，叫全名也不對勁，只喊名也怪奇特的，若用學長學姊稱呼，一不小心談話又忘了是哪個學長學姊。

清醒之日

「我在你家，可以借用你房間嗎？」
「？」原來是高中同學來訪，如此地閒情逸致 … 到我的房間打電動！
「借你桌子打電動」

十一點多離開實驗室，一回到租屋的地方，從沒想過接下來才是那無法描述的黑歷史，為什麼連回家都要被酷刑拷問，這些不是我能決定的啊，學校那樣就算了，連好友都能從動態中發現個什麼，為什麼你們要處處肉搜呢？

次次日，收到通知要買電池回去弄熱水器，到處敲個房門請求工具支援，才發現自己有多孤單。不僅是電池沒電，連瓦斯都沒有，看來今天是冷水澡的清醒日！

「孤單的時候 就可以相互依靠」－《人渣本願》

根據以往的經驗，想認識的女生都會馬上有男朋友。我想這次也不意外，祝福各位，這個小劇場給大家帶來一陣子的歡笑。

學期的最後，授課老師說研究生只要有幾門課拿到 A+ 就好，研究方向就會比較明確。然而，在這當下把碩士課程八門都修完，回過頭來，發現到擅長的研究領域無法判定，對於研究能力更加迷惘。感謝小夥伴們的協助，全 A+ 成就取得。

「誰能溫暖我 誰能快點溫暖我」－《為美好的世界獻上祝福》

Read More +
解題區/出題解題

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

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

Task Description

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

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

Subtask

$N \le 100$

Input Format

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

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

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

Output Format

輸出一行最小距離和。

Sample Input 1

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

Sample Output 1

1
80

Solution

初學第一步－窮舉

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

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

踏出第一步－動態規劃

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

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

更進一步

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

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

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

計算基礎花費

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

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

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

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

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

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

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

優化一夜

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

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

優化二夜

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

優化三夜

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

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

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

優化四夜

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

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
#include <bits/stdc++.h>
using namespace std;
#define MAXNODE 5000
#define MAXCOL 32
#define MAXN 32
#define INF 0x3f3f3f
struct DancingNode {
    int left, right, up, down;
    int ch, rh, w;
} DL[MAXNODE];
struct HelperNode {
    int head, size, next, prev;
} HN[MAXNODE];
int help_head;
int s[MAXCOL];
int head, size, Ans;
void Remove(int c) {
    DL[DL[c].right].left = DL[c].left;
    DL[DL[c].left].right = DL[c].right;
    for (int i = DL[c].down; i != c; i = DL[i].down) {
    	if (HN[DL[i].rh].head == i)
    		HN[DL[i].rh].head = DL[i].left;
    	HN[DL[i].rh].size--;
        DL[DL[i].right].left = DL[i].left;
        DL[DL[i].left].right = DL[i].right;
        s[DL[i].ch]--;
    }
}
void Resume(int c) {
    for (int i = DL[c].down; i != c; i = DL[i].down) {
    	if (HN[DL[i].rh].head == i)
    		HN[DL[i].rh].head = DL[i].right;
    	HN[DL[i].rh].size++;
        DL[DL[i].right].left = i;
        DL[DL[i].left].right = i;
        s[DL[i].ch]++;
    }
    DL[DL[c].right].left = c;
    DL[DL[c].left].right = c;
}
void Reduce(int i) {
    int t = DL[i].rh;
    HN[HN[t].next].prev = HN[t].prev;
    HN[HN[t].prev].next = HN[t].next;
    s[DL[i].ch]--;
    
    DL[DL[i].down].up = DL[i].up;
    DL[DL[i].up].down = DL[i].down;
    for (int k = DL[i].right; k != i; k = DL[k].right) {
        DL[DL[k].down].up = DL[k].up;
        DL[DL[k].up].down = DL[k].down;
        s[DL[k].ch]--;
    }
}
void Recover(int i) {
    int t = DL[i].rh;
    HN[HN[t].next].prev = t;
    HN[HN[t].prev].next = t;
    s[DL[i].ch]++;
    
    DL[DL[i].down].up = i;
    DL[DL[i].up].down = i;
    for (int k = DL[i].right; k != i; k = DL[k].right) {
        DL[DL[k].down].up = k;
        DL[DL[k].up].down = k;
        s[DL[k].ch]++;
    }
}
void Select(int i, int j) {
    int s = DL[i].rh;
    Reduce(i);
    Remove(j);
    HN[HN[s].next].prev = HN[s].prev;
    HN[HN[s].prev].next = HN[s].next;
}
void Cancel(int i, int j) {
    int s = DL[i].rh;
    Resume(j);
    Recover(i);
    HN[HN[s].next].prev = s;
    HN[HN[s].prev].next = s;
}
int new_node(int up, int down, int left, int right) {
    DL[size].up = up, DL[size].down = down;
    DL[size].left = left, DL[size].right = right;
    DL[up].down = DL[down].up = DL[left].right = DL[right].left = size;
    return size++;
}
void new_row(int n, int R[][2], int rh) {
    int r, row = -1, k;
    int h = size;
    for (int i = 0; i < n; i++) {
        r = R[i][0];
        DL[size].ch = r, s[r]++;
        DL[size].rh = rh;
        DL[size].w = R[i][1];
        if (row == -1) {
            row = new_node(DL[DL[r].ch].up, DL[r].ch, size, size);
        } else {
            k = new_node(DL[DL[r].ch].up, DL[r].ch, DL[row].left, row);
        }
    }
    HN[rh].size = n;
    HN[rh].head = h;
    HN[rh].next = HN[help_head].next;
    HN[rh].prev = help_head;
    HN[HN[help_head].next].prev = rh;
    HN[help_head].next = rh;
}
void init(int m) {
    size = 0;
    help_head = 0, HN[help_head].next = HN[help_head].prev = help_head;
    head = new_node(0, 0, 0, 0);
    for (int i = 1; i <= m; i++) {
        new_node(i, i, DL[head].left, head);
        DL[i].ch = i, s[i] = 0;
        DL[i].rh = 0;	// important pointer (fail pointer)
    }
}
int Dx[MAXN], Dy[MAXN];
int markRStk[MAXNODE], markRidx = -1;
int markCStk[MAXNODE], markCidx = -1;
void printDL(int n) {
    for (int i = HN[help_head].prev; i != help_head; i = HN[i].prev) {
        int j = HN[i].head, prev = 1;
        printf("%2d ", DL[j].rh);
        do {
            while (prev < DL[j].ch)
                prev++, printf(" -- ");
            prev = DL[j].ch+1;
            int rh = DL[j].rh, ch = DL[j].ch;
            printf("[%2d]", DL[j].w - Dx[rh] - Dy[ch]);
            j = DL[j].right;
        } while (j != HN[i].head);
        while (prev <= n)
            prev++, printf(" -- ");
        puts("");
    }
    puts("----------------------");
}
int cH(int limit = INF) {
    int ins = 0;
    for (int i = HN[help_head].prev; i != help_head; i = HN[i].prev) {
        if (HN[i].size == 0)
            return INF;
        int j = HN[i].head, v = INF;
        int rh = DL[j].rh, ch;
        do {
            ch = DL[j].ch;
            v = min(v, DL[j].w - Dx[rh] - Dy[ch]);
            j = DL[j].right;
        } while (j != HN[i].head && v);
        if (v == INF || ins+v >= limit)
            return INF;
        if (v) {
            ins += v;
            Dx[rh] += v;
            markRStk[++markRidx] = v;
            markRStk[++markRidx] = rh;
        }
    }
/*
	printf("cH first\n");
	printDL(4);
*/
    for (int i = DL[head].right; i != head; i = DL[i].right) {
        if (DL[i].down == i)
            return INF;
        int j = DL[i].down, v = INF;
        int ch = DL[i].ch, rh;
        do {
            rh = DL[j].rh;
            v = min(v, DL[j].w - Dx[rh] - Dy[ch]);
            j = DL[j].down;
        } while (j != i && v);
        if (v == INF || ins+v >= limit)
            return INF;
        if (v) {
            ins += v;
            Dy[ch] += v;
            markCStk[++markCidx] = v;
            markCStk[++markCidx] = ch;
        }
    }
    return ins;
}
void rH(int backRidx, int backCidx) {
    while (markRidx > backRidx) {
        int a = markRStk[markRidx--];
        int b = markRStk[markRidx--];
        Dx[a] -= b;
    }
    while (markCidx > backCidx) {
        int a = markCStk[markCidx--];
        int b = markCStk[markCidx--];
        Dy[a] -= b;
    }
}
#ifdef SAME_CUT
int isSame(int vx, int vy) {
    int i = HN[vx].head, j = HN[vy].head;
    do {
        if (DL[i].ch != DL[j].ch && (DL[i].ch != vy && DL[j].ch != vx))
            return 0;
        int e1 = DL[i].w - Dx[vx] - Dy[DL[i].ch];
        int e2 = DL[j].w - Dx[vy] - Dy[DL[j].ch];
        if (e1 != e2)
            return 0;
        i = DL[i].right;
        j = DL[j].right;	
    } while (i != HN[vx].head && j != HN[vy].head);
    return i == HN[vx].head && j == HN[vy].head;
}
#endif
#ifdef MST_CUT
int parent[MAXN];
int findp(int x) {
    return x == parent[x] ? x : (parent[x]=findp(parent[x]));
}
int joint(int x, int y) {
    x = findp(x), y = findp(y);
    if (x == y) return 0;
    parent[y] = x;
    return 1;
}
int MST(int st, int ed, int limit = INF) {
    static int used[MAXN] = {}, cases = 0;
    if (HN[HN[help_head].prev].prev == help_head)
        return 0;
    int i = HN[help_head].prev;
    int n = 0;
    int ret = 0, in = INF, out = INF;
    cases++;
    do {
        int j = HN[i].head;
        do {
            int x = DL[j].rh, y = DL[j].ch;
            int w = DL[j].w - Dx[x] - Dy[y];
            if (x == st) {
                in = min(in, w);
            } else if (y == ed) {
                out = min(out, w);
            } else {
                if (used[x] != cases)
                    used[x] = cases, parent[x] = x, n++;
                if (used[y] != cases)
                    used[y] = cases, parent[y] = y, n++;
                n -= joint(x, y);
            }
            j = DL[j].right;
        } while (j != HN[i].head);
        i = HN[i].prev;
    } while (i != help_head);
    ret = in + out;
    if (ret >= limit)
        return limit;
    if (n > 1) {
        fprintf(stderr, "cc cut\n");
        return INF;
    }
    return ret;
}
#endif 
void DFS(int u, int lower_bound) {
#ifdef MST_CUT
    if (lower_bound + MST(u, 1, Ans - lower_bound) >= Ans) {
        return ;
    }
#else
    if (lower_bound >= Ans) {
        return ;
    }
#endif
    if (HN[HN[help_head].prev].prev == help_head) {
        Ans = lower_bound;
        return ;
    }
    
    int i = HN[u].head;
    vector< pair<int, int> > pick;
    do {
        int v = DL[i].ch;
        int runnable = 1;
        if (v == 1)
            runnable = 0;
#ifdef SAME_CUT
        if (runnable) {
            int e1 = DL[i].w - Dx[u] - Dy[v];
            for (int j = HN[u].head; j != i && runnable; j = DL[j].right) {
                int e2 = DL[j].w - Dx[u] - Dy[DL[j].ch];
                if (e1 != e2)
                    continue;
                if (isSame(v, DL[j].ch))
                    runnable = 0;
            }
        }
#endif
        if (runnable) {
            int backRidx = markRidx, backCidx = markCidx;
            int tu = HN[u].head;
            Select(tu, v);
            int c = cH(Ans - lower_bound);
            if (c != INF) {
                pick.push_back(make_pair(lower_bound + c + DL[i].w-Dx[u]-Dy[v], v));
            //	DFS(v, lower_bound + c + DL[i].w-Dx[u]-Dy[v]);
            }
            rH(backRidx, backCidx);
            Cancel(tu, v);
        }
        i = DL[i].right;
    } while (i != HN[u].head);
    sort(pick.begin(), pick.end());
    for (int i = 0; i < pick.size(); i++) {
        int v = pick[i].second;	
/*
		printf("%d -> %d, %d\n", u, v, pick[i].first);
		printDL(4);
*/
        int backRidx = markRidx, backCidx = markCidx;
        int tu = HN[u].head;
        Select(tu, v);
/*
		printf("Select\n");
		printDL(4);
*/
        int c = cH(Ans - lower_bound);
/*
		printf("Reduce\n");
		printDL(4);
*/
        if (c != INF)
            DFS(v, pick[i].first);
        rH(backRidx, backCidx);
        Cancel(tu, v);
    }
}
int main() {
    int G[MAXN][MAXN];
    int n;
    while (scanf("%d", &n) == 1) {
        init(n);
        for (int i = 1; i <= n; i++) {
            int R[MAXN][2], Rcnt = 0;
            for (int j = 1; j <= n; j++) {
                assert(scanf("%d", &G[i][j]) == 1);
                if (G[i][j])
                    R[Rcnt][0] = j, R[Rcnt][1] = G[i][j], Rcnt++;
            }
            new_row(Rcnt, R, i);
        }
        int backRidx = markRidx, backCidx = markCidx;
        Ans = INF;
        DFS(1, cH());
        rH(backRidx, backCidx);
        printf("%d\n", Ans);
    }
    return 0;
}
Read More +
手札日記

優化小插曲 下篇

2016・冬・續

學校課業

畢業門檻共計二十四學分，每堂課原則上都是三學分，意即要修八堂課，其中一部分可以選修外系碩班課，而我卻傻傻地全選資工碩班的課。大部分的人都是每學期修三到四門課，最後一學期專心做論文，前提指導教授派發的助教工作夠輕鬆，而不是每週都來亂的情況。

而這學期一開始有打算修英文課，開給碩班的英文課不外乎有很多系來修，對於我這種渣渣而言還稍微困難了些，但沒想到第一堂課的互動，在互動過程中回答問題時，讓我的人生增添了 黑歷史 ，隨後的分組成了邊緣人，於是我，退了。

計算機圖學

這裡只描述過程中的心情，實作細節有機會再發文

之前在大學時，課程中 GameMaker 這套軟體上弄遊戲就快崩潰，靠著小夥伴編寫劇本，自己兜弄 UI 設計與美工，討論如何做出適當的效果，反覆地測試運行流程，那時用了一整個學期，請參考 Github，那門課可說是我唯三拿到低於九十分的系上課程 (其中兩個是專題，也許是我嗆老師不懂設計、太過古板)。

第二個作業使用 Unity3D 開發，拉個 GUI 弄得昏天暗地。在物理引擎設定上，質量大小、推力反應、空氣阻力、碰撞動量 … 等，這些估計數據的感覺全還給老師，瞬間覺得自己是個蠢貨，怎麼調都不對。搞了幾天都弄不定，還好有小夥伴古古教我怎麼用，不然還真的會崩潰一陣子。

《田中君總是如此慵懶》

期末專案跟小夥伴古古一組，著手渲染技術之一的 Radiosity Github，用各式優化技巧和 WebGL、OpenGL 呈現作品，有機會我們再來談談大量的優化技術吧，將老師給的程序利用編譯器優化到極致！隨後再平行它！

演算法設計與分析

雖然沒修這門課，看著學弟去修這門課的作業，又知道助教是之前修平行的古耕竹同學，追求效能極限的對手！有對手就會成長，永不滿足！

其中一個作業是計算 0/1 背包的最佳解，課程中使用 branch-and-bound 算法，可想而知這種搜索法還是會還是會退化的。不管哪種寫法都是 NP-Complete，因此只能做常數優化，如何加速背包問題請參考《淺談背包問題 (0/1 Knapsack Problem) 優化那些事》那篇。花了好幾天加速他們的作業，用優化後的 DP 計算就能完爆他們早期的剪枝版本，在數據範圍小的時候，剪枝版本仍然比較好。

為了加速應用，使用了 DP 達到理想最佳解，為了加速這個 DP 計算，額外使用了 DP 優化計算，又為了優化 … 不斷地遞迴下去，最後達到貪心收斂？

「但我的愛和努力都還不成熟」－ 斯特拉的魔法

還有另一個 TSP 銷售員旅遊問題的實作優化，有空我再來跟大家說說吧。

優化心情一覽

摸遍上百人的 Code，還原所有學過的知識，偷學技術並變成自己的招式，Gate of Codeylon！

「Gate of Bootylon」

連續幾週都在高喊加速！你怎麼改，我就怎麼追。當成功加速時，就擺上切嗣大叔

「固有時制御 二倍數 Time Alter Double Accel」－ Fate/Zero

隔日發現第一名被奪走時，高喊奪回誓言「妹子被暴搜搶走！不行，我一定要攻下她。」

斯特拉的魔法

反反覆覆地，每天心情有如三溫暖般地洗禮，時時刻刻都無法鬆懈，我們越來越逼近真理，內心相當雀躍。古助教看著我與學弟們的競爭，也加入討論，感覺相當地不錯。

M <<「我在追求真理」
K >>「好好栽培學弟，學弟就是你的真理」
M <<「？」
K >>「真理還分性別嗎？」
M <<「…」

學弟們最後加速剪枝到優化後的 DP 無法追到。深感可惜，在沒有滿足鴿籠原理的情況下，我找不太到卡暴搜的測資。相反地，若很容易做到鴿籠原理，莫名其妙可以快個十幾倍。他們加速剪枝，我利用減少計算量加速，彼此都獲得更進一步的討論，有人討論果然可以走得更遠更快。

「想要走得更遠」－ 吹響吧！上低音號

最後這些題目被我們放在批改娘系統上，測資有卡實作邊界處理不好程序。雖然我不是最強，但在這第一學府就要成為最強，任何亂來的解法見一個斬一個，動用腦中所有的智慧生出最強的測資！啊，好像把別門課的作業弄壞掉了。當然作業只求計算正確，速度不要太離譜，這一點我們還是無法違背課程需求，要求每位同學與我們一起加速和探討似乎都錯了什麼。

實驗室雜談

教授語錄

「實作細節一點也不重要，根本沒人想聽」
「別太相信論文」
「理論都正確，但你做的有什麼用？」

紀錄心中的困惑

這一年，我越來越不能理解這奇怪的要求。當不想講細節的時候，又被問要怎麼做，細節大多都沒有乾淨的理論描述，完全依賴在計算機架構上，不同的架構又是另一套參數，講起來沒完沒了。

當被質疑報告的論文正確性時，偏偏那一篇還是 tier-1 conference (曾經最高的國際會議) 來的論文，到底是信還是不信呢？那我們被要求投那個國際會議時，又要秉持什麼樣的態度才合適。從研究論文開始，我也明白很多論文的部分成果有瑕疵。無妨地，我們可以理論推測那是很簡單的，但礙於篇幅有限。

當想做出一篇很有用的論文時，我覺得這困難到想休學不讀，根本沒有機會遇到有用的且可以寫出論文的機緣。所以當時跟老師談論到想休學，老師卻說論文隨便寫就發一篇，這樣隨便發一篇真的有用嗎？這麼說起來，我的人生好沒用呢。

「不是都怪你太弱嗎」－ 3 月的獅子

白色聖誕前奏

公司的大哥大姐打算在聖誕節辦活動，要來個交換禮物又是卡拉歌唱比賽，櫃檯姊姊和其他 AE 姊姊談論這事情時，恰好就在旁邊聽著。

A >>「打算辦交換禮物活動，妳覺得如何」
B >>「我覺得交換禮物不錯呀」
Morris 思考中
A >>「妳看看 Morris 好像不太行耶，主管 X 是打算換成卡拉歌唱比賽」
A >>「他們是不是不太願意花時間去挑禮物，他們工程師都這樣」
Morris 思考中
B >>「Morris 好像撐不住了」
A >>「比賽這種就不是每個人都能參與，有輸有贏的，好好的聖誕節為什麼要這樣？」

咱只是在思考這些活動的複雜度而已，行不行也不是我能決定的，喜歡估算的心錯了嗎？對於我來說負擔太大，不僅不會唱歌，連買禮物的品味都成問題。最後耗費一個週末只思考禮物要買些什麼。

「是從什麼時候開始的呢 會感覺聖誕好痛苦」－ 3 月的獅子

這學期在學校大大小小的事物都跟古同學有關，學弟常常要跟我描述是古學長 (演算法設計助教) 還是古學姊 (物聯網助教和計算機圖學一起修課)，同時他們兩位又有在我當平行助教的時候一起修平行課，稱呼方法搞得我好混亂，但跟這群有趣的傢伙同一屆不容易也很幸運。

跟學弟們討論到底要買什麼禮物作為交換，並抱怨著自己聖誕節魯到不行。此時，學弟突然說「你不是在追古學姊嗎？」當下的我愣了一下，學弟馬上補一刀「難道是追古學長嗎？」等等，這聽起來越來越不對勁，看來這兒已經沒我說話的餘地，學弟們的這套想法讓我不知所措。

「嗯 我會努力的」－ 斯特拉的魔法

週末買公司要求的禮物時，學弟學長們再三交代「不買給『她』嗎」？心想「她」到底是誰，在旁人的催促下，內心相當煎熬，我不懂的事情太多，對於沒有錯過的問題，完全不知如何下手解決，行動與常識近乎零，類似當機地直接跳過無法執行的部分。於是在聖誕的前一週，跑了數間店買了數個禮物準備。

不管送給誰，送不出去就在聖誕節自己拆，自己買的禮物自己拆，做好邊緣人的心理建設。同時，週末也順便補實驗室另一個女生生日蛋糕，請學弟騎車帶我去買蛋糕，畢竟去年她們替我慶生，這回有學弟們一起在實驗室，買蛋糕也不尷尬，這週末真是忙碌呀。

從此之後，學弟問任何話題，一提到「學姊」總要特別小心地回答，他們總是在刺探我的心意，無意也變有意。在一旁的實驗室「學姊」看著我慌忙地回答與再三確認不時在一旁偷笑。「也許，他們把你的潛意思說出來了。」在這陣子，我想我回答什麼都不是。

隨著聖誕逼近，學長學弟一直在催促我何時要送「很多人週一就開始送了，你怎麼還不送」我只能再三地反問「到底是想要我送給誰啦」想著適合自己品味的禮物，若要與眾不同，肯定是朝著送記憶卡內附程序，並留言到「這是我加速三倍快的程序，拜託了，請回禮追我」的那種。

「要我故意輸掉嗎？」－ 3 月的獅子

最後在聖誕前幾天，也就去公司上班的前一天把禮物送出去，這下你們就不會再撈叨我了吧？至於在聖誕節前一日，因為被要求上台唱中文歌而內心受創好幾天，有必要把聖誕節搞得如此繁複嗎？

跨年尾牙

實驗室尾牙總會問跨年打算怎麼過

M「反正沒人約，應該在實驗室過吧！」
A「學長，不用擔心，我們都在實驗室」
B「實驗室由我們守護，你放膽出去約吧，記得遊戲角色開著幫我放輪」
聽起來有點感動又有點頹廢是怎麼回事 …
M「約誰？」
B「當然是學姊」
… … …
M「你們認真的嗎？」
B「當然」
一群在旁邊偷笑
實驗室學姊對 M 說「你是認真的嗎？」
… …. … …. .. …
M「嗯，失敗了，跟你們過吧」

「明天早上能陪我一下嗎」－ 吹響吧！上低音號

2016・冬・終

Read More +
手札日記

優化小插曲 上篇

2016・冬・始

想想論文

第一階段

在學長的催促下，一時興起想個論文題目，加上所學的一些技巧來實作資料結構，再搭配資料加速的方式，就能在某些約束大小下灑亂數，沒有平行的版本就能比序列化版本快上一點，但平行版本並沒有想像中成正比地加速，一不小心就被 data layout 打到殘廢，被 imbalance 的問題搞死，沒弄好 scheduling 根本玩不下去。

論文可比想像中的困難多了，大多有用的技巧沒有理論，實務上的技術卻不被看好，對這個世界感到失望啊！

「那可比想像中的困難多了」－ 斯特拉的魔法

找了一篇別人寫的論文拿來平行，接著在平行處理根據現在的計算機架構設計。若沒有設計好資料結構和算法，效能遲遲無法引出來，有時還牽涉到實作成面的技巧，各方面都要顧及。

然而，實作總是不被看好的那一塊，大部分的論文都沒有提及或提供代碼，其他的論文到底可不可信？這一點深感疑惑啊！說不定早就錯得離譜，還一直誤導未來的讀者，當然我自己也不例外。至於如何解決這些問題，也許從一開始就是無解的題目呢。

「從一開始就是無解的題目呢」 － 果然我的青春戀愛喜劇搞錯了。

有一次大學長回實驗室，無聊問有沒有題目可以寫，於是拿了想要做的論文題目過去，三十分鐘就想到兩年前那篇論文怎麼做，再過三十分鐘的話我就不用做了，或許我們競賽時的題目，就是別人花費好幾個月所寫的論文結果。開始膽怯自己能力所寫的論文是多麼地無用。

儘管如此，著手實作論文再加上自己的算法設計，開始進行大量的優化探討，快還要再更快！如果寫不出教授滿意的論文，至少自己滿意就好。

第二階段

從沒想過會因為論文關連，再次翻起去年因可能沒畢業而著手的研究項目，降時間常數、降空間常數 … 等各種版本的實作。當初研究一點也不透徹，需求才是推動發展的主因！

「我真是半吊子吶」 － 3 月的獅子

研究中使用位元壓縮技術，將二元樹的記憶體配置弄得更好，又不能失去原本的操作性質，又要支持高效平行操作。折騰了好幾天，終於在常理測資範圍下，空間壓縮近 4 倍，平行加持下改善 25% 效能。回過頭來看看自己做了什麼，怎麼想都覺得很變態。

「啊！」－ 斯特拉的魔法

第三階段

因為需求一週就看了三四篇論文，連續過了好幾週，嘗試抽絲剝繭找到能用的工具，論文也許就是炒大鍋菜，但沒人告訴你加什麼料最好，更沒有人告訴你炒出來的那盤菜是否已經有人做過，做過就不能發論文了，一切只能靠自己，感覺起來就像在賭博，指導教授不見得有空理你呢。

實作技巧發揮到極致時，讓我們結合所有理論吧！一切皆為 $O(1)$ ，在平行上也 $O(1)$ 給你看。在十一月時，這一切發展得並不順利，根據理論推導的實作效能並不好，是不是有哪個傢伙的論文寫壞掉了？

「一個個都這樣 氣死我了」－ WWW.WORKING!! 迷糊餐廳

實習工作

工作機會

實習工作在完全不懂的電路設計自動化 Cadence 公司，每一個名詞基本上只有電機系畢業的學生懂，而我卻要在設計工具上撰寫客製化的功能，聽起來是多麼不可思議。實習的日子裡，同時在思考自己適合哪一類型的工作環境，是否能在這種大公司工作，還是小公司下工作快樂些。

那些沒上研究所的小夥伴都跑去找工作，而我卻在論文迷惘下找了份實習，體驗一下未來的工作和生活型態是否為所需。聽者小夥伴如此地抱怨：

「快把我變成社畜！不，算我求你了」友人如此哀求道
「這種需求我從來沒見過！」

某天早上六點收到一封信，看著英文的信封標題，想說一定是沒有分類成功的廣告信，社畜就這麼一如往常地去上班！等待編譯如此地漫長，打混之際來看封信 … “Google 面試邀請，下周可談”，基於各種因素，最後還是辭退了這次邀請，畢竟沒學歷沒英文，實在無法讓自己滿意地去面試工作。

「噓！」－ 吹響吧！上低音號

有些人覺得我浪費面試機會，但我足夠認識自己無法滿足他們的需求，現在還不是時候。你們才是更有機會的不是嗎？你們都擁有我想要的部分。

他們當時是藉由 Github 和 Blog 找到我的，努力經營自己的成果吧！總有一天，有人會因此欣賞你的。若沒有，就作為自己的最後一份驕傲吧！

工作環境

主管總是很有耐心跟我解釋顧客的需求和工具使用。身為一個工程師，對產品完全不熟悉怎麼跑測試，每一個元件之間的數學關係，多對一還是一對多，英文名詞聽得霧煞煞，有些還牽涉到多形和壓縮索引的格式。

相比從家教學生那兒聽來的，在一般公司就得自己摸，摸不好還要被主管嫌，而主管只會加油打氣，一點方向和意見都不給，甚至還導向錯誤的發展方向。仔細想想，這兒的工作環境挺不錯的。更沒想到公司在月底會替我們新進員工買生日蛋糕慶生，被排到十一月員工的慶生，心中很納悶明明過了快一個月，遲遲無法反應過來。

「生日禮物」遞上
「？？」
「來一下，來吃生日蛋糕」
「？？？」

「不管怎麼看 都是有人想找我茬」－ 吹響吧！上低音號

各個 Team 的成員和主管都有來，要求我們說些願望什麼的，對於咱這個肥宅而言，一夥人吃蛋糕的場景可說是零經驗，心中的小鹿都不知道消失到哪了。

2016・冬・上・結束

Read More +