高速ゼータ変換、高速メビウス変換

Inclusion and Exclusion

包除原理で活躍する二つのアルゴリズムです。この2つの記事がわかりやすかったです。
高速メビウス変換について - 篠突く雨の日記
ゼータ変換とメビウス変換 - pekempeyのブログ

メビウス:∩→∪
ゼータ :∪→∩と覚えればいいです。たぶんメビウスの方がよく使うと思います。なぜなら包除原理は∩を∪に変換する公式と考えられるからです。

例えば2,3,7のうちどれかの倍数(集合でいうと∪の関係)である数が100までの自然数で何個あるか考えましょう。
2かつ3の倍数(つまり集合でいうと∩の関係)である数の個数などは簡単に求まります。これをメビウス変換に渡せば求めるべき値が得られます。
ゼータ変換はこの逆で∪で表される関係を∩にしてくれます。

アルゴリズムはbitdpで、dp[i][S]:(i-1番目までのbitは包含関係が満たされればなんでもいいけどi番目からはSと等しい要素の和)でやります。少しワーシャルフロイドに似てるかな…?
たとえばdp[2][1110]は要素0010,0110,1010,1110の和になっています。これを配列を使いまわして書いたのが下のコードになります。

メビウス変換とゼータ変換はコード見てもらえればわかる通り、違うのは+と-だけなので、どっちかを覚えればよいです。
集合Sに含まれるTの和ではなく、Sを含むTの和でも少し変えるだけでよいです。

FFTみたいな使い方もできます。
(f*g)(X)=Σ(S+T=X)f(S)g(T)について
ζ(f*g)(X)=ζf(S)*ζf(T)が成立します。ここでSとTが同じ要素を持っていてもよいことに注意です。
これはFFTの、
(f*g)(n)=Σ(k<=n)f(k)g(n-k)について
FFT(f*g)(n)=FFTf(n)*FFTg(n)が成立するのととても似ています。

記事を読んでるときにS=110(2),T=100(2)のときS⊃TなのかT⊃Sなのかこんがらがりましたが、正しいのは前者です。

vector<int> zeta(vector<int> f) { //get F(S)=Σ(T⊂S)G(T)
	int n = sz(f);
	for(int i = 0; (1 << i) < n; i++) {
		rep(j, 0, n) if(j & (1 << i)) f[j] += f[j ^ (1 << i)];
	}
	return f;
}
/*
vector<int> zeta(vector<int> f) { //get F(S)=Σ(S⊂T)G(T)
	int n = sz(f);
	for(int i = 0; (1 << i) < n; i++) {
		rep(j, 0, n) if(!(j & (1 << i))) f[j] += f[j | (1 << i)];
	}
	return f;
}
*/

vector<int> moebius(vector<int> f) { //inverse of the first zeta
	int n = sz(f);
	for(int i = 0; (1 << i) < n; i++) {
		rep(j, 0, n) if(j & (1 << i)) f[j] -= f[j ^ (1 << i)];
	}
	return f;
}

int N;
int d[4] = {2, 3, 7, 11};
int cnt[110];

void solve() {
	N = 4;
	int res = 0;
	rep(i, 1, 100 + 1) {
		rep(j, 0, N) {
			if(i % d[j] == 0) {
				res++;
				break;
			}
		}
	}
	debug(res);
	rep(bit, 0, (1 << N)) {
		int a = 1;
		rep(i, 0, N) {
			if(bit & (1 << i)) {
				a *= d[i];
			}
		}
		cnt[bit] = 100 / a;
	}
	cnt[0] = 0;
	debug(vi(cnt, cnt + (1 << N)));
	vi vec = moebius(vi(cnt, cnt + (1 << N)));
	vi vec2 = zeta(vec);
	rep(bit, 0, (1 << N)) {
		debug(bitset<4>(bit), vec[bit], vec2[bit]);
	}
}