KUPC 2011 F: ボ〜ル - komiyamの日記

問題概要

原点を含まない円がN(<1500)個ある。この中から最大K個選んで、半直線OA(Aは第1、第2象限にある)と交わる確率が最大になるようにしたい。その確率を求める問題。

考えたこと

幾何は問題じゃなく、極座標で考えると区間の被覆する区間を最大化する問題に落ちることが分かる。
フローのF？でもノード数多いしそもそもフローだとしたらdoubleな時点で自分には解けない。
残るは貪欲かDPか。入力サイズ的にDPだろうなあ。
状態を(何番目まで見たか、今までいくつ選んだか、直前に選んだ区間は何か)としたら状態数O(N^2*K)で解けない。
本番中はここまで考えて、計算量が落とせなかったので別の問題に移った。
解法を人から教えてもらって再挑戦。
幾何パートはEPSとか使わずに処理する。doubleだからいつでもEPSを使う、というのは危険な思想だ。EPSとか基本的に最後の手段だし…。
DPパートはキレイに書けた。送信すると2ケースだけ合わない。
2ケースだけだから何かコーナーケースにハマった？よく分からない。
後日ジャッジデータが公開されてからローカルで試してみた。DPの更新式の最も最後の処理にミスが判明。訂正して無事AC。

後から通したコード

#include <cstdio>
#include <vector>
#include <algorithm>
#include <cmath>
#include <complex>
#include <cassert>
using namespace std;

const int MAX_N = 1509;
const double PI = atan(1.0)*4.0;
const pair<double, double> nothing = make_pair(0,0);

int N, K;

double dp[MAX_N][MAX_N + 1][2];
double xs[MAX_N], ys[MAX_N], rs[MAX_N];

int next_idx[MAX_N];
int cover_next[MAX_N];
vector< pair<double, double> > ranges;

//角度の区間へデータを変換する
void convert(){

	ranges.assign(N, nothing);
	for(int i=0; i<N; i++)if(ys[i] + rs[i] > 0){
		complex<double> center(xs[i], ys[i]);
		double r = rs[i];
		ranges[i].first = arg(center) - asin(r/abs(center));
		ranges[i].second = arg(center) + asin(r/abs(center));

		if( ranges[i].first < -PI ){
			ranges[i].first += 2*PI;
		}
		if( ranges[i].second > PI){
			ranges[i].second -= 2*PI;
		}

		if( xs[i] > 0 && ys[i] - rs[i] <= 0){
			ranges[i].first = 0.0;
		}
		if( xs[i] < 0 && ys[i] - rs[i] <= 0){
			ranges[i].second = PI;
		}
	}

}

//無駄な区間を減らしたり、次のindexを求めたり
void init(){
	convert();

	for(int i=0; i<N; i++){
		assert(ranges[i].first <= ranges[i].second);
//		assert(0 <= ranges[i].first && ranges[i].second <= PI);
	}


	vector< pair<double, double> > tmp;
	sort(ranges.begin(), ranges.end());
	for(int i=0; i<N; i++)if(ranges[i] != nothing){
		bool ok = true;
		for(int j=0; j<N; j++)if(i!=j){
			if(ranges[j].first <= ranges[i].first && ranges[i].second <= ranges[j].second){
				ok = false;
			}
		}
		if(ok){
			tmp.push_back(ranges[i]);
		}
		else{
			ranges[i] = nothing;
		}
	}

	ranges = tmp;
	N = ranges.size();
	K = min(K, N);

	//左端でソート
	sort(ranges.begin(), ranges.end());

	//被覆する中で一番得なもの
	fill(cover_next, cover_next+N, -1);
	for(int i=0; i<N-1; i++){
		double most_right = ranges[i].second;
		for(int j=i+1; j<N; j++){
			if(ranges[j].first <= ranges[i].second && most_right < ranges[j].second){
				most_right = ranges[j].second;
				cover_next[i] = j;
			}
		}
	}

	//被覆しない最初のもの
	fill(next_idx, next_idx+N, N);
	for(int i=0; i<N-1; i++){
		if(cover_next[i] == -1){
			next_idx[i] = i + 1;
		}
		else if(cover_next[i] < N-1){
			next_idx[i] = cover_next[i]+1;
		}
	}


	if(0)for(int i=0; i<N; i++){
		printf("[%f, %f], c_n:%d, ncn:%d\n", ranges[i].first, ranges[i].second, cover_next[i],next_idx[i]);
	}

	for(int i=0; i<N; i++){
		assert( ranges[i].second - ranges[i].first > 0);
	}

}


double solve(){
	init();
	if(N==0){
		return 0.0;
	}
	//dp[i][j][k] = {i番目まで見た。これまでにj個使った。
	//                 k=0ならi番目は使ってない、k=1ならi番目を使っている
	//k=0についてはi番目の区間はそれまでのとオーバーラップしない
	for(int i=0; i<N; i++){
		for(int j=0; j<=K; j++){

			int next = cover_next[i];
			int non_cover_next = next_idx[i];

			//何もせずにパス
			dp[i+1][j][0] = max(dp[i+1][j][0], dp[i][j][0]);
			dp[non_cover_next][j][0] = max(dp[non_cover_next][j][0], dp[i][j][1]);

			//新しく採用する
			if(j < K){
				dp[i][j+1][1] = max(dp[i][j+1][1], dp[i][j][0] + (ranges[i].second - ranges[i].first));
			}

			//重なってる部分があるのを伸ばす
			if(next != -1 && j < K){
				dp[next][j+1][1] = max(dp[next][j+1][1], dp[i][j][1] + (ranges[next].second - ranges[i].second));
			}
		}
	}

	double ans = 0.0;
	for(int i=0; i<=K; i++)for(int k=0; k<2; k++){
		ans = max(ans, dp[N][i][k]);
	}

	return ans/PI;
}


int main(){

	scanf("%d%d",&N,&K);
	for(int i=0; i<N; i++){
		scanf("%lf%lf%lf", xs+i, ys+i, rs+i);
	}

	printf("%.7f\n", solve());

	return 0;
}