BZOJ 3199: [Sdoi2013]escape

$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;
}