Official

G - Divisible by 3 Editorial by Nyaan

この問題は、乗法的関数の prefix sum を計算するアルゴリズムを利用することで解くことが出来る問題です。

乗法的関数とは？

正整数を引数に取る関数 \(f(n)\) が次の条件を満たす時、\(f(n)\) は 乗法的 であると呼びます。

  • 任意の \(\gcd(a, b) = 1\) である正整数 \(a, b\) に対して \(f(ab) = f(a) f(b)\) が成り立つ。

言い換えると、正整数 \(n\) が素数 \(p_1, p_2, \dots, p_k\) を用いて \(n = p_1^{e_1} p_2^{e_2} \dots p_k^{e_k}\) と表せる時、乗法的関数 \(f\) について \(f(n)\) が \(f(p_1^{e_1}), f(p_2^{e_2}), \dots, f(p_k^{e_k})\) の積として表せる、という解釈が出来ます。

乗法的関数に関する問題としてしばしば登場するものとして prefix sum の計算が挙げられます。すなわち次の問題です。

乗法的関数の prefix sum

乗法的関数 \(f(n)\) が与えられる。(素数 \(p\) に対して \(f(p^e)\) は \(\mathrm{O}(1)\) で計算できると考えてよい。) 次の式を計算せよ。

\[\sum_{1 \leq n \leq N} f(n)\]

今回の問題は、実は乗法的関数の prefix sum の計算に帰着させることが出来ます。

はじめにいくつかの事実を確認します。

事実 1

\(\sigma(n)\) を \(n\) の約数の総和を返す関数とする。\(\sigma(n)\) は乗法的関数である。

この事実は、\(n = p_1^{e_1} p_2^{e_2} \dots p_k^{e_k}\) における \(\sigma(n)\) が

\[\sigma(n) = \prod_{i=1}^k (1 + p_i + \dots + p_i^{e_i})\]

と表せるという公式から従います。

事実 2

正整数 \(M\) が与えられる。\(g(n)\) を総積が \(n\) である長さ \(M\) の正整数列の個数を返す関数とする。\(g(n)\) は乗法的関数である。

この事実は、いくらかの組合せ的解釈により \(n = p_1^{e_1} p_2^{e_2} \dots p_k^{e_k}\) において

\[g(n) = \prod_{i=1}^k \binom{e_i + M - 1}{M - 1}\]

になることを証明することで確認できます。(証明は練習問題とします, 青~黄 diff 程度)

以上より、今回の問題は次のように言い換えられます。

今回の問題の言い換え

乗法的関数 \(g(n)\) が与えられる。(素数 \(p\) に対して \(g(p^e)\) は \(\mathrm{O}(1)\) で計算できると考えてよい。) 次の式を計算せよ。

\[\sum_{1 \leq n \leq N, \sigma(n) \bmod 3 = 0} g(n)\]

この問題内の \(\sigma(n) \bmod 3 = 0\) の部分を取り除くと乗法的関数の prefix sum を計算する問題となります。この部分を上手く取り除いてみましょう。

乗法的関数 \(h(n)\) を素数 \(p\) に対して次式が成り立つように置きます。

\[ h(p^e) = \begin{cases} 0 & \sigma(p^e) \bmod 3 = 0\\ g(p^e) & \text{otherwise} \\ \end{cases} \]

すると、 \(\sigma(n)\) の乗法性により次式が成り立ちます。

\[ h(n) = \begin{cases} 0 & \sigma(n) \bmod 3 = 0 \\ g(n) & \sigma(n) \bmod 3 \neq 0 \\ \end{cases} \]

よって、今回の問題は次のように言い換えられます。

今回の問題の言い換え(2)

乗法的関数 \(g(n), h(n)\) が与えられる。(素数 \(p\) に対して \(g(p^e), h(p^e)\) は \(\mathrm{O}(1)\) で計算できると考えてよい。) 次の式を計算せよ。

\[\sum_{n=1}^N \left( g(n) - h(n) \right)\]

もちろん上式は

\[\sum_{n=1}^N g(n) - \sum_{n=1}^N h(n)\]

