👊パンチ👊ゼータ変換👊パンチ👊

累積和はわかるけど、約数包除はややこしくてわからない！

競プロでgcd (最大公約数)・lcm (最小公倍数) をなんやかんやする問題が出たときに、何も考えずにゼータ変換のライブラリを貼って解く（「殴る」）ことができることがあるのを紹介します。ABC-EやABC-Fを速解きして暖色パフォが出せたらうれしいですね
※ ネタバレを多分に含むのでご注意ください

はじめに
方法・実装
実例
まとめ
参考

はじめに

約数集合での上位集合に対するゼータ変換 (以下ゼータ変換とします) とは、各素数が各添字に対応する多次元の累積和をとることで、添字gcd畳み込みなどの操作が高速にできたりする便利な変換です*1。愚直にやれば $O(N^2)$ になる色々な操作が $O(N\log\log N)$ になって速いです。これを使って、簡単な式変形で様々な問題を解くことができたりします。

方法・実装

ゼータ変換

具体的には、長さ $N$ の数列 $A$ に対して、
$\displaystyle B_i = \sum_{i|j} A_j$
で定義される長さ $N$ の数列 $B$ に変換する操作がゼータ変換です。実装はエラトステネスの篩同様に

vector<bool> sieve(a.size(), true);
for (size_t i =  2; i < a.size(); i++) {
    if (!sieve[i]) continue;
    for (size_t j = (a.size()-1)/i; j > 0; j--) {
        sieve[i*j] = false;
        a[j] += a[i*j];
    }   
}

のように実装すると、計算量 $O(N\log\log N)$ で変換することができます。
素数(≒次元)ごとに指数の大きい方から累積和をとるようなアルゴリズムになっていることが分かると思います。
参考: 高速ゼータ変換の約数版 O(N log(log(N))) - noshi91のメモ
参考: ゼータ変換・メビウス変換を理解する - Qiita

メビウス変換

メビウス変換 (ゼータ変換した配列を元に戻す変換) も累積和を元の配列に戻すように逆順に引いていくと

vector<bool> sieve(a.size(), true);
for (size_t i =  2; i < a.size(); i++) {
    if (!sieve[i]) continue;
    for (size_t j = 1; i*j < a.size(); j++) {
        sieve[i*j] = false;
        a[j] -= a[i*j];
    }   
}

のように実装できます。

GCD畳み込み

長さ $N$ の配列 $A,\ B$ があったときに、これをゼータ変換してそれぞれ $\mathcal{A},\ \mathcal{B}$ とします。 $\mathcal{A},\ \mathcal{B}$ の各成分をかけてメビウス変換すると添字gcdで畳み込みができます。具体的には、
$\displaystyle\begin{align} \mathcal{A}_i\mathcal{B}_i &= \left( \sum_{i|j} A_j \right) \left( \sum_{i|k} B_k \right) \\ &= \sum_{i|j} \sum_{i|k} A_j B_k \\ &= \sum_{i|\mathrm{gcd}(j,k)} A_j B_k \\ &= \sum_{i|g} \sum_{\mathrm{gcd}(j,k)=g} A_j B_k \end{align}$
となるため、 $\left\{\mathcal{A}_i\mathcal{B}_i\right\}_i$ をメビウス変換すると配列 $C$ が
$\displaystyle\begin{align} C_i = \sum_{\mathrm{gcd}(j,k)=i} A_j B_k \end{align}$
のように得られます。この記事ではこれをgcd畳み込みと呼ぶことにします。

LCM畳み込み

ゼータ/メビウス変換をする際の累積和をとる向きを逆にすると添字lcmでの畳み込みができます。先ほどの配列 $C$ が
$\displaystyle\begin{align} C_i = \sum_{\mathrm{lcm}(j,k)=i} A_j B_k \end{align}$
のように得られます。しかし、lcmの値は大きくなって配列に収まりきらない場合が多く、あまり使われません。

実例

ABC106-B 105

