考慮二分斜率 \(k\)，用斜率為 \(k\) 的直線切簡單多邊形，要計算在直線之下，在多邊形區域內的整點數量之和。
把每條邊向 \(x\) 軸做垂線，差分成求若干個直角梯形的整點數，對每個梯形用類歐計算即可。
你以為這就完了？那就想的太簡單了。
在簡單多邊形的邊界上，每條邊上的點、每個角上的點的貢獻是極難處理的，而且這題切完後會變成多個簡單多邊形，並不能簡單判斷應該取正還是取負貢獻。下面的做法來自 std 的實現方式（
考慮一種“畢克定理”式的貢獻計算方式：假設在每個點放了一個半徑可以看作 \(0\) 的圓盤，對於一個簡單多邊形，我們計算在多邊形內的圓盤的重量之和，其中嚴格在內部的點貢獻 \(1\)，邊上的點貢獻 \(\frac 12\)，角上的點貢獻 \(\frac {\alpha}{2\pi}\)。
這個“重量”的貢獻差分到每個梯形上，是可以簡單加減的。
我們看看在這種計算方式下會算出什麼：假設取的斜率 \(k\) 是極遠處的一個座標互質的點，這條“新線”切多邊形的任何線都不產生整點。
- 內部整點的貢獻是正確的。
- 邊上的點的貢獻為 \(\frac 12\)，需要列舉每條原始多邊形上的邊加上貢獻。
- 角上的點的貢獻之和為 \(\frac{n}{2}-1\)（多邊形內角和），其中 \(n\) 為這個多邊形的點數。而我們需要的點數是 \(n-2\)（新的線切出來的兩個點一定不是整點），把 \(\frac{n}{2}-1\) 乘上 \(2\) 就對了。
我們透過把貢獻拆到每個角度和邊上，就完美解決了邊界的問題！
最後是一些二分的細節，我們可以取大於 \(V^2\) 的一個大質數 \(P\)，取出 \((P,i),(i,P)\) 點列作為斜率進行二分。
最後對於斜率相等的點還要二分是第幾個，需要二分 \(x\) 並且實現數座標 \(\le x\) 的所有點，這個在數點的函式里加一個限制就行。
時間複雜度 \(O(n\log^2 V)\)。
// what is matter? never mind.
//#pragma GCC optimize("Ofast")
//#pragma GCC optimize("unroll-loops")
//#pragma GCC target("sse,sse2,sse3,sse4,popcnt,abm,mmx,avx,avx2")
#include<bits/stdc++.h>
#define For(i,a,b) for(int i=(a);i<=(b);++i)
#define Rep(i,a,b) for(int i=(a);i>=(b);--i)
#define ll long long
#define int long long
#define ull unsigned long long
#define i128 __int128
#define SZ(x) ((int)((x).size()))
#define ALL(x) (x).begin(),(x).end()
using namespace std;
#define fi first
#define se second
#define pb push_back
#define mkp make_pair
typedef pair<int,int>pii;
typedef vector<int>vi;
#define maxn 500005
#define inf 0x3f3f3f3f
// i=[0..n] (a*i+b)/c
i128 floor_sum(i128 n,i128 a,i128 b,i128 c){
i128 res=0;
// For(i,0,n)res+=(a*i+b)/c;return res;
res=(n+1)*(b/c)+n*(n+1)/2*(a/c);
a%=c,b%=c;
if(a<0)a+=c,res-=n*(n+1)/2;
if(b<0)b+=c,res-=n+1;
if(!a)return res+(n+1)*(b/c);
// i128 m=(a*n+b)/c;
// i128 t=floor_sum(m-1,c,c-b-1,a);
// res=n*m-t;
i128 t=a*(n+1)+b;
if(t>=c) res+=floor_sum(t/c-1,c,t%c,a);
return res;
}
i128 floor_sum(i128 a,i128 b,i128 c,i128 l,i128 r){
return floor_sum(r,a,b,c)-floor_sum(l-1,a,b,c);
}
#define double long double
typedef double db;
const db eps=1e-10,pi=3.14159265358979323846;
struct P{
int x,y;
P(int x=0,int y=0):x(x),y(y){}
P&operator +=(P o){return x+=o.x,y+=o.y,*this;}
P&operator -=(P o){return x-=o.x,y-=o.y,*this;}
P&operator *=(int o){return x*=o,y*=o,*this;}
P&operator /=(int o){return x/=o,y/=o,*this;}
friend P operator +(P a,P b){return a+=b;}
friend P operator -(P a,P b){return a-=b;}
friend P operator *(P a,int b){return a*=b;}
friend P operator /(P a,int b){return a/=b;}
friend bool operator <(P a,P b){return a.x==b.x?a.y<b.y:a.x<b.x;}
friend int operator %(P a,P b){return a.x*b.x+a.y*b.y;} // dot
friend int operator *(P a,P b){return a.x*b.y-a.y*b.x;} // cross
inline double ang(){return atan2(y,x);}
inline double l(){return sqrt((*this)*(*this));}
void read(){cin>>x>>y;}
void out(){cout<<"("<<x<<","<<y<<")"<<endl;}
};
const int W=1000000000,V=W*W+W+1;
pii getid(int x){
if(x<=V)return mkp(x,V);
return mkp(V,V*2-x);
}
int gcd(int x,int y){
x=abs(x),y=abs(y);
return __gcd(x,y);
}
int n,k;
P a[maxn];
db ang[maxn],ang2[maxn];
int tmp[maxn];
bool isin(P a,int k1,int k2,int xl){
return a.x<=xl && (i128)a.y*k2<=(i128)a.x*k1;
}
i128 Div(i128 a,i128 b){
if(a>=0) return a/b;
return (a-b+1)/b;
}
int gclamp(i128 x){
return (int)max(min(x,(i128)V),-(i128)V);
}
void clamp(int&l,int&r,i128 k,i128 b){
if(k==0){
if(b<0)l=1,r=0;
return;
}
if(k>0){
i128 t=Div(b,k);
r=min(r,gclamp(t));
}else{
i128 t=Div((-b)+(-k)-1,-k);
l=max(l,gclamp(t));
}
}
int seg(P u,P v,int k1,int k2,int xl){
if(u.x>v.x)swap(u,v);
int d=gcd(v.x-u.x,v.y-u.y);
int l=1,r=d-1;
i128 dx=(v.x-u.x)/d,dy=(v.y-u.y)/d;
clamp(l,r,dx,xl-u.x);
clamp(l,r,dy*k2-dx*k1,(i128)u.x*k1-(i128)u.y*k2);
return max(0ll,r-l+1);
}
pair<int,db>calc(P u,P v,int id,int k1,int k2,int xl){
int res1=0; db res2=0;
if(u.x==v.x)return mkp(res1,res2);
if(u.x>v.x)swap(u,v);
int l=u.x+1,r=min(v.x-1,xl);
int a=v.y-u.y,c=v.x-u.x,b=u.y*c-u.x*a-1;
clamp(l,r,(i128)k2*a-(i128)k1*c,-(i128)b*k2);
// cout<<"l,r "<<l<<" "<<r<<"\n";
if(l<=r) res1+=2*(floor_sum(a,b,c,l,r));//cout<<"a,b,c "<<a<<" "<<b<<" "<<c<<"\n";
// cout<<"res1 "<<res1<<"\n";
l=u.x+1,r=min(v.x-1,xl);
clamp(l,r,(i128)k1*c-(i128)k2*a,(i128)b*k2-1);
a=k1,b=0,c=k2;
// cout<<"l,r "<<l<<" "<<r<<"\n";
if(l<=r) res1+=2*(floor_sum(a,b,c,l,r));//cout<<"a,b,c "<<a<<" "<<b<<" "<<c<<"\n";
// cout<<"res1 "<<res1<<"\n";
// cout<<"res1 "<<res1<<"\n";
res1+=seg(P(u.x,0),u,k1,k2,xl);
res1+=seg(P(v.x,0),v,k1,k2,xl);
res1+=tmp[id];
// cout<<"res1 "<<res1<<"\n";
if(isin(u,k1,k2,xl)) res2+=ang[id];
if(isin(v,k1,k2,xl)) res2+=1-ang[id];
return mkp(res1,res2);
}
int counts(int k1,int k2,int xl){
// cout<<"count "<<k1<<" "<<k2<<" "<<xl<<"\n";
int res1=0; db res2=0;
For(i,0,n-1){
P u=a[i],v=a[(i+1)%n];
tmp[i]=seg(u,v,k1,k2,xl);
if(u.x==v.x)continue;
auto [c1,c2]=calc(u,v,i,k1,k2,xl);
if(u.x<v.x) res1-=c1,res2-=c2;
else res1+=c1,res2+=c2;
}
For(i,0,n-1){
if(isin(a[i],k1,k2,xl)) res2+=ang2[i];
// res1+=seg(a[i],a[(i+1)%n],k1,k2,xl);
res1+=tmp[i];
}
int res=res1+round(res2);
// cout<<"Res "<<res<<"\n";
res/=2;
return res;
}
void work()
{
n=read(),k=read();
For(i,0,n-1)a[i].read();
For(i,0,n-1){
P u=a[i]-a[(i+n-1)%n],v=a[i]-a[(i+1)%n];
ang2[i]=atan2(u*v,u%v);
if(ang2[i]<0) ang2[i]+=2*pi;
ang2[i]/=pi;
u=a[i],v=a[(i+1)%n];
if(u.x>v.x)swap(u,v);
if(u.x==v.x)ang[i]=0;
else ang[i]=((v-u).ang()+pi/2)/pi;
}
int l=0,r=2*V;
while(l+1!=r){
int mid=l+r>>1;
pii kk=getid(mid);
if(counts(kk.fi,kk.se,W)>=k)r=mid;
else l=mid;
}
pii k1=getid(l),k2=getid(r);
int rk=k-counts(k1.fi,k1.se,W);
l=0,r=W;
while(l+1!=r){
int mid=l+r>>1;
if(counts(k2.fi,k2.se,mid)-counts(k1.fi,k1.se,mid)>=rk)r=mid;
else l=mid;
}
int ry=(((i128)r*k2.fi)/k2.se);
cout<<r<<" "<<ry<<"\n";
}
signed main()
{
int T=read();
while(T--)work();
return 0;
}
/*
*/