Day 102-Q4 | Dirichletの双曲線法（約数総和関数）

問題

正の整数$N$が与えられる。$d(i)$を整数$i$の正の約数の個数とするとき、$D(N)=\sum_{i=1}^{N} d(i)$を求めよ。$N$は非常に大きい（最大$10^{12}$）ため、すべての$i$について$d(i)$を個別に計算してはいけない。

$D(N)$は約数ペア$(a,b)$で$ab\le N$を満たすものの個数と等しい。Dirichletの双曲線法を使うと、対称性$ab\le N \iff (a\le\sqrt N \text{かつ任意の}b)\text{または}(b\le\sqrt N\text{かつ任意の}a)$を利用して

$$D(N) = 2\sum_{a=1}^{\lfloor\sqrt N\rfloor}\left\lfloor\frac{N}{a}\right\rfloor - \lfloor\sqrt N\rfloor^2$$

という$O(\sqrt N)$の閉じた式に落とし込める。

入力形式

制約

$1 \le N \le 10^{12}$

入出力例

入力例1

出力例1

$d(1..10)=1,2,2,3,2,4,2,4,3,4$の総和が27。

入力例2

1000000000000

出力例2

27785452449086

$N=10^{12}$でも$r\approx10^6$回のループのみで即座に計算できる。

概念図

ヒント（段階的開示）

ヒント1: 方向性

$d(i)$を$i$ごとに素因数分解して求めるのは$N\le10^{12}$では到底間に合わない。「$D(N)$」を「$i$ごとの和」ではなく、「約数ペア$(a,b)$の個数」という別の見方に変換できないか考える。

ヒント2: アプローチ

$D(N)=\sum_{a=1}^{N}\lfloor N/a\rfloor$という形にまず変形する。この和は「$\lfloor N/a\rfloor$の値が変わらない区間ごとにまとめる」古典的な数論テクニック（商列挙）で$O(\sqrt N)$になるが、区間の右端の計算がやや煩雑。Dirichletの双曲線法を使うと、対称性から単純な1本のループだけで同じ答えが求まる、より簡潔な式に到達できる。

ヒント3: 誘導（コード骨格）

def D(N):
    r = int(N ** 0.5)
    while (r + 1) * (r + 1) <= N:
        r += 1
    while r * r > N:
        r -= 1
    total = 0
    for a in range(1, r + 1):
        total += N // a
    return 2 * total - r * r

$r=\lfloor\sqrt N\rfloor$を浮動小数点誤差なしに補正することと、`2*sum - r*r`という最終式さえ覚えていれば実装は数行で終わる。

模範解答 (Python)

import sys


def solve():
    N = int(sys.stdin.readline())

    r = int(N ** 0.5)
    while (r + 1) * (r + 1) <= N:
        r += 1
    while r * r > N:
        r -= 1

    total = 0
    for a in range(1, r + 1):
        total += N // a

    print(2 * total - r * r)


solve()

計算量: $O(\sqrt N)$（$N=10^{12}$でも$r\approx10^6$回のループのみ）。

Step-by-Step 解説

1「約数の個数の総和」を「格子点の個数」に翻訳
$D(N)$は、$1\le a,b$かつ$ab\le N$を満たす整数ペア$(a,b)$の総数と等しい。

2$\sqrt N$で領域を分割（双曲線法の核心）
$a\le\sqrt N$の範囲と$b\le\sqrt N$の範囲は正方形$a,b\le\sqrt N$の部分で重複するため、2倍してから正方形分を1回引く。

3整数演算で$r=\lfloor\sqrt N\rfloor$を正確に求める
浮動小数点誤差の可能性を`while`ループで補正し、必ず整数演算のみで確定させる。

4$O(\sqrt N)$ループで合計
`for a in range(1,r+1): total += N//a`で単純に足すだけ。最後に`2*total-r*r`が答え。

よくあるミス

ミス	原因	正しい書き方
`r=int(N**0.5)`をそのまま使い補正しない	浮動小数点誤差で$r$が実際の$\lfloor\sqrt N\rfloor$と1ずれることがある	補正ループで整数のみの比較により確定させる
$D(N)=\sum\lfloor N/a\rfloor$（片側だけ）を答えとして出力	双曲線法の「2倍して正方形を引く」変形をしていない	閉じた式$2\sum_{a=1}^{r}\lfloor N/a\rfloor - r^2$を使う
`r*r`をfloatにキャストしてオーバーフローを気にする	Pythonの整数は多倍長なのでオーバーフローの心配は不要	整数のまま`r*r`で計算すればよい
$N=1$などの境界値で式が壊れると誤解する	$N=1$なら$r=1$、`range(1,2)`で`a=1`は必ず回る	実際に$N=1$: $r=1$, total=1, 答え=$2\times1-1=1=d(1)$で正しい

次のステップ

発展: 約数総和$\sigma(i)=\sum_{d|i}d$の総和$\sum_{i=1}^N\sigma(i)$を同様に双曲線法で$O(\sqrt N)$で求める（$\sum_{a=1}^{r}(a\cdot\lfloor N/a\rfloor)$型の式に一般化）
発展: Min-25篩や杜教篩と組み合わせ、乗法的関数の前綴和を$O(N^{2/3})$でさらに高速化する
次回予告: Auction Algorithm（オークションアルゴリズム・価格上昇による分散型最適割当）

自己評価

完全理解おおむね理解要復習

Day 102-Q4 — Dirichletの双曲線法（約数総和関数）

問題