IOI2020集训队作业-4 (CF607E, AGC038E, AGC030C)

A - CF607E Cross Sum

题意

有$n$条直线。

设集合$D$为：这些直线中两两交点构成的可重集合，其中一个点的出现次数等于有多少对直线的交点是这个点。

给出一个点$o$和一个整数$m$，你需要求$D$中到$o$最近的$m$个点到$o$的距离的和。

n≤50000,m≤3×min⁡{107,∣D∣}n\le 50000, m\le 3\times \min \{ 10^7,|D|\}n≤50000,m≤3×min{107,∣D∣}。时限7s。

Sol

二分求出第$m$近的点到$o$的距离。检查的时候需要统计在圆内的直线交点个数，可以转化成求直线与圆交形成的弦的交点个数。假设弦的端点与圆心的连线与$x$轴的夹角分别为$l_i$和$r_i$，则两条弦有交当且仅当$l_i<l_j<r_i<r_j$，这是个二维偏序的形式。这一步的复杂度是$O(n\log n\log V)$的，$V$与坐标范围和精度有关。

接下来用与上面类似的思想扫出所有的在圆内的交点，复杂度$O(m\log n)$。

实现细节

1）由于在求$l_i<l_j<r_i<r_j$的时候是角度的比较，值域非常小（是$(-\pi ,\pi ]$），所以这一步的比较中一定要把$eps$设得非常非常小（或者干脆不设$eps$）。
2）为了卡常，在求二维偏序的时候要让每个交点只被数一次（也就是数无序对$(i,j)$）。
3）注意考虑$\pi$和$-\pi$实际上是同一个角，要特殊判断一下，保证计数不重不漏。
4）算上在圆上的点，得到的点数可能会大于$m$。所以最后算答案的时候，先算出严格在圆内的点的点数$tot$以及它们到圆心的距离的和，然后再加上$(m-tot)\cdot d$就是答案（其中$d$是圆的半径）。

Code

#include <cstdio>
#include <iostream>
#include <algorithm>
#include <cstring>
#include <cmath>
#include <vector>
#define PB push_back
#define ll long long
#define db long double
using namespace std;
template <class T>
inline void rd(T &x) {
	x=0; char c=getchar(); int f=1;
	while(!isdigit(c)) { if(c=='-') f=-1; c=getchar(); }
	while(isdigit(c)) x=x*10-'0'+c,c=getchar(); x*=f;
}
const int maxn=5e4+10;
const db eps=1e-16,Pi=acos(-1.0);
int dcmp(db x) { return x<-eps?-1:(x>eps); }
struct Point {
	db x,y;
	Point(db x=0,db y=0): x(x),y(y) {}
	friend Point operator +(Point A,Point B) { return Point(A.x+B.x,A.y+B.y); }
	friend Point operator -(Point A,Point B) { return Point(A.x-B.x,A.y-B.y); }
	friend Point operator *(Point A,db B) { return Point(A.x*B,A.y*B); }
};
db Cross(Point A,Point B) { return A.x*B.y-A.y*B.x; }
db len(Point A) { return sqrt(A.x*A.x+A.y*A.y); }
db get_ag(db x) {
	if(x>Pi) x-=2*Pi;
	if(x<-Pi) x+=2*Pi;
	if(dcmp(x-Pi)==0||dcmp(x+Pi)==0) return Pi;
	return x;
}
db get_ag(Point A) { return get_ag(atan2(A.y,A.x)); }
Point o;
db ki[maxn],bi[maxn]; int n;
Point per[maxn];

db xi[maxn*20]; int m;
int find(db x) {
	int l=1,r=m;
	while(l<r) {
		int mid=l+r>>1;
		if(dcmp(xi[mid]-x)==0) return mid;
		if(xi[mid]>x) r=mid-1;
		else l=mid+1;
	}
	return l;
}

