解題區/解題區 - Zerojudge

b434. 圖片的梯度

Problem

可以參考 中央大學影像處理實驗室 課程 講義第八章 邊擷取

邊緣擷取可以靠相鄰像素的差異來判斷是否是邊緣，若顏色差異很大則表示是邊。然而圖片像素是不連續的離散點，因此沒有辦法利用微積分學到的微分來得精準，而是要採用差分 (減法)，至於遮罩要取多大是個問題，一般而言還會搭配降噪的矩陣來完成。

如 Prewitt operator、Sobel operator、Kirsch operator … 等，矩陣都有搭配降噪處理。之所以用矩陣表示，是為了讓 GPU 快速地計算數值。

Sample Input

1
2
3
2 2
0 0 0 255 255 255
255 255 255 255 255 255

Sample Output

1
2
3
2 2
255 255 255 0 0 0
0 0 0 0 0 0

Solution

按照公式模擬即可。

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
#include <bits/stdc++.h> 
using namespace std;
class IMAGE {
public:
    struct Pixel {
        double r, g, b, a;
        Pixel(double x = 0, double y = 0, double z = 0, double w = 0):
            r(x), g(y), b(z), a(w) {}
        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, a-x.a);
    	}
    	Pixel operator+(const Pixel &x) const {
        	return Pixel(r+x.r, g+x.g, b+x.b, a+x.a);
    	}
    	Pixel operator*(const double x) const {
        	return Pixel(r*x, g*x, b*x, a*x);
    	}
        Pixel operator/(const double x) const {
        	return Pixel(r/x, g/x, b/x, a/x);
    	}
    	void print() {
    		printf("%d %d %d", (int)round(r), (int)round(g), (int)round(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();
    }
    Pixel pabs(Pixel x) {
        return Pixel(fabs(x.r), fabs(x.g), fabs(x.b));
    }
    void border() {
        for (int i = 0; i < H; i++) {
            for (int j = 0; j < W; j++) {
                Pixel dx(0, 0, 0), dy(0, 0, 0);
                if (i+1 < H)
                    dx = data[i+1][j] - data[i][j];
                if (j+1 < W)
                    dy = data[i][j+1] - data[i][j];
                dx = pabs(dx), dy = pabs(dy);
                data[i][j] = (dx + dy)/2.0;
            }
        }
    }
    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' : ' ');
    }
} test;
int main() {
    test.read();
    test.border();
    test.print();
    return 0;
}
Read More +
解題區/解題區 - Zerojudge

b427. 漸層色彩

Problem

給予兩個顏色的像素，在一張空白影像中使用漸層效果。

  • 水平由左到右的漸層
  • 從中心點擴散的漸層

Sample Input

1
3 3 0 0 0 0 255 255 255 255 255

Sample Output

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

Solution

線性內插，特別小心中心點 (x, y) 會是一個浮點數情況。

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
#include <bits/stdc++.h> 
using namespace std;
class IMAGE {
public:
    struct Pixel {
        double r, g, b, a;
        Pixel(double x = 0, double y = 0, double z = 0, double w = 0):
            r(x), g(y), b(z), a(w) {}
        void read() {
            scanf("%lf %lf %lf %lf", &r, &g, &b, &a);
        }
        Pixel operator-(const Pixel &x) const {
        	return Pixel(r-x.r, g-x.g, b-x.b, a-x.a);
    	}
    	Pixel operator+(const Pixel &x) const {
        	return Pixel(r+x.r, g+x.g, b+x.b, a+x.a);
    	}
    	Pixel operator*(const double x) const {
        	return Pixel(r*x, g*x, b*x, a*x);
    	}
        Pixel operator/(const double x) const {
        	return Pixel(r/x, g/x, b/x, a/x);
    	}
    	void print() {
    		printf("%d %d %d %d", (int)round(r), (int)round(g), (int)round(b), (int)round(a));
    	}
    };
    int W, H;
    static const int MAXN = 300;
    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 alloc() {
        int TYPE;
        Pixel st, ed;
        scanf("%d %d %d", &W, &H, &TYPE);
        st.read(), ed.read();
        gradient(st, ed, TYPE);
    }
    void gradient(Pixel st, Pixel ed, int TYPE) {
        if (TYPE != 0 && TYPE != 1)
            return;
        if (TYPE == 0) {
            for (int i = 0; i < H; i++)
                for (int j = 0; j < W; j++) 
                    data[i][j] = W-1 ? st + (ed - st) * j / (W-1) : st;
        } else {
            double cx = (H-1)/2.0, cy = (W-1)/2.0;
            double cr = hypot(cx, cy);
            for (int i = 0; i < H; i++)
                for (int j = 0; j < W; j++)
                    data[i][j] = cr ? st + (ed - st) * hypot(i-cx, j-cy) / cr : st;
        }
    }
    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' : ' ');
    }
} test;
int main() {
    test.alloc();
    test.print();
    return 0;
}
Read More +
解題區/解題區 - Zerojudge

b426. 宇宙光明體

Problem

模擬兩張圖片的合成。

將前景圖片疊在背景圖片之上，按照權重比例來融合重疊部分得像素。

Sample Input

1
2
3
4
5
0 0 0.5
1 1
255 255 255
1 1
0 0 0

Sample Output

1
2
1 1
128 128 128

Solution

原來是宇宙光明體，還以為是宇宙大覺者呢。

