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

· · 题解

题意简述

一条直线把平面分成陆地与水域，起点在陆地一侧，距海岸为 d。两种区域中的速度和单位面积价值不同。求从起点出发并在时间 t 内返回时，闭合曲线所围区域的最大价值。

解题思路

固定一段曲线的端点与长度时，圆弧围出的面积最大。把凹陷部分向外推也不会增加路程。因此，最优曲线可以取成凸曲线；若它进入水域，就只与海岸相交两次，并且每种介质中的边界都是圆弧。

全文使用圆心角的一半。设圆弧的弦长为 c，半圆心角为 x，定义：

F(x)=\frac{x}{\sin x}

该圆弧的长度为 cF(x)，弓形面积为：

A(c,x)=\frac{c^2(2x-\sin2x)}{8\sin^2x}

当 x 很小时，直接计算会出现消去误差。代码改用泰勒展开：

\begin{aligned} F(x) & =1+\frac{x^2}{6}+\frac{7x^4}{360} \\ A(c,x) & =c^2\left(\frac{x}{6}+\frac{x^3}{45}+\frac{x^5}{315}\right) \end{aligned}

完全不进入水域时，由等周不等式，最优曲线是周长为 vt 的圆。它可以放在陆地一侧，并让起点落在圆周上，价值为：

\frac{Vv^2t^2}{4\pi}

下面只考虑进入水域的曲线。两段陆地圆弧的半径必须相同，否则可以把少量弧长从曲率较小的一段移给另一段，从而增加总面积。它们在海岸处的切线方向也相同。

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

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

右侧交换 u 与 \lambda-u 后不变，且在合法区间内只有这两个可能值。两侧取不同值时，它们在起点处光滑相接，合成一段完整圆弧；取相同值时，曲线关于起点到海岸的垂线对称，并在起点处形成拐角。这两类已经覆盖全部最优形状。

先处理光滑相接的情形。设完整陆地圆弧与水域圆弧的半圆心角分别为 a,b，陆地圆弧半径为 r。公共弦长为 2r\sin a。用满全部时间可得：

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

两部分面积分别为：

\begin{aligned} A_L & =\frac{r^2(2a-\sin2a)}{2} \\ A_W & =A(2r\sin a,b) \end{aligned}

陆地圆弧能够经过起点的条件是：

r(1-\cos a)\ge d

内部最优点还满足两条一阶条件。第一条来自两种介质之间的弧长分配，第二条来自海岸交点的移动：

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

消去 b 后得到：

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

若右侧属于 [0,1]，只需检查 a=\arcsin\sqrt z 与 a=\pi-\arcsin\sqrt z。再由两条一阶条件用 atan2 唯一确定 b，最后检查经过起点的条件。光滑类的边界要么退化为全陆圆，要么恰好让起点成为陆地圆弧最深处，此时也属于下面的对称类。

再处理起点有拐角的情形。设海岸交点到起点垂足的距离为 x，两种圆弧的半圆心角仍记为 a,b。每段陆地弦长为 \sqrt{d^2+x^2}，于是时间方程为：

\frac{2F(a)}{v}\sqrt{d^2+x^2}+\frac{2F(b)}{w}x=t

记 p=2F(a)/v、q=2F(b)/w。方程的非负根可以稳定地写成：

x=\frac{(t-pd)(t+pd)}{tq+p\sqrt{(t-pd)(t+pd)+d^2q^2}}

它要求 pd\le t。为了保证陆地圆弧在抵达交点前不再次穿过海岸，还需要 0\le a\le\pi/2。两种区域的面积为：

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

当 b\to\pi 时，水域圆弧趋于与海岸相切的整圆，直接代入会形成 0\times\infty。剩余水域路程为 w(t-pd)，所以水域面积取极限：

\frac{w^2(t-pd)^2}{4\pi}

对称类的目标函数不保证单峰。代码在 [0,a_{\max}]\times[0,\pi] 上计算 200\times200 网格，收集所有八邻域局部最高点，再从每个候选点做二维缩步搜索。网格包含四条边界和所有退化端点；连续搜索把两个步长反复减半，足以满足三位小数的精度要求。

设网格边数为 N，局部候选数为 M，每个候选缩步 K 次。每组数据的时间复杂度为 O(N^2+MK)，空间复杂度为 O(N^2)。

正确性证明

凸化不会增加曲线长度，也不会减少任何一种区域的面积。固定端点和弧长后，圆弧围出的面积最大。因此至少存在一条最优曲线，只与海岸相交两次，并在每种介质内由圆弧组成。

同一介质中的两段陆地圆弧必须具有相同半径和相同接岸方向。距离方程关于两个半圆心角对称，所以两侧只能取同一个解，或分别取一对互换的解。前者正是对称拐角类，后者正是完整圆弧类，分类没有遗漏。

