2015-06-16

b405. 記憶中的序列

Problem

給一個序列，支持插入元素、刪除元素、反轉區間、區間極值、區間總和、區間修改以及最重要的 版本回溯 。

題目輸入會拿上一個答案進行 XOR 的加密，

Sample Input

10
9 81 28 6 40 68 74 50 82 62 
15
1 1 83
2 11
3 7 8
4 4 10 28
5 2 6
85 82 89
14 13 13
56 61
57 57 40
58 61
59 61 63
60 57 59 4
61 58 49
137 139 136
37 38 42

decrypt input

10
9 81 28 6 40 68 74 50 82 62
15
1 1 83
2 11
3 7 8
4 4 10 28
5 2 6
6 1 10
7 4 4
0 5
1 1 16
2 5
3 5 7
4 1 3 60
5 2 9
6 4 7
7 4 8

Sample Output

Solution

持久化 treap，區間反轉，區間最大值，區間最小值，區間總和，插入元素，刪除元素，操作回溯。

這坑爹的題目描述擺明就是拖人下水，只好含淚地寫下去。持久化類型的題目應該是告一段落，就差一個 LCP + HASH 的維護了吧。麻煩的地方在於持久化的懶操作標記傳遞時都要新增加一個節點，這項舉動導致記憶體量非常多，1 秒處理的資訊吃掉 512 MB。二分記憶體上限，終於找到內存池的最多元素個數，中間不小心把多項資訊的 int 型態宣告成 long long，導致內存池一直不夠大。

如果不考慮回溯操作，速度上沒有比暴力算法來得快上很多，treap 基於隨機數的調整，而且每一次修改元素都要增加很多記憶體來補，導致速度上無法領先太多。但如果用暴力法，會造成版本回溯問題，空間會到 $O(
n^2)$，不允許的狀態，持久化造就的平均空間使用$O(Q log N \times sizeof(Node))$，看起來仍然是很巨大的。

#include <bits/stdc++.h>
using namespace std;
#define MAXN 7100000
#define MAXQ 60005
class PersistentTreap {
public:
    struct Node;
    static Node *null;
    struct Node {
        Node *lson, *rson;
        int key, size;
        long long val, mnval, mxval, sumval;
        long long adel; // delivery
        int rdel;
        
        Node(long long c = 0, int s = 1): size(s) {
            lson = rson = null;
            key = rand();
            val = mnval = mxval = sumval = c;
            adel = rdel = 0;
        }
        void update() {
            size = 1;
            size += lson->size + rson->size;
        }
        void pushUp() {
            mnval = mxval = sumval = val;
            if (lson != null) {
                mnval = min(mnval, lson->mnval + lson->adel);
                mxval = max(mxval, lson->mxval + lson->adel);
                sumval += lson->sumval + lson->adel * lson->size;
            }
            if (rson != null) {
                mnval = min(mnval, rson->mnval + rson->adel);
                mxval = max(mxval, rson->mxval + rson->adel);
                sumval += rson->sumval + rson->adel * rson->size;
            }
        }
        void pushDown(PersistentTreap &treap) {
            if (adel) {
                val += adel, mnval += adel, mxval += adel;
                if (lson != null) {
                	lson = treap.cloneNode(lson);
                    lson->adel += adel;
                }
                if (rson != null) {
                	rson = treap.cloneNode(rson);
                    rson->adel += adel;
                }
                adel = 0;
            }
            if (rdel&1) {
                if (lson == null)
                    lson = rson, rson = null;
                else if (rson == null)
                    rson = lson, lson = null;
                else
                    swap(lson, rson);
                
                if (lson != null) {
                	lson = treap.cloneNode(lson);
                    lson->rdel += rdel;
                } 
                if (rson != null) {
                    rson = treap.cloneNode(rson);
                    rson->rdel += rdel;
                }
                rdel = 0;
            }
            pushUp();
        }
    } nodes[MAXN], *root[MAXQ];
    int bufIdx, verIdx;
    
    Node* merge(Node* a, Node* b) {
        if (a == null)	return cloneNode(b);
        if (b == null)	return cloneNode(a);
        if (a != null)	a->pushDown(*this);
    	if (b != null)	b->pushDown(*this);
        Node *ret;
        if (a->key < b->key) {
            ret = cloneNode(a);
            ret->rson = merge(a->rson, b);
        } else {
            ret = cloneNode(b);
            ret->lson = merge(a, b->lson);
        }
        ret->update();
        ret->pushUp();
        return ret;
    }
    void split(Node* a, Node* &l, Node* &r, int n) {
        if (n == 0) {
            l = null, r = cloneNode(a);
            return ;
        } 
        if (a->size <= n) {
            l = cloneNode(a), r = null;
            return ;
        }         
        if (a != null)
            a->pushDown(*this);
        if (a->lson->size >= n) {
            r = cloneNode(a);
            split(a->lson, l, r->lson, n);
            r->update();
            r->pushUp();
        } else {
            l = cloneNode(a);
            split(a->rson, l->rson, r, n - (a->lson->size) - 1);
            l->update();
            l->pushUp();
        }
    }
    
    void insert(Node* &a, Node *ver, int pos, long long s[], int sn) {
        Node *p, *q, *r;
        int n = sn;
        split(ver, p, q, pos);
        build(r, 0, n - 1, s);
        p = merge(p, r);
        a = merge(p, q);
    }
    void erase(Node* &a, Node *ver, int pos, int n) {
        Node *p, *q, *r;
        split(ver, p, q, pos - 1);
        split(q, q, r, n);
        a = merge(p, r);
    }
    void reverse(Node* &a, Node *ver, int left, int right) {
        Node *p, *q, *r;
        split(ver, p, q, left - 1);
        split(q, q, r, right - left + 1);
        q->rdel++;
        a = merge(p, q);
        a = merge(a, r);
    }
    void add(Node* &a, Node *ver, int left, int right, long long val) {
        Node *p, *q, *r;
        split(ver, p, q, left - 1);
        split(q, q, r, right - left + 1);
        q->adel += val;
        a = merge(p, q);
        a = merge(a, r);
    }
    long long q_sum, q_max, q_min;
    void q_dfs(Node *u, int left, int right) {
        left = max(left, 1);
        right = min(right, u->size);
        if (left > right)
            return;        
        if (u != null)
            u->pushDown(*this);
        int lsz = u->lson != null ? u->lson->size : 0;
        int rsz = u->rson != null ? u->rson->size : 0;
        if (left == 1 && right == lsz + rsz + 1) {
            q_sum += u->sumval;
            q_max = max(q_max, u->mxval);
            q_min = min(q_min, u->mnval);
            return ;
        }
        if (left <= lsz+1 && lsz+1 <= right) {
            q_sum += u->val;
            q_max = max(q_max, u->val);
            q_min = min(q_min, u->val);
        }
        
        if (u->lson != null && left <= lsz)
            q_dfs(u->lson, left, right);
        if (u->rson != null && right > lsz+1)
            q_dfs(u->rson, left - lsz - 1, right - lsz - 1);
    }
    void q_init() {
        q_max = LONG_LONG_MIN;
        q_min = LONG_LONG_MAX;
        q_sum = 0;
    }
    long long findMax(Node *ver, int left, int right) {
        q_init();
        q_dfs(ver, left, right);
        return q_max;
    }
    long long findMin(Node *ver, int left, int right) {
        q_init();
        q_dfs(ver, left, right);
        return q_min;
    }
    long long findSum(Node *ver, int left, int right) {
        q_init();
        q_dfs(ver, left, right);
        return q_sum;
    }
    void print(Node *ver) {
        if (ver == null)	return;
        ver->pushDown(*this);
        print(ver->lson);
        printf("[%3lld]", ver->val);
        print(ver->rson);
    }
    void init() {
        bufIdx = verIdx = 0;
        root[verIdx] = null;
    }
private:
    Node* cloneNode(Node* u) {
        Node *ret;
        if (u == null) {
            return u;
        } else {
        	if (bufIdx >= MAXN)
        		exit(0);
            assert(bufIdx < MAXN);
            ret = &nodes[bufIdx++];
            *ret = *u;
            return ret;
        }
    }
    void build(Node* &a, int l, int r, long long s[]) {
        if (l > r)	return ;
        int m = (l + r) /2;
        Node u = Node(s[m]), *p = &u, *q;
        a = cloneNode(p), p = null, q = null;
        build(p, l, m-1, s);
        build(q, m+1, r, s);
        p = merge(p, a);
        a = merge(p, q);
        a->update();
        a->pushUp();
    }
} tree;
PersistentTreap::Node t(0, 0);
PersistentTreap::Node *PersistentTreap::null = &t;
int main() {
    int N, Q;
    long long A[65536];
    long long cmd, x, y, v;
    while (scanf("%d", &N) == 1) {
        for (int i = 1; i <= N; i++)
            scanf("%lld", &A[i]);
        tree.init();
        tree.insert(tree.root[0], tree.null, 0, A+1, N);
        scanf("%d", &Q);
        int verIdx = 0;
        long long encrypt = 0, ret = 0;
        for (int i = 1; i <= Q; i++) {
            scanf("%lld", &cmd);
            cmd ^= encrypt;
            if (cmd == 1) { 		// insert
                scanf("%lld %lld", &x, &v);
                x ^= encrypt, v ^= encrypt;
                long long B[] = {v};
                tree.insert(tree.root[i], tree.root[verIdx], (int) x, B, 1);
                verIdx = i;
            } else if (cmd == 2) { 	// erase
                scanf("%lld", &x);
                x ^= encrypt;
                tree.erase(tree.root[i], tree.root[verIdx], (int) x, 1);
                verIdx = i;
            } else if (cmd == 3) { 	// reverse
                scanf("%lld %lld", &x, &y);
                x ^= encrypt, y ^= encrypt;
                tree.reverse(tree.root[i], tree.root[verIdx], (int) x, (int) y);
                verIdx = i;
            } else if (cmd == 4) { 	// [x, y] += v
                scanf("%lld %lld %lld", &x, &y, &v);
                x ^= encrypt, y ^= encrypt, v ^= encrypt;
                tree.add(tree.root[i], tree.root[verIdx], (int) x, (int) y, v);
                verIdx = i;
            } else if (cmd == 5) {  // max(A[x:y])
                scanf("%lld %lld", &x, &y);
                x ^= encrypt, y ^= encrypt;
                ret = tree.findMax(tree.root[verIdx], (int) x, (int) y);
                printf("%lld\n", ret);
                encrypt = ret;
                tree.root[i] = tree.root[verIdx];
                verIdx = i;
            } else if (cmd == 6) {  // min(A[x:y])
                scanf("%lld %lld", &x, &y);
                x ^= encrypt, y ^= encrypt;
                ret = tree.findMin(tree.root[verIdx], (int) x, (int) y);
                printf("%lld\n", ret);
                encrypt = ret;
                tree.root[i] = tree.root[verIdx];
                verIdx = i;
            } else if (cmd == 7) {  // sum(A[x:y])
                scanf("%lld %lld", &x, &y);
                x ^= encrypt, y ^= encrypt; 
                ret = tree.findSum(tree.root[verIdx], (int) x, (int) y);
                printf("%lld\n", ret);
                encrypt = ret;
                tree.root[i] = tree.root[verIdx];
                verIdx = i;
            } else if (cmd == 0) {  // back Day x-th
                scanf("%lld", &x);
                x ^= encrypt;
                tree.root[i] = tree.root[x];
                verIdx = i;
            } else {
                puts("WTF");
                return 0;
            }
        }
    }
    return 0;
}

