解題區/解題區 - Zerojudge

b486. 變態史考古

Problem

題目描述

隨著考古學者近幾年的研究，逐漸地分析出變態人類文化族群的祖先。同時也知道，每一代人會將其變態基因傳遞給下一代，而下一代的變態指數會比上一代多一。若子代有數個，每個子代帶有的變態性向會有所歧異。

每當新聞報導某個變態父子關係時，人們總會熱議某兩個變態的相似程度。口耳相傳的消息是不可信任的，現在給予事情發生的時間軸，請驗證人們口中的變態相似程度。

變態相似程度定義

$\textit{sim}(A, B) = \frac{\textit{dist}(\textit{Original}(A), \textit{LCA}(A, B)) \times 2}{\textit{dist}(\textit{Original}(A), A) + \textit{dist}(\textit{Original}(B), B)}$

其中 $\textit{Original}(A)$ 表示根據事實推測得到的最早源頭祖先，$\textit{dist}(A, B)$ 表示在樹狀圖中兩點的距離，$\textit{LCA}(A, B)$ 表示 $A, \; B$ 最小公共祖先 (Lowest Common Ancestor，簡稱 LCA)。

例如當前知道的祖先關係如下

1
2
3
4
5
6
7
    P
   / 
  Q
 / \
A   R
   /
  B
  • $\textit{Original}(A) = \textit{Original}(B) = P$，
  • $\textit{dist}(P, A) = 3, \; \textit{dist}(P, B) = 4, \; \textit{LCA}(A, B) = Q, \; \textit{dist}(Q, P) = 2$。
  • 最後 $\textit{sim}(A, B) = \frac{4}{7}$

Sample Input

1
2
3
4
5
6
7
8
9
10
5 9
news 3 2
news 4 2
sim 3 4
sim 3 5
news 2 1
news 5 4
sim 3 4
sim 3 5
sim 4 4

Sample Output

1
2
3
4
5
1/2
-1
2/3
4/7
1/1

Solution

維護一個有根森林的建造，在建造過程中找尋 LCA。出這一題時，測資只有打算卡暴力解，意即使用慢慢往上爬找交會的算法。預計會有很多種解法，因為建樹結果是唯一的，於是柳州 101 果真給我一個離線的方案。

使用離線處理把建樹的資訊保留，利用倍增算法找攀升的祖先，隨後模擬時，使用並查集維護當前所在樹的樹根，由於要找到 LCA 到當前樹根的距離，其實就相當於找到最終樹根的距離去扣除即可。

另外的解法是採用 Link/Cut Tree，直接套模板找 LCA 就相當簡單。

離線處理 + 倍增算法

倍增算法實作時特別小心快取置換，例如陣列宣告 arr[MAXLOGN][MAXN] 和 arr[MAXN][MAXLOGN] 使用 dp 建表時，前者速度會比後者快上三倍之多 (當 MAXN = 200000 以上)。

但是在找 LCA 時，access memory 又會變慢，所以見仁見智。若是在 RMQ 的 sparse table 處理上會有很大的影響。

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
#include <bits/stdc++.h> 
using namespace std;
const int MAXN = 100005;
const int MAXLOGN = 17;
struct Event {
    char s[8];
    int x, y;
    void read() {
        scanf("%s %d %d", s, &x, &y);
    }
};
int fa[MAXLOGN][MAXN], visited[MAXN], dep[MAXN];
int parent[MAXN];
vector<int> g[MAXN];
void dfs(int u) {
    visited[u] = 1;
    for (auto v : g[u]) {
        fa[0][v] = u;
        dep[v] = dep[u] + 1;
        dfs(v);
    }
}
int findp(int x) {
    return x == parent[x] ? x : (parent[x]=findp(parent[x]));
}
int LCA(int x, int y) {
    if (dep[x] < dep[y])	swap(x, y);
    int d = dep[x] - dep[y];
    for (int i = MAXLOGN-1; i >= 0; i--) {
        if ((d>>i)&1) {
            d -= 1<<i;
            x = fa[i][x];
        }
    }
    if (x == y)	return x;
    for (int i = MAXLOGN-1; i >= 0; i--) {
        if (fa[i][x] != fa[i][y]) {
            x = fa[i][x], y = fa[i][y];
        }
    }
    return fa[0][x];
}
void sim(int x, int y) {
    int rx, ry, lca;
    rx = findp(x), ry = findp(y);
    if (rx != ry) {
        puts("-1");
        return ;
    }
    lca = LCA(x, y);
    int a = (dep[lca]-dep[rx]+1)*2, b = dep[x]+dep[y]-dep[rx]*2+2, g;
    g = __gcd(a, b), a /= g, b /= g;
    printf("%d/%d\n", a, b);
}
int main() {
    int n, m;
    while (scanf("%d %d", &n, &m) == 2) {
        for (int i = 1; i <= n; i++)
            g[i].clear(), visited[i] = 0, parent[i] = i;
        vector<Event> e(m);
        for (int i = 0; i < m; i++) {
            e[i].read();
            if (e[i].s[0] == 'n')
                g[e[i].y].push_back(e[i].x);
        }
        
        memset(visited, 0, (n+1) * sizeof(visited[0]));
        for (int i = 1; i <= n; i++) {
            if (visited[i] == 0)
                dfs(i);
        }
        
        for (int i = 1; i < MAXLOGN; i++)
            for (int j = 1; j <= n; j++)
                fa[i][j] = fa[i-1][fa[i-1][j]];
                
        for (int i = 0; i < m; i++) {
            if (e[i].s[0] == 'n') {
                parent[e[i].x] = e[i].y;
            } else {
                sim(e[i].x, e[i].y);
            }
        }
    }
    return 0;
}

LCT

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
#include <bits/stdc++.h> 
using namespace std;
class LCT { // Link-Cut Tree
public:
    static const int MAXN = 100005;
    struct Node {
        static Node *EMPTY;
        Node *ch[2], *fa;
        int rev;
        int size;
        Node() {
            ch[0] = ch[1] = fa = NULL;
            rev = 0; 
            size = 1;
        }
        bool is_root() {
            return fa->ch[0] != this && fa->ch[1] != this;
        }
        void pushdown() {
            if (rev) {
                ch[0]->rev ^= 1, ch[1]->rev ^= 1;
                swap(ch[0], ch[1]);
                rev ^= 1;
            }
        }
        void pushup() {
            if (this == EMPTY)	return ;
            size = 1;
            if (ch[0] != EMPTY)
                size += ch[0]->size;
            if (ch[1] != EMPTY)
                size += ch[1]->size;
        }
    } _mem[MAXN];
    
