Flight Safety (BZOJ1020)

题意

世界上有很多岛屿，可以看成是 C 个边形 ($C\leq20$，每个点数 M 不超过 30，可以非凸)。一架飞机沿一条 N 点折线飞行 ($N\leq20$)，求航线上一个点，最大化其到最近岛屿的距离。给出的所有点都是 $[-1000,10000]$内的整数坐标。要求最终答案精度在 $10^{-2}$以下 (UVA 精度 $10^{-3}$)。

分析

首先我想，既然航线被划分成了 “在岛上”-“在海上”-“在岛上”-… 的交替模式，那么我们对每个在海上的段求最远距离不是很容易吗？取个中点不就行了。

反例 (中点显然不是最优的)：

那么求截得这段海上航线的两条线段的角平分线与航线的交点不就是最远点吗？

反例：

额。。。这个里面的花，我们或许只能用更奇怪的方法解决这道题了。

爬山 (90pts)

爬山算法是求解怪异函数的最优化问题的常用方法。

如果将航线上每个点到最近岛屿的距离看成有关飞行距离的函数，那么我们会看到一个非常可爱的函数图像。这个图像有如下性质：斜率的绝对值不超过 1。原因是：从任意一处出发，向任意方向飞行，你与最近岛屿的距离变化速度不可能超过你飞行的速度。

因此，我们证明了这个函数是不会突变的。这就意味着我们可以让很多个小人从不同位置爬这个图像，他们不会因为函数突变而摔死或者爬不动 (即若函数不连续则爬山算法不可解极值)。这样，最终爬地最高的小人就代表着函数的极值。

这种算法的基本步骤是：

• 1. 按一定间隔放置小人
• 2. 每个小人有一个特定的步长 s
• 3. 小人向着函数增长的方向爬 s 的长度
• 4.s 变小

最终每个小人都会稳定在一座山峰处。

但是我们在准备实现这个算法之前需要注意如下几个问题：

• 要多少个小人可以保证几乎每个峰都被找到，或没找到的不会影响精度
• 小人的初始步长 s
• 小人步长的减少率 v
• 小人要爬多少步
• 以上的算法如果可以 100% 找到最优解，时间复杂度是否允许。

首先，我们可以证明，两个小人间距小于 $10^{-2}$才一定可以找到最优解，而这和暴力找点没有任何区别，时间复杂度都不足以通过此题。因此我们可以冒险将间距调大，爬山步数调小，减少率设为 2/3，并且面向数据编程，就可以在洛谷获得 90 分的好成绩。

二分 (100pts)

我并没有亲身尝试这个恐怖的算法，但是这是复杂度严格 $O(NMClogA)$的 (A 为答案大小)。这是一个非常优秀的复杂度，谁看了都想打，但是不巧，看看它是怎么做的你就不想打了。

• 二分距离 x
• 将每个多边形的边向外平移 x，形成多个由圆弧和线段组成的形状
• 求航线与这些形状的交点，并对交点离散化后逆差分，求出每一段被形状覆盖的次数
• 如果存在一段没有被覆盖，则 x 扩大，反之减小。

看到了那个 “由圆弧和线段组成的形状” 了没？这东西原本就及其的恐怖，还需要求和线段的交点，实在是太难打。

迭代 (100pts)

于是我们继续利用爬山的那种偷懒的思想。

分析爬山不可行的原因，是有很多没用的人在爬。国家供养不起那么多废物，于是就 TLE 了。所以我们要想办法清理掉更多的废物，保留最有可能爬上最高峰的人。

于是迭代算法就应运而生了。

对于一条线段，它到最近岛屿的距离一定是在 $[D-\frac{L}{2},D+\frac{L}{2}]$的区间内的 (D 为其中点到最近岛屿的距离，L 为其长度)。如果这个区间和 $[X,+\infty]$(X 为已知的最远距离) 没有交集，或者 L 已经非常短，可以认为这条线段就是个废物，否则，我们可以将其拆成前半段和后半段递归处理。

因此我们可以用一个队列 (甚至有限队列) 维护当前所有有用的线段，每次生成两个线段，将有用的保留。这样，我们只会用有用的线段进行类似爬山的操作，这样就可以很快获得答案。

