kmjp's blog

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

TopCoder SRM 684 Div1 Medium DivFree、Div2 Medium DivFree2

ゴリ押し解法でも定数倍高速化で通せたのか…。
https://community.topcoder.com/stat?c=problem_statement&pm=14170

問題

整数N,Kに対し、以下の特性を持つ数列の組み合わせ数を答えよ。

  • 素数はNで、各要素は1~Kのいずれかの整数である。
  • 連続する2要素を前からA,Bとする。A≦BかA mod B !=0の最低片方を満たす。

解法

O(NK)解法は愚直なDPで解ける。まずF(n,k)を以下のように置く。
F(n,k) := n要素の数列であり、最後の要素はkで、かつ問題分中の隣接要素の条件を満たしているものの数

隣接要素の条件は、Bに対し条件を満たすAと満たさないAは満たすAの方が多いため、満たさないAの分を減らした方が速い。
 F(n,k) = \sum_i F(n-1,i) - \sum_{k>i,i|k} F(n-1,k)
上記FをDPで求め、F(N,*)の総和を取ればよい。
Div2はこれで解ける。

Div1はO(NK)では解けないので、何かしらの改善が必要。
自分が見た範囲では3通りの方法を見つけた。

同じ要素をまとめて高速化

O(NK)でもN,Kが1000位までは試せるので、一度F(n,k)を見てみよう。
すると、F(n,k)はkが近い範囲では割と同じ値を取ることに気づく。
F(n,k)をkの値ごとに分類すると凡そ2√K個程度のグループになる。(K=50000の最大ケースだと、グループは446個)

そこで、F(n,k)の計算をするときに同じグループに属する複数のkに対し計算を省略することでO(N√K)で解くことができる。
結構時間がギリギリなので、mod10^9+7を取る回数を減らすと何とか通る。
自分の解法は1.7s程度かかった。

同じ要素をまとめて高速化+行列累乗

F(n,*)をK要素のベクトルの数列の第n項と見なす。
するとF(n,*)はF(n-1,*)の値で求まるので行列累乗テクで高速にF(N,*)を求めることができる。
とはいえ元のK要素をそのまま行列にすると50000次正方行列の乗算が必要になってTLEする。
しかし上のテクで446次正方行列に落とし込むと、ギリギリ行列乗算+バイナリ法で間に合うようだ。
こちらはO(K√K*logN)になるのかな。

包除原理

他人のコードが何をやっているかわからなかったが、以下でわかった。
Topcoder Single Round Match 684 - Codeforces

隣接要素の条件がなければ、数列の取り方はK^Nである。
ここから隣接要素の条件に反するケースを引こう。

ここで重要な事実として、連続して隣接要素の条件に反するのは、高々logK要素までである。
隣接要素の条件に反するには、BがAの約数かつBがA未満でなければならない。
よって連続で条件に反する場合、A,A/2,A/4,A/8…と並べるのが最長で、それでもlog_2(K)≦20要素である。

上記Editorialに則り、以下のように変数を取る。

  • Z[i] := i要素の数列で、隣接要素の条件を満たすもの。
  • W[i] := i要素の数列で、隣接要素の条件を満たさないもの。
    • すなわちZ[i]+W[i] = K^iである。
  • CS[i] := 数列の(i+1)要素における隣接関係i個において、全て隣接要素の条件を満たさないケース。

ここで上記サイトにある以下の関係が登場する。
 W_i = W_{i-1}*K + \sum_j (-1)^{j+1} CS_j Z_{i-(j+1)}
最初の項は簡単で、すでに最初(i-1)要素の時点で条件を満たさないなら、最後の要素が何であってもやはり条件は満たさない。
後ろの項は、(i-j)要素目までは条件を満たすのに、その後追加したj要素がいずれも直前の要素に対し隣接要素における条件を連続して違反するケースを示している。
上記のとおり連続して違反するのは高々20要素なので、jは20までループすれば十分である。

  • 1の累乗を取っているのは、包除原理の要素で同じ構成を多重カウントしないためである。

以下、上は同じ要素をまとめたDP、下は包除原理解である。

ll dp[2][50505];
ll* from=dp[0];
ll* to=dp[1];

int mapto[50505];
int rev[50505];
int num[50505];
vector<pair<int,int>> del[505];
ll mo=1000000007;

class DivFree {
	public:
	int N,K;
	int dfcount(int n, int k) {
		N=n;K=k;
		int i,j,x;

		ZERO(dp);
		for(j=1;j<=k;j++) {
			to[j]=k;
			for(x=j*2;x<=k;x+=j) to[j]--;
		}
		int md=0;
		ZERO(num);
		mapto[0]=0;
		for(j=1;j<=k;j++) {
			mapto[j]=mapto[j-1]+(to[j]!=to[j-1]);
			rev[mapto[j]]=j;
			md=max(md,mapto[j]);
			num[mapto[j]]++;
		}
		
		for(i=2;i<=md;i++) {
			del[i].clear();
			int x=rev[i];
			for(j=x*2;j<=k;j+=x) {
				if(del[i].size() && del[i].back().first==mapto[j]) {
					del[i].back().second++;
				}
				else {
					del[i].push_back({mapto[j],1});
				}
			}
		}
		
		
		ZERO(dp);
		for(j=1;j<=md;j++) from[j]++;
		N--;
		while(N--) {
			ll tot=0;
			for(i=1;i<=md;i++) (tot+=from[i]*num[i])%=mo;
			to[1]=1;
			for(i=2;i<=md;i++) {
				to[i]=tot;
				FORR(r,del[i]) to[i]-=from[r.first]*r.second;
				to[i]=to[i]%mo;
				if(to[i]<0) to[i]+=mo;
			}
			swap(from,to);
		}
		
		ll ret=0;
		for(i=1;i<=md;i++) (ret+=from[i]*num[i])%=mo;
		return ret;
		
	}
}
ll dp[22][50505];
ll ok[50505],fail[50505];
ll bad[22];
ll mo=1000000007;

class DivFree {
	public:
	int dfcount(int N, int K) {
		int i,j,x;
		
		ZERO(dp);
		ZERO(bad);
		for(i=1;i<=K;i++) dp[0][i]=1;
		for(i=1;i<=20;i++) {
			for(j=1;j<=K;j++) {
				for(x=j*2;x<=K;x+=j) (dp[i][j] += dp[i-1][x])%=mo;
				bad[i]+=dp[i][j];
			}
			bad[i]%=mo;
		}
		
		ll kp=1;
		ok[0]=1;
		for(i=1;i<=N;i++) {
			ok[i]=kp=kp*K%mo;
			fail[i]=fail[i-1]*K%mo;
			
			int mi=1;
			for(j=1;j<=min(20,i-1);j++) {
				fail[i]+=mi*bad[j]*ok[i-j-1];
				mi=-mi;
			}
			fail[i]=(fail[i]%mo+mo)%mo;
			ok[i]=(ok[i]+mo-fail[i])%mo;
		}
		
		return ok[N];
		
	}
}

まとめ

446要素にまとまることは本番中わかっていたが、最適化がいまいちで10s以上かかっていた。
446要素でも頑張れば通せたのか…本番もっとあがいてみればよかった。