解題區/解題區 - UVa

UVa 13032 - Marbles in Jars

Problem

有 $N$ 個罐子，其中有一個罐子放置彈珠時，每一個彈珠重量會是其他罐子的 1.1 倍，也就是說其中有一個魔法罐子，在魔法罐子中的彈珠都會特別重。現在限制每一個罐子放置彈珠的數量，請問有多少放置彈珠方案，可以找到魔法罐子。

Sample Input

1
2
3
4
5
6
7
3
3
1 2 3
3
1 2 2
2
2 2

Sample Output

1
2
3
Case 1: 1
Case 2: 0
Case 3: 2

Solution

從題目描述中理解到，要找到魔法罐子最簡單的策略就是每一個罐子都擁有不同數量的彈珠，由於每一個罐子限制彈珠數量不一致，做起來就稍微複雜。

假設先依照可仿製彈珠數量多寡排序，將放置少量的彈珠優先放置在前面，依序討論放入新的罐子又與之前放置不同個數的彈珠的方法數。

定義 dp[i][j] 為前 $i$ 個罐子，其中有一罐放置最多的彈珠 $j$ 個。由於已經由小到大排序好，放入新的罐子，分成兩種情況考慮，其中一種是比 $i$ 個罐子放置更多的彈珠，不然就是放置較少的彈珠數量，數學組合一下即可。時間複雜度 $\mathcal{O}(N M^2)$

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
#include <bits/stdc++.h>
using namespace std;
const long long MOD = 1000000007;
int main() {
    int testcase, cases = 0;
    scanf("%d", &testcase);
    while (testcase--) {
        int N, M[128];
        scanf("%d", &N);
        for (int i = 1; i <= N; i++)
            scanf("%d", &M[i]);
        sort(M+1, M+1+N);
        
        long long dp[128][128] = {};
        dp[0][0] = 1;
        for (int i = 1; i <= N; i++) {
            for (int j = 0; j <= M[i]; j++) {
                if (j - (i-1) >= 0) {
                    dp[i][j] = dp[i][j] + dp[i-1][j] * (j - (i-1));
                    dp[i][j] %= MOD;
                }
                for (int k = j+1; k <= M[i]; k++) {
                    dp[i][k] = dp[i][k] + dp[i-1][j];
                    dp[i][k] %= MOD;
                }
            }
        }
        
        long long ret = 0;
        for (int i = 0; i <= M[N]; i++)
            ret = (ret + dp[N][i])%MOD;
        printf("Case %d: %lld\n", ++cases, ret);
    }
    return 0;
}
Read More +
解題區/解題區 - UVa

UVa 13038 - Directed Forest

Problem

給予一張有向森林圖，要求任兩點若有一條路徑，則兩點不可著相同顏色，請問最少著色數為何？

Sample Input

1
2
3
4
5
6
7
8
9
10
11
12
3
4  2
1  2
3  4
4  3
1  2
2  3
4  1
7  3
3  6
3  7
1  5

Sample Output

1
2
3
Case 1: 2
Case 2: 4
Case 3: 2

Solution

明顯地，最長路徑長度若為 $L$，至少需要 $L$ 個顏色，因此針對每一個樹根進行最長路徑搜索即可，又由於題目給定森林，把每一個樹根抓出來取最大值即可。複雜度 $\mathcal{O}(N)$

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
#include <bits/stdc++.h>
using namespace std;
const int MAXN = 131072;
vector<int> g[MAXN];
int indeg[MAXN];
int dfs(int u) {
    int ret = 1;
    for (auto &v : g[u])
        ret = max(ret, dfs(v)+1);
    return ret;
}
int solve(int root, int N) {
    return dfs(root);
}
int main() {
    int testcase, cases = 0;
    int N, M, u, v;
    scanf("%d", &testcase);
    while (testcase--) {
        scanf("%d %d", &N, &M);
        for (int i = 0; i <= N; i++)
            g[i].clear(), indeg[i] = 0;
        for (int i = 0; i < M; i++) {
            scanf("%d %d", &u, &v);
            g[u].push_back(v), indeg[v]++;
        }
        
        int ret = 1;
        for (int i = 1; i <= N; i++) {
            if (indeg[i] == 0)
                ret = max(solve(i, N), ret);
        }
        printf("Case %d: %d\n", ++cases, ret);
    }
    return 0;
}
Read More +
解題區/解題區 - UVa

UVa 13005 - Blood groups

Problem

還記得孟德爾遺傳法則嗎？從血型分成 A, B, AB, O 型四種，由三個遺傳因子 A, B, i 兩兩組合而成。而外星人突破孟德爾遺傳法則，由 $N$ 個遺傳因子 (不包含隱性因子 i) 決定表徵，並且一個子代會從 $N$ 個父代分別繼承一個遺傳因子。

現在給予在場 $N$ 個父代分別有哪些遺傳因子 (若沒有告知，則剩餘都是隱性因子 i)，接下來會有 $Q$ 組詢問，問某一個遺傳因子組合是否可以被在場 $N$ 個父代組合而成。

Sample Input

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

Sample Output

1
2
3
4
5
6
7
8
N
Y
N
Y
N
Y
Y
N

Solution

每一個父代都提供一個遺傳因子，問最後能不能滿足所有匹配，顯而易見地是一題完美二分匹配。建造 source - 父代 - 遺傳因子 - sink。若某一組詢問需要父代提供遺傳因子 $x$，就去找尋有哪些父代可以提供遺傳因子，並且拉一條邊 父代 - 遺傳因子 x 流量為 1。隱性提供者特別判斷一下即可。

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
#include <bits/stdc++.h>
using namespace std;
const int MAXV = 40010;
const int MAXE = MAXV * 200 * 2;
const int INF = 1<<29;
typedef struct Edge {
    int v, cap, flow;
    Edge *next, *re;
} Edge;
class MaxFlow {
    public:
    Edge edge[MAXE], *adj[MAXV], *pre[MAXV], *arc[MAXV];
    int e, n, level[MAXV], lvCnt[MAXV], Q[MAXV];
    void Init(int x) {
        n = x, e = 0;
        for (int i = 0; i < n; ++i) adj[i] = NULL;
    }
    void Addedge(int x, int y, int flow){
        edge[e].v = y, edge[e].cap = flow, edge[e].next = adj[x];
        edge[e].re = &edge[e+1], adj[x] = &edge[e++];
        edge[e].v = x, edge[e].cap = 0, edge[e].next = adj[y];
        edge[e].re = &edge[e-1], adj[y] = &edge[e++];
    }
    void Bfs(int v){
        int front = 0, rear = 0, r = 0, dis = 0;
        for (int i = 0; i < n; ++i) level[i] = n, lvCnt[i] = 0;
        level[v] = 0, ++lvCnt[0];
        Q[rear++] = v;
        while (front != rear){
            if (front == r) ++dis, r = rear;
            v = Q[front++];
            for (Edge *i = adj[v]; i != NULL; i = i->next) {
                int t = i->v;
                if (level[t] == n) level[t] = dis, Q[rear++] = t, ++lvCnt[dis];
            }
        }
    }
    int Maxflow(int s, int t){
        int ret = 0, i, j;
        Bfs(t);
        for (i = 0; i < n; ++i) pre[i] = NULL, arc[i] = adj[i];
        for (i = 0; i < e; ++i) edge[i].flow = edge[i].cap;
        i = s;
        while (level[s] < n){
            while (arc[i] && (level[i] != level[arc[i]->v]+1 || !arc[i]->flow)) 
                arc[i] = arc[i]->next;
            if (arc[i]){
                j = arc[i]->v;
                pre[j] = arc[i];
                i = j;
                if (i == t){
                    int update = INF;
                    for (Edge *p = pre[t]; p != NULL; p = pre[p->re->v])
                        if (update > p->flow) update = p->flow;
                    ret += update;
                    for (Edge *p = pre[t]; p != NULL; p = pre[p->re->v])
                        p->flow -= update, p->re->flow += update;
                    i = s;
                }
            }
            else{
                int depth = n-1;
                for (Edge *p = adj[i]; p != NULL; p = p->next)
                    if (p->flow && depth > level[p->v]) depth = level[p->v];
                if (--lvCnt[level[i]] == 0) return ret;
                level[i] = depth+1;
                ++lvCnt[level[i]];
                arc[i] = adj[i];
                if (i != s) i = pre[i]->re->v;
            }
        }
        return ret;
    }
} g;
int main() {
    int N, Q, B, x;
    while (scanf("%d %d", &N, &Q) == 2) {
        set< set<int> > S;
        set<int> A[128];
        for (int i = 0; i < N; i++) {
            scanf("%d", &B);
            for (int j = 0; j < B; j++) {
                scanf("%d", &x);
                A[i].insert(x);
            }
            if (B != N)	A[i].insert(0);
        }
        
        for (int i = 0; i < Q; i++) {			
            int source = 2*N+2, sink = 2*N+3;
            g.Init(2*N+5);
            
            int used[128] = {};
            for (int j = 0; j < N; j++)	// parent
                g.Addedge(source, j, 1);
                
            scanf("%d", &B);
            for (int j = 0; j < B; j++) {
                scanf("%d", &x);
                g.Addedge(N+x, sink, 1);
                for (int k = 0; k < N; k++) {
                    if (A[k].count(x)) {
                        g.Addedge(k, N+x, 1);
                        used[k] = 1;
                    }
                }
            }
            
            int allused = 1;
            for (int j = 0; j < N; j++) {
                if (A[j].count(0) && B != N)
                    used[j] = 1;
                allused &= used[j];
            }
            
            if (!allused) {
                puts("N");
                continue;
            }
            
            int flow = g.Maxflow(source, sink);
            puts(flow == B ? "Y" : "N");
        }
    }
    return 0;
}
Read More +
解題區/解題區 - 其他題目

