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

高速ゼータ変換・高速メビウス変換が気になっていたところにタイムリーに先日のARC 100のE問題が来て（ちゃんと解けなかったので）折角なのでこの機会にマスターしようと思い競プロ有識者の皆さんのブログやツイートを漁ってみました。

・・・・・・

ゼータ変換／メビウス変換が互いに逆変換の関係なのはわかったんだけど…どっちがどっち？

というか、4種類あるんだけど…どれがどれ？

$O(N2^N)$ なループ rep(i,N) rep(j,1<<N) *1の中の実装が4種類

if (!(j & (1<<i))) f[j] += f[j|(1<<i)];
if (j & (1<<i)) f[j] += f[j^(1<<i)];
if (!(j & (1<<i))) f[j] -= f[j|(1<<i)];
if (j & (1<<i)) f[j] -= f[j^(1<<i)];

どれがゼータ変換で、どれがメビウス変換？

ゼータ変換／メビウス変換それぞれにS⊆T, T⊆Sの2種類ずつあるということみたいで混乱するけれど、
（けんちょん先生が全ての謎を解いてくれるのを待ちつつ）ネットの集合知やWikipediaの助けを借りて答えを探してみよう…

7/3 20:26 誰もがこれはゼータ変換だと言うタイプが(1)に来るように (1)(2),(3)(4)の順序を入れ替えました。

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

ゼータ変換(1)。（これはみんなゼータ変換って言ってる）
g(S)=「Sを部分集合として持つ全ての集合Tについてのf(T)の総和」
$\displaystyle g(s)=\sum_{s\subseteq T}{f(T)}$
この逆変換がメビウス変換(1)。

 f[0]: 00000001 => 11111111  g[000]=f[000]+f[001]+f[010]+f[011]+f[100]+f[101]+f[110]+f[111]
 f[1]: 00000010 => 10101010  g[001]=f[001]+f[011]+f[101]+f[111]
 f[2]: 00000100 => 11001100  g[010]=f[010]+f[011]+f[110]+f[111]
 f[3]: 00001000 => 10001000  g[011]=f[011]+f[111]
 f[4]: 00010000 => 11110000  g[100]=f[100]+f[101]+f[110]+f[111]
 f[5]: 00100000 => 10100000  g[101]=f[101]+f[111]
 f[6]: 01000000 => 11000000  g[110]=f[110]+f[111]
 f[7]: 10000000 => 10000000  g[111]=f[111]

iwiwi 2012-04-22, 高速ゼータ変換
todo314 2012-06-14, 高速ゼータ変換*2
pekempey 2016-10-30, （"上位集合に対する"ゼータ変換、に相当か）
prime 2016-12-01, 高速ゼータ変換
hamayanhamayan 2017-05-20, 高速ゼータ変換「状態xを含む状態yを全列挙」
omochan 2017-09-26, "正しくない"ゼータ変換
kyopro_friends 2018-07-01, 高速ゼータ変換
sugarknri 2018-07-03, ゼータ変換

(2) `rep(i,N) rep(j,1<<N) if (j&(1<<i)) f[j]+=f[j^(1<<i)];`