と変形できるため、\(g(n)\) の prefix sum から \(h(n)\) の prefix sum を引いたものを計算できればよいです。
以上より今回の問題を乗法的関数の prefix sum を計算する問題に帰着させることが出来ました。

乗法的関数の prefix sum の計算方法はいくつかありますが、ここでは

  • Lucy DP, および
  • black algorithm あるいは Min_25 篩

と呼ばれるアルゴリズムを利用した解法を説明します。

Lucy DP

以降では乗法的関数 \(f(n)\) の \(1 \leq n \leq N\) における prefix sum を考えます。また、関数 \(F(n), F_{\mathrm{prime}}(n)\) と整数集合 \(Q_n\) を次のように定義します。

\[F(n) = \sum_{1 \leq m \leq n} f(m)\]

\[F_{\mathrm{prime}}(n) = \sum_{1 \leq p \leq n, p : \mathrm{prime}} f(p)\]

\[Q_n = \left\lbrace m : m = \left \lfloor \frac{n}{k} \right \rfloor \text{を満たす正整数 } k \text{ が存在する} \right\rbrace\]

乗法的関数の prefix sum を 2 つのステップに分けて計算することにします。

  1. \(F_{\mathrm{prime}}(n)\) を \(n \in Q_N\) において全て計算する
  2. \(F_{\mathrm{prime}}(n)\) を用いて \(F(N)\) を計算する

1 番目のステップとして Lucy DP と呼ばれる DP が、\(2\) 番目のステップとして black algorithm あるいは Min_25 篩が登場します。

まずは Lucy DP の説明をします。Lucy DP は Project Euler: Problem 10 のスレッド(AC すると見られます) にある Lucy_Hedgehog 氏の投稿に由来する DP です。

簡単のため \(p\) が素数の場合に \(f(p) = 1\) が成り立つ場合を例にとって説明します。この時、\(F_{\mathrm{prime}}(n)\) は素数計数関数 \(\pi(n)\) と一致するので、ここでは \(\pi(n)\) の計算方法として説明します。

\(1 \leq x \leq \lfloor \sqrt{N} \rfloor,n \in Q_N\) の範囲において \(S(x, n)\) を次のように定義します。

  • \(S(x, n)\) : \(2\) 以上 \(n\) 以下の整数のうち、エラトステネスの篩のアルゴリズムにおいて \(x\) 以下の素数では篩われなかった整数の個数

\(\pi(n) = S(\lfloor \sqrt{N} \rfloor, n)\) が成り立つので、\(S(\lfloor \sqrt{N} \rfloor, \ast)\) を全て計算できればよいです。

DP の遷移として、\(S(x, n)\) を \(S(x-1, \ast)\) を用いて表すことを考えます。

  • \(x\) が素数でない場合

\(x\) において篩う処理が発生しないので \(S(x,n) = S(x-1,n)\) となります。

  • \(x\) が素数であり、\(n \lt x^2\) の場合

\(n\) 以下の \(x\) の倍数はすでに篩われています。よってこちらも \(S(x, n) = S(x-1, n)\) となります。

  • \(x\) が素数であり \(n \geq x^2\) の場合

\(x\) で篩われる数を \(mx\) と表すと、\(m\) のなす集合は「\(x-1\) 以下の素数で篩われなかった数の集合」から「\(x-1\) 以下の素数の集合」を除いたものになります。よって

\[S(x, n) = S(x-1,n) - \left(S(x-1, \lfloor n/x \rfloor) - \pi(x-1) \right)\]

であり、\(\pi(x-1) = S(x-1,x-1)\) に注意すると

\[S(x, n) = S(x-1,n) - S(x-1, \lfloor n/x \rfloor) + S(x-1, x-1)\]

が成り立ちます。以上より次式が成り立ちます。

\[ S(x, n) = \begin{cases} S(x-1, n) - S(x-1, \lfloor n/x \rfloor) + S(x-1, x-1) & (x \text{ is prime}) \wedge x^2 \leq n \\ S(x-1, n) & \mathrm{otherwise} \end{cases} \]

