Problem Description
Alice is interesting in computation geometry problem recently. She found a interesting problem and solved it easily. Now she will give this problem to you :
You are given N distinct points (Xi,Yi) on the two-dimensional plane. Your task is to find a point P and a real number R, such that for at least ⌈N2⌉ given points, their distance to point P is equal to R.
Input
The first line is the number of test cases.
For each test case, the first line contains one positive number N(1≤N≤105).
The following N lines describe the points. Each line contains two real numbers Xi and Yi (0≤|Xi|,|Yi|≤103) indicating one give point. It’s guaranteed that N points are distinct.
Output
For each test case, output a single line with three real numbers XP,YP,R, where (XP,YP) is the coordinate of required point P. Three real numbers you output should satisfy 0≤|XP|,|YP|,R≤109.
It is guaranteed that there exists at least one solution satisfying all conditions. And if there are different solutions, print any one of them. The judge will regard two point’s distance as R if it is within an absolute error of 10−3 of R.
Sample Input
1
7
1 1
1 0
1 -1
0 1
-1 1
0 -1
-1 0
Sample Output
0 0 1
题目保证有解,若n>=3,则解一定是不共线3点的外接圆
随机选择3点,判断是否符合条件
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=5e4+7;
const int maxp=1e5+7;
const double inf=1e100;
const double eps=1e-8;
const double pi =acos(-1.0);
int sgn(double x){
if(fabs(x)<eps)return 0;
return x>0?1:-1;
}
//double hypot(double x,double y){return sqrt(x*x+y*y);}
struct Point{
double x,y;
Point(double x=0,double y=0):x(x),y(y){}
bool operator ==(const Point& b){return !sgn(x-b.x)&&!sgn(y-b.y);}
bool operator <(const Point& b)const{return !sgn(x-b.x)?y<b.y:x<b.x;}
double operator^(const Point& b){return x*b.y-y*b.x;}//叉积
double operator*(const Point& b){return x*b.x+y*b.y;}//点积
Point operator +(const Point &b){return Point(x+b.x,y+b.y);}
Point operator -(const Point& b){return Point(x-b.x,y-b.y);}
Point operator *(const double &k){return Point(x*k,y*k);}
Point operator /(const double &k){return Point(x/k,y/k);}
double len(){return hypot(x,y);}
double len2(){return x*x+y*y;}
double rad(Point a,Point b){Point p=*this;return fabs(atan2(fabs((a-p)^(b-p)),(a-p)*(b-p)));}//pa和pb夹角
double angle(){return atan2(y,x);}//倾斜角
double angle(Point B){//夹角
Point A=(*this);
return acos(A*B/A.len()/B.len());
}
double distance(Point b){return hypot(x-b.x,y-b.y);}
Point trunc(double r){//化为长度为r的向量
double l=len();
if(!sgn(l))return *this;
r/=l;return Point(x*r,y*r);
}
Point rotleft(){return Point(-y,x);}//逆时针旋转90 度
Point rotright(){return Point(y,-x);}//顺时针旋转90 度
Point rotate(double rad){//rad为弧度 逆时针旋转
Point A=(*this);
double c=cos(rad),s=sin(rad);
return Point(A.x*c-A.y*s,A.x*s+A.y*c);
}
Point rotate(Point p,double angle){//绕p逆时针旋转后点
Point v=((*this)-p).rotate(angle);
return Point(p.x+v.x,p.y+v.y);
}
void input(){scanf("%lf%lf",&x,&y);}
void show(){printf("(%.2lf %.2lf)\n",x,y);}
};
struct Line{
Point s,e;
Line(){};
Line(Point s,Point e):s(s),e(e){};
void adjust(){swap(s,e);}
int relation(Point p){//点和直线关系1在左侧 2在右侧 3在直线上
int c=sgn((p-s)^(e-s));
if(c<0)return 1;
if(c>0)return 2;
return 3;
}
double len(){return s.distance(e);}
Point crosspoint(Line v){//求两直线的交点//要保证两直线不平行或重合
double a1=(v.e-v.s)^(s-v.s);
double a2=(v.e-v.s)^(e-v.s);
return Point((s.x*a2-e.x*a1)/(a2-a1),(s.y*a2-e.y*a1)/(a2-a1));
}
void input(){s.input(),e.input();}
void show(){s.show(),e.show();}
};
struct circle{
Point p;//圆心
double r;//半径
circle(){}
circle(Point p,double r):p(p),r(r){}
circle(Point a,Point b,Point c){//三角形的外接圆//利用两条边的中垂线得到圆心
Line u=Line((a+b)/2,((a+b)/2)+((b-a).rotleft()));
Line v=Line((b+c)/2,((b+c)/2)+((c-b).rotleft()));
p=u.crosspoint(v),r=p.distance(a);
}
Point point(double a){return Point(p.x+cos(a)*r,p.y+sin(a)*r);}//根据圆心角算圆上的点
void input(){p.input();scanf("%lf",&r);}
void show(){printf("%lf %lf %lf\n",p.x,p.y,r);}
double area(){return pi*r*r;}
double circumference(){return 2*pi*r;}//周长
int relation(Point b){//点和圆的关系//0 圆外//1 圆上//2 圆内
double dst=b.distance(p);
if(sgn(dst-r)<0)return 2;
else if(!sgn(dst-r))return 1;
return 0;
}
};
Point P[maxp];
int n;
void solve(){
if(n==1){
printf("%lf %lf 1\n",P[0].x+1,P[0].y);
return;
}
if(n<5){
printf("%lf %lf %lf\n",(P[0].x+P[1].x)/2.0,(P[0].y+P[1].y)/2.0,P[0].distance(P[1])/2.0);
return;
}
while(1){
int i=rand()%n,j=rand()%n,k=rand()%n;
if(i==j||i==k||j==k)continue;
if(!sgn((P[j]-P[i])^(P[k]-P[i])))continue;
circle O(P[i],P[j],P[k]);
if(fabs(O.p.x)>1e9||fabs(O.p.y)>1e9||fabs(O.r)>1e9)continue;
int cnt=0;
for(int i=0;i<n;++i)cnt+=O.relation(P[i])==1;
if(cnt<<1>=n){
O.show();return;
}
}
}
int main(){
int t;
scanf("%d",&t);
while(t--){
srand(time(NULL));
scanf("%d",&n);
for(int i=0;i<n;++i)P[i].input();
solve();
}
return 0;
}