计算几何（三）

LOJ #2401. 「JOISC 2017 Day 4」Dragon 2

首先这题有个条件就是每次询问的$A,B$不会和之前的询问相同，那么我们对于一个询问$A,B$，我们找出龙数更少的那个集合假设是$A$，那么我们只要能在$O(C|A|)$的时间复杂度内解决这个问题，我们的总复杂度就是$O((n+q)\sqrt{n}C)$的。
因为对于$|A|\le \sqrt{n}$，一次询问复杂度就是$O(\sqrt{n}C)$的。
对于$|A|>\sqrt{n}$，只会有$O(\frac{|B|}{|A|})$次询问复杂度涉及到$A$，
也就是说$|A|$对于总复杂度的贡献是$O(nc)$的，有因为这样的$A$只有$O(\sqrt{n})$个，
所以总复杂度就是$O((n+q)\sqrt{n}C)$的。

那么我们可以发现，我们把每个点关于我们的线段$(D_1,E_1,D_2,E_2)$，求出两个极角，把每个点看做二元组$(x,y)$，那么:开始盗图

看完这两张图应该就可以马上想到怎么用极角转化为二维偏序问题然后可持久化线段树解决了。
$\mathcal{A}CCode$
1A真是太感动了。

#include<bits/stdc++.h>
#define maxn 30005
#define rep(i,j,k) for(int i=(j),LIM=(k);i<=LIM;i++)
#define per(i,j,k) for(int i=(j),LIM=(k);i>=LIM;i--)
#define db long double
#define pdd pair<db,db>
#define pb push_back
#define mp make_pair
#define Pi 3.1415926535897932384626433832795
using namespace std;

char cb[1<<16],*cs=cb,*ct=cb;
#define getc() (cs==ct&&(ct=(cs=cb)+fread(cb,1,1<<16,stdin),cs==ct)?0:*cs++)
template<class T>void read(T &res){
	char ch;bool f=0;
	for(;!isdigit(ch=getc());) if(ch == '-') f = 1;
	for(res=ch-'0';isdigit(ch=getc());res=res*10+ch-'0');
	(f) && (res=-res); 
}
int n,m;
vector<pdd>G[maxn][2];
vector<db>sb[maxn][2][2];
vector<int>rt[maxn][2];
db RD(db a){ 
	for(;a >= 2 * Pi; a-= 2 *Pi);
	for(;a < 0; a += 2 * Pi);
	return a;
}
#define maxp maxn * 100
int lc[maxp],rc[maxp],sm[maxp],tot;
int newnode(int u){
	lc[++tot] = lc[u] , rc[tot] = rc[u] , sm[tot] = sm[u];
	return tot;
}
void ins(int &u,int l,int r,int p){
	u = newnode(u);
	sm[u] ++;
	if(l == r) return;
	int m = l+r>>1;
	p <= m ? ins(lc[u],l,m,p) : ins(rc[u],m+1,r,p);
}
int qry(int u,int l,int r,int ql,int qr){
	if(l>qr||ql>r||ql>qr||!u) return 0;
	if(ql<=l&&r<=qr) return sm[u];
	int m = l+r>>1;
	return qry(lc[u],l,m,ql,qr) + qry(rc[u],m+1,r,ql,qr); 
}

