JZOJ 5267. 费马点问题

Description

给n个点，每个点的横纵坐标为[0,10000]内的整数，找出一个点，使这个点到其他所有点的距离之和最小，也就是求费马点。
找到的点可能是小数。

Input

第一行一个整数n表示点数，后面n行每行两个整数，表示点的坐标。

Output

第一行一个整数n表示点数，后面n行每行两个整数，表示点的坐标。

Sample Input

5
8 8
10 8
7 10
8 6
5 7

Sample Output

9.40

Data Constraint

对于30%数据n<=10
​对于100%数据n≤100

Solution

• 经典的费马点问题，直接使用模拟退火算法——一个完完全全的贪心求近似最优解的算法。

• 引用一下度娘：

模拟退火算法（Simulate Anneal，SA）是一种通用概率演算法，用来在一个大的搜寻空间内找寻命题的最优解。模拟退火是由S.Kirkpatrick, C.D.Gelatt和M.P.Vecchi在1983年所发明的。V.Černý在1985年也独立发明此演算法。模拟退火算法是解决TSP问题的有效方法之一。

模拟退火的出发点是基于物理中固体物质的退火过程与一般组合优化问题之间的相似性。模拟退火算法是一种通用的优化算法，其物理退火过程由加温过程、等温过程、冷却过程这三部分组成。

• 进入正题，我们在范围内随机 $M=30$ 个可能点，

• 并设置限制 $lim=10^4$ ，降温率 $delta=0.99$ ，温度下限阈值 $T_{min}=0.01$。

• 每次处理一个可能点 $(x_i,y_i)$ ，那么在矩形 $(x±lim,y±lim)$ 中随机出一个点 $(x_i',y_i')$ 。

• $O(N)$ 计算一次 $(x_i',y_i')$ 的答案，若能更新答案，就将 $(x_i,y_i)$ 移到 $(x_i',y_i')$ 。

• 最后把 $lim=lim*delta$ ，则继续处理每个可能点，直到 $lim<=T_{min}$ 时停止。

Code

#include<cstdio>
#include<cstdlib>
#include<cmath>
#include<ctime>
using namespace std;
const int N=10001,M=31,way[4][2]={{1,1},{-1,1},{-1,-1},{1,-1}};
struct data
{
    double x,y;
}a[N],b[M],t;
int n,mnx,mxx,mny,mxy;
double ans=1e9;
inline int read()
{
    int X=0,w=1; char ch=0;
    while(ch<'0' || ch>'9') {if(ch=='-') w=-1;ch=getchar();}
    while(ch>='0' && ch<='9') X=(X<<3)+(X<<1)+ch-'0',ch=getchar();
    return X*w;
}
inline double get()
{
    double s=0;
    for(int i=1;i<=n;i++) s+=sqrt((t.x-a[i].x)*(t.x-a[i].x)+(t.y-a[i].y)*(t.y-a[i].y));
    return s;
}
inline int random(int x,int y)
{
    return x+rand()%(y-x+1);
}
int main()
{
    srand(time(NULL));
    n=read();
    mnx=mny=N;
    for(int i=1;i<=n;i++)
    {
        a[i].x=read(),a[i].y=read();
        if(a[i].x<mnx) mnx=a[i].x;
        if(a[i].x>mxx) mxx=a[i].x;
        if(a[i].y<mny) mny=a[i].y;
        if(a[i].y>mxy) mxy=a[i].y;
    }
    for(int i=1;i<M;i++)
    {
        b[i].x=random(mnx*N,mxx*N)*1.0/N;
        b[i].y=random(mny*N,mxy*N)*1.0/N;
        t=b[i];
        double p=get();
        if(p<ans) ans=p;
    }
    for(double lim=N;lim>0.01;lim*=0.99)
        for(int i=1;i<M;i++)
        {
            t.x=random((b[i].x-lim)*N,(b[i].x+lim)*N)*1.0/N;
            t.y=random((b[i].y-lim)*N,(b[i].y+lim)*N)*1.0/N;
            double p=get();
            if(p<ans) ans=p,b[i]=t;
        }
    printf("%.2lf",ans);
    return 0;
}