感觉这种类似于让你找第 \(k\) 大距离的计算几何题其实都挺套路的。
二分一个答案 \(t\)，然后思考一下什么样的点对满足原点到它们的连线的距离 \(\le t\)。以原点为圆心 \(t\) 为半径画个圆，然后分情况讨论：
- 如果 \((x_i,y_i)\) 或 \((x_j,y_j)\) 在圆内，那么必然符合条件。
- 如果 \((x_i,y_i)\) 和 \((x_j,y_j)\) 都在圆外，那么考虑做这两个点到圆的切线，那么我们惊奇地发现，原点到这两点连线的距离 \(\le t\)，当且仅当这两点对应的两条切线的辐角组成的区间只有相交和包含关系。
于是就做完了，离散化一下然后树状数组做即可。注意由于包含和无交都合法，因此如果 \(l>r\) 直接交换 \(l,r\) 即可，不用拆成 \([l,2\pi]\) 和 \([0,r]\) 两个区间。
时间复杂度 \(n\log n\log V\)。
const int MAXN=1e5;
const double EPS=1e-9;
const double Pi=acos(-1);
int n;ll k;
int sgn(double x){return ((x<-EPS)?-1:((x<EPS)?0:1));}
struct point{
double x,y;
point(){x=y=0;}
point(double _x,double _y){x=_x;y=_y;}
friend point operator+(const point &X,const point &Y){return point(X.x+Y.x,X.y+Y.y);}
friend point operator-(const point &X,const point &Y){return point(X.x-Y.x,X.y-Y.y);}
friend point operator*(const point &X,const double &Y){return point(X.x*Y,X.y*Y);}
friend double operator |(const point &X,const point &Y){return X.x*Y.y-X.y*Y.x;}
friend double operator ^(const point &X,const point &Y){return X.x*Y.x+X.y*Y.y;}
friend bool operator ==(const point &X,const point &Y){return sgn(X.x-Y.x)==0&&sgn(X.y-Y.y)==0;}
friend bool operator !=(const point &X,const point &Y){return sgn(X.x-Y.x)!=0||sgn(X.y-Y.y)!=0;}
double operator ~()const{return sqrt(x*x+y*y);}
double operator !()const{return atan2(y,x);}
}p[MAXN+5];
void norm(double &d){
while(d<-EPS)d+=2*Pi;
while(d>2*Pi-EPS)d-=2*Pi;
}
int t[MAXN*5+5];
void add(int x,int v){for(int i=x;i<=MAXN*5;i+=(i&(-i)))t[i]+=v;}
int query(int x){int ret=0;for(int i=x;i;i&=(i-1))ret+=t[i];return ret;}
ll calc(double d){
ll res=0;
vector<pair<double,double> >vec;
vector<pii>nvec;vector<double>key;
int cnt=0;
for(int i=1;i<=n;i++){
if(~p[i]<(d+EPS));
else{
double ang=acos(d/(~p[i]));
double l=(!p[i])-ang,r=(!p[i])+ang;
norm(l);norm(r);
if(l>r)swap(l,r);
vec.pb(mp(l,r));
}
}
res=1ll*n*(n-1)/2;
for(pair<double,double>p:vec)key.pb(p.fi),key.pb(p.se);
sort(key.begin(),key.end());
for(pair<double,double>p:vec){
int L=lower_bound(key.begin(),key.end(),p.fi)-key.begin();
int R=lower_bound(key.begin(),key.end(),p.se)-key.begin();
nvec.pb(mp(L+1,R+1));
}
memset(t,0,sizeof(t));
sort(nvec.begin(),nvec.end());
for(int i=0;i<nvec.size();i++){
int l=nvec[i].fi,r=nvec[i].se;
res-=query(r)-query(l-1);add(r,1);
}
return res;
}
int main(){
scanf("%d%lld",&n,&k);
for(int i=1;i<=n;i++)scanf("%lf%lf",&p[i].x,&p[i].y);
// calc(0.8);
double l=EPS,r=20000,p=0;
while(fabs(r-l)>EPS){
double mid=(l+r)/2;
if(calc(mid)>=k)p=r=mid;
else l=mid;
}printf("%.10lf\n",p);
return 0;
}