Test Code

暴力法驗證用，不支持加密的輸入運行，用一個 vector 完成。

#include <bits/stdc++.h>
using namespace std;
#define MAXN (1<<19)
#define MAXQ 65536
int main() {
    int N, Q;
    long long A[65536];
    long long cmd, x, y, v;
    while (scanf("%d", &N) == 1) {
        for (int i = 1; i <= N; i++)
            scanf("%lld", &A[i]);
        vector< vector<long long> > L;
        L.push_back(vector<long long>(A+1, A+1+N));
        scanf("%d", &Q);
        int verIdx = 0;
        long long encrypt = 0, ret = 0;
        for (int i = 1; i <= Q; i++) {
            scanf("%lld", &cmd);
            vector<long long> TA = L[verIdx];
            if (cmd == 1) {         // insert
                scanf("%lld %lld", &x, &v);
                TA.insert(TA.begin()+x, v);
            } else if (cmd == 2) {  // erase
                scanf("%lld", &x);
                TA.erase(TA.begin()+x-1);
            } else if (cmd == 3) {  // reverse
                scanf("%lld %lld", &x, &y);
                reverse(TA.begin()+x-1, TA.begin()+y);
            } else if (cmd == 4) {  // [x, y] += v
                scanf("%lld %lld %lld", &x, &y, &v);
                for (int i = (int) x-1; i < y; i++)
                    TA[i] += v;
            } else if (cmd == 5) {  // max(A[x:y])
                scanf("%lld %lld", &x, &y);
                long long ret = LONG_LONG_MIN;
                for (int i = (int) x-1; i < y; i++)
                    ret = max(ret, TA[i]);
                printf("%lld\n", ret);
            } else if (cmd == 6) {  // min(A[x:y])
                scanf("%lld %lld", &x, &y);
                long long ret = LONG_LONG_MAX;
                for (int i = (int) x-1; i < y; i++)
                    ret = min(ret, TA[i]);
                printf("%lld\n", ret);
            } else if (cmd == 7) {  // sum(A[x:y])
                scanf("%lld %lld", &x, &y);
                long long ret = 0;
                for (int i = (int) x-1; i < y; i++)
                    ret += TA[i];
                printf("%lld\n", ret);
            } else if (cmd == 0) {  // back Day x-th
                scanf("%lld", &x);
                TA = L[x];
            } else {
                puts("WTF");
                return 0;
            }
            L.push_back(TA), verIdx = i;
            printf("version %d\n", i);
            for (auto &x: TA)
                printf("[%3lld]", x);
            puts("");
        }
    }
    return 0;
}

Test Generator

小測資檢驗用

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <algorithm>
#include <set>
using namespace std;
int main() {
     srand (time(NULL));
    int testcase, n, m, x, y;
    int A[100];
    testcase = 1;
    freopen("in.txt", "w+t", stdout);
//    printf("%d\n", testcase);
    while (testcase--) {
        n = 32;
        printf("%d\n", n);
        for (int i = 0; i < n; i++) {
            printf("%d%c", rand(), i == n-1 ? '\n' : ' ');
        }
        m = 1000;
        printf("%d\n", m);
        vector<int> L;
        L.push_back(n);
        for (int i = 1; i <= m; i++) {
            while (true) {
                int cmd = rand()%8;
                if (n == 1 && cmd == 2)
                    continue;
                if (cmd == 0) {
                	int x = rand()%i;
                	printf("%d %d\n", cmd, x);
                	n = L[x];
                }
                if (cmd == 1)
                    printf("%d %d %d\n", cmd, rand()%(n+1), rand()), n++;
                if (cmd == 2)
                    printf("%d %d\n", cmd, rand()%n+1), n--;
                if (cmd == 3) {
                    x = rand()%n+1, y = rand()%n+1;
                    if (x > y)  swap(x, y);
                    printf("%d %d %d\n", cmd, x, y);
                }
                if (cmd == 4) {
                    x = rand()%n+1, y = rand()%n+1;
                    if (x > y)  swap(x, y);
                    printf("%d %d %d %d\n", cmd, x, y, rand());
                }
                if (cmd == 5 || cmd == 6 || cmd == 7) {
                    x = rand()%n+1, y = rand()%n+1;
                    if (x > y)  swap(x, y);
                    printf("%d %d %d\n", cmd, x, y);
                }
                break;
            }
            L.push_back(n);
        }
    }
    return 0;
}

Read More +

2015-06-13

學校課程/虛擬實境

虛擬實境建模 SVD 筆記

由於在選讀論文時，和上課都遇到相關問題，但一直對這 SVD 的實際感觸沒很深，了解數學計算，卻不知道特性一直很困惑我。在這裡收集一下曾經學過的內容，當初在線性代數時沒有學到這一塊，至於 $LU$ 分解有沒有上到過，我自己也深感困惑，不過看一下網路資源勉強還能理解。

至於 SVD 則在線性代數非常重要，在許多領域時常會遇到，如巨量資料 (massive data)？機器學習？影像處理？就我所知理解看看吧！

SVD

定義從一般式

$$A_{m, n} = U_{m, m} \Sigma_{m, r} V_{n, r}^T$$

得到縮減過後的相似矩陣

$$A_{m, n} \cong U_{m, r} \Sigma_{r, r} V_{n, r}^T$$

