bzoj 2961: 共点圆 cdq分治

       这道题目数据很弱不保证我的程序完全正确QAQ。

       另外这道题目在2013年集训队论文中有提到。以下是窝的口胡>.<：

       对于一个点(x0,y0)和一个圆心为(x,y)的圆，显然当(x0-x)^2+(y0-y)^2<=x^2+y^2时点在圆内，化简得到：2y·y0>=-2x·x0+x0^2+y0^2，然后可以把2y0除到右边去，当y0>0时有：

       y>=(-x0/y0)x+(x0^2+y0^2)/2y0

       显然这是一个半平面交的形式，因此圆心(x,y)需要在直线y=(-x0/y0)x+(x0^2+y0^2)/2y0的上面。要让所有的圆心都满足这个条件，我们可以对所有的圆心求一个下凸壳，然后到这个下凸壳上面去二分查找斜率最接近-x0/y0的然后只要这个圆心满足条件则所有的圆心满足条件。

       当y0<0时同理，维护一个上凸壳即可。

       因此得到一个在线的用平衡树维护全凸壳的做法，时间复杂度O(NlogN)。

       由于没有强制在线，因此采用cdq分治，[l,r]中得到[l,mid]的上下凸壳更新[mid+1,r]中的答案即可。可以采用线性扫描维护凸壳的方法。因此为O(NlogN)。

AC代码如下：

#include<iostream>
#include<cstdio>
#include<cmath>
#include<algorithm>
#define sqr(x) (x)*(x)
#define inf 1e40
#define eps 1e-10
#define N 500005
using namespace std;

int n,cnt,q1[N],q2[N]; struct node{ int op,id,p; double x,y,k; }a[N],b[N]; bool ok[N];
bool cmp(node u,node v){ return u.k<v.k; }
double getk(int x,int y){
	if (fabs(a[x].x-a[y].x)<eps) return inf;
	return (a[y].y-a[x].y)/(a[y].x-a[x].x);
}
double dis(int x,int y){
	return sqr(a[x].x-a[y].x)+sqr(a[x].y-a[y].y);
}
void solve(int l,int r){
	if (l==r) return;
	int mid=(l+r)>>1,i,j=l,k=mid+1,tp1=0,tp2=0;
	for (i=l; i<=r; i++)
		if (a[i].id<=mid) b[j++]=a[i]; else b[k++]=a[i];
	for (i=l; i<=r; i++) a[i]=b[i];
	solve(l,mid); j=1;
	for (i=l; i<=mid; i++) if (!a[i].op){
		while (tp1>1 && getk(q1[tp1-1],i)+eps>getk(q1[tp1-1],q1[tp1])) tp1--;
		q1[++tp1]=i;
		while (tp2>1 && getk(q2[tp2-1],i)<getk(q2[tp2-1],q2[tp2])+eps) tp2--;
		q2[++tp2]=i;
	}
	for (i=mid+1; i<=r; i++) if (a[i].op){
		if (a[i].y<0){
			while (tp1>1 && getk(q1[tp1-1],q1[tp1])<a[i].k) tp1--;
			if (tp1>0 && dis(q1[tp1],0)<dis(q1[tp1],i)) ok[a[i].p]=0;
		} else{
			while (j<tp2 && getk(q2[j],q2[j+1])<a[i].k) j++;
			if (j<=tp2 && dis(q2[j],0)<dis(q2[j],i)) ok[a[i].p]=0;
		}
	}
	solve(mid+1,r);
	for (i=j=l,k=mid+1; i<=r; i++)
		if (j<=mid && a[j].x<a[k].x || k>r) b[i]=a[j++]; else b[i]=a[k++];
	for (i=l; i<=r; i++) a[i]=b[i];
}
int main(){
	scanf("%d",&n); int i; bool flag=0;
	for (i=1; i<=n; i++){
		scanf("%d%lf%lf",&a[i].op,&a[i].x,&a[i].y);
		if (a[i].op){ a[i].p=++cnt; ok[cnt]=flag; } else flag=1;
		if (a[i].y) a[i].k=-a[i].x/a[i].y; else a[i].k=inf;
		a[i].id=i;
	}
	sort(a+1,a+n+1,cmp); solve(1,n);
	for (i=1; i<=cnt; i++) puts(ok[i]?"Yes":"No");
	return 0;
}

by lych

2016.4.16