- 参考这里最下面
- 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 <unordered_map>
#include <vector>
#include <set>
#include <iostream>
#include <algorithm>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <sys/time.h>
#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;
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);
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]);
}
}
}
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;
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;
}
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 << 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 << 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 << 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;
}
- 以下是旧版本,期望线性法存在问题:s为零会出问题,和参考不一致
#include <unordered_map>
#include <vector>
#include <iostream>
#include <algorithm>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <sys/time.h>
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;
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);
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]);
}
}
}
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;
}
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;
}