int main(){
	read(n),read(m);
	static db A[maxn],B[maxn],D1,E1,D2,E2,A1,A2;
	static int C[maxn];
	rep(i,1,n) read(A[i]),read(B[i]),read(C[i]);
	read(D1),read(E1),read(D2),read(E2);
	A1 = atan2(E2-E1,D2-D1) , A2 = atan2(E1-E2,D1-D2);
	rep(i,1,n){
		db a = atan2(B[i]-E1,A[i]-D1) , b = atan2(B[i]-E2,A[i]-D2);
		if(RD(a - A1) < Pi)
			G[C[i]][0].pb(mp(RD(a-A1),RD(A2-b)));
		else 
			G[C[i]][1].pb(mp(RD(A1-a),RD(b-A2)));
	}
	rep(i,1,m) rep(t,0,1){
		vector<db>&v0 = sb[i][t][0] , &v1 = sb[i][t][1];
		for(auto v:G[i][t]) v0.pb(v.first),v1.pb(v.second);
		sort(v0.begin(),v0.end());
		sort(v1.begin(),v1.end());
		vector<vector<int> >ad;
		ad.resize(v0.size()+1);
		for(auto v:G[i][t])
			ad[lower_bound(v0.begin(),v0.end(),v.first)-v0.begin()+1].pb(lower_bound(v1.begin(),v1.end(),v.second)-v1.begin()); 
		rt[i][t].resize(v0.size()+1);
		rep(j,1,v0.size()){
			rt[i][t][j]=rt[i][t][j-1];
			for(int v:ad[j])
				ins(rt[i][t][j],0,v1.size()-1,v);
		}
	}
	int q;read(q);
	for(;q--;){
		int a,b;bool fg = 0;
		read(a),read(b);
		if(G[a][0].size() + G[a][1].size() > G[b][0].size() + G[b][1].size()){
			swap(a,b);
			fg = 1;
		}
		int ans = 0;
		rep(t,0,1){
			vector<db>&v0 = sb[b][t^1][0] , &v1 = sb[b][t^1][1];
			for(auto u:G[a][t])
				ans += qry(rt[b][t^1][upper_bound(v0.begin(),v0.end(),RD(Pi-u.first))-v0.begin()] , 0 , v1.size()-1 , 0 , upper_bound(v1.begin(),v1.end(),RD(Pi-u.second))-v1.begin()-1);
		}
		if(fg){
			rep(t,0,1){
				vector<db>&v0 = sb[b][t][0] , &v1 = sb[b][t][1];
				for(auto u:G[a][t])
					ans += 
						qry(rt[b][t][rt[b][t].size()-1],0,v1.size()-1,upper_bound(v1.begin(),v1.end(),RD(u.second))-v1.begin(),v1.size()-1) - 
						qry(rt[b][t][upper_bound(v0.begin(),v0.end(),RD(u.first))-v0.begin()] , 0 , v1.size()-1 , upper_bound(v1.begin(),v1.end(),RD(u.second))-v1.begin(),v1.size()-1);
			}
		}	
		else{
			rep(t,0,1){
				vector<db>&v0 = sb[b][t][0] , &v1 = sb[b][t][1];
				for(auto u:G[a][t])
					ans += qry(rt[b][t][upper_bound(v0.begin(),v0.end(),u.first)-v0.begin()] , 0 , v1.size()-1 , 0 , upper_bound(v1.begin(),v1.end(),u.second)-v1.begin()-1);
			}
		}
		printf("%d\n",ans);
	}
}

JSOI 2018 战争

闵可夫斯基和+判断点是否在凸包内。

$\mathcal{A}CCode$

#include<bits/stdc++.h>
#define maxn 200005
#define rep(i,j,k) for(int i=(j),LIM=(k);i<=LIM;i++)
#define per(i,j,k) for(int i=(j),LIM=(k);i>=LIM;i--)
#define LL long long
#define Ct const
#define pb push_back
using namespace std;

int n,m,q;

