ganmo::cout

競技プログラミング始めました

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

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

競プロで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):
熨斗袋さんからさらにパンチ力の高い解法を頂きました

問題の式は
\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に間に合います。

まとめ

包除がややこしくてわからなくても累積和の要領で殴れました。
いかがでしたか?かえってややこしくなってない?