Sasha Circle (CF549E)

题目大意

平面上有两个点集，点集大小分别为 $n,m(n,m\leq 10^4)$，所有的点都在整点处，且坐标的绝对值不超过 $C(C\leq 10^4)$。求一个圆将两个点集隔开。

引理

• 1. 整点集的凸包大小上限为 $O(C^\frac{2}{3})$。证明：我不会。
• 2. 抛物面 $z=x^2+y^2$与平面 $ax+by+cz=1$的交为一个椭圆。这个椭圆在平面 $z=0$上的投影是一个圆。证明：联立两个方程不难发现这个结论。

记号

• 令 $convex(A)$为点集 $A$的凸包。

• 令 $|A|$为集合 $A$的大小 (所含元素的个数)。

方法一

假定点集 $A$被划分到了圆内，点集 $B$在圆外。一个非常直观的想法是：只保留 $A$的凸包上的点即可。这么做一定是对的，正确性非常显然。利用这一步操作，$|A|$被缩减到了 $|convex(A)|=O(C^{\frac{2}{3}})$。

但是对于 $B$就不一样了。$B$中的每一个点似乎都有可能限制最终圆的大小和位置。所以，我们可能需要保留 $B$中所有的点。

为了让划分圆包含所有 $A$中的点，不包含 $B$中的点，我们可以想到的一个充要条件是：圆心到任意一个 $A$中的点的距离小于到任意一个 $B$中的点的距离。根据这个条件，我们枚举 $a\in A,b\in B$，则圆心必然在 $a,b$中垂线划分出的包含 $a$的半平面内。

我们总共会得到 $O(|convex(A)||B|)$个半平面，再通过半平面交判定合法性。总时间复杂度为 $O(|convex(A)||B|\log(|convex(A)||B|))$。

这个时间复杂度有点大，不能通过本题。

方法二

本题还有一个性质，我们一定可以把最终的划分圆尽量缩小，使得它经过至少两个 $A$中的点。

如果枚举 $A$中的每一个点，以这个点为反演中心做圆的反演，那么问题就转化成了：求平面上一条直线，使直线的两侧各只有一种颜色的点。

这种做法同样只能做到 $O(|convex(A)||B|\log(|convex(A)||B|))$的复杂度，不能通过本题。

方法三

我们可以枚举最终的划分圆经过的两个点，并判定是否存在一个合法的圆心，使得这两个点和这个圆心确定的圆满足题意。

实际上，最终的划分圆可能经过的 $A$中的点对的个数是 $O(|convex(A)|)$的。为什么呢？

做一个二维到三维的变换：$f(x,y)\rightarrow (x,y,x^2+y^2)$。这些点全部都会落在一个抛物面 $P:z=x^2+y^2$上。而任何一个划分圆都可以被表示为某一个平面 $F:ax+by+cz=1$和抛物面 $P$的交在平面 $z=0$上的投影，而且圆内的点即为 $F$下方的点，圆外的点即为 $F$上方的点。

一个合法的划分圆必然将两种点完全分开。因而，我们可以求出变换后的 $A_3,B_3$的凸包 (用下标表示维度)，并判定是否存在一个平面同时与 $convex(A_3),convex(B_3)$相切。如果不存在，即意为 $A_3,B_3$的凸包相交因而不存在合法的划分圆，反之意为存在一个合法的划分圆。

我们注意到，我们只需要求 $B_3$的上凸壳和 $A_3$的下凸壳。而 $B_3$的上凸壳上的点一定是 $B_2$的凸包上的点，而 $A_3$的下凸壳上的点就是整个 $A_2$。

为了找到一个合法的划分平面，强制让划分平面与 $convex(A_3)$的上凸壳的棱相切即可。由于 $convex(A_3)$的上凸壳的棱数只有 $O(|convex(A)|)$，因此最终的划分圆可能经过的 $A$中的点对的个数是 $O(|convex(A)|)$的。

枚举 $convex(A_3)$的每一条棱，以及 $convex(B_3)$的每一个点，并判断它们确定的平面是否是合法的划分平面即可。这样做略微有些困难，因而我们可以将这个问题再转化成原模型下的问题去解决：

枚举 $convex(A_3)$的每一条棱对应的 $A_2$中的点对，用以确定划分圆圆心所在的直线，再枚举 $A_2,B_2$中的每一个点，用来确定圆心所在直线上的区间。如果最后这些区间有交，即判定为有解。

问题是如何求出 $convex(A_3)$呢？我们注意到，$A_3$的上凸壳的每一个面对应的平面都满足：该平面对应的划分圆包含所有 $A_2$中的点。因此，我们 $A_3$的上凸壳又对应了 $convex(A_2)$的一个三角剖分，这个三角剖分满足：任意一个三角形的外接圆都包含 $A_2$的每一个点。

这个三角剖分可以用分治来求解。给定一个凸包，任选一条边 $(a,b)$，在凸包上找到一个 $c$，使得 $(a,b,c)$构成的三角形包含整个凸包，然后递归处理 $a\cdots c$和 $c\cdots b$中的点即可。该算法的正确性不难证明。

凭借这个算法，我们就可以在 $O(|convex(A_2)|^2)$的时间内求出 $convex(A_3)$的上凸壳，然后花费 $O(|convex(A_2)|(|convex(A_2)|+|B|))$的时间求解原问题。