struct Pt{
	LL x,y;
	Pt (Ct LL &x=0,Ct LL &y=0):x(x),y(y){}
	Pt operator +(Ct Pt &B)Ct{ return Pt(x+B.x,y+B.y); }
	Pt operator -(Ct Pt &B)Ct{ return Pt(x-B.x,y-B.y); }
	bool operator <(Ct Pt &B)Ct{ return x == B.x ? y < B.y : x < B.x; }
	LL operator *(Ct Pt &B)Ct{ return x * B.y - y * B.x; }
};
vector<Pt>P,Q;
vector<Pt> Convex(Ct vector<Pt>&A){
	vector<Pt>B = A , R;
	sort(B.begin(),B.end());
	for(auto u:B){
		for(;R.size() >= 2 && (u-R[R.size()-2]) * (R.back()-R[R.size()-2]) <= 0;R.pop_back());
		R.pb(u);
	}
	int L = R.size();
	per(i,B.size()-1,0){
		auto u = B[i];
		for(;R.size() > L && (u-R[R.size()-2]) * (R.back()-R[R.size()-2]) <= 0;R.pop_back());
		R.pb(u);
	}
	return R; // the last point is same as the first point
}
vector<Pt>Minkowski(Ct vector<Pt>&A,Ct vector<Pt>&B){
	vector<Pt>s1,s2;
	rep(i,1,A.size()-1) s1.pb(A[i] - A[i-1]);
	rep(i,1,B.size()-1) s2.pb(B[i] - B[i-1]);
	vector<Pt>R;
	R.pb(A[0] + B[0]);
	for(int i=0,j=0;i<s1.size() || j<s2.size();){
		static Pt P;
		if(j >= s2.size() || (i < s1.size() && s1[i] * s2[j] <= 0))P = R.back() + s1[i++];
		else P = R.back() + s2[j++];
		for(;R.size() >= 2 && (P-R.back()) * (R.back() - R[R.size()-2]) == 0;R.pop_back());
		R.pb(P);
	}
	if(R.size() > 3 && (R.back() - R[R.size() - 2]) * (R[1] - R[0]) == 0){
		rep(i,0,R.size()-2) 
			R[i] = R[i+1];
		R.resize(R.size()-2);
		R.pb(R[0]);
	}
	return R;
}
bool cmp(Ct Pt &A,Ct Pt &B){ return A * B < 0; }
LL sqr(LL a){ return a * a; }
int main(){
	scanf("%d%d%d",&n,&m,&q);
	P.resize(n) , Q.resize(m);
	rep(i,0,n-1) scanf("%lld%lld",&P[i].x,&P[i].y);
	rep(i,0,m-1) scanf("%lld%lld",&Q[i].x,&Q[i].y),Q[i].x*=-1,Q[i].y*=-1;
	P = Convex(P) , Q = Convex(Q);
	P = Minkowski(P,Q);
	rep(i,1,P.size()-2) P[i] = P[i] - P[0];
	for(;q--;){
		static Pt A;
		scanf("%lld%lld",&A.x,&A.y);
		A = A - P[0];
		if(A * P[1] < 0 || A * P[P.size()-2] > 0) puts("0");
		else{
			int t = lower_bound(P.begin()+1,P.end()-1,A,cmp) - P.begin();
			if(A * P[t] == 0){
				if(sqr(A.x) + sqr(A.y) <= sqr(P[t].x) + sqr(P[t].y)) puts("1");
				else puts("0");
			}
			else{
				if((A - P[t-1]) * (P[t] - P[t-1]) >= 0) puts("1");
				else puts("0"); 
			}
		}
	}

}

CF1195F Geometers Anonymous Club

闵可夫斯基和的性质。
$\mathcal{A}CCode$

#include<bits/stdc++.h>
#define maxn 100005
#define rep(i,j,k) for(int i=(j),LIM=(k);i<=LIM;i++)
#define per(i,j,k) for(int i=(j),LIM=(k);i>=LIM;i--)
using namespace std;

int n,q,A[maxn*3],B[maxn*3],ans[maxn];
int gcd(int a,int b){ return !b ? a : gcd(b,a%b); }
vector<pair<int,int> >G[maxn];
vector<int>E[maxn];

pair<int,int>SL(pair<int,int>a){
	int t = gcd(abs(a.first) , abs(a.second));
	a.first /= t , a.second /= t;
	return a;
}

int tr[maxn];
void upd(int u,int v){ for(;u<=n;u+=u&-u) tr[u] +=v ;}
int qry(int u){ int r=0;for(;u;u-=u&-u) r += tr[u];return r; }
map<pair<int,int>,int>mAp;