給定 $A$ 矩陣，要變成三個矩陣相乘，可以壓縮儲存空間，一般來說都是拿來得到近似的 $A$ 矩陣，藉由刪除過小的 eigenvalue 來完成，eigenvalue 就是所謂的矩陣特徵值。

在學線性代數時，沒有仔細想過 eigenvalue 的意義所在，每一個 eigenvalue 也會對應一個 eigenvector，這些 eigenvalue 大小對應的分佈的相關性大小，而 eigenvector 則表示新的轉換軸的方向。類似最小平方法得到的結果。

但在巨量資料課程中我們可以知道一件事情，若是稀疏矩陣 $A$，經過 SVD 運作後，得到 $U, V$ 矩陣都是稠密的，因此空間儲存反而划不來，因此會有另外一種解法 CUR，可以上網搜尋 massive data mining coursera 在這裏可以知道關於 CUR 的資訊。

從課程中的例子來看，假設在電影評論系統中，則 A 矩陣表示電影對人的喜好度，則 SVD 相當於如下描述 (假想，單純的矩陣概念)

$A(i, j)$ 電影編號 i 對用戶 j 的好感度
$U(i, j)$ 電影編號 i 對電影屬性 j 的成分程度
$\Sigma(i, j)$ 電影屬性 i 跟電影屬性 j 的之間重要性。
$V(i, j)$ 表示用戶 i 對電影屬性 j 的偏愛程度。

也就是說有一些電影屬性不是很重要，可以移除掉，這樣的表示法得到一個近似的矩陣 $A$。

在巨量資料中一個有趣的運用，可以用在購物推薦系統中，來達到預測相關性，跟某人對其他產品的喜好程度或者是需求性，但是在處理之前會遇到矩陣過大，但又非常稀疏，為了解決紀錄問題，就有奇奇怪怪的 Dimensionality Reduction 來完成，SVD 就是其中一種方法。

其他應用部分，將不重要的特徵值消除，來求得哪些元素是需要的，看到有一篇文章在做 PCA 主成份分析。

$$\begin{align} & A_{m, n} V_{n, r} = U_{m, r} \Sigma_{r, r} V_{n, r}^T V_{n, r} \\ & \because V_{n, r}^T V_{n, r} = I \\ & \therefore A_{m, n} V_{n, r} = U_{m, r} \Sigma_{r, r} \\ \end{align}$$

如此一來就可以知道要保留哪幾個主要資訊。

在虛擬實境論文中也是有這樣的處理，做的是 CT 建模，將模型的每個頂點，要計算 QEM (最小二次誤差) 的重要性，可以將其表示成一個矩陣，利用 SVD 方式得到，來決策要保留前 k 個重要性的節點，再將其他節點移除，來讓新得到節點盡可能地可以貼近真實模型。

之所以會這麼做，是要利用內插法來強迫得到邊緣偵測，但這樣會造成很多新頂點，為了要把過多頂點消除，使用 QEM 作為評價方法，利用 SVD 來完成大規模刪除的操作。

Map-Reduce

都明白 SVD 是一個很繁複的做法，有一個 $O(N^3)$ 的算法來求得所有 eigenvalue，若要用一般電腦跑，光是得到一個大矩陣所有 eigenvalue 複雜度相當高，其複雜程度沒有細讀過。有一種迭代的方法，類似馬可夫過程，每一次迭代可以逼近一個 eigenvalue，為了得到 n 個 eigenvalue，想必要迭代非常多次。運用 Map-Reduce 分散到很多台群落電腦加快矩陣乘法，這種方法的複雜度是可以期待的。

詳細操作請參考

eigenvalue $\lambda$ 滿足

$Av = \lambda v$

現在開始迭代 $v$ 讓其等式趨近於穩定

$A x_{i} = \lambda x_{i+1}$

初始化$x_{i} = (1, 1, ..., 1, 1)$ ，其中$|x_{i}| = 1$ 藉由長度為 1 的向量，來求得$\lambda$，之所以要對 $x$ 進行正規化，其一就是要求出$\lambda$，另一個原因是要看何時穩定，當$|x_{i} - x_{i+1}|$ 差距夠小時，則表示進入穩定狀態。

那就可以得到矩陣 $A$ 的其中一個 eigenvalue，為了求得另外一組解，得到新的$A^* = A - \lambda x x^T$，再對$A^*$ 找到 eigenvalue，不斷地迭代下去。

Read More +

2015-06-12

手札日記

閱讀人生『閱讀』社群外出紀錄

前處理

在這個畢業季節，考量著那做不到服務小時的畢業門檻，近乎危急之下，學校邀請閱讀社群創辦人來講這個社群的心路歷程與目標，被家裡人催著參加活動，剛好參加了這個講座。

在進入這場講座之前，用很久的 Facebook 都不曉得有『閱讀』社群的存在，活動標榜「百萬人都愛上的華人最大臉書閱讀社群」的確有朋友也喜歡這個粉絲頁，從未見過演算法會讓我知道這個訊息，稍微看了一下內容，也許就跟長輩們常常分享的心靈雞湯差不多吧？既然參加活動，事前準備按個讚也不為過，但一個星期中也從來都沒見過相關訊息跳出來。

好吧！參加活動就能開始理解到底是賣什麼。

《吹響吧！低音號》誰知道呢那個人比較特別

運行

入場

在這炙熱的畢業季，走到學校後門的白色貨櫃屋，比起以往只從租屋處到學校的田中小徑，這一天走了一條不尋常的道路，不習慣的環境，停個車也擔憂打擾到其他住在附近的人。曾經聽同學說這是個茶鋪，對我來說是高檔了點，只要有網路電腦的地方，在哪都無所謂，這種高檔店一年去的次數一隻手就能數出來。

吱吱悟悟走向協辦單位拿取資料填寫、簽名報到，很久沒聽到認識人以外的聲音，聽不懂活動人員給予我的下一個指令是什麼，傻里傻氣地按照當初在網路上看的活動訊息行事，心想找知道就找個朋友跟我一起來，但又想到朋友常回應的那一句話「我時數滿了，為什麼要去？」看來是沒望了。

也許在別人眼中舉止是多麼笨拙，選擇一杯飲料作為低銷入場、再找個位子坐下，久違的一步，對我來說可是比寫代碼還辛苦。為了防止看起來突兀，盡可能地選擇一個可以躲藏的位子，所思又想「到底哪個位子才適合我，這一切好困難」在這小茶鋪，擺著如電視劇一樣的一堆小椅子、小投影幕及少數的桌子，既不是可以躲藏自我的專業講廳，也不是階梯式的椅子擺放，彼此跟彼此之間的距離少了點隔閡，心中的畏懼又增添一層。

《約會、戀愛究竟是什麼》可能無法接受自己連普通人都做不了的現實吧！

記憶讀取

活動一開始，從主持人 鄭俊德 口中開始道出自己的經驗、家庭背景，與其款款地道出自己有多成功、怎麼成功的，那些不具有說服別人跟著一起做的舉動，介紹自己的一些小事開始，來龍去脈講的是自己怎麼從家庭圈圈中跳出來，又說著自己如何從學歷下走出自己。

插曲

中間有個小插曲活動「 如何說出自己的故事？ 」先別急著講長篇大論，能不能用一分鐘介紹給旁邊的人認識，現在的我們越來越少用文字對談，資訊相當多導致大部分的時間都在聽別人怎麼做，如何說給別人聽變得相當困難。

要我們介紹自己給自己旁邊的人，「等等，我身旁坐得是 … ？」此時心中不斷地迴盪，既然都走出門總得再跨一步『我是 …』我想已經在別人眼中死了吧，自覺說出來會被當作活在平行世界似的。跟妹子聊天什麼的，不是奇蹟可以形容的，但想到自身的負面情緒會影響他人，說話又更加的怯步。

「我叫翔雲，資工系，來這邊的機緣只是單純畢業門檻，平常是個碼農。中興轉學來到這裡，在別人眼中只是離家鄉花蓮更近一點 …」

話語就此打斷，對方是個法文系女生，跟我仇視的英文有點血緣關係，大一生不刷 Facebook，覺得在上面會消耗時間太多，大一生活多采多姿又忙碌，這一點我是很清楚的，曾經也過著自以為忙碌的代碼生活。原來她來這裡的原因是因為主持人是位 專欄作家 ，看來是一位社會中現實真正的生活者，看著雜誌、與朋友交談、或者在圖書館中靜靜讀書度過一日的生活的吧，在我心中假想的形象不斷地浮現。不知道能不能說個一分鐘，更多情況無法進行即時的對談，想說的話很多，很想讓別人記住我，但又怕這一段成為 黑歷史 ，躊躇之際又選擇放棄說話。

