crazy_cloud's Blog
0%

Codeforces 744D Hongcow Draws a Circle

Posted on Edited on In

题目大意

在二维平面上有 \(n\) 个红点和 \(m\) 个蓝点。

你要画一个圆，使得这个圆内部不含任何蓝点，且至少包含一个红点。边界上的点可以算也可以不算。

你要计算圆的半径最大能是多少。

\(1\leq n,m\leq10^3,1\leq x_i,y_i\leq10^4\)

没有两点重合。

题目分析

首先这个圆一定会经过至少一个蓝点。接下来有二分答案和圆的反演两种思路。

算法一

考虑二分答案。

你可能会注意到这题直接看似乎并不满足二分性。但是如果我们先枚举一个红点作为一定在圆上的点二分答案（这个肯定满足二分性），然后再枚举蓝点二分答案，我们可以证明后面枚举的蓝点答案要么小于当前最优值，要么满足二分性，也就是说我们可以先判断一遍当前答案在这个点是否可行，如果可行就二分答案。

至于怎么判定呢？确定了半径，圆心其实就是在这个半径的一个圆上旋转，我们可以计算出其它所有点进入和退出圆的角度，然后极角排序扫描线就好了。

时间复杂度 \(O(n^2\log n\log X)\)。

算法二

上面的算法复杂度比较大，这里我们可以使用经典套路使复杂度降下来。我们将点集 random_shuffle 一下，然后期望二分次数就变成了 \(O(\log n)\) 次。

时间复杂度变为了 \(O(n^2\log n+n\log n\log X)\)。

算法三

考虑圆的反演，我们枚举圆经过的蓝点并将其作为反演中心，选一个合适的反演半径，然后所有经过这个点的圆都变成了一条直线。

显然，直线离中心越近，原平面中圆的半径越大。也就是说我们要在反演过后的平面中找到一条离反演中心最近的直线使其远离反演中心的一侧半平面内没有蓝点，且至少包含一个红点。

对蓝点做凸包，然后求出红点在上面的切线，扫描线即可。

时间复杂度 \(O(n^2\log n)\)。

算法四？

出题人给出了一个未经证明的观察结论：反演后的蓝凸包上的点，一定是原平面中蓝点集 Delaunay 三角剖分中和反演中心蓝点有边相连的点。

也就是说，整个算法的过程中，其实凸包上点集的大小和是线性的。那么算法三的最后面就可以暴力计算了。

时间复杂度 \(O(n^2\log n)\)。

代码实现

我辣么菜当然是写算法二啦！

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
#include <algorithm>
#include <iostream>
#include <cstdio>
#include <cctype>
#include <cmath>

using namespace std;

typedef long double db;

inline int read()
{
	int x=0,f=1;
	char ch=getchar();
	while (!isdigit(ch)) f=ch=='-'?-1:f,ch=getchar();
	while (isdigit(ch)) x=x*10+ch-'0',ch=getchar();
	return x*f;
}

const int N=1005;
const int V=N<<1;
const db EPS=1e-4;
const db INF=1e+9;
const db pi=acosl(-1);
const int ERROR=-100000;

inline bool equ(db x,db y){return fabs(x-y)<=EPS;}
inline int sgn(db x){return equ(x,0.)?0:(x<0?-1:1);}

struct P
{
	db x,y;

	inline P (db x_=0.,db y_=0.){x=x_,y=y_;}

	inline P operator+(P const p)const{return P(x+p.x,y+p.y);}
	inline P operator-(P const p)const{return P(x-p.x,y-p.y);}
	inline P operator*(db const k)const{return P(x*k,y*k);}

	inline db operator*(P const p)const{return x*p.x+y*p.y;}
	inline db operator^(P const p)const{return x*p.y-y*p.x;}

	inline db mod2()const{return (*this)*(*this);}
	inline db mod()const{return sqrt(fabs(mod2()));}

	inline db arg(){return atan2(y,x);}
}pts[V];

int n,m;
db ans;

typedef pair<db,int> event;
#define mkp(a,b) make_pair(a,b)
#define ft first
#define sd second

event evt[V<<1];
db dist[V][V];

inline bool cmp(event p,event q){return p.ft<q.ft;}

inline void adjust(db &x){x=x>pi?x-pi*2:x<-pi?x+pi*2:x;}

inline bool judge(int id,db r)
{
	int cur=id<=n,cnt=0,ret=0;
	for (int i=1;i<=n+m;++i)
		if (i!=id)
		{
			db d=dist[i][id];
			if (sgn(d-2.*r)>0) continue;
			db alpha=acosl(d/(2.*r))-1e-30*(i>n),mid=(pts[i]-pts[id]).arg();
			db st=mid-alpha,en=mid+alpha;
			adjust(st),adjust(en);
			int cst=i<=n?1:ERROR;
			cur+=(sgn(en-st)<0)*cst,evt[++cnt]=mkp(st,cst),evt[++cnt]=mkp(en,-cst);
		}
	ret=cur,sort(evt+1,evt+1+cnt,cmp);
	for (int i=1;i<=cnt;++i) ret=max(ret,cur+=evt[i].sd);
	return ret>=1;
}

inline void binary_search(int id)
{
	db l=ans,r=ans;
	for (;judge(id,r)&&r<=INF;r*=2);
	for (db mid;sgn(r-l)>=0;)
	{
		mid=(l+r)*.5;
		if (judge(id,mid)) l=(ans=mid)+EPS;
		else r=mid-EPS;
	}
}

int main()
{
	freopen("circle.in","r",stdin),freopen("circle.out","w",stdout);
	n=read(),m=read();
	for (int i=1,x,y;i<=n+m;++i) x=read(),y=read(),pts[i]=P(x,y);
	random_shuffle(pts+1,pts+1+n),random_shuffle(pts+1+n,pts+1+n+m);
	for (int i=1;i<=n+m;++i)
		for (int j=1;j<=n+m;++j)
			dist[i][j]=(pts[i]-pts[j]).mod();
	ans=1.;
	for (int i=1;i<=n;++i)
	{
		if (!judge(i,ans)) continue;
		binary_search(i);
	}
	for (int i=n+1;i<=n+m;++i)
	{
		if (!judge(i,ans)) continue;
		binary_search(i);
	}
	if (ans<INF*.5) printf("%.10lf\n",(double)ans);
	else printf("-1\n");
	fclose(stdin),fclose(stdout);
	return 0;
}