手札日記

轉折序曲

2016・秋・始

雜務生活

上一屆的學長被拖到七月最後一天才口試結束，快八月底才送論文離開，實驗室雜務頓時落到身上來，心情非常複雜。離開的學長留了了一句話「不用擔心，你永遠有學長！」的確，當兵休學的學長們還不確定何時回來，目前還有數不出來幾屆的學長，前人所經歷的驚濤駭浪無法想像。

不過對我來說，在這一年內轉變好多，已經可以自己騎腳踏車去光華採購實驗室器材，從一個不出門買東西的肥宅是相當挑戰的。然而，決定要買什麼仍然非常非常困難的，目標求在各種約束條件下，能購入物品價值在未來可能所需的期望值最高，這問題在競賽中的約束相當多，用到現實生活中，連輸入都不給！

助教生涯結束

起初，因為沒有人願意接 PF 課程助教，實驗室成員首先發難接下，當下也不確定其他助教的素質，還是暫時性地接下來。就算是簡單的 C 語言，也不會每個大學畢業都很熟悉，連我也是半吊子的，要教台灣第一學府的莘莘學子們，沒有能力也要有責任學會，壓力並不小。

開學初助教名單尚未決定前，打開信箱十封中有八封都是助教工作相關的，沒意外地，第三輪噩夢即將開始，當然，助教這裡沒有必過法寶，寄信詢問是沒用的。

「我給你加油」－ 這個美術社大有問題！

當助教聽起來不錯，但自己的未來呢？心想論文連個底都沒有，覺得自己做什麼都是在逃避。按照普遍的流程來看，不少人題目決定好就陸續在剩餘兩學期慢慢做。

也許，很多朋友會說「去找指導教授討論看看」當然，我也這麼嘗試過，教授總是回「自己回去想想」我想研究是一個人的戰爭，不然就沒有研究能力吧？別以為老師總會給你方向。如果你遇到貴人給你指點，不妨多留點感恩的心吧。

「完全無法溝通」－ 這個美術社大有問題！

先不談論文內容要寫些什麼，對我來說用英文寫才是最大的難關，於是老師今年連課程小考題目描述都不寫，叫助教輪番用英文寫，每一週都要替學生寫小考題目，為了領區區的七千元，寫題目描述、生測資、提供解答一手包辦，處理學生雜務外，還得處理額外的系統維護和新功能加入。

用「學生本業是學習，不是賺錢！」安慰自己，花了兩三天寫好的英文描述，上陣前被 PF 全改掉，嫌我們描述複雜且沒邏輯順序，條列式描述全被改成好幾段的故事描述。我想學生們也很困惑，到底是來練英文的，還是來學寫程式的，條列式不是很容易編寫程式嗎？我自己認為，初學者在英文單詞對應邏輯判斷的 邊界條件 非常難適應。這些不是我追尋的理想，也不是我喜歡的交流學習方式，再加上有實習工作來找我，於是推辭掉了助教工作。

「我要推辭掉」－ 這個美術社大有問題！

推辭掉助教工作後，在職班博班學長向我問道「要不要去他們公司實習，目前論文也沒有方向吧？」仔細想想，換點環境做事也許能找到些什麼，於是帶著零基礎的英文去了外商公司實習，對方不嫌棄我這一部分讓我好感動。每週花個兩天到公司上班，其餘日子就在實驗室處理雜務接電話，教教學弟如何寫 C。

然而，去公司實習的過程並不順利，指導教授 PF 向我要求走學校正常管道，心裡這麼想

「什麼是正常管道呢？產學合作嗎？」
「透過產學合作計畫不知道要申請多久，而產學合作當作論文發展的一部分嗎？」
「去年看著產學合作的情況，又不給數據進行實驗，猜想數據情況寫論文，他們也不採用得到的理論」
「系辦轉發的徵學碩的實習生的公司徵人信都有，為什麼我就不能去外頭看看？」

老師只回了我一句「這可能是違法的，你還有領實驗室計畫的薪水」當下的我非常不能理解，這是嫌我工作時數不夠？別的實驗室研究生不一定有薪水可領，那我也沒有一定要領的使命。與其擔心我不領實驗室薪水會餓死，先想想一個月六千連吃飯都不夠，還不能外出打工的道理都不能理解，又不是在台北有家住。都來台北讀書了，那一點點錢的約束還不如自由點更省錢。

結果，老師擔心我助教時數不夠 「助教簽退時間不能造假」，更問我 「當初找你來做什麼的？」，畢竟系統大部分是我寫的，不然學弟們又要花一個學期從零開始撰寫系統。回顧當時問助教意願時，我只回了「能不當嗎？」心情更加地複雜，為什麼用這樣的口吻反問我這些？

「明明『活著』才是最痛苦的」－ 橘色奇蹟

最後，相信自己所做的事情沒有錯，向系辦申請助教離職手續。

指導學弟

雖說沒有接任助教一職，看著碩一學弟們弄小考題目、出測試資料到編寫程式風格，經驗和見解沒有預期中的多，於是夜間指導碩一學弟寫程式、操作環境，講著怎麼寫比較好、比較不容易出 bug、要怎麼抓 bug，多設計某個算法可以更快。回頭一想，是不是太過嘮叨了？

可是我好想從討論中挖掘更多更好的想法啊，在上研究所的這一年裡，能討論的對象大部分都不是現實生活身邊的人，跟著網友們學習。現在，請帶著愚蠢的我前進吧！

「真說起來我自己都不清楚自己對什麼感興趣」－ 斯特拉的魔法

本實驗室管理機器是最麻煩的，大部分的經歷中都不會有管理伺服器，出事解決的經驗更少，而我也是學長慢慢指導和自己投入大量時間才會一些。當然，剩餘這一年中要把伺服器交給學弟！

在某一夜中，看學弟順利裝了兩台 server，卻莫名其妙卡在第三台無法安裝，於是一個一個排除可能性。

  1. 發現無法讀取 USB 開機，重新燒錄後的結果，仍無法抓取，連燒錄軟體都讀不到。
  2. 拿 Live 檢查的 USB 隨身碟，BIOS 自動偵測 bootable 裝置一時還沒反應，開機在失敗好幾次才進入 Live，隨後又偵測不到。
  3. 拿新的 USB 隨身碟燒錄最新版 iso，停在特定檔案燒錄十幾分鐘未果。(看來是壞掉)
  4. 回頭拿舊的 USB 隨身碟燒相同的 iso，順利進入並成功安裝！

不玩啦！一直重複開機一個小時過去、兩個小時過去 … 壞了兩個 USB，開新的還有問題，人品差到不行。

「看吧 這個世界中哪裡都沒有你的容身之處」－ 3 月的獅子

與同學相聚

在十月的某個周末，我們來了個高中數學組的小小聚會。不知不覺經歷了五、六年，同學們各個厲害，有當警察在賺錢、目標當老師的、更有著決定往更高學府鑽研的，每個人說著自己的目標和生活的那份自信，讓我好羨慕！只有我找不到其他要做的事情，就只是寫覺得有趣的程式，卻不知道如何賺錢以及人生的下一步在哪。

「找不到其他要做的事情」－ 3 月的獅子

網友軼事

在 2016 秋這個時候上映動畫《你的名字》，靠小夥伴們訂票，這是我們圖片幕後團隊久違的聚會。在去看之前，那些情人氣氛早就在各處醞釀，而我走在人群中會被嘲笑吧？

「我就是覺得 自己在人群會被嘲笑」－ 斯特拉的魔法

平時我有玩線上遊戲〈楓之谷〉，有一位好友在國中時認識，在研究所的時候再次在網路上被搭話，才發現早在好幾年前就認識，談著談著，好像只有我記得那段日子的事情。常常會發現某些事情，我記得太過清楚，甚至連自己在娃娃車上誰推我去哪裡玩，買些什麼都記在腦海裡。

雖說是女網友，說特別也倒沒什麼，不過有男朋友的人還傳封訊息「走 我們一起去看」給我，當下的我還正在同一建築物參加表親的婚禮，而她還說著自己一個人去看電影的情況。後來還發現到，因碩二搬到學校附近住，她居然還住在附近，能推薦些吃的給我知道，生病時還能幫我找間診所。

機緣可說是巧妙呢

2016・秋・結束

Read More +
手札日記

轉折前奏曲

2016・夏・始

前言

距離上次寫日記已經超過六個月，大學時期還能一個月寫一篇，一到研究所就沒辦法隨心所欲發揮。這也收到不少朋友抱怨到部落格怎停了，其一原因在於寫一篇耗費數個小時，想做的事情太多，不想做的事情更多，在這些事情交替之下，寫一篇看似有內容的文章變得相當困難，也許這就是逐漸步入社會吧！怪不得有些人一進公司工作，接著就音訊全無，不久的我也會這樣吧，能繼續耍廢打混的日子不多了。

雖說沒辦法耗費心思寫文章，但還會在 Facebook 上發動態，不用擔心，接著就讓我整理一下這幾個月的動態吧！

想要追到

拿著在高等編譯器學到的知識，針對以前運行效果不佳的資料結構進行優化，有著更快更好的目標，想要追求到妳的那份心情無法抑制。研究算法與寫題目，就像追妹子一樣困難，也可以選擇兩方，但咱的目標不多，先優化下去追，追不到再嘗試開平行追！別說我犯規，有時候需要點不折手段，你說是吧？

七月時 —— 加速 Delaunay 三角剖分計算，實作平行版本的前置作業，優化弄得萌萌哒停不下來。

「現在的我 不想讓這個後悔留到 10 年後一直殘留」- 橘色奇蹟

學長畢業潮

才隔了一年，學長接二連三地畢業、休學、換實驗室，不像當初大學可以擁有三年的學長，來一年就要準備接管實驗室的一些雜事，當一屆只有一個男生時，管機器之類的顯得孤單。有一次一個人扛著 2U 伺服器回來，裡頭硬碟和四張顯卡都插滿笨重處理器，才發現一旁的警告標語要求兩人協力安裝，頓時才發現大家都要離開這兒，剩下我這一個不堪用的傢伙。

在這難熬的碩一升碩二暑假，還能每天看得到大二同學認真刷題，維護機器不能停，又不時想著自己的未來在哪，論文又要寫些、做些什麼。儘管如此，熬過一年也突破當初在這裡的誓言，早在一年前，休學與離開一事充斥在腦海裡，若沒有各位的支持，說實在難以在這環境下活著。

「謝謝你們 我現在又有勇氣了」－ ReLIFE

至強融核

想必每個人看到這中文翻譯名稱，都會覺得太酷了！