    int bufIdx;
    LCT() {
        Node::EMPTY = &_mem[0];
        Node::EMPTY->fa = Node::EMPTY->ch[0] = Node::EMPTY->ch[1] = Node::EMPTY;
        Node::EMPTY->size = 0;
        bufIdx = 1; 
    }
    void init() {
        bufIdx = 1;
    }
    Node* newNode() {
        Node *u = &_mem[bufIdx++];
        *u = Node();
        u->fa = u->ch[0] = u->ch[1] = Node::EMPTY;
        return u;
    }
    void rotate(Node *x) {
        Node *y;
        int d;
        y = x->fa, d = y->ch[1] == x ? 1 : 0;
        x->ch[d^1]->fa = y, y->ch[d] = x->ch[d^1];
        x->ch[d^1] = y;
        if (!y->is_root())
            y->fa->ch[y->fa->ch[1] == y] = x;
        x->fa = y->fa, y->fa = x;
        y->pushup(), x->pushup();
    }
    void deal(Node *x) {
        if (!x->is_root())	deal(x->fa);
        x->pushdown();
    }
    void splay(Node *x) {
        Node *y, *z;
        deal(x);
        while (!x->is_root()) {
            y = x->fa, z = y->fa;
            if (!y->is_root()) {
                if (y->ch[0] == x ^ z->ch[0] == y)
                    rotate(x);
                else
                    rotate(y);
            }
            rotate(x);
        }
        x->pushup();
    }
    Node* access(Node *u) {
        Node *v = Node::EMPTY;
        for (; u != Node::EMPTY; u = u->fa) {
            splay(u);
            u->ch[1] = v;
            v = u;
            v->pushup();
        }
        return v;
    }
    void mk_root(Node *u) {
        access(u)->rev ^= 1, splay(u);
    }
    void cut(Node *x, Node *y) {
        mk_root(x);
        access(y), splay(y);
        y->ch[0] = x->fa = Node::EMPTY;
    }
    Node* cut(Node *x) {
        Node *u, *v, *rt = find(x);
        mk_root(rt);
        access(x), splay(x);
        for (v = x->ch[0]; v->ch[1] != Node::EMPTY; v = v->ch[1]);
        x->ch[0]->fa = x->fa;
        x->fa = x->ch[0] = Node::EMPTY;
        return v;
    }
    void link(Node *x, Node *y) {
        mk_root(x);
        x->fa = y;
    }
    Node* find(Node *x) {
        for (x = access(x); x->pushdown(), x->ch[0] != Node::EMPTY; x = x->ch[0]);
        return x;
    }
    //
    Node* lca(Node *x, Node *y) {
        if (x == y)	return x;
        if (find(x) != find(y))
            return Node::EMPTY;
        Node *u, *v = Node::EMPTY, *LCA = Node::EMPTY;
        access(y), splay(y);
        for (u = x; u != Node::EMPTY; u = u->fa) {
            splay(u);
            if (u->fa == Node::EMPTY) {
                LCA = u;
//				u->ch[1]->push_add(c), v->push_add(c);
            }
            u->ch[1] = v;
            v = u;
            v->pushup();
        }
        return LCA;
    }
    int path(Node *x, Node *y) {
        int ret;
        Node *u, *v = Node::EMPTY;
        access(y), splay(y);
        for (u = x; u != Node::EMPTY; u = u->fa) {
            splay(u), u->pushdown();
            if (u->fa == Node::EMPTY) {
                ret = u->ch[1]->size + v->size;
            }
            u->ch[1] = v;
            v = u;
            v->pushup();
        }
        return ret;
    }
    inline int label(Node *x) {
        return x - _mem;
    }
} tree;
LCT::Node *LCT::Node::EMPTY;
LCT::Node *A[100005], *node_x, *node_rt;
int main() {
    int n, m, x, y, c, u, v;
    char cmd[8];
    while (scanf("%d %d", &n, &m) == 2) {
        tree.init();
        for (int i = 1; i <= n; i++)
            A[i] = tree.newNode();
        
        for (int i = 0; i < m; i++) {
            scanf("%s", cmd);
            if (cmd[0] == 'n') {		// link
                scanf("%d %d", &x, &y);
                tree.link(A[x], A[y]);
            } else if (cmd[0] == 's') {	// sim
                scanf("%d %d", &x, &y);
                node_x = tree.lca(A[x], A[y]);
                int lca = tree.label(node_x);
                if (lca == 0) { 
                    puts("-1");
                } else {
                    node_rt = tree.find(A[lca]);
                    int p1, p2;
                    p1 = tree.path(A[x], A[y]) + 1;
                    p2 = tree.path(A[lca], node_rt) + 1;
                    int a = p2*2, b = p1 + p2*2 - 1, g;
                    g = __gcd(a, b), a /= g, b /= g;
                    printf("%d/%d\n", a, b);
                } 
            }
        }
    }
    return 0;
}
Read More +
解題區/解題區 - Zerojudge

b483. 史蒂芙的觀察日記

Problem

背景

史蒂芙在學院讀書時，教授曾經教過最小生成樹 (Minimum Spanning Tree，簡稱 MST)，MST 最常見到的兩個算法 Kruskal’s Algorithm、Prim’s Algorithm，還有其他的算法來得到 MST。

在某一次作業中，史蒂芙被特別指派的某一道題目，題目描述「給定好一棵 MST，接著新增加一條邊，請問新的最小生成樹成本為何。」聰明的史蒂芙想到一個樸素的解法「增加新的一條邊若造成 MST 上存在一個環，只要把環上最權重最大的邊移除，最小生成樹就大功告成。」

萬萬沒想到，正高興地要拿著 $O(V)$ 算法交差時，卻沒想到一連串的詢問，於是史蒂芙被命令回去觀察觀察一番再交出報告。「看來史蒂芙還是只有被欺負的份呢。」

題目描述

給定 $N$ 個點、$M$ 條雙向邊，依序給定每一條邊，每增加一條邊，輸出當前的最小生成樹成本，若圖不連通，則輸出最小生成森林成本。

Sample Input

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

Sample Output

1
2
3
4
5
6
7
8
9
10
11
12
Case #1
1
3
3
Case #2
9
17
24
30
27
22
19

Solution

用 Link/Cut tree 解 MST 相當瘋狂的舉動，可以將找環上最重邊移除的操作降到 $O(\log V)$。最小生成樹是一個無向的樹，而邊的權重可以轉化成一個點，例如 ei(x, y, v) 相當於 x <-> ei <-> y 其中節點 x, y 帶有的權重為 0，而 ei 是一個虛設點帶權重為 v。

假設增加一條邊 ej(x, y, v) 會產生一個環，先檢查節點 x, y 是否在同一個集合內，若不在同一個集合直接建起這一條邊。反之，找到路徑 path(x, y) 上的點權最大值，接著映射回去，移除掉虛設點的兩個連接邊 x <-> ei 和 ei <-> y。

由於是一個無向樹，在尋找路徑時可以使用以下寫法 (否則要用另外一種，但這種寫法比較簡潔)

1
2
3
4
5
Node* maxPath(Node *x, Node *y) {
    mk_root(x);
    access(y), splay(y);
    return y->mxv;
}

如果 x 本身不是根，則這個操作會改變樹的上下層關係，對於有向樹會有差別，如 Node* find(Node *x) 中反映出來並非有向樹的唯一根，只能用來計算並查集的關係 find(x) == find(y)。

  • 能不能減少虛設點的使用？
    目前我沒辦法解決這問題，若要減少虛設點，則每個節點的權重代表維護到父節點的邊權重 (可以參考 UVa 11994 - Happy Painting)，仍然可以找到路徑上的最重邊，問題在於加入新的一條邊時，兩個有根樹要合併，由於其中一端可能不是其中一個的樹根，其上下層關係需要反轉，則邊權代表邊權的能力就會失效。也許可以利用標記傳遞來重新配置，寫起來恐怕會非常複雜。而虛設點會讓點數增加兩倍，耗時多兩倍是可以接受的方案。

小記

