はじめに

畳み込みを知っていますか？僕はよく知りません

Wikipedia 先生に聞いてみると

畳み込み（たたみこみ、英: convolution）とは、関数 $g$ を平行移動しながら関数 $f$ に重ね足し合わせる二項演算である。

と言われましたが, 正直よくわかりません...

競技プログラミングの文脈においては, 畳み込みというと大抵の場合

二つの数列 $a = (a_0, \ldots ,a_n), b = (b_0, \ldots ,b_m)$ から

$$
c_i = \sum_{j=0}^i a_jb_{i-j}
$$

となる数列 $c = (c_0, \ldots ,c_{n+m})$ を求めるような演算を指すのではないでしょうか

これは, 歴史上もっとも天才なアルゴリズムとの呼び声高い Fast Fourier Transform を用いることで $\Theta((n+m)\log(n+m))$ ほどの計算量オーダーで求めることができます (すごい！)

さて, 数学ではどうなのかは知りませんが競技プログラミングでは

上記のような添え字の和での畳み込み

$$
c_k = \sum_{i+j = k} a_ib_j
$$

のほかにも,

添え字の bitwise AND での畳み込み

$$
c_k = \sum_{i \& j = k} a_ib_j
$$

添え字の bitwise OR での畳み込み

$$
c_k = \sum_{i | j = k} a_ib_j
$$

そして今回の本題である添え字の bitwise XOR での畳み込み

$$
c_k = \sum_{i \oplus j = k} a_ib_j \quad
$$

など, ほかにもさまざまなものがあるようです

なんと上記の畳み込みはすべて $\Theta((n+m)\log(n+m))$ でできることが知られています (すごすぎ！)

bitwise AND, bitwise XOR については, 高速ゼータ変換/高速メビウス変換について理解するとかなり直感的 (個人差があります) だったのですが, bitwise XOR だけ毛色が違くてあまりにもわからなかったので自分なりにまとめることにしました

本文

定義とアイデア

改めてになりますが, bitwise XOR convolution とは

$$
c_k = \sum_{i \oplus j = k} a_ib_j
$$

というような二項演算のことです ($\oplus$ はbitwise XORを表しています)

簡単のため, $n = 2^k$ ($k$ が被ってしまいましたが別物です)として,

$a = (a_0, \ldots ,a_{n-1}), b = (b_0, \ldots ,b_{n-1})$ から

$$
c_k = \sum_{i \oplus j = k} a_ib_j
$$

としたときに

$$c = a \ast b$$

と書くことにします

さて, これは各項を愚直に計算するとこで $\Theta(n^2)$ 回の演算で求めることができますが, どうすれば高速化できるでしょうか

ここで登場するのは伝家の宝刀 †分割統治† です

任意の高速化は分割統治であると考えることができる(異論は認めます)ので, 今回も使っていきましょう

どのように分割するかというと, わかりやすく真ん中で半分に分けます

つまり,

$$a^l = (a_0, \ldots ,a_{n/2-1}),\quad a^r = (a_{n/2}, \ldots ,a_{n-1}),$$

$$b^l = (b_0, \ldots ,b_{n/2-1}),\quad b^r = (b_{n/2}, \ldots ,b_{n-1})$$

から

$$c^l = (c_0, \ldots ,c_{n/2-1}),\quad c^r = (c_{n/2}, \ldots ,c_{n-1})$$

を求めます

$n = 2^k$ より, 数列の前半部分の添え字は下から $k$ 番目の bit が立っておらず, 後半部分は立っています

よって, その bit にだけ注目することで

$c^l = a^l \ast b^l + a^r \ast b^r,\quad c^r = a^l \ast b^r + a^r \ast b^l$

と計算することができました！

計算量はどうなったかというと, 長さ $n$ の数列同士の畳み込みに $T(n)$ かかるとすると半分の長さの畳み込み $4$ 回になったので

$T(n) = 4T(n/2) + \Theta(n)$

よって, Master theorem より $T(n) = \Theta(n^2)$ です！

... あれ？

分割統治をした結果, 計算量は変わりませんでした！いかがでしたか？

Karatsuba method との関係

困りました

ここで, 同じようなアイデアで高速化を実現している Karatsuba method を利用します

Karatsuba method とは, 今回で言うところの添え字の和での畳み込みを $\Theta(n^{\log_2 3})$ ほどの計算量オーダーで行うことができるアルゴリズムです

$n = 2^k$ 桁の整数が $X, Y$ あったときに, 上位 $n/2$ 桁と下位 $n/2$ 桁に分けて

$X = aR + b, Y = cR + d \quad(R = 10^{n/2})$ とすると

$$\begin{align}
XY &= (aR + b)(cR + d)\cr
&= acR^2 + (ad + bc)R + bd
\end{align}$$

