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

PTC 201506 六月小結

雖然有在 Facebook 上面討論，遲遲沒有貼出來。由於代碼有點長，所以就放在 Github 上進行瀏覽，在此只提供文字說明，題解代碼 點我。

Problem A Magic Box

給定最多 n = 30 個數字，任意挑選至少一個數字相乘得到 P，求 P mod 1000000007 = K 的組合數有多少。

背景知識

  • 中間相遇法
  • 乘法反元素

藉由中間相遇法，拆成 A, B 兩堆，每一堆最多有 n/2 個數字，分別計算組合數的乘積結果的組合數 P mod 1000000007。窮舉 A 那一堆其中一個組合乘積為 PA，方法數為 CNT_PA，而為了使 PA * PB mod 1000000007 = K，利用乘法反元素，找到 PB = K * PA^-1 mod 1000000007，利用映射找到 CNT_PB，因此累計方法數就會增加 CNT_PA * CNT_PB。

特別小心 K = 0 要例外處理。

Problem B Hidden Lines

平面上給定一堆線段，請問由上而下看時，有多少被藏住的線段。

背景知識

  • 分而治之
  • 計算幾何
  • 掃描線算法

將線段按照左端點的 x 座標值排序，接著對 x 軸座標進行分治，分治時維護區間內的天際線，合併兩個天際線的算法最慘需要 $O(n \log n)$，若不使用掃描線算法可以降到 $O(n)$，這類似合併排序，將兩個曲折的線段合併成一個。

最後遞歸 $T(n) = 2T(n/2) + O(n) = O(n \log n)$ 時間複雜度可以接受，但是代碼複雜度就很高。

Problem C Bombs

有兩種炸彈分別是單核和雙核炸彈，每一個核心都有兩個感測器，感測器又有兩種圓形、方形。若炸彈不會引爆，則至少一個核心被破壞。一個核心被破壞，則兩個感測器處於不運作。

圓形感測器不運作，則表示通過此感測器的線路都必須被剪斷。方形感測器不運作，則表示通過此感測器的線路保持正常。給予線路連接數個炸彈的感測器情況，請問是否存在一種剪斷線路的方案，滿足炸彈都不會被引爆。

背景知識

  • 2-SAT

變形的 2-SAT 問題，對於每一個線段剪或不剪都是一個未知數，求方程式是否有解。對於單核炸彈而言，那個核心必須被破壞，保證通過兩個感測器的線路都必須符合不被運作的情況，一開始就指派線路剪斷或不能被剪斷。對於雙核炸彈而言，當一個核心正常，則另一個核心必然要被破壞，反之亦然，則可以 核心 A 的某條線路若正常運作導致 B 的某條線路是壞的 意即 A -> NOT B 的邏輯方程。

最後套用 2-SAT 問題進行矛盾檢查即可。

Problem D Quick Reversal

數個區間反轉操作，求奇數位數字總和

背景知識

  • 模擬
  • 離散處理

區間離散 + 暴力 $O(Q^2)$。區間反轉可以用 splay tree 或者是塊狀表加速，根據 SGU 187 的經驗，不會卡到 $O(Q \log Q)$

Problem E Milk Delivery

求數條 s-t 路徑恰好經過 K 條路，路上權值和最大。

背景知識

  • 矩陣乘法
  • 圖論

建表 $O(N^3 \log K)$，詢問 $O(1)$。

Read More +
解題區/解題區 - Zerojudge

a481. 樹的維護

Problem

操作有三種

  • 1 x y 詢問樹路徑 x - y 經過的點權重和
  • 2 x w 修改 x 的點權為 w
  • 3 x y 將 x 父節點修改成 y

Sample Input

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

Sample Output

1
2
3
11
12
10

Solution

題目看起來是一個有根樹，由於權重落在點上，又由於父節點修改可以用額外數組紀錄，因此可以當作無向樹操作。

假設知道一條邊的兩端點 e(u, v)，進行 cut(u, v) 操作，不知道誰父誰子的切割時，可以轉兩次通到樹根，讓另一端點落在左子樹。

1
2
3
4
5
void cut(Node *x, Node *y) {
    mk_root(x);
    access(y), splay(y);
    y->ch[0] = x->fa = Node::EMPTY;
}

由於是無向樹，詢問路徑則藉由將其中一端變成根，接著藉由 access() 將路徑接到樹根，再藉由 splay() 轉到 splay tree 的根，此時帶權值就是路徑總和。

1
2
3
4
5
int sumPath(Node *x, Node *y) {
    mk_root(x);
    access(y), splay(y);
    return y->sum;
}
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
#include <bits/stdc++.h> 
using namespace std;
class LCT { // Link-Cut Tree
public:
    static const int MAXN = 131072;
    struct Node {
        static Node *EMPTY;
        Node *ch[2], *fa;
        int rev;
        int val, sum, size;
        Node() {
            ch[0] = ch[1] = fa = NULL;
            rev = 0;
            val = sum = 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() {
            sum = val, size = 1;
            if (ch[0] != EMPTY)
                sum += ch[0]->sum, size += ch[0]->size;
            if (ch[1] != EMPTY)
                sum += ch[1]->sum, size += ch[1]->size;
        }
        inline void push_deal(int c) {
            val = c;
            pushup();
        }
    } _mem[MAXN];
    