atcoder.jp
$B_i=1$ を下位集合に対してゼータ変換することで $N$ 以下のすべての正整数の約数の個数を $O(N\log\log N)$ で求められます。

ただ、制約が $N\le 200$ で全探索 $O(N^2)$ が通るのでゼータ変換を使うのは完全に👊パンチ👊です。制約が $N\le 10^6$ でも多分通ります。

ABC172-D Sum of Divisors

atcoder.jp
式を変形すると

$\displaystyle\begin{align} \sum_{K=1}^N K \times f(N) =& \sum_{K=1}^N K \times \#\left\{ i \in \mathbb{Z} \cap [1,N] \mid i|K \right\} \\ =& \sum_{K=1}^N K \sum_{i|K} 1 \\ \end{align}$

となります。 $\sum_{i|K}1$ は下位集合に対するゼータ変換を使えば $O(N\log\log N)$ で求まるのでACできます。ただし、空間 $O(N)$ なので定数倍は重いです。

ABC170-D Not Divisible

atcoder.jp
長さ $V:=\max A$ の配列 $B$ を $B_i=\#\{j \mid A_j = i \}$ と構築すると、求める性質は

$\displaystyle\begin{align} 「\sum_{j|i} B_j = 1」かつ「B_i = 1」 \end{align}$

と言い換えることができ、これも下位集合に対するゼータ変換を使えば $O(N\log\log N)$ で求まるのでACできます。

ABC162-E Sum of gcd of Tuples (Hard)

atcoder.jp
$K\le10^5$ なので、長さ $K$ の配列 $\mathcal{B}$ を

$\displaystyle\begin{align} \mathcal{B}_i &= \#\left\{ (A_1,\cdots,A_N) \middle| \mathrm{gcd}(A_1,\cdots,A_N)\text{は}i\text{の倍数} \right\} \\ &= \sum_{i|g} \#\left\{ (A_1,\cdots,A_N) \middle| \mathrm{gcd}(A_1,\cdots,A_N)=g \right\} \\ \end{align}$

のようにすると、 $\mathcal{B}$ をメビウス変換すればgcdが $i$ になる数列の数が求まります。そこでgcdを $i$ の倍数にするには各要素に $i$ の倍数を使えばいいので

$\displaystyle\begin{align} \mathcal{B}_i &= \left\lfloor \frac{K}{i} \right\rfloor^N \end{align}$

が求まります。最後に $\mathcal{B}$ をメビウス変換で $B$ に変換して

$\displaystyle\begin{align} \sum_{i=1}^K i \cdot B_i \end{align}$

を求めて出力すればACです。 $B$ の構築に繰り返し二乗法で $O(K\log N)$ 、メビウス変換で $O(K\log\log K)$ 、最後の答えが $O(K)$ なので $O(K(\log N+\log\log K))$ で間に合います。

参考: [AtCoder 参加感想] 2020/04/12:ABC 162 | maspyのHP
参考: AtCoder頻出？貼るだけ高速ゼータ・メビウス変換＆約数拡張 - ゆきノート（仮）

追記(2020/7/25)：
熨斗袋さんからさらにパンチ力の高い解法を頂きました

Sum of gcd of Tuples は更に👊出来ますね

まず 3 つの列 a, b, c を gcd 畳み込みすると、ゼータ変換を Z として
Z⁻¹(Z( Z⁻¹(Z(a)⚪︎Z(b)) )⚪︎Z(c))
になるんですがこれは
Z⁻¹(Z(a)⚪︎Z(b)⚪︎Z(c))
なので、結局複数の列を一気に畳めることが分かります
— 応用微積 (@noshi91) 2020年7月24日

列 S を S[1] = S[2] = … S[K] = 1 とすると、
Z⁻¹(Z(S)⚪︎Z(S)⚪︎ … (K 個) … ⚪︎Z(S))
を計算できればいいので、S をゼータ変換して各要素を k 乗してメビウス変換すると OK です
— 応用微積 (@noshi91) 2020年7月24日