按照公式純模擬即可。

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
#include <bits/stdc++.h> 
using namespace std;
const double eps = 0;
class IMAGE {
public:
    struct Pixel {
        long double r, g, b, a;
        Pixel(long double x = 0, long double y = 0, long double z = 0, long double w = 0):
            r(x), g(y), b(z), a(w) {}
        void read() {
            int p1, p2, p3;
            scanf("%d %d %d", &p1, &p2, &p3);
            r = p1, g = p2, b = p3;
        }
        Pixel operator-(const Pixel &x) const {
        	return Pixel(r-x.r, g-x.g, b-x.b, a-x.a);
    	}
    	Pixel operator+(const Pixel &x) const {
        	return Pixel(r+x.r, g+x.g, b+x.b, a+x.a);
    	}
    	Pixel operator*(const long double x) const {
        	return Pixel(r*x, g*x, b*x, a*x);
    	}
        Pixel operator/(const long double x) const {
        	return Pixel(r/x, g/x, b/x, a/x);
    	}
    	void print() {
    		double p1 = r, p2 = g, p3 = b, p4 = a;
    		printf("%.0lf %.0lf %.0lf", p1 + eps, p2 + eps, p3 + eps);
    	}
    };
    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 add(IMAGE &other, int X, int Y, double OPACITY) {
        for (int i = 0; i < other.H && i+X < H; i++) {
            for (int j = 0; j < other.W && j+Y < W; j++) {
                if (i+X >= 0 && j+Y >= 0)
                    data[i+X][j+Y] = other.data[i][j] * OPACITY + data[i+X][j+Y] * (1 - OPACITY);
            }
        }
    }
    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' : ' ');
    }
} back, fore;
int main() {
    double OPACITY;
    int X, Y;
    scanf("%d %d %lf", &X, &Y, &OPACITY);
    fore.read();
    back.read();
    back.add(fore, Y, X, OPACITY);
    back.print();
    return 0;
}
Read More +
解題區/解題區 - Zerojudge

b425. 抽菸動作 請勿模仿

Problem

替影像打上馬賽克效果，馬賽克效果是把區域內部的每一個像素改變，改變的方法為從舊的圖像中，搜索以其為中心的 11 x 11 方形內部所有像素的平均。

這一題要馬賽克區域為圓形，接著找到方形內部的像素平均即可。

Sample Input

1
2
3
4
0 1 2
1 2
1 2 3
4 5 6

Sample Output

1
2
3
1 2
3 4 5
3 4 5

Solution

馬賽克的方案也許還有其他，從這裡可以知道平均造就出新的像素顏色。那麼可以列出方程式，所謂的解馬賽克估計也是去找線性聯立方程組，利用高斯消去法就能求解，可惜的是會有多組解的情況吧。

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;
class IMAGE {
public:
    struct Pixel {
        double r, g, b, a;
        Pixel(double x = 0, double y = 0, double z = 0, double w = 0):
            r(x), g(y), b(z), a(w) {}
        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, a-x.a);
    	}
    	Pixel operator+(const Pixel &x) const {
        	return Pixel(r+x.r, g+x.g, b+x.b, a+x.a);
    	}
    	Pixel operator*(const double x) const {
        	return Pixel(r*x, g*x, b*x, a*x);
    	}
        Pixel operator/(const double x) const {
        	return Pixel(r/x, g/x, b/x, a/x);
    	}
    	void print() {
    		printf("%d %d %d", (int)round(r), (int)round(g), (int)round(b));
    	}
    };
    int W, H;
    static const int MAXN = 300;
    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 blur(long long X, long long Y, long long R) {
        if (R < 0)	return ;
        for (int i = 0; i < H; i++) {
            for (int j = 0; j < W; j++) {
                if ((i-X)*(i-X) + (j-Y)*(j-Y) <= R*R) {
                    Pixel B(0, 0, 0);
                    int cnt = 0;
                    for (int a = -5; a <= 5; a++) {
                        for (int b = -5; b <= 5; b++) {
                            if (i+a < H && j+b < W && i+a >= 0 && j+b >= 0) {
                                B = B + data[i+a][j+b];
                                cnt++;
                            }
                        }
                    }
                    tmp[i][j] = B / cnt;
                } else {
                    tmp[i][j] = data[i][j];
                }
            }
        }
        
        for (int i = 0; i < H; i++)
            for (int j = 0; j < W; j++)
                data[i][j] = tmp[i][j];
    }
    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' : ' ');
    }
} test;
int main() {
    long long X, Y, R;
    scanf("%lld %lld %lld", &X, &Y, &R);
    test.read();
    test.blur(Y, X, R);
    test.print();
    return 0;
}
Read More +
解題區/解題區 - Zerojudge

b424. 圖片縮放

Problem

此題有精準度誤差，未能提交成功。

進行圖片縮放，有兩種方案來補足缺失或者是簡化的像素計算。不管縮放如何，在新的坐標系的點可以對應到舊坐標系中，並且會有四個像素點 (正方形) 包含這個對應點，有可能會有所缺失，但至少包含一個點。

  • 最近鄰點法 (nearest-neighbor approach)
    找正方形中哪個點最靠近新的對應點，就令其像素顏色等於最靠近的點。
  • 雙線性內插 (bi-linear interpolation)
    顏色由相鄰四個點決定，使用雙性內插來得到新的像素顏色。

特別注意，雖然邊長根據縮放四捨五入，但在計算時先不考慮最終結果的長寬，否則縮放倍率也會跟著四捨五入。

可以參考 中央大學影像處理實驗室 課程 講義第二章

Sample Input

1
2
3
2 2 0
2 1 
1 2 3 255 4 5 6 255

Sample Output

1
2
3
4 2
1 2 3 255 3 4 5 255 4 5 6 255 4 5 6 255
1 2 3 255 3 4 5 255 4 5 6 255 4 5 6 255

