素数の個数や和を高速に計算できるLucy DPをやさしく説明

Project Euler, Problem 10: Summation of Primes のForumで投稿されていたLucy Hedgehog氏の解答（Lucy DPとも呼ばれる）がなかなか理解できなかったので、自分用のメモとして分かりやすさを最優先で説明しました。

素数の個数を求めるアルゴリズムの考え方

このアルゴリズムはMeissel–Lehmer algorithm (Wikipedia)の考えを元にしているということで、Wikiでは以下のように述べられています。

実際にすべての素数をリストすることなく$x$以下の素数の数を求める

一瞬そんなことができるのかと思いますが、これを一歩一歩説明して行きたいと思います

エラトステネスのふるい

プログラマであればそんなのは知っているよと言われるかもしれませんが、これをよ～く見るとこのアルゴリズムの肝が見えてきます。

例としてN=100の場合を以下に示します。

1. ふるいに使うのは$\sqrt{N}=10$以下の素数なので$[2,3,5,7]$
2. それらの素数の倍数を取り除いて言って残ったのが（#4: で白い四角の数）が素数
3. でも、ここではその素数ではなくて、ふるいで取り除かれた数に注目します。

ふるいで取り除かれた数

各回で、取り除かれた数の箱にはそれぞれ色が着いていますが、これを表にしてその数をそれぞれ$p$で割ります

p	取り除かれた数	pで割る	個数
$2$	$4,6,8,...98,100$	$\color{red}{2,3,4...,49,50}$	$49$
$3$	$9,15,21,...,93,99$	$\color{red}{3,5,7,...,31,33}$	$16$
$5$	$25,35,...,85,95$	$\color{red}{5,7,...,17,19}$	$6$
$7$	$49,77,91$	$\color{red}{7,11,13}$	$3$

これをよく見ると、

• 例えば$p=3$で取り除かれた$9,15,21,...,93,99$を3で割った$\color{red}{3,5,7,...,31,33}$は$p=2$で残った数の最初の方（図で赤字の数字）と一致している
• 同様に$p=5$で取り除かれた$25,35,...,85,95$を5で割った$\color{red}{5,7,...,17,19}$は$p=3$で残った数の最初の方（図で赤字の数字）と一致している

ここでその残った数の最初の方を表記するために、$S(v,p)$を$p$まで ふるいを行って残った数の$v$以下の個数（素数も含む）とします

\begin{align}
&S(100,1) = 99　(1を除く)　\\
&S(100,2) = 99 - 49 = 50 　(2を含む) \\
&S(33,2) = 17 　(2を含む) \\
&S(20,3) = 8  (2,3を含む) \\
&S(3,3) = 2 (2,3のみ、すなわち3以下の素数) \\
\end{align}

先程の素数$p$で取り除かれた数の個数を$R(p)$とすると