暑假都在忙著裝 Xeon Phi (中文翻譯: 至強融核)，參考各個網站提供的資訊，在一些奇奇怪怪的環境下安裝，弄了好幾週才開始寫在別於 GPU 的環境上寫平行，最後看到 200 個核心同時運行，內心相當激動。至於，這些研究對於論文有沒有幫助又是另一回事，過度期待是不好的，畢竟買來的 Xeon Phi 版本屬於 Knights Corner (騎士號角)，直到 Knights Landing (騎士登陸) 才是比較能抗衡 GPU 的版本。感謝各方的幫助，才完成環境設定，這中間充滿了 BUG 啊！

「各位看到、聽到、體驗到的 BUG 只不過是冰山一角，真的。」
「沒有 BUG 的改版就像我的女朋友，沒有存在過。」

「今天要向大神告白」－ ReLIFE

接著，當然要在還有學長的時候，嘗試把 Xeon Phi 弄上 OJ 上，讓我們好好地挑戰平行極限，看看不同的平行環境的極致效能。然而，評測方法不常被討論，這變得是只能由自己定義，除了這些外，運行模式也不好處理，因為屬於協處理器，也可以提供特殊的運行流程，心中總沒有一個底。至今，由於尚未需要實驗，又暫時擱淺一陣子，只有出出例題來玩玩，但編譯器本身就有所不同，能加入的暗示也不一定有相同行為，這些都是很困難的。

「就算遇到困難 只要有你在 我一定能努力活下去」－ Re: 從零開始的異世界

暑期外出

在這暑期的最後一波，小夥伴約我參加 COSCUP 2016，看到一堆 LLVM、一堆 js 以及非常非常多的 Pokemon Go，不少講者先讓附近站點灑花後才開始講，好讓大家不會太無聊。最大收穫是能在前端跑接近原生的 C/C++，會議提到的 emscripten 也許能來個前端 OJ，來個壯烈的分散架構吧！

這次會議應用周邊有一堆 hubot，這不得聯想先前亂搞的 AIML 和 live2D 初階應用，感覺有很多有趣的玩法。在大神們提到一堆網頁架構，需要將之前寫的整份砍掉才行，先破壞再創造，追求效能極致的道路。最後，活動收尾還是來個板橋高中桌遊團，在 Inker 的指引下，得知有一本《あなたの知らない超絶技巧プログラミングの世界》神書，讓你不只會 C，還會讓你成為魔法師！

(✪ω✪) － 這個美術社大有問題！

2016・夏・結束

Read More +
解題區/解題區 - 其他題目

2017 Facebook Hacker Cup Round 1

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

Facebook 2017 Hacker Cup Round 1

A. Pie Progress

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

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

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

B. Fighting the Zombies

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

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

C. Manic Moving

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

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

D. Beach Umbrellas

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

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

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

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

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

Solution A

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
#include <bits/stdc++.h>
using namespace std;
int main() {
    int testcase, cases = 0;
    scanf("%d", &testcase);
    while (testcase--) {
        int N, M;
        scanf("%d %d", &N, &M);
        int dp[305][305] = {};
        const int INF = 0x3f3f3f3f;
        for (int i = 0; i <= N; i++) {
            for (int j = 0; j <= N; j++)
                dp[i][j] = INF;
        }
        
        dp[0][0] = 0;
        for (int i = 0; i < N; i++) {
            int A[305];
            for (int j = 0; j < M; j++)
                scanf("%d", &A[j]);
            sort(A, A+M);
            for (int j = 0, sum = 0; j < M; j++) {
                sum += A[j];
                A[j] = sum;
            }
            for (int j = i; j < N; j++) {
                if (dp[i][j] == INF)
                    continue;
                for (int k = 0; k < M && k+j <= N; k++) {
                    dp[i+1][j+k+1] = min(dp[i+1][j+k+1], dp[i][j] + A[k] + (k+1)*(k+1));
                }
            }
            for (int j = i+1; j <= N; j++)
                dp[i+1][j] = min(dp[i+1][j], dp[i][j]);
        }
        printf("Case #%d: %d\n", ++cases, dp[N][N]);
    }
    return 0;
}

Solution B

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
#include <bits/stdc++.h>
using namespace std;
int main() {
    int testcase, cases = 0;
    scanf("%d", &testcase);
    while (testcase--) {
        int N, R;
        scanf("%d %d", &N, &R);
        vector< pair<int, int> > A;
        set<int> SX, SY;
        for (int i = 0; i < N; i++) {
            int x, y;
            scanf("%d %d", &x, &y);
            A.push_back(make_pair(x, y));
            SX.insert(x), SY.insert(y);
        }
        sort(A.begin(), A.end());
        
        int ret = 0;
        for (auto LX : SX) {
            for (auto LY : SY) {
                int cnt = 0;
                vector<int> used(N, 0);
                for (int i = 0; i < N; i++) {
                    if (A[i].first >= LX && A[i].first <= LX+R
                        && A[i].second >= LY && A[i].second <= LY+R)
                    cnt++, used[i] = 1;
                }
                for (auto TX : SX) {
                    for (auto TY: SY) {
                        int dd = 0;
                        for (int i = 0; i < N; i++) {
                            if (used[i])
                                continue;
                            if (A[i].first >= TX && A[i].first <= TX+R
                                && A[i].second >= TY && A[i].second <= TY+R)
                                dd++;
                        }
                        ret = max(ret, dd+cnt);
                    }
                }
            }
        }
        printf("Case #%d: %d\n", ++cases, ret);
    }
    return 0;
}

Solution C

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
#include <bits/stdc++.h>
using namespace std;
int main() {
    int testcase, cases = 0;
    scanf("%d", &testcase);
    while (testcase--) {
        int N, M, K;
        scanf("%d %d %d", &N, &M, &K);
        
        long long g[105][105] = {};
        const long long INF = 1LL<<60;
        for (int i = 1; i <= N; i++) {
            for (int j = 1; j <= N; j++)
                g[i][j] = INF;
            g[i][i] = 0;
        }
        for (int i = 0; i < M; i++) {
            int A, B;
            long long G;
            scanf("%d %d %lld", &A, &B, &G);
            g[A][B] = min(g[A][B], G);
            g[B][A] = min(g[B][A], G);
        }
        
        for (int k = 1; k <= N; k++) {
            for (int i = 1; i <= N; i++) {
                for (int j = 1; j <= N; j++)
                    g[i][j] = min(g[i][j], g[i][k]+g[k][j]);
            }
        }
        
//		for (int i = 1; i <= N; i++) {
//			for (int j = 1; j <= N; j++)
//				printf("%lld ", g[i][j]);
//			puts("");
//		}
        
        int S[5005], D[5005];
        for (int i = 0; i < K; i++)
            scanf("%d %d", &S[i], &D[i]);
        static long long dp[5005][105][2] = {};
        for (int i = 0; i <= K; i++) {
            for (int j = 0; j <= N; j++)
                dp[i][j][0] = INF, dp[i][j][1] = INF;;
        }
        dp[0][1][0] = 0;
        
        for (int i = 0; i < K; i++) {
            int s1 = S[i], d1 = D[i];
            for (int j = 1; j <= N; j++) {
                long long cc;
                cc = g[j][s1]+g[s1][d1];
                dp[i+1][d1][0] = min(dp[i+1][d1][0], dp[i][j][0] + cc);
                cc = g[j][s1];
                dp[i+1][s1][1] = min(dp[i+1][s1][1], dp[i][j][0] + cc);
                
                if (dp[i][j][1] != INF && i > 0) {
                    int sP = S[i-1], dP = D[i-1];
                    cc = g[j][s1]+g[s1][dP];
                    dp[i+1][dP][1] = min(dp[i+1][dP][1], dp[i][j][1] + cc);
                    cc = g[j][s1]+g[s1][dP]+g[dP][d1];
                    dp[i+1][d1][0] = min(dp[i+1][d1][0], dp[i][j][1] + cc);
                }
            }
            
//			for (int j = 1; j <= N; j++)
//				printf("%lld ", dp[i+1][j][0]);
//			puts("");
        }
        
        long long ret = -1;
        for (int i = 1; i <= N; i++) {
            if (dp[K][i][0] != INF)
                ret = max(ret, dp[K][i][0]);
        }
        printf("Case #%d: %lld\n", ++cases, ret);
    }
    return 0;
}

Solution D

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
#include <bits/stdc++.h>
using namespace std;
const long long MOD = 1000000007LL;
void exgcd(long long x, long long y, long long &g, 
	long long &a, long long &b) {
    if (y == 0)
        g = x, a = 1, b = 0;
    else
        exgcd(y, x%y, g, b, a), b -= (x/y) * a;
}
long long inverse(long long x, long long p) {
    long long g, b, r;
    exgcd(x, p, g, r, b);
    if (g < 0)	r = -r;
    return (r%p + p)%p;
}
int main() {
    int testcase, cases = 0;
    scanf("%d", &testcase);
    while (testcase--) {
        int N, M, R[2048], S = 0, mxR = 0;
        scanf("%d %d", &N, &M);
        for (int i = 0; i < N; i++)
            scanf("%d", &R[i]), S += R[i], mxR = max(mxR, R[i]);
        if (N == 1) {
            printf("Case #%d: %lld\n", ++cases, M);
            continue;
        }
        fprintf(stderr, "%d %d %d\n", N, S, M);
        long long invNplus = 1;
        map<long long, long long> C;
        {
            long long f = 1;
            for (int i = 1; i <= N; i++)
                f = (f * i)%MOD;
            invNplus = inverse(f, MOD);
            
            int l = 1, r = 1;
            f = 1;
            for (int i = M - 2*S; i <= M - 2*S + 2*mxR; i++) {
                int tM = i+N-1;
                if (tM < N)
                    continue;
                int L = tM-N+1, R = tM;
                if (r < L)
                    l = r = f = L;
                while (l < L)
                    f = (f*inverse(l, MOD))%MOD, l++;
                while (r < R)
                    r++, f = (f * r)%MOD;
                C[tM] = (f * invNplus)%MOD;
//				printf("C(%lld %d) = %lld, %lld\n", tM, N, C[tM], f);
            }
        }
        long long ret = 0;
        for (int i = 0; i < N; i++) {
            for (int j = 0; j < N; j++) {
                if (i == j)
                    continue;
                int tM = M + R[i] + R[j] - 2*S;
                if (tM+N-1 < N)
                    continue;
//				printf("add C(%d %d)\n", tM+N-1, N);
                ret += C[tM+N-1];
                ret %= MOD;
            }
        }
        
        long long f = 1;
        for (int i = 1; i <= N-2; i++)
            f = (f * i)%MOD;
        ret = ret * f;
        ret %= MOD;
        assert(ret >= 0);
        printf("Case #%d: %lld\n", ++cases, ret);
    }
    return 0;
}
Read More +
解題區/解題區 - 其他題目

2017 Facebook Hacker Cup 資格賽

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

Facebook Hacker Cup 2017 Qualification Round

A. Progress Pie

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

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

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

