Grid Points (Jinan Regional 2022. F)
考虑二分斜率 \(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;
}
/*
*/
