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

此為計算機圖學課程中的一環,自由挑選主題。而這個 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 +

轉折序曲

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 +