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

[Tmt514] Beverage Cup 2 E - Kirino in Google Code Jam

contents

  1. 1. Problem
  2. 2. Solution
    1. 2.1. atan2 誤差測試
    2. 2.2. 誤差數據產生
    3. 2.3. 對照組
    4. 2.4. Final

Problem

2015 Google Code Jam Round 1A 中，C - Logging 有點嚴重的問題！很多精準度不高的代碼都通過了！現在要求去找到測資 tricky cases，讓他們通通 WA！

原題是對於任意平面點，找到要移除最小的平面點，使得自己在凸包邊緣上。標準的極角排序，維護 180 度以內 (半平面) 的點個數。

atan2 esp = 1e-12 not enough for $x, y \in \left [0, 10^6 \right]$，極角排序替代方案！快來戳我的講法，看 POJ 類似的題目，精準度要壓到 1e-16 呢，所以誤差還要抓更小才行。

Solution

等等，challenge 這一題時，自己的代碼也 WA 了！

一開始誤以為是 index out of bound，assert() 下去都沒發生。後來估計是 atan2(y, x) 的誤差，對此誤差早有耳聞，解題撞壁的機會很低。咱很蠢，隨機生了 1000 萬個平面點，按照 atan2 排序，檢查相鄰兩個座標的外積 (叉積) 是否為逆時針，跑一次隨機生成大概能爆出一個誤差的相鄰點，十來次得到 n 個點，隨機組合這 n 個點於四個象限，再亂撒少數個點，再跟自己撰寫的代碼做比對 (中間又測了上萬組)。

atan2 誤差測試

atan2_eps_test
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 <stdio.h>
#include <math.h>
#include <vector>
#include <assert.h>
#include <time.h>
#include <algorithm>
using namespace std;
#define eps 1e-8
struct Pt {
    double x, y;
    Pt(double a = 0, double b = 0):
    	x(a), y(b) {}	
    Pt operator-(const Pt &a) const {
        return Pt(x - a.x, y - a.y);
    }
    Pt operator+(const Pt &a) const {
        return Pt(x + a.x, y + a.y);
    }
    Pt operator*(const double a) const {
        return Pt(x * a, y * a);
    }
    bool operator==(const Pt &a) const {
    	return fabs(x - a.x) < eps && fabs(y - a.y) < eps;
    }
    bool operator<(const Pt &a) const {
        if (fabs(x - a.x) > eps)
            return x < a.x;
        if (fabs(y - a.y) > eps)
            return y < a.y;
        return false;
    }
    double length() {
        return hypot(x, y);
    }
    void read() {
        scanf("%lf %lf", &x, &y);
    }
};
double dot(Pt a, Pt b) {
    return a.x * b.x + a.y * b.y;
}
double cross(Pt o, Pt a, Pt b) {
    return (a.x-o.x)*(b.y-o.y)-(a.y-o.y)*(b.x-o.x);
}
double cross2(Pt a, Pt b) {
    return a.x * b.y - a.y * b.x;
}
int between(Pt a, Pt b, Pt c) {
    return dot(c - a, b - a) >= -eps && dot(c - b, a - b) >= -eps;
}
int onSeg(Pt a, Pt b, Pt c) {
    return between(a, b, c) && fabs(cross(a, b, c)) < eps;
}
struct Seg {
    Pt s, e;
    double angle;
    int label;
    Seg(Pt a = Pt(), Pt b = Pt(), int l=0):s(a), e(b), label(l) {
//		angle = atan2(e.y - s.y, e.x - s.x);
    }
    bool operator<(const Seg &other) const {
        if (fabs(angle - other.angle) > eps)
            return angle > other.angle;
        if (cross(other.s, other.e, s) > -eps)
            return true;
        return false;
    }
    bool operator!=(const Seg &other) const {
        return !((s == other.s && e == other.e) || (e == other.s && s == other.e));
    }
};
Pt getIntersect(Seg a, Seg b) {
    Pt u = a.s - b.s;
    double t = cross2(b.e - b.s, u)/cross2(a.e - a.s, b.e - b.s);
    return a.s + (a.e - a.s) * t;
}
double getAngle(Pt va, Pt vb) { // segment, not vector
    return acos(dot(va, vb) / va.length() / vb.length());
}	
double distSeg2Point(Pt a, Pt b, Pt p) {
    Pt c, vab;
    vab = a - b;
    if (between(a, b, p)) {
        c = getIntersect(Seg(a, b), Seg(p, Pt(p.x+vab.y, p.y-vab.x)));
        return (p - c).length();
    }
    return min((p - a).length(), (p - b).length());
}
Pt rotateRadian(Pt a, double radian) {
    double x, y;
    x = a.x * cos(radian) - a.y * sin(radian);
    y = a.x * sin(radian) + a.y * cos(radian);
    return Pt(x, y);
}
int monotone(int n, Pt p[], Pt ch[]) {
    sort(p, p+n);
    int i, m = 0, t;
    for(i = 0; i < n; i++) {
        while(m >= 2 && cross(ch[m-2], ch[m-1], p[i]) <= 0)
            m--;
        ch[m++] = p[i];
    }
    for(i = n-1, t = m+1; i >= 0; i--) {
        while(m >= t && cross(ch[m-2], ch[m-1], p[i]) <= 0)
            m--;
        ch[m++] = p[i];
    }
    return m-1;
}
const double pi = acos(-1);
int main() {
    srand ( time(NULL));
    
    vector< pair<double, Pt> > V;
    for (int i = 0; i < 10000000; i++) {
        int X, Y;
        X = (rand()*rand())%2000000 - 1000000;
        Y = (rand()*rand())%2000000 - 1000000;
        V.push_back(make_pair(atan2(Y, X), Pt(X, Y)));
    }
    sort(V.begin(), V.end());
    for (int i = 0; i+1 < V.size(); i++) {
        if (cross2(V[i].second, V[i+1].second) != 0 && fabs(V[i].first - V[i+1].first) < 1e-12)
            printf("%.0lf %.0lf %.0lf %.0lf\n", V[i].second.x, V[i].second.y, V[i+1].second.x, V[i+1].second.y);
    }
//	long long a = 599430;
//	long long b = 1428722;
//	long long c = 648824;
//	long long d = 1546451;
//	long long a = 885828;
//	long long b = 737167;
//	long long c = 898961;
//	long long d = 748096;
//
//	printf("%lld\n", a * d - b * c);
//	printf("%lf %lf, %d\n", atan2(b, a), atan2(d, c), fabs(atan2(b, a) - atan2(d, c)) < 1e-12);
    return 0;
}
/*
885828 737167 898961 748096
854425 732894 706721 606199
-971645 988877 -838743 853618
993753 -652232 991850 -650983
620105 659001 916399 973880
*/

