[模板]计算几何

平面最近点对

将平面内点按\(x\)坐标排序,分治\(x\)坐标,设\(ret=min(f(l,mid),f(mid+1,r))\),将\(x\in[mid-ret,mid+ret]\)内的点按\(y\)坐标排序,算每个点与相邻的\(6\)个点的距离找最优解即可.

时间复杂度:\(O(nlogn)\).

#define N 100005
#define INF 1e15
struct point{
    double x,y;
}p[N];
inline double sqr(double k){
    return k*k;
}
inline double dis(point x,point y){
    return sqrt(sqr(x.x-y.x)+sqr(x.y-y.y));
}
inline bool cmpx(point x,point y){
    if(x.x!=y.x) return x.x<y.x;
    return x.y<y.y;
}
inline bool cmpy(point x,point y){
    if(x.y!=y.y) return x.y<y.y;
    return x.x<y.x;
}
inline double min_d(int l,int r){
    double ret=INF;
    if(r-l<=20){
        for(int i=l;i<r;++i)
            for(int j=i+1;j<=r;++j)
                ret=min(ret,dis(p[i],p[j]));
        return ret;
    }
    int mid=l+r>>1;
    ret=min(min_d(l,mid),min_d(mid+1,r)); 
    while(p[l].x+ret<p[mid].x) ++l;
    while(p[r].x-ret>p[mid].x) --r;
    sort(p+l,p+1+r,cmpy);
    for(int i=l;i<r;++i)
        for(int j=min(r,i+6);j>i;--j)
            ret=min(ret,dis(p[i],p[j])); 
    sort(p+l,p+1+r,cmpx);
    return ret;
}
inline double min_dis(){
    sort(p+1,p+1+n,cmpx);
    return min_d(1,n);
}

推荐

http://www.cnblogs.com/xdruid/archive/2012/05/27/CP.html

凸包

时间复杂度:\(O(nlogn)\).

#define N 100005
#define eps 1e-13
struct point{
    int x,y;double an;
}a[N],v[N];
int n,u,vn;
inline double sqr(int k){
    return (double)(k*k);
}
inline point dec(point x,point y){
    return (point){x.x-y.x,x.y-y.y,0.0};
}
inline int mul(point x,point y){
    return x.x*y.y-y.x*x.y;
}
inline double dis(point x,point y){
    return sqrt(sqr(abs(x.x-y.x))+sqr(abs(x.y-y.y)));
}
inline bool cmp(point x,point y){
    if(fabs(x.an-y.an)<eps)
        return dis(x,a[1])>dis(y,a[1]); 
    return x.an<y.an;
} 
inline void convex(){
    u=1;
    for(int i=2;i<=n;++i)
        if((a[i].x<a[u].x)||(a[i].x==a[u].x&&a[i].y<a[u].y)) u=i;
    int tmp=a[u];a[u]=a[1];a[1]=tmp;
    for(int i=2;i<=n;++i)
        a[i].an=atan2(a[i].y-a[1].y,a[i].x-a[1].x);
    sort(a+2,a+1+n,cmp);
    v[++vn]=a[1];v[++vn]=a[2];a[++n]=a[1];
    for(int i=3;i<=n;++i){
        if(fabs(a[i].an-a[i-1].an)<eps) continue;
        while(vn>1&&mul(dec(a[i],v[vn-1]),dec(v[vn],v[vn-1]))>0) --vn;
        v[++vn]=a[i];
    }
} 

旋转卡壳

时间复杂度:\(O(n)\).