B. Lazy Loading

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

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

C. Fighting the Zombie

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

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

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

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

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

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

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

Solution A

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
#include <bits/stdc++.h>
using namespace std;
int main() {
    const double eps = 1e-8;
    const double pi = acos(-1);
    int testcase, cases = 0;
    scanf("%d", &testcase);
    while (testcase--) {
        int P, X, Y;
        scanf("%d %d %d", &P, &X, &Y);
        int ret = 0;
        
        if (X == 50 && Y == 50) {
            ret = 1;
        } else if (hypot(X-50, Y-50) > 50) {
        } else if (P == 100) {
            ret = 1;
        } else if (P == 0) {
            
        } else {
            int vx = X-50, vy = Y-50;
            int tx = 0, ty = 50;
            double theta = acos((vx*tx+vy*ty)/hypot(vx, vy)/hypot(tx, ty));
            if (tx*vy - ty*vx > 0)
                theta = 2*pi-theta;
            double t = (double) P/100.0*2*pi;
            if (theta <= t)
                ret = 1;
        }
        printf("Case #%d: %s\n", ++cases, ret ? "black" : "white");
    }
    return 0;
}

Solution B

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
#include <bits/stdc++.h>
using namespace std;
int main() {
    int testcase, cases = 0;
    scanf("%d", &testcase);
    while (testcase--) {
        int n;
        vector<int> A;
        scanf("%d", &n);
        for (int i = 0; i < n; i++) {
            int x;
            scanf("%d", &x);
            A.push_back(x);
        }
        sort(A.begin(), A.end());
        int ret = 0;
        int r = n-1, l = 0;
        while (r >= l) {
            int need = ((50 + A[r]-1) / A[r]);
            if (l + need-1 > r)
                break;
            l += need-1, r--;
            ret++;
        }
        printf("Case #%d: %d\n", ++cases, ret);
    }
    return 0;
}

Solution C

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
#include <bits/stdc++.h>
using namespace std;
int main() {
    int testcase, cases = 0;
    scanf("%d", &testcase);
    while (testcase--) {
        int H, S;
        scanf("%d %d", &H, &S);
        vector< pair<int, int> > A;
        double ret = 0;
        for (int i = 0; i < S; i++) {
            char s[128];
            scanf("%s", s);
            int X = 0, Y = 0, Z = 0;
            for (int j = 0, x = 0, sign = 1, idx = 0; s[j]; j++) {
                if (isdigit(s[j]))
                    x = x * 10 + s[j] - '0';
                if (s[j+1] == '\0' || !isdigit(s[j])) {
                    x = x * sign;
                    if (idx == 0)
                        X = x;
                    else if (idx == 1)
                        Y = x;
                    else
                        Z = x;
                    if (s[j] == '-')
                        sign = -1;
                    else
                        sign = 1;
                    x = 0, idx++;
                }
            }
            int l = X+Z, r = X*Y+Z;
            static const int OFF = 1024;
            static double dp[2][OFF];
            for (int i = 0; i < OFF; i++)
                dp[0][i] = dp[1][i] = 0;
            dp[0][0] = 1;
            for (int i = 0; i < X; i++) {
                int p = i&1, q = i&1^1;
                for (int j = 0; j < OFF; j++)
                    dp[q][j] = 0;
                for (int j = 0; j < OFF; j++) {
                    if (dp[p][j] <= 0)	continue;
                    for (int k = 1; k <= Y; k++)
                        dp[q][j+k] += dp[p][j]*((double) 1.f/Y);
                }
            }
            double sum = 0;
            for (int j = 0; j <= X*Y; j++) {
                if (j+Z >= H)
                    sum += dp[(X-1)&1^1][j];
            }
            ret = max(ret, sum);
        }
        printf("Case #%d: %.6lf\n",++cases, ret);
    }
    return 0;
}
Read More +
解題區/出題解題

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

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

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

輸入格式

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

輸出格式

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

範例輸入 1

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

範例輸出 1

1
7178

Solution

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

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

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

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

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

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

優化初夜

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

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
#include <bits/stdc++.h>
using namespace std;
namespace {
    static const int MAXW = 1000005;
    static const int MAXN = 1005;
    struct BB {
        int w, v, c;
        BB(int w = 0, int v = 0, int c = 0):
            w(w), v(v), c(c) {}
        bool operator<(const BB &x) const {
            return w * c  < x.w * x.c;
        }
    };
    static int run(BB A[], int dp[], int W, int N) {
        static int MQ[MAXW][2];
        for (int i = 0, sum = 0; i < N; i++) {
            int w = A[i].w, v = A[i].v, c = A[i].c;
            sum = min(sum + w*c, W);
            for (int j = 0; j < w; j++) {
                int l = 0, r = 0;
                MQ[l][0] = 0, MQ[l][1] = dp[j];
                for (int k = 1; k*w+j <= sum; k++) {
                    if (k - MQ[l][0] > c)	
                        l++;
                    int dpv = dp[k*w+j] - k*v;
                    while (l <= r && MQ[r][1] <= dpv)
                        r--;
                    r++;
                    MQ[r][0] = k, MQ[r][1] = dpv;
                    dp[k*w+j] = max(dp[k*w+j], MQ[l][1] + k*v);
                }
            }
        }
    }	
    static int knapsack(int C[][3], int N, int W) {
        vector<BB> A;
        for (int i = 0; i < N; i++) {
            int w = C[i][0], v = C[i][1], c = C[i][2];
            A.push_back(BB(w, v, c));
        }
        assert(N < MAXN);
        static int dp1[MAXW+1], dp2[MAXW+1];
        BB Ar[2][MAXN];
        int ArN[2] = {};
        memset(dp1, 0, sizeof(dp1[0])*(W+1));
        memset(dp2, 0, sizeof(dp2[0])*(W+1));
        sort(A.begin(), A.end());
        int sum[2] = {};
        for (int i = 0; i < N; i++) {
            int ch = sum[1] < sum[0];
            Ar[ch][ArN[ch]] = A[i];
            ArN[ch]++;
            sum[ch] = min(sum[ch] + A[i].w*A[i].c, W);
        }
        run(Ar[0], dp1, W, ArN[0]);
        run(Ar[1], dp2, W, ArN[1]);
        int ret = 0;
        for (int i = 0, j = W, mx = 0; i <= W; i++, j--) {
            mx = max(mx, dp2[i]);
            ret = max(ret, dp1[j] + mx);
        }
        return ret;
    }
}
int main() {
    int W, N;
    assert(scanf("%d %d", &W, &N) == 2);
    int C[MAXN][3];
    for (int i = 0; i < N; i++)
        assert(scanf("%d %d %d", &C[i][1], &C[i][0], &C[i][2]) == 3);
    printf("%d\n", knapsack(C, N, W));
    return 0;
}

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

優化二夜

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

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
static int run(BB A[], int dp[], int W, int N) {
    static int MQ[MAXW][2];
    for (int i = 0, sum = 0; i < N; i++) {
        int w = A[i].w, v = A[i].v, c = A[i].c;
        sum = min(sum + w*c, W);
        for (int j = 0; j < w; j++) {
            int l = 0, r = 0;
            MQ[l][0] = 0, MQ[l][1] = dp[j];
            for (int k = 1, tw = w+j, tv = v; tw <= sum && k <= c; k++, tw += w, tv += v) {
                int dpv = dp[tw] - tv;
                while (l <= r && MQ[r][1] <= dpv)
                    r--;
                r++;
                MQ[r][0] = k, MQ[r][1] = dpv;
                dp[tw] = max(dp[tw], MQ[l][1] + tv);
            }
            for (int k = c+1, tw = (c+1)*w+j, tv = (c+1)*v; tw <= sum; k++, tw += w, tv += v) {
                if (k - MQ[l][0] > c)   
                    l++;
                int dpv = dp[tw] - tv;
                while (l <= r && MQ[r][1] <= dpv)
                    r--;
                r++;
                MQ[r][0] = k, MQ[r][1] = dpv;
                dp[tw] = max(dp[tw], MQ[l][1] + tv);
            }
        }
    }
}

優化三夜

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

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

改善了 10% 效能

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
 static int run(BB A[], int dp[], int W, int N) {
    static int MQ[MAXW][2];
    for (int i = 0, sum = 0; i < N; i++) {
        int w = A[i].w, v = A[i].v, c = A[i].c;
        sum = min(sum + w*c, W);
        for (int j = 0; j < w; j++) {
            int l = 0, r = 0;
            MQ[l][0] = 0, MQ[l][1] = dp[j];
            int cache_max = MQ[l][1], cache_idx = MQ[l][0];
            int k_bound;
            k_bound = min((sum-j)/w, c);
            for (int k = 1, tw = w+j, tv = v; k <= k_bound; k++, tw += w, tv += v) {
                // tw = k*w+j, tv = k*v;
                int dpv = dp[tw] - tv;
                while (l <= r && MQ[r][1] <= dpv)
                    r--;
                r++;
                MQ[r][0] = k, MQ[r][1] = dpv;
                if (r == l)    cache_max = dpv, cache_idx = k;
                dp[tw] = max(dp[tw], cache_max + tv);
            }
            k_bound = (sum-j)/w;
            for (int k = c+1, tw = (c+1)*w+j, tv = (c+1)*v; k <= k_bound; k++, tw += w, tv += v) {
                int dpv = dp[tw] - tv;
                while (l <= r && MQ[r][1] <= dpv)
                    r--;
                r++;
                MQ[r][0] = k, MQ[r][1] = dpv;
                if (r == l)   
                    cache_max = dpv, cache_idx = k;
                else if (k - cache_idx > c)   
                    l++, cache_idx = MQ[l][0], cache_max = MQ[l][1];
                dp[tw] = max(dp[tw], cache_max + tv);
            }
        }
    }
}