2016 Facebook Hacker Cup Round 2

這一場是凌晨的比賽，因此沒有興趣參與，還是來翻譯一下吧！

A. Boomerang Decoration

要製作回力鏢的雙翼，需要兩個對稱形狀的翼構成，簡單來說是兩個相同字串。給予左翼 $L$ 和右翼 $R$，現在 A 和另一名小夥伴 B 同時更改左翼 $L$，使得 $L = R$，請問最少時間為何。

每一時刻，$A$ 可以選擇 $[0, x]$ 都改變成顏色 $c_1$，而 $B$ 則可以選擇 $[y, n-1]$ 全部變成顏色 $c_2$，要求 $A$ 和 $B$ 的塗色區間不可重疊，意即 $y > x$。

由於每一次塗色是前綴或者後綴，那對於 A 而言，一定是選擇 $x$ 嚴格遞減，反之對於 $y$ 是嚴格遞增，否則之前的操作會白工。藉由上述推測得知，A 和 B 塗色的區間不會重疊，得到最後一定要求左區段和右區段工作時間最大值。

目標找到 A 完成左區段 dpL[0...x] 的最少時間，以及 B 完成 dpR[y...n-1] 的最少時間，最後答案為 min(max(dpL[0...x], dpR[x+1...n-1]))。計算 dpL[] 和 dpR[] 的方法是一樣的，左右相反即可，

由於只能更改左翼，右翼不能更改，那麼一旦改了左翼位置 $x$，事實上要全部變動成右翼，變數個數 $d$ 即是右翼不同的連續相同字符片段。根據貪心策略，若左翼位置 $x$ 和右翼位置 $x$ 相同，則變動費用為 dpL[x] = dpL[x-1]，否則就是 dpL[x] = d。時間複雜度 $\mathcal{O}(n)$。

B. Carnival Coins

參與一場遊戲，獎品無限多個，遊戲規則要求一次投擲 $n$ 硬幣，硬幣出現正面的機率 $P$，只要出現 $K$ 個以上正面即可獲得一份獎品，現在手裡握有 $N$ 個硬幣，你可以邀請很多個小夥伴幫忙參與遊戲，將這 $N$ 枚硬幣分給小夥伴玩。在最佳策略下，請問獲得獎品個數的期望值為何？

現在的目標要找到怎麼分，小夥伴個數不是問題，每一個小夥伴要拿多少枚硬幣是主要問題。定義 dp[i][j] 為投擲 $i$ 枚硬幣恰好 $j$ 枚正面的機率，遞推得到 dp[i][j] = dp[i-1][j-1]*P + dp[i-1][j]*(1-P)。由於大於等於 $K$ 枚只能獲得一份獎品，定義 dpw[i] 為投擲 $i$ 枚硬幣獲得一份獎品的期望值，遞推得到 dpw[i] = sum(dp[i][j]), j >= K。

最後，定義 dpw[i] 為分配 $i$ 枚硬幣給小夥伴的最大期望值個數，得到 dpw[i] = max(dpw[i-j]+dpv[j])。每一步都是 $\mathcal{O}(N^2)$，時間複雜度 $\mathcal{O}(N^2)$。

C. Snakes and Ladders

有一個奇怪的收藏家，收集很多梯子，每一個梯子位於水平位置 $x_i$ 高度為 $H_i$，然而當地特有蛇種喜歡懸掛在兩個等高 $H_p$ 的梯子之間，並且兩個梯子中間的其餘梯子都滿足 $H_r \le H_p$。若蛇懸掛在兩個距離 $L$ 的梯子之間，則需要餵養 $L^2$ 單位的飼料，請問飼主一天要餵多少單位的飼料。

搭配組合的數學計算，由於懸掛的組合保證中間不會有高於兩側，則可以用單調堆完成紀錄，由左而右依序加入梯子，在單調堆中維護梯子高度遞減的位置。由於等高的梯子個數需要合併，否則沒辦法計算 $L$，為了計算 $L^2$。當從左而右加入梯子時，往左找到合法的等高梯子，由於得知 $\sum (X - x_i)^2$，把式子展開得到 $nX^2 - 2X \sum x_i + \sum x_i^2$，因此需要對於每一個高度，要合併 $\sum x_i$ 和 $\sum x_i^2$。

由於每一個元素至多 push 和 pop 一次，時間複雜度 $\mathcal{O}(N)$。

D. Costly Labels

在一棵 $N$ 個節點的無根樹，每一個節點可填入 $K$ 種顏色，對於每一個節點填不同顏色都有不同花費，若一個節點 $u$ 的鄰居中有兩個以上節點顏色相同，則 $u$ 需要額外花費 $P$ 來維持平衡。請問著色的最少花費為何。

若不考慮額外花費 $P$ 這一點，這題就是再簡單不過的貪心。從鴿籠原理知道，若一個點的鄰居大於 $K$ 個，則他必然存在兩個相同顏色的鄰居，一定需要花費 $P$ 維持平衡。由於圖給定一棵樹，自然而然地就設想到樹形 dp，維護子樹的最少花費進行合併。

由於子樹的根著色不用考慮，但方案合併要考慮父節點要著色的方案，因此得到狀態 dp[u][c1][c2] 節點 $u$ 著色 $c_1$，$u$ 的父節點著色 $c_2$。分成兩種狀態轉移，是否要花費 $P$，若 $u$ 花費 $P$，直接合併子樹 $v$ 的最小值 dp[v][?][c1]。若不花費 $P$，勢必要讓子樹都填入不同顏色，由於每一個節點配對顏色花費都不同，最小化著色方案是標準的二分圖帶權匹配，可以用 KM 算法或者是最少費用流完成，只要排除預設 $u$ 的父節點顏色建圖即可，幸好不花費的數值都小於等於 $K$，KM 算法提供 $\mathcal{O}(K^3)$。整體的時間複雜度最慘 $\mathcal{O}(N \times K^2 \times K^3)$。

Solution

A. Boomerang Decoration

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
#include <bits/stdc++.h>
using namespace std;
const int MAXS = 131072;
char sa[MAXS], sb[MAXS];
int Ldp[MAXS], Rdp[MAXS];
int main() {
    int testcase, cases = 0;
    int n;
    scanf("%d", &testcase);
    while (testcase--) {
        scanf("%d", &n);
        scanf("%s %s", sa+1, sb+1);
        
        memset(Ldp, 0, sizeof(Ldp));
        memset(Rdp, 0, sizeof(Rdp));
        
        int c;
        c = 0;
        for (int i = 1; i <= n; i++) {
            if (sb[i] != sb[i-1])
                c++;
            if (sa[i] != sb[i])
                Ldp[i] = c;
            else
                Ldp[i] = Ldp[i-1];
        }
        c = 0;
        for (int i = n, j = 1; i >= 1; i--, j++) {
            if (sb[i] != sb[i+1])
                c++;
            if (sa[i] != sb[i])
                Rdp[j] = c;
            else
                Rdp[j] = Rdp[j-1];
        }
        int ret = n;
        for (int i = 1; i <= n; i++) {
            ret = min(ret, max(Ldp[i], Rdp[n-i]));
        }
        printf("Case #%d: %d\n", ++cases, ret);
    }
    return 0;
}

B. Carnival Coins

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
#include <bits/stdc++.h>
using namespace std;
const int MAXN = 4096;
double dp[MAXN][MAXN], dpw[MAXN], dpv[MAXN];
int main() {
    int testcase, cases = 0;
    scanf("%d", &testcase);
    while (testcase--) {
        int N, K;
        double P;
        scanf("%d %d %lf", &N, &K, &P);
        memset(dp, 0, sizeof(dp));
        dp[0][0] = 1;
        for (int i = 0; i <= N; i++) {
            for (int j = 0; j <= i; j++) {
                dp[i+1][j] += dp[i][j] * (1-P);
                dp[i+1][j+1] += dp[i][j] * P;
            }
            dpw[i] = 0;
            for (int j = K; j <= i; j++)
                dpw[i] += dp[i][j];
        }
        
        memset(dpv, 0, sizeof(dpv));
        for (int i = 1; i <= N; i++) {
            for (int j = 1; j <= i; j++) {
                dpv[i] = max(dpv[i], dpw[j] + dpv[i-j]);
            }
        }
        
        double ret = dpv[N];
        printf("Case #%d: %.9lf\n", ++cases, ret);
    }
    return 0;
}

