一篇教会你ACM模拟退火SA

最新推荐文章于 2020-04-25 23:49:41 发布

转载最新推荐文章于 2020-04-25 23:49:41 发布 · 置顶 · 2.5k 阅读

22 ·

CC 4.0 BY-SA版权

原文链接：https://cardiffle.github.io/2019/08/05/%E6%A8%A1%E6%8B%9F%E9%80%80%E7%81%ABSA/

本文深入解析模拟退火算法(SA)的原理及其在优化问题中的应用，通过对比爬山算法，阐述SA如何通过随机性跳出局部最优解，达到全局最优。文章详细介绍了SA的流程、参数设置及在各类问题中的应用案例。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

（原创文章，转载自自己博客，转载请标明来源。）

模拟退火SA全解

简介

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

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

以上内容来自度娘，所以说，SA 是一种非常玄学的随机化算法，非常适合在比赛中遇到非常不好用正常方法解决时强行A题，非常适合我在队里的位置（~~抱紧大腿QaQ~~）可能交一页就能过了（汪汪大笑）。

爬山算法

爬山算法的思路是，每次在当前找到的方案附件寻找一个新的方案（常用方式是随机一个差值），如果这个解更优就直接转移。

对于单峰函数来说，这样可以找到最优解了（但是对于单峰函数，为什么我们不三分呢）

对于多数我们求解的函数，导致麻烦的情况多时长成上图这个样子，如果简单地爬山可能爬到一个最优解里爬不出来了。

显然算法是否得出最优解与初始解的位置以及搜寻的附近解的区域大小有关。

但是如果我们的 $\Delta$ 设置的比较大（左右横跳一次的半径）的话有几率跳出去，但是太大了又可能跳来跳去从而找不到最优解。~~欧皇的一键AK算法~~

对于我这样的非酋无非太不友好了。

模拟退火

简单说，模拟退火是一种随机化算法，用于求函数的极值。当一个问题的方案数量极大（甚至是无穷的）而且不是一个单峰函数时，我们常使用模拟退火求解。

它与爬山算法最大的不同是，在寻找到一个局部最优解时，赋予了它一个跳出去的概率，也就有更大的机会能找到全局最优解。

原理

模拟退火的原理也和金属退火的原理近似：将热力学的理论套用到统计学上，将搜寻空间内每一点想像成空气内的分子；分子的能量，就是它本身的动能；而搜寻空间内的每一点，也像空气分子一样带有“能量”，以表示该点对命题的合适程度。演算法先以搜寻空间内一个任意点作起始：每一步先选择一个“邻居”，然后再计算从现有位置到达“邻居”的概率。

——百度百科

在 OI 领域，对应的，每次随机出一个新解，如果这个解更优，则接受它，否则以一个与温度和与最优解的差相关的概率接受它。

通过对热力学定律的模拟，结合计算机对离散数据的处理，我们定义：如果当前状态与最优状态的能量差为 $\Delta E$ ，当前温度为 $T$ ， $k$ 为一个随机数，则发生状态转移的概率是：

$P(\Delta E) = e^{\frac{\Delta E} {kT}}$

显然如果 $ \Delta E >0$ ，转移一定会成功，但是如果对于 $ \Delta E <0$ ，我们则会计算概率接受这个新解。

过程

模拟退火的过程中维护T，另包含三个参数：初始温度 $T_0$ ，降温系数 $\Delta$ ，中止温度 $T_k$ 。

$T_0$ 是一个比较大的数，$ \Delta$ 是一个略小于1的正数， $T_k$ 是一个略大于 0 的正数。

先让温度 $T=T_0$ ，然后每次降温时 $T=T·\Delta$ ，直到 $T≤T_k$ 为止。此时当前接为最优解。

Wiki 上的有名动图：

Attention

程序开始时，我们先要 srand(一个常数)。这个常数可以决定分数，可以使用 233333，2147483647 等等。

一遍模拟退火往往跑出最优解，这时候不要慌，常数一换，又是一条好汉，~~冲冲冲~~

可以用一个全局变量记录所有跑过的 SA 的最优解，每次从哪个最优解开始继续 SA，可以减少误差。

时间复杂度

