kmjp's blog

競技プログラミング参加記です

AtCoder ARC #022 : D - スプリンクラー

ARC

実装がかなり面倒な問題。
http://arc022.contest.atcoder.jp/tasks/arc022_4

問題

2次元座標上N個の異なる格子点(X[i],Y[i])が与えられる。
各点から、原点を通るような円を書く。
いずれかの円に含まれる格子点の数を求めよ。

解法

以後説明のため、Pをx,yの絶対値の上限とする。

ケース1 : N<=100、P==100

各点が愚直に各円に含まれるか判定すればよい。
O(N*P^2)なので余裕で間に合う。

ケース2 : N<=1000、P==1000

自分が本番解けたのはここまで。
累積和を用いて一次元分処理を減らす。
各円について、Y座標ごとにX座標の最左点および最右点を求めておく。
後で各Y座標においていずれかの円に含まれる点を列挙する。
これでO(NP + P^2)なので間に合う。

ケース3 : N<=10^5、P==1000

ケース2に追加して、明らかに処理不要な円＝他の円に含まれる円の処理を省きたい。
処理すべき円は、与えられる格子点において凸包を成す点だけである。
凸包を成す点の数はO(P)程度なので、ケース2でO(NP)かかる処理がO(P^2)で済むようになる。

以前凸包を求めるライブラリを作ってあったので、凸包を使うアイデアを聞いたらケース2→3はすぐコードを書けた。

int N;
int X[100001],Y[100001];
int HH[5001][5001];
int SQ[2000001];
int hoge[2001][2001];


ll veccross(pair<int,int> p1,pair<int,int> p2,pair<int,int> p3) {
	p3.first-=p1.first;p2.first-=p1.first;
	p3.second-=p1.second;p2.second-=p1.second;
	return p3.first*(ll)p2.second-p2.first*(ll)p3.second;
}

vector<int> convex_hull(vector< pair<int, int> >& vp) {
	vector<pair<pair<int, int>, int> > sorted;
	vector<int> res;
	int i,k=0,rb;
	
	if(vp.size()<=2) {
		if(vp.size()>=1) res.push_back(0);
		if(vp.size()>=2) res.push_back(1);
		return res;
	}
	
	FOR(i,vp.size()) sorted.push_back(make_pair(vp[i],i));
	sort(sorted.begin(),sorted.end());
	res.resize(vp.size()*2);
	/* bottom */
	FOR(i,vp.size()) {
		while(k>1 && veccross(vp[res[k-2]],vp[res[k-1]],sorted[i].first)<=0) k--;
		res[k++]=sorted[i].second;
	}
	/* top */
	for(rb=k, i=vp.size()-2;i>=0;i--) {
		while(k>rb && veccross(vp[res[k-2]],vp[res[k-1]],sorted[i].first)<=0) k--;
		res[k++]=sorted[i].second;
	}
	res.resize(k-1);
	return res;
}


void solve() {
	int f,i,j,k,l,x,y;
	
	cin>>N;
	
	FOR(x,5001) FOR(y,5001) HH[x][y]=-10000;
	for(i=1;i<=2000000;i++) {
		SQ[i]=SQ[i-1];
		if((SQ[i]+1)*(SQ[i]+1)<=i) SQ[i]++;
	}
	
	FOR(i,N) cin>>X[i]>>Y[i];
	vector<pair<int,int> > V;
	
FOR(i,N) V.push_back(make_pair(X[i],Y[i]));
	vector<int> C=convex_hull(V);
	
	FOR(i,C.size()) {
		int cx=V[C[i]].first;
		int cy=V[C[i]].second;
		if(cx<-1000 || cx>1000) return;
		if(cy<-1000 || cy>1000) return;
		
		ll r=cx*cx+cy*cy;
		for(y=cy-(SQ[r]+1);y<=cy+(SQ[r]+1);y++) {
			if(r-(cy-y)*(cy-y)>=0) {
				x=SQ[r-(cy-y)*(cy-y)];
				HH[y+2500][cx-x+2500]=max(HH[y+2500][cx-x+2500],cx+x+2500);
			}
		}
	}
	
	ll ret=0;
	FOR(y,5001) {
		int prey=-5000,my=-5001;
		FOR(x,5001) {
			if(HH[y][x]>-10000) {
				if(x<=my) {
					ret+=x-prey;
					prey=x;
					my=max(my,HH[y][x]);
				}
				else {
					ret+=my-prey+1;
					prey=x;
					my=HH[y][x];
				}
			}
		}
		ret += my-prey+1;
	}
	cout << ret << endl;
}