Solution

除了陷阱的邊長四捨五入外，精準度根據雙性內插影響。

雙線性內插可以由下列兩種方案，或者還有其他的

  • 兩次 lerp 插值，先水平後垂直或者先垂直後水平
  • 特別的 $f(x, y) = (1-a)(1-b) g(j, k) + a(1-b) g(j+1, k) + (1-a) b g(j, k+1) + a b g(j+1, k+1)$

特別發現到，第二種方案沒有用到除法運算，但乘法運算的優先順序會影響到誤差。然而我的代碼是混用造成難以評估的誤差。關於第二種方案的坐標系，可以參考講義中的圖示。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
#include <bits/stdc++.h> 
using namespace std;
const double eps = 1e-18;
class IMAGE {
public:
    struct Pixel {
        long double r, g, b, a;
        Pixel(long double x = 0, long double y = 0, long double z = 0, long double w = 0):
            r(x), g(y), b(z), a(w) {}
        void read() {
            int p1, p2, p3, p4;
            scanf("%d %d %d %d", &p1, &p2, &p3, &p4);
            r = p1, g = p2, b = p3, a = p4;
        }
        Pixel operator-(const Pixel &x) const {
        	return Pixel(r-x.r, g-x.g, b-x.b, a-x.a);
    	}
    	Pixel operator+(const Pixel &x) const {
        	return Pixel(r+x.r, g+x.g, b+x.b, a+x.a);
    	}
    	Pixel operator*(const long double x) const {
        	return Pixel(r*x, g*x, b*x, a*x);
    	}
        Pixel operator/(const long double x) const {
        	return Pixel(r/x, g/x, b/x, a/x);
    	}
    	void print() {
    		double p1 = r, p2 = g, p3 = b, p4 = a;
    		printf("%.0lf %.0lf %.0lf %.0lf", p1 + eps, p2 + eps, p3 + eps, p4 + eps);
    	}
    };
    int W, H;
    static const int MAXN = 512;
    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();
    }
    int isValid(int x, int y) {
        return x >= 0 && y >= 0 && x < H && y < W;
    }
    void resize(long double SX, long double SY, int trans = 0) {
        if (trans != 0 && trans != 1)
            return ;			
        int rW = (int) round(SX * W);
        int rH = (int) round(SY * H);
        if (rW > 256 || rH > 256 || rW <= 0 || rH <= 0)
            return ;
        for (int i = 0; i < rH; i++) {
            for (int j = 0; j < rW; j++) {
                long double x = i / SY;
                long double y = j / SX;
                int lx, rx, ly, ry;
                lx = floor(x), rx = ceil(x);
                ly = floor(y), ry = ceil(y);
                
                Pixel v[] = {data[lx][ly], data[lx][ry], data[rx][ly], data[rx][ry]};
                int px[] = {lx, lx, rx, rx};
                int py[] = {ly, ry, ly, ry};
                if (trans == 0) {
                    Pixel t1, t2;
                    if (isValid(px[0], py[0]) && isValid(px[1], py[1]) && isValid(px[2], py[2]) && isValid(px[3], py[3])) {
                        long double a = x - lx, b = y - ly;
                        tmp[i][j] = v[0] * (1 - a) * (1 - b) + v[2] * a * (1 - b)
                                    + v[1] * (1 - a) * b + v[3] * a * b;
                    } else if (isValid(px[0], py[0]) && isValid(px[2], py[2])) {
                        if (rx != lx)
                            tmp[i][j] = v[0] + (v[2] - v[0]) * ((x - lx) / (rx - lx));
                        else
                            tmp[i][j] = v[0];
                    } else if (isValid(px[0], py[0]) && isValid(px[1], py[1])) {
                        if (ry != ly)
                            tmp[i][j] = v[0] + (v[1] - v[0]) * ((y - ly) / (ry - ly));
                        else
                            tmp[i][j] = v[0];
                    } else if (isValid(px[0], py[0])) {
                        tmp[i][j] = v[0];
                    } else {
                        puts("WTF");
                    }
                } else {
                    int c = -1;
                    double mndist;
                    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;
                    }
                    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];
    }
    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' : ' ');
    }
} test;
int main() {
    double X, Y;
    int TYPE;
    scanf("%lf %lf %d", &X, &Y, &TYPE);
    test.read();
    test.resize(X, Y, TYPE);
    test.print();
    return 0;
}
Read More +
解題區/解題區 - Zerojudge

b423. 魔術橡皮擦

Problem

提供影像去背、相鄰去色兩種方案。

  • 去背對整張影像找像素差小於某一個範圍內，就將其 alpha 值設成 0。
  • 相鄰去色對鄰近四個方向上下左右進行搜索，將像素差低於某個閥值的都加入到集合中進行 alpha 值設成 0。

Sample Input

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

Sample Output

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

