与えられた整数 $N$ に対し、 $N$ 以下の素数の数を数えるアルゴリズムを理解します。

Library Checkerで高速な実装を調べてみると、提出107115が最速ですが、このソースコードは未定義動作を含んでいます（ランキング上位を見てみると、62015, 107115, 147826, 150290,150294, 150330, 160040 はすべて同じバグを含んでいます。つまり、62015で埋め込まれたバグのようです）。上位陣の他の提出もほとんど同じアルゴリズムとなっています（Library Checkerへの最古の提出6832もsmallsへの最適化がされていないことを除いて実質的に同じアルゴリズムに見えます。どうやら昔からあるアルゴリズムのようです。smallsへの最適化は提出10100で出現しています）。

このアルゴリズムの詳しい解説が以下にありますが、もう少し行間を埋めていきたいと思います。

何をしているのか

参考：素数の数え上げと乗法的関数の和

上記の非常に詳しいスライドを見ればわかりますが、素数計数関数は以下の再帰関数で計算可能です。

$S(1, \nu) = \nu - 1$
$S(i, \nu) = S(i-1, \nu) - S(i-1, \lfloor\nu/i\rfloor) + \pi(i-1) \text{ if } i \text{ is prime and } i^2 \le \nu$
$S(i, \nu) = S(i-1, \nu) \text{ if } i \text{ is not prime or } i^2 \gt \nu$
$\pi(\nu) = S(\nu, \nu)$

この再帰関数を配列を使いまわす動的計画法（DP）で計算するというアルゴリズムになっていますが、いくつか工夫がほどこされています。

基本的なアルゴリズム

$S(i, \nu)$ の意味

$S(i, \nu)$ は、2以上 $\nu$ 以下の自然数のうち、エラトステネスの篩で $i$ 以下の素数までを使って篩ったとき、残っているものの数です。残っているものとは、 $i$ までの素数と、最小素因数が $i$ を超えるもの（(i+1)-rough number）の数です。 ※後での説明の都合から、X-rough numberには1も含む定義を採用します（1を含まない定義もよく使われます）。それにかかわらず、 $S(i, \nu)$ の意味するところは「2以上 $\nu$ 以下」なので1はカウント対象外です。

例： $S(4, 100) = 34$ です。2以上100以下の自然数のうち4以下の素数までを使って篩ったとき残っているものは、以下の34個です。

素数（25個。2, 3も残っていることに注意）
合成数だが最小素因数が5以上のもの（25, 35, 49, 55, 65, 77, 85, 91, 95の9個）

再帰の様子

$S(31, 1000)$ を計算するときの再帰の様子を示しました。まず、 $S(i, \nu)$ は、 $i$ 以下の最大の素数を $p$ として $S(p, \nu)$ と同じなので、縦方向では素数のみを示しました。また、横方向では、 $\lfloor 1000/j\rfloor$ としてあり得る値のみを示しました。水色の部分は、 $i^2 \ge \nu$ であるため、 $\sqrt i$ 以下の最大の素数を $p$ として $S(p, \nu)$ と同じになる部分です。

第二引数 $\nu$ としてあり得る値

再帰的に呼び出された場合も含めた場合、第二引数としてあり得る値は自明ではありません。例えば、 $j_1, j_2$ を自然数としたとき、 $\lfloor\lfloor n / j_1\rfloor / j_2\rfloor$ はどのような値を取りえるでしょうか？

実は、 $\lfloor\lfloor n / j_1\rfloor / j_2\rfloor$ は $\lfloor n / j_1 / j_2\rfloor$ に等しいことが示せます。その理由は、次のように説明できます。

まず、 $\lfloor n / j_1\rfloor$ は $n / j_1$ と比べて0以上 $(j_1-1)/j_1$ 以下しか小さくありません。同様に考えれば、 $\lfloor \lfloor n / j_1\rfloor / j_2\rfloor$ は $\lfloor n / j_1\rfloor / j_2$ と比べて0以上 $(j_2-1)/j_2$ 以下しか小さくないとわかります。したがって、 $\lfloor \lfloor n / j_1\rfloor / j_2\rfloor$ は $n / j_1 / j_2$ よりも0以上 $(j_2-1)/j_2 + (j_1-1)/j_1/j_2 = 1 - 1/j_1/j_2$ 以下しか小さくありません。 $n / j_1 / j_2$ よりも0以上1未満しか小さくない非負整数は、唯一 $\lfloor n / j_1 / j_2\rfloor$ だけです。したがって、 $\lfloor \lfloor n / j_1\rfloor / j_2\rfloor$ は $\lfloor n / j_1 / j_2\rfloor$ に等しくなります。

余談：符号付きの場合であっても、C言語の除算（零への丸め）ならば符号と絶対値を独立に考えられるので同様の結果が成り立ちます。

よって、 $\lfloor\lfloor n / j_1\rfloor / j_2\rfloor$ がとりえる値の集合は $\lfloor n / j_1\rfloor$ の集合と同じです。この結果を再帰的に使えば、 $\lfloor\lfloor\lfloor n / j_1\rfloor / j_2\rfloor / j_3\rfloor$ がとりえる値の集合も $\lfloor n / j_1\rfloor$ の集合と同じということがわかります。