R(p)	S(v,p)で表す		個数
$R(2)$	$S(50,1)+0=S(100//2,1)+S(1,1)$		$49$
$R(3)$	$S(33,2)+1=S(100//3,2)+S(2,2)$		$16$
$R(5)$	$S(20,3)+2=S(100//5,3)+S(2,2)$		$6$
$R(p)$	$...=S(100//p,p-1)+S(p-1,p-1)$		$-$

となり$S(v,p)$は$S(v,p-1)$から$R(p)$を引けばよいので

\begin{align}
S(v,p) &= S(v,p-1) - R(p) \\
& = S(v,p-1) - (S(v//p, p-1)-S(p-1,p-1)) \cdots (1) \
\end{align}

という漸化式が得られました

S(v,p)の記憶領域

式$(1)$は単純な漸化式ではなくて右辺に$S(v//p, p-1)$があるので$v//p$の取りうる値すべてを記憶する必要があります
ただ少し試してみれば$N=100$のとき$v//p$は以下の19個しかないことが分かります

\begin{align}
& [100, 50, 33, 25, 20, 16, 14, 12, 11, 10, 9, 8, 7, 6, 5, 4, 3, 2, 1]\\
\end{align}

要は$1\le i \le \sqrt{N}$とその$i$で$N$を割った商$N//i$です。

S(v,p)がどう変わっていくか

アルゴリズムは$N=100$のときは$S(V)$の19個の領域を用意して、$(1)$の漸化式に従って更新すればよいことになります。これを図示すると、

矢印が$S(v//p, p-1)$の流れを示しています。これをよく見ると影響はvの小さいほうから大きい方にしかないので、Vを降順に更新していけば記憶領域は1行分あれば十分だと言うことが分かります。

この過程を終えて最後の行の$S(100,7)$が「100までの数に7までふるいをかけて残った数]=素数の個数25個となるわけです。

プログラミング

ここまでは結構複雑でしたがプログラムはシンプルです。ポイントは、

• 最初は$p$が素数かどうか分かっていないので$S(v,p)$の更新を行いながら$S[p]==S[p-1]$であれば素数でないので処理を行わない
• $S[v]$ の初期値は$v-1$（1を除くすべての数）
• $v<p^2$では$v//p < p$となり$S[v//p]-S[p-1]\le 0$となるので処理をする必要がない

N = 100
r = int(N**0.5)
V = [N//i for i in range(1,r+1)]
V += list(range(V[-1]-1,0,-1))
print(f"# V = {V}")
S = {i:i-1 for i in V}
print(f"# S = {S}")
for p in range(2,r+1):
    if S[p] == S[p-1]: continue # p is not prime
    sp = S[p-1]  # num. of primes smaller than p
    for v in V:
      if v < p**2: break  
      S[v] -= (S[v//p] - sp)  
    print(f"# p={p}, S(p-1)={sp}, S(v,p)={list(S.values())}")
print(f"# Number of primes less or equals to {N} is {S[N]}")
# V = [100, 50, 33, 25, 20, 16, 14, 12, 11, 10, 9, 8, 7, 6, 5, 4, 3, 2, 1]
# S = {100: 99, 50: 49, 33: 32, 25: 24, 20: 19, 16: 15, 14: 13, 12: 11, 11: 10, 10: 9, 9: 8, 8: 7, 7: 6, 6: 5, 5: 4, 4: 3, 3: 2, 2: 1, 1: 0}
# p=2, S(p-1)=0, S(v,p)=[50, 25, 17, 13, 10, 8, 7, 6, 6, 5, 5, 4, 4, 3, 3, 2, 2, 1, 0]
# p=3, S(p-1)=1, S(v,p)=[34, 18, 12, 10, 8, 6, 6, 5, 5, 4, 4, 4, 4, 3, 3, 2, 2, 1, 0]
# p=5, S(p-1)=2, S(v,p)=[28, 16, 11, 9, 8, 6, 6, 5, 5, 4, 4, 4, 4, 3, 3, 2, 2, 1, 0]
# p=7, S(p-1)=3, S(v,p)=[25, 15, 11, 9, 8, 6, 6, 5, 5, 4, 4, 4, 4, 3, 3, 2, 2, 1, 0]
# Number of primes less or equals to 100 is 25

本当に高速になっているのか

計算時間をsympyのprimerangeで素数を発生して数えた場合と比較しまが、その差は歴然としています。

N	$10^2$	$10^6$	$10^{8}$	$10^{10}$
primerange	0.01s	1.62s	192s	-
Lucy DP	0.01s	0.01s	0.42s	14.45s

素数の和を求める

個数を和に変えるのは驚くほど簡単で

• $S(v,p)$をふるいで取り除かれた数の和に変更する
• 式$(1)$を以下に変える

\begin{align}
&S(v,p) = S(v,p-1) - \color{red}{p \times}(S(v//p, p-1)-S(p-1,p-1)) \cdots (2) \
\end{align}

したがってプログラムは2行を以下のように変更するだけです

:
S = {i:i*(i+1)//2-1 for i in V}
:
      S[v] -= p*(S[v//p] - sp)
:

これはまさに【Project Euler】Problem 10: 素数の合計なので、こちらの投稿にプログラムと評価を追加しました。

【応用編】最小素因数の和

最小素因数とは以下のように$n$を素因数分解したときの最小の素因数のことです

n	2	3	4	5	6	7	8	9	10
最小素因数	2	3	2	5	2	7	2	3	2

上のアルゴリズムの素数$p$で取り除かれた数の個数を$R(p)$を使うと、

• 最小素因数が$p$になる個数は$R(p)+1$($p$自身)

となることが分かるので答えは$\sum_p^N p \times R(p)$と「素数の和」を足したものです

したがってプログラムは以下のように「素数の数」と「素数の和」を両方合わせたコードになります

N = 10**6
r = int(N**0.5)
V = [N//i for i in range(1,r+1)]
V += list(range(V[-1]-1,0,-1))
C, S = {i:i-1 for i in V}, {i:(i*(i+1)//2-1) for i in V}
lpsum = 0
for p in range(2,r+1):
    if S[p] == S[p-1]: continue # p is not prime
    cp, sp = C[p-1], S[p-1]     # sum,count of primes < p
    for v in V:
      if v < p**2: break
      S[v] -= (p*(S[v//p] - sp))
      C[v] -= (C[v//p] - cp)
      if v == N:
        lpsum += p*(C[v//p] - cp)

print(f"Answer: {lpsum+S[N]}")

この考え方はProject Euler、Problem 521: Smallest Prime Factorを解くのに役に立ちます

（開発環境：Google Colab）