POJ 1981 Circle and Points(n3 & n2logn)解法与疑点解析

Detail
题目大意：
给出平面上n个点求单位圆的最大覆盖点数

朴素做法O(n3)

n2枚举两个点+单位半径确定（2个）圆，任取其一方向构造唯一圆，由数学关系求出圆心坐标，最后枚举点到圆心距离是否大于1即可。代码如下：

#include <iostream>
#include <cmath>
#include <algorithm>
#include <cstring>
#include <cstdio>
#define In inline
#define pi atan(1.0)*4
#define enter puts("")
#define rep(i,a,b) for(int i=(a);i<=(b);i++)
#define eps 1e-6
using namespace std;

In ll read()
{
	ll ans = 0;
	char ch = getchar(), las = ' ';
	while(!isdigit(ch)) las = ch, ch = getchar();
	while(isdigit(ch)) ans = (ans << 1) + (ans << 3) + ch - '0', ch = getchar();
	if(las == '-')	ans = -ans;
	return ans;
}
In void write(ll x)
{
	if(x < 0) x = -x, putchar('-');
	if(x >= 10)	write(x / 10);
	putchar(x%10 + '0');
}
int sgn(double x){
	if(fabs(x) < eps) return 0;
	return x < 0 ? -1 : 1;
}

struct Point{
	double x, y;
	Point(){}
	Point(double x, double y):x(x), y(y){}
	void input(){scanf("%lf%lf", &x, &y);}
	Point operator += (Point b){return Point(x+b.x, y+b.y);}
	Point operator + (Point b){return Point(x+b.x, y+b.y);}
	Point operator - (Point b){return Point(x-b.x, y-b.y);}
	Point operator * (double k){return Point(x*k, y*k);}
	Point operator / (double k){return Point(x/k, y/k);}
	double operator ^ (Point b){return x*b.y-y*b.x;}
	bool operator == (Point b){
		return sgn(x-b.x) == 0 && sgn(y-b.y) == 0;
	}
	bool operator < (Point b){
		return sgn(x-b.x) == 0 ? sgn(y-b.y)<0 : x < b.x;
	}
	void norm(){
		double dis = sqrt(x*x+y*y);
		x/=dis, y/=dis;
	}
	double distance(Point b){
		return sqrt((x-b.x)*(x-b.x)+(y-b.y)*(y-b.y));
	}
	Point rotate(int k){// 1 :left  2: right
		if(k == 1)	return Point{y, -x};
		else return Point{-y, x};
	}	
};
typedef Point Vector;
const int N = 500;
Point p[N];
	 
void run_case(){
	int n;
	while(~scanf("%d", &n), n){
		int ans = 1;
	 	rep(i, 0, n-1) p[i].input();
	 	rep(i, 0, n-2) rep(j, i+1, n-1){
	 		double dis = p[i].distance(p[j]);
	 		if(sgn(dis-2)>0) continue;
	 		int cnt = 0;
	 		Vector ab = p[j]-p[i];
	 		Point m = (p[i]+p[j])/2, mid;
	 		ab.norm();
	 		Vector oR = ab.rotate(1);
	 		/*mid = m + oR*sqrt(1-dis*dis/4);
	 		rep(k, 0, n-1)	if(sgn(mid.distance(p[k])-1)<=0) cnt ++;
	 		ans = max(ans, cnt);*/
	 		cnt = 0;
	 		oR = ab.rotate(2);
	 		mid = m + oR*sqrt(1-dis*dis/4);
	 		rep(k, 0, n-1)	if(sgn(mid.distance(p[k])-1)<=0) cnt ++;
	 		ans = max(ans, cnt);
	 	}
	 	write(ans);
	 	enter;
	 }		
}

int main()
{
//	freopen("C:\\Users\\MARX HE\\Desktop\\input.txt","r",stdin);
//	freopen("C:\\Users\\MARX HE\\Desktop\\output.txt","w",stdout);
	int _ = 1;
//    int _ = read();
    while(_ --){
		run_case();
    }
}