namespace TREE {
	int c[maxn*10];
	void init() { for(int i=1;i<=m;++i) c[i]=0; }
	void add(int i) { for(;i<=m;i+=i&-i) c[i]++; }
	int query(int i) { int ans=0; for(;i;i-=i&-i) ans+=c[i]; return ans; }
	int query(int l,int r) { return query(r)-query(l-1); }
}
#include <iomanip>
db L[maxn],R[maxn];
int Li[maxn],Ri[maxn];
int tot;
int id[maxn];
int rnk[maxn];
bool cmprnk(int x,int y) { return Li[x]!=Li[y]?Li[x]<Li[y]:Ri[x]>Ri[y]; }
ll sol(db d) {
	m=tot=0;
	for(int i=1;i<=n;++i) {
		Point p=per[i];
		db x=len(p-o);
		if(dcmp(x-d)>=0) continue;
		db l,r;
		if(dcmp(x)==0) {
			l=get_ag(atan(ki[i]));
			r=get_ag(l+Pi);
		}
		else {
			db ag=get_ag(p-o),delt=acos(x/d);
			l=get_ag(ag-delt),r=get_ag(ag+delt);
		}
		if(l>r) swap(l,r);
		L[++tot]=l,R[tot]=r,id[tot]=i;
		rnk[tot]=tot;
		xi[++m]=l,xi[++m]=r;
	}
	{
		sort(xi+1,xi+m+1);
		int p=0;
		for(int l=1,r;l<=m;l=r+1) {
			r=l;
			while(r+1<=m&&dcmp(xi[r+1]-xi[l])==0) r++;
			xi[++p]=xi[l];
		}
		m=p;
	}
	for(int i=1;i<=tot;++i) Li[i]=find(L[i]),Ri[i]=find(R[i]);
	sort(rnk+1,rnk+tot+1,cmprnk);
	ll ans=0; TREE::init();
	for(int u=1;u<=tot;++u) {
		int i=rnk[u];
		if(Li[i]+1<=Ri[i]-1)
			ans+=TREE::query(Li[i]+1,Ri[i]-1);
		TREE::add(Ri[i]);
	}
	return ans;	
}
vector<int> qans;
int ql,qr;
namespace BST {
	const int M=maxn;
	int fa[M],ch[M][2],val[M],id[M];
	int rt,ncnt;
	int get(int x) { return ch[fa[x]][1]==x; }
	void rotate(int x) {
		int f=fa[x],ff=fa[f],d=get(x);
		fa[x]=ff; if(ff) ch[ff][ch[ff][1]==f]=x;
		fa[ch[x][d^1]]=f; ch[f][d]=ch[x][d^1];
		fa[f]=x,ch[x][d^1]=f;
	}
	void splay(int x,int gl) {
		for(int f=fa[x];f!=gl;rotate(x),f=fa[x])
			if(fa[f]!=gl) rotate(get(x)==get(f)?f:x);
		if(!gl) rt=x;
	}
	void insert(int v,int i) {
		int u=rt,f=0;
		while(u) f=u,u=ch[u][v>val[u]];
		val[u=++ncnt]=v,id[u]=i;
		fa[u]=f; if(f) ch[f][v>val[f]]=u; else rt=u;
		splay(u,0);
	}
	void query(int x) {
		if(!x) return;
		if(val[x]<ql) return query(ch[x][1]);
		if(val[x]>qr) return query(ch[x][0]);
		qans.PB(id[x]);
		query(ch[x][0]),query(ch[x][1]);
	}
}

