LOJ #2159. 「POI2011 R1」收缩点 Plot（最小圆覆盖，倍增代替二分）

最新推荐文章于 2025-03-16 17:18:06 发布

Freopen

最新推荐文章于 2025-03-16 17:18:06 发布

阅读量237

点赞数

分类专栏：计算几何倍增

本文链接：https://blog.csdn.net/qq_35950004/article/details/103706497

版权

计算几何同时被 2 个专栏收录

29 篇文章

订阅专栏

倍增

15 篇文章

订阅专栏

该博客详细介绍了如何解决LOJ #2159问题，即「POI2011 R1」收缩点 Plot。博主指出在寻找最小圆覆盖时，不能简单地使用二分查找，因为检查操作的复杂度较高。为了解决这个问题，博主提出了使用倍增算法来替代二分，这样总复杂度变为O(答案logn)。文章着重强调了这种方法在处理精度问题上的挑战，并提供了通过验证的代码。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

题目

首先可以二分答案 $a n s$ 。
然后再从左往右在最小圆半径 $< = a n s$ 的情况下尽量多分点在同一组。
因为最小圆覆盖需要随机，所以不能直接从左往右加。
这个过程也不能二分，因为 $c h e c k$ 的复杂度是 $O （二分答案）$ 的。
但是发现我们有总答案 $= n$
于是可以用倍增代替二分，那么我们 $c h e c k$ 的总复杂度是 $O(答案\log n)$ 的。
注意这个输出方案极度卡精。

$\mathrm{AC \ Code}$

#include<bits/stdc++.h>
#define maxn 100005
#define Pt Point
#define Ct const 
#define db double
#define rep(i,j,k) for(int i=(j);i<=(k);i++)
#define eps 1e-20
using namespace std;

int n,m;

inline db sqr(Ct db &x){ return x*x; }
struct Pt{
	db x,y;Pt(db x=0,db y=0):x(x),y(y){}
	Pt operator +(Ct Pt &B)Ct{ return Pt(x+B.x,y+B.y); }
	Pt operator -(Ct Pt &B)Ct{ return Pt(x-B.x,y-B.y); }
	Pt operator *(Ct db &B)Ct{ return Pt(x*B,y*B); }
	Pt operator /(Ct db &B)Ct{ return Pt(x/B,y/B); }
	db operator *(Ct Pt &B)Ct{ return x*B.y-y*B.x; }
	db Len()Ct{ return sqr(x)+sqr(y); }
	Pt Rotate90(){ return Pt(y,-x); }
}P[maxn],O[maxn],rP[maxn];
int cnt;db r[maxn];

Pt Intersec(Ct Pt &p1,Ct Pt &v1,Ct Pt &p2,Ct Pt &v2){
	db t = (p2 - p1) * v2 / (v1 * v2);
	return p1 + (v1 * t);
}

db cal(int L,int R,Pt &O,db &r){
	random_shuffle(&P[L],&P[R+1]);
	r=0;
	rep(i,L,R) if((P[i]-O).Len()>r){
		O=P[i],r=0;
		rep(j,L,i-1) if((P[j]-O).Len()>r){
			O=(P[i]+P[j])/2,r=(P[i]-P[j]).Len()/4;
			rep(k,L,j-1) if((P[k]-O).Len()>r){
				if(fabs((P[i]-P[j])*(P[i]-P[k]))>1e-8)
					O=Intersec((P[i]+P[j])/2,(P[i]-P[j]).Rotate90(),(P[i]+P[k])/2,(P[i]-P[k]).Rotate90());
				else{
					db dij = (P[i]-P[j]).Len() , dik = (P[i] - P[k]).Len() , djk = (P[j] - P[k]).Len();
					if(dij > dik && dij > djk) O=(P[i]+P[j])/2;
					if(dik > dij && dik > djk) O=(P[i]+P[k])/2;
					if(djk > dij && djk > dik) O=(P[j]+P[k])/2;
				}
				r=(P[i]-O).Len();
			}
		}
	}
	rep(i,L,R) P[i]=rP[i];
	return r;
}

bool check(db d){
	cnt = 0;	
	for(int i=1,j,p;i<=n;i=i+p){
		if(++cnt > m) return 0;
		j=2,p=0;
		for(;i+j-1<=n;j<<=1) if(cal(i,i+j-1,O[cnt],r[cnt]) > d) break;
		for(int k=j>>1;k;k>>=1) if(i+p+k-1<=n && cal(i,i+p+k-1,O[cnt],r[cnt]) <= d) p+=k;
		cal(i,i+p-1,O[cnt],r[cnt]);
	}
	return 1;
}

int main(){
	scanf("%d%d",&n,&m);
	rep(i,1,n) scanf("%lf%lf",&P[i].x,&P[i].y),rP[i]=P[i];
	db L=0,R=1e14,mid;
	for(int stp=0;stp<=100;stp++){
		mid=(L+R)*0.5;
		if(check(mid)) R=mid;
		else L=mid;
		if(R-L < eps) break;
	}
	mid=(L+R)*0.5;
	printf("%.15lf\n",sqrt(mid));
	check(mid+1e-3);
	printf("%d\n",cnt);
	rep(i,1,cnt) printf("%.15lf %.15lf\n",O[i].x,O[i].y);
}