優化四夜

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

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

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
#include <bits/stdc++.h>
using namespace std;
namespace {
    static const int MAXW = 1000005;
    static const int MAXN = 1005;
    static const int MAXC = 1<<12;
    struct BB {
        int w, v, c;
        BB(int w = 0, int v = 0, int c = 0):
            w(w), v(v), c(c) {}
        bool operator<(const BB &x) const {
            return w * c  < x.w * x.c;
        }
    };
    static bool cmpByWeight(BB x, BB y) {
        return x.w < y.w;
    }
    static int run(BB A[], int dp[], int W, int N) {
        static int MQ[MAXC][2];
        
        for (int i = 0, sum = 0; i < N; i++) {
            int w = A[i].w, v = A[i].v, c = A[i].c;
            assert(c < MAXC);
            sum = min(sum + w*c, W);
            if (c != 1) {
                for (int j = 0; j < w; j++) {
                    int l = 0, r = 0;
                    MQ[r][0] = 0, MQ[r][1] = dp[j];
                    int cache_max = MQ[r][1], cache_idx = MQ[r][0];
                    int k_bound;
                    r = (r+1)&(MAXC-1);
                    k_bound = min((sum-j)/w, c);
                    for (int k = 1, tw = w+j, tv = v; k <= k_bound; k++, tw += w, tv += v) {
                        // tw = k*w+j, tv = k*v;
                        int dpv = dp[tw] - tv;
                        while (l != r && MQ[(r-1+MAXC)&(MAXC-1)][1] <= dpv)
                            r = (r-1+MAXC)&(MAXC-1);
                        MQ[r][0] = k, MQ[r][1] = dpv;
                        if (l == r)	cache_max = dpv, cache_idx = k;
                        r = (r+1)&(MAXC-1);
                        dp[tw] = max(dp[tw], cache_max + tv);
                    }
                    k_bound = (sum-j)/w;
                    for (int k = c+1, tw = (c+1)*w+j, tv = (c+1)*v; k <= k_bound; k++, tw += w, tv += v) {
                        int dpv = dp[tw] - tv;
                        while (l != r && MQ[(r-1+MAXC)&(MAXC-1)][1] <= dpv)
                            r--;
                        MQ[r][0] = k, MQ[r][1] = dpv;
                        if (l == r)	cache_max = dpv, cache_idx = k;
                        else if (k - cache_idx > c)	
                            l = (l+1)&(MAXC-1), cache_idx = MQ[l][0], cache_max = MQ[l][1];
                        r = (r+1)&(MAXC-1);
                        dp[tw] = max(dp[tw], cache_max + tv);
                    }
                }
            } else if (c == 1) {
                for (int j = sum; j >= w; j--)
                    dp[j] = max(dp[j], dp[j-w]+v);
            }
        }
    }	
    static int greedy(int C[][3], int N, int W) {
        struct GB {
            int w, v, c;
            GB(int w = 0, int v = 0, int c = 0):
                w(w), v(v), c(c) {}
            bool operator<(const GB &x) const {
                if (v * x.w != x.v * w)
                    return v * x.w > x.v * w;
                return c > x.c;
            }
        };
        vector<GB> A;
        for (int i = 0; i < N; i++) {
            int w = C[i][0], v = C[i][1], c = C[i][2];
            A.push_back(GB(w, v, c));
        }
        sort(A.begin(), A.end());
        int ret = 0;
        for (int i = 0; i < N; i++) {
            int t = min(A[i].c, W/A[i].w);
            if (t == 0)
                return -1;
            W -= t*A[i].w;
            ret += t*A[i].v;
            if (W == 0)
                return ret;
        }
        return ret;
    }
    static int knapsack(int C[][3], int N, int W) {
        // filter
        {
            int filter = greedy(C, N, W);
            if (filter != -1)
                return filter;
        }
        vector<BB> A;
        for (int i = 0; i < N; i++) {
            int w = C[i][0], v = C[i][1], c = C[i][2];
            A.push_back(BB(w, v, c));
        }
        // reduce
        {
            sort(A.begin(), A.end(), cmpByWeight);
            map<pair<int, int>, int> R;
            for (int i = 0; i < N; i++)
                R[make_pair(A[i].w, A[i].v)] = i;
            for (int i = 0; i < N; i++) {
                int c = A[i].c;
                map<pair<int, int>, int>::iterator it;
                for (int k = 1; k <= c; k <<= 1) {
                    int w = A[i].w * k, v = A[i].v * k;
                    it = R.find(make_pair(w, v));
                    if (it != R.end() && i != it->second) {
                        int j = it->second;
                        A[j].c ++;
                        A[i].c -= k;
                    }
                    c -= k;
                }
                if (c > 0) {
                    int w = A[i].w * c, v = A[i].v * c;
                    it = R.find(make_pair(w, v));
                    if (it != R.end() && i != it->second) {
                        int j = it->second;
                        A[j].c ++;
                        A[i].c -= c;
                    }
                }
            }
        }
        static int dp1[MAXW+1], dp2[MAXW+1];
        BB Ar[2][MAXN];
        int ArN[2] = {};
        memset(dp1, 0, sizeof(dp1[0])*(W+1));
        memset(dp2, 0, sizeof(dp2[0])*(W+1));
        sort(A.begin(), A.end());
        int sum[2] = {};
        for (int i = 0; i < N; i++) {
            int ch = sum[1] < sum[0];
            Ar[ch][ArN[ch]] = A[i];
            ArN[ch]++;
            sum[ch] = min(sum[ch] + A[i].w*A[i].c, W);
        }
        run(Ar[0], dp1, W, ArN[0]);
        run(Ar[1], dp2, W, ArN[1]);
        int ret = 0;
        for (int i = 0, j = W, mx = 0; i <= W; i++, j--) {
            mx = max(mx, dp2[i]);
            ret = max(ret, dp1[j] + mx);
        }
        return ret;
    }
}
int main() {
    int W, N;
    assert(scanf("%d %d", &W, &N) == 2);
    int C[MAXN][3];
    for (int i = 0; i < N; i++)
        assert(scanf("%d %d %d", &C[i][1], &C[i][0], &C[i][2]) == 3);
    printf("%d\n", knapsack(C, N, W));
    return 0;
}
Read More +
解題區/出題解題

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

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

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

輸入格式

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

輸出格式

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

範例輸入 1

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

範例輸出 1

1
107

範例輸入 2

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

範例輸出 2

1
1735

Solution

Branch-and-bound

bound knapscak problem 古耕竹同學提供

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

想要更快嗎？多看論文且實作它吧！

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

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

參考代碼

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
#include <bits/stdc++.h>
using namespace std;
struct Item {
    int v, w;
    double cp() {
        return (double) v / w;
    }
};
int cmpCP(Item a, Item b) {
    double cpA = a.cp();
    double cpB = b.cp();
    if (cpA == cpB)
        return a.w > b.w;
    return cpA > cpB;
}
int n, m;
int ret = 0;
Item items[1005];
int prefixW[1005];
int prefixV[1005];
 
int h(int i, int w, int b, int &j, int &f) {
    int fw = i == 0 ? w : (w + prefixW[i-1]);
    for (j = max(b, i); j < n && prefixW[j] <= fw; j++);
    if (j >= n) {
        return prefixV[n-1];
    }
    int v = 0;
    f = 0;
    v += prefixV[j-1] - (i == 0 ? 0 : prefixV[i-1]);
    w -= prefixW[j-1] - (i == 0 ? 0 : prefixW[i-1]);
    if (w != 0) {
        f = 1;
        float k = (float) w / items[j].w;
        w -= k * items[j].w;
        v += k * items[j].v;
    }
    return v;
}
 
int greedy(int w) {
    int v = 0;
    for(int i = 0; i < n; i++)  {
        if (w >= items[i].w){
            w -= items[i].w;
            v += items[i].v;
        } else {
            return v;
        }
    }
    return v;
}
void dfs(int i, int w, int v, int b) {
    if (i >= n) {
        ret = max(ret, v);
        return;
    }
    int j = n, f;
    int hv = v + h(i, w, b, j, f);
    if (hv <= ret)
        return;
    if (f == 0) {
        ret = max(ret, hv);
        return;
    }
    if (w >= items[i].w)
        dfs(i + 1, w - items[i].w, v + items[i].v, j);
    dfs(i + 1, w, v, j);
    return;
}
 
int main() {
    while (scanf("%d %d", &m, &n) == 2) {
        for (int i = 0; i < n; i++)
            scanf("%d %d", &items[i].v, &items[i].w);
        stable_sort(items, items+n, cmpCP);
        for (int i = 0, w = 0, v = 0; i < n; i++) {
            w += items[i].w;
            v += items[i].v;
            prefixW[i] = w;
            prefixV[i] = v;
        }
        ret = greedy(m);
        dfs(0, m, 0, 0);
        printf("%d\n", ret);
    }
    return 0;
}

結論

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

Dynamic Programming

前言

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

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

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

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

優化初夜

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

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

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

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

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

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

優化二夜

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

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

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

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

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

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

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

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

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
#include <bits/stdc++.h>
using namespace std;
namespace {
    static const int MAXW = 5000005;
    static const int MAXN = 1005;
    static void run(int A[][2], int dp[], int W, int N) {
        for (int i = 0, sum = 0; i < N; i++) {
            int w = A[i][0], v = A[i][1];
            for (int j = min(W, sum + w); j >= w; j--)
                dp[j] = max(dp[j], dp[j-w]+v);
            sum += w;
        }
    }
    static int knapsack(int C[][2], int N, int W) {
        vector< pair<int, int> > A;
        for (int i = 0; i < N; i++)
            A.push_back(make_pair(C[i][0], C[i][1]));
        N = A.size();
        assert(N < MAXN);
        static int dp1[MAXW+1], dp2[MAXW+1];
        int Ar[2][MAXN][2], ArN[2] = {};
        memset(dp1, 0, sizeof(dp1[0])*(W+1));
        memset(dp2, 0, sizeof(dp2[0])*(W+1));
        sort(A.begin(), A.end());
        int sum[2] = {};
        for (int i = 0; i < N; i++) {
            int ch = sum[1] < sum[0];
            Ar[ch][ArN[ch]][0] = A[i].first;
            Ar[ch][ArN[ch]][1] = A[i].second;
            ArN[ch]++;
            sum[ch] += A[i].first;
        }
        run(Ar[0], dp1, W, ArN[0]);
        run(Ar[1], dp2, W, ArN[1]);
        int ret = 0;
        for (int i = 0, j = W, mx = 0; i <= W; i++, j--) {
            mx = max(mx, dp2[i]);
            ret = max(ret, dp1[j] + mx);
        }
        return ret;
    }
}
int main() {
    int W, N, C[MAXN][2];
    while (scanf("%d %d", &W, &N) == 2) {
        for (int i = 0; i < N; i++)
            assert(scanf("%d %d", &C[i][1], &C[i][0]) == 2);
        printf("%d\n", knapsack(C, N, W));
    }
    return 0;
}
Read More +
學校課程/平行程式

平行優化技巧－基礎篇

在撰寫平行程式時，仍會使用到編譯器的技巧，大幅度地減少指令個數，有時候這些優化甚至會把寫壞的分享記憶體 (shared memory) 存取給移除掉，如果能充分理解這些不好的設計，就不會造成實驗上發生的謬誤，還以為自己跑得很慢，平行根本沒加速到，這樣草草交出的報告，助教也收到煩了 (這裡不討論實驗方法本身就是錯的情況)。

內存篇

首先，來談談 cache line 是什麼吧？現在的機器都有著記憶體階層架構 (memory hierarchy)，按照大小排序，其中最小的就是 cache line，接下來才是 L1, L2, L3 cache，最後是最大的主記憶體 (main memory)，但某些超級電腦的設計就不是如此，而有些輔助處理器 (如 intel xeon phi) 缺少 L3 cache，而像 GPU 則沒有 cache 設計，但有提供有 bank 存取資料。這些都跟處理問題的類型有關，才導致有那些特殊設計。