接下来是这个解法疑惑的思考：
可以注意到代码注释部分其实是因为圆上两个点加已知半径只能确定两个圆(当然如果已知两点在直径上圆是唯一确定的）。
而加上两个圆一起判之后代码很不幸的T了

然鹅注释掉任一边的判断之后就过了，os我才跑了5e7就不行了，看看人家cf神机1s 1e8(
言归正传；直接开始考虑我们认为只判一边出现漏解的情况，
如图：

当考虑ab为已知点且选择顺时针方向圆心时，考虑c漏解，则当ac为已知点时b将被考虑，可以证明由此情况可以推知只考虑单方向圆心不会漏解。

求圆弧交（O(n2logn))

类似扫描线算法，枚举每一个点为圆心，再枚举每一个点与坐标圆的交点，这里选[-PI, PI]，规定小角度为1，大角度为-1.
排序后枚举最大值即可。

In ll read()
{
	ll ans = 0;
	char ch = getchar(), las = ' ';
	while(!isdigit(ch)) las = ch, ch = getchar();
	while(isdigit(ch)) ans = (ans << 1) + (ans << 3) + ch - '0', ch = getchar();
	if(las == '-')	ans = -ans;
	return ans;
}
In void write(ll x)
{
	if(x < 0) x = -x, putchar('-');
	if(x >= 10)	write(x / 10);
	putchar(x%10 + '0');
}
int sgn(double x){
	if(fabs(x) < eps) return 0;
	return x < 0 ? -1 : 1;
}

struct Point{
	double x, y;
	Point(){}
	Point(double x, double y):x(x), y(y){}
	void input(){scanf("%lf%lf", &x, &y);}
	Point operator += (Point b){return Point(x+b.x, y+b.y);}
	Point operator + (Point b){return Point(x+b.x, y+b.y);}
	Point operator - (Point b){return Point(x-b.x, y-b.y);}
	Point operator * (double k){return Point(x*k, y*k);}
	Point operator / (double k){return Point(x/k, y/k);}
	double operator ^ (Point b){return x*b.y-y*b.x;}
	bool operator == (Point b){
		return sgn(x-b.x) == 0 && sgn(y-b.y) == 0;
	}
	bool operator < (Point b){
		return sgn(x-b.x) == 0 ? sgn(y-b.y)<0 : x < b.x;
	}
	double distance(Point b){
		return sqrt((x-b.x)*(x-b.x)+(y-b.y)*(y-b.y));
	}	
};
const int N = 500;
Point p[N];
struct A{
	double ang;
	int ct;
	bool friend operator < (const A &a, const A &b){
		return a.ang < b.ang;
	}
}Ang[N*2];
void run_case(){
	int n;
	while(~scanf("%d", &n), n){
		int ans = 1;
		rep(i, 1, n) p[i].input();
		rep(i, 1, n){
			int cnt = 0;
			rep(j, 1, n){
				if(i == j) continue;
				if(p[i].distance(p[j]) > 2) continue;
				double angle = atan2(p[j].y-p[i].y,p[j].x-p[i].x);
				Point mid = (p[i] + p[j])/2;
				double ag = acos(mid.distance(p[i]));
				Ang[cnt].ang = angle - ag;
				Ang[cnt ++].ct = 1;
				Ang[cnt].ang = angle + ag;
				Ang[cnt ++].ct = -1;
			}
			sort(Ang, Ang+cnt);
			int temp = 1;
			rep(i, 0, cnt-1){
				temp += Ang[i].ct;
				ans = max(ans, temp);
			}
		}
		write(ans);
		enter;
	}					
}

int main()
{
//	freopen("C:\\Users\\MARX HE\\Desktop\\input.txt","r",stdin);
//	freopen("C:\\Users\\MARX HE\\Desktop\\output.txt","w",stdout);
	int _ = 1;
//    int _ = read();
    while(_ --){
		run_case();
    }
}

这里有一个细节问题可以思考一下，temp=1是圆心含有一个点。
可以造出样例使最后的Ang数组中的ang是在[-PI, Pi]之外的，关于这部分有两种处理方法：
1.
任其溢出，思考方式和上一种解法类似，简单说是对于a为圆心时如果b的作用弧度为[-189,-171], c的作用弧度为[172,188]。则此时枚举点a取得的极大值为2，然鹅他的最值不应该是3么？程序有错吗？显然没有（（（
因为如果b ,c 处于a的x轴分界线两侧，那么在枚举b/c时至少有一个点时其余两点不在分界线（负半轴）处，所以最大值枚举成立。
2.
但如果细心的小伙伴仍想精细处理该边界问题，则可以在枚举angle+ag时增加判断，如果越过PI则在Ang[cnt ++].PI -1, Ang[cnt].-PI.+1。这样排序后就可以做到不漏啦。
于是问题就这么愉快地解决了（梦中解决系列++

这么详细的解释就不需要附上代码了吧，打字真累…