Solution

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
#include <bits/stdc++.h>
using namespace std;
int X, Y, TYPE, DIFF;
namespace Image {
    const int MAXN = 256;
    int n, m, rgba[MAXN][MAXN][4], used[MAXN][MAXN] = {};
    int RGB_diff(int x, int y, int a, int b) {
        return (rgba[x][y][0]-rgba[a][b][0])*(rgba[x][y][0]-rgba[a][b][0])+
                (rgba[x][y][1]-rgba[a][b][1])*(rgba[x][y][1]-rgba[a][b][1])+
                (rgba[x][y][2]-rgba[a][b][2])*(rgba[x][y][2]-rgba[a][b][2]);
    }
    void gTransparent(int x, int y, long long diff) {
        if (x < 0 || y < 0 || x >= m || y >= n || diff < 0)
            return ;
        diff = diff * diff;
        for (int i = 0; i < m; i++)
            for (int j = 0; j < n; j++)
                if (RGB_diff(x, y, i, j) <= diff)
                    rgba[i][j][3] = 0;
    }
    void lTransparent(int x, int y, int ox, int oy, long long diff) {
        if (ox < 0 || oy < 0 || ox >= m || oy >= n || diff < 0)
            return ;
        static const int dx[] = {0, 0, 1, -1};
        static const int dy[] = {1, -1, 0, 0};
        if (x < 0 || y < 0 || x >= m || y >= n)
            return ;
        if (RGB_diff(x, y, ox, oy) > diff * diff || used[x][y])
            return ;
        rgba[x][y][3] = 0, used[x][y] = 1;
        for (int i = 0; i < 4; i++)
            lTransparent(x+dx[i], y+dy[i], ox, oy, diff);
    }
    void print() {
        printf("%d %d\n", n, m);
        for (int i = 0; i < m; i++) {
            for (int j = 0; j < n; j++)
                printf("%d %d %d %d%c", rgba[i][j][0], rgba[i][j][1], rgba[i][j][2], rgba[i][j][3], j == n-1 ? '\n' : ' ');
        }
    }
}
int main() {
    scanf("%d %d %d %d", &X, &Y, &TYPE, &DIFF);
    scanf("%d %d", &Image::n, &Image::m);
    for (int i = 0; i < Image::m; i++)
        for (int j = 0; j < Image::n; j++)
            for (int k = 0; k < 4; k++)
                scanf("%d", &Image::rgba[i][j][k]);
    swap(X, Y);
    if (TYPE == 0)
        Image::gTransparent(X, Y, DIFF);
    if (TYPE == 1)
        Image::lTransparent(X, Y, X, Y, DIFF);
    Image::print();
    return 0;
}
Read More +
解題區/解題區 - Zerojudge

b422. Colorful Life and Monochromatic Life

Problem

將一張 rgb 表示的彩色影像變成灰階影像。

套用公式 round((r + g + b)/3) 得到新的像素色彩。

Sample Input

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

Sample Output

1
2
3
3 2
2 2 2 5 5 5 8 8 8
11 11 11 14 14 14 17 17 17

Solution

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
#include <bits/stdc++.h> 
using namespace std;
int main() {
    int n, m;
    while (scanf("%d %d", &n, &m) == 2) {
        printf("%d %d\n", n, m);
        for (int i = 0; i < m; i++) {
            for (int j = 0; j < n; j++) {
                int r, g, b, w;
                scanf("%d %d %d", &r, &g, &b);
                w = (int) round((r + g + b)/3.0);
                printf("%d %d %d%c", w, w, w, j == n-1 ? '\n' : ' ');
            }
        }
    }
    return 0;
}
Read More +
解題區/解題區 - Zerojudge

b419. 公平的硬幣

Problem

一天，liouzhou_101 去打印店裏打印了幾張複習資料，花了 4 毛錢，他給了店主1張塊錢的鈔票，結果店主補了他兩枚硬幣，一個5毛錢，一個 1 毛錢。

在走回宿舍的途中，liouzhou_101 一直在想一個問題，5 毛錢的硬幣和 1 毛錢的硬幣哪個更公平？所謂公平，就是指將一枚硬幣拋擲 1 次，他正面朝上的概率是 12，反面朝上的概率也是 12。

結果他一會去就把1毛錢的硬幣拋了 1000 次， 記下來有 531 次正面朝上。然後他突然說道：“這硬幣不公平！！！”

真的不公平嗎？他希望你計算一下出現這種情況的概率。

Sample Input

1
2
3
4
5
10 5 5
20 6 8
100 40 60
1000 531 531
5000 2345 2456

Sample Output

1
2
3
4
5
0.2461
0.2310
0.9648
0.0037
0.1093

Solution

跟柳柳州互相出題目，已經來到了第 n 天，儘管我出的都是垃圾跟搬運題目，一出題被電得好愉悅，出在多的垃圾題都電不到柳柳州。

在擲硬幣找區間次數 $[a, b]$ 的機率，由於是離散的統計上就只能推組合數學，但又不是曲棍棒的模型 (Hockey Stick Pattern)，算了先手動窮舉，用 Stirling’s formula 計算組合，後來發現有幾組數據有問題。

當 n 太大的時候，機率分布就相當稀疏，直接套用數值積分，直接去積分自然分布得到公式稍微困難，利用 simpson 或者是 romberg 去解決，但這樣還會有一個問題，只有在平均值的部分機率非常高 (類似離散的 00000100000)，導致自適應辛普森積分會判斷錯誤，提醒剖半積分後測資就正確了不少。似乎合作出了一個可怕的題目，儘管如此，還是得用 mathlab 或者是 wolframalpha 來協助測資的檢驗，剩下的部分就交給柳柳州了。

關於自然分布 (normal distribution)：

  • 骰子擲 $n$ 次，命中的機率 $p = 1/2$
  • 期望值 $\mu = np = n/2$
  • 標準差為 $\sigma = \sqrt{np(1-p)} = \sqrt{n/4}$。

期望值 $E[x] = \sum_{k=1}^{n} k \binom{n}{k} p^k (1-p)^{n-k}$，消階層之後提出，計算後得到 $E[x] = np$。

