求点集中最近点对--期望线性法

• 参考这里最下面
• n是1200左右时，和暴力解法耗时差不多
• n越大，期望线性法越高效
• n=99999时，暴力需要58秒，期望线性法需要小于1秒
• n=307200时，暴力需要548秒，期望线性法需要5秒左右
• CPU：Intel® Core™ i7-9700 CPU @ 3.00GHz
• 期望线性法的最坏情况，即每次都重新划分，在n=9999时需要39秒
• 20221123：加了分治，需要注意的是函数rec中的pt数组，n=99999时，如果使用vector，需要4.94秒，如果使用普通数组，需要0.07秒。因为vector的数据在堆上，普通数组的数据在栈上，访问栈比堆快得多，而且普通数组数据连续，cache命中率高。同时因为普通数组的数据在栈上导致数组太大会段错误，所以将普通数组放在静态区，访问静态区同样比堆快得多。
• 20221123新版本：增加了保存所有最近点对

#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#define EPS 1e-3

using namespace std;

struct lwd {
    int i1;
    int i2;
};
struct cmp {
    bool operator() (const lwd &l1, const lwd &l2) {
        if(l1.i1 == l2.i1) return l1.i2 < l2.i2;
        return l1.i1 < l2.i1;
    }
};
vector<lwd> res,  res2;
set<lwd, cmp> res1;

void simple(vector<int> p1, vector<int> p2, int &a, int &b)
{
    int n = p1.size(), index = 0;
    double p10 = p1[0], p11 = p1[1], p20 = p2[0], p21 = p2[1];
    double dis = (p10 - p11) * (p10 - p11) + (p20 - p21) * (p20 - p21);
    dis = sqrt(dis);
    a = 0, b = 1;
    for (int i = 0; i < n; ++i)
    {
        for (int j = i + 1; j < n; ++j)
        {
            double p1i = p1[i], p1j = p1[j], p2i = p2[i], p2j = p2[j];
            double tmp = (p1i - p1j) * (p1i - p1j) + (p2i - p2j) * (p2i - p2j);
            tmp = sqrt(tmp);
            if(tmp > dis) continue;
            if (tmp < dis)
            {
                dis = tmp;
                a = i;
                b = j;
                res.clear();
            }
            lwd lwl;
            lwl.i1 = i;
            lwl.i2 = j;
            res.emplace_back(lwl);
        }
    }
}

/***************************************  分治  *******************************************/

struct pt
{
    int x, y, id;
};

struct cmpx
{
    bool operator()(const pt &a, const pt &b) const
    {
        return a.x < b.x || (a.x == b.x && a.y < b.y);
    }
} cmp_x;

struct cmpy
{
    bool operator()(const pt &a, const pt &b) const { return a.y < b.y; }
} cmp_y;

// vector vp;
double mindist;
int ansa, ansb;

inline void upd_ans(const pt &a, const pt &b)
{
    double ax = a.x, bx = b.x, ay = a.y, by = b.y;
    double dist = sqrt((ax - bx) * (ax - bx) + (ay - by) * (ay - by));
    if(dist > mindist) return;
    if (dist < mindist){
        mindist = dist, ansa = a.id, ansb = b.id;
        res1.clear();
    }
    lwd lwl;
    lwl.i1 = a.id;
    lwl.i2 = b.id;
    if(a.id > b.id) {
        lwl.i2 = a.id;
        lwl.i1 = b.id;
    }
    res1.insert(lwl);
}

void rec(vector<pt> &vp, int l, int r)
{
    if (r - l <= 3)
    {
        for (int i = l; i <= r; ++i)
            for (int j = i + 1; j <= r; ++j)
                upd_ans(vp[i], vp[j]);
        sort(vp.begin() + l, vp.begin() + r + 1, cmp_y);
        return;
    }

    int m = (l + r) >> 1;
    int midx = vp[m].x;
    rec(vp, l, m), rec(vp, m + 1, r);
    inplace_merge(vp.begin() + l, vp.begin() + m + 1, vp.begin() + r + 1, cmp_y);

    // vector t(vp.size());
    static pt t[100000];
    int tsz = 0;
    for (int i = l; i <= r; ++i)
        if (abs(vp[i].x - midx) <= mindist)
        {
            t[tsz++] = vp[i];
        }
    for(int i = 0; i < tsz; ++i) {
        for(int j = i+1; j < tsz && t[i].y-t[j].y <= mindist; ++j) {
            upd_ans(t[i], t[j]);
        }
    }
}

