TopCoderOpen 2020 North Asia : Hard EllysIncinerator

GCJを思い出す。
https://community.topcoder.com/stat?c=problem_statement&pm=16217&rd=18197

問題

2次元空間上にいくつか交わらない円が置かれている。
ここに、新たな円を1つ置きたい。
円の中心は(0,0)-(1000,1000)の矩形内になければならず、かつ既存の円と交わってはならない。
置ける最大の最大半径を求めよ。

解法

円の置き方の候補は以下のいずれかである。

• 2つの円に接する
• 1の円に接し、かつ中心が矩形の外周上

それぞれのケースにおいて置く円の半径を二分探索すると中心の位置が決まるので、既存の円と混じらわないかチェックしよう。

int N;
vector <double> X,Y,R;

class EllysIncinerator {
	public:
	
	int ok2(double cx,double cy,double r) {
		int i;
		if(0>cx||0>cy||cx>1000||cy>1000) return 0;
		
		FOR(i,N) {
			double r2=hypot(X[i]-cx,Y[i]-cy);
			if(r2<R[i]+r-(1e-9)) return 0;
		}
		return 1;
	}
	
	int ok(int a,int b,double D) {
		double r=hypot(X[a]-X[b],Y[a]-Y[b]);
		
		double tx,ty;
		
		if(r>R[a]+R[b]+2*D) {
			double dx=X[b]-X[a];
			double dy=Y[b]-Y[a];
			double dl=hypot(dx,dy);
			dx/=dl;
			dy/=dl;
			dx=X[a]+dx*(R[a]+D);
			dy=Y[a]+dy*(R[a]+D);
			if(ok2(dx,dy,D)) return 1;
		}
		else {
			double dx=X[b]-X[a];
			double dy=Y[b]-Y[a];
			double r1=R[a]+D;
			double r2=R[b]+D;
			double A=(dx*dx+dy*dy+r1*r1-r2*r2)/2;
			double d=(dx*dx+dy*dy)*r1*r1-A*A+1e-9;
			if(d<0) d=0;
			
			tx=X[a]+(A*dx+dy*sqrt(d))/(dx*dx+dy*dy);
			ty=Y[a]+(A*dy-dx*sqrt(d))/(dx*dx+dy*dy);
			if(ok2(tx,ty,D)) return 1;
			tx=X[a]+(A*dx-dy*sqrt(d))/(dx*dx+dy*dy);
			ty=Y[a]+(A*dy+dx*sqrt(d))/(dx*dx+dy*dy);
			if(ok2(tx,ty,D)) return 1;
			
			
			
		}
		return 0;
	}
	
	double ok3(int a,double D) {
		D+=R[a];
		if(X[a]<D) {
			double dy=sqrt(D*D-X[a]*X[a]);
			if(ok2(0,Y[a]+dy,D-R[a])) return 1;
			if(ok2(0,Y[a]-dy,D-R[a])) return 1;
		}
		if(X[a]+D>1000) {
			double dy=sqrt(D*D-(1000-X[a])*(1000-X[a]));
			if(ok2(1000,Y[a]+dy,D-R[a])) return 1;
			if(ok2(1000,Y[a]-dy,D-R[a])) return 1;
		}
		if(Y[a]<D) {
			double dx=sqrt(D*D-Y[a]*Y[a]);
			if(ok2(X[a]+dx,0,D-R[a])) return 1;
			if(ok2(X[a]-dx,0,D-R[a])) return 1;
		}
		if(Y[a]+D>1000) {
			double dx=sqrt(D*D-(1000-Y[a])*(1000-Y[a]));
			if(ok2(X[a]+dx,1000,D-R[a])) return 1;
			if(ok2(X[a]-dx,1000,D-R[a])) return 1;
		}
		return 0;
	}
	
	double getMax(vector <double> X, vector <double> Y, vector <double> R) {
		::X=X;
		::Y=Y;
		::R=R;
		N=X.size();
		
		int x,y,i;
		double ma=0;
		FOR(y,N) FOR(x,y) {
			double L=0,R=2000;
			FOR(i,60) {
				double M=(L+R)/2;
				if(ok(x,y,M)) L=M;
				else R=M;
			}
			ma=max(ma,L);
		}
		FOR(x,N) {
			double L=0,R=2000;
			FOR(i,60) {
				double M=(L+R)/2;
				if(ok3(x,M)) L=M;
				else R=M;
			}
			ma=max(ma,L);
		}

		return ma;
	}
}

まとめ

最初2つ目のケースを見落としてたので、出てたらレート落ちてたかな。