而標準差 $\sigma$ 先從變異數下手

  • 定義： $Var[x] = E[(x - \mu)^2] = \sum (x - \mu)^2 P[X = x]$
  • 拆解步驟 1： $Var[x] = \sum (x^2 - 2 \mu x + \mu^{2}) P[X = x]$
  • 拆解步驟 2： $Var[x] = E[x^2] - 2 E[x] E[x] + E[x] E[x] = E[x^2] - E[x]^2$
  • 拆解步驟 3： $E[x^2] = E[x(x-1)] + E[x] = n(n-1)p^2 + np$
  • 拆解步驟 4： $Var[x] = np(1-p)$

接下來直接帶入自然分布的公式，參數 $(\mu, \sigma)$，特別小心使用辛普森積分 $[a, b]$ 時，由於問題是離散的，要改變成 $[a - 0.5, b + 0.5]$，防止 $a = b$ 時造成積分錯誤。

特別感謝學弟 firejox 的數學指導。

不久之後，firejox 說道內建函數 erf() 可以解決積分問題，詳見 wiki Gauss error function，因此這一題也可以不寫辛普森積分，套用內建函數在不到 10 行內解決此題。

暴力檢測

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
#include <bits/stdc++.h>
using namespace std;
const long double PI = acos(-1), e = exp(1);
const int MAXN = 100001;
double f[MAXN];
double stirling(int n) {
    if (n < MAXN)
        return f[n];
    return log2(sqrt(2*PI*n)) + n*log2(n/e);
}
double logC(int n, int m) {
    return stirling(n) - stirling(n - m) - stirling(m);
}
int main() {
    f[0] = 0;
    for (int i = 1; i < MAXN; i++)
        f[i] = log2(i) + f[i-1];
    int n, a, b;
    int cases = 0;
    while (scanf("%d %d %d", &n, &a, &b) == 3) {
        double ret = 0;
        for (int i = a; i <= b; i++)
            ret += pow(2, logC(n, i) - log2(2) * n);
        printf("%.4lf\n", ret);
        if (++cases > 100)
            break;
    }
    return 0;
}

數值積分

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
#include <bits/stdc++.h>
using namespace std;
const long double PI = acos(-1), e = exp(1);
const int MAXN = 1001;
double f[MAXN];
double stirling(int n) {
    if (n < MAXN)
        return f[n];
    return log2(sqrt(2*PI*n)) + n*log2(n/e);
}
double logC(int n, int m) {
    return stirling(n) - stirling(n - m) - stirling(m);
}
class SimpsonIntegral {
public:
    const double eps = 1e-6;
    const double pi = acos(-1);
    const double sqr2pi = sqrt(2 * pi);
    double mu, sigma; // function coefficient
    double f(double x) { // usually replace
        return exp(-(pow(x - mu, 2)/2.0/sigma/sigma))/ sigma / sqr2pi;
    }
    void setCoefficient(double m, double s) {
        this->mu = m, this->sigma = s;
    } 
    double simpson(double a, double b, double fa, double fb, double fab) {
        return (fa + 4 * fab + fb) * (b - a) / 6.0;
    }
    double integral(double l, double r) {
        return integral(l, r, eps);
    }
private:
    double integral(double a, double b, double c, double eps, double A, double fa, double fb, double fc) {
        double ab = (a + b)/2, bc = (b + c)/2;
        double fab = f(ab), fbc = f(bc);
        double L = simpson(a, b, fa, fb, fab), R = simpson(b, c, fb, fc, fbc);
        if (fabs(L + R - A) <= 15 * eps)
            return L + R + (L + R - A) / 15.0;
        return integral(a, ab, b, eps/2, L, fa, fab, fb) + integral(b, bc, c, eps/2, R, fb, fbc, fc);
    }
    double integral(double a, double b, double eps) {
        double ab = (a + b) /2;
        double fa = f(a), fb = f(b), fab = f(ab);
        return integral(a, ab, b, eps, simpson(a, b, fa, fb, fab), fa, fab, fb);
    }
} tool;
double binom(double n, double a, double b) {
    double a1, b1;
    double d = sqrt(n/2.0);
    a1 = 0.5 * (1.0 + erf((a - n/2.0 - 0.5) / d));
    b1 = 0.5 * (1.0 + erf((b - n/2.0 + 0.5) / d));
    return b1 - a1;
}
int main() {
    f[0] = 0;
    for (int i = 1; i < MAXN; i++)
        f[i] = log2(i) + f[i-1];
    int n, a, b;
    int cases = 0;
    while (scanf("%d %d %d", &n, &a, &b) == 3) {
        cases++;
        if (n < MAXN) {
            double ret = 0;
            for (int i = a; i <= b; i++)
                ret += pow(2, logC(n, i) - log2(2) * n);
            printf("%.4lf\n", ret);
        } else {
//			tool.setCoefficient(n/2.0, sqrt(n/4.0));
//			double ret = 0, l = a - 0.5, r = b + 0.5;
//			if (r < n/2.0 || l > n/2.0)
//				ret = tool.integral(l, r);
//			else
//				ret = tool.integral(l, n/2.0) + tool.integral(n/2.0, r);
//			printf("%.4lf\n", ret);
            printf("%.4lf\n", binom(n, a, b));
        } 
    }
    return 0;
}
Read More +
解題區/解題區 - Zerojudge

b417. 區間眾數

Problem

背景

大學上課是怎麼回事的呢？談論到學期成績計算，很多奇聞軼事可以說的，例如教授任選三題批改、同學任選三題作答、滿分三百、沒有考試、 … 等。筆試類型的考試只是最常見的一種，還有所謂的多才華加分，例如上台報告、舉手發問、參與展覽、參加競賽、 … 等。