《果然我的青春戀愛搞錯了》這世界上有著只要在場就會讓氣氛變差的人

圈圈

主持人講到，對於父母給予的環境到底要怎麼跨出圈圈，圈圈中的 熱情、專業、經濟 如何永恆運作？在圈圈中發現世界，還是在圈圈中活在自我 (用他的弟弟來講例子，弟弟好像被出賣了！)，如何將圈圈放大？理解的越多，圈圈越大就會越複雜。用熱情推向專業再來獲得金錢，空有熱情造就不出專業，將自己的熱情變成專業，但更多的是沒有熱情，先用專業帶來的經濟，將獲得的金錢培養出熱情，再利用熱情跳到另一個新的圈圈。

從某些人身上的確發現這一點，不少人在大學習得專業卻沒有熱情，偶爾會罵著他們「到底是不是來學習，還是來玩的？」在知道自己讀的專業也不過是生活的一小塊時，那種想法也隨之消去，在大學累積人脈才是重要的，學業只要維持在最低限度之上就好。

一個問題被創造出來通常是有解法的，即便那解法不夠快速，講者一講到這個部分，身為計算機專業可說是反對再三，有些問題連解法都沒有，這裡的沒有可以說是人一生的時間也算不出答案的問題。這裡笑笑就好，畢竟現實生活中這種問題鮮少。

既然自己無法解決，肯定要問問朋友，朋友也可以問問朋友，相互幫忙是不錯的選擇。「 世界那麼大，理解得卻很少。 」腦袋記憶的時間有限，學習專研也是有所限制，什麼都要理解是不可能的。既然現在流行分散式系統架構，單機運作什麼的，隨著學習曲線不斷地陡增，越來越明白這困境。「 教練，我想會一個交朋友的技能啊！ 」

《果然我的青春戀愛搞錯了》想要完全理解是一個非常自以為是、獨裁又傲慢的願望

習得性無助

關於習得性無助可以參考維基百科。

講者用實際實驗的例子告訴我們，實驗中把一條狗彼得關進一個地板接有電擊房間中，當狗尋找躲避電擊的方法，若發生找不到解決方案，彼得會進入放棄狀態，直接接受地板上的電擊，毫無抵抗地趴下。

接著把彼得帶到另一個房間與其他兩隻狗放在一起，這個房間與剛才只差了一點，跨過某一個衡桿後，另一邊是沒有電擊的，另外兩隻狗開始掙扎逃脫，其中一隻狗恰好跨過衡桿躲避電擊，另一隻狗也就跟著那隻狗一起在衡桿的另一側，特別的是彼得 選擇放棄 ，仍然在電擊區承受電擊，這就是所謂的 習得性無助 。

講者描述到很多人接觸數學是多麼地恐懼，也許我面對的是對 英文的恐懼 ，因為聽不懂發音，好像也沒有努力的必要，於是我選擇放棄，一般人認為理所當然的事情，做不到那不是更容易放棄嗎？其實一開始英文成績都還算不錯，國中基測還是英文滿分，隨著聽力成績的壓迫，近乎猜測期望值的分數，讓我越來越不行。

《吹響吧！低音號》不想和其他傢伙一樣

於是我走向另一邊，每天窩著打程序爬到現在，解題數上千、文章篇數上千、累積點閱次數一百多萬。現在的我相當無助，畢不了業，英文不好甚至沒有工作、未來，即便從他人口中道出「Morris 很強」，對我來說不具有一絲的鼓勵，更明白前面有更多無法理解的內容， 世界都用英文封藏起來 。

自殺

這社群有一點很有趣，關於文章分享的部分，提供不只有大眾愛看的內容，原以為只是喜愛讚數的社群，對於負面情緒的文章從來不在意，講者口中道出「 用正向文章鼓勵負面情緒的人們是辦不到的 」因此他們仍然會放一些讚數不高的文章，不總是討好大眾愛看的內容，比起以往粉絲頁的經營方式這一點令我相當訝異。

這一點我很喜歡，也許正是處於負面那一群的人們。

甚至會因為有人看了那些正面文章，向他們哭訴要自殺了呢，要是我頂多到憎恨的地步罷了，「你們口口說聲簡單、正向的事物，這些我都沒有，是不是該結束人生呢？」被這樣的話語絆住雙腳。

講者描述「他們想要自殺，目標不是結束生命，而是解決問題」聽到 解決問題 四個字耳目一亮，在平常用程序解決算法問題正是解決問題，對於解決問題的熱情，自殺真是不錯的方法，時間複雜度低、速度快速、實作方法簡單，不外乎是這世界快好的 $O(1)$ 作法。

彩蛋

《吹響吧！低音號》我會努力的，所以你也要好好努力

《吹響吧！低音號》那時突然想到必須和前輩們競爭的現實，令我感到害怕

《吹響吧！低音號》不過，還想做得更好

《獻給阿爾吉儂的花束》大概重複一百萬次吧，因為我有大把的時間，反正我沒有女朋友、朋友

《果然我的青春戀愛搞錯了》我也想「真實」

《吹響吧！低音號》快住口

Read More +

2015-06-12

學校課程/資訊安全

近代加密數位簽章 DSA 筆記

DSA

Digital Signature Standard (DSS) 數位憑證標準，而有一個 DSA 演算法來完成憑證。

每一次憑證傳送三個參數 $(r, s, M)$，利用私用金鑰進行加密。讓對方使用公鑰進行驗證。

$r = (g^k \mod p) \mod q$
$s = k^{-1}(SHA(M) + x \times r) \mod q$

give receiver $(r, s, M)$, we have random (one-time) $k$ & private key $x$

每一次簽名時，最好使用隨機亂數 $k$ 作為簽章用途，隨用即棄，若發生重複使用，將會發生被攻擊，被攻擊者拿到私鑰後，攻擊者就能進行偽裝。

Attack

if reused $k$

如果重複使用隨機變數 $k$，可以藉由夾擊的方式得到，雖然不知道加密的隨機亂數 $k$ 究竟是什麼，藉由夾擊可以把 $k$ 消除，利用乘法反元素反求 $x$ 是多少，公式如下所述。

$k = s1^{-1} (SHA(M1) + x \times r) \mod q$ - (1)
$k = s2^{-1} (SHA(M2) + x \times r) \mod q$ - (2)

subtract (1) - (2)

$0 = s1^{-1} SHA(M1) - s2^{-1} SHA(M2) + (s1^{-1} - s2^{-1}) x \times r \mod q$
$s1^{-1} SHA(M1) - s2^{-1} SHA(M2) = (s2^{-1} - s1^{-1}) x \times r \mod q$
$x = (s1^{-1} SHA(M1) - s2^{-1} SHA(M2)) \times (s2^{-1} - s1^{-1})^{-1} \times r^{-1} \mod q$

get private key. Dangerous !

Verify

驗證的正確性與否，接收者需要計算以下四個步驟來完成驗證。

$w = s^{-1} \mod q$
$u1 = (SHA(M) \times w) \mod q$
$u2 = r \times w \mod q$
$v = g^{u1} y^{u2} \mod q$

test $v = r$, why correct ?

理論的正確性可以參考以下證明，把公開的 $y = g^x \mod p$ 拿出來，這個安全性取決於離散對數的算法難度。

known $y = g^x \mod p$, $x$ is a private key
because $s = k^{-1}(SHA(M) + x \times r) \mod q$
then $k = s^{-1}(SHA(M) + x \times r) \mod q$

want $r = (g^k \mod p) \mod q$

$r = (g^k \mod p) \mod q = g^{s^{-1}(SHA(M) + x \times r)} \mod q$
$r = g^{s^{-1} SHA(M)} \times g^{x (s^{-1} \times r)} \mod q$
$r = g^{s^{-1} SHA(M)} \times g^{x^{(s^{-1} \times r)}} \mod q$
$r = g^{s^{-1} SHA(M)} \times y^{(s^{-1} \times r)} \mod q$

而我們事實上就是在處理指數部分，

離散對數

解決問題 $y = g^x \mod p$，當已知 $y, g, p$，要解出 $x$ 的難度大為提升，不像國高中學的指數計算，可以藉由 log() 運算來完成，離散對數可以的複雜度相當高，當 $p$ 是一個相當大的整數時，通常會取用到 256 bits 以上，複雜度則會在 $O(2^{100})$ 以上。