#define N 100005
struct point{
    int x,y;
}a[N];
int n;
inline point dec(point x,point y){
    return (point){x.x-y.x,x.y-y.y};
}
inline int mul(point x,point y){
    return x.x*y.y-x.y*y.x;
}
inline double sqr(int k){
    return (double)(k*k);
}
inline double dis(point x){
    return sqrt(sqr(x.x)+sqr(x.y));
}
inline int Next(int k){
    if(++k>n) return 1;
    return k;
}
inline double rorate(){
	point p;
    double di,dia=0.0;
    if(n==1) return dia;
    for(int i=1,j=2;i<=n;++i){
    	p=dec(a[Next(i)],a[i]);
        while(abs(mul(p,dec(a[j],a[i])))<abs(mul(p,dec(a[Next(j)],a[i]))))
            j=Next(j);
        dia=max(dia,max(dis(dec(a[i],a[j])),dis(dec(a[Next(i)],a[Next(j)]))));
    }
    return dia;
}

判断点是否在凸多边形内

inline point dec(point x,point y){
    return (point){x.x-y.x,x.y-y.y};
}
inline double mul(point x,point y){
    return x.x*y.y-y.x*x.y;
}
inline bool onseg(point p,point a,point b){  
    if(cmp(mul(dec(a,p),dec(b,p)))) return false;  
    return cmp(a.x-p.x)*cmp(b.x-p.x)<=0&&cmp(a.y-p.y)*cmp(b.y-p.y)<=0;
} 
inline int chk(point p){
    int cnt=0,k,d1,d2;
    for(int i=1;i<=n;++i){
        if(onseg(p,a[i],a[i+1])) return -1;
        k=cmp(mul(dec(a[i+1],a[i]),dec(p,a[i])));
        d1=cmp(a[i].y-p.y);d2=cmp(a[i+1].y-p.y);
        if(k>0&&d1<=0&&d2>0) ++cnt;
        if(k<0&&d2<=0&&d1>0) --cnt;
    }
    if(cnt) return 1;
    return 0;
}

半平面交

将半平面按极角排序(靠近法向量方向的半平面更小)
双端队列维护半平面交:顺序加入半平面,当队尾的2个半平面的交点在当前半平面外，删除队尾的半平面;当队头的2个半平面的交点在当前半平面外，删除队头的半平面.
\(P.S.\)形如\(ax+by+c\;\geq\;0\)的方程组的解\(A(x,y)\)满足条件\(\overrightarrow{BA}\;\times\;\overrightarrow{BC}\;\geq\;0(B,C\;\in\;ax+by+c=0)\)(\(ax+by+c<0\)反之)

typedef long double ld;
struct point{
    ld x,y;
}p[N];
struct line{
    point s,e;ld an;
}a[N],q[N];
int n,m,h,t;
inline point dec(point x,point y){
    return (point){x.x-y.x,x.y-y.y};
}
inline ld mult(point x,point y){
    return x.x*y.y-x.y*y.x;
}
inline bool cmp(line x,line y){
    if(x.an==y.an) return mult(dec(x.e,x.s),dec(y.s,x.s))>0;
    return x.an<y.an;
}
inline point inter(line a,line b){
    ld s1,s2,t;point ret;
    s1=mult(dec(b.e,a.s),dec(a.e,a.s));
    s2=mult(dec(a.e,a.s),dec(b.s,a.s));
    t=s2/(s1+s2);
    ret.x=b.s.x+(b.e.x-b.s.x)*t;
    ret.y=b.s.y+(b.e.y-b.s.y)*t;
    return ret;
}
inline bool chk(line x,line y,line z){
    point a=inter(x,y);
    return mult(dec(a,z.s),dec(z.e,z.s))>0;
}
inline bool hpi(int k){
	for(int i=1;i<=n;++i) 
	 	a[i].an=atan2(a[i].e.y-a[i].s.y,a[i].e.x-a[i].s.x);
	sort(a+1,a+1+n,cmp);
    h=1;t=0;
    for(int i=1;i<=n;++i){
        while(h<t&&chk(q[t],q[t-1],a[i])) --t;
        while(h<t&&chk(q[h],q[h+1],a[i])) ++h;
        q[++t]=a[i];
    }
    while(h<t&&chk(q[t],q[t-1],q[h])) --t;
    while(h<t&&chk(q[h],q[h+1],q[t])) ++h;
    return t-h+1>=3;
}