この DP は \(n\) のみを添え字に持って in-place に更新していくことができます。

  • Python による \(\pi(N)\) を計算する関数の実装例
import math
def Pi(N):
    sqrtN = math.isqrt(N)
    Q = [N // i for i in range(1, sqrtN + 1)]
    Q += list(range(Q[-1] - 1, 0, -1))
    S = {i: i - 1 for i in Q}
    for x in range(2, sqrtN + 1):
        if S[x] > S[x - 1]:
            for n in Q:
                if n < x * x: break
                S[n] -= S[n // x] - S[x - 1]
    return S[N]

DP の計算量を考えます。まず、ABC239Ex 解説 の観察 1, 2, 3 を認めるものとします。例えば \(|Q_N| = \mathrm{O}(\sqrt{N})\) です。
さて、計算量は S(v) -= S[v // p] - S[p - 1]の行が呼び出される回数を考えればよいです。\(x\) と \(N^{\frac{1}{4}}\) の大小で場合分けして考えます。

  • \(2 \leq x \leq N^{\frac{1}{4}}\) のとき

素数は素数定理より \(\mathrm{O}\left(\frac{N^{\frac{1}{4}}}{\log N}\right)\) 個あり、DP テーブルのサイズは \(\mathrm{O}(\sqrt{N})\) なので合計 \(\mathrm{O}\left(\frac{N^{\frac{3}{4}}}{\log N}\right)\) 回になります。

  • \(N^{\frac{1}{4}} \leq x \leq \sqrt{x}\) のとき

各 \(x\) について \(x^2 \leq n \leq N\) を満たす部分を更新しますが、\(\sqrt{N} \leq x^2\) より更新する個数は \(\mathrm{O}\left(\frac{N}{x^2}\right)\) となるのがわかります。よって素数定理と合わせると \(\mathrm{O}\left(\frac{1}{\log N} \sum_{x=N^{1/4}}^{\sqrt{N}} \frac{N}{x^2} \right) = \mathrm{O}\left(\frac{N^{\frac{3}{4}}}{\log N}\right)\) 回であることがわかります。

以上より全体の計算量は \(\mathrm{O}\left(\frac{N^{\frac{3}{4}}}{\log N}\right)\) になります。

ここでは \(\pi(n)\) を用いて説明しましたが、このアルゴリズムは他の乗法的関数における \(F_{\mathrm{prime}}(n)\) の列挙にも適用することが出来ます。
今回の問題で登場する \(g, h\) の素数 \(p\) における振る舞いは次の通りです。

\[g(p) = M\]

\[h(p) = \begin{cases} 0 & p \bmod 3 = 2 \\ M & \mathrm{otherwise} \end{cases} \]

です。\(g(p)\) の場合は定数なので \(\pi(n)\) を列挙した後に \(M\) 倍すれば所望のものを得られます。\(h(p)\) の場合は少し難しいですが、\(3\) で割って \(1\) 余る素数の個数と \(2\) 余る素数の個数を Lucy DP で求めることを考えれば所望のものを計算できます。(詳細は略します)

black algorithm

Lucy DP を用いて \(F_{\mathrm{prime}}(n)\) を \(n \in Q_N\) において列挙することができました。これらを用いて \(F(N)\) を計算するアルゴリズムである black algorithm の概要を説明します。

  • 考案者による記事 にも書かれているように、black algorithm は平易なアルゴリズムである点が特長です。

まず、以下の条件を満たす \(N\) 頂点の根付き木を考えます。ここで、\(n \geq 2\) に対して \(n\) の素因数のうち最大のものを \(\mathrm{gpf}(n)\) と表します。

  • 頂点に \(1\) から \(N\) までの番号がついている。
  • 頂点 \(1\) は根であり、頂点 \(n\) \((n \geq 2)\) の親は \(\frac{n}{\mathrm{gpf}(n)}\) である。

例えば \(N = 16\) における木は次の通りです。(葉の頂点を丸で、葉でない頂点を四角で囲っている。また、頂点 \(n\) \((n \geq 2)\) とその親を結ぶ辺を \(\mathrm{gpf}(n)\) と対応させた色で塗っている。)

image

この木の特徴として、頂点 \(n\) から頂点 \(1\) までのパスに \(n\) の素因数に対応する辺が素因数の大きい順に並んでいるという点があります。よって、現在の \(f(n)\) の値を持ちながらこの木の上を適切に DFS することで \(f(1), f(2), \dots, f(N)\) を \(\mathrm{O}(N)\) で列挙することが出来ます。そのため prefix sum も \(\mathrm{O}(N)\) で求められますが、これは低速です。

この DFS アルゴリズムを \(F_{\mathrm{prime}}(n)\) を用いて高速化してみましょう。今、葉でない頂点 \(n\) にいて、\(f(n)\) がわかっているとします。ここで、

  • \(k\) 番目に小さい素数を \(p_k\),
  • \(p_i = \mathrm{gpf}(n)\),
  • \(p_j\) を \(N/n\) 以下の最大の素数

とすると、\(n\) の子の頂点を小さい順に並べた列は \(n p_i, n p_{i+1}, \dots, n p_j\) のように表せます。そこで、\(f(n p_i), f(n p_{i+1}), \dots, f(n p_j)\) をまとめて足し上げることを考えます。\(n p_i\) は特殊処理するとして、\(p_{i+1}, \dots, p_j\) は全て \(n\) と互いに素なので、\(f\) の乗法性を利用すると

\[\sum_{k=i+1}^j f(n p_k) = f(n) \sum_{k=i+1}^j f(p_k) = f(n) \left( F_{\mathrm{prime}} \left( \left \lfloor \frac{N}{n} \right \rfloor \right) - F_{\mathrm{prime}}(p_i) \right)\]

となり、総和が \(F_{\mathrm{prime}}(n)\) を用いて高速に計算できます。この式を利用して、「ある頂点 \(n\) に移動するたびに、その頂点の子の \(f\) を足し上げる」という操作を行うことにすれば、DFS において葉の頂点に移動する必要が無くなります。(例えば \(N = 16\) の例だと、\(n=1,2,3,4,8\) 以外の頂点に訪問する必要が無くなります。)

以上の内容を適切に実装することで、\(F(N)\) を \(\mathrm{O}(根付き木の葉でない頂点の個数)\) で計算することが出来ます。具体的に計算量を式として表すと \(\mathrm{O}(N^{1 - \varepsilon})\) (\(\mathrm{O}(N)\) かつ、任意の \(\varepsilon \gt 0\) に対して \(\Omega(N^{1 - \varepsilon})\) が成り立つ、の意) になりますが、実際に葉でない頂点の個数を数えてみると \(N = 10^{10}\) で \(6298637 \simeq 6.3 \times 10^6\) 個と十分少ないため定数倍が非常に軽いことが確認できます。(black algorithm は競技プログラミングにおいて出題される \(N \leq 10^{12}\) 程度の制約においては十分高速に動作します。)

black algorithm は以上のように平明なアルゴリズムですが、計算量が \(\mathrm{O}(N^{1 - \varepsilon})\) である欠点もあります。そこで別解として Min_25 篩と呼ばれるアルゴリズムも紹介します。

Min_25 篩

Min_25 篩とは Min_25 氏によって考案されたアルゴリズムです。Min_25 篩は black algorithm と異なり \(\mathrm{O}(N^{\frac{2}{3}})\) という \(\mathrm{o}(N)\) の計算量評価が得られているアルゴリズムです。
また、black algorithm が \(F(N)\) のみを計算するのに対して、Min_25 篩は \(n \in Q_N\) に対して \(F(n)\) の計算結果が得られるという副産物もあります。
ただし、本来の Min_25 篩のアルゴリズムはかなり複雑なので、ここでは簡略化した \(\mathrm{O}\left(\frac{N^{\frac{3}{4}}}{\log N}\right)\) のアルゴリズムを説明します。

直感的に説明すると、Min_25 篩は Lucy DP の逆をするアルゴリズムです。
\(1 \leq x \leq \lfloor \sqrt{N} \rfloor,n \in Q_N\) の範囲において \(S(x, n)\) を次のように定義します。

  • \(S(x, n)\) : \(2\) 以上 \(n\) 以下の整数のうち、エラトステネスの篩のアルゴリズムにおいて \(x\) 以下の素数では篩われなかった整数に対する \(f(n)\) の総和

すると、明らかに \(S(\lfloor \sqrt{N} \rfloor, n) = F_{\mathrm{prime}}(n)\) であり、求めたいものは \(S(1, N) + 1\) です。よって、Lucy DP と逆方向に、すなわち \(S(\lfloor \sqrt{N} \rfloor, \ast)\) から始めて第 1 添え字が小さくなる方向に計算していくことを考えます。

Lucy DP と同様に \(S(x, n)\) と \(S(x-1, n)\) の差分を考えてみましょう。

  • \(x\) が素数でない場合

Lucy DP と同様 \(S(x, n) = S(x-1, n)\) です。

  • \(x\) が素数であり、かつ \(n \lt x^2\) である場合

Lucy DP と同様 \(S(x, n) = S(x-1, n)\) です。

  • \(x\) が素数であり、かつ \(n \geq x^2\) である場合

\(x\) において篩われる数は \(x^c m\) (\(m\) は \(x\) より大きい素数の積で表せる整数) という形をしています。\(c\) の値を固定して考えると、\(m\) としてあり得るものは

  • 『\(n/x^c\) 以下の「\(x\) 以下の素数で篩ったときに残った整数」』の集合から
  • 「\(x\) 以下の素数」の集合を取り除き、
  • さらに \(c \geq 2\) の場合は \(1\) を加えたもの

になります。この事実を元に立式して整理すると、結局次の式を手に入れられます。

\[S(x-1, n) = S(x, n) + \sum_{1 \leq c, x^{c+1} \leq n} f(x^c)\left( S(x,\lfloor n/x^c \rfloor) - F_{\mathrm{prime}}(x) \right) + f(x^{c+1})\]

上記の式を元に漸化式を計算すればよく、これは Lucy DP と同様に in-place な DP として高速化できます。計算量もまた Lucy DP と同様の方法で評価できて \(\mathrm{O}\left(\frac{N^{\frac{3}{4}}}{\log N}\right)\) になります。

また、このアルゴリズムおよび Lucy DP の両方を Fenwick tree を用いて高速化することで \(\mathrm{O}\left(N^{\frac{2}{3}}\right)\) に計算量を落とすことができて、高速化を施したアルゴリズムが本来の Min_25 篩です。(詳細は筆者も知りません…)

以上に説明したアルゴリズム、すなわち Lucy DP および (black algorithm あるいは Min_25 篩) を実装することで \(g(n), h(n)\) の prefix sum を計算することができるので今回の問題を解けます。計算量は \(\mathrm{O}\left(\frac{N^{\frac{3}{4}}}{\log N}\right)\) あるいは定数倍の軽い \(\mathrm{O}(N^{1 - \varepsilon})\) で十分高速です。

発展的話題: より計算量が良いアルゴリズム

本記事では乗法的関数の prefix sum を計算するアルゴリズムとして black algorithm と Min_25 篩を説明しましたが、乗法的関数の prefix sum を計算するアルゴリズムは他にも存在します。特に、中国の競技プログラミング界隈においては研究が盛んな分野で、様々なアルゴリズムが発明されています。
計算量の上では、 zhoukangyang 氏が 2023 年に発表した \(\mathrm{\tilde{O}}(\sqrt{N})\) のアルゴリズムが最速となっています。(ブログ, 論文) また、競技プログラミングでの実用上の最速は朱震霆氏の計算量 \(\Theta \left( \frac{N^{\frac{2}{3}} }{\left(\log N\right)^{\frac{4}{3}}}\right)\) のアルゴリズムであるようです。
また、乗法的関数が特定の条件を満たす場合は杜教篩や PowerfulNumber 篩と呼ばれるアルゴリズムを適用できます。
このように、乗法的関数の prefix sum は様々な人がアルゴリズムを提案しているので、興味がある方は調べてみると面白いと思います。

参考文献

posted:
last update: