NEFU709(第K个圆的半径)

题目：The Kth Circle

 

题意：平面上给定n个点的坐标，我们可以选其中3点构成一个圆，当然也就有3段弧，我们定义第K个圆满足条件：

(1)除了已选的3个点，剩下的n-3个点中有K个点在圆内，有n-K-1个点在圆外

(2)保证所有的点在圆的三条弧中其中一条弧的一侧

解析：本题实际上就是得到这样一个结论，满足条件的圆其中有两点一定在这些点的凸包上，并且这两点相邻，这样我们就可以先

求凸包，然后枚举凸包的每条边，设边的端点为A，B，再枚举点B，然后我们把角APB按照从小到大的顺序排序，其中第K+1个就

是我们要求的，对于枚举凸包的每一条边都求出第K+1个点与凸包边的端点形成的圆的半径，记录最小的。

由于点的数目2000，比较大，直接枚举会超时。所以采用堆优化。

 

#include <iostream>
#include <algorithm>
#include <stdio.h>
#include <cmath>

using namespace std;
const int N=2005;

struct Point
{
	int x,y;
};

struct HEAP
{
	double var;
	int id;
};

HEAP heap[N];//小顶堆
Point p[N];

int stack[N];

bool readIn(int &N,int &K)
{
	if(scanf("%d%d",&N,&K)==EOF)
		return false;
	for(int i=1;i<=N;i++)
		scanf("%d%d",&p[i].x,&p[i].y);
	return true;
}

int cross(Point A,Point B,Point C)
{
    return (B.x-A.x)*(C.y-A.y)-(B.y-A.y)*(C.x-A.x);
}

bool cmp(const Point &a,const Point &b)
{
	return cross(p[1],a,b)>0 ;
}

double dist(Point A,Point B)
{
    return sqrt((A.x-B.x)*(A.x-B.x)+(A.y-B.y)*(A.y-B.y));
}

double getAng(const Point &a,const Point &b,const Point &c)
{
	double A,B,C;
	C=dist(b,c);
	A=dist(a,b);
	B=dist(a,c);
	double ang_a=acos((A*A+B*B-C*C)/(2*A*B));
	return ang_a;
}

/**根据角a和b对点排序*/
bool heapCmp(const HEAP &a,const HEAP &b)
{
	return a.var>b.var;
}

void insert(const HEAP &a,int &cnt)
{
	heap[++cnt]=a;
	int  j=cnt;
	while(j>1)
	{
		if(heapCmp(heap[j/2],heap[j]))
		{
			swap(heap[j/2],heap[j]);
			j/=2;
		}
		else break;
	}
}

void heapDown(const int &id,const int &cnt)
{
	int lc=id*2,rc=id*2+1,j=id;
	if(lc<=cnt&&heapCmp(heap[j],heap[lc])) j=lc;
	if(rc<=cnt&&heapCmp(heap[j],heap[rc])) j=rc;
	if(j!=id)
	{
		swap(heap[id],heap[j]);
		heapDown(j,cnt);
	}
	return;
}

/**已知三点求三角形的面积*/
double triangleArea(const Point &a,const Point &b,const Point &c)
{
	return fabs(cross(a,b,c))/2.0;
}

/**已知圆上三点求圆半径*/
double getRadius(const Point &a,const Point &b,const Point &c)
{
	double A=dist(a,b);
	double B=dist(a,c);
	double C=dist(b,c);
	return (A*B*C)/(triangleArea(a,b,c)*4.0);
}

void Solve(const int &N,const int &K)
{
	int id=1;
	for(int j=1;j<=N;j++)
		if(p[id].x<p[j].x||(p[id].x==p[j].x&&p[id].y<p[j].y))
			id=j;
	swap(p[1],p[id]);
	sort(p+2,p+N+1,cmp);
	int top=0;
	stack[++top]=1,stack[++top]=2;
	for(int i=3;i<=N;i++)
	{
		while(top>1&&cross(p[stack[top-1]],p[stack[top]],p[i])<0) top--;
		stack[++top]=i;
	}
	stack[++top]=stack[1];
	double ans=1000000000;
	HEAP tmp;
	for(int i=1;i<top;i++)
	{
		int cnt=0;
		for(int j=1;j<=N;j++)
		{
			if(j!=stack[i]&&j!=stack[i+1])
			{
				tmp.id=j;
				tmp.var=getAng(p[j],p[stack[i]],p[stack[i+1]]);
				if(cnt<K) insert(tmp,cnt);
				else
				{
					if(heapCmp(tmp,heap[1]))
					{
						heap[1]=tmp;
						heapDown(1,cnt);
					}
				}
			}
		}
		double tmp=getRadius(p[stack[i]],p[stack[i+1]],p[heap[1].id]);
		ans=ans>tmp? tmp:ans;
	}
	printf("%.2lf\n",ans);
}

int main()
{
	int N,K;
	while(readIn(N,K))
	{
		K++;
		Solve(N,K);
	}
	return 0;
}