C. Snakes and Ladders

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
#include <bits/stdc++.h>
using namespace std;
const long long MOD = 1000000007;
struct State {
    int H;
    long long sum, sqsum, n;
    State(int H, long long sum = 0, 
        long long sqsum = 0, long long n = 0):
        H(H), sum(sum), sqsum(sqsum), n(n) {}
};
int main() {
    int testcase, cases = 0;
    scanf("%d", &testcase);
    while (testcase--) {
        int n, x, h;
        scanf("%d", &n);
        vector<std::pair<int, int>> A;
        for (int i = 0; i < n; i++) {
            scanf("%d %d", &x, &h);
            A.push_back(make_pair(x, h));
        }
        sort(A.begin(), A.end());
        long long ret = 0;
        stack<State> stk;
        for (int i = 0; i < n; i++) {
            while (!stk.empty() && stk.top().H < A[i].second)
                stk.pop();
            if (!stk.empty() && stk.top().H == A[i].second) {
                State e = stk.top();
                long long X = A[i].first;
                long long N = e.n;
                long long S = N*X%MOD * X%MOD;
                S = (S - X*2%MOD*e.sum%MOD + e.sqsum)%MOD;
                e.sum = (e.sum + X)%MOD;
                e.sqsum = (e.sqsum + X*X%MOD)%MOD;
                e.n++;
                ret = (ret + S)%MOD;
                stk.pop(), stk.push(e);
            } else {
                long long X = A[i].first;
                stk.push(State(A[i].second, X, X*X%MOD, 1));
            }
        }
        ret = (ret + MOD)%MOD;
        printf("Case #%d: %lld\n", ++cases, ret);
    }
    return 0;
}

D. Costly Labels

版本 1

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
#include <bits/stdc++.h>
using namespace std;
const int MAXN = 128;
const int MAXM = 1048576;
struct Node {
    int x, y, cap;
    int cost;// x->y, v
    int next;
} edge[MAXM];
class MinCost {
public:
    int e, head[MAXN], pre[MAXN], record[MAXN], inq[MAXN];
    int dis[MAXN];
    int n;
    const int INF = 0x3f3f3f3f;
    void Addedge(int x, int y, int cap, int cost) {
        edge[e].x = x, edge[e].y = y, edge[e].cap = cap, edge[e].cost = cost;
        edge[e].next = head[x], head[x] = e++;
        edge[e].x = y, edge[e].y = x, edge[e].cap = 0, edge[e].cost = -cost;
        edge[e].next = head[y], head[y] = e++;
    }
    pair<int, int> mincost(int s, int t) {
        int mncost = 0;
        int flow, totflow = 0;
        int i, x, y;
        while(1) {
        	for (int i = 0; i < n; i++)
        		dis[i] = INF;
            int oo = dis[0];
            dis[s] = 0;
            deque<int> Q;
            Q.push_front(s);
            while(!Q.empty()) {
                x = Q.front(), Q.pop_front();
                inq[x] = 0;
                for(i = head[x]; i != -1; i = edge[i].next) {
                    y = edge[i].y;
                    if(edge[i].cap > 0 && dis[y] > dis[x] + edge[i].cost) {
                        dis[y] = dis[x] + edge[i].cost;
                        pre[y] = x, record[y] = i;
                        if(inq[y] == 0) {
                            inq[y] = 1;
                            if(Q.size() && dis[Q.front()] > dis[y])
                                Q.push_front(y);
                            else
                                Q.push_back(y);
                        }
                    }
                }
            }
            if(dis[t] == oo)
                break;
            flow = 0x3f3f3f3f;
            for(x = t; x != s; x = pre[x]) {
                int ri = record[x];
                flow = min(flow, edge[ri].cap);
            }
            for(x = t; x != s; x = pre[x]) {
                int ri = record[x];
                edge[ri].cap -= flow;
                edge[ri^1].cap += flow;
                edge[ri^1].cost = -edge[ri].cost;
            }
            totflow += flow;
            mncost += dis[t] * flow;
        }
        return make_pair(mncost, totflow);
    }
    void init(int n) {
    	this->n = n;
        e = 0;
        for (int i = 0; i <= n; i++)
            head[i] = -1;
    }
} g;
int N, K, P;
int costG[1024][32];
vector<int> treeG[1024];
long long dp[1024][32];	// dp[subtree: u][color of u: c] = mincost, don't care penalty P
long long dpv[1024];
void dfs(int u, int p) {
    for (int i = 0; i < K; i++)
        dp[u][i] = costG[u][i];
    for (auto &v : treeG[u]) {
        if (v == p)	continue;
        dfs(v, u);
        
        int cost = P;
        for (auto &adj_v : treeG[v]) {
            if (adj_v == u)	continue;
            cost += dpv[adj_v];
        }
        
        if (treeG[v].size() > K) {
            for (int i = 0; i < K; i++)
                dp[u][i] += cost;
        } else {
            // find disjoint-color by maximum weighted matching
            for (int i = 0; i < K; i++) {	// u-color
                int source = 0, sink = treeG[v].size()+K+1;
                g.init(treeG[v].size()+K+2);
                for (int j = 0; j < K; j++) {
                    if (j == i)	continue;
                    g.Addedge(treeG[v].size()+j+1, sink, 1, 0);
                }
                for (int j = 0; j < treeG[v].size(); j++) {
                    int adj_v = treeG[v][j];
                    if (adj_v == u)	continue;
                    for (int k = 0; k < K; k++) {
                        if (k == i)	continue;
                        g.Addedge(j+1, treeG[v].size()+k+1, 1, dp[adj_v][k]);
                    }
                    g.Addedge(source, j+1, 1, 0);
                }
                pair<int, int> mm = g.mincost(source, sink);
                dp[u][i] += min(cost, mm.first);
            }
        }
    }
    dpv[u] = INT_MAX;
    for (int i = 0; i < K; i++)
        dpv[u] = min(dpv[u], dp[u][i]);
}
long long final(int v) {
    int cost = P;
    for (auto &adj_v : treeG[v])
        cost += dpv[adj_v];
    
    if (treeG[v].size() > K)
        return cost;
    // find disjoint-color by maximum weighted matching
    for (int i = 0; i < K; i++) {	// u-color
        int source = 0, sink = treeG[v].size() + K + 1;
        g.init(treeG[v].size()+K+2);
        for (int j = 0; j < K; j++) {
            g.Addedge(j+treeG[v].size()+1, sink, 1, 0);
        }
        for (int j = 0; j < treeG[v].size(); j++) {
            int adj_v = treeG[v][j];
            for (int k = 0; k < K; k++)
                g.Addedge(j+1, treeG[v].size()+k+1, 1, dp[adj_v][k]);
            g.Addedge(source, j+1, 1, 0);
        }
        
        pair<int, int> mm = g.mincost(source, sink);
        cost = min(cost, mm.first);
    }
    return cost;
}
int main() {
    int testcase, cases = 0;	
    scanf("%d", &testcase);
    while (testcase--) {
        scanf("%d %d %d", &N, &K, &P);
        for (int i = 1; i <= N; i++) {
            for (int j = 0; j < K; j++) {
                scanf("%d", &costG[i][j]);
            }
        }
        
        for (int i = 1; i <= N; i++)
            treeG[i].clear();
        for (int i = 1; i < N; i++) {
            int x, y;
            scanf("%d %d", &x, &y);
            treeG[x].push_back(y);
            treeG[y].push_back(x);
        }
        
        dfs(1, -1);
        
        long long ret = INT_MAX;
        for (int i = 0; i < K; i++)
            ret = min(ret, dp[1][i]);
        ret += final(1);
        printf("Case #%d: %lld\n", ++cases, ret);
    }
    return 0;
}