不幸地，某 M 修了一堂怪課，教授每一堂課結束後總是會大聲宣揚「葛萊分多加 5 分」對於需要不斷複習演練的某 M 而言，這門課的即時評分方式讓其非常挫敗。神奇的是，有一次教授這麼問同學「給一序列，請計算眾數、平均數、中位數，答對加分。」某 M 剎那間傻住，大學的學期成績可以用這種問題獲得，當他在懷疑這個問題時，問題早已被同學回答走。

某 M 非常不甘心，計算眾數、平均數、中位數就能夠加分，要求只是 n = 10 的序列。既然加分有理，出題吧！

題目描述

平均數、中位數可以藉由區間總和、區間 K 小問題來完成，現在來解決眾數。

給定一個長度為 $n$ ( $1 \le n \le 100000$ ) 的正整數序列 $s$ ($1 \le s_i \le n$)，對於 $m$ ( $1 \le m \le 1000000$) 次詢問 $l, r$，每次輸出 $[s_l, \text{... },s_r]$ 中，出現最多次的次數以及有多少種數字出現最多次。

Sample Input

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

Sample Output

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

Solution

眾數區間查找 Range mode query 性質可以參考維基百科。

在這一題中，這裡提供三種解法，但我的寫法中只有莫隊算法可以通過，其餘算法受到時間跟空間上的限制，需要額外的優化，細節就不多做說明：

  • 莫隊算法
    只能離線處理，無法強制在線詢問，維護眾數最高頻率指針，來統計眾數的對應次數和等價眾數的個數，由於指針移動都是 $O(1)$，整體上時間效能為 $O(N^{1.5})$。
  • 分塊在線 1
    離線預處理 $O(N^{0.5} \times N^{0.5} \times \log{N})$
    提供在線詢問 $O(\text{ans} \times \log{N})$，由於知道每一個塊中的代表眾數都可能是答案，因此預處理每一個塊的眾數，其餘的非眾數捨棄掉，最慘情況是每個數字都出現一次，在線詢問時複雜度會掉到 $O(N \log N)$。
  • 分塊在線 2
    • 離線預處理時間 $O(N \times L^2)$，提供在線詢問 $O(N/L)$，其中 $L$ 表示有多少個分塊。
    • 預處理分塊編號 $PILE[i, j]$ 的眾數情況 (類似莫隊維護指針)，因此記憶體消耗 $O(N \times L^2)$。
    • 對於每一處詢問，會對應預處理中的區域塊的紀錄 $PILE[L, R]$，接著頭尾兩端的殘存數字再依序加入，意即 $A[a, L-1]$ 和 $A[R+1, b]$ 的數字加入的複雜度 $O(1)$，隨後再進行撤銷 $O(1)$。
    • 假設詢問次數 $Q$ 與 $N$ 呈現性關係，整體需要 $O(N \times L^2 + N^2 / L)$，則 $L = N^{1/3}$ 是最好的情況，意即每一個塊具有 $O(N^{2/3})$ 個數字，最後得到整體時間效能為 $O(N^{1.66})$。

由於這一題目標不是找到眾數為何 (這牽涉到最小字典順序)，而是找眾數的出現個數和有幾種眾數可能，莫隊算法會快於其他算法。若是要求最小眾數莫隊算法需要 $O(N^{1.5} \log{N})$，分塊在線 1 也是 $O(N^{1.5} \log{N})$。

分塊在線 2 在此題不僅僅遇到記憶體過大，只好鬆弛 $L$ 來壓低記憶體用量，但增加時間需求，時間上慢個 10 倍以上。

莫隊算法

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
#include <bits/stdc++.h> 
using namespace std;
const int MAXN = 100005;
const int MAXV = 100005;
const int MAXQ = 1000005;
class Offline {
public:
    struct Event {
        static int PILE, belong[MAXV];
        int l, r, qid;
        Event(int a = 0, int b = 0, int e = 0):
            l(a), r(b), qid(e) {}
        bool operator<(const Event &x) const {
            if (belong[l] != belong[x.l])
                return l < x.l;
            return r < x.r; 
        }
    } event[MAXQ];
    int A[MAXN], N, qsize;
    int ret[MAXQ][2];
    void init(int N) {
        this->N = N, qsize = 0;
        memset(val_count, 0, sizeof(val_count));
        memset(mode_freq, 0, sizeof(mode_freq));
        mode_idx = 0;
        Event::PILE = sqrt(N);
        for (int i = N; i >= 1; i--)
            Event::belong[i] = i / Event::PILE;
    }
    void q_add(int l, int r) {
        event[qsize] = Event(l, r, qsize);
        qsize++;
    }
    void run() {
        sort(event, event+qsize);
        
        int l = 1, r = 0;
        for (int i = 0; i < qsize; i++) {
            while (r < event[i].r)	r++, update(A[r], 1);
            while (r > event[i].r)	update(A[r], -1), r--;
            while (l > event[i].l)	l--, update(A[l], 1);
            while (l < event[i].l)	update(A[l], -1), l++;
            ret[event[i].qid][0] = mode_idx;
            ret[event[i].qid][1] = mode_freq[mode_idx];
        }
    }
private:
    // unrolled
    int val_count[MAXV], mode_freq[MAXN], mode_idx;
    inline void update(int x, int val) {
        if (val == 1) {
            mode_freq[val_count[x]]--;
            val_count[x]++;
            mode_freq[val_count[x]]++;
            if (mode_freq[mode_idx+1] != 0)
                mode_idx++;
        } else if (val == -1) {
            mode_freq[val_count[x]]--;
            val_count[x]--;
            mode_freq[val_count[x]]++;
            if (mode_freq[mode_idx] == 0)
                mode_idx--;
        }
    }
} off;
int Offline::Event::PILE = 16, Offline::Event::belong[MAXV];
namespace mLocalStream {
    class M {
    public:
        static const int N = 1048576;
        char buf[N], *p, *end;
        M() {
            p = buf, end = buf + N - 1;
        }
        void printInt(int x) {
            static char stk[16];
            int idx = 0;
            stk[idx++] = '\n';
            if (!x)	
                stk[idx++] = '0';
            while (x)
                stk[idx++] = x%10 + '0', x /= 10;
            while (idx) {
                if (p == end) {
                    *p = '\0';
                    printf("%s", buf), p = buf;
                }
                *p = stk[--idx], p++;
            }
        }	
        static inline int readchar() {
            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++;
        }
        static 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;
        }
        ~M() {
            *p = '\0';
            printf("%s", buf);
        }		
    } bprint;
}
int main() {
    int n, m, l, r, a, b;
    while (mLocalStream::M::ReadInt(&n)) {
        mLocalStream::M::ReadInt(&m);
        off.init(n);
        for (int i = 1; i <= n; i++)
            mLocalStream::M::ReadInt(&off.A[i]);
        for (int i = 0; i < m; i++) {
            mLocalStream::M::ReadInt(&l);
            mLocalStream::M::ReadInt(&r);
            off.q_add(l, r);
        }
        
        off.run();
        for (int i = 0; i < m; i++)
            printf("%d %d\n", off.ret[i][0], off.ret[i][1]);
    }
    return 0;
}