Multicore

Cache Size

為什麼 cache line 很重要？因為存取速度最快，但能存放的資料也很小，通常是 64 bytes 或者 32 bytes。在快取設計中，一次把位址連續的資料一起搬靠近 CPU，這樣簡單的設計，導致我們在撰寫程式時，若能將結構定義得好便能讓使用率高，這非常仰賴編成者對於系統架構的了解。例如色彩 RGB，可以宣告成 struct RGB {flaot R, G, B}; 或者 float R[], G[], B[]，如果 RGB 三者會共同計算，例如色彩轉換通常是三者交互作用的表達式，因此前者構造方式較為妥當。

例題：批改娘 10087. Sparse Matrix Multiplication (OpenMP)

Cache Hierarchy

然而，有些情況充分使用反而不好，沒錯，就是平行程式，通常平行程式不外乎會有分享記憶體，這時候充分利用反而是一件壞事，因為我們需要讓數個核心 (core) 同時運行，因為每一個 core 都有各自的 cache line，若在不同 core 上同時修改相近的位址時，很容易產生 dirty bit，這導致要掉回 L3 cache (同一個 CPU，但分享在不同的 core 上) 進行同步。通常編譯器能做到 independent 分析，那麼就會利用暫存器配置來避開這問題，如果不行，就要靠自己手動來調整。

例題：批改娘 10085. Parallel Count (debug)

GPU

GPU 的記憶體設計分成四種，經常分成 on-chip 和 off-chip，差別在於是不是內建焊死，因此 on-chip 內建記憶體有比較好的效能，而 off-chip 則是外接記憶體，存取速度相對於 on-chip 慢個三倍以上。特別注意，有好就有壞，提供越快存取速度的記憶體能存放的容量越小，因此非常不容易配置。

在編寫 CUDA 或者是 OpenCL 時，提供 shared memory 的設計，如果請求數量過大，將會自動轉入 global memory，這一點沒有明確告知，只有在執行時期才會發現突然變慢。若採用存取速度較快的寫法，有時候不見得比較快，

例題：批改娘 10101. Fast Game of Life (CUDA)

Read More +
手札日記

走在平行道路上－後篇

到頭來，最好的平行就是不交流，最好的進展就是不溝通。－Antisocial Parallelism

崩壞季節

在學期後半段，「四月是你的謊言，五月是否要換題目，六月是否要重新來過」聽到這些是否能明白我在說什麼呢？沒錯，就是傳奇的碩二崩壞日子，若是以畢業導向來讀研究所，花個一年或一年半修課，剩餘時間還要產出論文，想而知壓力是如此地大，有些領域甚至要實作一個系統，那麼忙碌程度就明顯不同，學長們崩壞的日子越來越近。

平行助教

追求優化極致

課堂作業的平行題目大多都是純數值計算，相當枯燥乏味，但優化技術屬於編譯器範疇，對於平行課程而言，只需要知道概念就能交差，對於這一點我也無可奈何，數次詢問老師是否要以速度作為評分標準，得到的回答都是「不」，心中帶點傷心，很想鼓勵那些寫得不錯又願意花時間鑽研的同學。看到沒有鑽研的同學嫌這門課簡單、好混，真的好想拉著他們一起來研究。

與同學互動

相較於往年的平行課程，今天很多題目都被設計可以 Online Judge，但有批改娘後，同學在半夜用不錯的想法寫出較快的程序，剛睡醒的我馬上發現自己寫的程序不斷地被幹掉，又花了一個下午加上現有知識再追回去。反反覆覆地過著這種日子，每天醒來都會畏懼最好的寫法會使用什麼樣的概念。心裡過得很痛苦，即使如此我也挺開心的，同學不嫌棄廢廢的我球員兼裁判，還願意在我出的題目下花時間，在下非常感動。

「即使如此我也挺開心的」－《Re: 從零開始的異世界生活》

當我拚命優化 local memory 存取，卻在替同學 debug 時發現意外地加速，於是新境界到來，順便跟同學交流一下加速部份，甚至連開檔時間都要省！一起追尋神乎其技的感覺非常不賴。

效能進展 3571 ms (24-core CPU) -> 2567 ms (GPU, partial local memory) -> 2472 ms (GPU, full local memory) -> 1675 ms (GPU, full local memory + work group opt) -> 967 ms (GPU, global memory + I/O opt + embedded kernel code)

「來吧 快和我一起來修行吧」－《田中君總是如此慵懶》

「看著秒速被拉近，從 profiler 中看出當前問題卡在 memory bandwidth & usage，要著手這一塊優化了對吧。」

「就算是 1 μs，也省給你看！」

「又有人刷新速度極限，睡覺時收到訊息很想裝死。要是再強一點就好了，又要等價交換什麼了嗎？只剩下生命了啊」

「相互學習競爭，別的實驗室好厲害啊，一不留意效能就輸了，得快點想辦法。」

「遇到第一次執行時慢 800 ms 左右，接下來都執行相同輸入又不見這消失的時間，難道是所謂的 library cache miss ? 可是執行檔不到 1 MB 耶。」

「照理來說，演算法迭代到後期的收斂是可預測的，所以硬體能猜測的效能會加快，從平行的結果卻是越跑越慢嗎？明明一開始贏了，後來卻輸了，咱不明白啊，一定是寫壞掉了」

「一早起床就發現相當疲勞，一打開批改娘就發現排名被刷掉。如果贏不了，起碼不能輸，不然這樣有辱實驗室之名！M 之神啊，請賜予我力量活著。不然沒有信心活下去 …」

「動用腦中所有的智慧考慮這個問題」－《熊巫女》

事實上，這名跟我拚搏的學生只有一位 R04922075 古耕竹同學，就與 tmt514 談論的結果，一門課只需要有一位同學認真與助教一起學習就好，「一名足矣！」的精神已經深植我心。

後來才知道他之所以能在長時間拚搏都是因為住在實驗室裡，醒來就可以寫程式做實驗，累了就在實驗室睡，聽起來可是研究最高生活境界，當然聽起來是個廢人似的，但能省下龐大的住宿費相當吸引我，再加上在台北這種鬼地方，還得忍受每天一個多小時的通勤，先不談住宿有多昂貴，學校研究生宿舍不夠住真的不方便，不像以前中央大學，也許地理位置偏僻，研究生宿舍還一堆空位等著人去住呢！除非指導教授收不到學生，否則進住實驗室這件事情我還做不到啊。

「我做不到啊」－《田中君總是如此慵懶》

出題考驗

老師每週只要求同學們撰寫幾支小程序，套用平行有快就好，但時候發生平行反而慢，或者不如預期效果時，通常是因為平行撰寫後，導致編譯器的優化等級起不來，可是老師又不想讓同學學習優化，在兩頭難的情況下，只好進行自主研究，包裝成線上評測只是單純方便實驗，為了減少 IO 導致實驗不準確，可是費了相當大的心力在思考如何出題。

「你就隨便研究一下，然後自學吧！」－《田中君總是如此慵懶》

在沒有任何經驗傳承下，大部分的程式碼都被學長摧毀，拼湊零碎的知識弄成一題，再寫各種版本進行測試，調校參數進行實驗，簡直是在做大規模的研究似的，每天都如此反覆地度過好幾個星期。

「現在我的心裡已經裝不下你了」－《田中君總是如此慵懶》

思考如何解題、設計題目時，總會瞬間 CPU 全滿，誰也無法搶資源！這樣的日子過久了，覺得身心相當疲勞。

「現在的我需要休息」－《熊巫女》

到了學期後半，步入分散式計算的領域，但不知道能不能用單元測試的想法，建構出一套測試 Hadoop 題目的系統，雖然沒辦法完全模擬群集計算的環境，作為驗證程序正確性也有一定可靠性。採用 Hadoop Streaming 的方式進行 Judge。

曾經未完成的夢，是否在這裡能找到答案？

超艱難修課進行式

期中考

迎來在台大第二次期中考，遲遲無法適應這邊的考試方式，一部分原因也許是在以前學校大多都有考古題撐腰，而現在屬於無依無靠的情況下，又加上題目是一坨英文，看錯題目意思的機會太高，而在答題速度方面考試時間上會來不及，期中考變得一項大挑戰。其中一門課提供開書考，甚至允許同學上網搜尋答案，但是一打開試卷，突然冒出了一題「請看以下這一篇論文，然後回答下列問題」頓時的反射動作是直接跳過，經常看的英文論文通常是十幾頁在跳，考試時間兩個半小時，沒時間耗在這。

然而，在考試結束後，老師特別講到「作為一個研究生，看論文的能力很重要的」於是特地被老師點出來「你怎麼沒寫那一題呢？」只能苦笑地回應，在心中暗自地述說『因為我怕英文啊』那一題佔了 10% 成績，連看都沒看就放棄了，而那張試卷寫到最後一刻才把所有題目寫完，到底該不該慶幸沒去看那一題呢？

「我一點兒也不知道」－《熊巫女》

心靈打擊

這陣子受到強烈打擊，原因可以歸納出以下幾點：

  1. 大學長火速交了女朋友並且擺脫魯蛇行列
  2. 弄了半天的環境，被橫插一行突然 PASS
  3. 寫得程序在不知名的情況下被女同學說很醜
  4. 程序一直寫不好、不快
  5. 愚蠢至極

碩二學長則一直在暗示要抓交替，如報帳和財務管理、伺服器管理 … 等，準備進入三者皆無的狀態，徹底用到廢掉。

「你真的被詛咒了嗎？」－《Re: 從零開始的異世界生活》

新生活挑戰

由於助教薪水大概一個月七千，但發送方式不是學校支付，拖了好幾個月才會到，當然一個月七千在台北吃飯，若要吃點好的，一個月七千一點也不夠，都到了這把歲數，還是要想辦法養活自己的一小部分吧。於是小夥伴表示「薪水到現在都還沒發出來，一個月六千哪能過生活。」最後在小夥伴的催促下，Morris 打開家教網找零工啦。但打電話是多麼令人害怕得一件事情 …

「我絕對不會這麼失敗的」－《田中君總是如此慵懶》

被學長推去面試家教，明明都已經鼓起勇氣出發，卻又過程中潑冷水。去也不是，不去也不是，這搞得我好混亂，到底還能相信這世界什麼地方？

「妳這麼關心我，真不好意思，其實你不用太在乎我的，我喜歡一個人待著。」－《田中君總是如此慵懶》

白活系列

加速 NPC

「同樣都是人類，為什麼差距這麼大！」－《線上遊戲的老婆不可能是女生》