    int bufIdx;
    LCT() {
        Node::EMPTY = &_mem[0];
        Node::EMPTY->fa = Node::EMPTY->ch[0] = Node::EMPTY->ch[1] = Node::EMPTY;
        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();
    }
    void access(Node *u) {
        Node *v = Node::EMPTY;
        for (; u != Node::EMPTY; u = u->fa) {
            splay(u);
            u->ch[1] = v;
            v = u;
            v->pushup();
        }
    }
    void mk_root(Node *u) {
        access(u), splay(u);
        u->rev ^= 1;
    }
    void cut(Node *x, Node *y) {
        mk_root(x);
        access(y), splay(y);
        y->ch[0] = x->fa = Node::EMPTY;
    }
    void link(Node *x, Node *y) {
        mk_root(x);
        x->fa = y;
    }
    Node* find(Node *x) {
        access(x), splay(x);
        for (; x->ch[0] != Node::EMPTY; x = x->ch[0]);
        return x;
    }
    //
    void dealNode(Node *x, int c) {
        mk_root(x);
        x->push_deal(c);
    }
    int sumPath(Node *x, Node *y) {
        mk_root(x);
        access(y), splay(y);
        return y->sum;
    }
} tree;
LCT::Node *LCT::Node::EMPTY;
LCT::Node *A[131072], *node_x, *node_y;
int p[131072];
int main() {
    int n, m, x, y, c, u, v, cmd;
    while (scanf("%d", &n) == 1) {
        tree.init();
        for (int i = 1; i <= n; i++)
            A[i] = tree.newNode();
        
        for (int i = 1; i <= n; i++) {
            scanf("%d %d", &p[i], &y);
            if (p[i])
                tree.link(A[i], A[p[i]]);
            tree.dealNode(A[i], y);
        }
        
        scanf("%d", &m);
        for (int i = 0; i < m; i++) {
            scanf("%d %d %d", &cmd, &x, &y);
            if (cmd == 1) {
                printf("%d\n", tree.sumPath(A[x], A[y]));
            } else if (cmd == 2) {
                tree.dealNode(A[x], y);
            } else if (cmd == 3) {
                tree.cut(A[x], A[p[x]]);
                tree.link(A[x], A[y]);
                p[x] = y;
            }
        }
    }
    return 0;
}
Read More +
解題區/解題區 - 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 +
解題區/解題區 - 其他題目

BZOJ 1047 [HAOI2007] 理想的正方形

Problem

題目連結

在一個 $N \times M$ 的矩形中，找到一個 $K \times K$ 的正方形區域，求所有正方形區域中最大值減去最小值的差最小為何。

Sample Input

1
2
3
4
5
6
5 4 2
1 2 5 6
0 17 16 0
16 17 2 1
2 10 2 1
1 2 2 2

Sample Output

1
1

Solution

明顯地是一個滑動窗口，對於一維情況很明顯地直接套用單調隊列即可。

對於二維的滑動窗口，則可以對每一列維護一個 列 的單調隊列，對於某一行掃描時，再維護一個 行 的單調隊列，往右移動窗口時，把移動到 列 的單調隊列的極值加入，並且把超出窗口的列的極值移除。

由上至下、由左而右掃描每一行上的元素，找到以此元素為正方形右下角的正方形區域最大最小值。時間複雜度 $O(N^2)$，怕重複的代碼太多，直接寫成模板搭配 std::deque<> 來使用，但這樣的寫法會比純數組來得慢 (method interface 和 cache 的影響 …)。

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
#include <bits/stdc++.h> 
using namespace std;
const int MAXN = 1024;
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 cmp1(int a, int b) { return a > b; }
bool cmp2(int a, int b) { return a < b; }
inline int readchar() {
    const int N = 1048576;
    static char buf[N];
    static char *p = buf, *end = buf;
    if(p == end) {
        if((end = buf + fread(buf, 1, N, stdin)) == buf) return EOF;
        p = buf;
    }
    return *p++;
}
inline int ReadInt(int *x) {
    static char c, neg;
    while((c = readchar()) < '-')    {if(c == EOF) return 0;}
    neg = (c == '-') ? -1 : 1;
    *x = (neg == 1) ? c-'0' : 0;
    while((c = readchar()) >= '0')
        *x = (*x << 3) + (*x << 1) + c-'0';
    *x *= neg;
    return 1;
}
int main() {
    int N, M, K, val;
    while (ReadInt(&N)) {
        ReadInt(&M);
        ReadInt(&K);
                
        MonoQueue<int, cmp1> mxC[MAXN];
        MonoQueue<int, cmp2> mnC[MAXN];
        for (int i = 0; i < M; i++)
            mxC[i].setWindow(K), mnC[i].setWindow(K);
        int ret = INT_MAX;
        
        for (int i = 0; i < N; i++) {
            MonoQueue<int, cmp1> mxR;
            MonoQueue<int, cmp2> mnR;
            mxR.setWindow(K), mnR.setWindow(K);
            for (int j = 0; j < M; j++) {
                ReadInt(&val);				
                mxC[j].add(i, val);
                mxR.add(j, mxC[j].top().second);
                
                mnC[j].add(i, val);
                mnR.add(j, mnC[j].top().second);
                
                if (i >= K-1 && j >= K-1)
                    ret = min(ret, mxR.top().second - mnR.top().second);
            }
        }
        printf("%d\n", ret);
    }
    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 +