つまり、 $\nu$ としてとりえる値は、 $\lfloor N / j_1\rfloor$ がとりえる値です。結局、 $\nu$ としてとりえる値の集合は $Q = \{ 1, 2, 3, \dots, \lfloor\sqrt N\rfloor, N/\lfloor\sqrt N\rfloor, \dots, N/3, N/2, N/1\}$ となります（※ $\lfloor\sqrt N\rfloor$ と $N/\lfloor\sqrt N\rfloor$ は同じである可能性があります）。

DPで求める

これを動的計画法で求める方法を考えます。非常に単純には、 $\lfloor\sqrt N\rfloor$ までの素数 $p$ （縦方向）と $q \in Q$ （横方向）のすべての組み合わせについて、配列を順番に埋めていけばよいですが、これではメモリを使いすぎです。幸い、 $S(p, q_1)$ の計算には、 $p$ の直前の素数を $prev\_p$ として、 $S(prev\_p, q_2), q_2 \in Q, q_2 \ge q_1$ の結果のみが必要なことが利用できます。つまり、直前の結果しか必要ないので、配列を使いまわすことができます。また、 $q_1$ の大きい方から順に更新していけば、（表裏の二つの配列を行ったり来たりして使いまわす方式ではなく）単一の配列を上書き更新することができます。

つまり、以下のように実装できます。

function prime_count( N )
  v = int(sqrt(N))
  S[q] = q - 1 for q ∈ Q # init

  for c = 0, 1, 2, ... do
    p = (c+1)-th prime # 2, 3, 5, ...
    for q ∈ Q, descending order
      if q >= p * p then
        S[q] -= S[q/p] - c
      endif
    endfor
  endfor
  return S[N]
endfunction

なお、S[q]を配列でとってしまうと大きすぎです（ $\Omega(N)$ のメモリ要求が発生してしまいます）。実際に使う部分は高々2 * v個しかないので、うまく圧縮したいです。連想配列を使ってもいいですが、高速化のためには配列で頑張りたいです。そこで、 $q \in \{ 1, 2, 3, \dots, \lfloor\sqrt N\rfloor\}$ の場合を格納するSmall[v+1]と、 $q \in \{N/\lfloor\sqrt N\rfloor, \dots, N/3, N/2, N/1\}$ の場合を格納するLarge[v+1]を用意します。 Smallの方では、インデクスはqをそのまま使えば良さそうです。 Largeの方では、インデクスはN/qを使えば良さそうです。このようにすることで、必要なメモリ量を削減することができます。

このアルゴリズムで最も時間のかかる工程はS[q]を更新する手順で、これは $\Theta(N^\frac34 / \log N)$ 回発生します。最も時間を要するのは、 $p \sim N^\frac14$ 付近のループです。 $p \le N^\frac14$ の範囲に素数は $\Theta(N^\frac14 / \log N)$ 個存在しますが、それらについて毎回少なくとも、大きいqのS[q]（Large[N/q]に相当する部分。v個ある）の更新が必要、という部分が計算量の大半を占めています。それより大きい部分では、途中のif文でスキップされる割合が急速に増加するため、 $p$ としてとりえる数の候補が多いにもかかわらず、計算量を支配しません。特に、 $p \ge N^\frac13$ の範囲では、S[N]だけが更新されます。

三種類の更新

上の疑似コードでS[q] -= S[q/p] - cと単純に書いていたところは、インデクス変換により、以下の三種類の更新に場合分けすることができます。

Large ← Large の更新

// j ∈ {1}∪(p, v/p] if p^4 < N
// j ∈ {1}           if p^4 >= N

// S(p, N/j) ← S(p, N/j/p)
Large[j] -= Large[j * p] - c

Large ← Small の更新

// j ∈ (v/p, v]   if p^4 < N
// j ∈ (p, N/p/p] if p^4 >= N

// S(p,N/j) ← S(p,N/j/p)
Large[j] -= Small[N/j/p] - c

Small ← Small の更新

// j ∈ [p*p, v) if p^4 < N
// j ∈ ∅        if p^4 >= N

// S(p, j) ← S(p, j/p)
Small[j] ← Small[j/p] - c

`roughs`を使うアルゴリズム

ここまでは再帰関数を素直にDPに変換したものであり、それほど難しくはありません。上記疑似コードの例も10行程度と非常にわかりやすいコードです（実際にはインデクス変換の計算が面倒ですが……）。しかし、あのアルゴリズムに出てきたroughsがまだ出てきていません。

ここからは、本題となるroughsを使うアルゴリズムの動作を説明します。 roughsは(p+1)-rough numberを管理している配列ですが、なぜそれを管理しているのかを説明していきます。

`roughs`が出てくる理由

先ほどの図を見れば、 $S(31, 1000)$ を計算する際に、 $S(p, \nu)$ のすべてが使われるわけではないということがわかります。例えば、 $S(p, 1000/4), S(p, 1000/8), S(p, 1000/9), S(p, 1000/12), S(p, 1000/16), S(p, 1000/18), ...$ などは使われていません。 $S(p, 1000/6)$ は使われているため、合成数 $j$ について $S(p, 1000/j)$ が使われない、というわけでもなさそうです。これをどう考えたらよいでしょうか。

