Hatena::ブログ(Diary)

simezi_tanの日記

2013-05-22

AOJ 2446 Enumeration + 高速メビウス変換まとめ

| 15:11 | AOJ 2446 Enumeration + 高速メビウス変換まとめを含むブックマーク

問題

日本語(http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=2446

制約条件

n≦20

m≦10^18


方針

包除原理。

高速メビウス変換というものを使うと、集合に対する関数f, gがあって、

f(S) = Σ[T⊆S](-1)^|S - T| g(T)

が成り立っていて、g(S)がわかっているとき、全てのSに対してf(S)を、

O(n2^n)で求められる。


この問題にこの式を適用するとき、Sは{i | 1≦i≦n}の部分集合で、

g(S) = | ∩[i∈S] Ai | とすれば、

f(S) = -(-1)^|S| | ∪[i∈S] Ai |となる。


なので、絶対値をとれば、| ∪[i∈S] Ai |が簡単に求まる。


高速ゼータ変換は結構情報があるのだけれど、その逆変換であるこれはあまり情報がなかった。

高速ゼータ変換については

のあたりを参照に。

rep(i, n) rep(j, 1 << n) if(!(j & 1 << i)) f[j] += f[j | 1 << i];

でできる。


で、高速メビウス変換は、

rep(i, n) rep(j, 1 << n) if(j & 1 << i) f[j] -= f[j ^ 1 << i];

でできて。

証明は以下の通り。


f0(S) = g(S)

fk(S) =

k∈Sでない fk-1(S)

k∈Sである fk-1(S-{K}) + fk-1(S)

とすると、

fn(S) = f(S)であることを示せばよい。


fk(S) = Σ[T⊆S, S-T⊆{1, 2, ..., k}] (-1)^|S-T| g(T)であることを示せばよい。

k = 0のときは定義から成り立つ。

k∈Sでないとき、S-T⊆{1, 2, ..., k}は{1, 2, ..., k-1}を満たすから、fk(S) = fk-1(S)

k∈Sのとき、総和の部分はS-Tがkを含むか含まないかの二つに分けられて、

Σ[T⊆S, S-T⊆{1, 2, ..., k-1}] (-1)^|S-T| g(T)

Σ[T⊆S, k∈S-T, S-T⊆{1, 2, ..., k}] (-1)^|S-T| g(T)


前者はfk-1(S)そのもの。

後者は、Sにkが含まれ、Tにkが含まれないので、

Σ[T⊆S-{K}, S-{k}-T⊆{1, 2, ..., k-1}] -(-1)^|S-{k}-T| g(T)とかけて、

これは-fk-1(S-{K})に相当する。

すなわちfk = -fk-1(S-{K}) + fk-1(S)

ソースコード

ll n, m, a[20], p[20], f[1 << 20];

int main(){
	cin >> n >> m;
	rep(i, n) cin >> a[i];
	rep(i, n) cin >> p[i];
	rep(i, 1 << n){
		ll t = 1;
		rep(j, n) if(i & 1 << j){
			ll u = a[j] / __gcd(t, a[j]);
			if(1.0 * t * u < 1.5e18) t *= u;
			else t = m + 1;
		}
		f[i] = m / t;
	}
	rep(i, n) rep(j, 1 << n) if(j & 1 << i) f[j] -= f[j ^ 1 << i];
	
	long double ans = 0;
	rep(i, 1 << n){
		long double e = 1;
		rep(j, n) e *= (i & 1 << j ? p[j] : 100 - p[j]) / 100.0l;
		ans += e * (m - abs(f[i]));
	}
	cout << fixed << setprecision(20) << ans << endl;
	
	return 0;
}

tokoharu-sakuratokoharu-sakura 2013/06/15 16:15 メビウス変換とゼータ変換の説明の式がおかしいと思います。(jを見るとすごく大きくなってしまう)

simezi_tansimezi_tan 2013/07/26 05:02 あ、本当だ、一部i, jが逆になっていますね。
下のソースコードのほうのは正しいです。

修正しておきました、ありがとうございます。

トラックバック - http://d.hatena.ne.jp/simezi_tan/20130522/1369203086
Connection: close