實際上有一個有趣的算法 body-step, giant-step algorithm，中文翻譯為 小步大步算法 ，在 ACM-ICPC 競賽中也可以找到一些題目來玩玩，算法的時間複雜度為 $O(\sqrt{p})$，空間複雜度也是 $O(\sqrt{p})$。相信除了這個外，還有更好的算法可以完成。

小步大步算法其實很類似塊狀表的算法，分塊處理，每一塊的大小為 $\sqrt{p}$，為了找尋答案計算其中一塊的所有資訊，每一塊就是一小步，接著就是利用數學運算，拉動數線，把這一塊往前推動 (或者反過來把目標搜尋結果相對往塊的地方推動)。因此需要走 $\sqrt{p}$ 大步完成。

Read More +

2015-06-06

學校課程/虛擬實境

虛擬實境考古筆記

虛擬實境這一門課可以讓你體驗有如在上另一門課的虛擬實境，前半段講生物結構，來了解軟硬體需要那些特性、效能，後半段在講視覺建模算法，牽涉到計算幾何、線性代數 … 等課程。

在幾何課程中介紹的 Delaunay triangulation 對於建模算法是相當重要，接著在多重解析度部分牽涉到傳輸跟處理速度，在資料結構的儲存上很有趣。為了加速運算效能、漸進式需求，各種轉換壓縮都會跑出來。

複習

(15%)[Definition]
詳細說明一個虛擬實境 (virtual reality) 系統需要具備那些特性，並說明一個虛擬實境系統要用到哪些軟體技術。

提供人們能感知的動態環境，視覺、嗅覺、觸覺，並且能與環境互動中得到反饋，得到逼真感覺。需要以下幾個軟體
1. 影像處理與繪製
2. 物理引擎
3. 擷取資料的整合
(15%)[Human perception]
說明人類 1. 單眼立體視覺的原理，2. 雙眼立體視覺的原理；並3. 比較其優缺點。
1. 單眼立體視覺靠單一眼球調整焦距，將數張影像交給大腦來整合，這必須藉由對物體認知來補足。
2. 雙眼立體視覺的原理主要靠左右兩個眼球不同的相位差，來計算立體影像的距離差距。
3. 視野範圍以及反應時間，單眼必須藉由調整焦距來補足，雙眼可以同時獲取不同相位差的影像。
(15%)[Hardware]
說明感應手套 (sensing glove) 與力回饋手套 (force-feedback glove) 的原理，並比較他們在應用上的差別(可舉例說明)。

感應手套是由操作者的反應傳輸給感測器，藉由手套上感應器之間的變化數據得到操作者給予的狀態。力回饋手套則是將物理特性反饋給使用者，藉由手套上的機械，施予使用者力反饋。

感應手套為使用者操作，力反饋手套則相反由環境操作。
(15%)[Modeling]
說明五種 3D 模型的表面表示式 (surface representation)，並比較這些表示式的優缺點。

目前共計主要分成五種表示法：
- regular square grid meshes (RSG)
- triangulated regular meshes (TRM)
- triangulated irregular networks (TINs)
- polygon meshes
- parametric patches
  
  RSG 建造簡單，但會有四點不共平面，影響到貼圖繪製上的問題。TRM 解決四點不共平面問題，但與 RSG 問題無法凸顯邊緣。TINs 解決無法凸顯邊緣問題，但需要經過最佳化算法來得到不規則三角網格，儲存方式必須記錄所有的點座標，由於要經過最佳化算法，處理速度會慢上很多。動態更動座標的效果較差。polygon meshes 建造相對困難，多邊形的最小單位是三角形，頂點儲存量比 TINs 少很多，parametric patches 藉由取樣點套入參數，讓表面接近平滑，針對部分訊息的改造，平滑效果造成的反饋比較逼真。
(15%)[Multiresolution modeling]
說明並比較以小波轉換為基礎的多重解析度模塑 (wavelet-based multiresolution modeling) 技術與 Lindstrom 等人的對稱三角網格 (symmetric triangulated networks) 多重解析度模塑技術的原理及優缺點。
- wavelet-based multiresolution modeling
  解析度資料量與面積成正比，考慮到與觀察者之間的距離，用另外三張影響轉換到更高解析度，每一次解析度增加兩倍。
- symmetric triangulated networks
  考慮到地形高地、平坦性對於觀察者的影響，投影到使用者的平面，判斷是否差異過大需要提供更高解析，解析度資料量所需節點數成正比
  
  相較之下，小波轉換給予的傳輸資料反應會比較慢，同時也會在使用者感受不到差異的地方提供更高解析度，很容易浪費資源，對稱三角網格在資料鏈結會複雜許多，沒有小波轉換的儲存架構來得容易理解，但對稱三角網格反應速度快，支持快速的解析度變換。
(15%)[Multiresolution modeling]
說明漸進網格 (progressive mesh) 多重解析度模塑與二次誤差 (quadric-error-metric) 多重解析度模塑的相同與差異處，並比較兩者的優缺點。

漸進網格藉由拔點、縮邊，提供不同解析度，依序拔點，考慮四個因素來決策拔點順序，這四個因素分別為影響到其他點之間的距離、縮邊的長度總和、頂點材質的差異、對邊界影響程度。

考慮的是從精細模型轉換到粗糙模型的差異變化，挑選變成粗糙的影響最小的先拔除，這樣的計算量會比較大，拔除一點後，鄰近的點的誤差值計算也要跟著重新計算。

二次誤差類似漸進網格的方式進行拔點縮點，但縮邊完的的結果不會只有三種可能，縮到兩端點的其中一個、或者是中點，提供一個數學計算來找到縮到哪一個地方更好。

二次誤差考慮的是從端點誤差值得相加，並且產生一個新的頂點，有點霍夫曼編碼的意味在。計算從粗糙模型上的點到精細模型平面的距離平方，與之前的模型 (Progressive mesh - Hoppe) 恰恰相反，並不是考慮精細模型上的點到粗糙模型平面的距離平方。因此對於 QEM 的誤差可以累加。問題是粗糙模型上的頂點是無法得知，屬於未知數，用內插法。QEM 誤差精確度只是約略，速度快但品質沒有 Progressive mesh 來得好。

(15%)[Multiresolution modeling]
說明 (a) 動態載入 (dynamic loading) 的意義。(b) 動態載入與多重解析度模塑技術整合會產生什麼問題。(c)解決上述問題的三種方法。

動態載入為在記憶體不足時，將部分資料存放在磁碟中，當程式運行到需要那些資料進行運算時再從磁碟中讀取，這樣的過程稱為動態載入。

動態載入對於多重解析度塑模來說，主要是記憶體問題以及所需要的解析度跟數據大小關係，解析度度越高，所需要的資料量越多。但是畫面呈現的解析度不一下，要怎麼進行資料讀取，動態載入的數據結構定義則需要特別設計，例如是要 response time 最小、計算複雜度最低，讀取數量最小 … 等。

假設需要不同解析度可以靠 Discrete layered multiresolution models，每一種解析度彼此獨立儲存，這樣的方法可以提供基礎的不同解析度，但讀取時間相對長，儲存在磁碟的空間需求大，效果呈現是最快的，但是離散結果會有不連續的解析度變化。

為了解決解析度連續變化問題，可以使用 progressive mesh 來達到，提供以點、線與解析度呈現線性變化。

若需要不同解析度，可以利用 hierarchical models 來完成畫面中每一處的解析度不同，這是一種資料結構，當需要更高解析度時，則走訪更深的樹節點來印出所需要的畫面。

評價效果好壞利用 Delaunay triangulation 計算平滑程度。
(15%)[Rendering]
說明shading的意義，敘述 Gouraud shading 與 Phong shading 的作法，最後比較他們的優缺點。

Read More +

2015-06-05

解題區/解題區 - POJ

POJ 1151 - Atlantis

Problem

給一堆平行兩軸的矩形，請問面積的聯集為何。

Sample Input

2
10 10 20 20
15 15 25 25.5
0

Sample Output

1 2	Test case #1 Total explored area: 180.00

Solution

利用掃描線算法，從 x 軸小掃描到大，維度 y 軸上的覆蓋情況。

若單純用離散方法在 2D 網格上填充，算法複雜度是 $O(n^2)$，最慘情況還會再增加到 $O(n^4)$，平均上來說離散方法是可行，但是記憶體是 $O(n^2)$。

由於 ITSA 的那題 n 會高達 30000，於是拿這一題來練手，掃描線算法複雜度 $O(n \log n)$，利用線段樹操作時，維護某個區間的完全覆蓋數，當完全覆蓋數等於 0 時，則計算子樹的覆蓋長度總和，雖然沒辦法完全覆蓋，則部分覆蓋可以從子樹計算得到。