DPの更新手順が、S[q] -= S[q/p] - cとなっていることに注目します。つまり、参照されるのはq/pとしてあり得る数だけです。 qをN/jと表せば、q/pはN/j/pとして表せます。逆に言うと、j*pと表せないkについてはS[N/k]の形で読み出されることがありません。よって、「(p+1)-rough number k について、S[N/k]の形では読まれることがない」と結論付けることができます。

このことを生かして更新回数を削減しているのがあのアルゴリズムです。 roughsは、v以下の(p+1)-rough numberを管理しており、sはその個数を示しています。今回のpで割り切れる、つまり(p+1)-rough numberでなくなる数について、std::erase_ifのようなことをやって前に詰めています。 roughsだけでなく、largesも詰める必要があるため、更新手順が-=を使ったものでなくなっています。この詰める手順も、qの大きい方からやっていくというのと合致していて、単一の配列を上書き更新していくことができます。

細かい注意ですが、N/kがN/k'と一致している場合、k'が(p+1)-rough numberでなければ、k'経由でS[N/k']を読むことはあります。 k>vの場合はこのようなk'が存在する可能性があり、それがゆえにsmallsの方はこういう更新回数削減が使えません。

配列のインデクス変換の効率よい実装

インデクス変換のアイデアは単純ですが、効率の良いコードに落とし込むのはなかなか工夫が必要です。あのコードが何をやっているかを詳しく見てみます。

`smalls`の方

smallsのインデクスはなんだか右シフトされている感じがあります。これは、2以上の $p$ について $S(p, 2i) = S(p, 2i-1)$ が成り立つ（偶数は(p+1)-rough numberにならない）ことを利用した圧縮です。引数が奇数とわかっていればそのまま右シフトすればよく、偶奇がわからない場合は1引いてから右シフトしています。

`larges`の方

smalls[d >> 1] - pcがインデクスになっていますが、かなり謎です。これは解説（眠れない夜は素数の個数でも数えましょう - えびちゃんの日記）を読まないとわかりませんでした。これはどうやらdがroughsの中で何番目に位置するかを計算しているようです。

読み出すデータは全て更新前のデータ、つまりprev_pでの値であることに注意して、以下のように考えればそれがわかります。まず、roughsには(prev_p+1)-rough numberが入っています。次に、S[d]は2以上d以下の数で、 $prev\_p$ までの素数と最小素因数が $prev\_p$ を超えるもの（(prev_p+1)-rough number）の数が入っているはずです。また、pcは、prev_pが2の時に0であることからわかるように、 $prev\_p$ までの素数の個数より一つ少ないです。一方、実はsmall[idx]はS[2*idx+1] - 1を管理しているので、一つ少ない数が返ってきます。 dは奇数なので、small[d>>1]がS[d]-1を表しています。よって、この二つの-1が打ち消しあって、small[d>>1] - pcが、2以上d以下の(prev_p+1)-rough numberの数です。つまり、dはsmall[d>>1] - pc番目の、2以上d以下の(prev_p+1)-rough numberです。 0-indexとの差分は、roughs[0]に1が居座っていることで打ち消され、結局dはroughs[small[d>>1] - pc]にあることがわかりました。

計算例： $p = 5$ の周回で、 $d = 5\times11 = 55$ を使う場合を考えます。まず5-rough numberを列挙すると、1, 5, 7, 11, 13, 17, 19, 23, 25, 29, 31, 35, 37, 41, 43, 47, 49, 53, 55, ...となります（要するに6で割って1余るか5余る数です）。 55が何番目かを考えると、0-indexで18番目であることがわかります。この18を求めているのがsmall[d>>1] - pcの部分です。 $S(3, 55) = 20$ であり、small[d>>1]は19を返します。 pcは1である（0で始まり、p = 3のループで一回だけインクリメントされた）ので、引き算すると18が得られます。

$p > N^\frac14$ の計算簡略化

あのアルゴリズムを見ると、DPをやり終えてからreturn S[N]とするのではなく、なんだか最後の方でよくわからない計算をしています。そもそもよく見ると、DPで更新すべき $p$ の範囲は $N^\frac12$ までだったはずが、 $N^\frac14$ までで脱出してしまっています。これは、先ほど説明した「三種類の更新」のうち、さらにp^4 < Nかで範囲が変わることに対応しています。

三種類の更新

Large ← Large の更新

Large[1] -= Large[p] - c

だけが実行されます。 Large[p]を更新するチャンスはまだありますが、その寄与は後で考えます。また、roughsには $N^\frac14$ を超えて $N^\frac12$ 以下の素数が全て格納されています。よって、

ret = larges[0]; // Large[1]
for( k = 1; k < s; ++k ) {
    ret += larges[k] - ( pc + k - 1 ); // Large[p] 
}

とやればこの部分の計算を一気に片付けることができます。ここのくくりだしは、特に計算量のオーダーを変えているわけではなさそうですが、シンプルになっていて高速に動作するようです。

なお、あの実装では... += (s + 2 * (pc - 1)) * (s - 1) / 2;などとしていますが、これは上記コードの( pc + k - 1 )部分の寄与を誘導変数最適化したものです。これにより、ループ内が極限までシンプルになり、高速に動作するようです。

Large ← Small の更新

後回しにした以下の部分の寄与を考える必要があります。