誤差數據產生

testdata
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
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <set>
#include <algorithm>
using namespace std;
double random() {
    int r = rand();
    return (double) r / RAND_MAX;
}
int vvv[20] = {885828, 737167, 898961, 748096,
854425, 732894, 706721, 606199,
-971645, 988877, -838743, 853618,
993753, -652232, 991850, -650983,
620105, 659001, 916399, 973880};
main() {
    int n, m, t, a, b, c, tmp;
    int z, i, j, k, l, p, q;
    srand ( time(NULL));
    freopen("in.txt", "w+t", stdout);
    
    int T = 10000;
    printf("%d\n", T);
    while(T--) {
        int n = 30, m = rand()%n + 1;
        printf("%d\n", n);
        set< pair<int, int> > S;
        S.insert(make_pair(0, 0));
        printf("%d %d\n", 0, 0);
        for (int i = 1; i < n; i++) {
            int x, y;
            do {
                if (rand()%10 < 8) {
                    x = vvv[rand()%20], y = vvv[rand()%20];
                    if (rand()%2)	x = -x;
                    if (rand()%2)	y = -y;
                } else {
                    x = rand()%100 - 50, y = rand()%100 - 50;
                }
            } while (S.count(make_pair(x, y)));
            S.insert(make_pair(x, y));
            printf("%d %d\n", x, y);
        }		
    }
    
    return 0;
}

對照組

利用外積的計算方案，誤差情況會更小。又因為是整數座標，基本上是不產誤差。