當年在 NCPC 搞不出來的 Problem I Christmas Gifts (NP-hard)，在賽後用 DLX 運行效果不錯，在啟發函數加上延遲標記更是屹立排名前數位已久。最近又因為平行把題目挖回來討論，在去年釣到大一學弟來解，便以飛快的速度擊破測資，最後達到加速 20x。再把當初需要跑 30 秒的測資來運行，現在只需要短短的 50ms。而都這麼快了 … 應該不用平行吧。

動態規劃

為了出題目，根據研究論文後的結果，出了一題常見的動態規劃題目，也就是 matrix chain multiplication，儘管算法最終可以加速到 $O(N^{\log_2 3})$，但在 30 年前就存在 matrix chain multiplication 的 $O(N \log N)$ 作法，那在 $O(N^3)$ 套了一堆剖分方便平行到底在窮忙什麼？

「當面對諸多選項時，我實在不曉得到底該選哪一個才不算錯。」當一個 Flag 打開之後，瞬間掉了不少 Performance …

近期活動

最近才把原本要出題目的測資補起來，論文誇張的加速似乎是跟非常非常蠢的版本比較，但快個兩三倍仍相當值得，編程複雜度也是增加兩三倍。各種層面的題目都已經全數實作好放上批改娘上提供線上評測。心想當研究一直發展下去，難道都不會窮盡？清單上的靶子長得越來越奇特，研究學者果真是一群變態啊。

「我已經失去活下去的信心」－《田中君總是如此慵懶》

學期最後一堂 Workshop 結束，心裡總算舒坦了些，不用每周在廢物的羞恥 play 下度過。感謝上課同學們賜教予我！

「謝謝」－《田中君總是如此慵懶》

Read More +
手札日記

走在平行道路上－前篇

前提摘要

這學期修課過得相當恐懼，期初上課就受到多方刺激，因為許多課第一堂都問說以前有沒有學過什麼，易想而知地，對學店逃出的我而言，只能默默在心中說「有聽過這個詞嗎？難道就沒辦法修課，這樣子根本沒有課程可以修。」幸好地，雖然沒有受過台大大學部課程的指導，大部分的內容以前都自己研究和玩過。

回頭看看，那些英文授課以及期末實作論文之類的，內心恐懼越來越深，而在下學期要報英文論文，學長因為連發音不正確或者用詞錯誤受到老師的深入指導，每一次開會都久戰一個多小時。而連基礎英文都只能靠 google 翻譯，為了確信實驗結果，都必須親自實作一番才行。腦子不斷地去喘測未來，估量這學期不間斷地抗戰的生活，然後還要受到老師嚴厲地批判。

接下來，就以圖文的方式回顧這一學期吧！

解題相關

逐漸地放下在 UVa Online Judge 的刷題日子，更沒有參與各種線上例行賽，由於曾經寫非常多的題目，仍然有不少人經常來詢問題目，不管是學校課程或者是 UVa 上的舊題目都會被拿來問。然而，有一些特殊案例是拿著世界總決賽等級的題目來問，一看就知道要寫個天昏地暗，說不定還不會的題目呢。

「來，請吧」-《為美好的世界獻上祝福》

只是比別人稍微努力一點，對於高難度的題目反應是相當慢的。當然，我還是盡力回答，不過那陣子還要趕作業，閒暇時間一點也不想再開題目，一開下去不知道會不會一整天就過去，這樣可就沒辦法好好寫作業，滿腦子都是揮之不去的題目。所以有忽略一些人的提問，在此向大家說聲抱歉。

實驗生活

燃燒經費

從碩一剛進來時，碩班學長已經擔任採購財務管理已經兩年，為什麼是兩年呢？這些就要留給本人來說，事實上這裡很多人都充滿兩年以上的回憶，若要成為這裡的一分子，負擔是非常沉重的。每學期都要消耗實驗室經費，經過半年的我仍然幫不上學長，沒有研究目標就不知道要怎麼樣的實驗環境，能提出的採購項目原則上都不會通過。

「對不起，我太沒用了」-《蒼之彼方的四重奏》

那一陣子，由於要架設實驗室的群集計算，花了好個星期都沒辦法把網路架設好，一部分原因都是因為想要在虛擬機器如 Xen 上面維護，這時候網路設定和應用程式之間的權限開放變得困難，這些牽涉到軟體設計，參數總是不如預期地運作。

為了解決管理介面每天都睡不好，每到實驗室看著旁邊吵到不行的伺服器，不想接觸那無法運行的廢鐵們。這時候就提案買個 安眠藥 之類的，可想而知地被老師打槍，因不久之後因為課堂要使用，建造不起來的壓力非常大，再加上實驗室已經沒人有架設經驗，對系統比較熟的蕭大帥還要忙著畢業論文，只好用毅力嘗試架設。這故事的最後，採用原生的方式完成，暫時先別為了防止系統掛掉而採用虛擬化技術保護。

巧遇英文

蕭大帥學長看我破爛的英文，設想推薦我學習英文的方法，例如去學校語言中心聽力練習，或者到圖書館借閱英文文法書籍，又或者參與大一英文課程，然而都因為時間點不對，而且修課過重而沒有動力參加，英文帶來的恐懼對我來說不是一時造成的。但是看英文字幕的動畫不是問題，可惜地歐美翻譯速度慢，再加上翻譯味道對不上，一部分是文化上的不同，這導致要找到對應詞彙困難，有機會再來採用英文字幕看新番吧，擷圖一定非常有趣。

「You were asleep for 15 years but still have a silver tongue」-《只有我不存在的街道》

路人搭訕

每天過著早上九點左右到實驗室，晚上十一點左右回住的地方，其一原因是要避免上班族群和補習班下課學生。每天過著像上班族往常的通勤生活。

在某次下雨的夜晚，看到窗外的雨停了，便拿著來滴著雨水的傘提早離開，快速地走在人行道上，突然一旁有聲細語傳來「你正要回家嗎？」轉過頭來看，原來是名剛打完工妹子，心想『還好不是什麼怪人，但主動跟我搭話也有一點怪人成份吧？』心中充滿地慌張回道「是的，剛從實驗室離開。」

「哦，你是研究生嗎？我剛從打工那裡下班。今天都沒有少收錢，超開心的，之前經常少收錢」帶著愉悅的語調回道
「辛苦你啦，你在活動中心打工嗎？」因為沒見過面的面孔，為了確定身分還是問些資訊來吧。
「嗯嗯，在麥當勞打工。」
『…』帶有疑惑地，可是我一點也沒印象，去過幾次麥當勞，但店員的年齡對不上的。
「不過我都在內場忙，所以可能沒有見過吧。」

聊著聊著都走到捷運站上，一路上只應答一些簡單的句子，突然來個人靠在身邊聊天，難免不受到驚嚇。到了捷運站，一般都直接走樓梯而沒有特地走到有電扶梯的入口，也許就因為要跟著我走

「可以借我攙著嗎？」她攙著我的胳膊走下去，這時候幼小的心靈受到沉重打擊。
「你有在練肌肉嗎？感覺是有肌肉的觸感 …」
「沒啦！」聽到內心不斷地尖叫、嘶吼，這一切都來得太突然。

「確實在啊，腦子不正常的孩子。」-《為美好的世界獻上祝福》

看來經常到實驗室打混的日子，導致逐漸地看到幻覺，總算可以達到擴充實境的地步，由衷地敬佩自己。

算法數學分析

經由學長的推薦，選修陳文進老師的「演算法數學分析」課程，這門課可說是各種演算法常見的數論技巧，在競賽中也非常容易見到，但從看過數學家是如何定義這些符號以及運算性質，這門課可說是增廣嚴謹數學定義與工具的好課。

然而，每週寫起作業來非常刺激，雖然水泥數學課本〈Concrete Mathematics〉大部分都有提供解答，但是解法就不太明確，有提供驗證答案正確性的參考，於是乎有各種神妙的解法，甚至不用套用課程所講的一些技法，老師也非常鼓勵同學用一些已知的知識來推論答案，有時還比教科書來得簡單扼要。

寫起作業每次大概花了一天到兩天，每週一個章節，老師沒講到的章節要自行閱讀，難度不會太高，花點時間坑一下基本上都能完成。作業推論則是經過一個星期的哀嚎，有時候六日推不出來，放著過好幾天才突然想到，數學的美妙近在眼前卻總是差一步就能推出。

「大家就愉快地化為塵土吧」-《為美好的世界獻上祝福》

由於這門課已經很久沒開，沒法找到合適的人選擔任助教，而常見的球員兼裁判的人選卻不選這門課，於是每週作業大家輪著改。有一次我改大家的作業，看到許多神奇的寫法，但也非常地痛苦，腦補了好幾十分鐘仍猜不出這傢伙在寫什麼，而我們寫程式的人，又經常失去一些數學常識，如 程式的等號和數學的等號是不相同的。若這門課程有助教，這助教一定會崩潰的！

奴工打雜

老闆又叫我把批改娘架設給 Data Structrue and Algorithm (簡稱 DSA) 課程使用，雖然跟他們實驗室要了一台 伺服器安裝，經過三番測試終於架設起來，只有稍微跟他們說明如何建置題目和加入使用者。畢竟是做免錢的，剩下就看他們造化。

沒想到最後還是採用老方法，繼承他們 DSA 去年的腳本改作業，畢竟要設置批改系統非常不容易。而上次交給電機系使用，他們成功地運作這到讓我感到意外，資工系反而沒有花時間架起來玩。

過了不久收到 DSA 課程來信要徵助教，想到已經擔任一門課助教，再一門會往生，有人還認為我會想去當 DSA 的助教呢？若是當了，早就來一場助教與學生之間的效能拚搏，想到若要協助架設批改系統，多麼令人感到戰慄。

《無彩限的幻影世界》

超艱難修課挑戰

清明連假都在寫作業，每天都看似好像有點進展，程序仍然沒辦法跑。覺得程式越寫越退化，退化到與組語奮戰。覺得各課程作業非常噁心，一年抵四年所學，都覺得快被後輩看起來認定是個廢物。如果寫不出來，是不是都是我的錯呢？

「不全都是你的錯嗎？」《好想大聲說出心底的話》

看著作業需要的 LLVM API，看不懂文件只好猜呀猜，針對文件窮舉各種英文姿勢，歷經崩潰的驗證思路後，才察覺流程果然有點詭異，確定描述上的瑕疵後，換個思路打掉重來吧！也許還需要寫點測試，才能知道到底有沒有寫對。於是作業又再次進入輪迴，進度從零開始！

平行助教

提前研究

在學期開始前，預想而知地會受到老闆指派當平行程式助教，從上學期 C 語言程式助教那時開始，平行設計的課堂題目就開始實驗，把平行題目放上批改娘系統上，沒想到大一新生也有嘗試去解決那些題目，台大新生的學習能力遠遠超出我的想像。

