题解：P17095 [ICPC 2017 Qingdao R] Enclosure of Land and Water

· · 题解

题意简述

起点到直线海岸的距离为 d。陆地与水域中的速度和单位面积价值分别不同。求在时间 t 内从起点出发并返回时，闭合路径所围区域的最大总价值。

解题思路

先定义两个圆弧函数。设弦长为 c，圆心角为 a：

\begin{aligned} G(a) & =\frac{a}{\sin(a/2)} \\ S(c,a) & =\frac{c^2(a-\sin a)}{8\sin^2(a/2)} \end{aligned}

弧长为 cG(a)/2，弓形面积为 S(c,a)。当 a\to0 时，代码改用：

\begin{aligned} G(a) & =2+\frac{a^2}{12}+\frac{7a^4}{2880} \\ S(c,a) & =c^2\left(\frac{a}{12}+\frac{a^3}{360}+\frac{a^5}{10080}\right) \end{aligned}

若曲线完全位于陆地，等周不等式给出一个周长为 vt 的圆。其价值为 Vv^2t^2/(4\pi)。

若曲线经过水域，可先将凹陷部分向外凸化。该操作不增加路程，却不会减小面积，因此最优曲线只与海岸相交两次。

固定同一介质内一段边界的端点与长度，圆弧围出的面积最大。若两段陆地圆弧半径不同，还能在两段之间重新分配少量长度，使总面积增大。因此两段陆地圆弧半径相同，接岸角也相同。

设半径为 r，一段陆地圆弧的半圆心角为 u，海岸处的公共方向角为 \lambda。起点到海岸的距离满足：

d=2r\sin u\sin(\lambda-u)

等式右侧关于 u 与 \lambda-u 对称。固定 r 和 \lambda 后，合法范围内至多产生这两个解。两侧选择不同解时，两段圆弧在起点光滑连接。它们合成一段完整陆地圆弧。两侧选择相同解时，图形关于垂线对称，起点处形成角。因此下面两类覆盖全部最优曲线。

光滑类中，设陆地、水域圆弧的圆心角为 a,b，陆地圆弧半径为 r。两段圆弧共用弦长 c=2r\sin(a/2)。用满全部时间可得：

r=\frac{t}{a/v+\sin(a/2)G(b)/w}

对应价值为：

V\frac{r^2(a-\sin a)}{2}+WS(c,b)

陆地圆弧能经过起点，当且仅当 r(1-\cos(a/2))\ge d。

光滑类的内部极值可以直接求出。令 A=a/2、B=b/2。分别改变两种介质中的弧长，可得加权曲率条件。移动海岸交点，则得到接岸方向条件。两式为：

\begin{aligned} w\cos A+v\cos B & =0 \\ \frac{\sin B}{\sin A} & =\frac{Ww}{Vv} \end{aligned}

消去 B 后得到：

\sin^2 A=\frac{V^2(v^2-w^2)}{w^2(W^2-V^2)}

记等式右侧为 z。若 z\in[0,1]，枚举 A=\arcsin\sqrt{z} 和 \pi-A。再用两条一阶条件与 atan2 求 B，最后检查高度条件。

高度取等号时，起点是完整陆地圆弧的最深点。将该圆弧从起点切开，恰好得到下一类中的两段对称圆弧。其余退化边界也会落入全陆圆或下一类的边界，因此不会遗漏。

对称类中，设海岸交点到起点垂足的距离为 x。每段陆地圆弧的弦长为 \sqrt{d^2+x^2}。设每段陆地圆弧的圆心角为 a，水域圆弧的圆心角为 b，并定义：

\begin{aligned} p & =\frac{G(a)}{v} \\ q & =\frac{G(b)}{w} \end{aligned}

时间条件为 p\sqrt{d^2+x^2}+qx=t。令 e=(t-pd)(t+pd)，其非负解可稳定地写成：

x=\frac{e}{tq+p\sqrt{e+d^2q^2}}

陆地与水域面积分别为：

\begin{aligned} A_L & =dx+2S\left(\sqrt{d^2+x^2},a\right) \\ A_W & =S(2x,b) \end{aligned}

这里 0\le a\le\pi，且 pd\le t。当 b\to2\pi 时，直接计算会出现 0\times\infty。此时改用极限：

\begin{aligned} A_L & =2S(d,a) \\ A_W & =\frac{w^2(t-pd)^2}{4\pi} \end{aligned}

为了说明内部极值的约束，令 A=a/2、B=b/2，并定义：

\begin{aligned} \rho & =\frac{Vv}{Ww} \\ D & =4\sin^2A-\rho^2\sin^2B \end{aligned}

对 A、B 和 x 分别作一阶变分。消去拉格朗日乘子后，任意内部极值都满足：

\begin{aligned} \sqrt{D}-\rho\sin B\cot A-2\frac{v}{w}\cos B & =0 \\ \frac{t}{d}\sqrt{D}-\frac{4A}{v}-\frac{2\rho B}{w} & =0 \end{aligned}

第一式来自弧长分配与交点移动，第二式就是时间条件。边界则是 A=0、A=A_{\max}、B=0 或 B=\pi。因此数值优化只需覆盖这些边界与内部驻点，无需假设目标函数单峰。