void getans(db d,int lim) {
	m=tot=0;
	for(int i=1;i<=n;++i) {
		Point p=per[i];
		db x=len(p-o);
		if(dcmp(x-d)>=0) continue;
		db l,r;
		if(dcmp(x)==0) {
			l=get_ag(atan(ki[i]));
			r=get_ag(l+Pi);
		}
		else {
			db ag=get_ag(p-o),delt=acos(x/d);
			l=get_ag(ag-delt),r=get_ag(ag+delt);
		}
		if(l>r) swap(l,r);
		L[++tot]=l,R[tot]=r,id[tot]=i;
		rnk[tot]=tot;
		xi[++m]=l,xi[++m]=r;
	}
	{
		sort(xi+1,xi+m+1);
		int p=0;
		for(int l=1,r;l<=m;l=r+1) {
			r=l;
			while(r+1<=m&&dcmp(xi[r+1]-xi[l])==0) r++;
			xi[++p]=xi[l];
		}
		m=p;
	}
	for(int i=1;i<=tot;++i) Li[i]=find(L[i]),Ri[i]=find(R[i]);
	sort(rnk+1,rnk+tot+1,cmprnk);
	db ans=0;
	for(int u=1;u<=tot;++u) {
		int i=rnk[u];
		if(Li[i]+1<=Ri[i]-1) {
			ql=Li[i]+1,qr=Ri[i]-1,qans.clear();
			BST::query(BST::rt);
			for(int j=0;j<qans.size();++j) {
				int x=id[qans[j]],y=id[i];
				Point p;
				p.x=(bi[y]-bi[x])/(ki[x]-ki[y]);
				p.y=p.x*ki[x]+bi[x];
				ans+=len(p-o),lim--;
			}
		}
		BST::insert(Ri[i],i);
	}
	ans+=lim*d;
	printf("%.10Lf",ans);
}	
int main() {
	rd(n),rd(o.x),rd(o.y); o.x/=1000,o.y/=1000;
	int lim; rd(lim);
	for(int i=1;i<=n;++i) rd(ki[i]),rd(bi[i]),ki[i]/=1000,bi[i]/=1000;
	for(int i=1;i<=n;++i)
		if(ki[i]==0) per[i]=Point(o.x,ki[i]*o.x+bi[i]);
		else {
			db k2=-1/ki[i];
			db b2=o.y-o.x*k2;
			per[i].x=(b2-bi[i])/(ki[i]-k2);
			per[i].y=ki[i]*per[i].x+bi[i];
		}
//	db lb=eps,rb=1e10;
//	int tcnt=0;
//	while(rb-lb>1e-7&&(++tcnt)<=200) {
//		db mid=(lb+rb)/2;
//		if(sol(mid)>=lim) rb=mid;
//		else lb=mid;
//	}
	db lb=0;
	for(db delt=1e12;delt>1e-9;delt/=2)
		if(sol(lb+delt)<lim) lb+=delt;

	getans(lb,lim);

//	for(int i=1;i<=n;++i) if(dcmp(len(per[i]-o)-lb)<0) id.PB(i);
//	db ans=0;
//	for(int i=0;i<id.size();++i)
//		for(int j=i+1;j<id.size();++j) {
//			if(dcmp(ki[id[i]]-ki[id[j]])==0) continue;
//			Point p;
//			p.x=(bi[id[j]]-bi[id[i]])/(ki[id[i]]-ki[id[j]]);
//			p.y=p.x*ki[id[i]]+bi[id[i]];
//			db tmp=len(p-o);
//			if(dcmp(tmp-lb)<0) lim--,ans+=tmp;
////			cout<<id[i]<<' '<<id[j]<<':'<<tmp<<endl;
//			if(lim==0) break;
//		}
//	ans+=lim*lb;
//	printf("%.10Lf\n",ans);
	return 0;
}

---

B - AGC038E Gachapon

题意

有一个随机数生成器，每次独立生成一个$[0,n)$的数，生成第$i$个数的概率是$p_i$。

给一个长度为$n$的数组$b_i$，问期望至少需要运行这个随机数生成器多少次，才能使对于每个$i$都满足第$i$个数的出现次数大于等于$b_i$。

$p_i$的给出方式是，给出$n$个整数$a_0,a_1,\cdots a_{n-1}$，则$p_i=\frac{a_i}{{\sum }_{j=0}^{n-1}{a}_j}$。