ケース4 : N<=10^5、P==10^5

今度はO(P^2)も使えない。
ケース3では各円について全Y座標を対象としたが、実際は円の大半は他の円に覆いかぶさる。
よって、全部の円の和を取って生成される図形の辺に沿って処理すればよい。

そのため、隣接する円との交点を求め、他の円に覆いかぶさられないY座標のみ処理する。
図形の外周は各整数のY座標と数回ずつぐらいしか交差しないので、O(P)程度で処理できる。
凸包を求めるO(N*logN)の方が重いかな。

ll veccross(pair<int,int> p1,pair<int,int> p2,pair<int,int> p3) {
	p3.first-=p1.first;p2.first-=p1.first;
	p3.second-=p1.second;p2.second-=p1.second;
	return p3.first*(ll)p2.second-p2.first*(ll)p3.second;
}

vector<int> convex_hull(vector< pair<int, int> >& vp) {
	vector<pair<pair<int, int>, int> > sorted;
	vector<int> res;
	int i,k=0,rb;
	
	if(vp.size()<=2) {
		if(vp.size()>=1) res.push_back(0);
		if(vp.size()>=2) res.push_back(1);
		return res;
	}
	
	FOR(i,vp.size()) sorted.push_back(make_pair(vp[i],i));
	sort(sorted.begin(),sorted.end());
	res.resize(vp.size()*2);
	/* bottom */
	FOR(i,vp.size()) {
		while(k>1 && veccross(vp[res[k-2]],vp[res[k-1]],sorted[i].first)<=0) k--;
		res[k++]=sorted[i].second;
	}
	/* top */
	for(rb=k, i=vp.size()-2;i>=0;i--) {
		while(k>rb && veccross(vp[res[k-2]],vp[res[k-1]],sorted[i].first)<=0) k--;
		res[k++]=sorted[i].second;
	}
	res.resize(k-1);
	return res;
}

int N;
vector<pair<int,int> > V;
vector<pair<int,int> > cond[500002];