問題の式は
$\displaystyle\begin{align} &\sum_{1\le a_1, a_2, \ldots, a_N \le K} \mathrm{gcd}(a_1, a_2, \ldots, a_N) \\ =& \sum_{g=1}^K \sum_{\mathrm{gcd}(a_1, a_2, \ldots, a_N)=g} g \\ =& \sum_{g=1}^K g \sum_{\mathrm{gcd}(a_1, a_2, \ldots, a_N)=g} 1 \end{align}$
と変形できます。正整数のgcdは可換モノイドなので、gcd畳み込みも同様です。そうするとすべての項が1の配列 $B$ をN個畳み込めば
$\displaystyle\begin{align} (\underbrace{B\ast B \ast \cdots \ast B}_{N \text{times}})_g &= \sum_{\mathrm{gcd}(a_1, a_2, \ldots, a_N)=g} B_{a_1}B_{a_2}\cdots B_{a_N} \\ &= \sum_{\mathrm{gcd}(a_1, a_2, \ldots, a_N)=g} 1 \end{align}$
が得られ、最初の変形した式に代入すれば値が求まることが分かります。
$\underbrace{B\ast B \ast \cdots \ast B}_{N \text{times}}$ は、 $B$ をゼータ変換したのち各項それぞれ $N$ 乗してメビウス変換をすれば得られるので、解けました。計算量は $O(K(\log N + \log\log K))$ です。
（追記おわり）

AGC038-C LCMs

atcoder.jp
まず、与えられた式を計算しやすいように変形すると

$\displaystyle\begin{align} &\sum_{i=0}^{N-2} \sum_{j=i+1}^{N-1} \textrm{lcm}(A_i, A_j) \\ =& \sum_{0\le i \lt j \lt N} \textrm{lcm}(A_i, A_j) \\ =& \frac12 \left[ \sum_{i=0}^{N-1} \sum_{j=0}^{N-1} \mathrm{lcm}(A_i,A_j) - \sum_{i=0}^{N-1} \mathrm{lcm}(A_i,A_i) \right] \\ =& \frac12 \left[ \sum_{i=0}^{N-1} \sum_{j=0}^{N-1} \mathrm{lcm}(A_i,A_j) - \sum_{i=0}^{N-1} A_i \right] \end{align}$

となります。二重 $\sum$ のところが $O(N^2)$ かかってしまうのでさらに変形します。

$\displaystyle\begin{align} &\sum_{i=0}^{N-1} \sum_{j=0}^{N-1} \mathrm{lcm}(A_i,A_j) \\ =& \sum_{i=0}^{N-1} \sum_{j=0}^{N-1} \frac{A_iA_j}{\mathrm{gcd}(A_i,A_j)} \\ \end{align}$

ここで、 $V:=\max_i\{A_i\}$ とすると、 $1\le\mathrm{gcd}(A_i,A_j)\le V$ より

$\displaystyle\begin{align} &\sum_{i=0}^{N-1} \sum_{j=0}^{N-1} \frac{A_iA_j}{\mathrm{gcd}(A_i,A_j)} \\ =& \sum_{g=1}^{V} \sum_{\mathrm{gcd}(A_i,A_j)=g} \frac{A_i A_j}{g} \\ =& \sum_{g=1}^{V} \frac1g \sum_{\mathrm{gcd}(A_i,A_j)=g} A_i A_j \end{align}$

となります。ここで長さ $V$ の配列 $B$ を $B_i = \sum_{A_j=i} A_j = i\cdot \#\left\{ j \middle| A_j = i \right\}$ としておき、 $B$ 同士のgcd畳み込みをすれば

$\displaystyle\begin{align} \left\{ B \ast B \right\}_i &= \sum_{\mathrm{gcd}(j,k)=i} B_j B_k \\ &= \sum_{\mathrm{gcd}(A_j,A_k)=i} A_j A_k \end{align}$