/***************************************  期望线性法  *******************************************/
// p1,p2是坐标数组，a,b是最近点对的索引
void linear(vector<int> p1, vector<int> p2, int &a, int &b)
{
    int n = p1.size(), mi1 = p1[0], ma1 = p1[0], mi2 = p2[0], ma2 = p2[0];
    // 求点的范围
    for (int i = 1; i < n; ++i)
    {
        if (mi1 > p1[i])
            mi1 = p1[i];
        if (ma1 < p1[i])
            ma1 = p1[i];
        if (mi2 > p2[i])
            mi2 = p2[i];
        if (ma2 < p2[i])
            ma2 = p2[i];
    }
    int d1 = ma1 - mi1, d2 = ma2 - mi2;
    // 以前两点的距离划分网格
    a = 0, b = 1;
    double p10 = p1[0], p11 = p1[1], p20 = p2[0], p21 = p2[1];
    double s = (p10 - p11) * (p10 - p11) + (p20 - p21) * (p20 - p21);
    s = sqrt(s);
    double ss = s + EPS;
    int s1 = d1 / ss + 1, s2 = d2 / ss + 1;
    // 将前两点放入对应网格
    unordered_map<int, vector<int>> um;
    for (int i = 0; i < 2; ++i)
    {
        int x = (p1[i] - mi1) / ss, y = (p2[i] - mi2) / ss;
        um[y * s1 + x].emplace_back(i);
    }
    lwd lwl;
    lwl.i1 = 0, lwl.i2 = 1;
    res2.emplace_back(lwl);
    // 依次处理剩余点：计算该点与周围九个网格点的最小距离，小于当前划分值则重新划分
    for (int i = 2; i < n; ++i)
    {
        // 当前点所在网格
        int x = (p1[i] - mi1) / ss, y = (p2[i] - mi2) / ss;
        // 周围九个网格与当前点的最小距离
        double tmp = s;
        // 周围九个网格与当前点最近的点的索引
        int aa = -1;
        // 周围九个网格
        for (int ii = y - 1; ii <= y + 1; ++ii)
        {
            if (ii < 0 || ii >= s2)
                continue;
            for (int jj = x - 1; jj <= x + 1; ++jj)
            {
                if (jj < 0 || jj >= s1)
                    continue;
                // 当前网格的点集,都是i之前的点
                vector<int> vv = um[ii * s1 + jj];
                for (int kk = 0; kk < vv.size(); ++kk)
                {
                    double p1i = p1[i], p1v = p1[vv[kk]], p2i = p2[i], p2v = p2[vv[kk]];
                    double dis = (p1i - p1v) * (p1i - p1v) + (p2i - p2v) * (p2i - p2v);
                    dis = sqrt(dis);
                    if(dis > tmp) continue;
                    if (dis < tmp)
                    {
                        tmp = dis;
                        aa = vv[kk];
                        res2.clear();
                    }
                    lwd lwl;
                    lwl.i1 = vv[kk];
                    lwl.i2 = i;
                    res2.emplace_back(lwl);
                }
            }
        }
        // 重新划分：更新划分值，最近点对，划分结果
        if (tmp < s)
        {
            s = tmp;
            ss = s + EPS;
            a = aa, b = i;
            s1 = d1 / ss + 1, s2 = d2 / ss + 1;
            um.clear();
            for (int ii = 0; ii < i; ++ii)
            {
                int xx = (p1[ii] - mi1) / ss, yy = (p2[ii] - mi2) / ss;
                um[yy * s1 + xx].emplace_back(ii);
            }
        }
        um[y * s1 + x].emplace_back(i);
    }
}