代码

爬山

#pragma GCC optimize(3)
#include <iostream>
#include <cstring>
#include <cstdio>
#include <cmath>
#define eps 1e-8
#define oo 1e5
#define MX 102

using namespace std;

struct vec
{
    double x,y;
    void input(){scanf("%lf%lf",&x,&y);}
    void output(){printf("%.3lf %.3lf\n",x,y);}
    vec (){}
    vec (const double& a,const double& b){x=a,y=b;}
    vec operator + (const vec& t)const{return vec(x+t.x,y+t.y);}
    vec operator - (const vec& t)const{return vec(x-t.x,y-t.y);}
    vec operator * (const double& t)const{return vec(x*t,y*t);}
    vec operator / (const double& t)const{return vec(x/t,y/t);}
    double len(){return sqrt(x*x+y*y);}
    vec normal(){return *this/(len()+eps);}
};
double dot(const vec& a,const vec& b){return a.x*b.x+a.y*b.y;}
double crs(const vec& a,const vec& b){return a.x*b.y-a.y*b.x;}
bool seg_seg_int(const vec& ap,const vec& aq,const vec& bp,const vec& bq)
{
    if( min(ap.x,aq.x)>max(bp.x,bq.x)||
        min(bp.x,bq.x)>max(ap.x,aq.x)||
        min(ap.y,aq.y)>max(bp.y,bq.y)||
        min(bp.y,bq.y)>max(ap.y,aq.y))return 0;
    if( crs(bq-aq,ap-aq)*crs(bp-aq,ap-aq)>0||
        crs(ap-bq,bp-bq)*crs(aq-bq,bp-bq)>0)return 0;
    return 1;
}
double p_to_seg(const vec& a,const vec& p,const vec& q)
{
    vec v=q-p;
    double t=dot(a-p,v)/dot(v,v);
    if(t<0||t>1)return min((a-p).len(),(a-q).len());
    vec b=p+(v*t);
    return (a-b).len();
}
bool p_in_poly(const vec& a,const vec* p,int n)
{
    int cnt=0;
    vec far=a+vec(1,+oo);
    for(int i=1;i<=n;i++)cnt+=seg_seg_int(a,far,p[i],p[i%n+1]);
    return (cnt&1);
}
int C,N,M[MX];
vec line[MX];
vec poly[MX][MX];

double mindis(const vec& a)
{
    double ans=+oo;
    for(int i=1;i<=N;i++)
    {
        if(M[i]>=3&&p_in_poly(a,poly[i],M[i]))return 0;
        if(M[i]>=2)for(int j=1;j<=M[i];j++)ans=min(ans,p_to_seg(a,poly[i][j],poly[i][j%M[i]+1]));
        else ans=min(ans,(poly[i][1]-a).len());
    }
    return ans;
}

void input()
{
    scanf("%d%d",&N,&C);
    for(int i=1;i<=C;i++)line[i].input();
    for(int i=1;i<=N;i++)
    {
        scanf("%d",&M[i]);
        for(int j=1;j<=M[i];j++)poly[i][j].input();
    }
}

double cal_max(const vec& p,const vec& q)
{
    vec v=q-p;
    double len=v.len(),ret=0,di=max(1e-2/len,2e-3);
    for(double i=0;i<=1;i+=di)
    {
        double now=i,stp=di/2,lf=mindis(p+v*now);
        for(int j=1;j<=21;j++)
        {
            double tr=min(1.0,now+stp);
            double nowf=lf,nxtf=mindis(p+v*tr);
            if(nxtf>nowf)now=tr,lf=nxtf;
            else now=max(0.0,now-stp),lf=mindis(p+v*now);
            if(nowf+stp*1.857142857<=ret)break;
            stp*=0.65;
        }
        ret=max(ret,mindis(p+v*now));
    }
    return ret;
}

void work()
{
    double ans=0;
    for(int i=1;i<C;i++)ans=max(ans,cal_max(line[i],line[i+1]));
    printf("%.2lf\n",ans);
}

int main()
{
    input();
    work();
    return 0;
}

迭代