となり、元の式が求まります。
以上を実装すると、 $B$ の構築に $O(N+V)$ 、畳み込みに $O(V\log\log V)$ 、最終的な値の計算に $O(N+V)$ かかるので時間計算量は $O(N+V\log\log V)$ となります。問題の制約は $N\le 2\times 10^5,\ V\le 10^6$ なので制限時間の2秒に間に合います。

参考: 添え字 gcd での畳み込みで AGC038-C を解く - noshi91のメモ
参考: agc038_c - wiki.kimiyuki.net

ABC020-D LCM Rush

atcoder.jp
これも少し工夫することでゼータ変換で解けます。求める値を式変形すると

$\displaystyle\begin{align} \sum_{i=1}^N \textrm{lcm}(i,K) &= \sum_{i=1}^N \frac{iK}{\textrm{gcd}(i,K)} \\ &= \sum_{g=1}^{\min(N,K)} \sum_{\textrm{gcd}(i,K)=g} \frac{iK}{g} \\ &= K \sum_{g=1}^{\min(N,K)} \frac{1}{g} \sum_{\textrm{gcd}(i,K)=g, \\ 1\le i \le N} i \\ \end{align}$

となり、ここで $\textrm{gcd}(i,K)=g$ が成立するような $g$ はKの約数に限られることに注意すると、

$\displaystyle\begin{align} & K \sum_{g=1}^{\min(N,K)} \frac{1}{g} \sum_{\textrm{gcd}(i,K)=g, \\ 1\le i \le N} i \\ =& K \sum_{g|K} \frac{1}{g} \sum_{\textrm{gcd}(i,K)=g, \\ 1\le i \le N} i \\ \end{align}$

ここで $\sum_{\textrm{gcd}(i,K)=g, 1\le i \le N} i$ の部分をゼータ変換して考えると

$\displaystyle\begin{align} & \sum_{g|\textrm{gcd}(i,K), \\ 1\le i \le N} i \\ =& \sum_{g|i, 1\le i \le N} i ~~~~ (すでにg|Kに限っているため) \\ =& \frac12 \left\lfloor \frac{N}{g} \right\rfloor \left( \left\lfloor \frac{N}{g} \right\rfloor + 1 \right) g \end{align}$

これをすべての $g|K$ について求めておき、 $K$ の約数だけを添字としてメビウス変換を行うと $\sum_{\textrm{gcd}(i,K)=g, 1\le i \le N} i$ が求まります。これを元の式に代入すれば求まって、ACです。 $K$ に含まれる相異な素因数の個数を $p(K)$ 、 $K$ の約数の個数を $d(K)$ とすると計算量は $O(p(K)d(K)\log d(K))$ となり、十分高速で2sに間に合います。

「みんなのプロコン 2018」決勝-A Uncommon

atcoder.jp

$i=1,2,\ldots,M$ に対して

$\displaystyle\begin{align} \# \{ j \mid \mathrm{gcd}(a_j, i) = 1 \} \end{align}$

をそれぞれ求める問題です。これは、

$\displaystyle\begin{align} B_g &:= \# \{ j \mid a_j = g \}, \\ C_g &:= \delta_{i,g} \end{align}$

で定義される配列 $B,C$ のgcd畳み込みの $g=1$ の項をそれぞれ求めればいいです。ここでこれらをゼータ変換することを考えると、 $B$ は $i$ に依存しないのでそのままゼータ変換して $\mathcal{B}$ とすればいいのですが、 $C$ が $i$ によって変わるのでこれをどうにかする必要があります。しかし、 $C$ をゼータ変換して $\mathcal{C}$ とすると、 $g|i$ を満たすとき $\mathcal{C}_g=1$ 、そうでなければ $\mathcal{C}_g=0$ なので、先ほどのLCM Rush同様に添字の範囲を $i$ の約数に限定して $\mathcal{B}$ をメビウス変換すると $B\ast C$ が得られます。 $M \le 10^5$ ならば約数の個数 $d(M)$ は128を超えることがなく、また相異の素因数の個数 $p(M)$ は6を超えることがないため、 $O(Mp(M)d(M)\log d(M))$ でも2sに間に合います。