ゼータ変換(2)。（~~これをメビウス変換と呼ぶ人もいるが恐らく誤り~~？？これをメビウス変換としている論文 (thanks to nishioさん）
g(S)=「Sの全ての部分集合Tについてのf(T)の総和」
$\displaystyle g(S)=\sum_{t\subseteq S}^{}{f(t)}$
この逆変換がメビウス変換(2)。

 f[0]: 00000001 => 00000001  g[000]=f[000]
 f[1]: 00000010 => 00000011  g[001]=f[001]+f[000]
 f[2]: 00000100 => 00000101  g[010]=f[010]+f[000]
 f[3]: 00001000 => 00001111  g[011]=f[011]+f[010]+f[001]+f[000]
 f[4]: 00010000 => 00010001  g[100]=f[100]+f[000]
 f[5]: 00100000 => 00110011  g[101]=f[101]+f[100]+f[001]+f[000]
 f[6]: 01000000 => 01010101  g[110]=f[110]+f[100]+f[010]+f[000]
 f[7]: 10000000 => 11111111  g[111]=f[111]+f[110]+f[101]+f[100]+f[011]+f[010]+f[001]+f[000]

todo314 2012-06-14, (2)のSとTの関係を逆にした場合
pekempey 2016-10-30, ”下位集合に対する”ゼータ変換
prime 2016-12-01, 高速メビウス変換（←多分誤り）
omochan 2017-09-26, "正しい" ゼータ変換
kyopro_friends 2018-07-01, 高速メビウス変換（primeさんのブログを引用してる）
sugarknri 2018-07-03, メビウス変換

(3) `rep(i,N) rep(j,1<<N) if (!(j&(1<<i))) f[j]-=f[j|(1<<i)];`

メビウス変換(1)。ゼータ変換(1)の逆変換。包除原理。
$\displaystyle f(s)=\sum_{s\subseteq T}{(-1)^{|T\setminus s|}g(T)}$

 f[0]: 11111111 => 00000001  f[000]=(-1)^(3-0)*g[111] + (-1)^(2-0)*(g[110]+g[101]+g[011]) + (-1)^(1-0)*(g[100]+g[010]+g[001]) + (-1)^(0-0)*g[000]
 f[1]: 10101010 => 00000010  [f001]=(-1)^(3-1)*g[111] + (-1)^(2-1)*(g[101]+g[011]) + (-1)^2*g[001]
 f[2]: 11001100 => 00000100  f[010]=(-1)^(3-1)*g[111] + (-1)^(2-1)*(g[110]+g[011]) + (-1)^2*g[010]
 f[3]: 10001000 => 00001000  f[011]=(-1)^(3-2)*g[111] + (-1)^(2-2)*g[011]
 f[4]: 11110000 => 00010000  f[100]=(-1)^(3-1)*g[111] + (-1)^(2-1)*(g[110]+g[101]) + (-1)^(1-1)*g[100]
 f[5]: 10100000 => 00100000  f[101]=(-1)^(3-2)*g[111] + (-1)^(2-2)*g[101]
 f[6]: 11000000 => 01000000  f[110]=(-1)^(3-2)*g[111] + (-1)^(2-2)*g[110]
 f[7]: 10000000 => 10000000  f[111]=(-1)^(3-3)*g[111]

pekempey 2016-10-30, （"上位集合に対する"メビウス変換、に相当か）

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

メビウス変換(2)。ゼータ変換(2)の逆変換。包除原理。（ゼータ変換(2)をメビウス変換と呼ぶ人はこれを「メビウス逆変換」と呼んでいるが恐らく誤り）
$\displaystyle f(S)=\sum_{t\subseteq S}{(-1)^{|S\setminus t|}g(t)}$

 f[0]: 00000001 => 00000001  f[000]=(-1)^(1-1)*g[000]
 f[1]: 00000011 => 00000010  f[001]=(-1)^(1-1)*g[001] + (-1)^(1-0)*g[000]
 f[2]: 00000101 => 00000100  f[010]=(-1)^(1-1)*g[010] + (-1)^(1-0)*g[000]
 f[3]: 00001111 => 00001000  f[011]=(-1)^(2-2)*g[011] + (-1)^(2-1)*(g[010]+g[001]) + (-1)^(2-0)*g[000]
 f[4]: 00010001 => 00010000  f[100]=(-1)^(1-1)*g[100] + (-1)^(1-0)*g[000]
 f[5]: 00110011 => 00100000  f[101]=(-1)^(2-2)*g[101] + (-1)^(2-1)*(g[100]+g[001]) + (-1)^(2-0)*g[000]
 f[6]: 01010101 => 01000000  f[110]=(-1)^(2-2)*g[110] + (-1)^(2-1)*(g[100]+g[010]) + (-1)^(2-0)*g[000]
 f[7]: 11111111 => 10000000  f[111]=(-1)^(3-3)*g[111] + (-1)^(3-2)*(g[110]+g[101]+g[011]) + (-1)^(3-1)*(g[100]+g[010]+g[001]) + (-1)^(3-0)*g[000]

lan496 2015-05-14, 高速メビウス変換
pekempey 2016-10-30, "下位集合に対する"メビウス変換
prime 2016-12-01, 高速メビウス逆変換（←多分誤り）
hamayanhamayan 2017-05-20, メビウス変換「状態xの部分状態yを全列挙」
omochan 2017-09-26, "正しい"ゼータ変換の逆関数としてのメビウス変換
sugarknri 2018-07-03, メビウス逆変換

メモ

(1)-(3), (2)-(4) はそれぞれ逆変換の関係。
(1)-(2), (3)-(4) はそれぞれ、ビットを逆順にしてからもう片方の変換を行いその結果を逆順にしたものと等しい。

ゼータ関数とメビウス関数

Wikipediaを見てると

「（幾何学における平面上の）メビウス変換」と「数論的メビウス変換」があったり
「数論的メビウス関数」と「組合せ論的メビウス関数」があったり

でややこしいのだけれど、今ゼータ変換とかメビウス変換とか言ってるのは

隣接代数 (順序理論) - Wikipedia

とか

の辺りの話らしい。隣接代数の項目にある

ゼータ関数を掛ける→積分に相当
メビウス関数を掛ける→微分に相当（メビウス反転と呼ばれる。"このときのメビウス反転は包含と切除の原理と呼ばれるものである"とある）

というのがそれっぽい。
この観点で言うなら前掲の(1),(2)はゼータ変換だし(3),(4)はメビウス変換ということで良さそう。
（どちらも半順序集合を対象としていて、順序を決める包含関係 ⊆ がどちらに向いていようが成り立つ、ということで良いのかな）

追記

@naoya_t https://t.co/9taQEeaOq2
>ゼータ変換(2)。（これをメビウス変換と呼ぶ人もいるが恐らく誤り）

に関して、https://t.co/QGl7orwgdI
で紹介した論文ではこれをメビウス変換と呼んでます
— nishio hirokazu (@nishio) 2020年10月23日

「メビウス変換」多義語なので要注意。
弊記事でゼータ変換(2)としているものをメビウス変換と呼んでいる論文もある。
"この概念を指すもっと誤解の余地のない言葉として「Sum over Subsets(SOS)」がある" https://codeforces.com/blog/entry/45223

デモ実装

#include <bits/stdc++.h>
using namespace std;

#define rep(var,n)  for(int var=0;var<(n);++var)

ostream& operator<<(ostream &s, vector<int> v)
{
  int cnt = v.size();
  s << "[ ";
  for (int i=0; i<cnt; i++) {
        if (i > 0) s << ", ";
        s << bitset<16>(v[i]);
  }
  return s << " ]";
}

inline int ilog2(int x) {
    int y = 0;
    while (x > 1) {
        x /= 2;
        ++y;
    }
    return y;
}

template <typename T>
vector<T> zeta(vector<T> f, int type) {
    int N = f.size(), K = ilog2(N);
    for (int i=0,b=1; i<K; ++i,b<<=1) { // m = (1<<i)
        for (int j=0; j<N; ++j)
            switch (type) {
            case 1:
                if (!(j & b)) f[j] += f[j | b];
                break;
            case 2:
                if (j & b) f[j] += f[j ^ b];
                break;
            case 3:
                if (!(j & b)) f[j] -= f[j | b];
                break;
            case 4:
                if (j & b) f[j] -= f[j ^ b];
                break;
            }
    }
    return f;
}

void pp(vector<int>& a, vector<int>& b) {
    rep(i, a.size()) {
        cout << " f[" << i << "]: "
             << bitset<8>(a[i])
             << " => " << bitset<8>(b[i])
             << endl;
    }
}

int main() {
    vector<int> a(1<<3);
    rep(i,1<<3) a[i] = 1<<i;

    rep(type, 4) {
        vector<int> b = zeta(a, 1+type);
        cout << "type " << (1+type) << ":" << endl;
        pp(a, b);
        cout << endl;

        vector<int> c = zeta(b, 1+(type^2));
        pp(b, c);
        cout << endl;
    }

    return 0;
}

ゼータ変換を使う問題の例

ARC 100 E - Or Plus Max (700)

ここではゼータ変換のtype2のやつ
intの足し算の代わりに、2要素を持たせて上位2件をマージする感じで
→AC
https://beta.atcoder.jp/contests/arc100/submissions/2781917
↑これを書いた時点では混同していて、ゼータ変換(1)をゼータ変換、ゼータ変換(2)をメビウス変換だと思ってる

メビウス変換と包除原理

例として、100までの自然数に含まれる2,3,5の倍数について考えてみる。

a[1] = 2の倍数の個数 = 50
a[2] = 3の倍数の個数 = 33
a[3] = 6(=2*3)の倍数の個数 = 16
a[4] = 5の倍数の個数 = 20
a[5] = 10(=2*5)の倍数の個数 = 10
a[6] = 15(=3*5)の倍数の個数 = 6
a[7] = 30(=2*3*5)の倍数の個数 = 3

ということで

vector<int> a { 0, 50, 33, 16, 20, 10, 6, 3 };

を用意して、これにメビウス変換(2)をかけると

 f[0]: 0 -> 0
 f[1]: 50 -> 50
 f[2]: 33 -> 33
 f[3]: 16 -> -67
 f[4]: 20 -> 20
 f[5]: 10 -> -60
 f[6]: 6 -> -47
 f[7]: 3 -> 74

が得られる。それぞれの絶対値を見ると

|f[1]| = 2の倍数の個数 = 50
|f[2]| = 3の倍数の個数 = 33
|f[3]| = 2または3（少なくとも一方）の倍数の個数 = 67
|f[4]| = 5の倍数の個数 = 20
|f[5]| = 2または5（少なくとも一方）の倍数の個数 = 60
|f[6]| = 3または5（少なくとも一方）の倍数の個数 = 47
|f[7]| = 2,3,5（少なくともどれか1つ）の倍数の個数 = 74

になっている。やっていることは包除原理を用いて積集合のサイズから和集合のサイズを求める計算である。
（厳密には絶対値というか $(-1)^{|S|+1}$ を掛けたもの。←メビウス変換の式と包除原理の式を並べて考えれば分かる）

*1:配列のサイズを2^Nとする

*2:高速ゼータ変換の中で起こっていることはtodo314さんの記事が分かりやすかった。

naoya_t@hatenablog

いわゆるチラシノウラであります

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

というか、4種類あるんだけど…どれがどれ？

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

(2) `rep(i,N) rep(j,1<<N) if (j&(1<<i)) f[j]+=f[j^(1<<i)];`

(3) `rep(i,N) rep(j,1<<N) if (!(j&(1<<i))) f[j]-=f[j|(1<<i)];`

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

メモ

ゼータ関数とメビウス関数

追記

デモ実装

ゼータ変換を使う問題の例

ARC 100 E - Or Plus Max (700)

メビウス変換と包除原理

というか、4種類あるんだけど…どれがどれ？

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

(2) rep(i,N) rep(j,1<<N) if (j&(1<<i)) f[j]+=f[j^(1<<i)];

(3) rep(i,N) rep(j,1<<N) if (!(j&(1<<i))) f[j]-=f[j|(1<<i)];

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

メモ

ゼータ関数とメビウス関数

追記

デモ実装

ゼータ変換を使う問題の例

ARC 100 E - Or Plus Max (700)

メビウス変換と包除原理

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

(2) `rep(i,N) rep(j,1<<N) if (j&(1<<i)) f[j]+=f[j^(1<<i)];`

(3) `rep(i,N) rep(j,1<<N) if (!(j&(1<<i))) f[j]-=f[j|(1<<i)];`

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