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

loj

后面令(p_i=r_i+g_i)

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

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

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

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

如果所有的(p_i)能够满足只有最多一种质因数,那么就可以对于每种质因数(p)维护一个(mod p^{lfloorlog_p (max p_i) floor})的同余方程,然后把每个方程的解集大小相乘就得到要求的东西.这里再设一个(X),然后枚举(ain[0,x)),求所有开始时刻为(kX+a(kin Z))的贡献.注意到固定(X,a)之后,(kX+a+x_j ge r_j mod p_j)可以改写成一个(kin S mod frac{p_j}{gcd(p_j,X)})的方程.如果取(X=2520),那么会有(forall qin[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
",sm[i]-sm[i+1]);
	return 0;
}