代码先计算 200\times200 网格，并保留八邻域局部最高点。网格同时包含全部边界和退化端点。随后从每个候选点做 80 轮二维缩步搜索。网格负责定位不同驻点附近的候选区域，连续缩步负责提高局部极值的计算精度。

设网格边数为 N，缩步轮数为 K，局部最高点数量为 M。每组数据的计算量为 O(N^2+KM)，空间复杂度为 O(N^2)。

参考代码

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

using ld=long double;
const int M=200;
const int N=M+5;
const double pi=acos(-1);
const double neg=-1e100;
double d,v,w,vl,vw,t,lim;
double val[N][N];
double coef(double a)
{
    return a<1e-4?2+a*a/12+7*a*a*a*a/2880:a/sin(a/2);
}
double segment(double c,double a)
{
    if(a<1e-4)return c*c*(a/12+a*a*a/360+a*a*a*a*a/10080);
    double s=sin(a/2);
    return c*c*(a-sin(a))/(8*s*s);
}
double corner(double a,double b)
{
    double p=coef(a)/v;
    if(p*d>t+1e-12)return neg;
    double al,aw;
    if(2*pi-b<1e-10)
    {
        double r=max(0.0,t-p*d);
        al=2*segment(d,a);
        aw=w*w*r*r/(4*pi);
    }
    else
    {
        double q=coef(b)/w,e=max(0.0,(t-p*d)*(t+p*d)),x=e/(t*q+p*sqrt(e+d*d*q*q));
        al=d*x+2*segment(hypot(d,x),a);
        aw=segment(2*x,b);
    }
    return vl*al+vw*aw;
}
double smooth(double a,double b)
{
    if(a<1e-10||2*pi-b<1e-10)return neg;
    double s=sin(a/2),r=t/(a/v+s*coef(b)/w);
    if(r*(1-cos(a/2))+1e-12<d)return neg;
    double al=r*r*(a-sin(a))/2,aw=segment(2*r*s,b);
    return vl*al+vw*aw;
}
double smooth_best()
{
    ld lv=v,lw=w,lvl=vl,lvw=vw,den=lw*lw*(lvw*lvw-lvl*lvl);
    if(den==0)return neg;
    ld num=lvl*lvl*(lv*lv-lw*lw),z=num/den;
    if(z<=0||z>1+1e-12)return neg;
    z=max(0.0L,min(z,1.0L));
    ld x=asinl(sqrtl(z)),ans=neg;
    for(int i=0;i<2;i++)
    {
        ld a=i?pi-x:x;
        if(i&&fabsl(a-x)<1e-15)continue;
        ld sb=lvw*lw/(lvl*lv)*sinl(a),cb=-lw/lv*cosl(a);
        if(fabsl(sb*sb+cb*cb-1)>1e-8)continue;
        ld b=atan2l(sb,cb);
        ans=max(ans,ld(smooth(2*a,2*b)));
    }
    return ans;
}
double opt(double la,double lb,double (*fun)(double,double))
{
    double da=la/M,db=lb/M,ans=0;
    for(int i=0;i<=M;i++)
    {
        for(int j=0;j<=M;j++)
        {
            val[i][j]=fun(i*da,j*db);
            ans=max(ans,val[i][j]);
        }
    }
    for(int i=0;i<=M;i++)
    {
        for(int j=0;j<=M;j++)
        {
            if(val[i][j]<neg/2)continue;
            bool ok=1;
            for(int k=0;k<9;k++)
            {
                int x=i+k/3-1,y=j+k%3-1;
                if(x<0||x>M||y<0||y>M)continue;
                if(val[x][y]>val[i][j]||(val[x][y]==val[i][j]&&make_pair(x,y)<make_pair(i,j)))ok=0;
            }
            if(!ok)continue;
            double a=i*da,b=j*db,sa=da,sb=db,cur=val[i][j];
            int cnt=80;
            while(cnt--)
            {
                double na=a,nb=b,bst=cur;
                for(int k=0;k<9;k++)
                {
                    int x=k/3-1,y=k%3-1;
                    double aa=max(0.0,min(a+x*sa,la)),bb=max(0.0,min(b+y*sb,lb)),res=fun(aa,bb);
                    if(res>bst)
                    {
                        bst=res;
                        na=aa;
                        nb=bb;
                    }
                }
                if(bst>cur)
                {
                    a=na;
                    b=nb;
                    cur=bst;
                }
                else
                {
                    sa/=2;
                    sb/=2;
                }
            }
            ans=max(ans,cur);
        }
    }
    return ans;
}
double solve()
{
    double ans=vl*v*v*t*t/(4*pi);
    ans=max(ans,smooth_best());
    if(2*d/v>t)return ans;
    if(pi*d/v<=t)lim=pi;
    else
    {
        double l=0,r=pi;
        for(int i=0;i<80;i++)
        {
            double mid=(l+r)/2;
            if(coef(mid)*d/v<=t)l=mid;
            else r=mid;
        }
        lim=l;
    }
    return max(ans,opt(lim,2*pi,corner));
}
int main()
{
    ios::sync_with_stdio(false);
    cin.tie(nullptr);
    int T;
    cin>>T;
    cout<<fixed<<setprecision(3);
    while(T--)
    {
        cin>>d>>v>>w>>vl>>vw>>t;
        cout<<solve()<<'\n';
    }
    return 0;
}