完整圆弧类的内部最优点必须满足弧长分配与交点移动的一阶条件。代码枚举消元后至多两个 a，据此恢复 b 并检查几何可行性。它的全部边界已经包含在全陆圆或对称类中，因此该类最优值必被检查。

对称类中，给定 a,b 后，时间方程左侧随 x 严格递增，所以公式求出的非负根是唯一合法的 x。面积公式分别准确计算矩形、两段陆地弓形和水域弓形；b=\pi 的极限也被单独处理。网格覆盖完整参数闭区域，随后对每个局部最高点细化，因而数值搜索覆盖该类的全部极值与边界。

程序最终取全陆圆、完整圆弧类和对称拐角类三者的最大值。三类覆盖所有最优曲线，所以输出即为最大总价值。

参考代码

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

using ld=long double;
const int M=200;
const int N=M+5;
const ld pi=acosl(-1);
const ld neg=-1e100L;
ld d,v,w,vl,vw,t,lim;
ld f[N][N];
ld arc(ld x)
{
    if(x<1e-5)return 1+x*x/6+7*x*x*x*x/360;
    return x/sinl(x);
}
ld segment(ld c,ld x)
{
    if(x<1e-5)return c*c*(x/6+x*x*x/45+x*x*x*x*x/315);
    ld s=sinl(x);
    return c*c*(2*x-sinl(2*x))/(8*s*s);
}
ld calc(ld a,ld b)
{
    ld p=2*arc(a)/v;
    if(p*d>t+1e-15)return neg;
    ld al,aw;
    if(pi-b<1e-12)
    {
        ld r=max(0.0L,t-p*d);
        al=2*segment(d,a);
        aw=w*w*r*r/(4*pi);
    }
    else
    {
        ld q=2*arc(b)/w;
        ld e=max(0.0L,(t-p*d)*(t+p*d));
        ld x=e/(t*q+p*sqrtl(e+d*d*q*q));
        al=d*x+2*segment(hypotl(d,x),a);
        aw=segment(2*x,b);
    }
    return vl*al+vw*aw;
}
ld smooth(ld a,ld b)
{
    if(a<1e-12||pi-b<1e-12)return neg;
    ld s=sinl(a);
    ld r=t/(2*a/v+2*s*arc(b)/w);
    if(r*(1-cosl(a))+1e-15<d)return neg;
    ld al=r*r*(2*a-sinl(2*a))/2;
    ld aw=segment(2*r*s,b);
    return vl*al+vw*aw;
}
ld smooth_best()
{
    ld den=w*w*(vw*vw-vl*vl);
    if(fabsl(den)<1e-30)return neg;
    ld z=vl*vl*(v*v-w*w)/den;
    if(z<=0||z>1+1e-15)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-18)continue;
        ld sb=vw*w*sinl(a)/(vl*v);
        ld cb=-w*cosl(a)/v;
        if(fabsl(sb*sb+cb*cb-1)>1e-10)continue;
        ld b=atan2l(sb,cb);
        ans=max(ans,smooth(a,b));
    }
    return ans;
}
ld corner_best()
{
    ld da=lim/M,db=pi/M,ans=0;
    for(int i=0;i<=M;i++)
    {
        for(int j=0;j<=M;j++)
        {
            f[i][j]=calc(i*da,j*db);
            ans=max(ans,f[i][j]);
        }
    }
    vector<pair<int,int>> q;
    for(int i=0;i<=M;i++)
    {
        for(int j=0;j<=M;j++)
        {
            if(f[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(f[x][y]>f[i][j]||(f[x][y]==f[i][j]&&make_pair(x,y)<make_pair(i,j)))ok=0;
            }
            if(ok)q.push_back({i,j});
        }
    }
    for(auto [x,y]:q)
    {
        ld a=x*da,b=y*db,sa=da,sb=db,cur=f[x][y];
        for(int k=0;k<90;k++)
        {
            ld na=a,nb=b,nxt=cur;
            for(int i=-1;i<=1;i++)
            {
                for(int j=-1;j<=1;j++)
                {
                    ld aa=max(0.0L,min(lim,a+i*sa));
                    ld bb=max(0.0L,min(pi,b+j*sb));
                    ld res=calc(aa,bb);
                    if(res>nxt)
                    {
                        nxt=res;
                        na=aa;
                        nb=bb;
                    }
                }
            }
            if(nxt>cur)
            {
                a=na;
                b=nb;
                cur=nxt;
            }
            else
            {
                sa/=2;
                sb/=2;
            }
        }
        ans=max(ans,cur);
    }
    return ans;
}
ld solve()
{
    ld 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/2;
    else
    {
        ld l=0,r=pi/2;
        for(int i=0;i<90;i++)
        {
            ld mid=(l+r)/2;
            if(2*arc(mid)*d/v<=t)l=mid;
            else r=mid;
        }
        lim=l;
    }
    return max(ans,corner_best());
}
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<<double(solve())<<'\n';
    }
    return 0;
}