void solve() {
	int f,i,j,k,l;
	ll x,y;
	
	cin>>N;
	
	FOR(i,N) {
		cin>>x>>y;
		V.push_back(make_pair(x,y));
	}
	vector<int> C=convex_hull(V);
	vector<complex<double> > CC;
	
	FOR(i,C.size()) {
		int cx1=V[C[i]].first;
		int cy1=V[C[i]].second;
		int cx2=V[C[(i+1)%C.size()]].first;
		int cy2=V[C[(i+1)%C.size()]].second;
		complex<double> c1(cx1,cy1), c2(cx2,cy2);
		
		complex<double> v1=c2-c1,v2=-c1,v3;
		double ab=abs(v1);
		v1.real()/=ab;
		v1.imag()/=ab;
		double dot=v1.real()*v2.real()+v1.imag()*v2.imag();
		v3=v2-dot*v1;
		if(v1.real()*v2.imag()-v2.real()*v1.imag()>0) v3=-v3;
		v2=dot*v1-v3+c1;
		CC.push_back(v2);
	}
	
	FOR(i,C.size()) {
		ll cx1=V[C[i]].first;
		ll cy1=V[C[i]].second;
		complex<double> c(cx1,cy1);
		complex<double> p1=CC[(i+C.size()-1)%C.size()]-c;
		complex<double> p2=CC[i]-c;
		
		ll r=cx1*cx1+cy1*cy1;
		if(p1==p2) {
			for(y=-sqrt(r)-1;y<=sqrt(r)+1;y++) if(r-y*y>=0) {
				x=sqrt(r-y*y);
				if((x+1)*(x+1)+y*y<=r) x++;
				cond[250000+cy1+y].push_back(make_pair(cx1-x,1000000+i));
				cond[250000+cy1+y].push_back(make_pair(cx1+x+1,i));
			}
		}
		else if(p1.real()<=0 && p2.real()<=0) {
			if(p1.imag()<=p2.imag()) {
				for(y=ceil(p1.imag());y<=p2.imag();y++) if(r-y*y>=0) {
					x=sqrt(r-y*y);
					if((x+1)*(x+1)+y*y<=r) x++;
					cond[250000+cy1+y].push_back(make_pair(cx1-x,1000000+i));
					if(x==0) cond[250000+cy1+y].push_back(make_pair(cx1+x+1,0+i));
				}
			}
			else {
				for(y=-sqrt(r)-1;y<=sqrt(r)+1;y++) if(r-y*y>=0) {
					x=sqrt(r-y*y);
					if((x+1)*(x+1)+y*y<=r) x++;
					cond[250000+cy1+y].push_back(make_pair(cx1+x+1,i));
					if(y>=p2.imag() || y<=p1.imag())
						cond[250000+cy1+y].push_back(make_pair(cx1-x,1000000+i));
				}
			}
		}
		else if(p1.real()>0 && p2.real()<0) {
			for(y=floor(p1.imag());y>=-sqrt(r)-1;y--) if(r-y*y>=0) {
				x=sqrt(r-y*y);
				if((x+1)*(x+1)+y*y<=r) x++;
				cond[250000+cy1+y].push_back(make_pair(cx1+x+1,i));
			}
			for(y=-sqrt(r)-1;y<=p2.imag();y++) if(r-y*y>=0) {
				x=sqrt(r-y*y);
				if((x+1)*(x+1)+y*y<=r) x++;
				cond[250000+cy1+y].push_back(make_pair(cx1-x,1000000+i));
			}
		}
		else if(p1.real()<0 && p2.real()>0) {
			for(y=ceil(p1.imag());y<=sqrt(r)+1;y++) if(r-y*y>=0) {
				x=sqrt(r-y*y);
				if((x+1)*(x+1)+y*y<=r) x++;
				cond[250000+cy1+y].push_back(make_pair(cx1-x,1000000+i));
			}
			for(y=sqrt(r)+1;y>=p2.imag();y--) if(r-y*y>=0) {
				x=sqrt(r-y*y);
				if((x+1)*(x+1)+y*y<=r) x++;
				cond[250000+cy1+y].push_back(make_pair(cx1+x+1,i));
			}
		}
		else if(p1.real()>0 && p2.real()>0) {
			if(p1.imag()>p2.imag()) {
				for(y=floor(p1.imag());y>=p2.imag();y--) if(r-y*y>=0) {
					x=sqrt(r-y*y);
					if((x+1)*(x+1)+y*y<=r) x++;
					cond[250000+cy1+y].push_back(make_pair(cx1+x+1,i));
				}
			}
			else {
				for(y=-sqrt(r)-1;y<=sqrt(r)+1;y++) if(r-y*y>=0) {
					x=sqrt(r-y*y);
					if((x+1)*(x+1)+y*y<=r) x++;
					cond[250000+cy1+y].push_back(make_pair(cx1-x,1000000+i));
					if(y<=p2.imag() || y>=p1.imag())
						cond[250000+cy1+y].push_back(make_pair(cx1+x+1,i));
				}
			}
		}
	}
	
	ll ret=0;
	FOR(i,500001) if(!cond[i].empty()) {
		sort(cond[i].begin(),cond[i].end());
		int pre=-1000000,st=0;
		FOR(j,cond[i].size()) {
			if(j>0 && (cond[i][j].first==cond[i][j-1].first) && (cond[i][j].second/1000000==cond[i][j-1].second/1000000)) continue;
			if(cond[i][j].second>=1000000) {
				if(st==0) pre=cond[i][j].first;
				st++;
			}
			else {
				if(--st==0) ret+=cond[i][j].first-pre;
			}
		}
	}
	
	cout << ret << endl;
}

まとめ

各円のうち、外周を形成する部分だけ抜き出す処理が条件分岐だらけで汚い。
もう少しスマートに書けないものか。