N以下の素数の個数

N以下の素数の個数π(N)をO(N^(3/4))で求めるアルゴリズムを説明する。
（真面目にやるとO(N^(2/3+ε))になったり、O(N^(2/3)/(logN)^2)になったりするらしいですが、それらはまだ理解していないので……）

○DPテーブルの設計
DP[i][n]=(n以下の数のうちi番目の素数まで篩にかけて生き残る数の個数)
とする。例えば
DP[0][10]=9、（2,3,4,5,6,7,8,9,10）
DP[1][10]=5、（2,3,5,7,9）
DP[2][10]=4、（2,3,5,7）
DP[3][10]=4、（2,3,5,7）
となる。
定義よりDP[π(√N)][N]=π(N)となることがわかるので、これを求めることを目標とする

○DPの遷移
i番目の素数をpi$p_i$で表すことにする。(1-indexed)
明らかに
DP[i][n]=DP[i-1][n]-(n以下の数のうち、pi$p_i$による篩で新たに落とされる数の個数)
となる。
「n以下の数のうち、pi$p_i$による篩で新たに落とされる数」
⇔「n以下の数のうち、最小素因数がpi$p_i$であるような数（pi$p_i$自身を除く）」
⇔「n以下の数のうち、pi$p_i$未満の素因数をもたず、pi$p_i$を素因数にもつ数（pi$p_i$自身を除く）」
であり、これはpi$p_i$で割ることで次と1対1対応する。
「n/pi$n/p_i$以下の数のうち、pi$p_i$未満の素因数をもたない数（1を除く）」
⇔「n/pi$n/p_i$以下の数のうち、pi−1$p_{i-1}$以下の素因数をもたない数（1を除く）」★
定義より
DP[i][n]=(n以下の数のうちi番目の素数まで篩にかけて生き残る数の個数)
=(n以下の数のうち、pi$p_i$以下の素因数をもたない数（1を除く）の個数)+(pi$p_i$以下の素数の個数)
だったので、★を満たす数の個数はDP[i-1][ n/pi$n/p_i$ ]-(i-1)となる。
よって
DP[i][n]=DP[i-1][n]-DP[i-1][ n/pi$n/p_i$ ]+(i-1)
という漸化式を得ることができた。

○計算の省略
DPテーブルの第2添字に入る値は1～NまでのN通りがあり得るが、漸化式をよく睨むとDP[*][N]を求めるための計算にはN/kの形で書ける数しか必要ないことがわかる
よって平方分割により実際には高々O(√N)個の部分しか計算しなくてよいことがわかる

○計算の省略2
DPの漸化式を見ると、DP[i][*]はDP[i-1][*]のみから計算できることがわかる。
また漸化式の右辺に登場するものの第2添字は、左辺に登場するものの第2添え字以下であるから、第2添字の大きいほうから順に計算することでDPテーブルを使いまわすことができる。（01ナップサックを思い出そう）
よって以下DP[n]の形で表すことにする。

○計算の省略3
pi$p_i$による篩の影響を計算する際には、k<p2i$k<p_i^2$であるような範囲についてはDP[k]の値は変化しない。このことを考慮することで計算量を削減することができる。具体的には次の通り
・pi<N1/4$p_i<N^{1/4}$の範囲
毎回O(√N)箇所のDPテーブルすべてを計算する。
N^(1/4)以下の素数の個数は自明にO(N^(1/4))なので全体で高々O(N^(3/4))
・N1/4<pi≤N1/2$N^{1/4}<p_i\le N^{1/2}$の範囲
更新が必要なのはk≥p2i$k\ge p_i^2$範囲だが、前述の平方分割テクによりこの範囲で計算する必要があるのはN/p2i$N/p_i^2$個
したがって全体では∑N1/4<pi≤N1/2Np2i${\displaystyle \sum _{N^{1/4}<p_i\le N^{1/2}}\frac{N}{p_i^2}}$箇所の更新が必要になる。
∑N1/4<pi≤N1/2Np2i<∑k=N1/4+1N1/2Nk2<∫N1/2N1/4Nx2dx=N3/4−N1/2=O(N3/4)${\displaystyle \sum _{N^{1/4}<p_i\le N^{1/2}}\frac{N}{p_i^2}<\sum _{k=N^{1/4}+1}^{N^{1/2}}\frac{N}{k^2}<{\int }_{N^{1/4}}^{N^{1/2}}\frac{N}{x^2}dx=N^{3/4}-N^{1/2}=O(N^{3/4})}$
なので全体で高々O(N^(3/4))

以上のDPは、あらかじめ√N以下の素数を列挙しておけば行えるため、時間計算量O(N^(3/4))、空間計算量(N^(1/2))で計算できることがわかった。

見てわかるとおり、この評価は極めておおざっぱなので正確に評価すると多分O(N^(3/4)/logN)とかくらいになる気がする。

出典はproject eulerのthread。ただし、この問題に正解していないとみることができない。
冒頭で紹介したMeissel–Lehmer algorithmの簡略版（のはず）だが、threadに書き込んだユーザーの名前を借り、project eulerではしばしば"Lucy's Algorithm"と呼ばれる。