Run Away HDU - 1109

题意: 一个左下角在(0,0)的矩阵,长为X,宽为Y,在矩阵中有N个点.现在要求输出一个点,满足到这N个点距离最大值最小.

#include<cstdio>
#include<vector>
#include<cmath>
#include<time.h>
#include<string>
#include<string.h>
#include<iostream>
#include<algorithm>
#define PI acos(-1.0)
#define pb push_back
#define F first
#define S second
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
const int N=1e3+3;
const int MOD=1e9+7;
const double eps=1e-10;
template <class T>
bool sf(T &ret){ //Faster Input
    char c; int sgn; T bit=0.1;
    if(c=getchar(),c==EOF) return 0;
    while(c!='-'&&c!='.'&&(c<'0'||c>'9')) c=getchar();
    sgn=(c=='-')?-1:1;
    ret=(c=='-')?0:(c-'0');
    while(c=getchar(),c>='0'&&c<='9') ret=ret*10+(c-'0');
    if(c==' '||c=='\n'){ ret*=sgn; return 1; }
    while(c=getchar(),c>='0'&&c<='9') ret+=(c-'0')*bit,bit/=10;
    ret*=sgn;
    return 1;
}
struct Point { //点的定义
    double x, y;
    Point(double x=0.0, double y=0.0) : x(x), y(y) {}
    Point operator + (const Point &rhs) const {  //向量加法
        return Point(x + rhs.x, y + rhs.y);
    }
    Point operator - (const Point &rhs) const {  //向量减法
        return Point(x - rhs.x, y - rhs.y);
    }
}p[N];
double X,Y;
int n;
int sign(double x){
    return fabs(x)<x?0:x<0?-1:1;
}
double length(Point a,Point b){
    return sqrt( (a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y) );
}
double cul(Point t){
    double ans=1e20;
    for(int i=0;i<n;i++){
        ans=min(ans,length(t,p[i]));
    }
    return ans;
}
double myrand(){
    return rand()%10000/10000.0; ///rand_MAX=2^15-1;
}
void SA(Point &ans){
    ans={X/2,Y/2};///把矩阵中心当作是最初答案
    double T=sqrt(X*X+Y*Y)/2;///步长是半对角线
    double E=cul(ans);///最初的最大值
    while(T>eps){/// 循环条件大于eps
        Point next; 
        double mmax=0.0; ///保存100个点中最小距离最大的点
        for(int i=0;i<100;++i){/// 随机100个点去扩展
            double x,y;
            double ang=myrand()*2*PI;///随机一个角度
            x=ans.x+cos(ang)*T;
            y=ans.y+sin(ang)*T;
            if(sign(x)<0)   x=0.0;///不超过矩阵范围
            if(sign(x-X)>0)   x=X;
            if(sign(y)<0)   y=0.0;
            if(sign(y-Y)>0)   y=Y;
            Point t(x,y);
            double nE=cul(t); 
            if(sign(nE-mmax)>=0){
                next=t;
                mmax=nE;
            }
        }
        if(sign(mmax-E)>0 || exp((mmax-E)/T>myrand())){ /// 如果比当前值要大 或者 比当前值小,但是满足放弃当前局部最优解的条件.exp中的函数一定要满足和T成反比例
            E=mmax;
            ans=next;
        }
        T*=0.88;/// 乘的数越大,循环次数越多,理论上越逼近全局最优解,但时间会剧增,
    }
}
void mian(){
    scanf("%lf%lf%d",&X,&Y,&n);
    for(int i=0;i<n;i++){
        double x,y;
        scanf("%lf%lf",&x,&y);
        p[i]={x,y};
    }
    Point ans;
    SA(ans);///SA是模拟退火的英文简称
    printf("The safest point is (%.1f, %.1f).\n",ans.x,ans.y);
}
int main(void){
    srand(time(0));
    int T;
    sf(T);
    while(T--){
        mian();
    }

    return 0;
}