int main(){
	scanf("%d",&n);
	rep(i,1,n){
		int K;scanf("%d",&K);
		rep(j,0,K-1){
			scanf("%d%d",&A[j],&B[j]);
			if(j) G[i].push_back(SL(make_pair(A[j]-A[j-1],B[j]-B[j-1])));
		}
		G[i].push_back(SL(make_pair(A[0] - A[K-1] , B[0] - B[K-1])));
	}
	scanf("%d",&q);
	rep(i,1,q){
		scanf("%d%d",&A[i],&B[i]);
		E[B[i]].push_back(i);
	}
	rep(i,1,n){
		for(auto u : G[i]){
			int t = mAp[u];
			upd(t+1,1) , upd(i+1,-1);
			mAp[u] = i;
		}
		for(int u : E[i])
			ans [u] = qry(A[u]);
	}
	rep(i,1,q) printf("%d\n",ans[i]);
}

CF958E3 Guard Duty (hard)

直接分治构造即可。
$\mathcal{A}CCode$

#include<bits/stdc++.h>
#define maxn 20005
#define rep(i,j,k) for(int i=(j),LIM=(k);i<=LIM;i++)
#define per(i,j,k) for(int i=(j),LIM=(k);i>=LIM;i--)
using namespace std;

int n;
int x[maxn],y[maxn],c[maxn],ans[maxn],LOC;
bool cmp(const int &u,const int &v){ return (x[u] - x[LOC]) * (y[v] - y[LOC]) - (y[u] - y[LOC]) * (x[v] - x[LOC]) < 0; }

void solve(int l,int r){
	if(l>r) return;
	int loc = l;
	rep(i,l,r) if(x[c[i]] < x[c[loc]]) loc = i;
	swap(c[loc] , c[l]);
	LOC = c[l];
	sort(c+l+1,c+r+1,cmp);
	int s0=0,s1=0,k;
	rep(i,l+1,r) if(s0 == s1 && ((c[i] <= n) ^ (c[l] <= n))){
		k = i;
		break;
	}
	else if(c[i] <= n) s0 ++;
	else s1++;
	ans[c[l]] = c[k] , ans[c[k]] = c[l];
	solve(l+1,k-1)  ,solve(k+1,r);
}

int main(){
	scanf("%d",&n);
	rep(i,1,2*n) scanf("%d%d",&x[i],&y[i]) , c[i] = i;
	solve(1,2*n);
	rep(i,1,n) printf("%d\n",ans[i] - n);
}

#2586. 「APIO2018」选圆圈

$KDT$优化暴力，复杂度还是$O(n^2)$的。
随手旋转个角度就没人能卡你了。
$\mathcal{A}CCode$

#include<bits/stdc++.h>
#define maxn 300005
#define eps 1e-3
using namespace std;

int n,c[maxn];
struct Point{ 
	double mn[2],mx[2],s[2];int r,ch[2],ans;
}a[maxn];
void upd(int u){
	a[u].mn[0]=a[u].s[0]-a[u].r;
	a[u].mn[1]=a[u].s[1]-a[u].r;
	a[u].mx[0]=a[u].s[0]+a[u].r;
	a[u].mx[1]=a[u].s[1]+a[u].r;
	for(int i=0;i<2;i++) if(a[u].ch[i]) for(int j=0;j<2;j++)
		a[u].mn[j]=min(a[u].mn[j],a[a[u].ch[i]].mn[j]),
		a[u].mx[j]=max(a[u].mx[j],a[a[u].ch[i]].mx[j]);
}
int D,rt;
bool cmp1(const int &u,const int &v){ return a[u].s[D]<a[v].s[D]; }
bool cmp2(const int &u,const int &v){ return a[u].r==a[v].r?u<v:a[u].r>a[v].r; }
void Build(int &u,int l,int r,int d){
	if(l>r) return;
	int mid=(l+r)>>1;D=d,nth_element(c+l,c+mid,c+r+1,cmp1);
	u=c[mid];
	Build(a[u].ch[0],l,mid-1,d^1),Build(a[u].ch[1],mid+1,r,d^1),upd(u);
}
double si , co;
double sqr(double a){ return a * a; }
void del(int u,Point p){
	if(!u || p.s[0]+p.r<a[u].mn[0] - eps || p.s[1]+p.r<a[u].mn[1] - eps || 
		p.s[0]-p.r>a[u].mx[0] + eps || p.s[1]-p.r>a[u].mx[1] + eps) return;
	if(!a[u].ans && sqr(a[u].s[0]-p.s[0])+sqr(a[u].s[1]-p.s[1])<=sqr(p.r+a[u].r)+eps)
		a[u].ans = p.ans , a[u].s[0]=a[u].s[1]=0,a[u].r=-2e9;
	del(a[u].ch[0],p),del(a[u].ch[1],p),upd(u);
}

