loj 6587「ICPC World Finals 2019」交通堵塞

loj

后面令\(p_i=r_i+g_i\)

可以发现经过\(P=\mathrm{lcm}(1,2\dots100)\)时间,红绿灯状态会恢复到0时刻的状态,并且由于红绿灯出现的时间区间的端点都是整数,所以可以原问题可改为随机一个整数开始时刻\(k\in[0,P)\),求在第\(i\)个位置碰到红灯的概率

由于在第\(i\)个位置碰到红灯的概率等价于在前\(i-1\)个位置是绿灯且在第\(i\)个位置是红灯,所以可以改为求在前\(i-1\)个位置是绿灯概率减去在前\(i\)个位置是绿灯概率,这里记为\(f_i\).求\(f_i\)的朴素想法是列出\(i\)个同余方程,第\(j\)个形如\(k+x_j\ge r_j \bmod p_j\),然后算出解集大小,除掉\(\mathrm{lcm}_{j=1}^{i}p_j\)就是\(f_i\)

现在不妨考虑一些特殊情况下的同余方程组

• 如果两个同余方程\(x\in S_1 \bmod p,x\in S_2 \bmod q\),满足\(\gcd(p,q)=1\),那么有\(|\{x|x\in [0,pq),x\in S_1 \bmod p,x\in S_2 \bmod q\}|=|S_1||S_2|\)
• 如果两个同余方程\(x\in S_1 \bmod p,x\in S_2 \bmod q\),满足\(p|q\),那么可以以\(O(q)\)的复杂度合并两个方程,得到\(x\in S \bmod q\)

如果所有的\(p_i\)能够满足只有最多一种质因数,那么就可以对于每种质因数\(p\)维护一个\(\mod p^{\lfloor\log_p (\max p_i)\rfloor}\)的同余方程,然后把每个方程的解集大小相乘就得到要求的东西.这里再设一个\(X\),然后枚举\(a\in[0,x)\),求所有开始时刻为\(kX+a(k\in Z)\)的贡献.注意到固定\(X,a\)之后,\(kX+a+x_j \ge r_j \bmod p_j\)可以改写成一个\(k\in S \bmod \frac{p_j}{\gcd(p_j,X)}\)的方程.如果取\(X=2520\),那么会有\(\forall q\in[1,100],\frac{q}{\gcd(q,X)}\)最多只有一种质因数,这样就可以转化成上述特殊情况去做了.朴素实现可以做到\(O(nX(\max p_j))\)

#include<bits/stdc++.h>
#define LL long long
#define db double

using namespace std;
const int N=500+10,M=100+10;
int rd()
{
    int x=0,w=1;char ch=0;
    while(ch<'0'||ch>'9'){if(ch=='-') w=-1;ch=getchar();}
    while(ch>='0'&&ch<='9'){x=x*10+(ch^48);ch=getchar();}
    return x*w;
}
LL md=3623169640452065472ll,yyb=2520;
LL gcd(LL a,LL b){return b?gcd(b,a%b):a;}
LL n,a[N][3],bt[M][M],sz[M];
db sm[N];
int sp[M],tp,pm[M];

int main()
{
	freopen("1.in","r",stdin);
	freopen("1.out","w",stdout);
	pm[1]=1,sp[++tp]=1;
	for(int i=2;i<=100;++i)
	{
		if(pm[i]) continue;
		int sx=i;
		while(sx*i<=100) sx*=i;
		sp[++tp]=sx;
		for(int j=i;j<=100;j+=i) pm[j]=sx;
	}
	n=rd();
	for(int i=1;i<=n;++i)
		for(int j=0;j<=2;++j)
			a[i][j]=rd();
	for(int h=0;h<yyb;++h)
	{
		for(int i=1;i<=100;++i)
		{
			if(pm[i]!=i) continue;
			sz[i]=i;
			for(int j=0;j<i;++j)
				bt[i][j]=1;
		}
		for(int i=1;i<=n;++i)
		{
			int np=a[i][1]+a[i][2],g=gcd(yyb,np),q=np/g;
			for(int j=0,k=0;j<pm[q];++j,k=(k+1)%q)
				if(bt[pm[q]][j]&&(k*yyb+h+a[i][0])%np<a[i][1]) --sz[pm[q]],bt[pm[q]][j]=0;
			db na=1.0/2520;
			for(int j=1;j<=tp;++j)
				na*=(db)sz[sp[j]]/(db)sp[j];
			sm[i]+=na;
		}
	}
	sm[0]=1;
	for(int i=0;i<=n;++i) printf("%.10lf\n",sm[i]-sm[i+1]);
	return 0;
}