bitwise XOR convolution と Karatsuba method

# はじめに

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

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)$ を達成することができました

## 実装

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

```c++
#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](https://judge.yosupo.jp/problem/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
```

[提出結果](https://judge.yosupo.jp/submission/82312)


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

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

# おわりに

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

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

# 参考文献

https://ja.wikipedia.org/wiki/%E7%95%B3%E3%81%BF%E8%BE%BC%E3%81%BF

https://kazuma8128.hatenablog.com/entry/2018/05/31/144519

https://qiita.com/Euglenese/items/260f9ddf513f772d7e42

https://cympfh.cc/aiura/hadamard-xor-convolution

https://dic.nicovideo.jp/a/%E3%81%84%E3%81%8B%E3%81%8C%E3%81%A7%E3%81%97%E3%81%9F%E3%81%8B%3F

http://www.finetune.co.jp/~lyuka/technote/fract/bmul.html