bzoj4520 CQOI2016 K远点对

题意：给定$N$个点，求欧几里得距离第$K$大的点对的欧几里得距离的平方。
数据范围：$1\le N\le 10^5$，$1\le K\le 100$，$K\le \frac{N\times (N-1)}{2}$，$0\le X_i,Y_i<2^{31}$，时限$3s$

做法：网上大多都是随机化，要么是不太能给出复杂度证明的k-d树。这里我写了一个确定性的做法，且有稳定的复杂度。
我们考虑这样做，每次跑出一个凸包，然后在这个凸包上跑最远点对，然后把最远点对的一个点删掉，继续这样跑，跑$K$次，我们发现第$K$远的点对一定有一个点在这$K$个点里了，那我们拿一个堆维护，$O(NKlog(K))$的扫一遍就可以得出答案了，我们的总复杂度是$O(NKlog(N))$，然后我们就得到了不知道几分的好成绩……
考虑优化这个东西，我们发现这个凸包每次只需要$O(N)$的时间复杂度，因为我们可以从$(-1,-1)$对这些点极角排序，然后跑两遍凸包算法一拼就求出凸包了。然后关键就在于这个求最远点对，这里有一个算法能把它从$Nlog(N)$优化到$O(N)$。
这个算法利用了四边形不等式的性质，四边形不等式的定义是如果一个矩阵，满足对于任意一个$2\times 2$的子矩阵，假设是
$$\left[\begin{array}{cc}a& b\\ c& d\end{array}\right]$$
都有$a+d\ge b+c$，就称这个矩阵满足四边形不等式。这个性质可以推出对于任意一个子矩阵，左上角加右下角一定大于等于左下角加右上角，证明就是将这个子矩阵包含的所有$2\times 2$的矩阵的$a+d-b-c$加起来，刚好就能得到这个结果。顺便提一句，四边形不等式其实更多在dp中起优化作用，类似的有决策单调性，在矩阵里有完全单调性与之对应，满足的是任意两行两列的交点若有左上小于右上，则有左下小于右下。这个有什么用呢？很明显，满足四边形不等式的矩阵肯定满足这个完全单调性。
我们发现这里凸包的距离矩阵似乎也满足这个性质，假设我们对于一个$2\times 2$的矩阵，左上角在$(i,j)$，满足的就是$dis(p_i,p_j)+dis(p_{i+1},p_{j+1})\ge dis(p_i,p_{j+1})+dis(p_j,p_{i+1})$，画个图就很直观了：

这里我们可以利用三角形两边之和大于第三边，有：
$AE+BE\le AB$，$CE+DE\le CD$
相加可得$AE+BE+CE+DE\le AB+CD$
刚好就是上面这个式子，但是这样就完了吗，直觉告诉你并没有，这里有一个反例就是$i=j$的时候，如图：

其实这个也很好解决，我们可以把这个矩阵的主对角线以下的部分平移到另一部分的右边，剩下的地方用$-\infty$填充，如图：