研究如何平行程式不算難，但是寫得好與壞差別很多，大部分的情況都會快上一些，慘的時候甚至會慢上數倍，經過一層層地解析，研究為什麼會變慢，慢是因為什麼因素所導致？是架構嗎？還是運行流程？又或者是算法重複計算？當遇到瓶頸的時候，只能埋頭苦幹地研究，深刻地體會到研究並不會提升智力，只是增加知識層面而已，像我這種笨蛋卻過著研究生活，怎麼想都奇怪呢！

「也不會提升智力了」-《為美好的世界獻上祝福》

優化挑戰

平行不只有平行，牽涉到硬體架構，隨著平行的需求所進行的算法設計，若整個算法對硬體友不友善，將會影響程式的快慢。最常見的就是快取問題，第一個遇到的是 Structure Of Array / Array Of Structure 的不同，針對使用層面，沒想到他們居然會差異到 10 ~ 20% 效能差異，也就是遠本跑 50 秒的程式，居然可以提升到 40 秒內跑完，若要發這種論文，想必需要知道很多硬體設計。

「Study Hard」－《無彩限的幻影世界》

製作批改系統

如果沒有正常執行完畢，莫名其妙伺服器的內存量增加，也找不到地方砍，重啟顯卡也沒辦法解決，眼睜睜地看著內存爆炸。OpenCL 要製作成 Judge 題目還有一段路要走，每次 Judge 就重新啟動 Server 是最慘的抉擇。

「這跟說好的不一樣啊？」－《Re: 從零開始的異世界生活》

大致上知道 Memory Leak 的出現點，不管是程序有錯或被強制中斷都來不及執行釋放，少釋放一個物件造成內存少 100MB，經過兩三百個案例就得重開機。暫時解決方案是寫腳本偵測記憶體用量，指定超過用量自動重開，一旦沒寫好將變成不斷地重開。

好不容易架設好批改系統，但是測速度有嚴重的啟動損耗，有時候在沙盒裡面運行還特別慢，當只裝 nVidia driver，啟動 OpenCL 基本消耗一秒多。為了使用 Profiler，意外地安裝 CUDA 相關插件，頓時啟動消耗變成毫秒等級。為了解決這奇怪的現象，還一度把伺服器搞壞，整台重灌才順利解決。

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

Read More +
學校課程/平行程式

批改娘 10116. Fast Dynamic Programming Computing I (OpenMP)

題目描述

給定一序列矩陣，期望求出相乘這些矩陣的最有效方法的乘法次數。

sequence.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
#include <stdio.h>
#define MAXN 4096
#define INF (1LL<<60)
int N;
long long dp[MAXN*MAXN], SZ[MAXN+5];
int main() {
    while (scanf("%d", &N) == 1) {
        for (int i = 0; i <= N; i++)
            scanf("%lld", &SZ[i]);
        for (int i = 1; i <= N; i++) {
            for (int j = 0; j < N-i; j++) {
                int l = j, r = j+i;
                long long local = INF;
                for (int k = l; k < r; k++) {
                    long long t = dp[l*N+k] + dp[(k+1)*N+r] + SZ[l] * SZ[k+1] * SZ[r+1];
                    if (t < local)
                        local = t;
                }
                dp[l*N+r] = local;
            }
        }
        printf("%lld\n", dp[0*N+N-1]);
    }
    return 0;
}

輸入格式

有多組測資，每組第一行會有一個整數 $N$ 表示矩陣鏈上有 $N$ 個矩陣，第二行上會有 $N+1$ 個整數 $Z_i$，表示矩陣鏈的每一個行列大小，例如當 $N = 3$ 時，輸入 10 30 5 60 表示矩陣 $A_{10, 30} B_{30, 5} C_{5, 60}$ 相乘。

  • $1 \le N \le 4096$
  • $1 \le Z_i \le 4096$

輸出格式

對於每組測資輸出一行一個整數 $M$ 為最少乘法次數。

範例輸入

1
2
3
4
5
6
7
8
9
10
11
12
13
14
2
2 2 2
 
3
10 30 5 60
 
3
1 5 20 1
 
3
5 10 20 35
 
6
30 35 15 5 10 20 25

範例輸出

1
2
3
4
5
8
4500
105
4500
15125

Solution

假設有 $p$ 個處理單元，矩陣大小為 $n \times n$，分析一般平行運算時間 $T_p$ 如下所示：

$$\begin{align*} T_p &= \sum\limits_{k=1}^{n} \left\lceil \frac{k}{p} \right\rceil \times (n - k) \end{align*}$$

針對地 $k$ 次迭代，將第二層迴圈平行化，每一個執行緒處理 $\left\lceil \frac{k}{p} \right\rceil$ 個狀態計算，每個狀態繼續需要 $n-k$ 次計算。

$$\begin{align*} \sum\limits_{k=1}^{n} \left\lceil \frac{k}{p} \right\rceil \times k &= 1 \times (1 + 2 + 3 + \cdots + p) + \\ & \qquad 2 \times ((p+1) + (p+2) + (p+3) + \cdots + 2p) + \cdots + \\ & \qquad (n \bmod{p}) \times \left\lceil \frac{n}{p}\right\rceil (\lfloor n/p \rfloor \times p + 1 + \cdots + n) \\ &= \sum_{k=1}^{\lfloor n/p \rfloor} k \cdot \frac{p (2k+1) p}{2} + (n \bmod{p}) \times \left\lceil \frac{n}{p}\right\rceil \frac{(n \bmod{p})(n + \left\lfloor \frac{n}{p}\right\rfloor p + 1)}{2} \\ &= p^2 \left[ \frac{\left\lfloor n/p \right\rfloor(\left\lfloor n/p \right\rfloor+1)(2 \left\lfloor n/p \right\rfloor + 1)}{6} + \frac{\left\lfloor n/p \right\rfloor(\left\lfloor n/p \right\rfloor+1)}{4} \right] + (n \bmod{p})^2 \left\lceil n/p \right\rceil \frac{n + \left\lfloor n/p \right\rfloor p + 1}{2} \\ \sum\limits_{k=1}^{n} \left\lceil \frac{k}{p} \right\rceil \times n &= n \cdot p \cdot \frac{\left\lfloor n/p \right\rfloor (\left\lfloor n/p \right\rfloor + 1)}{2} + n \cdot (n - p \cdot \left\lfloor n/p \right\rfloor) \left\lceil n/p \right\rceil \end{align*}$$

總結一下 $T_p = \sum\nolimits_{k=1}^{n} \left\lceil \frac{k}{p} \right\rceil \times n - \sum\nolimits_{k=1}^{n} \left\lceil \frac{k}{p} \right\rceil \times k = O(n^3 / p)$

針對前半段 $\sum\limits_{k=1}^{n} \left\lceil \frac{k}{p} \right\rceil$ 有以下兩種推法可供參考。

方法一

$$\begin{align*} \sum\limits_{k=1}^{n} \left\lceil \frac{k}{p} \right\rceil &= p \cdot 1 + p \cdot 2 + p \cdot 3 + \cdots + p \cdot \left\lfloor n/p \right\rfloor + (n \bmod{p}) \cdot \left\lceil n/p \right\rceil \phantom{0\over0}\\ &= \sum\limits_{k=1}^{\left\lfloor n/p \right\rfloor} k \cdot p + (n - p \cdot \left\lfloor n/p \right\rfloor) \left\lceil n/p \right\rceil \\ &= p \cdot \frac{\left\lfloor n/p \right\rfloor (\left\lfloor n/p \right\rfloor + 1)}{2} + (n - p \cdot \left\lfloor n/p \right\rfloor) \left\lceil n/p \right\rceil && \blacksquare \end{align*}$$

方法二

套用定理如下：(解釋： $n$ 個人分成 $p$，每組人數最多差一。)

$$\begin{align} n = \left\lceil \frac{n}{p} \right\rceil + \left\lceil \frac{n-1}{p} \right\rceil + \cdots + \left\lceil \frac{n-p+1}{p} \right\rceil \end{align}$$

套用上式從 $k = n$ 往回推。

$$\begin{align*} \sum\limits_{k=1}^{n} \left\lceil \frac{k}{p} \right\rceil &= \left\lceil \frac{n}{p} \right\rceil + \left\lceil \frac{n-1}{p} \right\rceil + \cdots + \left\lceil \frac{n-p+1}{p} \right\rceil + \sum_{k=1}^{n-p} \left\lceil \frac{k}{p} \right\rceil \\ &= n + (n - p) + (n - 2p) + \cdots + (n - \lfloor n/p \rfloor \cdot p) + \sum_{k=1}^{n \bmod{p}} \left\lceil \frac{k}{p} \right\rceil \\ &= \sum\limits_{k=0}^{\left\lfloor n/p \right\rfloor - 1} (n - k p) + (n \bmod{p}) \\ &= n \left\lfloor n / p \right\rfloor - \frac{\left\lfloor n/p \right\rfloor (\left\lfloor n/p \right\rfloor - 1)}{2} \times p + (n \bmod{p}) && \blacksquare \end{align*}$$

基礎平行

  • Accepted (50883 ms, 132224 KB)

平行地方呼之欲出，針對第二層迴圈直接平行化即可，下述代碼會帶入一點累贅，其原因在於做了各種實驗所致，但不影響正確性。只做了減少 thread 建立時的 overhead 的處理。

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 <stdio.h>
#include <omp.h>
#define MAXN 4096
#define INF (1LL<<60)
int N;
long long dp[MAXN*MAXN], SZ[MAXN+5];
int main() {
    while (scanf("%d", &N) == 1) {
        for (int i = 0; i <= N; i++)
            scanf("%lld", &SZ[i]);
        for (int i = 0; i < N; i++)
            dp[i*N+i] = 0;
        #pragma omp parallel firstprivate(N)
        {
            for (int i = 1; i <= N; i++) {
                #pragma omp for
                for (int j = 0; j < N-i; j++) {
                    int l = j, r = j+i;
                    long long local = INF;
                    dp[l*N+r] = INF;
                    for (int k = l; k < r; k++) {
                        long long t = dp[l*N+k] + dp[(k+1)*N+r] + SZ[l] * SZ[k+1] * SZ[r+1];
                        if (t < local)
                            local = t;
                    }
                    if (local < dp[l*N+r])
                        dp[l*N+r] = local;
                }
            }
        }
        printf("%lld\n", dp[0*N+N-1]);
    }
    return 0;
}

進階平行

  • Accepted (9890 ms, 136320 KB)

從一般平行實驗結果中，發現明顯地加速倍率不對，明明使用 24 個核心，跑起來不到兩倍加速的原因到底在哪？是的，在第三層存取順序需要 dp[l][k] 和 dp[k+1][r] 隨著 $k$ 變大，在前者存取順序是連續的，而在後者存取順序每一次跳躍長度為 $N$，導致 dp[k+1][r] 常常會發生 cache miss，一旦發生 cache miss，需要數百的 cycle 等待資料被帶進記憶體中。