代码实现

这里给出第三种方法的实现：

#include <bits/stdc++.h>

using namespace std;

const double eps = 1e-5;

template <typename T> void cmin(T &x, const T &y)
{
    if(y < x) x = y;
}

template <typename T> void cmax(T &x, const T &y)
{
    if(y > x) x = y;
}

template <typename T> void read(T &x)
{
    x = 0; char c = getchar(); bool f = 0;
    while(!isdigit(c) && c!='-') c = getchar();
    if(c == '-') f = 1, c = getchar();
    while(isdigit(c)) x = x*10+c-'0', c = getchar();
    if(f) x = -x;
}

struct vec
{
    double x, y;

    vec (const double &x0 = 0, const double &y0 = 0) : x(x0), y(y0) {}
    vec operator + (const vec &t) const {return vec(x+t.x, y+t.y);}
    vec operator - (const vec &t) const {return vec(x-t.x, y-t.y);}
    vec operator * (const double &t) const {return vec(x*t, y*t);}
    double len2 () const {return x*x + y*y;}
    bool operator < (const vec &t) const {return (fabs(x-t.x)>eps) ? (x<t.x) : (y<t.y);}
    vec rot_90_a () const {return vec(-y, x);}
};

double crs(const vec &a, const vec &b)
{
    return a.x*b.y - a.y*b.x;
}

double dot(const vec &a, const vec &b)
{
    return a.x*b.x + a.y*b.y;
}

vector<vec> convex(vector<vec> src)
{
    vector<vec> stk;
    sort(src.begin(), src.end());
    for(auto i=src.begin(); i!=src.end(); i++)
    {
        while(stk.size()>=2 && crs(*i-stk.end()[-1], stk.end()[-1]-stk.end()[-2])>=-eps) stk.pop_back();
        stk.emplace_back(*i);
    }
    int top = stk.size();
    for(auto i=src.rbegin(); i!=src.rend(); i++)
    {
        while(stk.size()>top && crs(*i-stk.end()[-1], stk.end()[-1]-stk.end()[-2])>=-eps) stk.pop_back();
        stk.emplace_back(*i);
    }
    stk.pop_back();
    return stk;
}

int n, m;
vector<vec> A, B;

void input()
{
    read(n), read(m);
    for(int i=1; i<=n; i++)
    {
        vec tmp;
        read(tmp.x), read(tmp.y);
        A.emplace_back(tmp);
    }
    for(int i=1; i<=m; i++)
    {
        vec tmp;
        read(tmp.x), read(tmp.y);
        B.emplace_back(tmp);
    }
}

vec lin_lin_int(vec p1, vec v1, vec p2, vec v2)
{
    double k = crs(p2-p1, v2) / crs(v1, v2);
    return p1 + v1*k;
}

bool try_plain(vec a, vec b, const vector<vec> &in, const vector<vec> &out)
{
    double mxc = +1e30, mnc = -1e30;
    vec mid = (a+b)*0.5, per = (b-a).rot_90_a();
    for(auto i : in)
        if(crs(b-a, i-a) > eps)
            cmax(mnc, crs(b-a, lin_lin_int((i+a)*0.5, (i-a).rot_90_a(), mid, per)-a));
        else if(crs(b-a, i-a) < -eps)
            cmin(mxc, crs(b-a, lin_lin_int((i+a)*0.5, (i-a).rot_90_a(), mid, per)-a));
        else if(dot(b-a, i-a)>dot(b-a, b-a) || dot(b-a, i-a)<0) return false;
    for(auto i : out)
        if(crs(b-a, i-a) > eps)
            cmin(mxc, crs(b-a, lin_lin_int((i+a)*0.5, (i-a).rot_90_a(), mid, per)-a));
        else if(crs(b-a, i-a) < -eps)
            cmax(mnc, crs(b-a, lin_lin_int((i+a)*0.5, (i-a).rot_90_a(), mid, per)-a));
        else if(dot(b-a, i-a)<dot(b-a, b-a) && dot(b-a, i-a)>0) return false;
    if(mxc > mnc+eps) return true;
    else return false;
}

bool check(vector<vec>::iterator in_beg, vector<vec>::iterator in_end, const vector<vec> &in, const vector<vec> &out)
{
    if(try_plain(*in_beg, *(in_end-1), in, out)) return true;
    if(in_end-in_beg < 3) return 0;
    vec lft = *in_beg, rgt = *(in_end-1);
    double mxd = -1e18, tmp;
    vector<vec>::iterator mxp;
    for(auto i=in_beg+1; i!=in_end-1; i++)
    {
        double val = - dot(*i-lft, *i-rgt) / crs(*i-lft, *i-rgt);
        if(val > mxd) mxd = val, mxp = i;
    }
    if(check(in_beg, mxp+1, in, out) || check(mxp, in_end, in, out)) return true;
    else return false;
}

int main()
{
    input();
    vector<vec> cA = convex(A), cB = convex(B);
    if(check(cA.begin(), cA.end(), A, B)) puts("YES");
    else if(check(cB.begin(), cB.end(), B, A)) puts("YES");
    else puts("NO");
    return 0;
}