所有的运算在模$998244353$的意义下进行。

$n\le 400,1\le a_i,b_i$，且$\sum a_i\le 400,\sum b_i\le 400$

Sol

首先对题目要求的东西进行min-max容斥。假如枚举了一个某一个{1,2,3⋯n}\{1,2,3\cdots n\}{1,2,3⋯n}的子集$S$，我们现在要求期望至少进行多少轮，$S$中有一个数的出现次数达到了$b_i$。

我们对$S$中的元素记一个{xj}\{ x_j \}{xj​}，表示某个时刻每个数的出现次数。则期望操作的轮数等于（所有满足$x_j<b_j$的{xj}\{ x_j \}{xj​}出现的概率 * 这个情况期望会持续多少轮操作）。设$P$为操作到$S$集合中的元素的概率，那么对于任意一个情况它的期望持续轮数都是$\frac{1}{P}$。所以我们可以先不管每种情况的期望持续轮数，算出每种情况出现的概率的和之后乘上$\frac{1}{P}$就可以了。

现在考虑怎么求某种情况出现的概率：此时应恰好有$X={\sum }_{i\in S}{x}_i$次操作我们操作了$S$中的元素，并且每个元素被操作的次数也是确定了的。所以概率就应该是：

$$(X!)\cdot {\Pi }_{i\in S}(\frac{a_i}{{\sum }_{j\in S}{a}_j}{)}^{x_i}\frac{1}{x_i!}$$

这里不需要考虑$S$之外的元素是因为：1）$S$之外的元素被生成对{xi}\{ x_i \}{xi​}没有影响。2）由于$a_i\ge 1$，所以总是存在一个时刻，$S$中的元素被选的次数是$X$。

接下来考虑用$dp$优化上面的过程。观察发现，我们只需要知道$S$中所有元素的$a_i$之和与$x_i$之和就可以算$X!,(\frac{1}{\sum a_j}{)}^X$和$\frac{1}{P}$，而对于其它的部分以及容斥系数，每个元素的贡献都是独立的。所以把$\sum a_i$和$\sum b_i$作为$dp$状态，其它在加入一个元素的时候乘上就可以了。

时间复杂度$O((\sum a_i)(\sum b_i{)}^2)$。

Code

#include <cstdio>
#include <iostream>
#include <algorithm>
#include <cstring>
#define ll long long
using namespace std;
template <class T>
inline void rd(T &x) {
	x=0; char c=getchar(); int f=1;
	while(!isdigit(c)) { if(c=='-') f=-1; c=getchar(); }
	while(isdigit(c)) x=x*10-'0'+c,c=getchar(); x*=f;
}
const int mod=998244353;
const int N=410;
void Add(int &x,int y) { x+=y; if(x>=mod) x-=mod; }
int Pow(int x,int y) {
	int res=1;
	while(y) {
		if(y&1) res=res*(ll)x%mod;
		x=x*(ll)x%mod,y>>=1;
	}
	return res;
}
int f[2][N][N];
int a[N],b[N],n,s;
int inv[N],fac[N];
void getfac(int n) {
	fac[0]=1; for(int i=1;i<=n;++i) fac[i]=fac[i-1]*(ll)i%mod;
	inv[n]=Pow(fac[n],mod-2); for(int i=n;i>=1;--i) inv[i-1]=inv[i]*(ll)i%mod;
}
int main() {
	getfac(400);
	rd(n);
	for(int i=1;i<=n;++i) rd(a[i]),rd(b[i]),s+=a[i];
	int cur=0; f[cur][0][0]=mod-1;
	for(int i=1;i<=n;++i) {
		int nxt=cur^1; memcpy(f[nxt],f[cur],sizeof(f[cur]));
		for(int j=0;j<=400;++j)
			for(int k=0;k<=400;++k) if(f[cur][j][k]) {
				for(int x=0,t=1;x<b[i];++x,t=t*(ll)a[i]%mod)
					Add(f[nxt][j+a[i]][k+x],f[cur][j][k]*(ll)inv[x]%mod*t%mod*(ll)(mod-1)%mod);
			}
		cur=nxt;
	}
	int ans=0;
	for(int i=1;i<=400;++i) {
		int ii=Pow(i,mod-2);
		for(int j=0,t=ii;j<=400;++j,t=t*(ll)ii%mod) if(f[cur][i][j])
			Add(ans,f[cur][i][j]*(ll)t%mod*fac[j]%mod);
	}
	ans=ans*(ll)s%mod;
	printf("%d",ans);
	return 0;
}