由於矩陣鍊乘積計算時只用了矩陣的上三角，那複製一份到下三角矩陣去吧！dp[l][r] = dp[r][l]，如此一來在第三層迴圈發生 cache miss 的機會就大幅度下降。

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 <stdio.h>
#include <omp.h>
#define MAXN 4096
#define INF (1LL<<60)
int N, SZ[MAXN+5];
long long dp[MAXN*MAXN];
int main() {
    while (scanf("%d", &N) == 1) {
        for (int i = 0; i <= N; i++)
            scanf("%d", &SZ[i]);
        for (int i = 0; i < N; i++)
            dp[i*N+i] = 0;
        #pragma omp parallel firstprivate(N)
        {
            for (int i = 1; i <= N; i++) {
                #pragma omp for
                for (int j = 0; j < N-i; j++) {
                    int l = j, r = j+i;
                    long long local = INF, base = 1LL * SZ[l] * SZ[r+1];
                    long long *dp1 = dp + l*N, *dp2 = dp + r*N;
                    for (int k = l; k < r; k++) {
                        long long t = dp1[k] + dp2[k+1] + SZ[k+1] * base;
                        if (t < local)
                            local = t;
                    }
                    dp1[r] = dp2[l] = local;
                }
            }
        }
        printf("%lld\n", dp[0*N+N-1]);
    }
    return 0;
}

忘卻快取平行

  • Accepted (4264 ms, 134400 KB)

  • 參考論文 High-perofrmance Energy-efficient Recursive Dynamic Programming with Matrix-multiplication-like Flexible Kernels

最主要的精神 cache-oblivious algorithm 設計，利用大方陣切割成數個小方陣，矩陣跟矩陣之間在足夠小的情況下進行合併計算。算法概述如下：

cache-oblivious algorithm

每一個函數參數中的藍色區塊答案算出，而其他參數則是要合併的左矩陣和下矩陣，直到計算大小足夠小時，直接跑序列版本的程式，此時所有陣列可以全部帶進快取，這時候計算變得更加有利。再加上很多層的快取，每一個 CPU 的 cache 重複使用率相當高。

一般的平行度最高只有 $N$，因為我們只針對其中一個迴圈平行，而在 cache-oblivious algorithm 中，平行度最高為 $N^{1.415}$。這些都只是理論分析，實作時為了考慮硬體架構是沒辦法達到理想狀態的。

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
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <assert.h>
#define MAXN 4096
#define MAXD 6
#define MAXC 64
#define INF (1LL<<62)
int N;
long long dp[MAXN*MAXN], SZ[MAXN+5];
typedef struct MtxHead {
    long long *A;
    int bx, by, n;
} MtxHead;
MtxHead subMtx(MtxHead X, int i, int j) {
    MtxHead T = X;
    T.n >>= 1;
    if (i == 2)	T.bx += T.n;
    if (j == 2)	T.by += T.n;
    return T;
}
static inline int min(int x, int y) {
    return x < y ? x : y;	
}
static inline int max(int x, int y) {
    return x > y ? x : y;	
}
static inline long long loop(MtxHead X, MtxHead U, MtxHead V, int tl, int tr, int l, int r) {
    long long v = INF, t;
    long long comSZ = SZ[l] * SZ[r+1];
    for (int k = tl; k < tr; k++) {
        t = U.A[l*N+k] + V.A[(k+1)+r*N] + SZ[k+1] * comSZ;
        if (t < v)	v = t;
    }
    return v;
}
void Cloop(MtxHead X, MtxHead U, MtxHead V) {
    int n = X.n, l, r, tl, tr;
    long long v, t;
    for (int i = n-1; i >= 0; i--) {
        for (int j = 0; j < n; j++) {
            l = X.bx + i, r = X.by + j;
            v = X.A[l*N+r];
            
            tl = max(U.by, X.bx+i), tr = min(U.by+n, X.by+j);
            t = loop(X, U, V, tl, tr, l, r);
            if (t < v)	v = t;
            
            tl = max(V.bx, X.bx+i), tr = min(V.bx+n, X.by+j);
            t = loop(X, U, V, tl, tr, l, r);
            if (t < v)	v = t;
            
            X.A[l*N+r] = X.A[r*N+l] = v;
        }
    }
}
void Bloop(MtxHead X, MtxHead U, MtxHead V) {
    int n = X.n, l, r, tl, tr;
    long long v, t;
    for (int i = n-1; i >= 0; i--) {
        for (int j = 0; j < n; j++) {
            l = X.bx + i, r = X.by + j;
            v = X.A[l*N+r];
                        
            tl = max(U.by+i, X.bx+i), tr = min(U.by+n, X.by+j);
            t = loop(X, U, V, tl, tr, l, r);
            if (t < v)	v = t;
            
            tl = max(V.bx, X.bx+i), tr = min(V.bx+n, X.by+j);
            t = loop(X, U, V, tl, tr, l, r);
            if (t < v)	v = t;
            
            X.A[l*N+r] = X.A[r*N+l] = v;
        }
    }
}
void Aloop(MtxHead X) {
    int n = X.n;
    for (int i = 1; i <= n; i++) {
        for (int j = 0; j+i < n; j++) {
            int l = X.bx + j, r = X.by + j+i;
            long long v = X.A[l*N+r], t;
            long long comSZ = SZ[l] * SZ[r+1];
            for (int k = l; k < r; k++) {
                t = X.A[l*N+k] + X.A[(k+1)+r*N] + SZ[k+1] * comSZ;
                if (t < v)
                    v = t;
            }
            X.A[l*N+r] = X.A[r*N+l] = v;
        }
    }
}
void Cpar(MtxHead X, MtxHead U, MtxHead V, int dep) {
    if (X.n <= MAXC) {
        Cloop(X, U, V);
        return ;
    }
    #pragma omp task if (dep < MAXD)
    Cpar(subMtx(X, 1, 1), subMtx(U, 1, 1), subMtx(V, 1, 1), dep+1);
    #pragma omp task if (dep < MAXD)
    Cpar(subMtx(X, 1, 2), subMtx(U, 1, 1), subMtx(V, 1, 2), dep+1);
    #pragma omp task if (dep < MAXD)
    Cpar(subMtx(X, 2, 1), subMtx(U, 2, 1), subMtx(V, 1, 1), dep+1);
    #pragma omp task if (dep < MAXD)
    Cpar(subMtx(X, 2, 2), subMtx(U, 2, 1), subMtx(V, 1, 2), dep+1);
    #pragma omp taskwait
    //
    #pragma omp task if (dep < MAXD)
    Cpar(subMtx(X, 1, 1), subMtx(U, 1, 2), subMtx(V, 2, 1), dep+1);
    #pragma omp task if (dep < MAXD)
    Cpar(subMtx(X, 1, 2), subMtx(U, 1, 2), subMtx(V, 2, 2), dep+1);
    #pragma omp task if (dep < MAXD)
    Cpar(subMtx(X, 2, 1), subMtx(U, 2, 2), subMtx(V, 2, 1), dep+1);
    #pragma omp task if (dep < MAXD)
    Cpar(subMtx(X, 2, 2), subMtx(U, 2, 2), subMtx(V, 2, 2), dep+1);
    #pragma omp taskwait
}
void Bpar(MtxHead X, MtxHead U, MtxHead V, int dep) {
    if (X.n <= MAXC) {
        Bloop(X, U, V);
        return ;
    }
    Bpar(subMtx(X, 2, 1), subMtx(U, 2, 2), subMtx(V, 1, 1), dep+1);
    //
    #pragma omp task if (dep < MAXD)
    Cpar(subMtx(X, 1, 1), subMtx(U, 1, 2), subMtx(X, 2, 1), dep+1);
    #pragma omp task if (dep < MAXD)
    Cpar(subMtx(X, 2, 2), subMtx(X, 2, 1), subMtx(V, 1, 2), dep+1);
    #pragma omp taskwait
    //
    #pragma omp task if (dep < MAXD)
    Bpar(subMtx(X, 1, 1), subMtx(U, 1, 1), subMtx(V, 1, 1), dep+1);
    #pragma omp task if (dep < MAXD)
    Bpar(subMtx(X, 2, 2), subMtx(U, 2, 2), subMtx(V, 2, 2), dep+1);
    #pragma omp taskwait
    //
    Cpar(subMtx(X, 1, 2), subMtx(U, 1, 2), subMtx(X, 2, 2), dep+1);
    Cpar(subMtx(X, 1, 2), subMtx(X, 1, 1), subMtx(V, 1, 2), dep+1);
    Bpar(subMtx(X, 1, 2), subMtx(U, 1, 1), subMtx(V, 2, 2), dep+1);
}
void Apar(MtxHead X, int dep) {
    if (X.n <= MAXC) {
        Aloop(X);
        return ;
    }
    #pragma omp task if (dep < MAXD)
    Apar(subMtx(X, 1, 1), dep+1);
    #pragma omp task if (dep < MAXD)
    Apar(subMtx(X, 2, 2), dep+1);
    #pragma omp taskwait
    
    Bpar(subMtx(X, 1, 2), subMtx(X, 1, 1), subMtx(X, 2, 2), dep+1);
}
void Psmall(int N) {
    for (int i = 0; i < N; i++)
        dp[i*N+i] = 0;
    #pragma omp parallel
    {
        for (int i = N-1; i > 0; i--) {
            int comN = N-i;
            #pragma omp for
            for (int j = 0; j < i; j++) {
                int l = j, r = comN+j;
                long long local = INF;
                long long *dp1 = dp+l*N, *dp2 = dp+r*N;
                long long comSZ = SZ[l] * SZ[r+1];
                for (int k = l; k < r; k++) {
                    long long t = dp1[k] + dp2[(k+1)] + comSZ * SZ[k+1];
                    if (t < local)
                        local = t;
                }
                dp[l*N+r] = dp[r*N+l] = local;
            }
        }
    }
    printf("%lld\n", dp[0*N+N-1]);
}
int main() {
    while (scanf("%d", &N) == 1) {
        for (int i = 0; i <= N; i++)
            scanf("%lld", &SZ[i]);
        if (N <= 2048) {
        	Psmall(N);
        	continue;
        }
        int ON = N;
        while ((N&(-N)) != N)
            N++;
        for (int i = 0; i < N; i++) {
        	for (int j = i+1; j < N; j++)
        		dp[i*N+j] = INF;
            dp[i*N+i] = 0;
        }
        MtxHead X;
        X.n = N, X.bx = X.by = 0, X.A = dp;
        #pragma omp parallel
        {
        	#pragma omp single
        	Apar(X, 0);
    	}
        printf("%lld\n", dp[0*N+ON-1]);
    }
    return 0;
}
Read More +