最近在撸 Link/Cut Tree (LCT)，於是拿去解決最小生成樹 MST，複雜度仍然是 $O(n \log 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
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
#include <bits/stdc++.h> 
using namespace std;
class LCT { // Link-Cut Tree
public:
    static const int MAXN = 262144;
    struct Node {
        static Node *EMPTY;
        Node *ch[2], *fa;
        int rev;
        int val, size;
        Node *mxv;
        Node() {
            ch[0] = ch[1] = fa = NULL;
            rev = 0;
            val = 0;
            mxv = this, size = 1;
        }
        bool is_root() {
            return fa->ch[0] != this && fa->ch[1] != this;
        }
        void pushdown() {
            if (rev) {
                ch[0]->rev ^= 1, ch[1]->rev ^= 1;
                swap(ch[0], ch[1]);
                rev ^= 1;
            }
        }
        void pushup() {
            if (this == EMPTY)	return ;
            mxv = this, size = 1;
            if (ch[0] != EMPTY) {
                if (ch[0]->mxv->val > mxv->val)	mxv = ch[0]->mxv;
                size += ch[0]->size;
            }
            if (ch[1] != EMPTY) { 
                if (ch[1]->mxv->val > mxv->val)	mxv = ch[1]->mxv;
                size += ch[1]->size;
            } 
        }
    } _mem[MAXN];
    
    int bufIdx;
    LCT() {
        Node::EMPTY = &_mem[0];
        Node::EMPTY->fa = Node::EMPTY->ch[0] = Node::EMPTY->ch[1] = Node::EMPTY;
        Node::EMPTY->size = 0;
        bufIdx = 1; 
    }
    void init() {
        bufIdx = 1;
    }
    Node* newNode() {
        Node *u = &_mem[bufIdx++];
        *u = Node();
        u->fa = u->ch[0] = u->ch[1] = Node::EMPTY;
        return u;
    }
    void rotate(Node *x) {
        Node *y;
        int d;
        y = x->fa, d = y->ch[1] == x ? 1 : 0;
        x->ch[d^1]->fa = y, y->ch[d] = x->ch[d^1];
        x->ch[d^1] = y;
        if (!y->is_root())
            y->fa->ch[y->fa->ch[1] == y] = x;
        x->fa = y->fa, y->fa = x;
        y->pushup(), x->pushup();
    }
    void deal(Node *x) {
        if (!x->is_root())	deal(x->fa);
        x->pushdown();
    }
    void splay(Node *x) {
        Node *y, *z;
        deal(x);
        while (!x->is_root()) {
            y = x->fa, z = y->fa;
            if (!y->is_root()) {
                if (y->ch[0] == x ^ z->ch[0] == y)
                    rotate(x);
                else
                    rotate(y);
            }
            rotate(x);
        }
        x->pushup();
    }
    Node* access(Node *u) {
        Node *v = Node::EMPTY;
        for (; u != Node::EMPTY; u = u->fa) {
            splay(u);
            u->ch[1] = v;
            v = u;
            v->pushup();
        }
        return v;
    }
    void mk_root(Node *u) {
        access(u)->rev ^= 1, splay(u);
    }
    void cut(Node *x, Node *y) {
        mk_root(x);
        access(y), splay(y);
        y->ch[0] = x->fa = Node::EMPTY;
    }
    Node* _cut(Node *rt, Node *x) {
        Node *u, *v;
        mk_root(rt);
        access(x), splay(x);
        for (v = x->ch[0]; v->ch[1] != Node::EMPTY; v = v->ch[1]);
        x->ch[0]->fa = x->fa;
        x->fa = x->ch[0] = Node::EMPTY;
        return v;
    }
    void link(Node *x, Node *y) {
        mk_root(x);
        x->fa = y;
    }
    Node* find(Node *x) {
        for (x = access(x); x->pushdown(), x->ch[0] != Node::EMPTY; x = x->ch[0]);
        return x;
    }
    //
    Node* maxPath(Node *x, Node *y) {
        mk_root(x);
        access(y), splay(y);
        return y->mxv;
    }
    void debug(int n) {
    }
} tree;
LCT::Node *LCT::Node::EMPTY;
LCT::Node *A[262144];
int x[131072], y[131072];
int main() {
    int n, m, w, u, v, cmd;
    int cases = 0;
    while (scanf("%d %d", &n, &m) == 2) {
        printf("Case #%d\n", ++cases);
        tree.init();
        for (int i = 1; i <= n; i++)
            A[i] = tree.newNode();
        
        int vcnt = n+1, mst = 0;
        for (int i = 1; i <= m; i++) {
            scanf("%d %d %d", &x[i], &y[i], &w);
            A[i+n] = tree.newNode();
            A[i+n]->val = w;
            if (tree.find(A[x[i]]) != tree.find(A[y[i]])) {
                tree.link(A[x[i]], A[i+n]);
                tree.link(A[y[i]], A[i+n]);
                mst += w;
            } else {
                LCT::Node *e = tree.maxPath(A[x[i]], A[y[i]]);
                if (e->val > w) {
                    mst += w - e->val;
                    int eid = e - tree._mem - n;
//					printf("eeeee %d %d %d\n", x[eid], y[eid], e->val);
                    tree.cut(A[x[eid]], e);
                    tree.cut(A[y[eid]], e);
                    tree.link(A[x[i]], A[i+n]);
                    tree.link(A[y[i]], A[i+n]);
                }
            }
            printf("%d\n", mst);
        }
    }
    return 0;
}
/*
3 3
1 2 1
1 3 2
2 3 4
*/
Read More +
解題區/解題區 - Zerojudge

b478. 有限間距最長共同子序列

Problem

背景

有限間距最長共同子序列 (Gapped Longest Common Subsequence, 簡稱 GLCS) 是一種 LCS 的變形，而 LCS 想必早已成了經典題型。

回顧一下最長共同子序列 LCS 的解法，可以利用最簡單的 $O(NM)$ 動態規劃完成或者在特殊不重複情況下轉化成 LIS 在 $O(N \log N)$ 時間完成，轉換成 LIS 是有風險的，有機會退化成 $O(NM \log N)$。說不定還有比 $O(NM)$ 更快的算法去計算 LCS。

GLCS 的特點在於挑選子序列時，挑選到的相鄰字符位置之間最多有 $K$ 個字符不選。可強迫挑選出的相似序列盡可能是有意義的。

例如兩個字符串 $S, \; T$ 要求 GLCS，當 $K = 2$ 時：

$$S = \text{ACBDCAA} \\ T = \text{ADDBCDBAC}$$

則 ABCA 是一組合法的 $K = 2$ 時的 GLCS，挑選方案如下

1
2
3
4
5
    0123456789
S = ACBDCAA
    ^ ^ ^^
T = ADDBCDBAC
    ^  ^^  ^

若 $K = 1$ 時，由於字符串 $T$ 挑選 ABCA 時，前兩個字符之間的不選字符數為 2，不符合小於等於 $K$，則 ABCA 不能是 $K = 1$ 的合法 GLCS，而 CDA 是一組合法 $K = 1$ 時的 GLCS。

1
2
3
4
5
    0123456789
S = ACBDCAA
     ^ ^ ^
T = ADDBCDBAC
        ^^ ^

題目描述

給予兩個字符串 $S, \; T$ 以及限制間距 $K$，求出 GLCS 的長度為何。

Sample Input

1
2
3
2 ACBDCAA ADDBCDBAC
1 ACBDCAA ADDBCDBAC
0 ACBDCAA ADDBCDBAC

Sample Output

1
2
3
4
3
2

Solution

這一題是 Zerojudge a374. 5. 股票趨勢 的簡化，股票趨勢那題要求的 $K$ 會隨著不同字符變動。如果固定 $K$，則可以利用單調隊列降到 $O(N^2)$。

GLCS 的遞迴公式如下：

$$T[i, j] := \left\{\begin{matrix} \textit{undefined} & A[i] \neq B[j] \\ max \;(T[a, b])+1 & A[i] = B[j], i - K - 1 \le a < i, j - K - 1 \le b < j \end{matrix}\right.$$

可以轉換成二維的 RMQ 問題，套用樹套樹結構可以單筆詢問 $O(\log^2 n)$ 完成。由於 DP 的掃描線性質，若用單調隊列單一詢問可以在 $O(1)$ 完成。比較裸的二維單調隊列求正方形窗口極值可以參考 BZOJ 1047 [HAOI2007] 理想的正方形。

對於二維的滑動窗口，則可以對每一列維護一個 列 的單調隊列，對於某一行掃描時，再維護一個 行 的單調隊列，往右移動窗口時，把移動到 列 的單調隊列的極值加入，並且把超出窗口的列的極值移除。
由上至下、由左而右掃描每一行上的元素，找到以此元素為正方形右下角的正方形區域最大最小值。

時間複雜度從 $O(N^2 + R \log^2 N)$ 降到 $O(N^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
#include <bits/stdc++.h> 
using namespace std;
const int MAXN = 2048;
template<typename WTYPE, bool FUNC(WTYPE, WTYPE)>
class MonoQueue {
public:
    deque< pair<int, WTYPE> > Q;
    int K;
    void setWindow(int K) {
        this->K = K;
    }
    inline void add(int idx, WTYPE val) {
        while (!Q.empty() && Q.front().first <= idx - K)
            Q.pop_front();
        while (!Q.empty() && FUNC(val, Q.back().second))
            Q.pop_back();
        Q.push_back({idx, val});
    }
    inline pair<int, WTYPE> top() {
        return Q.front();
    }
};
bool cmp(int a, int b) { return a >= b; }
char s1[MAXN], s2[MAXN];
int dp[MAXN][MAXN];
int main() {
    int N, M, K, val;
    while (scanf("%d %s %s", &K, s1+1, s2+1) == 3) {
        N = strlen(s1+1), M = strlen(s2+1), K++;
        int GLCS = 0;
        MonoQueue<int, cmp> mxC[MAXN];
                
        for (int i = 0; i < M; i++)
            mxC[i].setWindow(K);
        for (int i = 0; i <= N; i++)
            for (int j = 0; j <= M; j++)
                dp[i][j] = 0;
                
        for (int i = 0; i <= N; i++) {
            MonoQueue<int, cmp> mxR;
            mxR.setWindow(K);
            for (int j = 0; j <= M; j++) {
                    
                val = dp[i][j];
                mxC[j].add(i, val);
                mxR.add(j, mxC[j].top().second);
                
                if (s1[i+1] == s2[j+1])
                    dp[i+1][j+1] = mxR.top().second+1;
                GLCS = max(GLCS, val);
            }
        }
        printf("%d\n", GLCS);
    }
    return 0;
}
Read More +
解題區/解題區 - Zerojudge

a374. 5. 股票趨勢

Problem

在限制若選用某個字母下，與前一個所選字母的間距不能大於某個特定值，請問 GLCS 長度為何？

Sample Input

1
2
3
4
5
6
7
8
ACBDCAA ADDBCDBAC
2
$
A 2 B 0 C 3 D 0 $
ACBDCAA ADDBCDBAC
1
C 4 A 6 $

Sample Output

1
2
4 3
4

Solution

明顯地，從一般的 LCS 變化，假設間距都是無窮大，則退化成一般的 LCS，那 LCS 也可以寫成以下方案

$$T[i, j] := \left\{\begin{matrix} \textit{undefined} & A[i] \neq B[j] \\ max \;(T[a, b])+1 & A[i] = B[j], a < i, b < j \end{matrix}\right.$$

從改寫過的 LCS 方案可以看出算法複雜度最慘會到 $O(N^4)$，事實上只有 $O(R N^2)$，其中 $R$ 表示字符匹配的數量。

由於要改寫成每一個字元都有其限制，則 GLCS 會變成以下遞迴公式

$$T[i, j] := \left\{\begin{matrix} \textit{undefined} & A[i] \neq B[j] \\ max \;(T[a, b])+1 & A[i] = B[j], i - K - 1 \le a < i, j - K - 1 \le b < j \end{matrix}\right.$$

可以標準的 RMQ 區間極值查找，但麻煩地方在於處於二維上的極值查找，可以套用線段樹來完成樹套樹達到二維效果。使用 zkw 式只是查找常數比較小，但記憶體會比較大，而基本線段樹的空間開銷是 $4 N^2$。接下來就要完成單點更新、區域極值查找。

時間複雜度 $O(N^2 + R \log^2 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
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
#include <stdio.h>
#include <string.h>
#include <algorithm>
using namespace std;
class D2RMQ {
public:
    struct Node {
        short mxv;
    };
    struct D2Tree {
        Node subtree[2048];
    } d2tree[2048];
    int Mmain, Msub;
    short ansMax;
    void init(int n, int m) {
        for (Mmain = 1; Mmain < n+2; Mmain <<= 1);
        for (Msub = 1; Msub < m+2; Msub <<= 1);
        build();
    }
    void build() {
        for (int i = (Mmain<<1)-1; i > 0; i--)
            sub_build(i);
    }
    void modify(int x, int y, int val) {
        static int s, t;
        d2tree[x+Mmain].subtree[y+Msub].mxv = val;
        t = x+Mmain;
        for (s = (y+Msub)>>1; s > 0; s >>= 1)
            d2tree[t].subtree[s].mxv = max(d2tree[t].subtree[s<<1].mxv, d2tree[t].subtree[s<<1|1].mxv);
        for (s = (x+Mmain)>>1; s > 0; s >>= 1)
            sub_modify(s, y);
    }
    void query(int lr, int rr, int lc, int rc) {
        static int s, t;
        ansMax = 0;
        if (lr <= rr && lc <= rc)
        for (s = lr+Mmain-1, t = rr+Mmain+1; (s^t) != 1;) {
            if (~s&1)
                sub_query(s^1, lc, rc);
            if (t&1)
                sub_query(t^1, lc, rc);
            s >>= 1, t >>= 1;
        }
    }	
private:
    void sub_build(int k) {
        memset(d2tree[k].subtree, 0, sizeof(Node) * Msub * 2);
//		for (int i = (Msub<<1)-1; i > 0; i--)
//			d2tree[k].subtree[i].mxv = 0;
    }
    void sub_modify(int k, int y) {
        for (int s = y+Msub; s > 0; s >>= 1)
            d2tree[k].subtree[s].mxv = max(d2tree[k<<1].subtree[s].mxv, d2tree[k<<1|1].subtree[s].mxv);
    }
    void sub_query(int k, int lc, int rc) {
        static int s, t;
        for (s = lc+Msub-1, t = rc+Msub+1; (s^t) != 1;) {
            if (~s&1)
                ansMax = max(ansMax, d2tree[k].subtree[s^1].mxv);
            if (t&1)
                ansMax = max(ansMax, d2tree[k].subtree[t^1].mxv);
            s >>= 1, t >>= 1;
        }
    }
} tool; 
char s1[1024], s2[1024];
char ss[32767];
int main() {
    int testcase, cases = 0;
    while (scanf("%s %s", s1+1, s2+1) == 2) {
        scanf("%d", &testcase);
        int n = strlen(s1+1), m = strlen(s2+1);
        int gap[256] = {};
        
        while (testcase--) {
            int GLCS = 0;
            for (int i = 0; i < 256; i++)
                gap[i] = 0x3f3f3f3f;
            while (scanf("%s", ss) == 1 && ss[0] != '$')
                scanf("%d", &gap[ss[0]]);
                
            tool.init(n, m);
            
            int K, lx, ly, rx, ry, dp;
            for (int i = 1; i <= n; i++) {
                for (int j = 1; j <= m; j++) {
                    if (s1[i] == s2[j]) {
                        K = gap[s1[i]];
                        lx = max(i-K-1, 1), ly = max(j-K-1, 1);
                        rx = i-1, ry = j-1;
                        
                        tool.query(lx, rx, ly, ry);
                        dp = tool.ansMax+1;
                        tool.modify(i, j, dp);
                        GLCS = max(GLCS, dp);
                    }
                }
            }
            printf("%d", GLCS);
            if (testcase)
                printf(" ");
        }
        puts("");
    }
    return 0;
}
Read More +
解題區/解題區 - Zerojudge

a243. 第四題：點燈遊戲

Problem

點燈遊戲，在二維地圖上，點擊一個燈泡後，亮會變暗、暗會變亮，還有周圍的燈泡也會改變狀態。滿足 $| x1 − x2 | + | y1 − y2 | = d$ 的位置的燈泡也會隨之變化，注意到是曼哈頓距離等於，而非小於等於。

給一個初始盤面，請問點燈結果有沒有解。

Sample Input

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
1 1 1
1
2 2 1
1 1
1 1
3 2 1
1 0 1
0 1 0
3 3 1
1 0 1
0 1 0
1 0 1
4 4 2
1 1 0 1
0 0 0 1
1 0 1 1
1 0 0 0
5 5 1
1 1 1 0 1
0 1 0 1 0
1 0 1 0 1
0 1 0 1 0
1 0 1 0 1
5 5 2
0 0 0 0 0
0 0 0 0 0
0 0 1 0 0
0 0 0 0 0
0 0 0 0 0
11 11 3
0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 1 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0
11 11 3
0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 1 1 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0
13 13 7
0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 1 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0

Sample Output

1
2
3
4
5
6
7
8
9
10
1
1
0
1
0
0
1
1
0
1

Solution

轉換成聯立方程組，當盤面是 $N \times M$，則用 $N \times M$ 個變數和 $N \times M$ 等式進行解方程。燈泡變換只有 1 -> 0 和 0 -> 1，利用高斯消去法解聯立時需要使用 XOR 取代原本的加減乘除計算。變數 $x(i, j) = 1$ 表示要 $(i, j)$ 按下，反之 0 表示不按，最後再模擬一次看能不能翻回全暗。

關於列方程式，下面是一個 $(2, 2)$ 位置被點開的影響周圍情況

1
2
3
4
0100
1110
0100
0000

相同地，周圍也可以影響到自己，那麼方程式就是 $x(1, 2) + x(2, 1) + x(2, 2) + x(2, 3) + x(3, 2) = L(2, 2)$。其中 + 號相當於 XOR 計算，而 L(i, j) 表示初始盤面狀態，當左式等於右式相當於 L(i, 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
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
#include <bits/stdc++.h> 
using namespace std;
class XOR_GAUSSIAN { // XOR Gaussian Elimination
public:
    static const int MAXN = 700;
    char mtx[MAXN][MAXN+1], varX[MAXN];
    int n;
    void compute() {
    	int row = 0, col = 0, arb = 0;
    	int equ = n, var = n;
    	while (row < equ && col < var) {
    		int c = row;
    		for (int i = row; i < equ; i++)
    			if (mtx[i][col])
    				c = i;
    		for (int i = 0; i <= var; i++)
                swap(mtx[c][i], mtx[row][i]);
    		if (mtx[row][col] == 0) {
    			col++, arb++;
    			continue;
    		}    		
    		for (int i = 0; i < equ; i++) {
    			if (i == row || mtx[i][col] == 0) continue;
    			for (int j = var; j >= 0; j--)
    				mtx[i][j] ^= mtx[row][j];
    		}
    		row++, col++;
    	}
    	memset(varX, 0, sizeof(varX));
    	for (int i = 0, j; i < equ; i++) {
    		if (mtx[i][var] == 0)
    			continue;
    		for (j = 0; j < var && mtx[i][j] == 0; j++);
    		varX[j] = mtx[i][var];
    	}
    }
    void init(int n) {
    	this->n = n;
    	for (int i = 0; i < n; i++)
    		for (int j = 0; j <= n; j++)
    			mtx[i][j] = 0;
    }
} gauss;
int main() {
    int n, m, d;
    int g[32][32];
    int cases = 0;
    while (scanf("%d %d %d", &m, &n, &d) == 3 && n) {
        for (int i = 0; i < n; i++)
            for (int j = 0; j < m; j++)
                scanf("%d", &g[i][j]);
        int N = n*m;
        gauss.init(N);
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < m; j++) {
                int x = i*m + j, tx, ty;
                gauss.mtx[x][N] = g[i][j];
                gauss.mtx[x][x] = 1;
                for (int dx = -d; dx <= d; dx++) {
                    int v = d - abs(dx);
                    int dy[2] = {-v, v}, dn = v == 0 ? 1 : 2;
                    for (int k = 0; k < dn; k++) {
                        tx = i + dx, ty = j + dy[k];
                        if (tx < 0 || ty < 0 || tx >= n || ty >= m)
                            continue;
                        gauss.mtx[x][tx*m+ty] = 1;
                    }
                }
            }
        }
        gauss.compute();
        
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < m; j++) {
                if (gauss.varX[i*m+j] == 0)
                    continue;
                g[i][j] ^= 1;
                int tx, ty;
                for (int dx = -d; dx <= d; dx++) {
                    int v = d - abs(dx);
                    int dy[2] = {-v, v}, dn = v == 0 ? 1 : 2;
                    for (int k = 0; k < dn; k++) {
                        tx = i + dx, ty = j + dy[k];
                        if (tx < 0 || ty < 0 || tx >= n || ty >= m)
                            continue;
                        g[tx][ty] ^= 1;
                    }
                }
            }
        }
        
        int ok = 1;
        for (int i = 0; i < n; i++)  {
            for (int j = 0; j < m; j++) { 
                ok &= g[i][j] == 0;
            }
        }
        printf("%d\n", ok);
    }
    return 0;
}
/*
2 3 1
1 0 1 0 0 0
3 3 2
1 0 0
0 0 0
0 0 0
*/
Read More +
解題區/解題區 - Zerojudge

b463. 用曲線調整色彩

Problem

給定三個曲線，每個曲線通過 n 個點，分別表示 r, g, b 的轉換圖示，請將影像的 rgb 替換成曲線給定的結果。

特別注意到曲線的最左、最右端點之外都是水平線。

Sample Input

1
2
3
4
5
2 0 0 255 255
2 0 0 255 255
2 0 0 255 255
1 1
128 128 128

Sample Output

1
2
1 1
128 128 128

Solution

處理最煩的就是通過 n 個點的曲線怎麼求，套用高斯消去法，假設多項式有 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
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
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <vector>
#include <algorithm>
using namespace std;
class GAUSSIAN { // Gaussian Elimination
public:
    static const int MAXN = 16;
    double mtx[MAXN][MAXN+1], var[MAXN];
    int exist[MAXN], n;
    vector< pair<double, double> > A;
    double f(double x) {
        if (x < A.front().first)
            return A.front().second;
        if (x > A.back().first)
            return A.back().second;
        double y = 0;
        for (int i = 0; i < n; i++)
            y += var[i] * pow(x, i);
        return y;
    }
    void add(pair<double, double> a) {
        A.push_back(a);
        sort(A.begin(), A.end());
    }
    void compute(int n) {
        const double eps = 1e-12;
        this->n = n;
        for (int i = 0, c; i < n; i++) {
            c = i;
            for (int j = i; j < n; j++)
                if (fabs(mtx[c][i]) < fabs(mtx[j][i]))
                    c = j;
            if (fabs(mtx[c][i]) < eps)
                continue;
            if (c != i) {
                for (int j = 0; j <= n; j++)
                    swap(mtx[c][j], mtx[i][j]);
            }
            for (int j = 0; j < n; j++) {
                if (i == j) continue;
                for (int k = n; k >= i; k--) {
                    mtx[j][k] -= mtx[j][i]/mtx[i][i]*mtx[i][k];
                }
            }
        }
        for (int i = 0; i < n; i++)
            exist[i] = 1;
        for (int i = n-1; i >= 0; i--) {
            if (fabs(mtx[i][i]) < eps) {
                exist[i] = 0;
                continue;
            }
            if (fabs(mtx[i][n]) < eps)
                var[i] = 0;
            else
                var[i] = mtx[i][n]/mtx[i][i];
            for (int j = i+1; j < n; j++)
                if (fabs(mtx[i][j]) > eps && exist[j])
                    exist[i] = 0;
        }
    }
};
class IMAGE {
public:
    struct Pixel {
        double r, g, b;
        Pixel(double x = 0, double y = 0, double z = 0):
        r(x), g(y), b(z) {}
        void read() {
            scanf("%lf %lf %lf", &r, &g, &b);
        }
        Pixel operator-(const Pixel &x) const {
            return Pixel(r-x.r, g-x.g, b-x.b);
        }
        Pixel operator+(const Pixel &x) const {
            return Pixel(r+x.r, g+x.g, b+x.b);
        }
        Pixel operator*(const double x) const {
            return Pixel(r*x, g*x, b*x);
        }
        Pixel operator*(const Pixel &x) const {
            return Pixel(r*x.r, g*x.g, b*x.b);
        }
        Pixel operator/(const double x) const {
            return Pixel(r/x, g/x, b/x);
        }
        bool operator==(const Pixel &x) const {
            return r == x.r && g == x.g && b == x.b;
        }
        void print() {
            printf("%d %d %d", (int)round(r), (int)round(g), (int)round(b));
        }
        void normal() {
            r = max(min(r, 255.0), 0.0);
            g = max(min(g, 255.0), 0.0);
            b = max(min(b, 255.0), 0.0);
        }
    };
    int W, H;
    static const int MAXN = 256;
    Pixel data[MAXN][MAXN];
    void read() {
        scanf("%d %d", &W, &H);
        for (int i = 0; i < H; i++)
            for (int j = 0; j < W; j++)
                data[i][j].read();
    }
    void print() {
        printf("%d %d\n", W, H);
        for (int i = 0; i < H; i++)
            for (int j = 0; j < W; j++)
                data[i][j].print(), printf("%c", j == W-1 ? '\n' : ' ');
    }
    void blind(GAUSSIAN func[]) {
        for (int i = 0; i < H; i++) {
            for (int j = 0; j < W; j++) {
                data[i][j].r = func[0].f(data[i][j].r);
                data[i][j].g = func[1].f(data[i][j].g);
                data[i][j].b = func[2].f(data[i][j].b);
                data[i][j].normal();
            }
        }
    }
} image;
int main() {
    GAUSSIAN func[3];
    for (int i = 0; i < 3; i++) {
        int gn;
        double gx, gfx;
        scanf("%d", &gn);
        for (int j = 0; j < gn; j++) {
            scanf("%lf %lf", &gx, &gfx);
            for (int k = 0; k < gn; k++)
                func[i].mtx[j][k] = pow(gx, k);
            func[i].mtx[j][gn] = gfx;
            func[i].add({gx, gfx});
        }
        func[i].compute(gn);
    }
    image.read();
    image.blind(func);
    image.print();
    return 0;
}
Read More +
解題區/解題區 - Zerojudge

b461. Fibonacci 之夢

Problem

給斐波那契数十進制的最後 $K$ 位$F_n$，反求該數字符合的最小項 $n$ 為何。

Sample Input

1
2
3
4
3 3
001
144
004

Sample Output

1
2
3
2
13
You've slept foolish!

Solution

前續 b443: 我愛 Fibonacci 的解法，得到 $\mod 10^K$ 的週期長度，分別求出 $K = 1, \; 2, \; 3, \cdots 17$，採用 b443 的作法，可以在 7 秒左右完成週期計算。

也可以選擇本地建表，得到週期長度分別為以下，直接打表使用。

1
UINT64 C[20] = {1, 60, 300, 1500, 15000, 150000, 1500000, 15000000, 150000000, 1500000000, 15000000000, 150000000000, 1500000000000, 15000000000000, 150000000000000, 1500000000000000, 15000000000000000, 150000000000000000, 1500000000000000000};

當要求尾 K 位的結果，要從尾 1 位的符合位置開始窮舉，隨後增加一個週期。例如現在已知尾 i 位符合的項數位置 `pos，那麼在位置 pos + C[i] 檢查是否有符合，並且滿足 pos + C[i] <= C[i+1]，之後就是循環週期不必窮舉。

由於週期之間的大小倍率不大，窮舉情況非常少，大約可以在數萬次內解決，為了加速數列計算，矩陣乘法需要 $O(\log n)$ 的時間，可以預處理相乘的矩陣，由於每一步跨越的項數是固定的，那麼計算量就會降到 $O(1)$。

關於模乘法

1
return (a*b - (long long)(a/(long double)mod*b+1e-3)*mod+mod)%mod;

由於 long double 有效位數 106-bits，用在 64-bits 模乘法非常夠用，取代模擬的加法計算，速度快個四到五倍。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <algorithm>
#include <limits.h>
#include <iostream>
#include <vector>
#include <map>
using namespace std;
#define MILLER_BABIN 10
typedef unsigned long long UINT64;
#define UNSTABLE_FAST
UINT64 mul(UINT64 a, UINT64 b, UINT64 mod) {
#ifdef UNSTABLE_FAST
    return (a*b - (long long)(a/(long double)mod*b+1e-3)*mod+mod)%mod;
#endif
    UINT64 ret = 0;
    for (a = a >= mod ? a%mod : a, b = b >= mod ? b%mod : b; b != 0; b>>=1, a <<= 1, a = a >= mod ? a - mod : a) {
        if (b&1) {
            ret += a;
            if (ret >= mod)
                ret -= mod;
        }
    }
    return ret;
}
struct Matrix {
    UINT64 v[2][2];
    int row, col; // row x col
    Matrix(int n, int m, int a = 0) {
        memset(v, 0, sizeof(v));
        row = n, col = m;
        for(int i = 0; i < row && i < col; i++)
            v[i][i] = a;
    }
    Matrix multiply(const Matrix& x, const long long mod) const {
        Matrix ret(row, x.col);
        for(int i = 0; i < row; i++) {
            for(int k = 0; k < col; k++) {
                if (!v[i][k])
                    continue;
                for(int j = 0; j < x.col; j++) {
                    ret.v[i][j] += mul(v[i][k], x.v[k][j], mod);
                    if (ret.v[i][j] >= mod)
                        ret.v[i][j] -= mod;
                }
            }
        }
        return ret;
    }
    Matrix pow(const long long& n, const long long mod) const {
        Matrix ret(row, col, 1), x = *this;
        long long y = n;
        while(y) {
            if(y&1)	ret = ret.multiply(x, mod);
            y = y>>1, x = x.multiply(x, mod);
        }
        return ret;
    }
} FibA(2, 2, 0);
#define MAXL (50000>>5)+1
#define GET(x) (mark[x>>5]>>(x&31)&1)
#define SET(x) (mark[x>>5] |= 1<<(x&31))
int mark[MAXL], P[50000], Pt = 0;
void sieve() {
    register int i, j, k;
    SET(1);
    int n = 46340;
    for (i = 2; i <= n; i++) {
        if (!GET(i)) {
            for (k = n/i, j = i*k; k >= i; k--, j -= i)
                SET(j);
            P[Pt++] = i;
        }
    }
}
UINT64 mpow(UINT64 x, UINT64 y, UINT64 mod) { // mod < 2^32
    UINT64 ret = 1;
    while (y) {
        if (y&1)
            ret = (ret * x)%mod;
        y >>= 1, x = (x * x)%mod;
    }
    return ret % mod;
}
UINT64 mpow2(UINT64 x, UINT64 y, UINT64 mod) {
    UINT64 ret = 1;
    while (y) {
        if (y&1)
            ret = mul(ret, x, mod);
        y >>= 1, x = mul(x, x, mod);
    }
    return ret;
}
void exgcd(long long x, long long y, long long &g, long long &a, long long &b) {
    if (y == 0)
        g = x, a = 1, b = 0;
    else
        exgcd(y, x%y, g, b, a), b -= (x/y) * a;
}
long long llgcd(long long x, long long y) {
    if (x < 0)    x = -x;
    if (y < 0)    y = -y;
    if (!x || !y)    return x + y;
    long long t;
    while (x%y)
        t = x, x = y, y = t%y;
    return y;
}
long long inverse(long long x, long long p) {
    long long g, b, r;
    exgcd(x, p, g, r, b);
    if (g < 0)	r = -r;
    return (r%p + p)%p;
}
int isPrime(long long p) { // implements by miller-babin
    if (p < 2 || !(p&1))	return 0;
    if (p == 2)				return 1;
    long long q = p-1, a, t;
    int k = 0, b = 0;
    while (!(q&1))	q >>= 1, k++;
    for (int it = 0; it < MILLER_BABIN; it++) {
        a = rand()%(p-4) + 2;
        t = mpow2(a, q, p);
        b = (t == 1) || (t == p-1);
        for (int i = 1; i < k && !b; i++) {
            t = mul(t, t, p);
            if (t == p-1)
                b = 1;
        }
        if (b == 0)
            return 0;
    }
    return 1;
}
long long pollard_rho(long long n, long long c) {
    long long x = 2, y = 2, i = 1, k = 2, d;
    while (true) {
        x = (mul(x, x, n) + c);
        if (x >= n)	x -= n;
        d = llgcd(x - y, n);
        if (d > 1) return d;
        if (++i == k) y = x, k <<= 1;
    }
    return n;
}
void factorize(int n, vector<long long> &f) {
    for (int i = 0; i < Pt && P[i]*P[i] <= n; i++) {
        if (n%P[i] == 0) {
            while (n%P[i] == 0)
                f.push_back(P[i]), n /= P[i];
        }
    }
    if (n != 1)	f.push_back(n);
}
void llfactorize(long long n, vector<long long> &f) {
    if (n == 1)
        return ;
    if (n < 1e+9) {
        factorize(n, f);
        return ;
    }
    if (isPrime(n)) {
        f.push_back(n);
        return ;
    }
    long long d = n;
    for (int i = 2; d == n; i++)
        d = pollard_rho(n, i);
    llfactorize(d, f);
    llfactorize(n/d, f);
}
// above largest factor
// ---------------------- //
int legendre_symbol(UINT64 d, UINT64 p) {
    if (d%p == 0)	return 0;
    return mpow2(d, (p-1)>>1, p) == 1 ? 1 : -1;
}
void factor_gen(int idx, long long x, vector< pair<long long, int> > &f, vector<long long> &ret) {
    if (idx == f.size()) {
        ret.push_back(x);
        return ;
    }
    for (long long i = 0, a = 1; i <= f[idx].second; i++, a *= f[idx].first)
        factor_gen(idx+1, x*a, f, ret);
}
void factor_gen(long long n, vector<long long> &ret) {
    vector<long long> f;
    vector< pair<long long, int> > f2;
    llfactorize(n, f);
    sort(f.begin(), f.end());
    int cnt = 1;
    for (int i = 1; i <= f.size(); i++) {
        if (i == f.size() || f[i] != f[i-1])
            f2.push_back(make_pair(f[i-1], cnt)), cnt = 1;
        else
            cnt ++;
    }
    factor_gen(0, 1, f2, ret);
    sort(ret.begin(), ret.end());
}
UINT64 cycleInFib(UINT64 p) {
    if (p == 2)	return 3;
    if (p == 3)	return 8;
    if (p == 5)	return 20;
    vector<long long> f;
    if (legendre_symbol(5, p) == 1)
        factor_gen(p-1, f);
    else
        factor_gen(2*(p+1), f);
    long long f1, f2;
    for (int i = 0; i < f.size(); i++) {
        Matrix t = FibA.pow(f[i]-1, p);
        f1 = (t.v[0][0] + t.v[0][1])%p;
        f2 = (t.v[1][0] + t.v[1][1])%p;
        if (f1 == 1 && f2 == 0)
            return f[i];
    }
    return 0;
}
UINT64 cycleInFib(UINT64 p, int k) {
    UINT64 s = cycleInFib(p);
    for (int i = 1; i < k; i++)
        s = s * p;
    return s;
}
// inverse fibonacci function
int main() {
    sieve();
    FibA.v[0][0] = 1, FibA.v[0][1] = 1;
    FibA.v[1][0] = 1, FibA.v[1][1] = 0;
    
    int K, testcase;
    UINT64 C[20] = {1, 60, 300, 1500, 15000, 150000, 1500000, 15000000, 150000000, 1500000000, 15000000000, 150000000000, 1500000000000, 15000000000000, 150000000000000, 1500000000000000, 15000000000000000, 150000000000000000, 1500000000000000000};
    UINT64 M[20];
    scanf("%d %d", &K, &testcase);
    for (int i = 1; i <= K; i++) {
        long long m = mpow2(10, i, LONG_LONG_MAX);
        //        vector<long long> f;
        //        map<long long, int> r;
        //
        //        llfactorize(m, f);
        //        for (auto &x : f)
        //            r[x]++;
        //
        //        UINT64 cycle = 1;
        //        for (auto &x : r) {
        //            UINT64 t = cycleInFib(x.first, x.second);
        //            cycle = cycle / llgcd(t, cycle) * t;
        //        }
        //        C[i] = cycle;
        M[i] = m;
    }
    
    while (testcase--) {
        char s[32];
        scanf("%s", s);
        long long fib;
        vector<UINT64> pos;
        sscanf(s, "%lld", &fib);
        if (fib == 0) {
            puts("1");
            continue;
        }
        sscanf(s+K-1, "%lld", &fib);
        for (int i = 0; i < C[1]; i++) {
            Matrix t = FibA.pow(i, 10);
            long long fn = t.v[1][0];
            if (fn == fib)
                pos.push_back(i);
        }
        for (int i = 2; i <= K; i++) {
            sscanf(s+K-i, "%lld", &fib);
            vector<UINT64> next;
            for (auto &x : pos) {
                Matrix t = FibA.pow(x, M[i]), tp = FibA.pow(C[i-1], M[i]);
                for (UINT64 c = x; c <= C[i]; c += C[i-1]) {
                    long long fn = t.v[1][0];
                    if (fn == fib)
                        next.push_back(c);
                    t = t.multiply(tp, M[i]);
                }
            }
            sort(next.begin(), next.end());
            next.resize(unique(next.begin(), next.end()) - next.begin());
            pos = next;
        }
        if (pos.size() == 0)
            puts("You've slept foolish!");
        else
            printf("%llu\n", pos[0]+1);
    }
    return 0;
}
/*
 17 1
 00012779675210023
 16 1
 0012569221817747
 */
Read More +
解題區/解題區 - Zerojudge

b452. 傻傻地幫人數錢錢

Problem

給一個彩色影像上面有數個相同大小的硬幣，硬幣之間不會重疊，但會有部分碰觸，請問影像中有多少個硬幣。

Sample Input

1
2
3
4
5
6
7
8
9
10
10 9
159 138 133 150 129 124 153 132 127 154 133 128 151 132 126 151 132 126 152 133 129 152 133 129 143 125 125 136 118 118
147 126 121 163 142 137 146 125 120 86 65 60 60 41 35 57 38 32 68 49 45 111 92 88 143 125 123 138 120 118
145 125 118 152 132 125 63 44 37 53 34 27 54 35 28 55 36 29 56 39 32 50 32 28 105 87 83 132 114 110
150 130 123 106 86 79 53 34 27 57 38 31 53 34 27 64 45 38 55 38 31 56 39 32 63 45 41 131 113 109
151 132 125 87 68 61 52 33 26 56 37 30 58 41 33 49 32 24 57 40 33 53 36 29 43 28 23 131 116 111
136 117 111 103 84 78 56 37 31 59 40 34 54 37 30 57 40 33 47 30 23 53 35 31 54 39 34 141 126 121
142 124 120 140 122 118 52 34 30 54 36 32 51 33 29 53 35 31 53 38 33 46 31 28 82 67 64 136 121 118
145 127 127 125 107 107 116 98 98 54 36 36 53 35 33 54 36 34 46 30 30 65 49 49 138 122 122 145 129 129
135 119 120 137 121 122 137 121 122 133 117 118 103 87 88 95 79 80 111 95 96 138 122 123 132 118 118 132 118 118

Sample Output

1
1

Solution

若圓形彼此之間不相交，可以用二值化 + 灌水法 (flood fill) + 分團大小檢測。現在圓形有相連可能，出題者給我附前測代碼啊，咱們兩個比一下速度 … 總之步驟是

  1. 根據亮度二值化
  2. 搭配索貝爾運算 (sobel) 選定閥值後找到邊緣
  3. 接著窮舉圓半徑，使用霍夫轉換將邊緣點推到同一個圓心
  4. 掃描一個 $7 \times 7$ 的矩形，內部點個數要出現 56 個，同時要滿足 49 格至少出現 36 格或者其中一格出現大於 4 次。

關於霍夫轉換，其中$x_c, \; y_c$ 表示推向圓心的座標，而 $r$ 表示窮舉的圓半徑，$gx$ 表示 x 方向的 sobel 差分，同理 $gy$，而 $g = \sqrt{gx^2 + gy^2}$。

$x_c = x - r \times (gx / g)$ $y_c = y - r \times (gy / g)$

當然閥值判定都還不到位，手動測試好幾個版本，OpenCV 的寫法值得去 trace 一下。代碼僅供玩玩，不代表在其他情況也能使用。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <stdlib.h>
#include <limits.h>
#include <vector>
#include <algorithm>
using namespace std;
class IMAGE {
public:
    struct Pixel {
        int r, g, b;
        Pixel(int x = 0, int y = 0, int z = 0):
        r(x), g(y), b(z) {}
        void read() {
            scanf("%d %d %d", &r, &g, &b);
        }
        Pixel operator-(const Pixel &x) const {
            return Pixel(r-x.r, g-x.g, b-x.b);
        }
        Pixel operator+(const Pixel &x) const {
            return Pixel(r+x.r, g+x.g, b+x.b);
        }
        Pixel operator*(const double x) const {
            return Pixel(r*x, g*x, b*x);
        }
        Pixel operator/(const double x) const {
            return Pixel(r/x, g/x, b/x);
        }
        bool operator==(const Pixel &x) const {
            return r == x.r && g == x.g && b == x.b;
        }
        void print() {
            printf("%3d", length());
        }
        int sum() {
            return r + g + b;
        }
        int length() {
            return abs(r) + abs(g) + abs(b);
        }
        int dist(Pixel x) {
            return abs((r + g + b) - (x.r + x.g + x.b));
        }
    };
    int W, H;
    static const int MAXN = 256;
    Pixel data[MAXN][MAXN];
    void read() {
        scanf("%d %d", &W, &H);
        for (int i = 0; i < H; i++)
            for (int j = 0; j < W; j++)
                data[i][j].read();
    }
    void print() {
        printf("%d %d\n", W, H);
        for (int i = 0; i < H; i++)
            for (int j = 0; j < W; j++)
                data[i][j].print(), printf("%c", j == W-1 ? '\n' : ' ');
    }
    inline Pixel getPixel(int x, int y) {
        if (x >= 0 && y >= 0 && x < H && y < W)
            return data[x][y];
        if (y < 0)  return data[min(max(x, 0), H-1)][0];
        if (y >= W) return data[min(max(x, 0), H-1)][W-1];
        if (x < 0)  return data[0][min(max(y, 0), W-1)];
        if (x >= H) return data[H-1][min(max(y, 0), W-1)];
        return Pixel(0, 0, 0);
    }
    void sobel(int i, int j, int &gx, int &gy) {
        static int dx[] = {-1, -1, -1, 0, 0, 0, 1, 1, 1};
        static int dy[] = {-1, 0, 1, -1, 0, 1, -1, 0, 1};
        static int yw[] = {-1, 0, 1, -2, 0, 2, -1, 0, 1};
        static int xw[] = {-1, -2, -1, 0, 0, 0, 1, 2, 1};
        Pixel Dx(0, 0, 0), Dy(0, 0, 0);
        for (int k = 0; k < 9; k++) {
            if (xw[k])
                Dx = Dx + getPixel(i+dx[k], j+dy[k]) * xw[k];
            if (yw[k])
                Dy = Dy + getPixel(i+dx[k], j+dy[k]) * yw[k];
        }
        gx = Dx.sum(), gy = Dy.sum();
    }
    int used[MAXN][MAXN];
    int gx[MAXN][MAXN], gy[MAXN][MAXN];
    double gxy[MAXN][MAXN];
    int isValid(int x, int y) {
        return x >= 0 && y >= 0 && x < H && y < W;
    }
    int hough_circle() {
        int mxb = INT_MIN, mnb = INT_MAX;
        for (int i = 0; i < H; i++) {
            for (int j = 0; j < W; j++) {
                int b = data[i][j].length();
                mxb = max(mxb, b), mnb = min(mnb, b);
            }
        }
        if (mxb - mnb < 300)
            return 0;
        int threshold = 250;
        for (int i = 0; i < H; i++) {
            for (int j = 0; j < W; j++) {
                if (data[i][j].length() >= threshold)
                    data[i][j] = Pixel(1, 1, 1);
                else
                    data[i][j] = Pixel(0, 0, 0);
            }
        }
        
        for (int i = 0; i < H; i++) {
            for (int j = 0; j < W; j++) {
                sobel(i, j, gx[i][j], gy[i][j]), gxy[i][j] = hypot(gx[i][j], gy[i][j]);
            }
        }
        
        int ret = 0;
        for (double r = 4; r <= min(H, W)/2; r += 1) {
            for (int i = 0; i < H; i++)
                for (int j = 0; j < W; j++)
                    used[i][j] = 0;
            for (int i = 0; i < H; i++) {
                for (int j = 0; j < W; j++) {
                    if (gxy[i][j] > 4) {
                        int xc, yc;
                        xc = round(i - r * (gx[i][j] / gxy[i][j]));
                        yc = round(j - r * (gy[i][j] / gxy[i][j]));
                        if (isValid(xc, yc))
                            used[xc][yc]++;
                    }
                }
            }
            
            int coins = 0;
            const int C = 3;
            for (int i = 0; i < H; i++) {
                for (int j = 0; j < W; j++) {
                    if (used[i][j] < 3)
                        continue;
                    int sum = 0, mx = 0, has = 0;
                    for (int p = -C; p <= C; p++) {
                        for (int q = -C; q <= C; q++) {
                            if (isValid(i+p, j+q))
                                sum += used[i+p][j+q], mx = max(mx, used[i+p][j+q]), has += used[i+p][j+q] > 0;
                        }
                    }
                    if (sum > 56 && (has > 36 || mx > 4)) {
                        coins++;
                        int cx = i, cy = j;
                        for (int p = -r-1; p <= r+1; p++) {
                            for (int q = -r-1; q <= r+1; q++) {
                                if (isValid(cx+p, cy+q))
                                    used[cx+p][cy+q] = 0;
                            }
                        }
                    }
                }
            }
            if (coins < ret - 2)
                break;
            ret = max(ret, coins);
        }
        return ret;
    }
} image;
int main() {
    image.read();
    printf("%d\n", image.hough_circle());
    return 0;
}
Read More +
解題區/解題區 - Zerojudge

b448. 哈哈鏡

Problem

進行圖片變形，效果類似哈哈鏡的作用。

對圖片水平中線進行座標 y' = sqrt(y) 變換，將靠近中線的像素盡可能地拉往中間，形成延伸的變形效果。

Sample Input

1
2
1 1
1 2 3

Sample Output

1
2
1 1
1 2 3

Solution

計算公式

  • 中線之上 (H/2) - pow(i-H/2, 2)/ (H/2)
  • 中線之下 pow(i-H/2, 2)/ (H/2) + H/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
#include <bits/stdc++.h>
using namespace std;
class IMAGE {
public:
    struct Pixel {
        int r, g, b;
        Pixel(int x = 0, int y = 0, int z = 0):
        r(x), g(y), b(z) {}
        void read() {
            scanf("%d %d %d", &r, &g, &b);
        }
        Pixel operator-(const Pixel &x) const {
            return Pixel(r-x.r, g-x.g, b-x.b);
        }
        Pixel operator+(const Pixel &x) const {
            return Pixel(r+x.r, g+x.g, b+x.b);
        }
        Pixel operator*(const double x) const {
            return Pixel(r*x, g*x, b*x);
        }
        Pixel operator/(const double x) const {
            return Pixel(r/x, g/x, b/x);
        }
        bool operator==(const Pixel &x) const {
            return r == x.r && g == x.g && b == x.b;
        }
        void print() {
            printf("%d %d %d", r, g, b);
        }
        int sum() {
            return r + g + b;
        }
        int length() {
            return abs(r) + abs(g) + abs(b);
        }
        int dist(Pixel x) {
            return abs((r + g + b) - (x.r + x.g + x.b));
        }
    };
    int W, H;
    static const int MAXN = 256;
    Pixel data[MAXN][MAXN], tmp[MAXN][MAXN];
    void read() {
        scanf("%d %d", &W, &H);
        for (int i = 0; i < H; i++)
            for (int j = 0; j < W; j++)
                data[i][j].read();
    }
    void print() {
        printf("%d %d\n", W, H);
        for (int i = 0; i < H; i++)
            for (int j = 0; j < W; j++)
                data[i][j].print(), printf("%c", j == W-1 ? '\n' : ' ');
    }
    int isValid(int x, int y) {
        return x >= 0 && y >= 0 && x < H && y < W;
    }
    void distorting_mirror() {
        int rW = W, rH = H;
        double ch = H / 2.0;
        for (int i = 0; i < rH; i++) {
            double x, y;
            if (i < ch)
                x = ch - pow(i-ch, 2)/ch;
            else
                x = pow(i-ch, 2)/ch + ch;
            for (int j = 0; j < rW; j++) {
                y = j;
                int lx, rx, ly, ry;
                lx = floor(x), rx = ceil(x);
                ly = floor(y), ry = ceil(y);
                int px[] = {lx, lx, rx, rx};
                int py[] = {ly, ry, ly, ry};
                int c = -1;
                double mndist = 1e+30;
                for (int k = 0; k < 4; k++) {
                    if (!isValid(px[k], py[k]))
                        continue;
                    double d = (x-px[k])*(x-px[k])+(y-py[k])*(y-py[k]);
                    if (c == -1 || mndist > d)
                        c = k, mndist = d;
                }
                assert (c >= 0);
                tmp[i][j] = data[px[c]][py[c]];
            }
        }
        W = rW, H = rH;
        for (int i = 0; i < H; i++) {
            for (int j = 0; j < W; j++) {
                data[i][j] = tmp[i][j];
            }
        }
    }
} image;
int main() {
    image.read();
    image.distorting_mirror();
    image.print();
    return 0;	
}
Read More +
解題區/解題區 - Zerojudge

a994. 10325 - The Lottery

Problem

模擬計算，刪除組合數字的倍數，請問剩下多少個數字可選。

Sample Input

1
2
3
4
5
6
10 2
2 3
20 2
2 4
100 3
3 5 7

Sample Output

1
2
3
3
10
45

Solution

這一題是非常容易的題目，利用排容原理可以在 $O(2^m)$ 的時間內完成，所以複雜度沒有太大的改善，若使用 bitmask 的方式撰寫，複雜度會落在 $O(m \times 2^m)$，中間會有大量的 gcd(a, b) 計算，歐基里德輾轉相除法的常數並不大，時間複雜度 $O(\log n)$。

為了加速運算，可以利用組合來完成，利用選用組合 1011，可以得到 lcm(1011) = lcm(lcm(1010), A[0]) 完成，因此只會有 $O(2^m)$ 次 gcd() 的成本，整整少了常數 m。

因此需要使用 lowbit = i&(-i) 的位元運算技巧，同時為了得到 2 的冪次的次方數，建議使用內置函數 __builtin 系列，編譯器最佳化計算。而 gcd 使用，也使用內建的 __gcd(a, b) 但內置函數通常只會設置在 unsigned int 意即 32-bits 無號整數，要防止運算 overflow。

有人會提案使用建表，這樣查找只需要 $O(1)$，但記憶體置換會造成負擔，有兩個 $O(2^{15})$ 的 cache page，而且還要前置建表的計算成本。根據實驗測試，建表的方案速度會比較慢。

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
#include <bits/stdc++.h>
using namespace std;
int main() {
    int n, m, A[16];
    unsigned int dp[1<<15], a, b;
    while (scanf("%d %d", &n, &m) == 2) {
        for (int i = 1; i <= m; i++)
            scanf("%d", &A[i]);
        int ret = n;
        dp[0] = 1;
        for (int i = 1; i < (1<<m); i++) {
            long long val;
            a = dp[i-(i&(-i))], b = A[__builtin_ffs(i&(-i))];
            if (a > n) {dp[i] = n+1; continue;}
            val = b / __gcd(a, b) * (long long) a;
            if (val <= n) {
                dp[i] = val;
                if (__builtin_popcount(i)&1)
                    ret -= n / val;
                else
                    ret += n / val;
            } else {
                dp[i] = n+1;
            }
        }
        
        printf("%d\n", ret);
    }
    return 0;
}
Read More +