// j ∈ (p, N/p/p]
// S(p,N/j) ← S(p,N/j/p)

Large[j] -= Small[N/j/p] - c

Large[j]はこの分だけ小さくなるはずなのに、上ではret -= Large[p] - cとしてしまっていました。つまり、引きすぎになっているので、その分だけ足し戻します。

特定のpの時のことを考えます。この時、いろいろなLarge[j]に対して「引きすぎ」の寄与が一回ずつ発生しますが、retへの寄与はどのj経由だろうと同じなので、どのj経由で寄与が発生したかはあまり興味がありません（だからと言って計算量が減るわけではなさそうですが）。重要なのは、pの次の素数からjとして選ばれるためそこから寄与が始まること、N/p/pまでで止めないといけないことです。 N/p/pが（std::lower_bound的な意味で）roughsの中のどの位置にいるかは、先ほど話したテクニックと同じで、smalls[(N/p/p-1) >> 1] - pcで求めることができます。 smallsはもう変化せず、またroughsはもうずらしたりしないので、pcは今の値を使えばよいことに注意します。

これを各pに対して行えばいいので、

for( k = 1; k < s; ++k ) {
  p = roughs[k];
  k2_max = smalls[(N/p/p-1)>>1] - pc
  for( k2 = k + 1; k2 <= k2_max; ++k2 ) {
    j = roughs[k2];
    ret += smalls[(N/p/j-1)>>1] - ( pc + k - 1 ); // Small[N/j/p];
  }
}

とやればよいです。あの実装で(e - l) * (pc + l - 1)などとしているのは、例によって誘導変数最適化です。ここのループはかなりの回数実行されるため、この追い出しは重要です。

Small ← Small の更新

何もやる必要がありません。

計算量の解析

あの実装の時間計算量は、軽い計算（除算以外の整数演算とメモリアクセス）が $\Theta(N^\frac34/\log N)$ 回、重い計算（除算）が $\Theta(N^\frac34/(\log N)^2)$ 回、の合計です。さすがに除算ともなると軽い計算の $\Theta(\log N)$ 倍の計算コストがかかるとみなすのが妥当で、その仮定の下で $\Theta(N^\frac34\log N)$ にちょうど収まっています。

余談：除算のコスト

普通のCPUではnビットの除算命令は $\Theta(n)$ 時間かけて行われます。基本的には筆算みたいな感じで上から商を決めていきます。最近だと1サイクルに6ビットくらい商を決められるらしいです（Radix-64 SRT dividerでググってみましょう）。これは倍精度浮動小数点数（精度53bit）の除算を十数サイクルでできることと整合性があります（商を6ビットずつ立てると9サイクルかかって、前処理とか後処理とかでもう少し時間がかかる）。

ニュートン法を使えば $\Theta(\log n)$ 時間でできそうに見えますが、倍精度(n=53)くらいでは実用的ではありません。ニュートン法よりも並列度の高いGoldschmidt法がIBMのメインフレームで使われていたらしいですが、これでもまだ実用的ではないそうです（参考：コンピュータアーキテクチャの話(103) Newton-Raphson法とGoldschmidt法(2) | TECH+（テックプラス）。四倍精度ならあるいは、との指摘も書かれていました）。これは、この手の方式は反復のたびに乗算が必要であり、1サイクルで反復を終わらせるのは無理であることが大きいです。仮にt=2サイクルかかるとすれば、7bitの精度で初期商を推定してから3反復で53bit精度に達しますが、正しく丸めるためにはこの精度では足りませんので、もう1反復必要で、合計で4反復必要です。上記記事の図27にあるように、反復回数がi回ならit+1サイクルかかるので、これでもう9サイクルです。巨大な乗算器回路を使い続けたのにもかかわらずRadix-64 SRT divider（これはかなり小さい）と同速度では、実用的ではありません。でもSIMD命令だとSRT除算器を作るよりは既存の乗算器を使いまわしたほうがいいという判断でGoldschmidt法が使われているかもしれません（IntelのCPUでdivsd命令は13サイクルなのにdivss命令は11サイクルもかかっているので、SRT除算器じゃなさそう？）。

整数除算が100サイクル近くかかるのは謎で、これはIntelが最適化をサボっていただけだと思われます。最近のIntelのプロセッサ（Sunny Cove以降？）は128bit整数÷64bit整数（のうち商が64bitに収まるもの）を18サイクルで完了させるそうです。これはきっとRadix-16 SRT dividerを使って1サイクルに4ビットずつ商を決めているのでしょう（ポート割り当てが3*p1なので、前処理と後処理がありそうです）。

参考：素数カウント( $\mathrm{O}(\frac{N^{\frac{3}{4}}}{\log N})$・高速化版) | Nyaan’s Library

高速化をしたい場合は整数同士の除算をdouble型同士の除算で計算するのが最も高速化への寄与が大きい。なぜなら、64bit整数同士の除算は最悪ケースで80クロックと低速なのに対してdouble型同士は11クロックとはるかに高速であり、またN=10¹²程度の大きさ同士の切り捨て除算ならばdouble型で計算しても誤差が発生しないことが規格で保証されているためである。

