题解：P16626 [GKS 2017 #D] Trash

· · 题解

题意简述

在 (0,0) 与 (P,0) 之间选择一条开口向下的抛物线。圆形垃圾的中心沿其运动。求不碰到天花板或障碍点时的最大半径。

解题思路

令 t=-a，则 t\ge0，轨迹可以写成：

f(x)=tx(P-x)

先二分答案 R。抛物线的最高点在 x=P/2，因此不碰天花板等价于：

0\le t<\frac{4(H-R)}{P^2}

下面判断这个范围内是否存在不会碰到任何障碍的 t。

固定障碍点 (X_i,Y_i)，考虑以它为圆心、R 为半径的开圆。对于圆内满足 0<x<P 的点 (x,y)，经过它的轨迹参数唯一确定为 t=y/(x(P-x))。

开圆与竖条 0<x<P 的交集连通，上式又是连续函数。因此，会穿过这个开圆的所有 t 构成一个区间。轨迹经过障碍点时的参数为：

t_i=\frac{Y_i}{X_i(P-X_i)}

这个区间一定包含 t_i。于是可以分别在 t_i 左右二分，求出该障碍对应的禁用区间。若 t_i 超出天花板给出的范围，只需先判断范围右端是否在禁用区间内。

还需判断给定的 t 是否穿过障碍圆。设轨迹上横坐标为 z 的点到障碍的距离平方为：

D_i(z)=(z-X_i)^2+(tz(P-z)-Y_i)^2

极值只可能在区间端点或驻点处。令 F_i(z)=D_i'(z)/2，展开得到三次方程：

F_i(z)=2t^2z^3-3t^2Pz^2+(t^2P^2+2tY_i+1)z-tPY_i-X_i

求出它在 [0,P] 内的所有实根，代回 D_i 取最小值即可。代码用三次方程求根公式处理一般情况。

当 tP<1 时，直接套公式会出现除以很小的 t。此时有：

F_i'(z)=1+2tY_i+t^2(P^2-6Pz+6z^2)>\frac{1}{2}

因此 F_i 严格递增且只有一个根，可以直接二分。

收集所有禁用区间的端点，并排序去重。这些端点把可选参数范围分成若干段。若存在合法参数，那么它要么是某个端点，要么位于某段内部；后一种情况检查该段中点即可。注意障碍允许相切，所以禁用区间是开区间，端点也必须检查。

设二分轮数为 I=70。每次判定需要 O(IN+N^2) 的时间，总时间复杂度为 O(I^2N+IN^2)，空间复杂度为 O(N)。

参考代码

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

const int N=15;
const int K=70;
const int M=60;
const long double PI=acosl(-1);
int n;
long double p,h,x[N],y[N];
long double get_value(int id,long double t,long double z)
{
    long double dy=t*z*(p-z)-y[id];
    return (z-x[id])*(z-x[id])+dy*dy;
}
long double get_derivative(int id,long double t,long double z)
{
    return z-x[id]+t*(p-2*z)*(t*z*(p-z)-y[id]);
}
long double get_second_derivative(int id,long double t,long double z)
{
    return 1+2*t*y[id]+t*t*(p*p-6*p*z+6*z*z);
}
long double get_distance(int id,long double t)
{
    vector<long double> root={0,p};
    if(t*p<1)
    {
        long double l=0,r=p;
        for(int i=0;i<M;i++)
        {
            long double mid=(l+r)/2;
            if(get_derivative(id,t,mid)<0)l=mid;
            else r=mid;
        }
        root.push_back((l+r)/2);
    }
    else
    {
        long double a=-1.5L*p;
        long double b=p*p/2+y[id]/t+1/(2*t*t);
        long double c=-p*y[id]/(2*t)-x[id]/(2*t*t);
        long double q=b-a*a/3;
        long double r=2*a*a*a/27-a*b/3+c;
        long double d=r*r/4+q*q*q/27;
        long double scale=max((long double)1,max(fabsl(r*r/4),fabsl(q*q*q/27)));
        if(d>=-1e-24L*scale)
        {
            d=max((long double)0,d);
            long double z=cbrtl(-r/2+sqrtl(d))+cbrtl(-r/2-sqrtl(d))-a/3;
            root.push_back(z);
        }
        else
        {
            long double angle=acosl(max((long double)-1,min((long double)1,-r/2/sqrtl(-q*q*q/27))));
            long double len=2*sqrtl(-q/3);
            for(int i=0;i<3;i++)root.push_back(len*cosl((angle+2*PI*i)/3)-a/3);
        }
        for(long double &z:root)
        {
            if(z<0||z>p)continue;
            for(int i=0;i<8;i++)
            {
                long double dd=get_second_derivative(id,t,z);
                if(fabsl(dd)<1e-24L)break;
                long double nz=z-get_derivative(id,t,z)/dd;
                if(nz<0||nz>p)break;
                z=nz;
            }
        }
    }
    long double ans=1e100L;
    for(long double z:root)
    {
        if(z>=0&&z<=p)ans=min(ans,get_value(id,t,z));
    }
    return ans;
}
bool hit(int id,long double t,long double radius)
{
    long double eps=1e-14L*max((long double)1,radius*radius);
    return get_distance(id,t)+eps<radius*radius;
}
bool valid(long double t,long double radius)
{
    if(t<0||t*p*p/4+radius>=h)return 0;
    for(int i=0;i<n;i++)
    {
        if(hit(i,t,radius))return 0;
    }
    return 1;
}
bool check(long double radius)
{
    if(radius>=h)return 0;
    long double max_t=4*(h-radius)/(p*p);
    vector<long double> cand={0,max_t};
    for(int i=0;i<n;i++)
    {
        long double center=y[i]/(x[i]*(p-x[i]));
        long double top=min(center,max_t);
        if(!hit(i,top,radius))continue;
        if(hit(i,0,radius))cand.push_back(0);
        else
        {
            long double l=0,r=top;
            for(int j=0;j<K;j++)
            {
                long double mid=(l+r)/2;
                if(hit(i,mid,radius))r=mid;
                else l=mid;
            }
            cand.push_back(l);
            cand.push_back(r);
        }
        if(center<=max_t)
        {
            if(hit(i,max_t,radius))cand.push_back(max_t);
            else
            {
                long double l=center,r=max_t;
                for(int j=0;j<K;j++)
                {
                    long double mid=(l+r)/2;
                    if(hit(i,mid,radius))l=mid;
                    else r=mid;
                }
                cand.push_back(l);
                cand.push_back(r);
            }
        }
        else cand.push_back(max_t);
    }
    sort(cand.begin(),cand.end());
    vector<long double> val;
    for(long double t:cand)
    {
        if(val.empty()||t-val.back()>1e-20L)val.push_back(t);
    }
    for(long double t:val)
    {
        if(valid(t,radius))return 1;
    }
    int siz=val.size();
    for(int i=1;i<siz;i++)
    {
        if(valid((val[i-1]+val[i])/2,radius))return 1;
    }
    return 0;
}
int main()
{
    ios::sync_with_stdio(false);
    cin.tie(nullptr);
    int t;
    cin>>t;
    for(int k=1;k<=t;k++)
    {
        cin>>n>>p>>h;
        for(int i=0;i<n;i++)cin>>x[i]>>y[i];
        long double l=0,r=h;
        for(int i=0;i<K;i++)
        {
            long double mid=(l+r)/2;
            if(check(mid))l=mid;
            else r=mid;
        }
        cout<<"Case #"<<k<<": "<<fixed<<setprecision(10)<<l<<'\n';
    }
    return 0;
}