#include <stdio.h> 
#include <string.h>
#include <map>
#include <vector>
#include <math.h>
#include <algorithm>
using namespace std;
const int MAXN = 262144;
class SegmentTree {
public:
    struct Node {
        int cover;				// #cover
        double L, R;			// value in real line, cover [L, R]
        double clen;			// cover length
        void init(int a = 0, double b = 0, double c = 0, double d = 0) {
            cover = a, L = b, R = c, clen = d;
        }
    } nodes[MAXN];
    double Y[MAXN];
    void pushDown(int k, int l, int r) {
    }
    void pushUp(int k, int l, int r) {
        if (nodes[k].cover > 0) {
            nodes[k].clen = nodes[k].R - nodes[k].L;
        } else {
            if (l == r)
                nodes[k].clen = 0;
            else
                nodes[k].clen = nodes[k<<1].clen + nodes[k<<1|1].clen;
        }
    }
    void build(int k, int l, int r) { 
        nodes[k].init(0, Y[l], Y[r+1], 0);
        if (l == r)
            return ;
        int mid = (l + r)/2;
        build(k<<1, l, mid);
        build(k<<1|1, mid+1, r);
        pushUp(k, l, r);
    } 
    // operator, add
    void addUpdate(int k, int l, int r, int val) {
        nodes[k].cover += val;
        if (nodes[k].cover > 0) {
            nodes[k].clen = nodes[k].R - nodes[k].L;
        } else {
            if (l == r)
                nodes[k].clen = 0;
            else
                nodes[k].clen = nodes[k<<1].clen + nodes[k<<1|1].clen;
        }
    }
    void add(int k, int l, int r, int x, int y, int val) {
        if (x <= l && r <= y) {
            addUpdate(k, l, r, val);
            return;
        }
        pushDown(k, l, r);
        int mid = (l + r)/2;
        if (x <= mid)
            add(k<<1, l, mid, x, y, val);
        if (y > mid)
            add(k<<1|1, mid+1, r, x, y, val);
        pushUp(k, l, r);
    }
    // query 
    double r_sum;
    void qinit() {
        r_sum = 0;
    }
    void query(int k, int l, int r, int x, int y) {
        if (x <= l && r <= y) {
            r_sum += nodes[k].clen;
            return ;
        }
        pushDown(k, l, r);
        int mid = (l + r)/2;
        if (x <= mid)
            query(k<<1, l, mid, x, y);
        if (y > mid)
            query(k<<1|1, mid+1, r, x, y);
        pushUp(k, l, r);
    }
} tree;
struct Event {
    double x, yl, yr;
    int val;
    Event(double a = 0, double b = 0, double c = 0, int d = 0):
    x(a), yl(b), yr(c), val(d) {}
    bool operator<(const Event &a) const {
        return x < a.x;
    }
};
double Lx[32767], Ly[32767], Rx[32767], Ry[32767], K[32767];
int N;
double areaCoverAll() {
    vector<Event> events;
    vector<double> VY;
    map<double, int> RY;
    for (int i = 0; i < N; i++) {
        double x1, x2, y1, y2;
        x1 = Lx[i], x2 = Rx[i];
        y1 = Ly[i], y2 = Ry[i];
        events.push_back(Event(x1, y1, y2,  1));
        events.push_back(Event(x2, y1, y2, -1));
        VY.push_back(y1), VY.push_back(y2);
    }
    
    sort(events.begin(), events.end());
    sort(VY.begin(), VY.end());	
    VY.resize(unique(VY.begin(), VY.end()) - VY.begin());
    
    for (int i = 0; i < VY.size(); i++) {
        tree.Y[i] = VY[i];
        RY[VY[i]] = i;
    }
    
    tree.build(1, 0, VY.size()-2);
    
    double area = 0, prevX = 0;
    for (int i = 0, j; i < events.size(); ) {		
        if (i > 0) {
            tree.qinit();
            tree.query(1, 0, VY.size()-2, 0, VY.size()-2);
            area += (events[i].x - prevX) * tree.r_sum;
        }
        j = i;
        while (j < events.size() && events[j].x <= events[i].x) {
            tree.add(1, 0, VY.size()-2, RY[events[j].yl], RY[events[j].yr] - 1, events[j].val);
            j++;
        }
        prevX = events[i].x;
        i = j;
    }
    return area;
}
int main() {
    int testcase, cases = 0;
    while (scanf("%d", &N) == 1 && N) {
        for (int i = 0; i < N; i++)
            scanf("%lf %lf %lf %lf", &Lx[i], &Ly[i], &Rx[i], &Ry[i]);
        printf("Test case #%d\n", ++cases);
        printf("Total explored area: %.2f\n", areaCoverAll());
        puts("");
    }
    return 0;
}
/*
2
10 10 20 20
15 15 25 25.5
*/

Read More +

2015-06-05

解題區/出題解題

[ACM 題目] 障礙物轉換

Problem

在計算幾何領域中，機器人規劃路線是一個實用的議題，由於機器人是一個有體積的物體，要避開其他的障礙物，障礙物碰撞計算會變得相當複雜。由於某 M 沒有學好這一部分，知道一種方法可以將地圖轉換，把機器人轉換成一個點，計算得到新的障礙物，如此一來就不用處理障礙物碰撞。

從上圖中，一個箭頭型機器人要從左下角移動到右上角，假設不考慮物體可以旋轉，請問是否能移動過去。若把機器人當作一點看待，則必須將物體邊緣增加，變成中間那張圖 configuration space，只剩下一個點和數個多邊形障礙物。接著可以轉換成 visibility graph，把剩下的白色空間進行梯形剖分或者三角剖分，將區域表示成一個點，相鄰區域拉一條邊，就成 dual graph，就能套用一般的最短路徑算法，來協助機器人行走。

處理這問題對於某 M 過於困難，於是將問題更簡化些，只需要把中間那一張圖其中一個障礙物變化求出即可。

現在給定一個凸多邊形障礙物以及移動凸多邊形，假設定位點在深黑色點方形部分，藉由貼著障礙物可以構出新的障礙物，請問新的障礙物為何？

Input

第一行會有一個整數 $T$，表示有多少組測資。

每一組測資會包含兩個凸多邊形。第一行有一個整數 $N$ 表示移動凸多邊形的頂點數量，接下來會有 $N$ 行，每一行會有兩個整數 $x, y$ 表示頂點座標。在 $N+1$ 行之後，會有一個整數 $M$ 表示障礙物凸多邊形的頂點數量，接下來會有 $M$ 行，每一行會有兩個整數 $x, y$ 表示頂點座標。

$2 < N, M < 512$
$-32767 < x, y < 32767$
假設定位點 (深黑色點方形部分) 都設在原點 (0, 0)。
輸入順序會按照順時針或逆時針給定。

Sample Input

Sample Output

1 2	Case #1: 89.0 Case #2: 115.5

Hint

Solution

可以用 $O(nm \log nm)$ 來完成此題，方法很簡單，先將機器人的圖形按照原點旋轉 180 度，接著將障礙物每一個頂點座標加上機器人的座標，總共會得到 $nm$ 的頂點，套用凸包算法求一次即可。

這一題應用 Minkowski Sum 來完成，對於兩個凸多邊形的 Minkowski Sum，可以在 $O(n+m)$ 時間內完成，將兩個凸包採用逆時針儲存，挑選最左，等價時抓下方 y 座標的點，可以知道新凸包的頂點一定是由兩個凸包的頂點$A_i, B_j$ 疊加而成，一開始抓到的左下方頂點疊加可以得到第一個頂點，接著要按照凸包順序求出，則比較$(A_i, A_{i+1})$ 和$(B_j, B_{j+1})$ 哪一個偏的角度小，就挑選哪一個往前推動 i = i+1 或者是 j = j+1。最後就能在 $O(n+m)$ 時間內完成。

假設不是凸多邊形，兩個簡單多邊形最慘會到 $O(n^2 m^2)$。

Read More +

2015-06-05

解題區/出題解題

[ACM 題目] 優勢商品

Problem

一份產品有很多屬性，把每一個屬性當作一個維度，當成 $D$ 個維度的向量，當維度越多時絕對支配 (各方面都是最好的) 的機率就越低。若把 $D = 2$，就是在二維空間中的支配點 (Zerojudge d555 平面上的極大點)。由於絕對優勢的產品隨著維度增加而變少，定義出一種相對優勢產品，給定一個 $K$，若 $p_i$ 支配 (k-dominate) $p_j$，必須滿足下列三條：

