LG15646 [ICPC 2022 Tehran R] Windcatchers
给定平面 \(n\) 个点（\(3\le n\le 4000\)），找一条直线作为中线，使中线上至少有两个点，且所有点到直线的距离 \(\ge\) 宽度（点在边界上允许）。最大化宽度，输出该宽度。坐标为整数且保证所有点不在一条直线上。
给一种 \(O(n^2\log n)\) 的做法（复杂度瓶颈为排序）。
一个显然的 \(O(n^3)\) 暴力是，枚举直线上的两个点，然后求其他点到这条至直线的距离。现在我们考虑固定直线经过一个点 \((x_0,y_0)\)，看看经过这个点的直线所能做到的最大宽度。
设当前的直线斜率为 \(k\)（暂时不考虑斜率不存在，即直线垂直 \(x\) 轴的情况），则直线方程为 \(y=k(x-x_0)+y_0\)。那么点 \((x_i,y_i)\) 到直线的距离可以表示为（直线距离公式）
这条直线对应的最大宽度为
可以发现绝对值内的形式可以看作一条自变量为 \(k\) 的折线，这条折线在 \(k=\frac{y_i-y_0}{x_i-x_0}\) 处的值为 \(0\)，右侧斜率为 \(|x_i-x_0|\)，左侧斜率为 \(-|x_i-x_0|\)。由于我们不考虑 \(k=\frac{y_i-y_0}{x_i-x_0}\) 时的值，于是拆分为两条左右射线（不含端点）即可。 问题就是求 \(k\) 时这些射线的纵坐标最小值。
不同的直线斜率 \(k\) 只有 \((n-1)\) 种，自然可以使用李超线段树维护这至多 \((n-1)\) 个离散坐标的纵坐标最小值，但是因为常数问题加上很容易写成 \(O(n^2\log^2n)\)，甚至都跑不到数据最大点。
如果我们只考虑斜率为正的射线（即折线的右半部分），可以变成如下的斜率优化问题
给 \(n\) 条线段和 \(n\) 个位置 \(x_i\)（\(x\) 递增），第 \(i\) 条线段斜率为 \(k_i\)，在 \(x_i\) 处纵坐标为 \(0\)，且对应的范围为 \([x_{i+1},x_n]\) ，求每个位置 \(x_i\) 的线段纵坐标最小值。
这个问题可以转化为单调栈维护上凸包，这里简述
- 维护一个从栈顶到栈底斜率递减的上凸包，意思是斜率越小的直线在 \(x\) 越大的位置越会成为最小值。
- 对于当前位置 \(x_i\)，如果栈顶对应线段的纵坐标值大于栈第二条线段的纵坐标值，则弹出栈顶，直到栈顶的元素纵坐标值最小。那么这个最小值就是该位置的答案。
- 现在要加入范围为 \([x_{i+1},x_n]\) 的线段，如果栈顶的线段的斜率大于该线段，则弹出，直到栈顶线段斜率小于该线段。
- 维持上凸包的正确形态（这条不是很好理解）：如果栈顶线段与该线段交点横坐标大于栈第二条线段与该线段横坐标，则弹出栈顶线段，直到不满足。可以自己画图理解。
对于斜率为负的射线，可以使用同样的方法。然后综合起来取最小值即可。
至此我们就可以 \(O(n)\) 地求这些位置各线段的最小值。
然后对于斜率不存在的直线，他们垂直于 \(x\) 轴，只需要一开始特判一下即可。
回到原题目，我们的算法就是
- 枚举直线上的一个端点 \((x_0,y_0)\)。
- 枚举直线上的另一个端点 \((x_i,y_i)\)，求出直线斜率 \(k=\frac{y_i-y_0}{x_i-x_0}\)。
- 对直线斜率排序并去重。
- 对于每个端点，找到他向左或向右线段所管理的区间，并将线段加入。
- 求每个斜率位置 \(k\) 这些线段的纵坐标最小值。这个最小值 \(\times \frac{1}{\sqrt{1+k^2}}\) 即为斜率为 \(k\) 的直线的宽度。
- 综合宽度的最大值。
可以发现复杂度瓶颈在对斜率排序。
实际写的过程中，使用 double 类型存数据我的写法会出现答案错误问题，故我不得不使用分数来记录数据，这就导致了常数可能会比较大（其实我换分数本地速度提升不明显）。可能极限数据点就差一点就能跑过去，于是不得不采用计时终止的方法，配合一开始乱序枚举端点才能通过此题。具体实现见代码。
#include <bits/stdc++.h>using namespace std;typedef long long ll;
typedef array<ll,2>ttfa;
const int N=16384;
const ll INF=0x7fffffff;int n,m;struct frac{ll x,y;frac(ll xx=0,ll yy=1){ll g=__gcd(abs(xx),abs(yy));x=abs(xx)/g,y=abs(yy)/g;if((xx<0)^(yy<0))x=-x;}bool operator < (const frac &B)const{return x*B.y<y*B.x;}bool operator > (const frac &B)const{return x*B.y>y*B.x;}bool operator == (const frac &B)const{return x*B.y==y*B.x;}
}pos[N];array<ll,2>a[N];
ll kk[N],bb[N];
int tot,num=0;map<ll,int>bot;inline bool cmp(int idx,int idy,int loc){if(kk[idx]==kk[idy])return bb[idx]<bb[idy];return pos[loc].x*(kk[idx]-kk[idy])<pos[loc].y*(bb[idy]-bb[idx]);
}
int bel[N];
int lef[N],rit[N];
int stk[N];//单调栈inline double calc(ll y,ll x){ll g=__gcd(abs(y),abs(x));y/=g,x/=g;return 1.0*y/x;
}pair<frac,int>lis[N];int main(){//freopen("P15646_17.in","r",stdin);//freopen("data3.out","w",stdout);scanf("%d",&n);for(int i=1;i<=n;++i){scanf("%lld%lld",&a[i][0],&a[i][1]);}random_shuffle(a+1,a+1+n);double ans=0;for(int i=1;i<=n;++i){pos[++m]=frac(a[i][0],1);++bot[a[i][0]];}sort(pos+1,pos+1+m);m=unique(pos+1,pos+1+m)-pos-1;for(int j=1;j<=m;++j){if(bot[pos[j].x]>=2){ll l=INF,r=INF;if(j>1)l=pos[j].x-pos[j-1].x;if(j<m)r=pos[j+1].x-pos[j].x;ans=max(ans,(double)min(l,r));}}for(int i=1;i<=n;++i){num=m=0;bb[tot=0]=INF;for(int j=1;j<=n;++j){if(i==j)continue;if(a[i][0]!=a[j][0]){lis[++num]=make_pair(frac(a[j][1]-a[i][1],a[j][0]-a[i][0]),j);}else{bb[0]=min(bb[0],abs(a[i][1]-a[j][1]));}}sort(lis+1,lis+1+num);m=num;for(int j=1;j<=m;++j)pos[j]=lis[j].first;m=unique(pos+1,pos+1+m)-pos-1;//pp();for(int j=1;j<=m;++j)bel[j]=0,lef[j]=rit[j]=0;for(int t=1,loc=1;t<=num;++t){int j=lis[t].second;frac k=frac(a[j][1]-a[i][1],a[j][0]-a[i][0]);ll dx=abs(a[j][0]-a[i][0]),dy=abs(a[i][1]-a[j][1]);while(pos[loc]<lis[t].first)++loc;if(k.x<0){if(!lef[loc]||(-kk[lef[loc]]>dx)){lef[loc]=++tot;kk[tot]=-dx,bb[tot]=-dy;}if(!rit[loc]||(kk[rit[loc]]>dx)){rit[loc]=++tot;kk[tot]=+dx,bb[tot]=+dy;}}else{if(!lef[loc]||(-kk[lef[loc]]>dx)){lef[loc]=++tot;kk[tot]=-dx,bb[tot]=+dy;}if(!rit[loc]||(kk[rit[loc]]>dx)){rit[loc]=++tot;kk[tot]=dx,bb[tot]=-dy;}}}int top=0;for(int j=1;j<=m;++j){while(top>=2&&!cmp(stk[top],stk[top-1],j))--top;if(top&&cmp(stk[top],bel[j],j))bel[j]=stk[top];if(rit[j]){while(top&&kk[stk[top]]>=kk[rit[j]])--top;while(top>=2){int r1=stk[top-1],r2=stk[top],r3=rit[j];if((bb[r2]-bb[r3])*(kk[r2]-kk[r1])>=(bb[r1]-bb[r2])*(kk[r3]-kk[r2]))--top;else break;}stk[++top]=rit[j];} }top=0;for(int j=m;j>=1;--j){while(top>=2&&!cmp(stk[top],stk[top-1],j))--top;if(top&&cmp(stk[top],bel[j],j))bel[j]=stk[top];if(lef[j]){while(top&&-kk[stk[top]]>=-kk[lef[j]])--top;while(top>=2){int r1=stk[top-1],r2=stk[top],r3=lef[j];if((bb[r3]-bb[r2])*(kk[r1]-kk[r2])<=(bb[r2]-bb[r1])*(kk[r2]-kk[r3]))--top;else break;}stk[++top]=lef[j];} }for(int k=1;k<=m;++k){int id=bel[k];double tmp=(double)(pos[k].x*kk[id]+pos[k].y*bb[id]);double kval=1.0*pos[k].x/pos[k].y;double tmp2=pos[k].y*sqrt(1.0+kval*kval);tmp/=tmp2;ans=max(ans,tmp);}if((double)clock()/CLOCKS_PER_SEC>=5.50)break;}printf("%.12lf\n",ans);return 0;
}
想做法花了十分钟，写过这题花了一整天（不止）。如果大佬有复杂度更小的做法或者优化建议也请多多指教啊！