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

AOJ 包除原理

問題

日本語（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;
}