ディリクレ双曲線法 (Dirichlet hyperbola method)はあまりなじみのない言葉かもしれませんが、競技プログラムでは有効なアルゴリズムだと思いますのでできるだけやさしく解説したいと思います。
ディリクレの畳み込み
ディリクレの畳み込みは以下の式で定義されます。いきなり難しいですが、これは単に演算の定義をしているだけです。
\boxed{(f*g)(n) = \sum\limits_{d|n} f(d) g\left(\frac{n}{d}\right) = \sum\limits_{ij=n} f(i) g(j) \ \ \cdots (1)}
一番単純な例として$f(n)$も$g(n)$も$n$の値に関係なく１を返す関数$1(n)$とします
\begin{align}
&n=12のとき \\
&ij=12は積が12になる全てのijのペアという意味なので \\
& (i,j) = [(1,12),(2,6),(3,4),(4,3),(6,2),(12,1)] \\ \\
& (f*g)(12) = \sum\limits_{ij=12} 1(i)\cdot 1(j)\\
& = 1(1)\cdot 1(12)+1(2)\cdot 1(6)+1(3)\cdot 1(4)+1(4)\cdot 1(3)+1(6)\cdot 1(2)+1(12)\cdot 1(1) \\
& = 1+1+1+1+1+1+1 = 6\\
& となり\\
& (f*g)(n) = (nの約数の数) = \sigma_0(n)\\
& という意味になります\\
& 従ってディリクレの畳み込みでは以下の関係式が成り立ちます\\
& \ \ \ \ (f*g)(n) = 1(n)*1(n) = \sigma_0(n)\\
& ディリクレ双曲線法では、求めたい関数に対してfとgのペアを見つけるのがポイントです\\
\end{align}
でも「それで何が面白いの？」と思った方はぜひ最後まで呼んでいただけたらと思います。あっと驚く結末が待ってますよ。この累積和を求める時に威力を発揮するのがディリクレ双曲線法です。
ディリクレ双曲線法
まず式（1）の累積和（１からｎまでの和）を取り、変形します。
\begin{align}
\sum_{t<n}(f*g)(n) &= \sum_{ \color{red} {t<n}}\sum_{ \color{red} {ij=t}} f(i) g(j) \\
& = \sum_{ \color{red} {ij \leq n}} f(i) g(j). \\
\end{align}
この式は例えば$n=15$のとき$xy = 15$の双曲線の下の格子点のすべての(i,j)に対して$f(i)g(j)$を求めて和を取ることを意味しています。
ここで下の図の赤で囲まれた$i=2$の部分の和を考えます。
\begin{align}
&\sum_{ j\leq n/2} f(2) g(j) \\
&=f(2)\sum_{ j\leq n/2}g(j) \\
&ここで\sum_{ j\leq n}g(j)が簡単に求まりG(n)とすると \\
&=f(2)G(n/2)\\
&となりシグマがなくなります。\\
\end{align}
\begin{align}
&したがって\color {blue}{青い点の部分}は \\
& f(1)G(n/1)+f(2)G(n/2)+f(3)G(n/3)\\
&= \color {blue}{\sum_{ i\leq \sqrt{n} } f(i)G(\lfloor n/i \rfloor) \cdots (2)} \\
&同様に\color {red}{赤い点の部分}は \\
&= \color {red}{\sum_{ j\leq \sqrt{n} } g(j)F(\lfloor n/j \rfloor) \cdots (3)} \\
&青と赤が重なっている\color {green}{緑の点の部分}は \\
&= \color {green}{ F(\lfloor \sqrt{n} \rfloor G(\lfloor \sqrt{n} \rfloor \cdots (4)}\\
& \\
&したがって求めたかった\\
&\sum_{ij \leq n}(f*g)(n) =(2)+(3)-(4) \\
& = \boxed{\sum_{ i\leq \sqrt{n} } f(i)G(\lfloor n/i \rfloor) + \sum_{ j\leq \sqrt{n} } g(j)F(\lfloor n/j \rfloor) - F(\lfloor \sqrt{n} \rfloor G(\lfloor \sqrt{n} \rfloor \cdots (5)} \\
\end{align}
これがディリクレ双曲線法の式です 一見複雑になっただけのように見えますがよく見ると シグマの上限が$\sqrt{n}$で抑えられている ので$O(\sqrt{n})$のアルゴリズムであることが分かります。
それでは具体的な例を使ってプログラムを書くところまで行きたいと思います。
求める関数f(i)とg(j)のディリクレ畳み込みで表す
最初に示した約数の数を求める関数$\sigma_0(n)=1(n)*1(n)$を使います
\begin{align}
&f(i)=g(j)=1(n) \\
&F(n)=G(n)= \sum_{i \le n}1(n) = n \\
\end{align}
ディリクレ双曲線法の式に代入
\begin{align}
&これを式(5)に代入します\\
&\sum_{ij \leq n}(f*g)(n) = \sum_{ i\leq \sqrt{n} } 1 \cdot \lfloor n/i \rfloor + \sum_{ j\leq \sqrt{n} } 1 \cdot \lfloor n/j \rfloor - \lfloor \sqrt{n} \rfloor \lfloor \sqrt{n} \rfloor \\
&= 2 \cdot \sum_{ i\leq \sqrt{n} } \lfloor n/i \rfloor - \lfloor \sqrt{n} \rfloor ^2\\
& (f*g)(n)はnの約数の数を求める関数\sigma_0(n) だったので\\
& その累積和をD(n)で表すと\\
&\boxed{D(n) = 2 \cdot \sum_{ i\leq \sqrt{n} } \lfloor n/i \rfloor - \lfloor \sqrt{n} \rfloor ^2 \cdots (6) }\\
\end{align}
これをプログラムにします。
式をプログラムにする
これは簡単ですね。sympyのdivisor_count関数と比較して答えが合っているか確認しました。
import sympy as sp
def D0n(n): # using divisor_count
return sum(sp.divisor_count(i) for i in range(n+1))
def D0h(n): # using Dirichlet hyperbola method
q = int(n**0.5)
return 2*sum(n//i for i in range(1,q+1)) - q**2
for n in [10, 100, 10**4]:
print(n, D0n(n), D0h(n))
# 10 27 27
# 100 482 482
# 10000 93668 93668
実行時間の比較
divisor_countとディリクレ双曲線法の比較をするとその差は歴然としてますね。
| n | $10^4$ | $10^6$ | $10^8$ | $10^{12}$ | $10^{16}$ |
|---|---|---|---|---|---|
| divisor_count | <1s | 15s | - | - | - |
| ディリクレ双曲線法 | <1s | <1s | <1s | <1s | 13s |
【応用編】以下の約数関数の累積和をディリクレ双曲線法を使って求めよ
\begin{align}
D_2(n) &= \large \sum_{i \le n} \sum_{d | n} d^2 \\
&= \large \sum_{i \le n} \sigma_2(n)
\end{align}
ポイントはペアとなる$f,g$を見つけることですが、それほど簡単ではないのでディリクレの畳み込み#具体例 (Wikipedia)を参照します。
\begin{align}
&\sigma_k = Id_k * 1 \\
& ここでId_k(n) = n^k\\
& ということなので\\
& f(n)=n^2, g=1(n)\\
& 級数の公式よりF(n)= \frac1{6}n(n+1)(2n+1), G(n)=n\\
& となります\\
\end{align}
これらを式(5)に代入します。
\begin{align}
&\sum_{ij \leq n}(f*g)(n) = \sum_{ i\leq \sqrt{n} } i^2 \cdot \lfloor n/i \rfloor + \sum_{ j\leq \sqrt{n} } F(\lfloor n/j \rfloor) - F(\lfloor \sqrt{n} \rfloor) \lfloor \sqrt{n} \rfloor \\
& = \sum_{ i\leq \sqrt{n} } (i^2 \cdot \lfloor n/i \rfloor + F(\lfloor n/i \rfloor)) - F(\lfloor \sqrt{n} \rfloor) \lfloor \sqrt{n} \rfloor\\
\end{align}
と求まりましたので、プログラムにします。
def sqsum(n):
return (n*(n+1)*(2*n+1))//6
def D2h(n):
q = int(n**0.5)
return sum((i**2)*(n//i) + sqsum(n//i) for i in range(1, q+1)) - sqsum(q) * q
for n in [10**4]:
print(n, D2h(n))
# 10000 400757638164
（開発環境：Google Colab）
この考え方はProject Euler、Problem 401: Sum of squares of divisorsを解くのに役に立ちます。