分塊在線 1

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
#include <bits/stdc++.h> 
using namespace std;
const int MAXN = 100005;
class ModeQuery {
public:
    int A[MAXN], N, PILE, belong[MAXN];
    vector< vector<int> > V, Vmode;
    vector<int> POS[MAXN];
    int q_dup[MAXN], q_cases;
    void make(int n) {
        N = n, q_cases = 0;
        for (PILE = 1; PILE * PILE < N; PILE <<= 1);
        
        V.clear(), Vmode.clear();
        for (int l = 0, r; l < N; l = r+1) {
            r = min(l + PILE - 1, N-1);
            V.push_back(vector<int>(A+l, A+r+1));
            Vmode.push_back(vector<int>());
        }
        
        for (int i = 0; i < V.size(); i++) {
            sort(V[i].begin(), V[i].end());
            int mx_cnt = 1, cnt = 1;
            for (int j = 1; j <= V[i].size(); j++) {
                if (j == V[i].size() || V[i][j] != V[i][j-1]) {
                    if (cnt > mx_cnt)	
                        mx_cnt = cnt, Vmode[i].clear();
                    if (cnt == mx_cnt)
                        Vmode[i].push_back(V[i][j-1]);
                    cnt = 1;
                } else {
                    cnt++;
                }
            }
        }
        
        for (int i = 1; i <= N; i++) // value range
            POS[i].clear(), q_dup[i] = 0;
        for (int i = 0; i < N; i++)
            POS[A[i]].push_back(i);
        for (int i = N-1; i >= 0; i--)
            belong[i] = i / PILE;
    }
    int compute_cnt(int x, int l, int r) {
        if (q_dup[x] == q_cases)
            return 0;
        q_dup[x] = q_cases;
        int va = (int) (lower_bound(POS[x].begin(), POS[x].end(), l) - POS[x].begin()) - 1;
        int vb = (int) (upper_bound(POS[x].begin(), POS[x].end(), r) - POS[x].begin()) - 1;
        return vb - va;
    }
    void query(int l, int r, int &a, int &b) {
        static int t;
        q_cases++;
        a = b = 0;
        if (belong[l] == belong[r]) {
            for (int i = l; i <= r; i++) {
                t = compute_cnt(A[i], l, r);
                if (t > a)	a = t, b = 0;
                if (t == a)	b++;
            }	
            return ;		
        }
        for (int i = belong[l]+1; i < belong[r]; i++) {
            for (int j = 0; j < Vmode[i].size(); j++) {
                t = compute_cnt(Vmode[i][j], l, r);
                if (t > a)	a = t, b = 0;
                if (t == a)	b++;
            }
        }
        for (int i = (belong[a]+1)*PILE-1; i >= l; i--) {
            t = compute_cnt(A[i], l, r);
            if (t > a)	a = t, b = 0;
            if (t == a)	b++;
        }
        for (int i = (belong[b])*PILE; i <= r; i++) {
            t = compute_cnt(A[i], l, r);
            if (t > a)	a = t, b = 0;
            if (t == a)	b++;
        }
    }
} dream;
int main() {
    int n, m, l, r, a, b;
    while (scanf("%d %d", &n, &m) == 2) {
        for (int i = 0; i < n; i++)
            scanf("%d", &dream.A[i]);
        dream.make(n);
        for (int i = 0; i < m; i++) {
            scanf("%d %d", &l, &r);
            l--, r--;
            dream.query(l, r, a, b);
            printf("%d %d\n", a, b);
        }
    }
    return 0;
}