版本 2

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
#include <bits/stdc++.h>
using namespace std;
const int MAXN = 128;
const int MAXM = 1048576;
struct Node {
    int x, y, cap;
    int cost;// x->y, v
    int next;
} edge[MAXM];
class MinCost {
public:
    int e, head[MAXN], pre[MAXN], record[MAXN], inq[MAXN];
    int dis[MAXN];
    int n;
    const int INF = 0x3f3f3f3f;
    void Addedge(int x, int y, int cap, int cost) {
        edge[e].x = x, edge[e].y = y, edge[e].cap = cap, edge[e].cost = cost;
        edge[e].next = head[x], head[x] = e++;
        edge[e].x = y, edge[e].y = x, edge[e].cap = 0, edge[e].cost = -cost;
        edge[e].next = head[y], head[y] = e++;
    }
    pair<int, int> mincost(int s, int t) {
        int mncost = 0;
        int flow, totflow = 0;
        int i, x, y;
        while(1) {
        	for (int i = 0; i < n; i++)
        		dis[i] = INF;
            int oo = dis[0];
            dis[s] = 0;
            deque<int> Q;
            Q.push_front(s);
            while(!Q.empty()) {
                x = Q.front(), Q.pop_front();
                inq[x] = 0;
                for(i = head[x]; i != -1; i = edge[i].next) {
                    y = edge[i].y;
                    if(edge[i].cap > 0 && dis[y] > dis[x] + edge[i].cost) {
                        dis[y] = dis[x] + edge[i].cost;
                        pre[y] = x, record[y] = i;
                        if(inq[y] == 0) {
                            inq[y] = 1;
                            if(Q.size() && dis[Q.front()] > dis[y])
                                Q.push_front(y);
                            else
                                Q.push_back(y);
                        }
                    }
                }
            }
            if(dis[t] == oo)
                break;
            flow = 0x3f3f3f3f;
            for(x = t; x != s; x = pre[x]) {
                int ri = record[x];
                flow = min(flow, edge[ri].cap);
            }
            for(x = t; x != s; x = pre[x]) {
                int ri = record[x];
                edge[ri].cap -= flow;
                edge[ri^1].cap += flow;
                edge[ri^1].cost = -edge[ri].cost;
            }
            totflow += flow;
            mncost += dis[t] * flow;
        }
        return make_pair(mncost, totflow);
    }
    void init(int n) {
    	this->n = n;
        e = 0;
        for (int i = 0; i <= n; i++)
            head[i] = -1;
    }
} g;
int N, K, P;
int costG[1024][32];
vector<int> treeG[1024];
long long dp[1024][32][32];	// dp[subtree: u][color of u: c][color of u'parent: c]
long long dpv[1024][32];
void dfs(int u, int p) {
    for (auto &v : treeG[u]) {
        if (v == p)	continue;
        dfs(v, u);
    }
    
    for (int i = 0; i < K; i++) {
        for (int j = 0; j < K; j++) {
            int cost = costG[u][i] + P;
            for (auto &v : treeG[u]) {
                if (v == p)	continue;
                cost += dpv[v][i];
            }
            dp[u][i][j] = cost;
            
            if (treeG[u].size() > K)
                continue;
                
            int source = 0, sink = treeG[u].size()+K+1;
            g.init(treeG[u].size()+K+2);
            for (int k = 0; k < K; k++) {
                if (k == j && p != -1)	continue;	// parent color
                g.Addedge(treeG[u].size()+k+1, sink, 1, 0);
            }
            int branch = 0;
            for (int it = 0; it < treeG[u].size(); it++) {
                int v = treeG[u][it];
                if (v == p)	continue;
                branch++;
                for (int k = 0; k < K; k++) {
                    if (k == j && p != -1)	continue;
                    g.Addedge(it+1, treeG[u].size()+k+1, 1, dp[v][k][i]);
                }
                g.Addedge(source, it+1, 1, 0);
            }
            
            pair<int, int> mm = g.mincost(source, sink);
            if (mm.second != branch)
                continue;
            dp[u][i][j] = min(dp[u][i][j], (long long) mm.first + costG[u][i]);
        }
    }
    
    for (int i = 0; i < K; i++) {
        dpv[u][i] = INT_MAX;
        for (int j = 0; j < K; j++)
            dpv[u][i] = min(dpv[u][i], dp[u][j][i]);
    }
}
int main() {
    int testcase, cases = 0;	
    scanf("%d", &testcase);
    while (testcase--) {
        scanf("%d %d %d", &N, &K, &P);
        for (int i = 1; i <= N; i++) {
            for (int j = 0; j < K; j++) {
                scanf("%d", &costG[i][j]);
            }
        }
        
        for (int i = 1; i <= N; i++)
            treeG[i].clear();
        for (int i = 1; i < N; i++) {
            int x, y;
            scanf("%d %d", &x, &y);
            treeG[x].push_back(y);
            treeG[y].push_back(x);
        }
        
        dfs(1, -1);
        
        long long ret = INT_MAX;
        for (int i = 0; i < K; i++)
            for (int j = 0; j < K; j++)
                ret = min(ret, dp[1][i][j]);
        printf("Case #%d: %lld\n", ++cases, ret);
    }
    return 0;
}
Read More +
解題區/解題區 - Zerojudge

b685. 高中組第五題-課堂抽籤

Problem

$N$ 個人抽籤決定跟誰同一組，請問在某些籤已經決定好的情況下，求恰好分成 $M$ 組的抽籤方法數。

轉換成 $N$ 個節點，恰好形成 $M$ 個環的方法數。

## Sample Input ##

1
2
3
4
5
2
5 1
2 3 1 5 4
5 2
2 3 1 5 4

Sample Output

1
2
0
1

Solution

在此特別感謝 lwc (Wei Chen Liao) 提供思路。

手動暴力打表後，上網蒐到 A125714 Alfred Moessner’s factorial triangle. 恰好是這個數列的答案。

由於每一個人的籤不重複，若把他抽到的籤當作指向出去的邊，那麼保證這一個點恰好一個出邊和一個入邊。當所有籤抽滿後，圖看起來是好幾個環。

遞迴定義 dp[i][j] 表示有 $i$ 個節點和 $j$ 個環。考慮新加入得點要抽到哪一個籤，加入到 $j$ 個環的情況，則是把某一個環的某一個位置打開加入，或者自身形成一個環，最後得到 dp[i][j] = dp[i-1][j-1] + (i-1)*dp[i-1][j]

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
#include <bits/stdc++.h>
using namespace std;
const int MAXN = 1024;
const long long MOD = 1000007;
long long dp[MAXN][MAXN];
int A[MAXN], used[MAXN];
int main() {
    dp[0][0] = dp[1][1] = 1;
    for (int i = 2; i < MAXN; i++) {
        for (int j = 1; j <= i; j++)
            dp[i][j] = (dp[i-1][j-1] + dp[i-1][j] * (i-1))%MOD;
    }
    int testcase, n, m;
    scanf("%d", &testcase);
    while (testcase--) {
        scanf("%d %d", &n, &m);
        for (int i = 1; i <= n; i++)
            scanf("%d", &A[i]);
        int a = n, b = 0;
        memset(used, 0, sizeof(used));
        for (int i = 1; i <= n; i++) {
            if (used[i])
                continue;
            int x = i, cnt = 0;
            while (x && used[x] == 0)
                used[x] = 1, x = A[x], cnt++;
            a -= cnt-1;
            if (x == i)	m--;
            if (x == 0)	b++;
        }
        if (m < 0)
            puts("0");
        else
            printf("%lld\n", dp[b][m]);
    }
    return 0;
}
Read More +
解題區/解題區 - UVa

UVa 12415 - Digit Patterns

Problem

讀入一個 regex，一個主字串 $s$，請問有不同的 $i$ 滿足 $s$ 的子字串 s[j...i] 匹配 regex。

正規表達式長度最多 500，主字串 $s$ 長度最多 $10^7$。

Sample Input

1
2
3
4
6 1(2+3)*4
012345
2 00*(10+100)*
00100

Sample Output

1
2
5
1 2 4 5

Solution

隔了一年才解決它，這一題很明顯地在編譯器的範疇，要將一個 regex 轉換成 NFA 甚至 DFA。一開始設想轉換成 DFA，但長度 500 轉換成 DFA，用狀態壓縮的方式狀態增長非常大，不知道是不是指數成長，上傳就得到 Runtime Error。

將 regex 轉換成 NFA (nondeterministic finite automata) 並不難，但節點會很多且存在很多 $a \overset{\varepsilon}{\rightarrow} b$ 這種類型的邊。在解析 regular expression 時，用線性 $\mathcal{O}(n)$ 或者是 $\mathcal{O}(n^2)$ 都沒關係。因為在計算答案至少 $\mathcal{O}(10^7)$。