int main(int argc, char **argv)
{
    srand((unsigned)time(NULL));
    int n = stoi(argv[1]);
    vector<int> p1(n), p2(n);
    for (int i = 0; i < n; ++i)
    {
        p1[i] = rand() % n;
        p2[i] = rand() % n;
    }
    // 期望线性法的最好情况
    // p1[0] = p1[1] = p1[2] = 0, p2[0] = p2[1] = p2[2] = 0;
    // for(int i = 3; i < n; ++i) {
    //     p1[i] = p1[i-1]+i;
    //     p2[i] = p2[i-1]+i;
    // }
    // 期望线性法的最坏情况
    // for(int l = 0, r = n-1; l < r; ++l, --r) {
    //     int tmp = p1[l];
    //     p1[l] = p1[r];
    //     p1[r] = tmp;
    //     tmp = p2[l];
    //     p2[l] = p2[r];
    //     p2[r] = tmp;
    // }
    // for(int i = 0; i < n; ++i) cout << p1[i] << "  " << p2[i] << endl;
    int a = -1, b = -1;
    // 暴力解法耗时和结果
    struct timeval start, end;
    gettimeofday(&start, NULL);
    simple(p1, p2, a, b);
    gettimeofday(&end, NULL);
    float cost_time = (end.tv_usec - start.tv_usec) / 1000000.0 + end.tv_sec - start.tv_sec;
    cout << cost_time << endl;
    //cout << a << "  " << b << endl;
    //cout << p1[a] << "," << p2[a] << "\t" << p1[b] << "," << p2[b] << endl;
    cout << res.size() << endl;
    for(int i = 0; i < res.size(); ++i) {
        cout << res[i].i1 << "  " << res[i].i2 << endl;
        cout << p1[res[i].i1] << "," << p2[res[i].i1] << "\t" << p1[res[i].i2] << "," << p2[res[i].i2] << endl;
    }
    cout << "******************************\n";
    // 期望线性解法
    gettimeofday(&start, NULL);
    linear(p1, p2, a, b);
    gettimeofday(&end, NULL);
    cost_time = (end.tv_usec - start.tv_usec) / 1000000.0 + end.tv_sec - start.tv_sec;
    cout << cost_time << endl;
    // cout << a << "  " << b << endl;
    // cout << p1[a] << "," << p2[a] << "\t" << p1[b] << "," << p2[b] << endl;
    cout << res2.size() << endl;
    for(int i = 0; i < res2.size(); ++i) {
        cout << res2[i].i1 << "  " << res2[i].i2 << endl;
        cout << p1[res2[i].i1] << "," << p2[res2[i].i1] << "\t" << p1[res2[i].i2] << "," << p2[res2[i].i2] << endl;
    }
    cout << "******************************\n";
    // 分治
    vector<pt> vp;
    pt ppt;
    for (int i = 0; i < n; ++i)
    {
        ppt.x = p1[i];
        ppt.y = p2[i];
        ppt.id = i;
        vp.emplace_back(ppt);
    }
    gettimeofday(&start, NULL);
    sort(vp.begin(), vp.begin() + n, cmp_x);
    mindist = 1E20;
    rec(vp, 0, n - 1);
    gettimeofday(&end, NULL);
    cost_time = (end.tv_usec - start.tv_usec) / 1000000.0 + end.tv_sec - start.tv_sec;
    cout << cost_time << endl;
    // cout << ansa << "  " << ansb << endl;
    // cout << p1[ansa] << "," << p2[ansa] << "\t" << p1[ansb] << "," << p2[ansb] << endl;
    cout << res1.size() << endl;
    for(lwd l : res1) {
        cout << l.i1 << "  " << l.i2 << endl;
        cout << p1[l.i1] << "," << p2[l.i1] << "\t" << p1[l.i2] << "," << p2[l.i2] << endl;
    }
    cout << "******************************\n";

    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
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
304
305
306

• 以下是旧版本，期望线性法存在问题：s为零会出问题，和参考不一致

#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 

using namespace std;

void simple(vector<int> p1, vector<int> p2, int &a, int &b)
{
    int n = p1.size();
    double dis = n * 2.2;
    for (int i = 0; i < n; ++i)
    {
        for (int j = i + 1; j < n; ++j)
        {
            double p1i = p1[i], p1j = p1[j], p2i = p2[i], p2j = p2[j];
            double tmp = (p1i - p1j) * (p1i - p1j) + (p2i - p2j) * (p2i - p2j);
            tmp = sqrt(tmp);
            if (tmp < dis)
            {
                dis = tmp;
                a = i;
                b = j;
            }
        }
    }
}

/***************************************  分治  *******************************************/

struct pt
{
    int x, y, id;
};

struct cmpx
{
    bool operator()(const pt &a, const pt &b) const
    {
        return a.x < b.x || (a.x == b.x && a.y < b.y);
    }
} cmp_x;

struct cmpy
{
    bool operator()(const pt &a, const pt &b) const { return a.y < b.y; }
} cmp_y;

// vector vp;
double mindist;
int ansa, ansb;

inline void upd_ans(const pt &a, const pt &b)
{
    double ax = a.x, bx = b.x, ay = a.y, by = b.y;
    double dist = sqrt((ax - bx) * (ax - bx) + (ay - by) * (ay - by));
    if (dist < mindist)
        mindist = dist, ansa = a.id, ansb = b.id;
}

void rec(vector<pt> &vp, int l, int r)
{
    if (r - l <= 3)
    {
        for (int i = l; i <= r; ++i)
            for (int j = i + 1; j <= r; ++j)
                upd_ans(vp[i], vp[j]);
        sort(vp.begin() + l, vp.begin() + r + 1, cmp_y);
        return;
    }

    int m = (l + r) >> 1;
    int midx = vp[m].x;
    rec(vp, l, m), rec(vp, m + 1, r);
    inplace_merge(vp.begin() + l, vp.begin() + m + 1, vp.begin() + r + 1, cmp_y);

    // vector t(vp.size());
    static pt t[100000];
    int tsz = 0;
    for (int i = l; i <= r; ++i)
        if (abs(vp[i].x - midx) < mindist)
        {
            t[tsz++] = vp[i];
        }
    for(int i = 0; i < tsz; ++i) {
        for(int j = i+1; j < tsz && t[i].y-t[j].y < mindist; ++j) {
            upd_ans(t[i], t[j]);
        }
    }
}

/***************************************  期望线性法  *******************************************/
// p1,p2是坐标数组，a,b是最近点对的索引
void linear(vector<int> p1, vector<int> p2, int &a, int &b)
{
    int n = p1.size(), mi1 = p1[0], ma1 = p1[0], mi2 = p2[0], ma2 = p2[0];
    // 求点的范围
    for (int i = 1; i < n; ++i)
    {
        if (mi1 > p1[i])
            mi1 = p1[i];
        if (ma1 < p1[i])
            ma1 = p1[i];
        if (mi2 > p2[i])
            mi2 = p2[i];
        if (ma2 < p2[i])
            ma2 = p2[i];
    }
    int d1 = ma1 - mi1, d2 = ma2 - mi2;
    // 以前两点的距离划分网格
    a = 0, b = 1;
    double p10 = p1[0], p11 = p1[1], p20 = p2[0], p21 = p2[1];
    double s = (p10 - p11) * (p10 - p11) + (p20 - p21) * (p20 - p21);
    s = sqrt(s);
    int s1 = d1 / s + 1, s2 = d2 / s + 1;
    // 将每个点放入对应网格
    unordered_map<int, vector<int>> um;
    for (int i = 0; i < n; ++i)
    {
        int x = (p1[i] - mi1) / s, y = (p2[i] - mi2) / s;
        um[y * s1 + x].emplace_back(i);
    }
    // 依次处理剩余点：计算该点与周围九个网格点的最小距离，小于当前划分值则重新划分
    for (int i = 2; i < n; ++i)
    {
        // 当前点所在网格
        int x = (p1[i] - mi1) / s, y = (p2[i] - mi2) / s;
        // 周围九个网格与当前点的最小距离
        double tmp = n * 2.2;
        // 周围九个网格与当前点最近的点的索引
        int aa = -1;
        // 周围九个网格
        for (int ii = y - 1; ii <= y + 1; ++ii)
        {
            if (ii < 0 || ii >= s2)
                continue;
            for (int jj = x - 1; jj <= x + 1; ++jj)
            {
                if (jj < 0 || jj >= s1)
                    continue;
                // 当前网格的点集
                vector<int> vv = um[ii * s1 + jj];
                for (int kk = 0; kk < vv.size(); ++kk)
                {
                    if (vv[kk] <= i && vv[kk] > 1)
                        continue;
                    double p1i = p1[i], p1v = p1[vv[kk]], p2i = p2[i], p2v = p2[vv[kk]];
                    double dis = (p1i - p1v) * (p1i - p1v) + (p2i - p2v) * (p2i - p2v);
                    dis = sqrt(dis);
                    if (dis < tmp)
                    {
                        tmp = dis;
                        aa = vv[kk];
                    }
                }
            }
        }
        // 重新划分：更新划分值，最近点对，划分结果
        if (tmp < s)
        {
            s = tmp;
            a = i, b = aa;
            s1 = d1 / s + 1, s2 = d2 / s + 1;
            um.clear();
            for (int ii = 0; ii < n; ++ii)
            {
                int xx = (p1[ii] - mi1) / s, yy = (p2[ii] - mi2) / s;
                um[yy * s1 + xx].emplace_back(ii);
            }
        }
    }
}

int main(int argc, char **argv)
{
    srand((unsigned)time(NULL));
    int n = stoi(argv[1]);
    vector<int> p1(n), p2(n);
    for (int i = 0; i < n; ++i)
    {
        p1[i] = rand() % n;
        p2[i] = rand() % n;
    }
    // 期望线性法的最好情况
    // p1[0] = 0, p2[0] = 0;
    // for(int i = 1; i < n; ++i) {
    //     p1[i] = p1[i-1]+i;
    //     p2[i] = p2[i-1]+i;
    // }
    // 期望线性法的最坏情况
    // for(int l = 0, r = n-1; l < r; ++l, --r) {
    //     int tmp = p1[l];
    //     p1[l] = p1[r];
    //     p1[r] = tmp;
    //     tmp = p2[l];
    //     p2[l] = p2[r];
    //     p2[r] = tmp;
    // }
    // for(int i = 0; i < n; ++i) cout << p1[i] << "  " << p2[i] << endl;
    int a = -1, b = -1;
    // 暴力解法耗时和结果
    struct timeval start, end;
    gettimeofday(&start, NULL);
    simple(p1, p2, a, b);
    gettimeofday(&end, NULL);
    float cost_time = (end.tv_usec - start.tv_usec) / 1000000.0 + end.tv_sec - start.tv_sec;
    cout << cost_time << endl;
    cout << a << "  " << b << endl;
    cout << p1[a] << "," << p2[a] << "\t" << p1[b] << "," << p2[b] << endl;
    // 期望线性解法
    gettimeofday(&start, NULL);
    linear(p1, p2, a, b);
    gettimeofday(&end, NULL);
    cost_time = (end.tv_usec - start.tv_usec) / 1000000.0 + end.tv_sec - start.tv_sec;
    cout << cost_time << endl;
    cout << a << "  " << b << endl;
    cout << p1[a] << "," << p2[a] << "\t" << p1[b] << "," << p2[b] << endl;
    // 分治
    vector<pt> vp;
    pt ppt;
    for (int i = 0; i < n; ++i)
    {
        ppt.x = p1[i];
        ppt.y = p2[i];
        ppt.id = i;
        vp.emplace_back(ppt);
    }
    gettimeofday(&start, NULL);
    sort(vp.begin(), vp.begin() + n, cmp_x);
    mindist = 1E20;
    rec(vp, 0, n - 1);
    gettimeofday(&end, NULL);
    cost_time = (end.tv_usec - start.tv_usec) / 1000000.0 + end.tv_sec - start.tv_sec;
    cout << cost_time << endl;
    cout << ansa << "  " << ansb << endl;
    cout << p1[ansa] << "," << p2[ansa] << "\t" << p1[ansb] << "," << p2[ansb] << endl;

    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
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