n以下の素数のsumを高速に計算する

目次

動機

Project Euler で、n以下の素数の合計を求める問題がありました。 Problem 10 - Project Euler

この問題は2*10⁶以下の素数の合計を求めるだけなので簡単なのですが、Lucy_Hedgehog さんという方が面白いアルゴリズムをpostしていました。

せっかくなので、なるべく初心者にも分かりやすいようにアルゴリズムの概要を説明したいと思います。 アルゴリズムとか計算量とかどうでもいい人は、

optimized_prime_sum.cpp · GitHub

のプログラムをコピペして使って下さい。

問題概要

競技プログラミングで以下のような問題が出た場合、どう実装すればいいでしょう？

\(n\)以下の素数の合計を出力して下さい。ただし、n <= 10⁶

エラトステネスの篩を使えばn以下の素数の列挙は \( O(n \log\log (n)) \) で出来るので、このアルゴリズムで十分だと思います。

では次の制約の場合はどうでしょうか。

\(n\)以下の素数の合計を出力して下さい。ただし、n <= 10¹²

答えは 64bitに収まらないので、10⁹ + 7で割った値を答えて下さい。

アトキンの篩を使えば 計算 \( O(n / \log\log(n)) \) , 空間 \( O(n^{1 / 2}) \) で計算出来る(らしい)ので、これで十分かも知れません。 が、今回は計算 \( O(n^{3/4}) \), 空間\( O(n^{1 / 2}) \) で計算出来る方法を紹介したいと思います。

準備

\( S(n) \) = (\(n\)以下の素数のsum)　とします。 また、

\( S(n,k) = \) (\(n\)以下の数について、\(k\)(0-indexed)番目の素数までエラトステネスの篩を行った後に残っている数の合計)

と定義します。 分かりにくいので具体例を挙げると、 \[ S(20,0) = 2+3+5+7+9+11+13+15+17+19 = 101 \] \[ S(20,1) = 2+3+5+7+11+13+17+19 = 77 \] \[ S(50,2) = 2+3+5+7+11+13+17+19+23+29+31+37+41+43+47+49 = 377 \] となります。 0-indexedで考えると、0番目の素数は2, 1番目の素数は3ですね。 図を見れば分かりやすいと思うので、図も付けます*1。背景が灰色の数 or 背景の色が濃い数 がふるい落とされていない数ですね。

\(S(20,0)\)

\(S(20,1)\)

\(S(50,2)\)

篩は\(n^{1/2}\)の数までやれば良いので、 \(S(n,k) = S(n)　(n < p_{k+1}^{2}) \) が成り立つ事が分かります。ただし\(p_k\)は\(k\)番目の素数です。

一方 \( S(n,k)　( p_{k}^{2} < n) \) については以下の式が成り立ちます。

\[ S(n,k) = S(n,k - 1) - p_{k} * ( S(n / p_{k}, k - 1) - S(p_{k - 1}) ) \]

この式は

(\(p_k\)まで篩をした結果の合計) = (\(p_{k - 1}\)まで篩をした結果の合計) \(-\) (\(p_{k}\)自身を除く、\(p_{k}\) の倍数)

を表しています。 第2項の(\(p_{k}\)自身を除く、\(p_{k}\) の倍数)以外は定義と同じなので、第2項の説明をしたいと思います。

第2項を言い換えると\(p_{k} * ( \) \(n/p_{k}\)以下の数で 全ての素因数が\( p_{k} \)以上の数の合計\()\) となります。 一方で\(S(n / p_{k}, k - 1)\)は ( \(n/p_{k}\)以下の数で \( p_{k - 1} \)までの篩で残っている数の合計) です。

ここで、\( p_{k - 1} \)までの篩で残っている数とは、

• \(p_{k - 1}\)以下の素数
• \(p_{k - 1}\) より大きく、\(p_{k}\)以上の素因数のみで構成される数

の2種類です。後者がまさに求めたい数なので、\(S(n,k - 1)\)から \(p_{k - 1}\) 以下の素数の合計を引いた数が、まさに求めたかった数となりますね。 また、\(p_k \le n^{1/2} \) なので、\(S(p_{k - 1})\)は普通の篩で求めれば十分でしょう。

計算量

※数式が嫌いな人は無視するといいと思います。

(暗黙のうちに有理数を整数にcastしている箇所がいくつかあります。)

メモ化しつつ計算すればよいので、 \(S (a,b) の (a,b) \)の組み合わせの個数のオーダー = O(計算量) になりそうです。この個数を数えます。 \(b\)を固定した時, aの候補は\(\{n/1,n/2, \cdots ,n/z, \cdots, n/n\} ,\) ただし \(p_{k}^{2} \le a\) です。

まず\(a\)の候補ですが、少なくとも\(\{n/1,n/2, \cdots ,n/z, \cdots , n/n\}\)のどれかなので、\(p_{k}^{2} \le a\)を無視すれば \(O(n^{1/2})\)です。 これは、 \(n^{1/2} \le a\) となるような aの個数が高々\(n^{1/2}\) 個しか無いことから分かります。

\(1 \le p_b \le n^{1/4} \) の時

aの候補は、制約を無視しても\(O(n^{1/2})\)なので、 \((a,b)\)の組み合わせは高々\(O(n^{3/4})\)となります。

\(n^{1/4} \le p_b \le n^{1/2} \) の時

aは、\( a \in \{n/1,n/2, ... n/n\}, p_{b}^{2} \le a \le n \)です。\(n^{1/2} \le p_{b}^2 \le a\)なので、\(a\)の値は飛び飛びで計算しにくいので、逆数を取り式変形します。すると

\[ 1 \le a' \le \frac{n}{p_{b}^{2}} \]

ここで、\(a' = n/a = \frac{n}{n/z} \) であり、\(1 \le a' \le n^{1/2}\)なる\(a'\)となるような\(z\)が必ず1つ存在します。よって\((a,b)\)の個数は

\[ \sum_{n^{1/4} \le p_b \le n^{1/2}}{\frac{n}{p_{b}^{2}}} \]

です。これは積分すれば上限がわかり、 \(O(n^{3/4})\) です。

以上から、\((a,b)\)の組み合わせの個数は\(O(n^{3/4})\) となることが分かりました。 上手にプログラムを実装すれば計算量も\(O(n^{3/4})\)になりそうです。

素直なプログラム

上記をメモ化を使った素直なプログラムに直してみます。

メモ化でunordered_mapを使っているので結構メモリを食べます。残念ながら10¹²までのsumは計算出来ないでしょう。

最適化したプログラム

1つ前の結果しか参照しない事、以前のindexしか使っていないことからテーブルの使い回しが出来ます。また、a(=n)の取りうる値の特性から、unordered_mapではなくvectorで十分です。

計算 \( O(n^{3/4}) \), 空間\( O(n^{1 / 2}) \)のプログラムです。

Intel® Core™ i5-5200U @ 2.20GHz で約16sで終了します。

アルゴリズムの改良や高速化は大歓迎です。

余談

実は\(O(n^{2/3} + \epsilon)\)のアルゴリズムがあります。 Meissel-Lehmer algorithm と呼ばれるprime counting functionがあり、これは計算 \( O(n^{2/3 + \epsilon}) \), 空間\( O(n^{1 / 3 + \epsilon}) \) なのですが、このアルゴリズムを拡張することで、同様の(計算、空間)計算量でsumも計算することが出来ます。

気になった人は調べて見るといいと思います。

*1:こちらの画像の一部を切り出してきました。エラトステネスの篩 - Wikipedia