得到最原生的 NFA 後，接下來要思考怎麼找到匹配的子字串 $s_{j, i}$，維護可轉移的狀態集合 $S$，若當前 $S$ 中存在 acceptable 狀態，則找到一個 $i$ 解。依序加入字符 $s_i$，維護 $S$ 的做法如下：

  • 將起始狀態丟入 $S$，計算閉包 (closure)，也就是所有 $a \overset{\varepsilon}{\rightarrow} b$ 能連到的狀態都應該被丟入 $S$。
  • 定義下一個狀態集合 $S'$，在 $S$ 中的每一個狀態 $q$，都嘗試藉由 $s_i$ 轉移到 $q'$，將所有 $q'$ 丟入 $S'$ 中。
  • 若 $\text{closure}(S')$ 存在 acceptable 狀態，則 $i$ 是一組解。
  • $S = \text{closure}(S' + \mathit{q}_{start})$。

上面的算法並沒有錯，但原生的 NFA 產生方式在計算閉包時很慢，過多的邊和重複狀態導致。因此需要重建邊，預處理每一個字符轉移，移除所有 $a \overset{\varepsilon}{\rightarrow} b$，這麼一來就能在時限內通過。時間複雜度 $\mathcal{O}(|s| \times (\text{state} + \text{edge}))$

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
#include <bits/stdc++.h>
using namespace std;
const int MAXSTATE = 32767;
const int epsilon = 0;
int toIndex(char c) {
    return c - '0' + 1;
}
struct State {
    vector<State*> trans[11];
    int ac, label;
    State() {
        ac = label = 0;
        for (int i = 0; i < 11; i++)
            trans[i].clear();
    }
} _mem[MAXSTATE];
int nodesize = 0;
State *newState() {
    assert(nodesize < MAXSTATE);
    State *p = &_mem[nodesize++];
    *p = State();
    return p;
}
struct NFA {
    State *st, *ac;
    NFA() {
        st = ac = NULL;
    }
    NFA(char c) {
        st = newState(), ac = newState();
        ac->ac = 1;
        st->trans[toIndex(c)].push_back(ac);
    }
    NFA(NFA L, char c) {	// (A)*
        st = ac = NULL;
        if (c != '*')
            return ;
        st = newState(), ac = newState();
        ac->ac = 1;
        st->trans[epsilon].push_back(L.st);
        L.st->trans[epsilon].push_back(L.ac);
        L.ac->trans[epsilon].push_back(L.st);
        L.ac->trans[epsilon].push_back(ac);
        L.ac->ac = 0;
    }
    NFA(NFA L, NFA R, char c) {
        if (R.st == NULL) {
            st = L.st, ac = L.ac;
            return;
        }
        if (c == '|') {
            st = newState(), ac = newState();
            ac->ac = 1;
            st->trans[epsilon].push_back(L.st);
            st->trans[epsilon].push_back(R.st);
            L.ac->trans[epsilon].push_back(ac);
            R.ac->trans[epsilon].push_back(ac);
            L.ac->ac = R.ac->ac = 0;
        } else if (c == '&') {
            st = newState(), ac = newState();
            ac->ac = 1;
            st->trans[epsilon].push_back(L.st);
            L.ac->trans[epsilon].push_back(R.st);
            R.ac->trans[epsilon].push_back(ac);
            L.ac->ac = R.ac->ac = 0;
        }
    }
};
NFA parser(string exp) {
    if (exp.length() == 0)
        return NFA();
    if (exp.length() == 1)
        return NFA(exp[0]);
        
    int l = 0, pos = -1;
    for (int i = 0; i < exp.length(); i++) {
        if (exp[i] == '(') {
            l++;
        } else if (exp[i] == ')') {
            l--;
        } else if (exp[i] == '+') {
            if (l == 0) {
                pos = i;
                break;
            }
        }
    }
    
    if (pos != -1) {
        NFA L = parser(exp.substr(0, pos));
        NFA R = parser(exp.substr(pos+1));
        return NFA(L, R, '|');
    }
    
    int hasStar = 0;
    string ls, rs;
    
    if (exp[0] == '(') {
        for (int i = 0; i < exp.length(); i++) {
            if (exp[i] == '(') {
                l++;
            } else if (exp[i] == ')') {
                l--;
                if (l == 0) {	// (...)...
                    if (i+1 < exp.length() && exp[i+1] == '*') {	// (...)*...
                        hasStar = 1;
                        ls = exp.substr(1, i-1), rs = exp.substr(i+2);
                    } else {	// (...)...
                        ls = exp.substr(1, i-1), rs = exp.substr(i+1);
                    }
                    break;
                }
            }
        }
    } else {	// ...(...) or ...*... or ......
        for (int i = 0; i < exp.length(); i++) {
            if (exp[i] == '(') {
                l++;
            } else if (exp[i] == ')') {
                l--;
            }
            if (l == 0) {
                if (i+1 < exp.length() && exp[i+1] == '*') {
                    hasStar = 1;
                    ls = exp.substr(0, i+1), rs = exp.substr(i+2);
                } else {
                    ls = exp.substr(0, i+1), rs = exp.substr(i+1);
                }
                break;
            }
        }
    }
    for (int i = 0; rs.length() > 0; ) {
        while (i < rs.length() && rs[i] == '*')
            i++;
        rs = rs.substr(i);
        break;
    }
    NFA L = parser(ls);
    NFA R = parser(rs);
    if (hasStar)
        L = NFA(L, '*');
    return NFA(L, R, '&');
}
State *gmap[MAXSTATE];
int relabel(NFA A) {
    int size = 0;
    State *u, *v;
    queue<State*> Q;
    Q.push(A.st), A.st->label = ++size;
    while (!Q.empty()) {
        u = Q.front(), Q.pop();
        gmap[u->label] = u;
        for (int it = 0; it < 11; it++) {
        	for (int i = 0; i < u->trans[it].size(); i++) {
        		v = u->trans[it][i];
                if (v->label == 0) {
                    v->label = ++size;
                    Q.push(v);
                }
            }
        }
    }
    return size;
}
char s[16777216];
int used[MAXSTATE], ACable[MAXSTATE];
int closure(vector<int> &A, int x, int cases) {
    queue<int> Q;
    State *u;
    int accept = 0;
    if (used[x] != cases)
        A.push_back(x), used[x] = cases;
    Q.push(x);
    while (!Q.empty()) {
        x = Q.front(), Q.pop();
        u = gmap[x];
        accept |= u->ac;
        if (u->trans[epsilon].size() == 0)
            continue;
        for (auto &y : u->trans[epsilon]) {
            if (used[y->label] != cases) {
                Q.push(y->label), used[y->label] = cases;
                A.push_back(y->label);
            }
        }
    }
    return accept;
}
void rebuild(int n, int ctype) {
    memset(ACable, 0, sizeof(ACable));
    memset(used, 0, sizeof(used));
    int cases = 0;
    for (int i = 1; i <= n; i++) {
        cases++;
        vector<int> cc;
        closure(cc, i, cases);
        for (int j = 1; j <= ctype; j++) {
            set<int> S;
            for (int k = 0; k < cc.size(); k++) {
                State *u = gmap[cc[k]];
                for (auto &p : u->trans[j])
                    S.insert(p->label);
                ACable[i] |= u->ac;
            }
            State *u = gmap[i];
            u->trans[j].clear();
            for (auto &x : S)
                u->trans[j].push_back(gmap[x]);
        }
    }
}
int main() {
    int n, m;
    char regex[512];
    while (scanf("%d %s", &n, regex) == 2) {
        nodesize = 0;	// global
        NFA nfa = parser(regex);
        m = relabel(nfa);
        rebuild(m, n);
        scanf("%s", s);
        
        memset(used, 0, sizeof(used));
        int cases = 0, flag = 0;
        
        vector<int> A;
        cases++;
        A.push_back(1);
        for (int i = 0; s[i]; i++) {
            vector<int> next;
            cases++;
            int accept = 0;
            for (int j = 0; j < A.size(); j++) {
                int x = A[j];
                State *u = gmap[x];
                if (u->trans[toIndex(s[i])].size() == 0)
                    continue;
                for (auto &y : u->trans[toIndex(s[i])]) {
                    if (used[y->label] != cases) {
                        used[y->label] = cases, next.push_back(y->label);
                        accept |= ACable[y->label];
                    }
                }
            }
            if (used[1] != cases)
                used[1] = cases, next.push_back(1);
            A = next;
            if (accept) {
                if (flag)	putchar(' ');
                flag = 1;
                printf("%d", i+1);
            }
        }
        puts("");
    }
    return 0;
}
/*
6 1(2+3)*4
012345
2 00*(10+100)*
00100
*/
Read More +
解題區/出題解題

[讀檔操作] 有限記憶體排序數組

背景

現在的系統架構中，內存 (Memory，台灣翻譯：記憶體，大陸翻譯：內存) 的大小仍無法完全像硬碟機一樣。如檔案大小 64GB，內存只有 4GB，處理檔案無法全部 In Memory，當然最近的硬體技術也逐漸朝著 All In Memory 的方式進行加速。回過頭來，在內存不足的情況下，來練習如何處理檔案大小遠大於內存的情況吧。

題目描述

給予一個 binary file 的檔案名稱，裡面用 4 個位元組為一個 signed 32-bit 整數，有數個直到 EOF，請排序後以 binary mode 導入 stdout 輸出。

輸入格式

標準輸入串流只有一行一個字串，表示檔案名稱 $F$。檔案大小最多 16 MB，而你的程序被限制最多使用 4 MB 的內存。

每個整數 $x$，限制條件 $0 \le x \le 10^9$

輸出格式

輸出在標準串流以 binary mode 下，請避開使用 console 輸出，會因為特殊字元 (如 「嗶」一聲) 導致系統嚴重當機。

範例輸入

1
1.dat

範例輸出

1
... binary mode 無法顯示

提示

download p10068-in.dat

1
2
3
4
5
6
7
8
9
morris821028@morris821028-PC MINGW64 ~/Desktop/10068
$ ./sort2.exe >test.out
1.dat
morris821028@morris821028-PC MINGW64 ~/Desktop/10068
$ xxd test.out
00000000: 0000 0000 0100 0000 0100 0000 0100 0000  ................
00000010: 0300 0000 0400 0000 0500 0000 0500 0000  ................
00000020: 0800 0000 0900 0000                      ........

1.dat 以 binary file 儲存 10 個 4 bytes 整數，依序為 5, 9, 3, 5, 0, 1, 1, 8, 4, 1，排序後即為 0, 1, 1, 1, 3, 4, 5, 5, 8, 9。

由於限制內存使用量，無法寫入暫存檔案，可利用多次讀檔完成。

Solution

乍看之下，這一題是最常見的 external sort (外部排序)，外部排序需要額外寫檔案，由於記憶體用量計算，一寫檔案會計算到使用的記憶體中，實際上這一題不寫檔也是能解決的。

  • 給定要排序的數據範圍
  • 二分搜尋可容納的排序範圍，利用平衡樹或者 hash 來完成計算 <某整數, 某整數出現個數>
  • 將可容納到 hash 的所有數字排序，再將其直接寫到輸出檔案。
  • 不斷地遞增二分搜尋的左邊邊界。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
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
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <limits.h>
#define swap(x, y) {int t; t = x, x = y, y = t;}
int fsize(FILE *fp) {
    int prev = ftell(fp);
    fseek(fp, 0L, SEEK_END);
    int size = ftell(fp);
    fseek(fp, prev, SEEK_SET);
    return size;
}
#define USAGEMEM (2<<20)
#define MAXELE ((2<<20)/8)
int A[MAXELE];
#define HSIZE 100007
#define HNODE 50000
typedef struct HashNode {
    int f, cnt;
    struct HashNode *next;
} HashNode;
HashNode *hhead[HSIZE], hnodes[HNODE];
int nodesize = 0;
int cmp(const void *x, const void *y) {
    return ((HashNode*) x)->f - ((HashNode*) y)->f;
}
unsigned int hash(int f) {
    unsigned int a = 63689, b = 378551;
    unsigned int value = 0;
    value = value * a + f, a = a * b;
    return value;
} 
void initHash() {
    memset(hhead, 0, sizeof(hhead));
    nodesize = 0;
}
int insertHash(int f) {
    unsigned int hidx = hash(f)%HSIZE;
    for (HashNode *p = hhead[hidx]; p != NULL; p = p->next) {
        if (f == p->f) {
        	p->cnt++;
            return 1;
        }
    }
    if (nodesize >= HNODE)	return 0;
    HashNode s = {.f = f, .cnt = 1, .next = NULL};
    hnodes[nodesize] = s;
    hnodes[nodesize].next = hhead[hidx];
    hhead[hidx] = &hnodes[nodesize];
    nodesize++;
    return 1;
}
int merge_int(FILE *fp, int l, int r) {
    initHash();
    int ret = 0, x, n = 0;
    fseek(fp, 0, SEEK_SET);
    while ((n = fread(A, sizeof(int), MAXELE, fp))) {
        for (int i = 0; i < n; i++) {
            if (A[i] >= l && A[i] <= r) {
                if (!insertHash(A[i]))
                    return 0;
            }
        }
    }
    return 1;
}
#ifdef _WIN32
#include <fcntl.h>
#endif
int block_process(FILE *fp) {
    int base = 0;
    int l, r, mid, ret;
    l = base, r = 1000000000, ret = 0;
#ifdef _WIN32
    if (setmode(fileno(stdout), O_BINARY) == -1) {
        perror("Cannot set stdout to binary mode");
        return 2;
    }
#endif
#ifdef __linux__
    FILE *const out = fdopen(dup(fileno(stdout)), "wb");
#endif
    while (l <= r) {
        mid = (l + r)/2;
        int status = merge_int(fp, base, mid);
        if (status == 1) {
            l = mid + 1, r = 1000000000;
            base = mid + 1;
            qsort(hnodes, nodesize, sizeof(HashNode), cmp);
            for (int i = 0; i < nodesize; i++) {
                int x = hnodes[i].f;
                for (int j = hnodes[i].cnt-1; j >= 0; j--) {
#ifdef _WIN32
                    fwrite(&x, sizeof(int), 1, stdout);
#endif
#ifdef __linux__
                    fwrite(&x, sizeof(int), 1, out);
#endif
                }
                ret += hnodes[i].cnt;
            }
        } else {
            r = mid - 1;
        }
    }
    return ret;
}
int main() {
    char fName[128];
    scanf("%s", fName);
    FILE *fin = fopen(fName, "rb");
    
    int n = fsize(fin) / sizeof(int);
    block_process(fin);
    return 0;
}
Read More +
解題區/出題解題

[讀檔操作] 有限記憶體找中位數

背景

現在的系統架構中，內存 (Memory，台灣翻譯：記憶體，大陸翻譯：內存) 的大小仍無法完全像硬碟機一樣。如檔案大小 64GB，內存只有 4GB，處理檔案無法全部 In Memory，當然最近的硬體技術也逐漸朝著 All In Memory 的方式進行加速。回過頭來，在內存不足的情況下，來練習如何處理檔案大小遠大於內存的情況吧。

題目描述

給予一個 binary file 的檔案名稱，裡面用 4 個位元組為一個 signed 32-bits 整數，有數個直到 EOF，保證恰好有奇數個，請輸出中位數為何？

輸入格式

標準輸入串流只有一行一個字串，表示檔案名稱 $F$。檔案大小最多 16 MB，而你的程序被限制最多使用 4 MB 的內存。

每個整數 $x$，限制條件 $0 \le x \le 10^9$

輸出格式

輸出一行整數 $Median$ 為 binary file 的中位數。

範例輸入

1
p10067-in.dat

範例輸出

1
97537111

提示

二分搜尋，分桶

download p10067-in.dat

p10067-in.dat

1
2
3
10825098
97537111
208681850

Solution

這一個問題很久以前見過，對岸通常叫做「海量資料找中位數」約束在記憶體不夠的情況下，如何高效率找到中位數，操作時允許重複讀檔。

遲遲沒有 Online Judge 進行記憶體用量檢測以及開放讀寫檔案，如今有機會擔任台灣大學計算機程式設計的課程助教，架了一個允許亂來的 Online Judge - Judge Girl (中文翻譯：批改娘系統)，讀寫檔案也不會被封鎖！洽詢 網站，目前只針對課堂學生開放。

回到問題，題目給訂 16MB 的整數序列，最多 $n = 4 \times 10^6$，由於記憶體不足沒辦法直接宣告陣列大小為 $n$ 的整數陣列，又由於數字可能會重複，就沒辦法利用 bitmask 進行壓縮，只能倚靠不斷地讀檔案進行計算。若以數值範圍 $[0, V]$，大致放分成兩種策略

  • 二分搜尋 (AC, 250ms)，效率 $\mathcal{O}(V \log V)$，開檔次數 $\mathcal{O}(\log V)$。
  • 分桶算法 (AC, 70ms)，效率 $\mathcal{O}(N)$，開檔次數 $\mathcal{O}(1)$。

讀取檔案時，在有限記憶體空間，配置一塊作為緩衝區，一次讀入一坨子序列，一個一個讀取效率非常差，別輕忽從磁碟讀取資料的速度。剩餘空間再進行紀錄用途。從效能比較 二分搜尋 慢於 分桶算法。

二分搜尋

二分答案 $x$，接著線性掃描序列，計算有多少個整數小於等於 $x$，藉此找到中位數。

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
#include <stdio.h>
int fsize(FILE *fp) {
    int prev = ftell(fp);
    fseek(fp, 0L, SEEK_END);
    int size = ftell(fp);
    fseek(fp, prev, SEEK_SET);
    return size;
}
#define MAXELE (3<<20)/4
int countLess(FILE *fp, int ans) {
    static int A[MAXELE];
    int ret = 0, x, n = 0;
    fseek(fp, 0, SEEK_SET);
    while ((n = fread(A, sizeof(int), MAXELE, fp))) {
        for (int i = 0; i < n; i++)
            ret += A[i] <= ans;
    }
    return ret;
}
int main() {
    char fName[128];
    scanf("%s", fName);
    FILE *fin = fopen(fName, "rb");
    
    int n = fsize(fin) / sizeof(int);
    
    int m = n/2 + 1;
    int l = 0, r = 1000000000, mid, ret = 0;
    
    while (l <= r) {
        mid = (l + r)/2;
        int tn = countLess(fin, mid);
        if (tn < m)
            l = mid + 1, ret = mid+1;
        else
            r = mid - 1;
    }
    printf("%d\n", ret);
    return 0;
}

分桶算法

分桶算法分成兩次掃瞄，假設內存分配最多 $I$ 個 32bits 整數 ($I \ll N$)。

  • 由於記憶體限制大小關係，將數值範圍分成 $I$ 堆，如 [0, V/I-1]、[V/I, 2V/I-1]、… 等，計算區間範圍內的個數，讀取一次檔案，可以得到中位數介於 [iV/I, (i+1)V/I-1]。
  • 接著再劃分一次，直到區間長度為 1，中位數便可得到。
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
#include <stdio.h>
#include <assert.h>
#include <string.h>
int fsize(FILE *fp) {
    int prev = ftell(fp);
    fseek(fp, 0L, SEEK_END);
    int size = ftell(fp);
    fseek(fp, prev, SEEK_SET);
    return size;
}
#define MAXELE ((1<<20)/4)
static int A[MAXELE];
static int B[MAXELE];
int countLess(FILE *fp, int ans) {
    int ret = 0, x, n = 0;
    fseek(fp, 0, SEEK_SET);
    while ((n = fread(A, sizeof(int), MAXELE, fp))) {
        for (int i = 0; i < n; i++)
            ret += A[i] <= ans;
    }
    return ret;
}
int main() {
    char fName[128];
    scanf("%s", fName);
    FILE *fin = fopen(fName, "rb");
    
    int n = fsize(fin) / sizeof(int);
    
    int m = n/2 + 1;
    int l = 0, r = 1000000000, mid, ret = 0;
    
    int WIDTH = r, SHIFT;
    
    for (SHIFT = 1; (1<<SHIFT) < MAXELE; SHIFT++);
    WIDTH = r / (1<<SHIFT);
    
    memset(B, 0, sizeof(B));
    fseek(fin, 0, SEEK_SET);
    while ((n = fread(A, sizeof(int), MAXELE, fin))) {
        for (int i = 0; i < n; i++) {
            assert((A[i]>>SHIFT) < MAXELE);
            B[A[i]>>SHIFT]++;
        }
    }
    
    int BASE = 0, SELECT;
    for (int i = 0, sum = m; i < MAXELE; i++) {
        sum -= B[i];
        if (sum <= 0) {
            BASE = i * (1<<SHIFT);
            SELECT = sum + B[i];
            break;
        }
    }
    
    memset(B, 0, sizeof(B));
    fseek(fin, 0, SEEK_SET);
    while ((n = fread(A, sizeof(int), MAXELE, fin))) {
        for (int i = 0; i < n; i++) {
            if (A[i] - BASE < MAXELE && A[i] >= BASE)
                B[A[i] - BASE]++;
        }
    }
    for (int i = 0, sum = SELECT; i < MAXELE; i++) {
        sum -= B[i];
        if (sum <= 0) {
            printf("%d\n", i + BASE);
            break;
        }
    }
    return 0;
}
Read More +
學校課程/數位影像生成

pbrt-v2 加速結構 BVH-contract 改寫

Description of implementation approach and comments

從一般 BVH 架構中，一般都是用 full binary tree，子節點要不 0 個要不 2 個。若有 $N$ 個 primitive 物件，則表示有 $N$ 個葉節點放置這些 primitive 物件，而有 $N-1$ 個內部節點紀錄 Bounding Box 的部分。在測試交點和遍歷走訪的使用上，最慘有一半是多餘計算和走訪，而另一半屬於加速結構。

在論文 Ray Specialized Contraction on Bounding Volume Hierarchies 中，企圖想要利用 generic tree 取代一般使用 full binary tree 實作，在不更動 BVH 的效能下，減少運行上較沒用的內部節點，用以提升遍歷走訪效能，以及提升內存快取的效率。

降低內部節點步驟如下所示：

  1. 利用原生方法建立 full binary tree 的 BVH (利用各種分割策略完成)
  2. 進行坍倒 (flatten)，將二元樹不連續的記憶體分布調整成線性分布，加速遍歷走訪的內存快取效率。
  3. 靜態調整成 generic tree 版本，藉由啟發式收縮 (Contract)，利用節點與其子節點的 Bounding Box 表面積比例，評估浪費的走訪節點。
  4. 動態調整部分，採用隨機取樣，根據取樣 ray，取樣走訪次數，將比較容易打到的節點盡可能收縮到接近根節點。

若要收縮節點 $N$，假設 Bounding box 計算交點花費為 $C_B$，穿過節點 $N$ 的 Bounding box 射線機率 $\alpha_N$，得到收縮前後的計算差值 $\delta(N)$，如下所示。

$$\begin{align} \delta(N) &= n_{N.\text{child}} \; C_B - (\alpha_N (1+n_{N.\text{child}}) +(1 - \alpha_N)) \; C_B \\ & = ((1 - \alpha_N) \; n_{N.\text{child}} - 1) \; C_B \end{align}$$

目標讓 $\delta(N) < 0$，得到

$\alpha(N) > 1 - \frac{1}{n_{N.\text{child}}}$

計算 $\delta(N)$ 顯得要有效率，但又沒辦法全局考量，需要提供一個猜測算法，藉由部分取樣以及步驟 2. 的表面積總和比例進行收縮。

實作部分

從實作中，在步驟 2. 約略可以減少 $25\%$ 的節點，在記憶體方面的影響量沒有太大影響，節點紀錄資訊也增加 (sizeof(struct Node) 相對地大上許多)。

在步驟 3. 中，根據 pbrt-v2 的架構，加速結構能取得的場景資訊並不容易竄改，大部分的類別函數都是 const function()，意即無法變動 object member 的值，但針對指向位址的值可以改變。這類寫法，猜想所有加速結構都是靜態決定，在多核心運行時較不會存在同步 overhead 的影響。

在此，換成另外一種修改方案，在 pbrt-v2/core/scene.h 的 bool Scene::Intersect(...) 函數中加入 aggregate->tick();。利用 aggregate->tick(); 這個函數，大部分呼叫都沒有更動樹狀結構。當呼叫次數 $T$ 達到一定次數時，加速結構會進行大規模的結構調整。

根據 pbrt rendering 的步驟，儘管不斷地測試或者是估計 $T$ 要設定的值，都無法接近全局的取樣評估，其中最大的原因是取樣順序和局部取樣調整，從理論上得知不可能會有比較好的結果。這樣的寫法提供簡便的方案去統計 pbrt 運算時較有可能的 ray 從哪裡射出，不用挖掘所有的不同類型 ray 進行取樣。

最後，修改檔案如下：

修改檔案清單

1
2
3
4
5
6
7
8
9
.
├── accelerators
│   ├── bvhcontract.cpp
│   └── bvhcontract.h
└── core
    ├── api.cpp
    ├── primitive.cpp
    ├── primitive.h
    └── scene.h

core/api.cpp

1
2
3
4
5
6
7
8
Primitive *MakeAccelerator(const string &name,
                           const vector<Reference<Primitive> > &prims,
                           const ParamSet &paramSet) {
    ...
    else if (name == "bvhcontract")
        accel = CreateBVHContractAccelerator(prims, paramSet);
    ...
}

core/primitive.h

1
2
3
4
5
6
7
8
9
class Primitive : public ReferenceCounted {
public:
    ...
    // MORRIS ADD
    virtual void tick();
    ...
protected:
    ...
};

core/scene.h

1
2
3
4
5
6
7
8
9
class Scene {
public:
    // Scene Public Methods
    bool Intersect(const Ray &ray, Intersection *isect) const {
        ...
        aggregate->tick();
        ...
    }
};

Generic Tree 設計與走訪測試

在 Full binary tree style - BVH 實作中，利用前序走訪分配節點編號範圍 $[0, 2N-1]$，因此節點編號 $u$ 的左子樹的編號為 $u+1$，只需要紀錄右子樹編號 secondChildOffset，這種寫法在空間和走訪時的快取都有效能改善。在標準程序中也單用迭代方式即可完成，不採用遞迴走訪，減少 push stack 的指令。

在 Generic tree 版本中，基礎節點紀錄架構如下：

1
2
3
4
5
6
7
8
9
10
11
struct LinearTreeNode {
    // node link information
    uint32_t parentOffset;
    uint32_t childOffsetHead, childOffsetTail;
    uint32_t siblingOffsetNext, siblingOffsetPrev;
    // faster record
    uint32_t numChild;
    // node data
    TreeData e;
    uint32_t visitCount;
};

在原生 BVH 求交點寫法中，根據節點的 Split-Axis 以及 Ray 的方向決定先左子樹還是先右子樹走訪，藉以獲得較高的剪枝機率。但全部先左子樹走訪的快取優勢比較高 (因為前序分配節點編號)，反之在 Split-Axis 有一半的機率走會走快取優勢不高的情況，在權衡兩者之間。

然而在 Generic Tree 實作上，若要提供 Split-Axis 則需要提供 childOffsetTail 和 siblingOffsetPrev 兩個指針，則多了兩個紀錄欄位，單一節點大小從 sizeof(LinearBVHNode) = 32拓展到 sizeof(LinearBVHContractNode) = 60，記憶體用量整整接近兩倍。從 Contract 算法中，節點個數保證無法減少一半，推論得到在 Contract 後記憶體用量會多一些。

走訪實作上分成遞迴和迭代兩種版本，遞迴在效能上會卡在 push esp, argument 等資訊上，而在迭代版本省了 call function overhead 和空間使用，但增加計算次數。而在我撰寫的版本中，還必須 access 父節點的資訊決定要往 siblingOffsetNext 還是 siblingOffsetprev，因此快取效能從理論上嚴重下滑。

遞迴和迭代走訪寫法如下：

遞迴版本走訪

1
2
3
4
5
6
7
8
9
10
void recursiveTraverse(LinearTreeNode *node, LinearTreeNode *_mem) {
    // proess
    uint32_t offset = node->childOffsetHead;
    if (offset == -1)
        return ;
    for (LinearTreeNode *u; offset != -1; offset = u->siblingOffsetNext) {
        u = &_mem[offset];
        recursiveTraverse(u, _mem);
    }
}

迭代版本走訪

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
void iteratorTraverse(uint32_t offset, LinearTreeNode *_mem) {
    bool process = true;
    while (offset != -1) {
        LinearTreeNode *node = &_mem[offset];
        if (process) {
            // process
        }
        if (node->childOffsetHead != -1 && process) {
            offset = node->childOffsetHead;
            process = true;
        } else if (node->siblingOffsetNext != -1) {
            offset = node->siblingOffsetNext;
            process = true;
        } else {
            offset = node->parentOffset;
            process = false;
        }
    }
}

節點大小對走訪效能影響

從實驗數據看來，遞迴版本比迭代版本在越大節點情況效能普遍性地好，預估是在遞迴版本造成的快取效能好上許多。

sizeof(LinearTreeNode) bytes \ Traversal Recursive Loop
32 6.049s 5.628s
44 6.651s 6.817s
60 7.460s 6.888s
92 9.361s 9.271s
156 16.844s 16.694s
220 25.294s 27.031s
284 28.181s 30.900s
540 28.560s 33.707s

實驗結果

由於 tick() 效果並不好，調整參數後與原生的作法差距仍然存在，單靠表面積提供的啟發式收縮，效率比加上動態結果還好。但從實驗結果中得知，實作方面存在一些尚未排除的效能障礙，在多核心部分效能差距就非常明顯，預計是在求交點時同步資料造成的 overhead 時間。

而在減少的節點數量，光是啟發是的表面收縮就減少 $25\%$ 節點量，而在動態收縮處理，儘管已經探訪 $10^7$ 收縮點數量近乎微乎其微 (不到一兩百個節點)。

sences \ BVH policy Native Contract(loop) Contract(recursive) Node Reduce
sibenik.pbrt 7.000s 10.502s 9.411s 99576 / 131457
yeahright.pbrt 12.297s 14.638s 14.210s 288707 / 376317
sponza-fog.pbrt 16m36.037s 21m09.960s 20m12.012s 91412 / 121155

Test sences

Reference paper

Ray Specialized Contraction on Bounding Volume Hierarchies, Yan Gu Yong He Guy E. Blelloch

Read More +
學校課程/數位影像生成

pbrt-v2 加速結構 Environment Lights 改寫

Median Cut Algorithm

根據論文 P. Debevec, A Median Cut Algorithm for Light Probe Sampling, SIGGRAPH 2006 Sketches and Applications. 中，預期要將 pbrt/lights/infinite.cpp 中的 class InfiniteAreaLight 用數個點光源取代 Infinite Area Light 的寫法，提升均勻取樣的效能，而 Median Cut Algorithm 在預處理非常快，根據用多少量的點光源將影響品質，若在品質不用太好的 rendering 環境下，這是一個不錯的降質提升速度的方案。

在 Median Cut Algorithm 取樣 64 equal-energy region 的速度在不同取樣個數差異比較

Env-Light.pbrt



Env-Light-new.pbrt



在 Median Cut Algorithm 不同數量 equal-energy region 在 256 samples 速度

Median Cut Algorithm equal-energy region



算法的步驟如下：

  1. 將入射光場影像 (light probe image) 切成好幾個矩形區域，每一個區域將取用一個點光源代替。將入射光場影像轉換成灰階亮度 $Y$，如論文中所提及的方案 $Y = 0.2125 R + 0.7154 G + 0.0721 B$ 這類型的轉換。
  2. 對於每一個區域將會根據最長維度切割成兩個子區域。切割成兩個子區域的亮度總和越接近越好。
  3. 若切割區域數量不到目標的數量，則重複步驟 2。
  4. 最後將每一個區域算出代表的點光源，並加總區域內的亮度和，隨後取樣根據亮度和作為取樣根據 (在 Spectrum MedianCutEnvironmentLight::Sample_L(const Point&, float, const LightSample&, float, Vector&, float*, VisibilityTester*) 中使用)，用每一個區域內部的 pixel 位置和亮度計算重心作為代表的點光源。

算法類似於 k-d Tree，但特別的是每次選擇區域維度最長的那一段進行切割，而不是像 k-d Tree 則是採用輪替維度進行。

Median Cut Algorithm 需要 $\mathcal{O}(HW)$ 時間 $\mathcal{O}(HW)$ 空間來預處理亮度資訊。若要切割成 $N$ 個區域，需要 $\mathcal{O}(\log N)$ 次迭代，每一次迭代增加兩倍區域數量。將一個區域切割時，針對維度最長的那一軸二分搜尋，二分搜尋計算其中一個區域的亮度和是否是整個區域亮度的一半，由於預處理完成的表格，可以在 $\mathcal{O}(1)$ 時間內計算任意平行兩軸的矩形內部元素總和。

維護 sumTable[i][j] 計算左上角 $(0, 0)$ 右下角 $(i, j)$ 的亮度和，計算任意平行兩軸矩形內部和只需要 $\mathcal{O}(1)$ 時間。

1
2
3
4
5
6
7
float getEnergy(float sumTable[], int width, int height) {
    float e = sumTable[VERT(ru, rv)];
    if (lu)	e -= sumTable[VERT(lu-1, rv)];
    if (lv)	e -= sumTable[VERT(ru, lv-1)];
    if (lu && lv)	e += sumTable[VERT(lu-1, lv-1)];
    return e;
}

另一種方案，可以從 pbrt 原生的 class MIPMap 取代上述寫法，MIPMap 的寫法採用分層式的架構，每一層將原圖長寬各縮小一半。計算某矩形元素總和時，藉由兩層的總和估計進行內插。理論上，這種寫法雖然不夠精準，但提供很好的快取優勢。

重心計算

Centroid Formula 1

一般重心計算採用以下公式：

$$X_c = \frac{\sum^{}_{(x, y) \in \mathit{region}} L_{(x, y)} \; x}{\sum^{}_{(x, y) \in \mathit{region}} L_{(x, y)}} \\ Y_c = \frac{\sum^{}_{(x, y) \in \mathit{region}} L_{(x, y)} \; y}{\sum^{}_{(x, y) \in \mathit{region}} L_{(x, y)}}$$

經由 Median-Cut Algorithm 在 Texmap 1 運行後，代表的點光源明顯地偏離亮區，因此為了讓代表的點光源更靠近亮區，我們將其重心公式修改成 Centroid Formula 2。

Centroid Formula 1

若以 $\mathit{Energy} \propto L^2_{(x, y)}$，能量與亮度二次成正比，則計算出的重心會更靠近亮區。

$$X_c = \frac{\sum^{}_{(x, y) \in \mathit{region}} L^2_{(x, y)} \; x}{\sum^{}_{(x, y) \in \mathit{region}} L^2_{(x, y)}} \\ Y_c = \frac{\sum^{}_{(x, y) \in \mathit{region}} L^2_{(x, y)} \; y}{\sum^{}_{(x, y) \in \mathit{region}} L^2_{(x, y)}}$$

比較結果如下：

由於 CSS 問題，詳細內容請到 http://morris821028.github.io/hw-rendering/homework3/submission/index.html 觀看。

Read More +