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
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
#define rep(i,a,b) for(int i=a;i<b;i++)
const double eps=1e-6;
const double pi=acos(-1.0);
//可以在Point里面用Vec了
class Point;
typedef Point Vec;
//三态函数比较;精度问题
int dcmp(double x){
if(fabs(x)<eps) return 0;
return x<0?-1:1;
}
struct Point{
double x,y;
Point(double _x=0,double _y=0):x(_x),y(_y){}
/*
*向量运算
*/
//向量与常数
Vec operator*(double p){
return Vec(x*p,y*p);
}
Vec operator/(double p){
return Vec(x/p,y/p);
}
//向量与向量
Vec operator-(Vec obj){
return Vec(x-obj.x,y-obj.y);
}
Vec operator+(Vec obj){
return Vec(x+obj.x,y+obj.y);
}
//点积
double operator*(Vec& obj){
return x*obj.x+y*obj.y;
}
//叉积
double operator^(Vec obj){
return x*obj.y-y*obj.x;
}
//两个向量的夹角 A*B=|A|*|B|*cos(th)
double Angle(Vec B){
return acos((*this)*B/(*this).len()/B.len());
}
//两条向量平行四边形的面积
double Area(Vec B){
return fabs((*this)^B);//
}
//向量旋转
//旋转公式
// Nx (cos -sin) x
// =
// Ny (sin cos) y
Vec Rotate(double rad){
return Vec(x*cos(rad)-y*sin(rad),x*sin(rad)+y*cos(rad));
}
//返回向量的法向量,即旋转pi/2
Vec Normal(){
//返回单位法向量,注意L不能为0
/*
double L=obj.len();
return Vec(-y/L,x/L);
*/
//返回法向量
return Vec(-y,x);
}
/*
* 向量的性质
*/
//返回向量的长度,或者点距离原点的距离
double len(){
return hypot(x,y);
}
//返回两点之间的距离
double dis(Point obj){
//return hypot(x-obj.x,y-obj.y); //hypot 给定直角三角形的两条直角边,返回斜边边长
return sqrt((x-obj.x)*(x-obj.x)+(y-obj.y)*(y-obj.y));
}
//向量的极角 atan2(y,x)
/*
*向量的关系
*/
bool operator==(Point obj){
return dcmp(x-obj.x)==0&&dcmp(y-obj.y)==0;
}
bool operator<(Point obj){
return x<obj.x||(x==obj.x&&y<obj.y);
}
};
const int maxn=1e5+10;
Point p[maxn];
//求a,b,c 垂直平分线的交点,好在也过了
/*
//法向量,即旋转90度
Vec Normal(Vec A){
return Vec(-A.y,A.x);
}
Point GetSec(Point P,Vec v,Point Q,Vec w){
Vec u=P-Q;
double t=(w^u)/(v^w);
return P+v*t;
}
Point get_point(Point a,Point b,Point c){
Point A,B;
A.x=(a.x+b.x)/2.0,A.y=(a.y+b.y)/2.0;
B.x=(b.x+c.x)/2.0,B.y=(b.y+c.y)/2.0;
Vec v=Normal(b-a),w=Normal(b-c);
return GetSec(A,v,B,w);
}
*/
//求三个点的外心圆
Point get_point(Point a,Point b,Point c){
double x=( (a.x*a.x-b.x*b.x+a.y*a.y-b.y*b.y)*(a.y-c.y)-(a.x*a.x-c.x*c.x+a.y*a.y-c.y*c.y)*(a.y-b.y) ) / (2*(a.y-c.y)*(a.x-b.x)-2*(a.y-b.y)*(a.x-c.x));
double y=( (a.x*a.x-b.x*b.x+a.y*a.y-b.y*b.y)*(a.x-c.x)-(a.x*a.x-c.x*c.x+a.y*a.y-c.y*c.y)*(a.x-b.x) ) / (2*(a.y-b.y)*(a.x-c.x)-2*(a.y-c.y)*(a.x-b.x));
return Point(x,y);
}
/*
随机的坑:
概率就是概率,尽管1/8的概率,但是100次就是找不出来结果,简直爆炸
下次最少直接100倍,如果必定存在的话,直接while(1)
*/
int main(){
int T;
scanf("%d",&T);
while(T--){
int n;
scanf("%d",&n);
rep(i,0,n){
double x,y;
scanf("%lf %lf",&x,&y);
p[i]=Point(x,y);
}
int wrong=1;
Point po;double r;
//有可能这四个点都在一条直线上
if(n<=4){
if(n==1||n==2){
po.x=p[0].x+1.0,po.y=p[0].y;
r=1.0;
}
if(n==3||n==4){
po=(p[0]+p[1])/2;
r=p[0].dis(p[1])/2.0;
}
}
else{
rep(i,0,1000){
int a=rand()%n,b=rand()%n,c=rand()%n;
if(fabs((p[a]-p[b])^(p[c]-p[b]))<=eps)continue;
po=get_point(p[a],p[b],p[c]);
r=po.dis(p[a]);
int num=0;
rep(j,0,n){
if(fabs(po.dis(p[j])-r)<eps)num++;
if(num*2>=n){
wrong=0;break;
}
}
if(!wrong)break;
}
}
printf("%f %f %f\n",po.x,po.y,r);
}
return 0;
}