这个矩阵就满足四边形不等式了。
然后这里有一个神仙算法把这个矩阵每一行的最大值用$O(n+m)$的时间复杂度求出，因为我们有刚才的结论，我们发现每行的最大值位置是单调上升的，这个可以用我们推出的完全单调性证明，因为对于第$i$行的最大值位置$x$前面都满足$A_{i,j}\ge A_{i,x}$，然后我们就发现$A_{i+1}$肯定也有这些式子，所以易证这是单调的。
我们发现$n\times m$的矩阵明显不好，我们想个办法把$m$减到$n$，这里我们可以把矩阵变成这样：我们找出一些行使得$A_{i,i}\ge A_{i,i+1}$，这样我们就是一个$n\times n$的矩阵了。
这里举个例子，比如矩阵
$$\left[\begin{array}{ccccccc}1.& 2.& 5& 3& 4& 2& 1\\ 2& 4& 8& 9& 11& 10& 9\\ 3& 6& 11& 13& 16& 15& 14\end{array}\right]$$
这个矩阵好像具备刚才的性质，我们来消一下，首先我们选择第一行的第一个和第二个位置，发现$1<2$，好第一行明显没用了，删。
$$\left[\begin{array}{cccccc}2.& 5.& 3& 4& 2& 1\\ 4& 8& 9& 11& 10& 9\\ 6& 11& 13& 16& 15& 14\end{array}\right]$$
然后我们继续，发现$2<5$，删
$$\left[\begin{array}{ccccc}5.& 3.& 4& 2& 1\\ 8& 9& 11& 10& 9\\ 11& 13& 16& 15& 14\end{array}\right]$$
然后我们继续，发现$5>3$了，$5$这一排可能是有用的，留着，进第二排。
$$\left[\begin{array}{ccccc}5& 3& 4& 2& 1\\ 8& 9.& 11.& 10& 9\\ 11& 13& 16& 15& 14\end{array}\right]$$
这里$9<11$，这时候我们就把$9$这列删了，然后我们发现前一列可能不满足条件，于是退后一列
$$\left[\begin{array}{cccc}5.& 4.& 2& 1\\ 8& 11& 10& 9\\ 11& 16& 15& 14\end{array}\right]$$
发现$5>4$，这里我们对$5$放心了，可以继续往下做，注意这里如果$5$比较小的话下场是和之前$1,2$下场一样的。
下一步就是
$$\left[\begin{array}{cccc}5& 4& 2& 1\\ 8& 11.& 10.& 9\\ 11& 16& 15& 14\end{array}\right]$$
发现$11>10$，过
$$\left[\begin{array}{cccc}5& 4& 2& 1\\ 8& 11& 10& 9\\ 11& 16& 15.& 14.\end{array}\right]$$
这里$15>14$，大没有留$14$之必要，删掉，
$$\left[\begin{array}{ccc}5& 4& 2\\ 8& 11& 10\\ 11& 16& 15\end{array}\right]$$
注意如果$15$较小也需要后退一排，然后消出来可能不足$n$列，找出一些补齐即可。
这里我们说一下这个做法的证明，不难发现我们消到第$m$列的时候，有${\forall }_{1\le i\le m}{\forall }_{i\le j\le m}{A}_{i,j}\ge A_{i,j+1}$，证明也很简单，消掉一些行的四边形不等式因为我们推出的第一个结论，肯定还是满足四边形不等式，然后因为我们构造的矩阵满足${\forall }_{1\le i\le m}{A}_{i,i}\ge A_{i,i+1}$，再将四边形不等式的定义式拿出来变化一下就有$a-b\ge c-d$，我们就能得到前几行也肯定是满足$A_{x,i}\ge A_{x,i+1}$，而一行反过来的话（$A_{x,i}\le A_{x,i+1}$），类比的我们可以得到后几行也一定是这样，这样的话我们就可以遇到第$A_{m,m+1}\ge A_{m,m}$时就可以大胆删了。复杂度的话每次是每选中一行，比较的右边的那个位置就$+1$，每个选中时的列产生的贡献只会最多退回一次，所以就是$O(n+m)$的。
既然有上面这个就好办多了，每次我们对$n\times n$的，把奇数行拿出来一跑，最后一块扫一遍得到偶数的答案。
时间复杂度的递归式：$T(n)=T(\frac{n}{2})+O(n)$，即$O(n)$
然后这题我们就做完了！
时间复杂度：$O(Nlog(N)+NKlog(K))$，瓶颈在最后的堆。
代码又慢又臭又长：

