ブログ名

競技プログラミングやお歌のお話をする高菜です。

トップ > バタフライ演算で高速化できるものまとめ

バタフライ演算で高速化できるものまとめ

TODO: バグっているし読みづらいので書き直す TODO: 問題例追加 Editorial - AtCoder Regular Contest 132

T を集合、 V を要素数 N の全順序有限集合 V = \left \lbrace 0, \dots, N - 1 \right \rbrace、 f: T ^ { \mathcal { P } ( V ) } \rightarrow T ^ { \mathcal { P } ( V ) } を写像とします。

この f が写像族 \left ( f _ { i }: T ^ { \mathcal { P } ( V ) } \rightarrow T ^ { \mathcal { P } ( V ) } \right ) _ { i \in V } の合成 f = \bigcirc _ { i \in V } ^ { \leftarrow } f _ i = ​f _ { N - 1 } \circ f _ { N - 2 } \circ \dots \circ f _ { 0 } で表し、さらに各 f _ i \ ( i \in V ) が写像族 \left ( f _ { i, j }: T ^ 2 \rightarrow T ^ 2 \right ) _ { i \not \in S \subseteq V } の直積（すなわち \left ( f _ i (x) (S), f _ i (x) \left (S \cup \lbrace i \rbrace \right ) \right ) = f _ { i, S } \left ( x ( S ), x \left ( S \cup \lbrace i \rbrace \right ) \right ) （for x \in T ^ { \mathcal { P } ( V ) }, i \not \in S \subseteq V）が成り立つ）で表すことを目指しましょう。

例

Hadamard 変換

T が線形空間（加群）とします。 f _ { i, S } が i, S に依存せず行列 \left ((1, 1), (1, -1) \right ) であるものが Hadamard 変換です。 このとき f は M (S, T) = ( -1 ) ^ { \lvert S \triangle T \rvert } （ただし \triangle は集合の対称差）で定まる行列 M になります。

Hadamard 変換は xor-convolution (symmetric-difference convolution) に利用することができます。 a, b \in T ^ { \mathcal { P } ( V ) } をベクトルとし、 ベクトル c \in T ^ { \mathcal { P } ( V ) } を c ( S ) = \sum _ { A \triangle B = S } a (A) \cdot b (B) で定義します。 このとき (M c) (S) = \sum _ { T \subseteq V } ( -1 ) ^ { \lvert S \triangle T \rvert } c ( T ) = \sum _ { A, B \subseteq V } ( -1 ) ^ { \lvert S \triangle A \triangle B \rvert } a ( A ) \cdot b ( B ) = (M a) ( S ) \cdot ( M b ) ( S ) が成り立ち、xor-convolution は Hadamard 変換により各点積に写ることがわかります。さらに T において 2 が可逆であれば 行列 \left ((1, 1), (1, -1) \right ) は可逆です。

ゼータ変換・メビウス変換

T が線形空間（加群）とします。 f _ { i, S } が i, S に依存せず行列 \left ((1, 0), (1, 1) \right ) であるものが、ゼータ変換です。 このとき f は M (S, T) = 1 \mathop{ \mathrm { if } } S \supseteq T \mathop { \mathrm { else } } 0 で定まる行列 M になります。

ゼータ変換は and-convolution に利用することができます。 a, b \in T ^ { \mathcal { P } ( V ) } をベクトルとし、 ベクトル c \in T ^ { \mathcal { P } ( V ) } を c ( S ) = \sum _ { A \cap B = S } a (A) \cdot b (B) で定義します。 このとき (M c) (S) = \sum _ { T \supseteq S } c ( T ) = \sum _ { A, B \supseteq S } a( A ) \cdot b ( B ) = (M a) ( S ) \cdot ( M b ) ( S ) が成り立ち、and-convolution はゼータ変換により各点積に写ることがわかります。さらに行列 \left ((1, 0), (1, 1) \right ) は可逆であり、 逆変換はメビウス変換です。

また \left ((1, 0), (1, 1) \right ) の代わりに \left ((1, 1), (0, 1) \right ) を用いたものは反転束におけるゼータ変換であり、or-convolution に対応します。

離散フーリエ変換

T が 1 の N 乗根を持つ整域（可換環に拡張したい場合は『コンピュータ代数ハンドブック』参照です。）であるとします。 このとき、 \mathcal { P } ( V ) の \lbrace 0, 1 \rbrace ^ V との全単射による定まる全順序でいい感じに対応付けていい感じに bit-reverse すると、 f _ {i, S} が行列 \left ( (1, w ^ k), ( 1, - w ^ k) \right ) （ただし k はいい感じ）みたいにいい感じに表せます。私は疲れました。

謎アダマール変換もどき

本題です。こちらのユーザー解説に「アダマール変換に用いられるアダマール行列 と同様の再帰構造を持ち」とありますが、それはなんでしょう。 見た目は行列ですが、書いてある行列がそのまま当記事の形のバタフライ演算に分解するのは私にはできませんでした。

atcoder.jp

結論、 f _ { i, S } ( x, y ) = ( x + y, x + 2 ^ i - y ) と定義すると、ユーザー解説の g, D _ 2 に対して、 g = f (D _ 2) が成り立ちます。

感想

f _ { i, S } が線形変換でないケースもあるのですね。泣きました。

［追記］意味がわからないよ！ という方向け

申し訳ないです。特に新しくもないよく知られたこと（少なくとも私はライブラリがほぼこの形でした。）を、私が理解しやすいように形式的にお話してしまっただけで、要はこの関数に渡せる形にしましょうということです。

fn butterfly<T>(a: &mut [T], f: impl FnMut(usize, usize, &mut T, &mut T)) {
    let n = a.len();
    assert!(n.is_power_of_two());
    let lg = n.trailing_zeros() as usize;
    for i in 0..lg {
        for bs in (0..n).filter(|&bs| bs >> i & 1 == 0) {
            let (left, right) = a.split_at_mut(bs | 1 << i);
            f(i, bs, &mut left[bs], &mut right[0])
        }
    }
}

たとえば、ゼータ変換：

butterfly(&mut a, |_, _, x, y| *y ^= *x);

謎アダマール変換もどき：

fn transform(a: &mut [u32]) {
    butterfly(a, |i, _, x, y| {
        let z = *x + *y;
        let w = *x + (1 << i) - *y;
        *x = z;
        *y = w;
    });
}

Rust でも (*x, *y) = (*x + *y, *x + (1 << i) - *y) がしたいなと言う気持ちになりますね。

離散フーリエ変換（bit-reverse 部分以外）とかアダマール変換もこの関数を呼び出すだけでできると思うのですがフーリエの方計算大変なのでまた気が向いたらで。