naoya_t@hatenablog

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

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

高速ゼータ変換・高速メビウス変換が気になっていたところにタイムリーに先日の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]

メモ

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

ゼータ関数メビウス関数

Wikipediaを見てると

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

とか

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

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

追記

デモ実装

#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さんの記事が分かりやすかった。