$\exists S' \subseteq S, |S'| = K$
$\forall s_r \in S', p_i.s_r \geq p_j.s_r$
$\exists s_t \in S', p_i.s_t > p_j.s_t$

用中文解釋這幾條數學公式，假設有 6 個維度的產品，目標 $K = 5$，則表示要在 6 的維度中任意挑選 5 個維度比較，若存在一種挑法可以支配，則可以說 $p_i$ k-dominate $p_j$

例如現在有 8 台筆記型電腦，共有 6 種屬性可以比較，數值越大越好。

Notebook	CPU	RAM	HDD	HDD S.	Q.	VRAM
$N_1$	3	3	5	6	6	8
$N_2$	9	4	9	7	7	9
$N_3$	8	4	7	9	2	7
$N_4$	5	6	8	9	5	9
$N_5$	9	7	9	6	2	4
$N_6$	6	6	6	5	3	5
$N_7$	5	7	3	8	4	6
$N_8$	4	4	8	6	6	4

若要判斷 $N_2$ 是否支配 $N_3$，從中挑取 (CPU, RAM, HDD, HDD S., Q.) 這 5 個屬性比較，得到 $N_2$ 勝過於 $N_3$。最後可以發現到 $N_2$ 這產品可以支配 $N_1, N_3, N_5, N_6, N_8$。

現在問題來了，要找到不被任何產品支配的產品 (k-dominant skyline)。如上表 8 項產品的情況，只會有 $N_2, N_4$。

Input

第一行會有一個整數 $T$ 表示有多少組測資。

每一組測資第一行會有三個整數 $N, D, K$，分別表示產品數量、產品屬性種類、以及支配比較的屬性種類。接下來會有 $N$ 行數據，每一行上會有 $D$ 個整數 $s_j$ 表示一個產品的屬性，產品按照輸入順序編號。

$N < 1024$
$0 < D, K < 8$
$0 < s_j < 1000$

Output

對於每組測資輸出一行，輸出該組測資所有的 k-dominant skyline 的產品編號，由小到大輸出。若不存在則輸出 NONE。

Sample Input

4
2 3 1
20 20 20
20 20 20
2 5 1
20 5 20 20 20
5 20 20 20 20
4 6 5
9 4 9 7 7 9
9 7 9 6 5 4
9 6 9 8 3 3
1 2 3 9 2 2
8 6 5
3 3 5 6 6 8
9 4 9 7 7 9
8 4 7 9 2 7
5 6 8 9 5 9
9 7 9 6 2 4
6 6 6 5 3 5
5 7 3 8 4 6
4 4 8 6 6 4

Sample Output

Case #1: 1 2
Case #2: NONE
Case #3: 1
Case #4: 2 4

Solution

這是一篇論文題目 k-dominant skyline set，網路上提出不少算法，看起來都是有容錯率。

這一題只是介紹題，沒有強求高效算法計算，直接 $O(N^2)$ 的比較，稍微優化一下就可以。若要快一點，按照總和由大排到小，從概念上總和越大，表示支配別人的機率越高，那麼根據這種貪心思路去篩選，可以減少很多不必要的比較，速度會更快些。

#include <stdio.h>
#include <stdlib.h>
#include <set>
#include <queue>
#include <string.h>
#include <algorithm>
using namespace std;
const int MAXN = 1024;
int N, D, K;
vector< vector<int> > C;
struct Product {
    int v[8], id;
    bool operator<(const Product &a) const {
        return sum() > a.sum();
    }
    int sum() const {
        int s = 0;
        for (int i = 0; i < D; i++)
            s += v[i];
        return s;
    }
    int domain(Product &u) {
        for (auto &x : C) {
            int h1 = 1, h2 = 0;
            for (auto &e : x) {
                h1 &= v[e] >= u.v[e];
                h2 |= v[e] >  u.v[e];
            }
            if (h1 && h2)   return 1;
        }
        return 0;
    }
} item[MAXN];
int used[MAXN];
int main() {
    int testcase, cases = 0;
    scanf("%d", &testcase);
    while (testcase--) {
        scanf("%d %d %d", &N, &D, &K);
        for (int i = 0; i < N; i++) {
            for (int j = 0; j < D; j++)
                scanf("%d", &item[i].v[j]);
            item[i].id = i;
        }
        sort(item, item+N);
        
        C.clear();
        for (int i = 0; i < (1<<D); i++) {
            if (__builtin_popcount(i) == K) {
                vector<int> A;
                for (int j = 0; j < D; j++) {
                    if ((i>>j)&1)
                        A.push_back(j);
                }
                C.push_back(A);
            }
        }
        
        vector<Product> filter;
        vector<int> ret;
        for (int i = 0; i < N; i++)
            used[i] = 0;
        for (int i = 0; i < N; i++) {
            if (used[i] == 0) {
                filter.push_back(item[i]);
                for (int j = i+1; j < N; j++) {
                    if (used[j] == 0 && item[i].domain(item[j]))
                        used[j] = 1;
                }
            }
        }
        for (int i = 0; i < filter.size(); i++) {
            int ok = 1;
            for (int j = 0; j < N && ok; j++) {
                if (item[j].domain(filter[i]))
                    ok = 0;
            }
            if (ok == 1)
                ret.push_back(filter[i].id);
        }
        
        sort(ret.begin(), ret.end());
        printf("Case #%d:", ++cases);
        for (int i = 0; i < ret.size(); i++)
            printf(" %d", ret[i]+1);
        puts(ret.size() ? "" : " NONE");
    }
    return 0;
}

Read More +

2015-06-02

解題區/解題區 - 其他題目

2015 Google Code Jam Round 2

超過凌晨的競賽，隔天起床有一種說不出的疲累感，坑了一個早上，仍然沒想到最後幾題。看來已經玩不下去！圍觀到這裡。

A. Pegman

在谷歌地圖上放置一個小人，小人會根據當前給予的方向指標持續地移動，直到碰到下一個方向指標，請問至少要改變幾個方向指標所指引的方向，滿足在任何一個小人位置都不會走到邊界外部。

由於放在非指標位置就不能移動，只需要考慮小人放在指標上，考慮每一個方向指標都指向另一個指標，這樣就能保證有數個循環，盡可能地挑選原方向。

#include <stdio.h>
#include <string.h>
#include <queue>
#include <functional>
#include <deque>
#include <algorithm>
using namespace std;
char sg[128][128];
int ig[128][128];
const int dx[] = {0, 0, 1, -1};
const int dy[] = {1, -1, 0, 0};
const char dd[] = "><v^";
int main() {
    int testcase, n, m, cases = 0;
    scanf("%d", &testcase);
    while (testcase--) {
        scanf("%d %d", &n, &m);
        for (int i = 0; i < n; i++)
            scanf("%s", sg[i]);
        
        int idx = 0;
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < m; j++) {
                if (sg[i][j] != '.') {
                    ig[i][j] = idx++;
                }
            }
        }
        
        int match = 0, ret = 0;
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < m; j++) {
                if (sg[i][j] != '.') {
                    int odd = 0;
                    for (int k = 0; k < 4; k++)
                        if (dd[k] == sg[i][j])
                            odd = k;
                    int cost = 1, mm = 0;
                    for (int k = 0; k < 4; k++) {
                        int x = i, y = j;
                        x += dx[k], y += dy[k];
                        while (x < n && x >= 0 && y < m && y >= 0) {
                            if (sg[x][y] != '.') {
                                if (k == odd)
                                    cost = 0;
                                mm = 1;
                                break;
                            }
                            x += dx[k], y += dy[k];
                        }
                    }
                    match += mm, ret += cost;
                }
            }
        }
        
        printf("Case #%d: ", ++cases);
        if (match != idx)
            puts("IMPOSSIBLE");
        else
            printf("%d\n", ret);
    }
    return 0;
}

B. Kiddie Pool

放滿游泳池的水，目標容積 $V$ 溫度 $X$，有 $N$ 種水源，每一種水源有各自的流速 $R$、溫度 $C$，水源可以同時放水，請問達到目標容積、溫度的最少時間為何。

當下直觀地沒考慮水源可以同時使用，卡了一陣子。

Linear Programming

發現可以同時運作，考慮二分最少時間，其中一種最簡單的應用線性規劃，判斷線性規劃方程式是否有解，python 內建 LP 也許很過分。自己曾經寫過 simplex，不過調校上相當痛苦，判斷有沒有解有困難。以下為 LP 模型，二分時間 $\text{limit_time}$。