时间复杂度： $O (玄学）$ 或是 $O (脸白程度）$

一般降温系数 $\Delta$ 与 1 的差减少一个数量级, 耗时大约多 10 倍； $T_0$ 和 $T_k$ 变化一个数量级, 耗时不会变化很大。

模拟退火的实际应用

洛谷P1337（JSOI2004）吊打XXX

传送门

Description

gty又虐了一场比赛，被虐的蒟蒻们决定吊打gty。gty见大势不好机智的分出了n个分身，但还是被人多势众的蒟蒻抓住了。蒟蒻们将n个gty吊在n根绳子上，每根绳子穿过天台的一个洞。这n根绳子有一个公共的绳结x。吊好gty后蒟蒻们发现由于每个gty重力不同，绳结x在移动。蒟蒻wangxz脑洞大开的决定计算出x最后停留处的坐标，由于他太弱了决定向你求助。不计摩擦，不计能量损失，由于gty足够矮所以不会掉到地上。

Input

输入第一行为一个正整数n(1<=n<=10000)，表示gty的数目。
接下来n行,每行三个整数xi，yi，wi，表示第i个gty的横坐标，纵坐标和重力。
对于20%的数据，gty排列成一条直线。
对于50%的数据，1<=n<=1000。
对于100%的数据，1<=n<=10000,-100000<=xi,yi<=100000

Output

输出1行两个浮点数（保留到小数点后3位），表示最终x的横、纵坐标。

Sample Input

3
0 0 1
0 2 1
1 1 1

Sample Output

0.577 1.000

Hint

Time Limit: 10 Sec

Memory Limit: 128 MB Sec Special Judge

首先，根据自然规律，一切自然变化进行的方向都是使能量降低，因为能量较低的状态比较稳定。

因为物重一定，绳子越短，重物越低，势能越小，势能又与物重成正比，所以，只要使得 $\sum_{i=1}^ndist[i]*weight[i]$ 也就是总的重力势能最小，就可以使系统平衡。

这个时候我们可以考虑模拟退火. 首先随机一个点 $T_0$ 作为初始解(为了加速收敛, 我们可以直接取各个点坐标的平均值所在的点). 然后随机两个值作为差值加到这个点的坐标上作为下一个解.

然后模拟退火直接往上套就可以了~~233~~。

具体实现就是一个 while 循环, 循环内有4步:

根据当前解找到下一个解
计算下一个解的 “能量” (也就是价值)
决定是否要接受这个新解
降温

找下一个解的时候有一个提高精度的小技巧: 根据当前温度决定差值的范围. 这样在降温即将结束接近最优解的时候可以有更大的概率更精确地命中最优解。

当然如果是解是离散的就不能这样搞了. 以及生成下一个解的时候万万不能全部重新随机生成, 那就和瞎随没区别了…要随机作出一些相对小的修改.

具体做法就是使用一个产生 $[0, 1]$ 随机实数的函数, 将随机区间转为 $[- 1, 1]$ 后乘上 $T$ 作为差值. (也就是生成一个 $[- T, T]$ 的随机值作为差值)

不过实际操作的时候我们较少直接输出最终解, 而是选择在模拟退火的过程中单独维护一个解, 只在遇到更优解的时候将其更新, 增加正确率.

另一个提高正确率的方法是: 多次进行模拟退火过程(或者说"重新烧热再退火一遍"), 每次取最优解.

还有就是最后烧完之后可以再在全局最优解的基础上进行爬山.

总结一下就是神仙方法集萃

这里要注意rand()和rand()*2-RAND_MAX的区别
rand()的范围是0~RAND_MAX-1
rand()**2-RAND_MAX 的范围是-RAND_MAX到RAND_MAX-1

AC代码：

#include <bits/stdc++.h>
#define re register
using namespace std;

inline int read() { //读入优化
    int X=0,w=1; char c=getchar();
    while (c<'0'||c>'9') { if (c=='-') w=-1; c=getchar(); }
    while (c>='0'&&c<='9') X=(X<<3)+(X<<1)+c-'0',c=getchar();
    return X*w;
}

struct node { int x,y,w; };

node a[1010];
int n,sx,sy;

double ansx,ansy; //全局最优解的坐标
double ans=1e18,t; //全局最优解、温度
const double delta=0.993;//0.993; //降温系数

inline double calc_energy(double x,double y) { //计算整个系统的能量
    double rt=0;
    for (re int i=1;i<=n;i++) {
        double deltax=x-a[i].x,deltay=y-a[i].y;
        rt+=sqrt(deltax*deltax+deltay*deltay)*a[i].w;
    }
    return rt;
}

inline void simulate_anneal() { //SA主过程
    double x=ansx,y=ansy;
    t=2000; //初始温度
    while (t>1e-14) {
        double X=x+((rand()<<1)-RAND_MAX)*t;
        double Y=y+((rand()<<1)-RAND_MAX)*t; //得出一个新的坐标
        double now=calc_energy(X,Y);
        double Delta=now-ans;
        if (Delta<0) { //接受
            x=X,y=Y;
            ansx=x,ansy=y,ans=now;
        }
        else if (exp(-Delta/t)*RAND_MAX>rand()) x=X,y=Y; //以一个概率接受
        t*=delta;
    }
}

inline void Solve() { //多跑几遍SA，减小误差
    ansx=(double)sx/n,ansy=(double)sy/n; //从平均值开始更容易接近最优解
    simulate_anneal();
    simulate_anneal();
    simulate_anneal();
    simulate_anneal();
    simulate_anneal();
}

int main() {
    srand(18253517); srand(rand());srand(rand()); //玄学srand
    //cout << RAND_MAX;//32767
    n=read();
    for (re int i=1;i<=n;i++) {
        a[i].x=read(),a[i].y=read(),a[i].w=read();
        sx+=a[i].x,sy+=a[i].y;
    }
    Solve();
    printf("%.3f %.3f\n",ansx,ansy);
    return 0;
}

$[外链图片转存失败(img-l9p7R9qw-1564972608934)(C:\Users\12\AppData\Roaming\Typora\typora-user-images\1564814223351.png)]$

$[外链图片转存失败(img-LZEnhqyg-1564972608934)(C:\Users\12\AppData\Roaming\Typora\typora-user-images\1564814907097.png)]$
交了五次就过了，今天有点小欧~~(汪汪大笑)~~

最终版本是Delta取了0.993，跑了5次SA， $T_0$ 取20001，srand(rand())*2。

总结一下，如果答案不是最优的怎么办？

首先看看自己的判定过程写的是否正确，确定正确了之后。

1.调大Delta

2.调大 $T_0$

3.调小 $T_k$

4.多跑几遍SA （居然看到带佬建议跑42次，155次）

5.srand(rand())* 2 ~~不要问我为什么~~

6.更换随机数种子SEED

下面是一张表供大家估算运行时间，左边是“降温系数”，上方是初温与末温的比值，表格内容是大致的迭代次数。

在这里插入图片描述
从上表可以看出， $\frac {T_0}{T_k}$ 涨一个数量级只会增加约25%的迭代次数，而1-Delta 涨一个数量级会增加十倍的运行时间。

不同类型题目如何生成新解

坐标系内：随机生成一个点，或者生成一个向量。
序列问题： $r a n d o m$ _ $s h u f f l e$ 或者随机交换两个数。
网格问题：可以看做二维序列，每次交换两个格子即可。

$r a n d o m$ _ $s h u f f l e$ 用法

函数原型定义：
template<class RandomAccessIterator>
   void random_shuffle(
      RandomAccessIterator _First, //指向序列首元素的迭代器
      RandomAccessIterator _Last  //指向序列最后一个元素的下一个位置的迭代器
);
字符串和数组中应用实例
#include <bits/stdc++.h>
using namespace std;
 
int main()
{
    char b[]="abcde";
    char a[10];
    cout<<"char:"<<endl;
    for(int i=0;i<10;i++)
    {
        strcpy(a,b);
        random_shuffle(a,a+5);
        cout<<a<<endl;
    }
 
    int c[]={1,2,3,4,5};
    int d[10];
    cout<<endl<<"int:"<<endl;
    for(int i=0;i<10;i++)
    {
        memcpy(d,c,sizeof(c));
        random_shuffle(d,d+5);
        for(int j=0;j<5;j++) cout<<d[j]<<" ";
        cout<<endl;
    }
}

简单应用

1.POJ-2420 A Star not a Tree?

求费马点：在一个确定的矩形里找一个点，最大化与其他已有点的距离的最小值。

技巧是开局确定一个为对角线一半的值为最大步长 $T_0$ ，每一次在已有ans点的位置随机角度 $\pi$ ，在用当前步长乘对应角度三角函数得到新位置。

#include <cstdio>
#include <cmath>
#include <cstring>
#include <string>
#include <map>
#include <queue>
#include <iostream>
#include <algorithm>
#include <time.h>
#include <iomanip>
using namespace std;
#define RE freopen("1.in","r",stdin);
#define SpeedUp std::cout.sync_with_stdio(false);

const int maxn=1e5+50;
const int inf = 0x7FFFFFFF;    //可改勿删，不知名bug，以前遇到过是越界导致，但我并没有找到越界？
const double PI = acos(-1.0);

int n;
double X,Y;
struct Point
{
    double x,y;
    Point(){}
    Point(double _,double __):x(_),y(__){}
    void read(){
        cin>>x>>y;
    }
}p[maxn];

double dis(Point a,Point b){
    return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}
double minDis(Point cur){
    double ans = 1e80;
    for (int i = 0; i < n; ++i){
        ans = min(ans,dis(cur,p[i]));
    }
    return ans;
}
double randNum(){       //rand()生成[0,32767),包装下生成[0,1)
    return rand()%10000/10000.0;
}

void solve(Point &ans,double T){
    double T_min = 1e-8;
    double E = minDis(ans);//当前高度，即估值
    int count = 30;//对当前点取周围count个点，再从里面挑最高的当下一个点
    while(T > T_min){ //温度，也即本题步长
        Point next;
        double nE = 0.0;
        for (int i = 0; i < count; ++i){
            Point tmp;
            double angle = randNum()*2*PI;     //[0,2π)
            tmp.x = ans.x + cos(angle)*T;      //ans.x + [-1,1)*T
            tmp.y = ans.y + sin(angle)*T;
            tmp.x = min(X,max(0.0,tmp.x));
            tmp.y = min(Y,max(0.0,tmp.y));}

int main(){
    // RE
    SpeedUp
    int t;
    srand(time(NULL));  //c++提交，g++会因为这个RE
    cin>>t;
    while(t--){
        cin>>X>>Y>>n;
        for (int i = 0; i < n; ++i){
            p[i].read();
        }
        Point ans = Point(X/2,Y/2); 
               //起始为中心点
        double T = sqrt(X*X+Y*Y)/2.0;     //步长为对角线长度，包含了整个域}
        cout <<"**"<< T;
        solve(ans,T);
        cout <<setiosflags(ios::fixed);
        cout<<"The safest point is ("<<setprecision(1)<<ans.x<<", "<<setprecision(1)<<ans.y<<")."<<endl;
    }
    return 0;
}

2.HDU5017 Ellipsoid

题意：求一个椭球面上的一个点到原点的最短距离。

$ax^2+by^2+cz^2+dyz+exz+fxy=1$

思路：这个题网上的题解写的乱七八糟的，着实误导人，看代码才明白。虽然是三维的，但是空间枚举坐标无法保证可成立行，所以在二维平面枚举 x , y ，计算出 z 的位置（把椭圆方程化成二次函数求根公式求z），八个方向扩展找领域解。正解是三分套三分（单调性实锤）

/*没有用到随机的退火，我觉得是个爬山
   某lyh大神的推论椭球二次函数都是单峰的，那根本就没必要退火
   坑点大概就是一开始看了一份假题解
   后来WA就是没考虑方程的负解
   还有就是参数类型错误，一定要用double
*/
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <algorithm>
#include <iostream>
using namespace std;
#define LL long long
#define CLR(x) memset(x, 0, sizeof(x))
const double INF = 1e10;
const double eps = 1e-9;

int dx[8] = {-1, -1, -1, 0, 0, 1, 1, 1}; //因为是三维坐标所有有八个方向走
int dy[8] = {-1, 0, 1, -1, 1, -1, 0, 1};

double a, b, c, d, e, f;

double dist(double x, double y, double z)
{
    return sqrt(x*x + y*y + z*z);
}

double GetVal(double x, double y) 
{ //啦啦解方程
    double A = c;
    double B = e*x + d*y;
    double C = a*x*x + b*y*y + f*x*y-1;
    double D  = B*B - 4*A*C; 
    if (D < 0) return INF;

    double det = sqrt(D);

    double res1 = (-B + det) / (2*A);
    double res2 = (-B - det) / (2*A);//求根公式

    return dist(x,y,res1) < dist(x,y,res2) ? res1 : res2;
}

double Search()
{
    double T = 1; //初始温度
    double delta = 0.99;

    double x = 0, y = 0, z = sqrt(1/c);

    while (T > eps)
    {
        for (int i = 0; i < 8; i++)
        {
            double r = x + dx[i]*T;
            double c = y + dy[i]*T;
            double k = GetVal(r, c);
            if (k > INF) continue;
            if (dist(r,c,k) < dist(x,y,z))
            {
                x = r; 
                y = c;
                z = k;
            }
        }
        T *= delta;
    }
    return dist(x, y, z);

}

int main()
{
    while(scanf("%lf%lf%lf%lf%lf%lf", &a, &b, &c, &d, &e, &f) != EOF)
    {
        printf("%.7lf\n", Search());
    }    
}

在这里插入图片描述

3.HDU 1109 Run Away 类似求费马点

题意：给n个点，在平面内选一个点A，使A到所有点的距离和最小，并求最小距离。

#include <stdio.h>
#include <string.h>
#include <math.h>
#include <algorithm>
#include <time.h>
#include <iostream>
using namespace std;
#define LL long long
#define CLR(x) memset(x, 0, sizeof(x))
const double INF = 1e8;
const double eps = 1e-3;

int dx[5] = {-1, 0, 1, 0}; //平面内
int dy[5] = {0, 1, 0, -1};

double X, Y;
int M;
double U[1005], V[1005];

struct node //数组和结构体都可用
{
    double x, y;
    double val;
} p[55];

double dist(double x, double y, double a, double b)
{
    return (x-a)*(x-a) + (y-b)*(y-b);
}

double GetSum(node h)
{
    double d = INF;
    for (int i = 0; i < M; i++)
        d = min(d, dist(h.x, h.y, U[i], V[i])); 
    return d;
}


double Search()
{
    double T = X+Y;
    double delta = 0.9;
    for (int i = 0; i <= 50; i++)
    {
        p[i].x = (rand()%1000+1) * 1.0/1000.00 * X; //队内dalao说可以暴力，不用随机
        p[i].y = (rand()%1000+1) * 1.0/1000.00 * Y;
        p[i].val = GetSum(p[i]);
    }
    node tmp;
    while (T > eps)
    {
        for (int i = 0; i <= 50; i++)
        {
            for (int j = 0; j <= 50; j++)
            {
                double r = (double)(rand()%1000+1)*1.0/1000.00*T;
                double c = sqrt(T*T - r*r);

                if (rand()&1) r *= -1;
                if (rand()&1) c *= -1;

                tmp.x = p[i].x + r;
                tmp.y = p[i].y + c;

                if (tmp.x<0 || tmp.x>X || tmp.y<0 || tmp.y>Y) continue;

                if (GetSum(tmp) > p[i].val)
                {
                    p[i] = tmp;
                    p[i].val = GetSum(tmp);
                }
            }
        }
        T *= delta;
    }
    int cur = 1;
    for (int i = 0; i <= 50; i++)
        if (p[i].val > p[cur].val)
            cur = i;
    printf("The safest point is (%.1lf, %.1lf).\n", p[cur].x, p[cur].y);
   //不要漏了句号！！！！血一般的教训
}

int main()
{
    srand(unsigned(time(0)));
    int t;
    scanf("%d", &t);
    while(t--)
    {
        CLR(U);
        CLR(V);
        scanf("%lf%lf%d", &X, &Y, &M);
        for (int i = 0; i < M; i++)
            scanf("%lf%lf", &U[i], &V[i]);
        Search();
    }
}

4.HDU 3932Groundhog Build Home 平面上最小园覆盖问题

在这题采用了分块随机数。

#include <stdio.h>
#include <string.h>
#include <math.h>
#include <time.h>
#include <algorithm>
#include <iostream>
using namespace std;
#define LL long long
#define CLR(x) memset(x, 0, sizeof(x))
#define MAX 100000000*1LL
#define PI acos(-1)
const double eps = 1e-2;
const double INF = 1e8;

double X, Y;
int N;
double px[1005], py[1005];

struct node
{
    double x, y, val;
} p[1005];

double dist(double x, double y, double a, double b)
{
    return sqrt((x-a)*(x-a) + (y-b)*(y-b));
}

double GetMax(node t)
{
    double d = 0;
    for (int i = 0; i < N; i++)
        d = max(d, dist(t.x, t.y, px[i], py[i])); //注意，这里是最大值
    return d;
}

void Search()
{
    double delta = 0.6; //！！！太接近1会TLE
    double T = max(X, Y);
    node tmp;

    while(T > eps)
    {
        for(int i = 0; i < N; i++)
        {
            for(int j = 0; j < 20; j++) //这里也是，过大会TLE
            {
                double d = (rand()%1000+1)/1000.0*2*PI;
                tmp.x = p[i].x + cos(d)*T;
                tmp.y = p[i].y + sin(d)*T;

                if (tmp.x<=0 || tmp.x>=X || tmp.y<=0 || tmp.y>=Y) continue;

                if (GetMax(tmp) < p[i].val)
                {
                    p[i] = tmp;
                    p[i].val = GetMax(tmp);
                }
            }
        }
        T *= delta;
    }

    node res;
    res.val = INF;;
    for (int i = 0; i < N; i++) // 求出全局最优解
        if (p[i].val < res.val) res = p[i];

    printf("(%.1lf,%.1lf).\n", res.x, res.y);
    printf("%.1lf\n", res.val);
}

int main()
{
    srand(unsigned(time(0)));
    while(scanf("%lf%lf%d", &X, &Y, &N) != EOF)
    {
        CLR(px);
        CLR(py);
        for (int i = 0; i < N; i++)
            scanf("%lf%lf", &px[i], &py[i]);

        for (int i = 0; i < N; i++) /*随机取N个点。。。笔者原本是最直接改上一题的代码。。。
                                      结果出事情了。。。事实告诉我们随机数要谨慎*/
        {
            p[i].x = (rand()%1000+1)/1000.0*X;
            p[i].y = (rand()%1000+1)/1000.0*Y;
            p[i].val = GetMax(p[i]);
        }
        Search();
    }
}

5.HDU - 3644：A Chocolate Manufacturer’s Problem

模拟退火，求多边形内最大圆半径

因为计算几何还不太会，所以找了篇题解参考一

**pro：**给定一个N边形，然后给半径为R的圆，问是否可以放进去。问题转化为多边形的最大内接圆半径。（N<50）；

**sol：**乍一看，不就是二分+半平面交验证是否有核的板子题吗。然而事情并没有那么简单。因为我们的多边形可能是凹多边形，而前面的方法只对凸多边形有效。

学习了下模拟退火的算法，这个随机算法只在最小圆覆盖的时候写过。这里再学一下，看起来更正宗一点的。每次在当前点的附近（R）找是否能优化，而这个R慢慢变小，使得趋紧答案的趋势更精细。

**判定点再多边形内：**同样，不能用检验是否在每条边的左边来判定，因为不是凸多边形；我们可以用射线法搞。

#include<bits/stdc++.h>
#define rep(i,a,b) for(int i=a;i<=b;i++)
using namespace std;
const int maxn=110;
const double pi=acos(-1.0);
const double inf=0x7fffffff;
struct point {
    double x,y;
    point(){}
    point(double xx,double yy):x(xx),y(yy){}
};
struct line{
    point s,p;
    line(){}
    line(point xx,point yy):s(xx),p(yy){}
};
double getdis(point w,point v){
    return sqrt((w.x-v.x)*(w.x-v.x)+(w.y-v.y)*(w.y-v.y));
}
point operator /(point a,double t){ return point(a.x/t,a.y/t);}
point operator *(point a,double t){ return point(t*a.x,t*a.y);}
point operator -(point w,point v){return point(w.x-v.x,w.y-v.y);}
point operator +(point w,point v){return point(w.x+v.x,w.y+v.y);}
double det(point w,point v){ return w.x*v.y-w.y*v.x;}
double dot(point w,point v){ return w.x*v.x+w.y*v.y;}
double ltoseg(point p,point a,point b){
    point t=p-a;
    if(dot(t,b-a)<=0) return getdis(p,a);
    else if(dot(p-b,a-b)<=0) return getdis(p,b);
    return fabs(det(t,b-a))/getdis(a,b);
}
point p[maxn],tp[maxn]; double dist[maxn]; int N; line L[maxn];
bool isinside(point a)
{
//算法描述：首先，对于多边形的水平边不做考虑，其次，
//对于多边形的顶点和射线相交的情况，如果该顶点时其所属的边上纵坐标较大的顶点，则计数，否则忽略该点，
//最后，对于Q在多边形上的情形，直接判断Q是否属于多边形。
    int ncross=0;
    rep(i,0,N-1) {
        point p1=p[i],p2=p[i+1];
        if(ltoseg(a,p[i],p[i+1])==0) return true; //在线段上
        if(p1.y==p2.y) continue; //默认做水平x轴的线，所以水平线不考虑
        if(a.y<min(p1.y,p2.y)) continue; //相离不考虑
        if(a.y>max(p1.y,p2.y)) continue;
        double t=det(a-p[i],a-p[i+1]);
        if((t>=0&&p[i].y<a.y&&p[i+1].y>=a.y)||(t<=0&&p[i+1].y<a.y&&p[i].y>=a.y)) ncross++;
    }
    return (ncross&1);
}
double getmindis(point a)
{
    double ans=inf;
    rep(i,0,N-1) ans=min(ans,ltoseg(a,p[i],p[i+1]));
    return ans;
}
int main()
{
    srand(unsigned(time(NULL)));
    while(~scanf("%d",&N)&&N){
        double X,Y,R; X=Y=0;
        rep(i,1,N) {
            scanf("%lf%lf",&p[i].x,&p[i].y);
            X=max(X,p[i].x);
            Y=max(Y,p[i].y);
        }
        p[0]=p[N];
        rep(i,0,N-1) L[i]=line(p[i],p[i+1]-p[i]);
        scanf("%lf",&R);
        int maxt=min(N,20);
        rep(i,0,maxt-1){
            tp[i]=(p[i]+p[i+1])/2;
            dist[i]=0;
        }
        double step=min(X,Y);
        const int maxd=10;
        const double rate=0.55;
        bool flag=0;
        const double EPS2=1e-6;
        while(step>EPS2&&!flag){
            rep(i,0,maxt-1){
                rep(j,0,maxd-1){
                    double d=rand()%360/360.0*2*pi;
                    point next=tp[i];
                    next.x+=step*sin(d);
                    next.y+=step*cos(d);
                    if(!isinside(next)) continue;
                    double tdis=getmindis(next);
                    if(tdis+EPS2>dist[i]){
                        dist[i]=tdis; tp[i]=next;
                    }
                    if(tdis+EPS2>=R){
                        flag=1; break;
                    }
                }
            }
            step*=rate;
        }
        if(flag) puts("Yes");
        else puts("No");
    }
    return 0;
}

6.Gym-101158J Cover the Polygon with Your Disk

模拟退火+圆与多边形的交

题意：一个多边形，给定圆的半径，求圆能和多边形相交的最大面积。

正解：把多边形整体区域三分，然后套圆和多边形相交的模板。

#include<cstdio>
#include<cmath>
#include<iostream>
#include<algorithm>
#include<iomanip>
#include<vector>
using namespace std;

const int MAXN = 200000 + 5;
const double PI = 3.1415926;
const double EPS = 1e-10;

double R;//反演圆半径

int judgeZero(double a)//a>0为1,<0为-1，=0为0
{
    return (a > EPS) - (a < -EPS);
}

struct Point
{
    double _x, _y;
    Point(double x = 0.0, double y = 0.0) :_x(x), _y(y) {}
    Point(const Point& p) { _x = p._x, _y = p._y; }
    bool operator<(const Point& p)const
    {
        if (judgeZero(_x - p._x) != 0) return _x < p._x;
        return _y < p._y;
    }
    bool operator==(const Point& p)const
    {
        return judgeZero(_x - p._x) == 0 && judgeZero(_y - p._y) == 0;
    }
    void toMove(Point a, double rad, double d)
    {
        _x = a._x + cos(rad)*d;
        _y = a._y + sin(rad)*d;
    }
    Point operator+(Point a) { return Point(_x + a._x, _y + a._y); }
    Point operator-(Point a) { return Point(_x - a._x, _y - a._y); }
    friend Point operator*(double a, Point p) { return Point(a*p._x, a*p._y); }
    friend istream& operator >> (istream& in, Point& point)
    {
        in >> point._x >> point._y;
        return in;
    }
    friend ostream& operator<<(ostream& out, const Point& point)
    {
        out << fixed << setprecision(8) << point._x << ' ' << point._y;
        return out;
    }
}crossp[MAXN];

double getDis(Point a, Point b)
{
    return sqrt((a._x - b._x)*(a._x - b._x) + (a._y - b._y)*(a._y - b._y));
}

typedef vector<Point> Polygon;

struct Circle
{
    Point _o;
    double _r;
    Circle(double x = 0.0, double y = 0.0, double r = 0.0) :_o(x, y), _r(r) {}
    Circle(const Point& o, double r) :_o(o), _r(r) {}
    Circle getAnti(const Point& point)
    {
        Circle antic;
        double dis = getDis(point, _o);
        double tmp = R*R / (dis*dis - _r*_r);
        antic._r = tmp*_r;
        antic._o._x = point._x + tmp*(_o._x - point._x);
        antic._o._y = point._y + tmp*(_o._y - point._y);
        return antic;
    }
    friend istream& operator >> (istream& in, Circle& circle)
    {
        in >> circle._o >> circle._r;
        return in;
    }
    friend ostream& operator<<(ostream& out, const Circle& circle)
    {
        out << circle._o << ' ' << circle._r;
        return out;
    }
};

double getCross(Point p1, Point p2, Point p)
{
    return (p1._x - p._x)*(p2._y - p._y) - (p2._x - p._x)*(p1._y - p._y);
}

bool onSegment(Point p, Point a, Point b)//点p与线段ab
{
    if (judgeZero(getCross(a, b, p)) != 0)
        return false;
    if (p._x<min(a._x, b._x) || p._x>max(a._x, b._x))//trick：垂直或平行坐标轴;
        return false;
    if (p._y<min(a._y, b._y) || p._y>max(a._y, b._y))
        return false;
    return true;
}

bool segmentIntersect(Point a, Point b, Point c, Point d)//线段ab与线段cd
{
    if (onSegment(c, a, b) || onSegment(d, a, b) || onSegment(a, c, d) || onSegment(b, c, d))
        return true;
    if (judgeZero(getCross(a, b, c)*getCross(a, b, d)) < 0 && judgeZero(getCross(c, d, a)*getCross(c, d, b)) < 0)//==0的特殊情况已经在前面排除
        return true;
    return false;
}

bool onLine(Point p, Point a, Point b)//点p与直线ab
{
    return judgeZero(getCross(a, b, p)) == 0;
}

bool lineIntersect(Point a, Point b, Point c, Point d)//直线ab与直线cd
{
    if (judgeZero((a._x - b._x)*(c._y - d._y) - (c._x - d._x)*(a._y - b._y)) != 0)
        return true;
    if (onLine(a, c, d))//两直线重合
        return true;
    return false;
}

bool segment_intersectLine(Point a, Point b, Point c, Point d)//线段ab与直线cd
{
    if (judgeZero(getCross(c, d, a)*getCross(c, d, b)) <= 0)
        return true;
    return false;
}

Point point_line_intersectLine(Point p1, Point p2, Point p3, Point p4)//求出交点,直线(线段)p1p2与直线(线段)p3p4
{//t=lamta/(lamta+1),必须用t取代lamta,不然算lamta可能分母为0
    double x1 = p1._x, y1 = p1._y;
    double x2 = p2._x, y2 = p2._y;
    double x3 = p3._x, y3 = p3._y;
    double x4 = p4._x, y4 = p4._y;
    double t = ((x2 - x1)*(y3 - y1) - (x3 - x1)*(y2 - y1)) / ((x2 - x1)*(y3 - y4) - (x3 - x4)*(y2 - y1));
    return Point(x3 + t*(x4 - x3), y3 + t*(y4 - y3));
}

bool inPolygon(Point a, Polygon polygon)//点是否含于多边形
{
    Point b(-1e15 + a._x, a._y);//向左无穷远的线段
    int count = 0;
    for (int i = 0;i < polygon.size();++i)
    {
        Point c = polygon[i], d = polygon[(i + 1) % polygon.size()];
        if (onSegment(a, c, d)) //如果点在线段上，那么一定在多边形内
            return true;
        if (judgeZero((a._x - b._x)*(c._y - d._y) - (c._x - d._x)*(a._y - b._y)) == 0)
            continue;//如果边与射线平行,trick1
        if (!segmentIntersect(a, b, c, d))
            continue;//如果边与射线没有交点
        Point lower;
        if (c._y < d._y) lower = c;
        else lower = d;
        if (onSegment(lower, a, b)) //如果纵坐标小的点与射线有交点,trick2
            continue;
        count++;
    }
    return count % 2 == 1;
}

bool segment_inPolygon(Point a, Point b, Polygon polygon)//线段是否含于多边形
{
    if (!inPolygon(a, polygon) || !inPolygon(b, polygon))//两个端点都不在多边形内
        return false;
    int tot = 0;
    for (int i = 0;i < polygon.size();++i)
    {
        Point c = polygon[i], d = polygon[(i + 1) % polygon.size()];
        if (onSegment(a, c, d))//以下只用记录一个交点，多的要么重复，要么一定在边上
            crossp[tot++] = a;
        else if (onSegment(b, c, d))
            crossp[tot++] = b;
        else if (onSegment(c, a, b))
            crossp[tot++] = c;
        else if (onSegment(d, a, b))
            crossp[tot++] = d;
        else if (segmentIntersect(a, b, c, d))//端点没有在线段上且相交
            return false;
    }
    sort(crossp, crossp + tot);//按x，y排序
    tot = unique(crossp, crossp + tot) - crossp;
    for (int i = 0;i < tot - 1;++i)
    {
        Point tmp = 0.5*(crossp[i] + crossp[i + 1]);
        if (!inPolygon(tmp, polygon))//中点不在多边形内
            return false;
    }
    return true;
}

double lineDis_inPolygon(Point a, Point b, Polygon polygon)//求出直线在多边形内的长度
{
    int tot = 0;
    for (int i = 0;i < polygon.size();++i)
    {
        Point c = polygon[i], d = polygon[(i + 1) % polygon.size()];
        if (onLine(c, a, b) && onLine(d, a, b))//如果边与直线重合，记录两个端点
        {
            crossp[tot++] = c;
            crossp[tot++] = d;
        }
        else if (segment_intersectLine(c, d, a, b))
            crossp[tot++] = point_line_intersectLine(a, b, c, d);
    }
    sort(crossp, crossp + tot);
    tot = unique(crossp, crossp + tot) - crossp;
    double ans = 0.0;
    for (int i = 0;i < tot - 1;++i)
    {
        Point tmp = 0.5*(crossp[i] + crossp[i + 1]);
        if (inPolygon(tmp, polygon))
            ans += getDis(crossp[i], crossp[i + 1]);
    }
    return ans;
}

bool inCircle(Point p, Circle a)//点是否在圆内
{
    return judgeZero(getDis(p, a._o) - a._r) <= 0;
}

double dis_toSegment(Point c, Point a, Point b)//点到线段的最短距离
{
    Point d = a - b;
    double tmp = d._y;//得到与ab垂直的向量
    d._y = d._x;
    d._x = -tmp;
    d.toMove(c, atan2(d._y, d._x), 10.0);
    Point intsect = point_line_intersectLine(a, b, c, d);
    if (onSegment(intsect, a, b))
        return getDis(c, intsect);
    return min(getDis(c, a), getDis(c, b));
}

Polygon dividePolygon(Polygon p, Point a, Point b)//直线ab划分多边形
{
    int n = p.size();
    Polygon newp;
    for (int i = 0;i < n;i++)
    {
        Point c = p[i], d = p[(i + 1) % n];
        if (judgeZero(getCross(b, c, a)) > 0)//只记录向量ab左侧的多边形
            newp.push_back(c);
        if (onLine(c, a, b))//防止重复记录点
            newp.push_back(c);
        else if (onLine(d, a, b))
            continue;
        else if (lineIntersect(a, b, c, d))
        {
            Point tmp = point_line_intersectLine(a, b, c, d);
            if (onSegment(tmp, c, d))
                newp.push_back(tmp);
        }
    }
    return newp;//要判断newp.size()>=3
}

double polygonArea(Polygon polygon)//多边形面积
{
    double ans = 0.0;
    for (int i = 1;i < polygon.size() - 1;i++)
        ans += getCross(polygon[i], polygon[i + 1], polygon[0]);
    return fabs(ans) / 2.0;
}

bool circle_inPolygon(Circle c, Polygon p)//圆含于多边形
{
    if (!inPolygon(c._o, p))
        return false;
    for (int i = 0;i < p.size();i++)
    {
        Point a = p[i], b = p[(i + 1) % p.size()];
        if (judgeZero(dis_toSegment(c._o, a, b) - c._r) < 0)
            return false;
    }
    return true;
}

bool polygon_inCircle(Circle c, Polygon p)//多边形含于圆
{
    for (int i = 0;i < p.size();i++)
        if (!inCircle(p[i], c))
            return false;
    return true;
}

bool polygon_intersectCircle(Circle c, Polygon p)//圆含于多边形，多边形含于圆，圆和多边形相交，相切没算进去
{
    if (polygon_inCircle(c, p) || circle_inPolygon(c, p))
        return true;
    for (int i = 0;i < p.size();i++)
    {
        Point a = p[i], b = p[(i + 1) % p.size()];
        if (judgeZero(dis_toSegment(c._o, a, b) - c._r) < 0)//算相切要<=0
            return true;
    }
    return false;
}

int point_line_intersectCircle(Circle c, Point a, Point b, Point p[])//线段ab与圆c的交点
{
    double A = (b._x - a._x)*(b._x - a._x) + (b._y - a._y)*(b._y - a._y);
    double B = 2.0 * ((b._x - a._x)*(a._x - c._o._x) + (b._y - a._y)*(a._y - c._o._y));
    double C = (a._x - c._o._x)*(a._x - c._o._x) + (a._y - c._o._y)*(a._y - c._o._y) - c._r*c._r;
    double deta = B*B - 4.0 * A*C;
    if (judgeZero(deta) < 0) return 0;
    double t1 = (-B - sqrt(deta)) / (2.0 * A);
    double t2 = (-B + sqrt(deta)) / (2.0 * A);
    int tot = 0;
    if (judgeZero(t1) >= 0 && judgeZero(1 - t1) >= 0)
        p[tot++] = Point(a._x + t1*(b._x - a._x), a._y + t1*(b._y - a._y));
    if (judgeZero(t2) >= 0 && judgeZero(1 - t2) >= 0)
        p[tot++] = Point(a._x + t2*(b._x - a._x), a._y + t2*(b._y - a._y));
    if (tot == 2 && p[0] == p[1]) return 1;
    return tot;
}

double sectorArea(Circle c, Point p1, Point p2)//扇形面积op1p2
{
    Point op1 = p1 - c._o, op2 = p2 - c._o;
    double sita = atan2(op2._y, op2._x) - atan2(op1._y, op1._x);
    while (judgeZero(sita) <= 0) sita += 2 * PI;
    while (judgeZero(sita - 2 * PI) > 0) sita -= 2 * PI;
    sita = min(sita, 2 * PI - sita);
    return c._r*c._r*sita / 2.0;
}

double area_triangle_fromCircle(Circle c, Point a, Point b)//圆点，a，b构成的三角形与圆相交的面积
{
    bool flag1 = inCircle(a, c);
    bool flag2 = inCircle(b, c);
    if (flag1 && flag2)//两点在圆内
        return fabs(getCross(a, b, c._o) / 2.0);
    Point p[2];
    int tot = point_line_intersectCircle(c, a, b, p);
    if (flag1) return fabs(getCross(a, p[0], c._o) / 2.0) + sectorArea(c, p[0], b);
    if (flag2) return fabs(getCross(b, p[0], c._o) / 2.0) + sectorArea(c, p[0], a);
    if (tot == 2) return sectorArea(c, p[0], a) + sectorArea(c, p[1], b) + fabs(getCross(p[0], p[1], c._o)) / 2;
    return sectorArea(c, a, b);
}

double area_polygon_intersectCircle(Circle c, Polygon p)//多边形和圆的相交面积
{
    double ans = 0.0;
    for (int i = 0;i < p.size();i++)//trick：只需每相邻的点与圆心相连求相交面积，有向面积可以相互抵消
    {
        Point a, b;
        a = p[i], b = p[(i + 1) % p.size()];
        int flag = judgeZero(getCross(a, b, c._o));
        ans += flag*area_triangle_fromCircle(c, a, b);
    }
    return fabs(ans);
}

Point geometryCenter(Polygon p)//几何中心
{
    double x = 0.0, y = 0.0;
    for (int i = 0;i < p.size();++i)
        x += p[i]._x, y += p[i]._y;
    return Point(x / p.size(), y / p.size());
}

Circle c1, c2;
Polygon p;
double x_l, x_r, y_l, y_r;
double RR, mxx;
int n;

double toCal(double x)
{
    double r = -150, l = 150;
    for (int i = 0;i < n;++i)
    {
        double q = (p[i]._x - x)*(p[(i + 1) % n]._x - x);
        if (judgeZero(q) <= 0)
        {
            Point tmp = point_line_intersectLine(p[i], p[(i + 1) % n], Point(x, -200), Point(x, 200));
            r = max(r, tmp._y);
            l = min(l, tmp._y);
        }
    }
    double ans = 0;
    while (l < r - EPS)//三分y轴
    {
        double p1 = (l + r) / 2, p2 = (p1 + r) / 2;
        double gt1 = area_polygon_intersectCircle(Circle(Point(x, p1), RR), p);
        double gt2 = area_polygon_intersectCircle(Circle(Point(x, p2), RR), p);
        if (gt1 > gt2)
            r = p2;
        else
            l = p1;
    }
    return max(area_polygon_intersectCircle(Circle(Point(x, l), RR), p), area_polygon_intersectCircle(Circle(Point(x, r), RR), p));
}

int main()
{
    mxx = 0;
    x_l = 150, x_r = -150;
    y_l = 150, y_r = -150;
    scanf("%d%lf", &n, &RR);
    for (int i = 0;i < n;++i)
    {
        double x, y;
        scanf("%lf%lf", &x, &y);
        p.push_back(Point(x, y));
        x_l = min(x_l, x);
        x_r = max(x_r, x);
        y_l = min(y_l, y);
        y_r = max(y_r, y);
    }
    while (x_l < x_r - EPS)//三分x轴
    {
        double p1 = (x_l + x_r) / 2;
        double p2 = (p1 + x_r) / 2;
        if (toCal(p1) > toCal(p2))
            x_r = p2;
        else
            x_l = p1;
    }
    printf("%.6f\n", max(toCal(x_r), toCal(x_l)));
    return 0;
}

450行的计算几何，我们考虑用模拟退火来写一写。

退火求得圆心，然后套板子去求最大面积即可。

等学会了计算几何再回来看

#include <bits/stdc++.h>
#define eps 1e-6
typedef double db;
using namespace std;
int n, r;
int sign(db k) {
    if (k > eps) return 1;
    else if (k < -eps) return -1;
    return 0;
}
int cmp(db k1, db k2) {
    return sign(k1 - k2);
}
int inmid(db k1, db k2, db k3) {
    return sign(k1 - k3) * sign(k2 - k3) <= 0;   // k3 在 [k1,k2] 内
}
struct point {
    db x, y;//db就是double
    point(db _x, db _y): x(_x), y(_y) {}
    point operator + (const point &k1) const {
        return (point) {
            k1.x + x, k1.y + y
        };
    }
    point operator - (const point &k1) const {
        return (point) {
            x - k1.x, y - k1.y
        };
    }
    point operator * (db k1) const {
        return (point) {
            x*k1, y*k1
        };
    }
    point operator / (db k1) const {
        return (point) {
            x / k1, y / k1
        };
    }
    int operator == (const point &k1) const {
        return cmp(x, k1.x) == 0 && cmp(y, k1.y) == 0;
    }
    // 逆时针旋转
    point turn(db k1) {
        return (point) {
            x*cos(k1) - y*sin(k1), x*sin(k1) + y*cos(k1)
        };
    }
    point turn90() {
        return (point) {
            -y, x
        };
    }
    bool operator < (const point k1) const {
        int a = cmp(x, k1.x);
        if (a == -1) return 1;
        else if (a == 1) return 0;
        else return cmp(y, k1.y) == -1;
    }
    db abs() {
        return sqrt(x * x + y * y);
    }
    db abs2() {
        return x * x + y * y;
    }
    db dis(point k1) {
        return ((*this) - k1).abs();
    }
    point unit() {
        db w = abs();
        return (point) {
            x / w, y / w
        };
    }
    void scan() {
        double k1, k2;
        scanf("%lf%lf", &k1, &k2);
        x = k1;
        y = k2;
    }
    void print() {
        printf("%.11lf %.11lf\n", x, y);
    }
    db getw() {
        return atan2(y, x);
    }
    point getdel() {
        if (sign(x) == -1 || (sign(x) == 0 && sign(y) == -1)) return (*this) * (-1);
        else return (*this);
    }
    int getP() const {
        return sign(y) == 1 || (sign(y) == 0 && sign(x) == -1);
    }
};
int inmid(point k1, point k2, point k3) {
    return inmid(k1.x, k2.x, k3.x) && inmid(k1.y, k2.y, k3.y);
}
db cross(point k1, point k2) {
    return k1.x * k2.y - k1.y * k2.x;
}
db dot(point k1, point k2) {
    return k1.x * k2.x + k1.y * k2.y;
}
db rad(point k1, point k2) {
    return atan2(cross(k1, k2), dot(k1, k2));
}
point proj(point k1, point k2, point q) { // q 到直线 k1,k2 的投影
    point k = k2 - k1;
    return k1 + k * (dot(q - k1, k) / k.abs2());
}
db disSP(point k1, point k2, point q) {
    point k3 = proj(k1, k2, q);
    if (inmid(k1, k2, k3)) return q.dis(k3);
    else return min(q.dis(k1), q.dis(k2));
}
struct circle {
    point o;
    db r;
    circle(db x, db y, db rr): o(point(x, y)), r(rr) {}
    void scan() {
        o.scan();
        scanf("%lf", &r);
    }
    int inside(point k) {
        return cmp(r, o.dis(k));
    }
};
vector<point> getCL(circle k1, point k2, point k3) { // 沿着 k2->k3 方向给出 , 相切给出两个
    point k = proj(k2, k3, k1.o);
    db d = k1.r * k1.r - (k - k1.o).abs2();
    if (sign(d) == -1) return {};
    point del = (k3 - k2).unit() * sqrt(max((db)0.0, d));
    return {k - del, k + del};
}
db getarea(circle k1, point k2, point k3) {
    // 圆 k1 与三角形 k2 k3 k1.o 的有向面积交
    point k = k1.o;
    k1.o = k1.o - k;
    k2 = k2 - k;
    k3 = k3 - k;
    int pd1 = k1.inside(k2), pd2 = k1.inside(k3);
    vector<point>A = getCL(k1, k2, k3);
    if (pd1 >= 0) {
        if (pd2 >= 0) return cross(k2, k3) / 2;
        return k1.r * k1.r * rad(A[1], k3) / 2 + cross(k2, A[1]) / 2;
    } else if (pd2 >= 0) {
        return k1.r * k1.r * rad(k2, A[0]) / 2 + cross(A[0], k3) / 2;
    } else {
        int pd = cmp(k1.r, disSP(k2, k3, k1.o));
        if (pd <= 0) return k1.r * k1.r * rad(k2, k3) / 2;
        return cross(A[0], A[1]) / 2 + k1.r * k1.r * (rad(k2, A[0]) + rad(A[1], k3)) / 2;
    }
}
vector<point>p;//多边形的点
db cal_area(point P, db R) {//圆心,半径
    circle O = circle(P.x, P.y, R);
    db res = 0;
    int n = p.size();
    for (int i = 0; i < n; i++) {
        res += getarea(O, p[i], p[(i + 1) % n]);
    }
    return res;
}
double dx[4] = {0, 0, 1, -1};
double dy[4] = {1, -1, 0, 0};
double sumx, sumy;
double serch(int n) {
    double T = 50.0;//不能开太大
    point O = point(sumx / n, sumy / n);
    double ans = cal_area(O, r);
    for (int k = 1; k <= 100000; k++) {
        double t = T / k;
        for (int i = 0; i < 4; i++) {
            double tx = O.x + dx[i] * t;
            double ty = O.y + dy[i] * t;
            point tmp = point(tx, ty);
            double area = cal_area(tmp, r);
            if (area > ans) {
                ans = area;
                O = tmp;
            }
        }
    }
    return ans;
}
int main() {
    scanf("%d%d", &n, &r);
    for (int i = 1; i <= n; i++) {
        double x, y;
        scanf("%lf%lf", &x, &y);
        p.push_back(point(x, y));
        sumx += x;
        sumy += y;
    }
    printf("%f\n", serch(n));
    return 0;
}