HDU - 6242：Geometry Problem（随机+几何）

Posted 2021-01-08 hua-dong

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 NN distinct points (Xi,Yi)(Xi,Yi) on the two-dimensional plane. Your task is to find a point PP and a real number RR, such that for at least ⌈N2⌉⌈N2⌉ given points, their distance to point PP is equal to RR. 

InputThe first line is the number of test cases. 

For each test case, the first line contains one positive number N(1≤N≤105)N(1≤N≤105). 

The following NN lines describe the points. Each line contains two real numbers XiXiand YiYi (0≤|Xi|,|Yi|≤103)(0≤|Xi|,|Yi|≤103) indicating one give point. It‘s guaranteed that NN points are distinct. 
OutputFor each test case, output a single line with three real numbers XP,YP,RXP,YP,R, where (XP,YP)(XP,YP) is the coordinate of required point PP. Three real numbers you output should satisfy 0≤|XP|,|YP|,R≤1090≤|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 RR if it is within an absolute error of 10−310−3 of RR. 
Sample Input

1
7
1 1
1 0
1 -1
0 1
-1 1
0 -1
-1 0

Sample Output

0 0 1

题意：给定N个点，求一个圆，使得圆上的点大于大于一半，保证有解。

思路：既然保证有解，我们就随机得到三角形，然后求外接圆取验证即可。

#include<bits/stdc++.h>
#define rep(i,a,b) for(int i=a;i<=b;i++)
using namespace std;
typedef long long ll;
const double eps=1e-3;
const double pi=acos(-1.0);
struct point{
    double x,y;
    point(double a=0,double b=0):x(a),y(b){}
};
int dcmp(double x){ return fabs(x)<eps?0:(x<0?-1:1);}
point operator +(point A,point B) { return point(A.x+B.x,A.y+B.y);}
point operator -(point A,point B) { return point(A.x-B.x,A.y-B.y);}
point operator *(point A,double p){ return point(A.x*p,A.y*p);}
point operator /(point A,double p){ return point(A.x/p,A.y/p);}
point rotate(point A,double rad){
    return point(A.x*cos(rad)-A.y*sin(rad), A.x*sin(rad)+A.y*cos(rad));
}
bool operator ==(const point& a,const point& b) {
     return dcmp(a.x-b.x)==0&&dcmp(a.y-b.y)==0;
}
double dot(point A,point B){ return A.x*B.x+A.y*B.y;}
double det(point A,point B){ return A.x*B.y-A.y*B.x;}
double dot(point O,point A,point B){ return dot(A-O,B-O);}
double det(point O,point A,point B){ return det(A-O,B-O);}
double length(point A){ return sqrt(dot(A,A));}
double angle(point A,point B){ return acos(dot(A,B)/length(A)/length(B));}
point jiaopoint(point p,point v,point q,point w)
{   //p+tv q+tw，点加向量表示直线，求直线交点
    point u=p-q;
    double t=det(w,u)/det(v,w);
    return p+v*t;
}
point GetCirPoint(point a,point b,point c)
{
    point p=(a+b)/2;    //ab中点
    point q=(a+c)/2;    //ac中点
    point v=rotate(b-a,pi/2.0),w=rotate(c-a,pi/2.0);   //中垂线的方向向量
    if (dcmp(length(det(v,w)))==0)    //平行
    {
        if(dcmp(length(a-b)+length(b-c)-length(a-c))==0) return (a+c)/2;
        if(dcmp(length(b-a)+length(a-c)-length(b-c))==0) return (b+c)/2;
        if(dcmp(length(a-c)+length(c-b)-length(a-b))==0) return (a+b)/2;
    }
    return jiaopoint(p,v,q,w);
}
const int maxn=100010;
point a[maxn]; int F[maxn];
bool check(point S,double R,int N){
    int num=0;
    rep(i,1,N){
        if(dcmp(length(a[i]-S)-R)==0) num++;
    }
    if(num>=(N+1)/2) return true; return false;
}
int main()
{
    int T,N,M;
    scanf("%d",&T);
    while(T--){
        scanf("%d",&N);
        rep(i,1,N) scanf("%lf%lf",&a[i].x,&a[i].y);
        if(N==1) printf("%.6lf %.6lf %.6lf
",a[1].x,a[1].y,0.0);
        else if(N==2) printf("%.6lf %.6lf %.6lf
",(a[1].x+a[2].x)/2,(a[1].y+a[2].y)/2,length(a[1]-a[2])/2);
        else {
            while(true){
                rep(i,1,N) F[i]=i;
                random_shuffle(F+1,F+N+1);
                point S=GetCirPoint(a[F[1]],a[F[2]],a[F[3]]);
                double R=length(S-a[F[1]]);
                if(check(S,R,N)) {
                    printf("%.6lf %.6lf %.6lf
",S.x,S.y,R);
                    break;
                }
            }
        }
    }
    return 0;
}