となりますこのままではさきほどと同じで計算量は落ちませんが, ここで天才をします

$$ (a-b)(c-d) = ac - ad - bc + bd $$

という式が天から降りてくると, 移項して

$$ ad + bc = ac + bd - (a-b)(c-d) $$

と書けますこれを代入すると

$$\begin{align}
XY = acR^2 + (ac + bd - (a-b)(c-d))R + bd
\end{align}$$

となります項が増えてしまったようにも見えますが, $ac$ と $bd$ は二回ずつ登場しているので, 結局掛け算が起こっているのは $ac$, $bd$, $(a-b)(c-d)$ の三回

よって, この計算量は

$T(n) = 3T(n/2) + \Theta(n)$ と書けます

よって, Master theorem より $T(n) = \Theta(n^{\log_2 3})$ です！

同じことが bitwise XOR convolution でもできないでしょうか？

$c^l = a^l \ast b^l + a^r \ast b^r,\quad c^r = a^l \ast b^r + a^r \ast b^l$

という分割をしていましたが, 同じように項数を減らしたいです

ここで, おもむろに $(a^l + a^r)\ast(b^l + b^r)$ と $(a^l - a^r)\ast(b^l - b^r)$ を計算してみます分配法則をみたすので

$$ (a^l + a^r)\ast(b^l + b^r) = a^l \ast b^l + a^l \ast b^r + a^r \ast b^l + a^r \ast b^r $$

$$ (a^l - a^r)\ast(b^l - b^r) = a^l \ast b^l - a^l \ast b^r - a^r \ast b^l + a^r \ast b^r $$

と書けますよって,

$$ (a^l + a^r)\ast(b^l + b^r) + (a^l - a^r)\ast(b^l - b^r) = 2(a^l \ast b^l + a^r \ast b^r) $$

$$ (a^l + a^r)\ast(b^l + b^r) - (a^l - a^r)\ast(b^l - b^r) = 2(a^l \ast b^r + a^r \ast b^l) $$

ですしたがって,

$$\begin{align}
c^l &= a^l \ast b^l + a^r \ast b^r \cr
&= \Big((a^l + a^r)\ast(b^l + b^r) + (a^l - a^r)\ast(b^l - b^r)\Big)/2
\end{align}$$

$$\begin{align}
c^r &= a^l \ast b^r + a^r \ast b^l \cr
&= \Big((a^l + a^r)\ast(b^l + b^r) - (a^l - a^r)\ast(b^l - b^r)\Big)/2
\end{align}$$

と書けます

これは半分の長さの畳み込み二回でできるので, 再帰的に行うことで計算量は

$T(n) = 2T(n/2) + \Theta(n)$ と書けます

よって, Master theorem より $T(n) = \Theta(n\log n)$ です！

ということで, $\Theta(n\log n)$ を達成することができました

実装

これらのことをそのまま実装するとこうなります

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

#include <atcoder/modint>
using namespace atcoder;
using mint = modint998244353;

vector<mint> xor_convolution(const vector<mint> &a, const vector<mint> &b){
  int n = int(a.size());
  if(n == 1) return {a[0] * b[0]};

  int m = n/2;
  vector<mint> ap(m), am(m), bp(m), bm(m);

  for(int i = 0; i < m; ++i){
    ap[i] = a[i] + a[i+m];
    am[i] = a[i] - a[i+m];

    bp[i] = b[i] + b[i+m];
    bm[i] = b[i] - b[i+m];
  }

  vector<mint> pc = xor_convolution(ap, bp);
  vector<mint> mc = xor_convolution(am, bm);

  vector<mint> c(n);

  for(int i = 0; i < m; ++i){
    c[i] = (pc[i] + mc[i])/2;
    c[i+m] = (pc[i] - mc[i])/2;
  }

  return c;
}

int main() {
  int n; cin >> n;
  vector<mint> a(1<<n), b(1<<n);
  for(auto&e : a){ int x; cin >> x; e = x; }
  for(auto&e : b){ int x; cin >> x; e = x; }

  vector<mint> c = xor_convolution(a, b);
  for(auto&e : c) cout << e.val() << " ";
  cout << "\n";
}

Library Checker さんの Bitwise Xor Convolution で verify しました

入力

3
1 2 3 4 5 6 7 8
9 10 11 12 13 14 15 16

出力

492 488 476 472 428 424 412 408

提出結果

この実装だとかなり無駄があるため遅いですが, 無事 AC することができました！

配列使いまわせるよ！とか非再帰でできるよ！とかはぜひいろいろ調べて実装してみてください

おわりに

ということで bitwise XOR convolution について勉強して, 書いてみましたいかがでしたでしょうか

自分としてはあまりにもわからなかったものが何とか飲み込めたので満足です

参考文献