int main(){//freopen("1","r",stdin);
	scanf("%d",&n);
	si = 0.1 , co = sqrt(1 - si * si);
	for(int i=1;i<=n;i++){
		scanf("%lf%lf%d",&a[i].s[0],&a[i].s[1],&a[i].r);
		double x = a[i].s[0] * co + a[i].s[1] * si , y = - a[i].s[0] * si + a[i].s[1] * co;
		a[i].s[0]=x,a[i].s[1]=y,c[i]=i;
	}
	Build(rt,1,n,0);
	sort(c+1,c+1+n,cmp2);
	for(int i=1;i<=n;i++) if(!a[c[i]].ans)
		a[c[i]].ans=c[i],del(rt,a[c[i]]);
	for(int i=1;i<=n;i++) printf("%d%c",a[i].ans," \n"[i==n]);
}

CF1019E Raining season

现在看来这题好像挺水的，至少比求树上直径的期望要水
边分治后建出两边的凸包再闵科夫斯基和即可。
$\mathcal{A}CCode$

#include<bits/stdc++.h>
#define maxn 400005
#define rep(i,j,k) for(int i=(j),LIM=(k);i<=LIM;i++)
#define per(i,j,k) for(int i=(j),LIM=(k);i>=LIM;i--)
#define pb push_back
#define vc vector
#define LL long long
#define Ct const
#define mp make_pair
#define db double
using namespace std;

char cb[1<<16],*cs=cb,*ct=cb;
#define getc() (cs==ct&&(ct=(cs=cb)+fread(cb,1,1<<16,stdin),cs==ct)?0:*cs++)
template<class T>void read(T &res){
	char ch;
	for(;!isdigit(ch=getc()););
	for(res=ch-'0';isdigit(ch=getc());res=res*10+ch-'0');
}
int n,m;
struct Pt{
	LL x,y;
	db ag;
	Pt (Ct LL &x=0,Ct LL &y=0):x(x),y(y){}
	Pt operator +(Ct Pt &B)Ct{ return Pt(x+B.x,y+B.y); }
	Pt operator -(Ct Pt &B)Ct{ return Pt(x-B.x,y-B.y); }
	LL operator *(Ct Pt &B)Ct{ return x * B.y - y * B.x; }
	bool operator <(Ct Pt &B)Ct{ return y == B.y ? x < B.x : y > B.y; }
	void ang(){ ag = atan2(y , x); }
}cst[maxn<<1];
vc<Pt>H;
int info[maxn],Prev[maxn<<1],to[maxn<<1],vis[maxn<<1],sz[maxn],cnt_e=1,tot;
void Node(int u,int v,Pt C){ Prev[++cnt_e]=info[u],info[u]=cnt_e,to[cnt_e]=v,cst[cnt_e]=C;  }
vector<pair<int,Pt> >G[maxn];

