Codeforces 549E Sasha Circle

problem

题意

• 给出两堆点，问是否存在一个圆，能把其中的一堆全部包含进去但另外一堆全部在圆外面。

思路

这题当时比赛的时候没有人写出来。赛后的题解也迟迟没有给出。给出的题解我只看懂了前半部分。。后半部分应该是比较屌的东西。我的代码还是参考了别人AC的代码写的。思路就是固定两个点，然后对于这一堆点，我们可以求出一个最大半径和一个最小半径，而对于另外的点，我们可以逐个对这个最大和最小半径（其实是圆心到两点中点的距离）进行更新。如果更新到最后，最大半径和最小半径之间仍有范围，那么说明这种情况有解。我能想到枚举所有两个点的组合，但是不知道为什么向代码里这样递归写是正确的。。这题细节颇多，其实就是对方向要有一个很好的认识。。感觉是一道不错的计算几何题目。。

题解翻译

• 看到半夜，终于看懂题解后半部分的意思。真的是太巧妙了，值得去写一个中文版。
• 首先，我们建立一个空间坐标系，也就是在原坐标系上加上z轴
  画出一个曲面$z = x^2+y^2$，然后我们把给定的所有点都投影到这个曲面上。即$(x,y,x^2+y^2)$
• 然后我们会发现一件事情，即任意一个不垂直于xoy平面的平面，和这个曲面相交得到的曲线，在xoy平面上的投影都是一个圆。
  这是可以证明的，联立$ax+by+z = 1$和$x^2+y^2 = z$
  得到 $(x+\frac{1}{2}a)^2+(y+\frac{1}{2}b)^2 = c+\frac{1}{4}(a^2+b^2)$
  对于xoy平面来说这是一个任意的圆。

• 所以我们可以得到一个结论，如果原题给定的两个集合中的点在曲面上投影，形成的两个新的集合，可以用一个不垂直于xoy平面的平面分开，那么在xoy平面上，肯定存在这样的一个圆能完成这件事。
  给出一张图，是题解里的，更加容易理解这件事情。
  

• 有了这个结论就好办了。我们对其中的一个集合求凸包，然后求出这些点构成的上凸平面。这可以通过分治的手段来做。求出这些平面之后，我们可以从平面上来判断，另外一个集合是否能和当前集合通过这个面分开（算一下到中垂线的距离）

• 这道题的巧妙之处在于把平面几何通过转换成立体几何，更加直观的证明了一个不太容易想到的结论。个人觉得是一道非常不错的题目。

AC代码

/*************************************************************************
  > File Name: pe.cpp
  > Author: znl1087
  > Mail: loveCJNforever@gmail.com 
  > Created Time: 日  6/ 7 11:19:33 2015
 ************************************************************************/

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cmath>
#include <algorithm>
#include <string>
#include <cstdlib>
#include <vector>
#include <set>
#include <map>
#include <cfloat>
using namespace std;
const double PI = acos(-1.0);
struct Point{
    double x,y;
    Point(double x = 0,double y = 0):x(x),y(y){}
};
typedef Point Vec;
Vec operator + (Vec a,Vec b){return Vec(a.x+b.x,a.y+b.y);}
Vec operator - (Point a,Point b){return Vec(a.x-b.x,a.y-b.y);}
Vec operator * (Vec a,double p){return Vec(a.x*p,a.y*p);}
Vec operator / (Vec a,double p){return a*(1/p);}
bool operator < (const Point& a,const Point& b){
    return a.x < b.x || (a.x == b.x && a.y < b.y);
}
const double EPS = 1e-10;
int dcmp(double x){
    if(fabs(x) < EPS) return 0;else return x < 0 ? -1 : 1;
}
bool operator == (const Point& a,const Point& b){
    return dcmp(a.x-b.x)==0 && dcmp(a.y-b.y)==0;
}
double ang(Vec v){return atan2(v.y,v.x);}

double Dot(Vec A,Vec B){return A.x * B.x + A.y * B.y;}
double Length(Vec A){return sqrt(Dot(A,A));}
double Angle(Vec A,Vec B){return acos(Dot(A,B)/Length(A)/Length(B));}

double Cross(Vec A,Vec B){return A.x * B.y - A.y * B.x;}

Point CirCenter(Point a,Point b,Point c){
    double a1 = b.x - a.x,b1 = b.y - a.y,c1 = (a1 * a1 + b1 * b1)/2;
    double a2 = c.x - a.x,b2 = c.y - a.y,c2 = (a2 * a2 + b2 * b2)/2;
    double d = a1 * b2 - a2 * b1;
    return Point(a.x + (c1 * b2 - c2 * b1)/d,a.y + (a1 * c2 - a2 * c1)/d);
}