倍精度浮動小数点数を使うと速くなるのは実際かなりそうなんですが、それでも依然として $\Theta(\log N)$ のコストがかかっていると考えたほうが良いです。いまどきのCPUでは、加算命令やロードストア命令は（スループットとして）1サイクルに複数回できるのに対して、除算命令は（スループットとして）十数サイクルに一回しかできず、何十倍もの差があります。また、短いビット幅の除算命令の方が有意に早く終わることなども、除算命令のコストが $\Theta(\log N)$ あることを示唆しています。実際、この実装では倍精度浮動小数点数除算を使っているにもかかわらず、それでもなお除算がボトルネックになっているようです。

詳細な解析

初期化

配列の初期化には、 $\Theta(N^\frac12)$ 回の除算が必要ですが、他の所に比べたら計算量は全然少ないです。

エラトステネスの篩部分

エラトステネスの篩を $N^\frac12$ までやっているので、 $\Theta(N^\frac12 \log\log N)$ 回のメモリアクセスが発生していますが、これも他の所に比べたら計算量は全然少ないです。

Large の更新部分

$N^\frac14$ までの素数 $p$ （ $\Theta(N^\frac14 / \log N)$ 個ある）について、 $N^\frac12$ 未満の(p+1)-rough numberは $\Theta(N^\frac12 / \log p)$ 個（あってる？）あり、その組み合わせの通り数分だけの更新が発生します。積分を計算すると、更新の回数は $\Theta(N^\frac34 / (\log N)^2)$ 回のようです。

Large ← Large の更新では除算が不要ですが、Large ← Small の更新では除算が必要です。計算量への寄与は、 $N^\frac14$ に近い素数での更新が大半を占めていて、その場合はほとんど Large ← Small の更新なので、ほぼ毎回除算が発生すると考えてよいです。

よって、計算コストは $\Theta(N^\frac34 / \log N)$ 程度です。

なお、係数は $\frac12$ と $e^{-\gamma} \sim 0.56$ の間だと思います。これはBuchstab関数 $\omega(u)$ の $u=2, \infty$ での値です。定性的には、次のように考えればこの係数を導くことができます。エラトステネスの篩で $p$ まで篩った残りが(p+1)-rough numberです。この時、 $p$ から $p^2$ までの(p+1)-rough numberは素数です。つまり、 $p$ から $p^2$ の区間では、平均素数密度と(p+1)-rough numberの密度は一致します。よって、 $p^2$ 以下の(p+1)-rough numberの密度は、 $1/\log p^2 = \frac12 / \log p$ 程度となります。

$p^2$ を超えた部分では、(p+1)-rough numberの密度はこれとは異なります。 (p+1)-rough numberの出現位置は、 $p\#$ （pの素数階乗）という周期を持ちます。その周期の中には、 $\phi(p\#)$ 個の(p+1)-rough numberが存在します（ $\phi$ はオイラーのトーシェント関数です）。その比は漸近的に $e^{-\gamma} / \log p$ となることが知られているようです（Mertens' third theoremというらしいです。https://arxiv.org/pdf/2308.07570.pdfから見つけました）。よって、(p+1)-rough numberが大きく偏っていないとすれば、 $p^2$ を超えたあたりでの密度は $e^{-\gamma} / \log p$ くらいということになります。

結局のところ、 $p$ としては $N^\frac14$ 付近のものが多いので、係数はおそらく $\frac12$ になると思います。ただ、 $N^\frac15$ くらいまでは係数は $e^{-\gamma}$ であり、 $N^\frac15$ までの素数は $N=10^{11}$ において37個と、 $N^\frac14$ までの素数102個の1/3以上を占めるので、実際には $\frac12$ より少し大きめと思ったほうが良さそうです。

ところで、 $p = 13$ までの場合、 $\phi(p\#)/p\# \times \log p$ は $\frac12$ よりも小さいです。小さい $p$ で実験していて、 $p$ から $p^2$ の区間こそが(p+1)-rough numberの密度が高い地点になっていて合わないなぁと悩んでいました。あくまで漸近挙動なので小さい $p$ では $\frac12$ との大小関係という定性的な性質ですらも成り立たないようです。

Large ← Large の更新

$N^\frac14$ までの素数 $p$ （ $\Theta(N^\frac14 / \log N)$ 個ある）について、 $p$ から $N^\frac12/p$ までの(p+1)-rough numberが何個あるかを求めます。まず、 $N^\frac16$ までの素数 $p$ については、上限である $N^\frac12/p$ は $p^2$ より大きいので、 $\Theta(N^\frac12 / p / \log p)$ 個あるとして良さそうです（あってる？）。それ以上の素数 $p$ に関しては、上限である $N^\frac12 / p$ は $p^2$ より小さいので、(p+1)-rough numberというのは全て素数です。したがって、 $\Theta(N^\frac12 / p / \log(N^\frac12 / p))$ 個くらいあります。

この組み合わせの通り数分だけの更新が発生します。積分を計算すると、 $N^\frac16$ までの素数 $p$ での更新が $\Theta(N^\frac12 / \log N)$ 回、そこから $N^\frac14$ までの素数 $p$ での更新が $\Theta(N^\frac12 / \log N)$ 回、それぞれ発生するようです。よって、合計の更新回数は $\Theta(N^\frac12 / \log N)$ 回のようです。

係数は、前半部分が $\omega(2) = \frac12$ くらい、後半部分が $\log\log(N^\frac13) - \log\log(N^\frac16) = \log2$ でしょうか？

Large ← Small の更新

Large ← Large の更新は、オーダー的にとるに足りないことを上で示したため、Largeの更新回数とオーダーは一致します。つまり、更新の回数は $\Theta(N^\frac34 / (\log N)^2)$ 回のようです。

Small の更新部分

$N^\frac14$ までの素数 $p$ （ $\Theta(N^\frac14 / \log N)$ 個ある）について、 $p^2$ から $N^\frac12$ までの自然数は $N^\frac12 - p^2$ 個あり、その組み合わせの通り数分だけの更新が発生します。積分を計算すると、更新の回数は $\Theta(N^\frac34 / \log N)$ 回のようです。 $N^\frac12$ まで計算しないといけない、という部分と、 $p^2$ までは更新しなくてよい、という部分は、いずれも同じオーダーです。前者の係数は $1$ であり、後者の係数は $\frac13$ です。よって、全体の係数は $\frac23$ になります。

ここには除算が出てこないので、計算コストはそのまま $N^\frac34 / \log N$ です。

Large の寄与を集める部分

$N^\frac14$ から $N^\frac12$ の間の素数 $p$ （ $\Theta(N^\frac12 / \log N)$ 個ある）につき、一回メモリアクセスが発生します。よって、計算コストは $\Theta(N^\frac12 / \log N)$ です。

Small の寄与を集める部分

$N^\frac14$ から $N^\frac13$ の間の素数 $p$ （ $\Theta(N^\frac13 / \log N)$ 個ある）につき、それよりも大きく $N/p/p$ までの素数 $q$ は $O(N/p/p/\log(N/p/p))$ 個あり、その組み合わせの通り数分だけのメモリアクセスと除算が発生します。

積分を計算すると、除算の回数は $\Theta(N^\frac34 / (\log N)^2)$ 回のようです。 $q$ は $p$ より大きくなければいけない、という制約は $\Theta((N^\frac14/\log N)^2)$ 回くらい計算コストが減るだけなので、オーダーに与える影響はありません。

ここでは毎回除算が発生するので、計算コストは $\Theta(N^\frac34 / \log N)$ 程度です。

計算量解析のまとめ

実際に $N = 10^{11}$ で呼び出したときの更新回数を参考に載せました。定数まで合わせた実行回数の見積もりも載せておきます。

素数計数関数由来で出てくる $\log x$ は、 $P(x) = \ln(x/e)$ で見積もっておくと、定数がほぼ一致します。例えば、 $P(N^\frac12) = 11.7$ を使うと、 $\pi(N^\frac12) = 27293 = 1.007 N^\frac12/P(N^\frac12)$ と近似できます。また、 $P(N^\frac14) = 5.33$ を使うと、 $\pi(N^\frac14) = 102 = 0.967 N^\frac14/P(N^\frac14)$ と近似できます。その他、 $P(N^\frac13) = 7.44$ を使うと、 $\pi(N^\frac13) = 626 = 1.003 N^\frac13/P(N^\frac13)$ と近似できます。

これらは、積分区間が最小の素数2から始まらない時は使わないほうがいいです。素数密度に由来する $\log x$ はそのままで正しいです。

エラトステネスの篩について、 $\log\log N$ 部分は $\Sigma_{p = 2, 3, 5, \dots}1/p$ に由来し、最初100個の素数で2.10（素数2の寄与を取り除くと1.60）くらいです。

場所	オーダー	実際の回数	見積もり
初期化	$\Theta(N^\frac12)$ _回除算	158114	$N^\frac12 / 2 = 158114$
エラトステネスの篩	$\Theta(N^\frac12 \log\log N)$	242093	$N^\frac12 / 2 \times 1.60 = 253000$
Large ← Largeの更新	$\Theta(N^\frac12/\log N)$	105604	$N^\frac12(\phi(3\#)/3\#/3 + \phi(5\#)/5\#/5)$ ${}+\frac12 N^\frac12(1/\log7 - 1/\log N^\frac16)$ ${}+(\log2)(N^\frac12/\log N^\frac12)$ $= 113000$ ※
Large ← Smallの更新	$\Theta(N^\frac34 / (\log N)^2)$ _回除算	3452671	$\frac12 N^\frac34 / 5.33 / 5.33 = 3130000$ 以上 $0.56 N^\frac34 / 5.33 / 5.33 = 3510000$ 以下
Small ← Smallの更新	$\Theta(N^\frac34 / \log N)$	11472474	$\frac23 N^\frac34 / 2 / 5.33 = 11200000$
Largeの寄与回収	$\Theta(N^\frac12 / \log N)$	27191	$N^\frac12 / 11.7 = 27000$
Smallの寄与回収	$\Theta(N^\frac34 / (\log N)^2)$ _回除算	2025525	$N^\frac34 / 5.33 / \log N^\frac14$ ${}- N^\frac23 / 7.44 / \log N^\frac13 = 2060000$

※小さい素数での反復の寄与がかなりある一方、積分による評価は漸近評価なので小さな素数については誤差が大きいです。そこで、素数3と素数5については積分せずに解析的にわかる値を使い、積分区間は7からとしました。

なんとか全部の計算回数の評価を数パーセント以内の誤差で抑えることができました。

こうしてみると際立つのですが、除算が発生するところではオーダーが真に $\Theta(N^\frac34 / \log N)$ より小さいことがわかります。除算のコストを $O(1)$ だと考えている場合はroughsを導入しても全体の計算量のオーダーは $\Theta(N^\frac34 / \log N)$ のままですが、除算のコストが $f(N) \subset \omega(1)$ かつ $f(N) \subset O(\log N)$ だとすると元々のアルゴリズムだと計算量のオーダーが $\Theta(f(N) N^\frac34 / \log N)$ だったのが $\Theta(N^\frac34 / \log N)$ に下がっているので、真にオーダーを改善したと言えそうです。なお、「Small ← Smallの更新」の部分は計算量が $\Theta(N^\frac34 / \log N)$ 時間で全体を支配していますが、非常に単純なループであるため多少はSIMD命令で実行できることも、オーダーの見た目から受ける印象よりも高速に動作することの要因のようです。

つまらない最適化

実は除算の回数は $\Theta(N^\frac12)$ 回に減らせます。作戦としては、除算が必要なところをなんとか逆数乗算にできないか、というものです。逆数は一度計算してしまえばそれ以降は使いまわせるため、除算回数を減らせます。逆数を使いまわして除算回数を減らすため、一部の除算においては割る順番を入れ替えます（順番を入れ替えても結果が変わらないことは上で示した通りです）。

この目的のため、invsという配列を新たに用意します。 invs[k]には、N / roughs[k]が入っています。 roughsを詰めるたびにinvsも詰めていかないといけない点に注意します。

あとは、浮動小数点数レジスタと整数レジスタを行ったり来たりするとかなり時間がかかるので、整数命令だけでやります。精度はよくわかりませんが、 $N=10^{11}$ （37ビット）として、それ＋64ビットの精度（102ビット）を確保しました。 (__uint128_t(1)<<102) / pとやって固定小数点数で逆数を求めます。これは $N^\frac12$ （19ビット）までの数を掛けても128ビットに収まるので大丈夫です。

これで精度が足りるかを確認します。うまく丸めた時に整数除算a / bと結果を一致させられるようにするためには、少なくとも商の $1/b$ の違いが明確に判別できないといけません。また、事前計算された1/bの誤差（たかだか1）は、そのa倍が商の誤差として現れます。つまり、商の固定小数点数表記の最下位ビットに対してa程度の誤差が生じうるということです。この誤差を加味しても $1/b$ の違いが判別できるためには、a * bが精度に対して十分小さければいいということです。今回の除算では、a * bが高々 $N=10^{11}$ なので、小数点以下は102ビットありますが、37ビットくらいは有効な精度ではないということです。逆に言えば、そのような誤差が生じてなお、除算結果に64ビットの精度を維持できるということです。全然余裕です。もしかすると64bit÷32bitの除算でもぎりぎり精度が足りるのかもしれません。

ここで、「うまく丸める」部分を正しく実装しないといけません。逆数をinvp = (__uint128_t(1)<<102) / pで作ってしまい、さらに切り捨て除算をinv_j * invp >> 102として実装してしまうとダメです。どちらも零への丸めになっているので、下方向への誤差が二回重なることにより、真の結果よりも1小さい結果が出てくることがあります。そこで、逆数をinvp = (__uint128_t(1)<<102) / p + 1と作ります。これは上向き丸めしているのとほとんど同じです（ぴったり割り切れた場合だけは上向き丸めよりも大きな値になります）。この計算手順において、途中結果を上向き丸めしてから最終計算結果を下向き丸めした場合、途中で丸めずに最終計算結果を得てから下向き丸めした値以上の値が得られることはよいでしょう（途中結果の最終結果への寄与が常に正であるためです。一般的な数式では、このようなことを言うためには注意深い解析が必要です）。逆に今度は真の値よりも1大きい結果が出てくることがないかが心配ですが、 $a/b$ を64ビットもの精度で表しているため、十分余裕があり、安全です。

この最適化により、Library Checkerにおいて単独最速の提出になっています（が、本当につまらない最適化なのであまりうれしくないです。計算量のオーダーも変わっていないですし）。

なお、手元で確認したところ、このコードで $N = 10^{11}$ の場合を計算するのにかかるサイクル数は20Mサイクルくらいでした。そもそも、上で確認したように、配列を読み出す回数が17.5M回はあるので、本当に1サイクルに1更新くらいの勢いでできているようです。そういえば、逆数乗算になっているところのスループットは乗算が二回必要なため $1/2 \text{ cycle}^{-1}$ しかないはずで、それを加味すると20Mサイクル以上かかりそうです。つまり、どこかの更新は1回あたり1サイクル未満で行われていることになります。機械語コードを見てみると、何をやっているかはよくわかりませんでしたが、Small ← Smallの更新部分がSIMD命令を使っている感じになっていたので、ここの更新スループットが高いのでしょう。 SIMD命令で強引に高速化しているだけというわけでもなく、プログラムの全体を通したIPCも3.9というおそろしく高い数字になっていました（無意味なマイクロベンチマーク以外ではまず見ることのない数字です）。よほどCPUを使い倒すコードになっていたようです。

アルゴリズムの改良の余地の検討

結論から言うと、上記提出のもととなったアルゴリズムが優秀すぎて改良する方法を思いつけませんでした。

Smallの更新部分？

Small[q]の更新部分が最も計算量が多くなっているので、ここを何とか改良できないか考えてみます。この部分を観察してみると、 $p$ 箇所をある数だけ減算する、という計算を繰り返しやっているだけです。ようするに区間一様加算・点取得ができればいいので、双対セグメント木なりFenwick木（Binary Indexed Tree, BIT）を工夫して使うなりすれば、この部分の計算量を削減することができそうです。しかしこれは、Small[q]を読むところに $\Theta(\log N)$ がついてしまうのでダメでした。 10倍くらい遅くなりました…………。

全然理解できていないですが、（上記 $\Theta(N^\frac34/\log N)$ よりもオーダーが改善された） $\Theta(N^\frac23)$ で計算するアルゴリズムもBITを使っているので、これと似たことをやっているのかなぁと思っています。いずれにせよ、 $\Theta(N^\frac1{12})$ の改善のために $\Theta(\log N)$ が消えるのは、 $N = 10^{11}$ 程度ではかなりの痛手で、時間計算量 $\Theta(N^\frac23)$ のアルゴリズムを単純に実装しただけでは高速化は難しそうです。 BITが遅い原因は二つあって、二分木なので木の高さが大きくてたくさんメモリを読まないといけない（多分木にすることで書きが遅くなる代わりに読みを改善可能）ことと、インデクス計算を簡単にした代わりに局所性がいまいち（根付近をまとめて持てば改善可能。参考：キャッシュフレンドリーな二分探索ーデータ構造を再考する | POSTD）なことです。 BITは競プロ界隈では定数倍が軽いと評判のようですが、これはポインタをたどったり（メモリレイテンシが見える）、二択の当てられない分岐が発生したりする（分岐予測ミスペナルティを受ける）タイプの $\Theta(\log N)$ がつくデータ構造・アルゴリズムと比べているのだと思われます（BITはこのどちらも発生しません）。

Small[q]は高々1ずつしか増えない単調非減少数列（ $N^\frac12/2$ エントリ）なので、表現するのに必要な情報量はたったの $N^\frac12/2$ ビット（ $N=10^{11}$ において、20KiB）です。実際にはその32倍ものサイズ（632KiB）で保持しており、簡潔からほど遠いデータ構造で持っているわけですが、それにもかかわらず読み書きが遅いのは何とかできないんでしょうか。というか、簡潔からほど遠いデータ構造で持っているがゆえにビット列操作命令などが活用できなくて遅くなっているとも考えたほうがいいのかもしれません。

とりあえず次は、時間計算量 $\Theta(N^\frac23)$ のアルゴリズムを理解することと、Small[q]をビット配列とかで表現できないか考えてみることをやってみたいかなと思っています。

まとめ

与えられた整数 $N$ に対し、 $N$ 以下の素数の数を数える、時間計算量 $\Theta(N^\frac34 / \log N)$ のアルゴリズム（実用上高速）を理解した
除算を逆数乗算にすることで除算回数を減らすつまらない最適化をすることで、Library Checkerで単独最速を取った
今後の課題：時間計算量 $\Theta(N^\frac23)$ のアルゴリズムの理解

よーる

素数を数えるアルゴリズムの理解

何をしているのか

基本的なアルゴリズム

$S(i, \nu)$ の意味

再帰の様子

第二引数 $\nu$ としてあり得る値

DPで求める

三種類の更新

Large ← Large の更新

Large ← Small の更新

Small ← Small の更新

`roughs`を使うアルゴリズム

`roughs`が出てくる理由

配列のインデクス変換の効率よい実装

`smalls`の方

`larges`の方

$p > N^\frac14$ の計算簡略化

三種類の更新

Large ← Large の更新

Large ← Small の更新

Small ← Small の更新

計算量の解析

余談：除算のコスト

詳細な解析

初期化

エラトステネスの篩部分

Large の更新部分

Large ← Large の更新

Large ← Small の更新

Small の更新部分

Large の寄与を集める部分

Small の寄与を集める部分

計算量解析のまとめ

つまらない最適化

アルゴリズムの改良の余地の検討

Smallの更新部分？

まとめ

何をしているのか

基本的なアルゴリズム

の意味

再帰の様子

第二引数としてあり得る値

DPで求める

三種類の更新

Large ← Large の更新

Large ← Small の更新

Small ← Small の更新

roughsを使うアルゴリズム

roughsが出てくる理由

配列のインデクス変換の効率よい実装

smallsの方

largesの方

の計算簡略化

三種類の更新

Large ← Large の更新

Large ← Small の更新

Small ← Small の更新

計算量の解析

余談：除算のコスト

詳細な解析

初期化

エラトステネスの篩部分

Large の更新部分

Large ← Large の更新

Large ← Small の更新

Small の更新部分

Large の寄与を集める部分

Small の寄与を集める部分

計算量解析のまとめ

つまらない最適化

アルゴリズムの改良の余地の検討

Smallの更新部分？

まとめ

$S(i, \nu)$ の意味

第二引数 $\nu$ としてあり得る値

`roughs`を使うアルゴリズム

`roughs`が出てくる理由

`smalls`の方

`larges`の方

$p > N^\frac14$ の計算簡略化