分塊在線 2

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
#include <bits/stdc++.h> 
using namespace std;
const int MAXN = 100005;
const int MAXS = 35;
class ModeQuery {
public:
    struct Kit {
        unsigned short cnt[MAXN], freq[MAXN];
        int mode_idx;
        Kit() {
            init();
        }
        void init() {
            mode_idx = 0;
        }
        void add(int x, int val) {
            freq[cnt[x]]--;
            cnt[x] += val;
            freq[cnt[x]]++;
            if (cnt[x] > mode_idx)
                mode_idx = cnt[x];
        }
        void resume(int x, int val) {
            freq[cnt[x]]--;
            cnt[x] -= val;
            freq[cnt[x]]++;
        } 
    } f[MAXS * MAXS/2], kit;
    int fr[MAXS][MAXS];
    int A[MAXN], N, PILE, belong[MAXN];
    vector< pair<int, int> > V;
    void heavy() {
        int ff = 0;
        for (int i = 0; i < V.size(); i++)
            for (int j = i; j < V.size(); j++)
                fr[i][j] = ff++;
        for (int i = 0; i < V.size(); i++)
            for (int j = 0; j < i; j++)
                fr[i][j] = ff;
        for (int i = 0; i < V.size(); i++) {
            f[fr[i][i]].init();
            for (int j = V[i].first; j <= V[i].second; j++)
                f[fr[i][i]].add(A[j], 1);
        }
        for (int i = 0; i < V.size(); i++) {
            for (int j = i+1; j < V.size(); j++) {
                f[fr[i][j]] = f[fr[i][j-1]];
                for (int k = V[j].first; k <= V[j].second; k++)
                    f[fr[i][j]].add(A[k], 1);
            }
        }
    }
    void make(int n) {
        N = n;
        PILE = (int) pow(n, 0.70) + 1;
                
        V.clear();
        for (int l = 0, r; l < N; l = r+1) {
            r = min(l + PILE - 1, N-1);
            V.push_back(make_pair(l, r));
        }
        for (int i = N-1; i >= 0; i--)
            belong[i] = i / PILE;
        heavy();
    }
    void query(int l, int r, int &a, int &b) {
        int bl = belong[l], br = belong[r];
        if (bl == br || r - l <= PILE + 10) {
            kit.init();
            for (int i = l; i <= r; i++)
                kit.add(A[i], 1);
            a = kit.freq[kit.mode_idx], b = kit.mode_idx;
            for (int i = l; i <= r; i++)
                kit.resume(A[i], 1);
            return ;
        }
        if (l >= 0 && belong[l-1] == bl)	
            bl++;
        if (r < N && belong[r+1] == br)
            br--;
        Kit &tmp = f[fr[bl][br]];
        int tmp_mode_idx = tmp.mode_idx;
        for (int i = bl*PILE - 1; i >= l; i--)
            tmp.add(A[i], 1);
        for (int i = (br+1)*PILE; i <= r; i++)
            tmp.add(A[i], 1);
        a = tmp.freq[tmp.mode_idx], b = tmp.mode_idx;
        for (int i = bl*PILE - 1; i >= l; i--)
            tmp.resume(A[i], 1);
        for (int i = (br+1)*PILE; i <= r; i++)
            tmp.resume(A[i], 1);
        tmp.mode_idx = tmp_mode_idx;
    }
} dream;
int main() {
    int n, m, l, r, a, b;
    while (scanf("%d %d", &n, &m) == 2) {
        for (int i = 0; i < n; i++)
            scanf("%d", &dream.A[i]);
        dream.make(n);
        for (int i = 0; i < m; i++) {
            scanf("%d %d", &l, &r);
            l--, r--;
            dream.query(l, r, a, b);
            printf("%d %d\n", b, a);
        }
    }
    return 0;
}
Read More +
解題區/解題區 - Zerojudge

b418. 八成真物

Problem

「偽物比真物更有價值」—貝木泥舟
「真物和偽物一樣有價值」—忍野咩咩
「偽娘比真娘更有價值」—精英島民

任兩個物品的相似度$sim(A, B) = \frac{|A \cap B|}{|A \cup B|}$，換句話說把 $A$ 具有的特徵和 $B$ 具有的特徵類型取交集、聯集個數，相除就能得到其相似度。例如有 5 個特徵，若 A 表示成 11001、B 表示成 01100$sim(A, B) = \frac{|{2}|}{|{1, 2, 3, 5}|} = 0.25$。

現在盤面上有 N 個物品、M 種特徵，請問相似度大於等於 0.8 的相似對數 $S$ 有多少種。為了讓這一題更有趣味，算法允許偽物，輸出$\frac{S}{N(N-1)/2} \times 100$。

Sample Input

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
8 10
0111101101
0111101100
0111101100
0111101101
1100101010
0111101100
0111101100
0011101111
  
8 10
1111101100
1011101100
1010101100
1010111010
0100110100
1011101100
1110110000
1011111100

Sample Output

1
2
53.57
25.00

Solution

利用 std::bitset 加速交集和聯集運算，將數個屬性壓縮到一個 32-bits 單位，然後藉 bitcount 盡可能在 $O(1)$ 的時間內完成，就有機會通過這一題。一般的 $O(N^2 M)$ 比 $O(N^2 M/32)$ 慢個五六倍。

此外，特別小心浮點數計算誤差，5.0/4.0 >= 0.8 == false，用乘法判定大於閥值的操作。

關於 LSH 的部分，還在進行測試，若 signature 計算太慢，還不如直接暴力法，若能分散找到 signature 或者是預先有辦法處理好，那分桶計算就會好一點。

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
#include <bits/stdc++.h> 
using namespace std;
int main() {
    int n, m;
    char s[1024];
    while (scanf("%d %d", &n, &m) == 2) {
        bitset<512> attr[n];
        for (int i = 0; i < n; i++) {
            scanf("%s", s);
            attr[i] = bitset<512>(s);
        }
        
        double ret = 0;
        for (int i = 0; i < n; i++) {
            for (int j = i+1; j < n; j++) {
                int a = (attr[i]&attr[j]).count();
                int b = (attr[i]|attr[j]).count();
                if (5 * a >= 4 * b)
                    ret++;
            }
        }
        printf("%.2lf\n", ret * 100 / (n * (n-1)/2));
    }
    return 0;
}
Read More +