vector<Point> ConvexHull(vector<Point> p){
    sort(p.begin(),p.end());
    p.erase(unique(p.begin(),p.end()),p.end());
    int n = p.size();
    vector<Point> ch(n+1);
    int m = 0;
    for(int i=0;i<n;i++){
        while(m>1 && dcmp(Cross(ch[m-1] - ch[m-2],p[i]-ch[m-2]))<=0)m--;
        ch[m++] = p[i];
    }
    int k = m;
    for(int i=n-2;i>=0;i--){
        while( m>k && dcmp(Cross(ch[m-1] - ch[m-2],p[i] - ch[m-2])) <= 0)m--;
        ch[m++] = p[i];
    }
    if( n > 1)m--;
    ch.resize(m);
    return ch;
}

bool judge(const vector<Point>& a,const vector<Point>& b,int i,int j){
    double minr =-DBL_MAX,maxr = DBL_MAX;
    int mid = -1;
    int na = a.size(),nb = b.size();
    bool flag = true;
    for(int k=i+1;k<j;k++){
        if(dcmp(Cross(a[k]-a[i],a[k]-a[j]))==0){
            flag = false;
            break;
        }
        Point center = CirCenter(a[i],a[j],a[k]);
        double r = Length(center - (a[i]+a[j])/2.0);
        if(dcmp(Cross(center-a[i],a[j]-a[i]))*dcmp(Cross(a[k]-a[i],a[j]-a[i]))<0)r = -r;
        if(dcmp(r-minr)>0)
            minr = r,mid = k;
    }
    if(flag){
        for(int k=(j+1)%na;;k = (k+1)%na){
            if(k == i)break;
            if(dcmp(Cross(a[k]-a[i],a[k]-a[j]))==0){
                flag = false;
                break;
            }
            Point center = CirCenter(a[i],a[j],a[k]);
            double r = Length(center - (a[i]+a[j])/2.0);
            if(dcmp(Cross(center-a[i],a[j]-a[i]))*dcmp(Cross(a[k]-a[i],a[j]-a[i]))>0)r = -r;
            if( dcmp(r-maxr)<0)maxr = r;
        }
    }
    //printf("maxr: %lf minr: %lf i: %d j:%d\n",maxr,minr,i,j);
    bool ok = true;
    if(flag){
        for(int k=0;k<nb;k++){
            if(dcmp(Cross(b[k]-a[i],a[j]-a[i]))==0 && dcmp(Dot(a[i]-b[k],a[j]-b[k]))<=0 ){
                ok = false;
                break;
            }
            Point center = CirCenter(a[i],a[j],b[k]);
            double r = Length(center - (a[i]+a[j])/2.0);
            if(dcmp(Cross(center-a[i],a[j]-a[i]))*dcmp(Cross(b[k]-a[i],a[j]-a[i]))<0)r = -r;
            //printf("r :%lf x: %lf y:%lf\n",r,b[k].x,b[k].y);
            if(dcmp(Cross(b[k]-a[i],a[j]-a[i]))>0){
                if(dcmp(r-maxr)<0)maxr = r;
            }
            else if(dcmp(Cross(b[k]-a[i],a[j]-a[i]))<0){
                r = -r;
                if(dcmp(minr-r)<0)minr = r;
            }
            if( dcmp(maxr-minr) <= 0){
                ok = false;
                break;
            }
        }
    }
    else ok = false;
    /*if(ok){
      printf("%lf %lf\n%lf %lf\n",a[i].x,a[i].y,a[j].x,a[j].y);
      printf("%lf %lf\n",minr,maxr);
      }*/
    return ok || (mid >=0 && (judge(a,b,i,mid) || judge(a,b,mid,j))); 
}
int n,m;
int main(){
    scanf("%d%d",&n,&m);
    vector<Point> ina(n),inb(m);
    int x,y;
    for(int i=0;i<n;i++)
        scanf("%d%d",&x,&y),ina[i].x = (double)x,ina[i].y = (double)y;
    for(int i=0;i<m;i++)
        scanf("%d%d",&x,&y),inb[i].x = (double)x,inb[i].y = (double)y;
    vector<Point> a = ConvexHull(ina),b = ConvexHull(inb);
    puts((judge(a,inb,0,a.size()-1) || judge(b,ina,0,b.size()-1)) ? "YES" : "NO");
    return 0;
}