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

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 R. 

InputThe first line is the number of test cases. 

For each test case, the first line contains one positive number 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个点，求一个圆，使得圆上的点大于大于一半，保证有解。

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

#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\n",a[1].x,a[1].y,0.0);
        else if(N==2) printf("%.6lf %.6lf %.6lf\n",(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\n",S.x,S.y,R);
                    break;
                }
            }
        }
    }
    return 0;
}