$x_i \le \text{limit_time} * R[i]$
$\sum{C_i x_i} = V X$
$\sum x_i = V$
$x_i \ge 0$

帶入 simplex algorithm，判斷四條方程式是否有解。代碼就不附，解不出來。

Greedy

另外一種方式，考慮限制 t 時間內能調配的最高溫 high_t、最低溫 low_t，這兩種溫度可以貪心法求得，接著考慮目標溫度 low_t <= X <= high_t，因為中間的溫度一定可解，具有連續性。

複雜度 $O(N \log X)$

#include <bits/stdc++.h> 
using namespace std;
const int MAXN = 128;
const double eps = 1e-10;
double R[MAXN], C[MAXN], V, X;
int N;
vector< pair<double, double> > A;
int checkTime(double limitT) {
    double v, t, low_t, high_t;
    double mxV, teC;
    
    v = t = 0;
    for (int i = 0; i < N; i++) {
        mxV = A[i].second * limitT; // flow in limit time
        teC = A[i].first;
        
        mxV = min(mxV, V - v);
        t = (t * v + mxV * teC) / (v + mxV);
        v += mxV;
        if (v >= V - eps)
            break;
    }
    if (v < V - eps)
        return 0;
    low_t = t;
    
    v = t = 0;
    for (int i = N-1; i >= 0; i--) {
        mxV = A[i].second * limitT; // flow in limit time
        teC = A[i].first;
        
        mxV = min(mxV, V - v);
        t = (t * v + mxV * teC) / (v + mxV);
        v += mxV;
        if (v >= V - eps)
            break;
    }
    if (v < V - eps)
        return 0;
    high_t = t;
    return X >= low_t - eps && X <= high_t + eps;
}
int main() {
    int testcase, cases = 0;
    scanf("%d", &testcase);
    while (testcase--) {
        scanf("%d %lf %lf", &N, &V, &X);
        
        A.clear();
        for (int i = 0; i < N; i++) {
            scanf("%lf %lf", &R[i], &C[i]);
            A.push_back(make_pair(C[i], R[i]));
        }
        sort(A.begin(), A.end());
        
        double mxT, mnT;
        mnT = A[0].first, mxT = A[N-1].first;
        
        printf("Case #%d: ", ++cases);
        if (X < mnT - eps || X > mxT + eps) {
            puts("IMPOSSIBLE");
        } else {
            double l = 0, r = 1e+9, mid; // MAXV / MINR = 10000.0 / 0.0001 
            double ret = 0;
            for (int it = 0; it < 256; it++) {
                mid = (l + r)/2;
                if (checkTime(mid))
                    r = mid, ret = mid;
                else
                    l = mid;
            }
            printf("%.8lf\n", ret);
        }
    }
    return 0;
}

Minkowski Sum

在此先將問題轉換成「求 1 秒內最多能湊出多少公升的 X 度水」，知道單位時間內的最多流量，就能貪心得到最短時間。

將每一種水源 $(R_i, C_i)$ 轉換成向量 $v_i = (R_i, R_i \times C_i)$，目標要在 $t$ 時間內，得到 $\sum{t_i v_i} = (V, V \times X)$，滿足所有的 $t_i \le t$

將問題壓縮成 $t_i \in \left \[ 0, 1 \right]$，將所有符合的 $\sum{t_i v_i}$ 疊加起來，得到的區域相當於在做 Minkowski Sum，區域是一個凸多邊形。接著在 $y = X x$ 尋找交集的最大 x 值。

Demo

上圖是一個 3 個水源的情況，假設要湊出溫度 $X = 0.5$，相當於找 $y = 0.5 x$，這三個向量在單位時間內的疊加為淡紅色區域。為了得到這淡紅色區域，使用極角排序，依序疊加可以得到下凸包，反過來疊加可以得到上凸包。

明顯地要找一個凸包跟一條線的交點，得到單位時間內的最大流量。此時的 x 軸表示流量、y 軸表示熱量，找到交集中最大的 x 座標值。

#include <stdio.h>
#include <math.h>
#include <algorithm>
#include <set>
#include <vector>
using namespace std;
#define eps 1e-6
struct Pt {
    double x, y;
    int label;
    Pt(double a = 0, double b = 0, int c = 0):
    x(a), y(b), label(c) {}
    Pt operator-(const Pt &a) const {
        return Pt(x - a.x, y - a.y);
    }
    Pt operator+(const Pt &a) const {
        return Pt(x + a.x, y + a.y);
    }
    Pt operator*(const double a) const {
        return Pt(x * a, y * a);
    }
    bool operator<(const Pt &a) const {
        if (fabs(x - a.x) > eps)
            return x < a.x;
        if (fabs(y - a.y) > eps)
            return y < a.y;
        return false;
    }
};
double dot(Pt a, Pt b) {
    return a.x * b.x + a.y * b.y;
}
double cross(Pt o, Pt a, Pt b) {
    return (a.x-o.x)*(b.y-o.y)-(a.y-o.y)*(b.x-o.x);
}
double cross2(Pt a, Pt b) {
    return a.x * b.y - a.y * b.x;
}
int between(Pt a, Pt b, Pt c) {
    return dot(c - a, b - a) >= -eps && dot(c - b, a - b) >= -eps;
}
int onSeg(Pt a, Pt b, Pt c) {
    return between(a, b, c) && fabs(cross(a, b, c)) < eps;
}
bool cmp(const Pt& p1, const Pt& p2)
{
    if (p1.y == 0 && p2.y == 0 && p1.x * p2.x <= 0) return p1.x > p2.x;
    if (p1.y == 0 && p1.x >= 0 && p2.y != 0) return true;
    if (p2.y == 0 && p2.x >= 0 && p1.y != 0) return false;
    if (p1.y * p2.y < 0) return p1.y > p2.y;
    double c = cross2(p1, p2);
    return c > 0 || (c == 0 && fabs(p1.x) < fabs(p2.x));
}
Pt getIntersect(Pt as, Pt ae, Pt bs, Pt be) {
    Pt u = as - bs;
    double t = cross2(be - bs, u)/cross2(ae - as, be - bs);
    return as + (ae - as) * t;
}
int cmpZero(double v) {
    if (fabs(v) > eps) return v > 0 ? 1 : -1;
    return 0;
}
//
const int MAXN = 128;
int N;
double R[MAXN], C[MAXN], V, X;
void solveWithMinkowskiSum() {
    vector<Pt> P;
    for (int i = 0; i < N; i++)
        P.push_back(Pt(R[i], R[i]*C[i]));
    sort(P.begin(), P.end(), cmp);  // solar sort
    
    vector<Pt> convex;
    double mxFlow = 0;              // in unit time with temperature X
    Pt u, v, o(0, 0), to(1, X);     // y = (X) x, maximum x
    
    u = v = Pt(0, 0);
    convex.push_back(u);
    for (int i = 0; i < N; i++) {
        v = u + P[i];
        u = v;
        convex.push_back(v);
    }
    
    reverse(convex.begin(), convex.end());
    u = v = Pt(0, 0);
    for (int i = N-1; i >= 0; i--) {
        v = u + P[i];
        u = v;
        convex.push_back(v);
    }
    
    for (int i = 0, j = (int) convex.size()-1; i < convex.size(); j = i++) {
        u = convex[j], v = convex[i];
        if (cmpZero(cross(o, to, u)) * cmpZero(cross(o, to, v)) < 0) {
            Pt in = getIntersect(o, to, u, v);
            mxFlow = max(mxFlow, in.x);
        }
        if (cmpZero(cross(o, to, v)) == 0)
            mxFlow = max(mxFlow, v.x);
    }
    
    if (fabs(V) < eps)
        printf("%.10lf\n", 0.0);
    else if (fabs(mxFlow) < eps)
        puts("IMPOSSIBLE");
    else
        printf("%.10lf\n", V / mxFlow);
}
int main() {
    int testcase, cases = 0;
    scanf("%d", &testcase);
    while (testcase--) {
        scanf("%d %lf %lf", &N, &V, &X);
        
        V *= 10000, X *= 10000;
        for (int i = 0; i < N; i++) {
            scanf("%lf %lf", &R[i], &C[i]);
            R[i] *= 10000, C[i] *= 10000;
        }
        
        printf("Case #%d: ", ++cases);
        solveWithMinkowskiSum();
    }
    return 0;
}