void dfs0(int u,int ff){
	int p = 0;
	for(auto v:G[u]) if(v.first != ff){
		if(!p) Node(u,v.first,v.second) , Node(v.first,u,v.second) , p = u;
		else{
			++tot;
			Node(tot,p,Pt()),Node(p,tot,Pt());
			p = tot;
			Node(p,v.first,v.second),Node(v.first,p,v.second);
		} 
		dfs0(v.first,u);
	}
}

void dfs1(int u,int ff,int tsz,int &mn,int &rt){
	sz[u] = 1;
	for(int i=info[u],v;i;i=Prev[i]) if(v=to[i],!vis[i] && v!=ff){
		dfs1(v,u,tsz,mn,rt);
		if(max(sz[v] , tsz - sz[v]) < mn)
			mn =  max(sz[v] , tsz - sz[v]) , rt = i;
		sz[u] += sz[v];
	}
}

int Gert(int u,int tsz){
	int mn = 0x3f3f3f3f , rt = -1;
	dfs1(u,0,tsz,mn,rt);
	return rt;
}

void ser(int u,int ff,Pt C,vc<Pt>&P){	
	P.pb(C);sz[u] = 1;
	for(int i=info[u],v;i;i=Prev[i]) if((v=to[i]) != ff && !vis[i])
		ser(v,u,C+cst[i],P) , sz[u] += sz[v];
}

vc<Pt> Convex(vc<Pt>P){
	sort(P.begin(),P.end());
	vc<Pt>R;
	#define Rb R.back()
	#define Rs R[R.size()-2]
	for(Pt u:P){
		for(;R.size() >= 2 && (u-Rs) * (Rb-Rs) <= 0;R.pop_back());
		R.pb(u);
	}
	return R;
}

bool cmp(const Pt &u,const Pt &v){ return u.ag < v.ag; }
vc<Pt> Minkowski(Ct vc<Pt>&P,Ct vc<Pt>&Q){
	vc<Pt>R,sp,sq;
	rep(i,1,P.size()-1){
		sp.pb(P[i] - P[i-1]);
		if(sp.back().x || sp.back().y) sp.back().ang();
		else sp.pop_back();
	}
	rep(i,1,Q.size()-1){
		sq.pb(Q[i] - Q[i-1]);
		if(sq.back().x || sq.back().y) sq.back().ang();
		else sq.pop_back();
	}
	R.pb(P[0]+Q[0]);
	for(int i=0,j=0;i<sp.size() || j<sq.size();)
		if(j >= sq.size() || (i < sp.size() && sp[i].ag > sq[j].ag)) R.pb(R.back() + sp[i++]);
		else R.pb(R.back() + sq[j++]);
	return R;
}

void Solve(int u){
	if(u == -1) return;
	vc<Pt>P,Q;
	vis[u] = vis[u^1] = 1;
	ser(to[u],to[u^1],cst[u],P);
	ser(to[u^1],to[u],Pt(0,0),Q);
	P = Convex(P) , Q = Convex(Q);
	P = Minkowski(P,Q);
	for(Pt v:P) H.pb(v);
	Solve(Gert(to[u],sz[to[u]]));
	Solve(Gert(to[u^1],sz[to[u^1]]));
}

int main(){freopen("1.in","r",stdin);
	read(n),read(m);
	rep(i,1,n-1){
		int u,v;Pt C;
		read(u),read(v),read(C.x),read(C.y);
		G[u].pb(mp(v,C)) , G[v].pb(mp(u,C));
	}
	tot = n;
	dfs0(1,0);
	Solve(Gert(1,tot));
	H = Convex(H);
	
	int j = 0;
	rep(i,0,m-1){
		if(j < (int)H.size() - 1)
			printf("%lld %lld\n",H[j].y,H[j+1].y);
		for(;j < (int)H.size() - 1 && (H[j+1] - H[j]) * Pt(1,-i) <= 0;j++);
		printf("%lld%c",j < H.size() ? H[j].y + H[j].x * i : 0," \n"[i==m-1]);
	}
}