#include<iostream>
#include<cstring>
#include<cassert>
#include<cmath>
#include<map>
#include<set>
#include<queue>
#include<stack>
#include<vector>
#include<cstdio>
#include<time.h>
#include<algorithm>
using namespace std;
#define REP(i,x,y) for(ll i=x;i<=y;i++)
#define rep(i,n) REP(i,1,n)
#define rep0(i,n) REP(i,0,n-1)
#define repG(i,x) for(ll i=pos[x];~i;i=e[i].next)
#define ll long long
#define db double
const ll N=1e6+7;
const ll INF=9223372036854775806;
ll n,k,sz;
struct pt{ll x,y;}p[N],g[N];
ll s[N],dw[N],z[N];
bool is[N],e[N];
ll Sq(ll x){return x*x;}
ll dis(pt a,pt b){return Sq(a.x-b.x)+Sq(a.y-b.y);}//算欧几里得距离
ll qdis(ll a,ll b){
 if(a==b)return 0;
 if(a>b||b>a+sz)return -INF;
 if(b>sz)b-=sz;
 return dis(g[a],g[b]);//矩阵查询
}
bool cmp2(ll x,ll y){return x>y;}
bool cmp(pt a,pt b){return (b.x+1)*(a.y+1)<(b.y+1)*(1+a.x);}
bool cross(pt a,pt b,pt c){return (b.y-a.y)*(c.x-b.x)<=(b.x-a.x)*(c.y-b.y);}
vector<ll>v1[50],v2[50],d[50],w[50];
void work(ll x){
 if(v1[x].size()<3){//递归边界
  rep0(i,v1[x].size()){
   ll mx=0,p;
   rep0(j,v2[x].size()){
    ll o=qdis(v1[x][i],v2[x][j]);
    if(o>mx){
     mx=o;
     p=j;
    }
   }
   d[x].push_back(p);
  }
  return;
 }
 w[x].clear();
 ll li=0,h=1;
 w[x].push_back(0);
 while(h<v2[x].size()){//将n*m矩阵转成n*n
  if(qdis(v1[x][li],v2[x][w[x][w[x].size()-1]])<qdis(v1[x][li],v2[x][h])){
   w[x].pop_back();
   if(li)li--;
   else{
    w[x].push_back(h);
    h++;
   }
  }
  else{
   if(li<v1[x].size()-1){
    w[x].push_back(h);
    li++;
   }
   h++;
  }
 }
 rep0(i,v2[x].size())e[i]=0;
 rep0(i,w[x].size())e[w[x][i]]=1;
 ll nw=v1[x].size()-w[x].size();
 w[x].clear();
 rep0(i,v2[x].size()){
  if(e[i])w[x].push_back(i);
  else if(nw)w[x].push_back(i),nw--;
 }
 v2[x+1].clear();
 v1[x+1].clear();
 rep0(i,w[x].size())v2[x+1].push_back(v2[x][w[x][i]]);
 rep0(i,v1[x].size())if(i&1)v1[x+1].push_back(v1[x][i]);//拿出奇数行
 work(x+1);//递归
 d[x].clear();
 ll tp=0;
 rep0(i,v1[x].size()){//确定答案
  if(i&1)tp=d[x+1][i/2],d[x].push_back(w[x][tp]);
  else{
   ll hp=(i==v1[x].size()-1)?v1[x].size()-1:d[x+1][i/2],mx=0,p;
   REP(j,tp,hp){
    ll o=qdis(v1[x][i],v2[x][w[x][j]]);
    if(o>mx)mx=o,p=j;
   }
   d[x].push_back(w[x][p]);
  }
 }
}
priority_queue<ll,vector<ll>,greater<ll> >Q;
int main(){
 scanf("%lld%lld",&n,&k);
 rep(i,n)scanf("%lld%lld",&p[i].x,&p[i].y);
 if(n<400){
  ll cnt=0;
  rep(i,n)rep(j,i-1)s[++cnt]=dis(p[j],p[i]);
  sort(s+1,s+cnt+1,cmp2);
  printf("%lld\n",s[k]);
  return 0;
 }
 sort(p+1,p+n+1,cmp);
 rep(u,k){
  ll tp=0,cnt=0;
  rep(i,n){
   if(is[i])continue;
   while(cnt>1&&cross(p[s[cnt-1]],p[s[cnt]],p[i]))cnt--;
   s[++cnt]=i;
  }
  rep(i,cnt)g[i]=p[s[i]],z[i]=s[i];
  for(ll i=s[cnt];i;i--){
   if(is[i])continue;
   while(tp>1&&cross(p[s[tp-1]],p[s[tp]],p[i]))tp--;
   s[++tp]=i;
  }
  REP(i,2,tp-1)g[++cnt]=p[s[i]],z[cnt]=s[i];
  sz=cnt;
  v1[1].clear();
  v2[1].clear();
  rep(i,sz)v1[1].push_back(i);
  rep(i,sz*2)v2[1].push_back(i);
  work(1);
  ll id,mx=0;
  rep(i,sz){
   ll o=qdis(i,v2[1][d[1][i-1]]);
   if(o>mx)mx=o,id=i;
  }
  is[z[id]]=1;
  dw[u]=z[id];
 }
 ll gg=0;
 memset(is,0,sizeof(is));
 rep(i,k){//堆计算答案
  rep(j,n){
   if(is[j])continue;
   if(gg<k)Q.push(dis(p[dw[i]],p[j])),gg++;
   else{
    ll o=dis(p[dw[i]],p[j]);
    if(o>Q.top()){
     Q.push(o);
     Q.pop();
    }
   }
  }
  is[dw[i]]=1;
 }
 printf("%lld\n",Q.top());
 return 0;
}

我才不会告诉你我是前几天集训的时候学的这个算法