#pragma GCC optimize(3)
#include <iostream>
#include <cstring>
#include <cstdio>
#include <cmath>
#define eps 1e-8
#define oo 1e5
#define MX 102

using namespace std;

struct vec
{
    double x,y;
    void input(){scanf("%lf%lf",&x,&y);}
    void output(){printf("%.3lf %.3lf\n",x,y);}
    vec (){}
    vec (const double& a,const double& b){x=a,y=b;}
    vec operator + (const vec& t)const{return vec(x+t.x,y+t.y);}
    vec operator - (const vec& t)const{return vec(x-t.x,y-t.y);}
    vec operator * (const double& t)const{return vec(x*t,y*t);}
    vec operator / (const double& t)const{return vec(x/t,y/t);}
    double len(){return sqrt(x*x+y*y);}
    vec normal(){return *this/(len()+eps);}
};
struct itv
{
    vec l,r;
    itv (){}
    itv (const vec& a,const vec& b){l=a,r=b;}
    vec mid(){return (l+r)/2.0;}
    double len(){return (r-l).len();}
};
double dot(const vec& a,const vec& b){return a.x*b.x+a.y*b.y;}
double crs(const vec& a,const vec& b){return a.x*b.y-a.y*b.x;}
bool seg_seg_int(const vec& ap,const vec& aq,const vec& bp,const vec& bq)
{
    if( min(ap.x,aq.x)>max(bp.x,bq.x)||
        min(bp.x,bq.x)>max(ap.x,aq.x)||
        min(ap.y,aq.y)>max(bp.y,bq.y)||
        min(bp.y,bq.y)>max(ap.y,aq.y))return 0;
    if( crs(bq-aq,ap-aq)*crs(bp-aq,ap-aq)>0||
        crs(ap-bq,bp-bq)*crs(aq-bq,bp-bq)>0)return 0;
    return 1;
}
double p_to_seg(const vec& a,const vec& p,const vec& q)
{
    vec v=q-p;
    double t=dot(a-p,v)/dot(v,v);
    if(t<0||t>1)return min((a-p).len(),(a-q).len());
    vec b=p+(v*t);
    return (a-b).len();
}
bool p_in_poly(const vec& a,const vec* p,int n)
{
    int cnt=0;
    vec far=a+vec(1,+oo);
    for(int i=1;i<=n;i++)cnt+=seg_seg_int(a,far,p[i],p[i%n+1]);
    return (cnt&1);
}
int C,N,M[MX];
vec line[MX];
vec poly[MX][MX];

double mindis(const vec& a)
{
    double ans=+oo;
    for(int i=1;i<=N;i++)
    {
        if(M[i]>=3&&p_in_poly(a,poly[i],M[i]))return 0;
        if(M[i]>=2)for(int j=1;j<=M[i];j++)ans=min(ans,p_to_seg(a,poly[i][j],poly[i][j%M[i]+1]));
        else ans=min(ans,(poly[i][1]-a).len());
    }
    return ans;
}

void input()
{
    scanf("%d%d",&N,&C);
    for(int i=1;i<=C;i++)line[i].input();
    for(int i=1;i<=N;i++)
    {
        scanf("%d",&M[i]);
        for(int j=1;j<=M[i];j++)poly[i][j].input();
    }
}

itv use[1000001];
void work()
{
    double best=0;
    int h=0,t=1;
    for(int i=1;i<C;i++)use[++h]=itv(line[i],line[i+1]);
    for(int i=1;i<=300000&&h>=t;i++)
    {
        itv now=use[t++],nxt1,nxt2;
        best=max(best,mindis(now.mid()));
        nxt1=itv(now.l,now.mid());
        nxt2=itv(now.mid(),now.r);
        if(nxt1.len()>=1e-4&&mindis(nxt1.mid())+nxt1.len()/2.0>best)use[++h]=nxt1;
        if(nxt2.len()>=1e-4&&mindis(nxt2.mid())+nxt2.len()/2.0>best)use[++h]=nxt2;
    }
    printf("%.6lf\n",best);
}

int main()
{
    int t;
    scanf("%d",&t);
    for(int i=1;i<=t;i++)
    {
        input();
        work();
    }
    return 0;
}