$O(n^2)$建Voronoi图，求对偶图后BFS即可
用Canvas写了个可视化
想写增量算法和Fortune算法，可是我好菜啊orz
point的cmp写错了，调试了很久，要一直记得精度啊，用sgn函数，否则不满足偏序
半平面交抄板子都能抄错orz
此外写代码最好一气呵成，别磨叽，这东西我写了2day你敢信？
#include <iostream> #include <cmath> #include <cstdio> #include <cstring> #include <vector> #include <algorithm> #include <map> using namespace std; typedef double DB; const DB EPS = 1e-8; int sgn(DB x) { return x < -EPS ? -1 : EPS < x; } struct Point { DB x, y; Point(DB x = 0, DB y = 0):x(x), y(y) { } DB dot(Point v) { return x * v.x + y * v.y; } DB cross(Point v) { return x * v.y - y * v.x; } DB norm() { return sqrt(this->dot(*this)); } DB length() { return norm(); } Point normalize(); bool operator < (const Point&rhs) const { return sgn(x - rhs.x) < 0 || (sgn(x - rhs.x) == 0 && sgn(y - rhs.y) < 0); } }; typedef Point Vector; Point operator + (Point a, Point b) { return Point(a.x + b.x, a.y + b.y); } Point operator - (Point a, Point b) { return Point(a.x - b.x, a.y - b.y); } Point operator * (Point a, DB p) { return Point(a.x * p, a.y * p); } Point operator / (Point a, DB p) { return Point(a.x / p, a.y / p); } DB dot(Vector a, Vector b) { return a.dot(b); } DB cross(Vector a, Vector b) { return a.cross(b); } Point Point::normalize() { return (*this) = (*this) / this->norm(); } bool operator == (Point a, Point b) { return sgn(a.x - b.x) == 0 && sgn(a.y - b.y) == 0; } struct Line { Point p; Vector v; DB ang; Line() {} Line(Point p, Vector v):p(p), v(v) { v.normalize(); ang = atan2(v.y, v.x); } Point getInterPoint(Line b) { Line a = *this; return a.p + a.v * (b.p.cross(b.v) - a.p.cross(b.v)) / a.v.cross(b.v); } }; bool cmpAng(const Line&a, const Line&b) { return a.ang < b.ang; } bool onLeft(Line l, Point p) { return l.v.cross(p - l.p) > 0; } vector<Point> halfplaneInter(vector<Line> l) { vector<Point> ret; sort(l.begin(), l.end(), cmpAng); int first, last, n = l.size(); Point *p = new Point[n]; Line *q = new Line[n]; q[first = last = 0] = l[0]; for (int i = 1; i < n; ++i) { while (first < last && !onLeft(l[i], p[last - 1])) --last; while (first < last && !onLeft(l[i], p[first])) ++first; q[++last] = l[i]; if (sgn(fabs(q[last].v.cross(q[last - 1].v))) == 0) { --last; if (onLeft(q[last], l[i].p)) q[last] = l[i]; } if (first < last) p[last - 1] = q[last - 1].getInterPoint(q[last]); } while (first < last && !onLeft(q[first], p[last - 1])) --last; if (last - first <= 1) return ret; p[last] = q[last].getInterPoint(q[first]); for (int i = first; i <= last; ++i) ret.push_back(p[i]); delete p; delete q; p = NULL, q = NULL; return ret; } namespace DCEL { struct HalfEdge; struct Face; struct Vertex { int id; HalfEdge *parEdge; Vertex() { parEdge = NULL; } }; struct HalfEdge { HalfEdge *twin, *succ; int orig; Face *parFace; HalfEdge() { twin = succ = NULL; orig = 0; parFace = NULL; } }; struct Face { int site; int id; HalfEdge *eHead; Face() { eHead = NULL; } }; }; using namespace DCEL; struct VoronoiN2 { vector<Point> sites; vector<Point> pointPool; map<Point, unsigned int> pointId; map<pair<int, int>, HalfEdge* > conj; vector<Face*> faces; int getId(Point x) { if (pointId.find(x) == pointId.end()) { pointId.insert(make_pair(x, pointPool.size())); pointPool.push_back(x); } return pointId.find(x)->second; } VoronoiN2(Point*p, int n, Point a, Point c) { for (int i = 0; i < n; ++i) sites.push_back(p[i]); vector<Line> base; vector<Line> ext; vector<Point> vertex; Point b = Point(a.x, c.y), d = Point(c.x, a.y); base.push_back(Line(a, d - a)); base.push_back(Line(d, c - d)); base.push_back(Line(c, b - c)); base.push_back(Line(b, a - b)); Face *outerFace = new Face; outerFace->id = 0; outerFace->site = -1; outerFace->eHead = new HalfEdge; outerFace->eHead->orig = -1; outerFace->eHead->parFace = outerFace; faces.push_back(outerFace); for (int i = 0; i < n; ++i) { ext.assign(base.begin(), base.end()); for (int j = 0; j < n; ++j) if (j != i) { Vector v = sites[j] - sites[i]; ext.push_back(Line(p[i] + v / 2, Point(-v.y, v.x))); } vertex = halfplaneInter(ext); int m = vertex.size(); Face *newFace = new Face; HalfEdge *ptn, *newPtn; newFace->site = getId(sites[i]); ptn = newFace->eHead = new HalfEdge; ptn->orig = getId(vertex[0]); ptn->parFace = newFace; for (int i = 1; i < m; ++i) { newPtn = new HalfEdge(); newPtn->orig = getId(vertex[i]); newPtn->parFace = newFace; ptn->succ = newPtn; conj.insert(make_pair(make_pair(getId(vertex[i - 1]), getId(vertex[i])), ptn)); ptn = newPtn; } ptn->succ = newFace->eHead; conj.insert(make_pair(make_pair(getId(vertex[m - 1]), getId(vertex[0])), ptn)); newFace->id = faces.size(); faces.push_back(newFace); } memset(head, -1, sizeof head); edgeCnt = 0; for (int k = 1; k < (int)faces.size(); ++k) { HalfEdge*i = faces[k]->eHead; do { HalfEdge*j = i->succ; if (i->twin == NULL) { map<pair<int, int>, HalfEdge* >::iterator jt = conj.find(make_pair(j->orig, i->orig)); if (jt == conj.end()) { i->twin = outerFace->eHead; } else { i->twin = jt->second; } } //printf("%d -> %d ", k, i->twin->parFace->id); to[edgeCnt] = i->twin->parFace->id; nxt[edgeCnt] = head[k]; head[k] = edgeCnt++; i = j; } while (i != faces[k]->eHead); } //puts(""); } static const int MAXN = 1000; static const int MAXE = 10000; int edgeCnt; int head[MAXN]; int to[MAXE], nxt[MAXE]; void buildDualGraph() { // memset(head, -1, sizeof head); // edgeCnt = 0; // for (int i = 1; i < (int)faces.size(); ++i) { // HalfEdge*j = faces[i]->eHead; // do { // to[edgeCnt] = j->twin->parFace->id; // nxt[edgeCnt] = head[i]; // printf("%d ", j->twin->parFace->id); // //printf("%d -> %d ", i, to[edgeCnt]); // head[i] = edgeCnt++; // j = j->succ; // } while (j != faces[i]->eHead); // } } int findLocation(Point a) { Point p, v; for (int i = 1; i < (int)faces.size(); ++i) { bool ok = true; HalfEdge*j = faces[i]->eHead; do { p = pointPool[j->orig]; v = pointPool[j->succ->orig]; if ((v - p).cross(a - p) < 0) { ok = false; break; } j = j->succ; } while (j != faces[i]->eHead); if (ok) { return faces[i]->id; } } return 0; } void print() { printf("sitePoints = new Array("); printf("{x: %lf, y: %lf}", sites[0].x, sites[0].y); for (int i = 1; i < (int)sites.size(); ++i) { printf(", {x: %lf, y: %lf}", sites[i].x, sites[i].y); } puts(");"); printf("faces = new Array( "); for (int i = 1; i < (int)faces.size(); ++i) { HalfEdge*j = faces[i]->eHead; printf("new Array({x: %lf, y: %lf}", pointPool[j->orig].x, pointPool[j->orig].y); j = j->succ; while (j != faces[i]->eHead) { printf(", {x: %lf, y: %lf}", pointPool[j->orig].x, pointPool[j->orig].y); j = j->succ; } if (i != (int)faces.size() - 1) puts("),"); else puts(")"); } puts(");"); } void solve(Point a) { buildDualGraph(); int x = findLocation(a); // printf("%d ", x); int *que = new int[10000], *dis = new int[1000]; int ql, qr; memset(dis, 0x3f, (sizeof (int)) * 1000); dis[x] = 0; que[ql = qr = 0] = x; while (ql <= qr) { int u = que[ql++]; for (int i = head[u]; i != -1; i = nxt[i]) { int v = to[i]; if (dis[u] + 1 <= dis[v]) { dis[v] = dis[u] + 1; que[++qr] = v; } } } printf("%d ", dis[0]); } }; int main() { #ifdef lol freopen("3199.in", "r", stdin); freopen("3199.out", "w", stdout); #endif int t; scanf("%d", &t); while (t--) { int n; scanf("%d", &n); Point a, b; cin >> a.x >> a.y; cin >> b.x >> b.y; Point *p = new Point[n]; for (int i = 0; i < n; ++i) cin >> p[i].x >> p[i].y; VoronoiN2 *v = new VoronoiN2(p, n, Point(0, 0), a); //v->print(); v->solve(b); } return 0; }