NKOJ 3966 （ZJOI 2008）瞭望塔（半平面交/模拟退火）

P3966【ZJOI2008】瞭望塔

问题描述

致力于建设全国示范和谐小村庄的H村村长dadzhi，决定在村中建立一个瞭望塔，以此加强村中的治安。

我们将H村抽象为一维的轮廓。如下图所示

我们可以用一条山的上方轮廓折线(x1, y1), (x2, y2), …. (xn, yn)来描述H村的形状，这里x1 < x2 < …< xn。瞭望塔可以建造在[x1, xn]间的任意位置, 但必须满足从瞭望塔的顶端可以看到H村的任意位置。可见在不同的位置建造瞭望塔，所需要建造的高度是不同的。为了节省开支，dadzhi村长希望建造的塔高度尽可能小。

请你写一个程序，帮助dadzhi村长计算塔的最小高度。

输入格式

第一行包含一个整数n，表示轮廓折线的节点数目。
接下来第一行n个整数, 为x1 ~ xn.
第三行n个整数，为y1 ~ yn。

输出格式

仅包含一个实数，为塔的最小高度，精确到小数点后三位。

样例输入 1

6
1 2 4 5 6 7
1 2 2 4 2 1

样例输出 1

1.000

样例输入 2

4
10 20 49 59
0 10 10 0

样例输出 2

14.500

提示

对于60%的数据， N ≤ 60；
对于100%的数据， N ≤ 300，输入坐标绝对值不超过106，注意考虑实数误差带来的问题。

---

此题显然是要把轮廓线上每一条直线拿出来做半平面交得到一个下凸壳，然后找两条轮廓上同一x坐标竖直距离最近的两个点，由于上下都是直线，那么显然只可能在端点取到，因此只需要枚举每一个端点即可。

另外，本题可以直接模拟退火+爬山，只需要调好参数即可，算上方凸壳的y坐标只需要取每条直线的y坐标最小值。
需要注意的是，最好在退火结束的位置多爬几次山，减小精度误差。

---

代码（半平面交）：

#include<stdio.h>
#include<iostream>
#include<algorithm>
#include<cstring>
#define N 1005
using namespace std;
struct node{double x,y;}Point[N],Line[N],S[N],T[N];
int n,top;
bool operator<(node a,node b)
{
    if(a.x==b.x)return a.y>b.y;
    return a.x<b.x;
}
node Cal(node a,node b)
{
    double k=(b.y-a.y)/(b.x-a.x);
    double t=b.y-k*b.x;
    return (node){k,t};
}
node Intersection(node a,node b)
{
    double x=(b.y-a.y)/(a.x-b.x);
    return (node){x,b.x*x+b.y};
}
bool judge(node p,node l)
{
    double y=l.x*p.x+l.y;
    return p.y<y;
}
void Halfplane()
{
    int i,j;
    for(i=1;i<n;i++)
    {
        if(Line[i].x==Line[i-1].x)continue;
        while(top>1&&judge(Intersection(S[top],S[top-1]),Line[i]))top--;
        S[++top]=Line[i];
    }
}
int main()
{
    int i,j;double ans=1e10;
    scanf("%d",&n);
    for(i=1;i<=n;i++)scanf("%lf",&Point[i].x);
    for(i=1;i<=n;i++)scanf("%lf",&Point[i].y);
    for(i=1;i<n;i++)Line[i]=Cal(Point[i],Point[i+1]);
    sort(Line+1,Line+n);
    Halfplane();
    for(i=1;i<top;i++)T[i]=Intersection(S[i],S[i+1]);
    for(i=1;i<n;i++)Line[i]=Cal(Point[i],Point[i+1]);
    Point[n+1].x=1e10;
    for(i=1,j=1;i<top;i++)
    {
        while(T[i].x>=Point[j].x)j++;
        ans=min(ans,T[i].y-(Line[j-1].x*T[i].x+Line[j-1].y));
    }
    T[top].x=1e10;
    for(i=1,j=1;i<=n;i++)
    {
        while(Point[i].x>=T[j].x)j++;
        ans=min(ans,S[j].x*Point[i].x+S[j].y-Point[i].y);
    }
    printf("%.3lf",ans);
}

---

代码（模拟退火）：

#include<stdio.h>
#include<iostream>
#include<algorithm>
#include<cstring>
#include<cmath>
#include<cstdlib>
#define N 505
using namespace std;
struct node{double x1,x2,k,b;}L[N];
int n,X[N],Y[N];
double ans=1e20,Now,New;
double Rand()
{return rand()%10000/10000.0;}
node GetLine(int x1,int y1,int x2,int y2)
{
    double lx=1.0*x1,rx=1.0*x2;
    double ly=1.0*y1,ry=1.0*y2;
    double k=(ry-ly)/(rx-lx);
    double b=ry-k*rx;
    return (node){lx,rx,k,b};
}
double f(int p,double x)
{return L[p].k*x+L[p].b;}
double F(double x)
{
    if(x<X[1]||x>X[n])return 1e20;
    double y,Max=-1e20;
    for(int i=1;i<n;i++)
    {
        y=f(i,x);
        Max=max(Max,y);
    }
    for(int i=1;i<n;i++)
    if(x>=L[i].x1&&x<=L[i].x2)
    {
        Max=Max-f(i,x);
        break;
    }
    if(ans>Max)ans=Max;
    return Max;
}
void SA(double T)
{
    Now=X[1]+rand()%(X[n]-X[1]);
    while(T>0.0001)
    {
        New=Now+T*(2*Rand()-1);
        double d=F(Now)-F(New);
        if(d>0||exp(d/T)>Rand())Now=New;
        T*=0.9;
    }
    for(int i=1;i<=500;i++)
    {
        New=Now+T*(2*Rand()-1);
        if(F(New)<F(Now))Now=New;
    }
}
int main()
{
    srand(19260817);
    int i,j,k,x,y;
    scanf("%d",&n);
    for(i=1;i<=n;i++)scanf("%d",&X[i]);
    for(i=1;i<=n;i++)scanf("%d",&Y[i]);
    for(i=1;i<n;i++)L[i]=GetLine(X[i],Y[i],X[i+1],Y[i+1]);
    for(i=1;i<=50;i++)SA(2000000);
    printf("%.3lf",ans);
}