fixed-bug-Problem-C-Logging
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
#include <stdio.h>
#include <math.h>
#include <algorithm>
#include <set>
#include <vector>
using namespace std;
#define eps 1e-6
#define MAXN 131072
struct Pt {
    double x, y;
    int label;
    Pt(double a = 0, double b = 0, int c = 0):
    x(a), y(b), label(c) {}
    Pt operator-(const Pt &a) const {
        return Pt(x - a.x, y - a.y);
    }
    Pt operator+(const Pt &a) const {
        return Pt(x + a.x, y + a.y);
    }
    Pt operator*(const double a) const {
        return Pt(x * a, y * a);
    }
    bool operator<(const Pt &a) const {
        if (fabs(x - a.x) > eps)
            return x < a.x;
        if (fabs(y - a.y) > eps)
            return y < a.y;
        return false;
    }
};
double dot(Pt a, Pt b) {
    return a.x * b.x + a.y * b.y;
}
double cross(Pt o, Pt a, Pt b) {
    return (a.x-o.x)*(b.y-o.y)-(a.y-o.y)*(b.x-o.x);
}
double cross2(Pt a, Pt b) {
    return a.x * b.y - a.y * b.x;
}
int between(Pt a, Pt b, Pt c) {
    return dot(c - a, b - a) >= -eps && dot(c - b, a - b) >= -eps;
}
int onSeg(Pt a, Pt b, Pt c) {
    return between(a, b, c) && fabs(cross(a, b, c)) < eps;
}
Pt D[4096];
bool cmp(const Pt& p1, const Pt& p2)
{
    if (p1.y == 0 && p2.y == 0 && p1.x * p2.x <= 0) return p1.x > p2.x;
    if (p1.y == 0 && p1.x >= 0 && p2.y != 0) return true;
    if (p2.y == 0 && p2.x >= 0 && p1.y != 0) return false;
    if (p1.y * p2.y < 0) return p1.y > p2.y;
    double c = cross2(p1, p2);
    return c > 0 || c == 0 && fabs(p1.x) < fabs(p2.x);
}
int main() {
    int N, testcase, cases = 0;
    double x, y;
    scanf("%d", &testcase);
    while (testcase--) {
        scanf("%d", &N);
        for (int i = 0; i < N; i++) {
            scanf("%lf %lf", &x, &y);
            D[i] = Pt(x, y);
        }
        printf("Case #%d:\n", ++cases);
        
        if (N == 1) {
            puts("0");
            continue;
        }
        
        for (int i = 0; i < N; i++) {
            vector< Pt > A;
            for (int j = 0; j < N; j++) {
                if (i == j)
                    continue;
                Pt p = D[j] - D[i];
                A.push_back(p);
            }
            sort(A.begin(), A.end(), cmp);
            int M = (int)A.size();
            int l = 0, r = 0, cnt = 1;
            int ret = 0;
            for (l = 0; l < M; l++) {
                if (l == r)
                    r = (r+1)%M, cnt++;
                while (l != r && cross2(A[l], A[r]) >= 0) {
                    r = (r+1)%M, cnt++;
                }
                ret = max(ret, cnt);
                cnt--;
            }
            printf("%d\n", N - ret);
        }
    }
    return 0;
}

Final

生了 10000 組，終於在裡面發現其中幾組。

final
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
#include <stdio.h>
int main() {
    puts("1\n30\n0 0\n-971645 838743\n748096 -988877\n-652232 -993753\n737167 -838743\n-48 27\n706721 -885828\n606199 854425\n659001 -993753\n898961 885828\n-659001 885828\n748096 -973880\n21 -13\n-748096 606199\n-732894 991850\n13 -12\n659001 -737167\n-32 -32\n737167 -748096\n650983 -971645\n650983 -732894\n854425 -606199\n-606199 885828\n916399 -988877\n652232 -838743\n-606199 988877\n-620105 -652232\n-748096 -737167\n24 -23\n916399 854425\n");
    return 0;
}
/*
1
30
0 0
-971645 838743
748096 -988877
-652232 -993753
737167 -838743
-48 27
706721 -885828
606199 854425
659001 -993753
898961 885828
-659001 885828
748096 -973880
21 -13
-748096 606199
-732894 991850
13 -12
659001 -737167
-32 -32
737167 -748096
650983 -971645
650983 -732894
854425 -606199
-606199 885828
916399 -988877
652232 -838743
-606199 988877
-620105 -652232
-748096 -737167
24 -23
916399 854425
*/