与えられた整数に対し、以下の素数の数を数えるアルゴリズムを理解します。
Library Checkerで高速な実装を調べてみると、提出107115が最速ですが、このソースコードは未定義動作を含んでいます(ランキング上位を見てみると、62015, 107115, 147826, 150290,150294, 150330, 160040 はすべて同じバグを含んでいます。つまり、62015で埋め込まれたバグのようです)。
上位陣の他の提出もほとんど同じアルゴリズムとなっています(Library Checkerへの最古の提出6832もsmalls
への最適化がされていないことを除いて実質的に同じアルゴリズムに見えます。どうやら昔からあるアルゴリズムのようです。smalls
への最適化は提出10100で出現しています)。
このアルゴリズムの詳しい解説が以下にありますが、もう少し行間を埋めていきたいと思います。
眠れない夜は素数の個数でも数えましょう - えびちゃんの日記
何をしているのか
上記の非常に詳しいスライドを見ればわかりますが、素数計数関数は以下の再帰関数で計算可能です。
この再帰関数を配列を使いまわす動的計画法(DP)で計算するというアルゴリズムになっていますが、いくつか工夫がほどこされています。
基本的なアルゴリズム
の意味
は、2以上以下の自然数のうち、エラトステネスの篩で以下の素数までを使って篩ったとき、残っているものの数です。 残っているものとは、までの素数と、最小素因数がを超えるもの((i+1)-rough number)の数です。 ※後での説明の都合から、X-rough numberには1も含む定義を採用します(1を含まない定義もよく使われます)。それにかかわらず、の意味するところは「2以上以下」なので1はカウント対象外です。
例:です。2以上100以下の自然数のうち4以下の素数までを使って篩ったとき残っているものは、以下の34個です。
再帰の様子
を計算するときの再帰の様子を示しました。 まず、は、以下の最大の素数をとしてと同じなので、縦方向では素数のみを示しました。 また、横方向では、としてあり得る値のみを示しました。 水色の部分は、であるため、以下の最大の素数をとしてと同じになる部分です。
第二引数としてあり得る値
再帰的に呼び出された場合も含めた場合、第二引数としてあり得る値は自明ではありません。 例えば、を自然数としたとき、はどのような値を取りえるでしょうか?
実は、はに等しいことが示せます。 その理由は、次のように説明できます。
まず、はと比べて0以上以下しか小さくありません。 同様に考えれば、 はと比べて0以上以下しか小さくないとわかります。 したがって、はよりも0以上以下しか小さくありません。 よりも0以上1未満しか小さくない非負整数は、唯一だけです。 したがって、はに等しくなります。
余談:符号付きの場合であっても、C言語の除算(零への丸め)ならば符号と絶対値を独立に考えられるので同様の結果が成り立ちます。
よって、がとりえる値の集合はの集合と同じです。 この結果を再帰的に使えば、がとりえる値の集合もの集合と同じということがわかります。
つまり、としてとりえる値は、がとりえる値です。 結局、としてとりえる値の集合は となります(※とは同じである可能性があります)。
DPで求める
これを動的計画法で求める方法を考えます。 非常に単純には、までの素数(縦方向)と (横方向)のすべての組み合わせについて、配列を順番に埋めていけばよいですが、これではメモリを使いすぎです。 幸い、の計算には、の直前の素数をとして、 の結果のみが必要なことが利用できます。 つまり、直前の結果しか必要ないので、配列を使いまわすことができます。 また、の大きい方から順に更新していけば、(表裏の二つの配列を行ったり来たりして使いまわす方式ではなく)単一の配列を上書き更新することができます。
つまり、以下のように実装できます。
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]
を配列でとってしまうと大きすぎです(のメモリ要求が発生してしまいます)。
実際に使う部分は高々2 * v
個しかないので、うまく圧縮したいです。
連想配列を使ってもいいですが、高速化のためには配列で頑張りたいです。
そこで、の場合を格納するSmall[v+1]
と、 の場合を格納するLarge[v+1]
を用意します。
Small
の方では、インデクスはq
をそのまま使えば良さそうです。
Large
の方では、インデクスはN/q
を使えば良さそうです。
このようにすることで、必要なメモリ量を削減することができます。
このアルゴリズムで最も時間のかかる工程はS[q]
を更新する手順で、これは回発生します。
最も時間を要するのは、付近のループです。
の範囲に素数は 個存在しますが、それらについて毎回少なくとも、大きいq
のS[q]
(Large[N/q]
に相当する部分。v
個ある)の更新が必要、という部分が計算量の大半を占めています。
それより大きい部分では、途中のif
文でスキップされる割合が急速に増加するため、としてとりえる数の候補が多いにもかかわらず、計算量を支配しません。
特に、の範囲では、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
が出てくる理由
先ほどの図を見れば、を計算する際に、のすべてが使われるわけではないということがわかります。 例えば、 などは使われていません。 は使われているため、合成数についてが使われない、というわけでもなさそうです。 これをどう考えたらよいでしょうか。
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+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+1)-rough number)の数が入っているはずです。
また、pc
は、prev_pが2の時に0であることからわかるように、までの素数の個数より一つ少ないです。
一方、実は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]
にあることがわかりました。
計算例:の周回で、を使う場合を考えます。まず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
の部分です。
であり、small[d>>1]
は19を返します。
pc
は1
である(0
で始まり、p = 3
のループで一回だけインクリメントされた)ので、引き算すると18
が得られます。
の計算簡略化
あのアルゴリズムを見ると、DPをやり終えてからreturn S[N]
とするのではなく、なんだか最後の方でよくわからない計算をしています。
そもそもよく見ると、DPで更新すべきの範囲はまでだったはずが、までで脱出してしまっています。
これは、先ほど説明した「三種類の更新」のうち、さらにp^4 < N
かで範囲が変わることに対応しています。
三種類の更新
Large ← Large の更新
Large[1] -= Large[p] - c
だけが実行されます。
Large[p]
を更新するチャンスはまだありますが、その寄与は後で考えます。
また、roughs
にはを超えて以下の素数が全て格納されています。
よって、
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 の更新
何もやる必要がありません。
計算量の解析
あの実装の時間計算量は、軽い計算(除算以外の整数演算とメモリアクセス)が 回、重い計算(除算)が回、の合計です。 さすがに除算ともなると軽い計算の倍の計算コストがかかるとみなすのが妥当で、その仮定の下で にちょうど収まっています。
余談:除算のコスト
普通のCPUではnビットの除算命令は時間かけて行われます。 基本的には筆算みたいな感じで上から商を決めていきます。 最近だと1サイクルに6ビットくらい商を決められるらしいです(Radix-64 SRT dividerでググってみましょう)。 これは倍精度浮動小数点数(精度53bit)の除算を十数サイクルでできることと整合性があります(商を6ビットずつ立てると9サイクルかかって、前処理とか後処理とかでもう少し時間がかかる)。
ニュートン法を使えば時間でできそうに見えますが、倍精度(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=1012程度の大きさ同士の切り捨て除算ならばdouble型で計算しても誤差が発生しないことが規格で保証されているためである。
倍精度浮動小数点数を使うと速くなるのは実際かなりそうなんですが、それでも依然として のコストがかかっていると考えたほうが良いです。 いまどきのCPUでは、加算命令やロードストア命令は(スループットとして)1サイクルに複数回できるのに対して、除算命令は(スループットとして)十数サイクルに一回しかできず、何十倍もの差があります。 また、短いビット幅の除算命令の方が有意に早く終わることなども、除算命令のコストがあることを示唆しています。 実際、この実装では倍精度浮動小数点数除算を使っているにもかかわらず、それでもなお除算がボトルネックになっているようです。
詳細な解析
初期化
配列の初期化には、回の除算が必要ですが、他の所に比べたら計算量は全然少ないです。
エラトステネスの篩部分
エラトステネスの篩をまでやっているので、回のメモリアクセスが発生していますが、これも他の所に比べたら計算量は全然少ないです。
Large の更新部分
までの素数(個ある)について、未満の(p+1)-rough numberは 個(あってる?)あり、その組み合わせの通り数分だけの更新が発生します。 積分を計算すると、更新の回数は回のようです。
Large ← Large の更新では除算が不要ですが、Large ← Small の更新では除算が必要です。 計算量への寄与は、に近い素数での更新が大半を占めていて、その場合はほとんど Large ← Small の更新なので、ほぼ毎回除算が発生すると考えてよいです。
よって、計算コストは程度です。
なお、係数はとの間だと思います。 これはBuchstab関数のでの値です。 定性的には、次のように考えればこの係数を導くことができます。 エラトステネスの篩でまで篩った残りが(p+1)-rough numberです。 この時、からまでの(p+1)-rough numberは素数です。 つまり、からの区間では、平均素数密度と(p+1)-rough numberの密度は一致します。 よって、以下の(p+1)-rough numberの密度は、程度となります。
を超えた部分では、(p+1)-rough numberの密度はこれとは異なります。 (p+1)-rough numberの出現位置は、(pの素数階乗)という周期を持ちます。 その周期の中には、個の(p+1)-rough numberが存在します(はオイラーのトーシェント関数です)。 その比は漸近的にとなることが知られているようです(Mertens' third theoremというらしいです。https://arxiv.org/pdf/2308.07570.pdfから見つけました)。 よって、(p+1)-rough numberが大きく偏っていないとすれば、を超えたあたりでの密度はくらいということになります。
結局のところ、としては付近のものが多いので、係数はおそらくになると思います。 ただ、 くらいまでは係数はであり、までの素数はにおいて37個と、までの素数102個の1/3以上を占めるので、実際にはより少し大きめと思ったほうが良さそうです。
ところで、までの場合、はよりも小さいです。 小さいで実験していて、からの区間こそが(p+1)-rough numberの密度が高い地点になっていて合わないなぁと悩んでいました。 あくまで漸近挙動なので小さいではとの大小関係という定性的な性質ですらも成り立たないようです。
Large ← Large の更新
までの素数(個ある)について、からまでの(p+1)-rough numberが何個あるかを求めます。 まず、までの素数については、上限であるはより大きいので、個あるとして良さそうです(あってる?)。 それ以上の素数に関しては、上限であるはより小さいので、(p+1)-rough numberというのは全て素数です。 したがって、個くらいあります。
この組み合わせの通り数分だけの更新が発生します。 積分を計算すると、までの素数での更新が回、そこからまでの素数での更新が回、それぞれ発生するようです。 よって、合計の更新回数は回のようです。
係数は、前半部分がくらい、後半部分がでしょうか?
Large ← Small の更新
Large ← Large の更新は、オーダー的にとるに足りないことを上で示したため、Largeの更新回数とオーダーは一致します。 つまり、更新の回数は回のようです。
Small の更新部分
までの素数(個ある)について、からまでの自然数は個あり、その組み合わせの通り数分だけの更新が発生します。 積分を計算すると、更新の回数は 回のようです。 まで計算しないといけない、という部分と、までは更新しなくてよい、という部分は、いずれも同じオーダーです。 前者の係数はであり、後者の係数はです。 よって、全体の係数はになります。
ここには除算が出てこないので、計算コストはそのままです。
Large の寄与を集める部分
からの間の素数(個ある)につき、一回メモリアクセスが発生します。 よって、計算コストはです。
Small の寄与を集める部分
からの間の素数(個ある)につき、それよりも大きくまでの素数は個あり、その組み合わせの通り数分だけのメモリアクセスと除算が発生します。
積分を計算すると、除算の回数は回のようです。 はより大きくなければいけない、という制約は回くらい計算コストが減るだけなので、オーダーに与える影響はありません。
ここでは毎回除算が発生するので、計算コストは程度です。
計算量解析のまとめ
実際にで呼び出したときの更新回数を参考に載せました。 定数まで合わせた実行回数の見積もりも載せておきます。
素数計数関数由来で出てくるは、で見積もっておくと、定数がほぼ一致します。 例えば、を使うと、と近似できます。 また、を使うと、と近似できます。 その他、を使うと、と近似できます。
これらは、積分区間が最小の素数2から始まらない時は使わないほうがいいです。 素数密度に由来するはそのままで正しいです。
エラトステネスの篩について、部分はに由来し、最初100個の素数で2.10(素数2の寄与を取り除くと1.60)くらいです。
場所 | オーダー | 実際の回数 | 見積もり |
---|---|---|---|
初期化 | 回除算 | 158114 | |
エラトステネスの篩 | 242093 | ||
Large ← Largeの更新 | 105604 | ※ |
|
Large ← Smallの更新 | 回除算 | 3452671 | 以上 以下 |
Small ← Smallの更新 | 11472474 | ||
Largeの寄与回収 | 27191 | ||
Smallの寄与回収 | 回除算 | 2025525 | |
※小さい素数での反復の寄与がかなりある一方、積分による評価は漸近評価なので小さな素数については誤差が大きいです。そこで、素数3と素数5については積分せずに解析的にわかる値を使い、積分区間は7からとしました。
なんとか全部の計算回数の評価を数パーセント以内の誤差で抑えることができました。
こうしてみると際立つのですが、除算が発生するところではオーダーが真により小さいことがわかります。
除算のコストをだと考えている場合はroughs
を導入しても全体の計算量のオーダーはのままですが、除算のコストがかつだとすると元々のアルゴリズムだと計算量のオーダーがだったのがに下がっているので、真にオーダーを改善したと言えそうです。
なお、「Small ← Smallの更新」の部分は計算量が時間で全体を支配していますが、非常に単純なループであるため多少はSIMD命令で実行できることも、オーダーの見た目から受ける印象よりも高速に動作することの要因のようです。
つまらない最適化
実は除算の回数は回に減らせます。 作戦としては、除算が必要なところをなんとか逆数乗算にできないか、というものです。 逆数は一度計算してしまえばそれ以降は使いまわせるため、除算回数を減らせます。 逆数を使いまわして除算回数を減らすため、一部の除算においては割る順番を入れ替えます(順番を入れ替えても結果が変わらないことは上で示した通りです)。
この目的のため、invs
という配列を新たに用意します。
invs[k]
には、N / roughs[k]
が入っています。
roughs
を詰めるたびにinvs
も詰めていかないといけない点に注意します。
あとは、浮動小数点数レジスタと整数レジスタを行ったり来たりするとかなり時間がかかるので、整数命令だけでやります。
精度はよくわかりませんが、(37ビット)として、それ+64ビットの精度(102ビット)を確保しました。
(__uint128_t(1)<<102) / p
とやって固定小数点数で逆数を求めます。
これは(19ビット)までの数を掛けても128ビットに収まるので大丈夫です。
これで精度が足りるかを確認します。
うまく丸めた時に整数除算a / b
と結果を一致させられるようにするためには、少なくとも商のの違いが明確に判別できないといけません。
また、事前計算された1/b
の誤差(たかだか1)は、そのa
倍が商の誤差として現れます。
つまり、商の固定小数点数表記の最下位ビットに対してa
程度の誤差が生じうるということです。
この誤差を加味してもの違いが判別できるためには、a * b
が精度に対して十分小さければいいということです。
今回の除算では、a * b
が高々なので、小数点以下は102ビットありますが、37ビットくらいは有効な精度ではないということです。
逆に言えば、そのような誤差が生じてなお、除算結果に64ビットの精度を維持できるということです。
全然余裕です。
もしかすると64bit÷32bitの除算でもぎりぎり精度が足りるのかもしれません。
ここで、「うまく丸める」部分を正しく実装しないといけません。
逆数をinvp = (__uint128_t(1)<<102) / p
で作ってしまい、さらに切り捨て除算をinv_j * invp >> 102
として実装してしまうとダメです。
どちらも零への丸めになっているので、下方向への誤差が二回重なることにより、真の結果よりも1小さい結果が出てくることがあります。
そこで、逆数をinvp = (__uint128_t(1)<<102) / p + 1
と作ります。これは上向き丸めしているのとほとんど同じです(ぴったり割り切れた場合だけは上向き丸めよりも大きな値になります)。
この計算手順において、途中結果を上向き丸めしてから最終計算結果を下向き丸めした場合、途中で丸めずに最終計算結果を得てから下向き丸めした値以上の値が得られることはよいでしょう(途中結果の最終結果への寄与が常に正であるためです。一般的な数式では、このようなことを言うためには注意深い解析が必要です)。
逆に今度は真の値よりも1大きい結果が出てくることがないかが心配ですが、を64ビットもの精度で表しているため、十分余裕があり、安全です。
この最適化により、Library Checkerにおいて単独最速の提出になっています(が、本当につまらない最適化なのであまりうれしくないです。計算量のオーダーも変わっていないですし)。
なお、手元で確認したところ、このコードでの場合を計算するのにかかるサイクル数は20Mサイクルくらいでした。 そもそも、上で確認したように、配列を読み出す回数が17.5M回はあるので、本当に1サイクルに1更新くらいの勢いでできているようです。 そういえば、逆数乗算になっているところのスループットは乗算が二回必要なためしかないはずで、それを加味すると20Mサイクル以上かかりそうです。 つまり、どこかの更新は1回あたり1サイクル未満で行われていることになります。 機械語コードを見てみると、何をやっているかはよくわかりませんでしたが、Small ← Smallの更新部分がSIMD命令を使っている感じになっていたので、ここの更新スループットが高いのでしょう。 SIMD命令で強引に高速化しているだけというわけでもなく、プログラムの全体を通したIPCも3.9というおそろしく高い数字になっていました(無意味なマイクロベンチマーク以外ではまず見ることのない数字です)。 よほどCPUを使い倒すコードになっていたようです。
アルゴリズムの改良の余地の検討
結論から言うと、上記提出のもととなったアルゴリズムが優秀すぎて改良する方法を思いつけませんでした。
Smallの更新部分?
Small[q]
の更新部分が最も計算量が多くなっているので、ここを何とか改良できないか考えてみます。
この部分を観察してみると、箇所をある数だけ減算する、という計算を繰り返しやっているだけです。
ようするに区間一様加算・点取得ができればいいので、双対セグメント木なりFenwick木(Binary Indexed Tree, BIT)を工夫して使うなりすれば、この部分の計算量を削減することができそうです。
しかしこれは、Small[q]
を読むところにがついてしまうのでダメでした。
10倍くらい遅くなりました…………。
全然理解できていないですが、(上記よりもオーダーが改善された)で計算するアルゴリズムもBITを使っているので、これと似たことをやっているのかなぁと思っています。 いずれにせよ、の改善のためにが消えるのは、程度ではかなりの痛手で、時間計算量のアルゴリズムを単純に実装しただけでは高速化は難しそうです。 BITが遅い原因は二つあって、二分木なので木の高さが大きくてたくさんメモリを読まないといけない(多分木にすることで書きが遅くなる代わりに読みを改善可能)ことと、インデクス計算を簡単にした代わりに局所性がいまいち(根付近をまとめて持てば改善可能。参考:キャッシュフレンドリーな二分探索 ー データ構造を再考する | POSTD)なことです。 BITは競プロ界隈では定数倍が軽いと評判のようですが、これはポインタをたどったり(メモリレイテンシが見える)、二択の当てられない分岐が発生したりする(分岐予測ミスペナルティを受ける)タイプのがつくデータ構造・アルゴリズムと比べているのだと思われます(BITはこのどちらも発生しません)。
Small[q]
は高々1ずつしか増えない単調非減少数列(エントリ)なので、表現するのに必要な情報量はたったのビット(において、20KiB)です。
実際にはその32倍ものサイズ(632KiB)で保持しており、簡潔からほど遠いデータ構造で持っているわけですが、それにもかかわらず読み書きが遅いのは何とかできないんでしょうか。
というか、簡潔からほど遠いデータ構造で持っているがゆえにビット列操作命令などが活用できなくて遅くなっているとも考えたほうがいいのかもしれません。
とりあえず次は、時間計算量のアルゴリズムを理解することと、Small[q]
をビット配列とかで表現できないか考えてみることをやってみたいかなと思っています。