kmjp's blog

競技プログラミング参加記です

AtCoder ARC #137 : F - Overlaps

ARC

うーむ、処理が重い。
https://atcoder.jp/contests/arc137/tasks/arc137_f

問題

以下の処理をN回繰り返す。
[0,1]の中から一様にランダムに2つの実数(x,y)を選ぶ。
その後、座標min(x,y)からmax(x,y)の間にシールを張るとする。

シールがK+1枚以上重なる箇所がない確率を求めよ。

解法

この問題においては、2N個の実数の相対的な大小だけが関係し、詳細な値は関係しない。
2N個の実数を昇順に並べ、2つずつN組のペアを作る場合、その組み方は(2N-1)!!通りである。
このうち(K+1)個以上のペアが共通部分を持つことがない組み合わせを数え上げることを考えよう。

2N個の実数のうち、区間の右端となる箇所をN個定めたとする。
そうすると各右端が対応する左端を何通り取ることができるかが一意に定まる。
この候補数を、左側から順に並べた数列Xとする。

Xが問題の条件を満たすには

  • 1≦X[i]≦K
  • X[i+1]≧X[i]-1
  • X[N]=1

でなければならない。
真ん中の観点に着目し、真ん中に違反する回数の偶奇をもとに包除原理を適用する。

1以上K以下の降下列で、要素毎に2以上減少するような数列Aを数え上げよう。

  • f(L,R,mask) := 多項式で、n次の係数は以下を満たすn要素の数列の個数に相当する
    • L以上R未満の降下列で、要素毎に2以上減少する
    • maskは、数列中にLを含むかどうか、および(R-1)を含むかどうかの2*2通り

まず分割統治法の要領で、f(1,K+1,*)を求める。
次に以下の多項式を考える。これらはf(1,K+1,*)から求められる。偶数次数の符号反転は、包除原理のための処理。

  • P(x) := 多項式で、n次の係数は以下を満たすn要素の数列の個数に相当する
    • 1以上K以下の降下列で、要素毎に2以上減少する。末尾の要素は1
    • ただし、偶数次数の項は符号反転しておく
  • Q(x) := 多項式で、n次の係数は以下を満たすn要素の数列の個数に相当する
    • 1以上K以下の降下列で、要素毎に2以上減少する。末尾の要素の縛りはない
    • ただし、偶数次数の項は符号反転しておく

求めるのは、Q(x)に対応する数列がいくつか続いたあと最後P(x)に対応する数列が続くような長さNの数列なので、P* (1+Q+Q^2+Q^3....)のN次の係数を求めればよい。
(1+Q+Q^2+Q^3....)のところはバイナリ法でもよいし等比数列の和の要領で解いても良い。
下記は前者で求めている。

const int mo=998244353;
ll modpow(ll a, ll n = mo-2) {
	ll r=1; a%=mo;
	while(n) r=r*((n%2)?a:1)%mo,a=a*a%mo,n>>=1;
	return r;
}


template<class T> vector<T> fft(vector<T> v, bool rev=false) {
	int n=v.size(),i,j,m;
	for(int m=n; m>=2; m/=2) {
		T wn=modpow(5,(mo-1)/m);
		if(rev) wn=modpow(wn);
		for(i=0;i<n;i+=m) {
			T w=1;
			for(int j1=i,j2=i+m/2;j2<i+m;j1++,j2++) {
				T t1=v[j1],t2=v[j2];
				v[j1]=t1+t2;
				v[j2]=ll(t1+mo-t2)*w%mo;
				while(v[j1]>=mo) v[j1]-=mo;
				w=(ll)w*wn%mo;
			}
		}
	}
	for(i=0,j=1;j<n-1;j++) {
		for(int k=n>>1;k>(i^=k);k>>=1);
		if(i>j) swap(v[i],v[j]);
	}
	if(rev) {
		ll rv = modpow(n);
		FOR(i,n) v[i]=(ll)v[i]*rv%mo;
	}
	return v;
}

template<class T> vector<T> MultPoly(vector<T> P,vector<T> Q,bool resize=false) {
	if(resize) {
		int maxind=0,pi=0,qi=0,i;
		int s=2;
		FOR(i,P.size()) if(norm(P[i])) pi=i;
		FOR(i,Q.size()) if(norm(Q[i])) qi=i;
		maxind=pi+qi+1;
		while(s*2<maxind) s*=2;
		P.resize(s*2);Q.resize(s*2);
		if(s<=64) { //fastpath
			vector<T> R(s*2);
			for(int x=0;x<2*s;x++) for(int y=0;x+y<2*s;y++) (R[x+y]+=P[x]*Q[y])%=mo;
			return R;
		}
	}
	P=fft(P), Q=fft(Q);
	for(int i=0;i<P.size();i++) P[i]=(ll)P[i]*Q[i]%mo;
	return fft(P,true);
}


array<vector<ll>,4> C;
int N,K;

array<vector<ll>,4> hoge(int L,int R) {
	array<vector<ll>,4> D;
	if(L+1==R) {
		D[0]={1LL}; // なし
		D[1]={0LL}; // Lのみ
		D[2]={0LL}; // Rのみ
		D[3]={0LL,(ll)L}; // LとR-1
		return D;
	}
	
	int M=(L+R)/2;
	auto A=hoge(L,M);
	auto B=hoge(M,R);
	int x,y,i;
	FOR(x,4) FOR(y,4) {
		if(x/2&&y%2) continue;
		int t=(x%2)+(y/2)*2;
		vector<ll> E=MultPoly(A[x],B[y],true);
		if(D[t].size()<E.size()) D[t].resize(E.size());
		FOR(i,E.size()) {
			(D[t][i]+=E[i])%=mo;
		}
		
	}
	return D;
	
}

void solve() {
	int i,j,k,l,r,x,y; string s;
	
	
	cin>>N>>K;
	C=hoge(1,K+1);
	
	
	vector<ll> X(N+1),Y(N+1);
	FOR(x,4) {
		FOR(y,N+1) if(y<C[x].size()) {
			if(x%2) (X[y]+=C[x][y])%=mo;
			(Y[y]+=C[x][y])%=mo;
		}
	}
	
	Y[0]=0;
	FOR(x,N+1) if(x%2==0) {
		X[x]=mo-X[x];
		Y[x]=mo-Y[x];
	}
	vector<ll> Z={1};
	FOR(i,18) {
		Y[0]++;
		Z=MultPoly(Z,Y,true);
		Y[0]--;
		Z.resize(1<<18);
		for(x=N+1;x<1<<18;x++) Z[x]=0;
		
		Y=MultPoly(Y,Y,true);
		Y.resize(1<<18);
		for(x=N+1;x<1<<18;x++) Y[x]=0;
	}
	Z=MultPoly(X,Z,true);
	Z.resize(1<<18);
	
	ll ret=Z[N]%mo;
	for(i=2*N-1;i>=1;i-=2) ret=ret*modpow(i)%mo;
	cout<<ret<<endl;
	
	
}

まとめ

最初実行に10秒ぐらいかかってしまい、NTTのコードを微修正してギリギリ6秒を切った。
自分のNTTのライブラリ遅いのかな…。