---

C - AGC030C Coloring Torus

题意

称对于一个$n\times n$的方格（格子的左上角是$(0,0)$，右下角是$(n-1,n-1)$）和颜色数$K$的涂色方案是好的，当且仅当：

• 每个方格都涂了$K$种颜色中的一个。
• 每种颜色都至少被涂了一次。
• 对于任意两种颜色$i,j$，满足对所有颜色为$i$的格子，与它相邻的格子中颜色为$j$的格子的数量都是一样的。这里我们定义两个格子$(a,b)$和$(c,d)$相邻，当且仅当$a=c$且$b-d=\pm 1\mathrm{m}\mathrm{o}\mathrm{d}n$或者$b=d$且$a-c=\pm 1\mathrm{m}\mathrm{o}\mathrm{d}n$。

现在给你颜色数$K$，你可以在$[1,500]$中选择一个$n$，构造一种好的涂色方案。

$K\le 1000$

Sol

当$K$时$4$的倍数的时候，一种合法的构造是：

我们发现如果对于某几个$i$，把$2i-1$和$2i$这两种颜色改成同一种颜色，构造方案仍然是合法的。这样就可以把这种构造扩展到$K$不是$4$的倍数的情况。

Code

#include <cstdio>
#include <iostream>
#include <algorithm>
#include <cstring>
#define ll long long
using namespace std;
template <class T>
inline void rd(T &x) {
	x=0; char c=getchar(); int f=1;
	while(!isdigit(c)) { if(c=='-') f=-1; c=getchar(); }
	while(isdigit(c)) x=x*10-'0'+c,c=getchar(); x*=f;
}
const int N=510;
int m;
int a[N][N];
int mp[N*2];
int main() {
	rd(m);
	if(m==1) {
		printf("1\n1\n");
		return 0;
	}
	int n=(m+3)/4*2;
	int tot=n*2;
	int tmp=(4-m%4)%4;
	for(int i=0;i<tot;++i) mp[i]=i;
	int p=2*(n-tmp);
	for(int i=1;i<=tmp;++i) mp[2*(n-i)+1]=mp[2*(n-i)]=p++;
	for(int i=0;i<tot;i+=2) {
		int y=i/2;
		for(int j=0;j<n;++j,y=(y+1)%n)
			a[j][y]=mp[i+(j&1)];
	}
	printf("%d\n",n);
	for(int i=0;i<n;++i,puts(""))
		for(int j=0;j<n;++j)
			printf("%d ",a[i][j]+1);
	return 0;
}
	printf("1\n1\n");
		return 0;
	}
	int n=(m+3)/4*2;
	int tot=n*2;
	int tmp=(4-m%4)%4;
	for(int i=0;i<tot;++i) mp[i]=i;
	int p=2*(n-tmp);
	for(int i=1;i<=tmp;++i) mp[2*(n-i)+1]=mp[2*(n-i)]=p++;
	for(int i=0;i<tot;i+=2) {
		int y=i/2;
		for(int j=0;j<n;++j,y=(y+1)%n)
			a[j][y]=mp[i+(j&1)];
	}
	printf("%d\n",n);
	for(int i=0;i<n;++i,puts(""))
		for(int j=0;j<n;++j)
			printf("%d ",a[i][j]+1);
	return 0;
}