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

与えられた整数 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個)

再帰の様子

図1: 再帰の様子

 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)個存在しますが、それらについて毎回少なくとも、大きいqS[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としてあり得る数だけです。 qN/jと表せば、q/pN/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がインデクスになっていますが、かなり謎です。 これは解説(眠れない夜は素数の個数でも数えましょう - えびちゃんの日記)を読まないとわかりませんでした。 これはどうやらdroughsの中で何番目に位置するかを計算しているようです。

読み出すデータは全て更新前のデータ、つまり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の数です。 つまり、dsmall[d>>1] - pc番目の、2以上d以下の(prev_p+1)-rough numberです。 0-indexとの差分は、roughs[0]1が居座っていることで打ち消され、結局droughs[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を返します。 pc1である(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=1012程度の大きさ同士の切り捨て除算ならば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)アルゴリズムの理解

MN-Core 2の機械語ゴルフに入門する

今週の水曜日にMN-Core勉強会があり、そこでMN-Core 2のエミュレータが公開されたようです。 そこで、MN-Core 2の機械語アセンブリ言語がvsmと呼ばれているらしい?)のコードゴルフをやってみます。 MN-Core 2の機械語なんて書いたことないので、丁寧に勉強していきます。

入門書を読む

MN-Core2 vsm ゴルフ · GitHub

この入門書に書いてある、hello.vsm(下記コード)を見ながらやっていきます。

imm ui"0" $n0
imm ui"0x48656c6c" $r2
imm ui"0x6f2c2077" $r3
imm ui"0x6f726c64" $r0
imm ui"0x21000000" $r1

sor $subpeid $n0 $omr1
nop
lpassa $lr2 $lr0/$imr1

nop/2

# d get $lr0n0c0b0m0 1

l1bmrsbor $llr0v $llb0
nop/3

# d get $lb0n0c0b0 1

l2bmrsbor $lb0 $lc0

# d get $lc0n0c0 1

nop/3
mvp/n64i01 $lc0@.0 $p0
nop; wait i01

読み解いていきます。

定数定義

imm ui"0" $n0
imm ui"0x48656c6c" $r2
imm ui"0x6f2c2077" $r3
imm ui"0x6f726c64" $r0
imm ui"0x21000000" $r1

ここは雰囲気でわかりますが、即値を代入する命令が並んでいるようです。 ちゃんとマニュアルを読んでいくと、

  • 即値はALUがimm命令により出力する値である。ダブルクオートが必須なので注意する。(3.2.2 即値)
  • imm命令のオペランドにできるのは)単精度および半精度整数(3.2.2 即値)
    • 単精度整数というのは32bit整数のことを指しているらしい(1.3 語長と演算精度)
  • ペイロードとなるリテラルの解釈は、<integer-number-literal>の場合……(3.2.2 即値)
    • <integer-number-literal>はどこにも定義されていなくて詳細不明
  • $[(l|ll)(r|s|m|n)<addr>によりSRAMが参照できる(?)(3.4.3 Debug get文)
    • より厳密な定義は3.6.1に書かれているようだけれど……
    • 基本的には以下を覚えておけば良さそう?

というわけで以下のような感じのようです。

imm ui"0" $n0          # LM1[0]  = (u32)0
imm ui"0x48656c6c" $r2 # GRF0[2] = (u32)"Hell"
imm ui"0x6f2c2077" $r3 # GRF0[3] = (u32)"o, w"
imm ui"0x6f726c64" $r0 # GRF0[0] = (u32)"orld"
imm ui"0x21000000" $r1 # GRF0[1] = (u32)"!\0\0\0"

PE番号判別

sor $subpeid $n0 $omr1
nop

早速よくわかりません。

subpeidは、自身のPE番号を示す2bit値だそうです(3.6.1.20 固定値入力オペランド)。 SIMDではなくスカラで頑張ろうということでしょうか。

sor命令は、半精度('s'hort)のビット論理和(or)命令のようです(3.6.12 ALU命令式・表3.9・表3.10)。 この命令は、output is all 0の時、フラグが1になるようです(表3.11)。 つまり、PE番号が0の時だけフラグが1となり、PE番号が1~3の場合にはフラグが0になるようです。

どうやら、PE番号に依存したマスクを作る命令だったようです。 PE番号が0のPEだけ動かしていくといった感じのようです。

ここにnop命令が必要な理由は分かりませんでした。

命令間で適切にステップ数を空けなければならない場合がある。(3.6.3 ハザードの回避)

とあるので、nop命令が必要な局面があることは分かりますが、マスクレジスタについてこのような待ちが必要な状況が存在するのか、明確な記述は見つけられませんでした。 ただ、

マスクレジスタは通常のデータではなくフラグを保持する特殊なPEメモリであり、(3.6.2 マスクレジスタ

PE メモリへの書き込みを伴う命令は完了まで6サイクルを要する。(3.6.3.8 PEメモリ書き込み ⇒ PEメモリ読み出し)

とあるので、マスクレジスタについてもGRFと同様6サイクルの待ち時間が必要なのかもしれません。 次節の場合と異なりnop命令が一つで済んでいるのは、「iサイクル目(i=0,1,2,3)用」のマスクレジスタ、のような仕組みになっているからだと思われます(参考:3.6.2 マスクレジスタ における「サイクル方向」の記述)。 つまり、「0サイクル目用」への書き込みは第0サイクル、読み出しはnop命令が一つあれば第8サイクル、と7サイクル空き、6サイクルの待ち時間を満たせるからと思われます(「1サイクル目用」「2サイクル目用」「3サイクル目用」もまったく同様です)。

(2024/02/26 20:13追記ここから) Jun Makinoさんから、Xにて以下の指摘をいただきました。

実際、このnopを取り除いても正しく動作しました。 なんかアセンブルエラーになる気がしていましたが、勘違いか何か他の要因でのアセンブルエラーだったのかもしれません。 したがって、正しくは「マスクレジスタはPEメモリではあるが、次のステップで読み出すことが可能である。」となるようです。 (2024/02/26 20:13追記ここまで)

VLIW命令セットの機械語を書いたことがないので、このようなnopを入れないと正しく動かないみたいなことが書かれているマニュアルを見ていると新鮮です。

sor $subpeid $n0 $omr1 # M[1] = PEID == 0 ? 1 : 0
nop # マスクレジスタへの書き込みは6サイクル待つ必要がある(?)ので1NOP

条件転送

lpassa $lr2 $lr0/$imr1

nop/2

lpassa命令は、倍精度('l'ong)のコピー命令のようです(3.6.12 ALU命令式)。 なぜpassaなのかはよくわかりませんが、GRAPE-DRの論文[1]に記載されているGRAPE-DRアセンブリ*1にもpassaが出現する(何の命令であるかは不明)ので、伝統のある命令名のようです。

(2024/02/26 20:13追記ここから) 百千万億 萬(つもいよろず)さんからXにて以下の指摘をいただきました。

passa命令は、pass a(a入力をそのまま転送)の意味のようです。 (2024/02/26 20:13追記ここまで)

/$imr1は、その直前に書かれているオペランド$lr0=GRF0[0])への書き込みをマスクレジスタ1番によりマスクする、といった意味のようです。 $omr1$imr1で別表現ですが、これは「マスクレジスタの書き込み時は、$omr<addr>」「マスクレジスタを使う時は$imr<addr>」という使い分けになっているだけで、同じマスクレジスタを指しているようです。

マスクはマスクレジスタ書き込み時にも適用できることに注意する。(3.5.1.13 omr - マスクレジスタへの書き込み)

と書いてあるので、この時に見間違えないように表記を変えているのでしょうか。

また、GRFのアドレスは32bit単位で指定するようですが(3.6.1.10 r - GRF0)、長語(64bit)で転送しているため、GRF0[2]を読み出しているように見えてGRF0[2~3]を読み出しているようです。 同じく、GRF[0]に書き込んでいるようで、GRF[0~1]に書き込んでいるようです。 これも明示的な記述がなくてわかりませんでしたが、普通のメモリと同じくそうなっているのだと思います。

結局、GRF0[2~3]の値をGRF0[0~1]に書き込む(ただし、PE番号が0の場合のみ)、ということのようです(出力オペランドが後ろのようです。これは出力オペランドが複数になりうる(3.6 PE命令文)ことへの対応でしょうか?)。 これにより、PE番号0の場合はGRF0[0~1]が"Hello, w"に、PE番号1~3の場合はGRF0[0~1]が(そのまま変わらず)"orld!\0\0\0"になった状態になりました。 なかなか技巧的ですね。

さて、この操作の後にnop/2があります。 このnopの後に続く/は、nop命令を何回繰り返すかを示すようです。

GRF0/GRF1に書き込んだ後は、データ競合により1ステップ空けてからでないと同じアドレスからは読み出しを開始できない。(3.6.3.8 PEメモリ書き込み ⇒ PEメモリ読み出し)

とのことなので、nop命令が一つ必要なのはわかりますが、二つ必要なのは不思議です。 これはどうやら、

PE命令は1命令でL2B以下全体を4サイクル制御する(1.2 機械語命令概観)

に由来するようです。 この文章は分かりづらいですが、要するにPE命令は長さ4のベクトル命令であるということだと思います。 これにより、GRF0[0]へ第28サイクル、第29サイクル、第30サイクル、第31サイクルに書き込みアクセスされることになります。 nop命令が一つだと、この後第36サイクルにGRF0[0]からの読み出しが発生します。 これは間が6サイクル空いていないので、データ競合となります(3.6.3.8 PEメモリ書き込み ⇒ PEメモリ読み出し)。 実際、lpassa $lr2 $lr0/0100に変えてみるとnop命令が一つでもアセンブルエラーとはなりませんでした(参考:3.6.3.8 PEメモリ書き込み ⇒ PEメモリ読み出し)。

結局、ステップの後半でGRF0[0]周辺に書き込みアクセスしなければよいわけです。 そこで、後述のストライドアクセス指定を使ってlpassa $lr2v4 $lr0v4/$imr1としてみたところ、nop命令が一つでもアセンブルエラーとなりませんでした。 やはり、本質的にはnop命令が一つでも大丈夫なようです。 これはMN-Core 2の機械語ゴルフでnop命令を減らすテクニックとなるかもしれません。

lpassa $lr2 $lr0/$imr1 # GRF0[2~3] -> GRF0[0~1] if M[1] 

nop/2 # 最低6サイクル空ける必要があるためnopを二つ

L1Bへの転送

l1bmrsbor $llr0v $llb0
nop/3

(後ろに$llで始まるオペランドが来る)l1bmrsbor命令は、L1B命令の縮約命令(3.6.8.13 l1bmr<rrn_opcode> - 2長語16x1縮約)で、半精度('s'hort)のbitwise logical or(bor)演算指定のもののようです(3.5.5 縮約演算指定)。

各L1BMについて、各MAB配下の4PEからサイクルあたり2長語ずつを読み出し、MAB方向には縮約、PE方向には結合してL1BMにサイクルあたり8長語で書き込む。(3.6.8.12 l1bmr<rrn_opcode> - 1長語16x1縮約)

と書いてありますが意味が分かりません。 まず、「縮約」というのは「総和を求める」とかそういう操作のことのようです。 「MAB方向」「PE方向」というのがどの方向のことを言うのかはよくわかりませんが、疑似コードを見る限り、次の操作を行うようです(3.3 疑似コードによる命令動作の記述 に記法が載っている)。

  • 以下を4サイクルくりかえす。以下は第iサイクルの動作を記述する(i = 0, 1, 2, 3)。ストライド指定などによりiに依存しうる。
    • L1B配下の16個のMABに対し、各MABの0番目のPEの<src>~<src+1>(64bit)に書かれた値の総和(など)を求め、L1Bメモリの<addr_b + 8 * i>番地に書き込む
    • L1B配下の16個のMABに対し、各MABの0番目のPEの<src+2>~<src+3>(64bit)に書かれた値の総和(など)を求め、L1Bメモリの<addr_b + 8 * i + 4>番地に書き込む
    • L1B配下の16個のMABに対し、各MABの1番目のPEの<src>~<src+1>(64bit)に書かれた値の総和(など)を求め、L1Bメモリの<addr_b + 8 * i + 1>番地に書き込む
    • L1B配下の16個のMABに対し、各MABの1番目のPEの<src+2>~<src+3>(64bit)に書かれた値の総和(など)を求め、L1Bメモリの<addr_b + 8 * i + 5>番地に書き込む
    • L1B配下の16個のMABに対し、各MABの2番目のPEの<src>~<src+1>(64bit)に書かれた値の総和(など)を求め、L1Bメモリの<addr_b + 8 * i + 2>番地に書き込む
    • L1B配下の16個のMABに対し、各MABの2番目のPEの<src+2>~<src+3>(64bit)に書かれた値の総和(など)を求め、L1Bメモリの<addr_b + 8 * i + 6>番地に書き込む
    • L1B配下の16個のMABに対し、各MABの3番目のPEの<src>~<src+1>(64bit)に書かれた値の総和(など)を求め、L1Bメモリの<addr_b + 8 * i + 3>番地に書き込む
    • L1B配下の16個のMABに対し、各MABの3番目のPEの<src+2>~<src+3>(64bit)に書かれた値の総和(など)を求め、L1Bメモリの<addr_b + 8 * i + 7>番地に書き込む

<src>には$llr0vが指定されています。 このうしろについているvは、ストライドアクセスを意味するようです。 インクリメント量が省略された表記であり、アクセス語長での1語分(3.6.1.10 r - GRF0)ずつのストライドアクセス(つまりストリームアクセス)を指定しているようです。

結局、

  • L1Bメモリの0番地には、PE番号0のPEのGRF0[0~1](64bit)から読み出した値が書かれる(16個全部のMABが同じ値をL1Bに送って、それが全部ビット毎論理和されるので結局値は変わらない)
    • つまり、"Hello, w"が書き込まれる
  • L1Bメモリの1番地には、PE番号1のPEのGRF0[0~1](64bit)から読み出した値が書かれる(同)
    • つまり、"orld!\0\0\0"が書き込まれる
  • L1Bメモリの2番地には、PE番号2のPEのGRF0[0~1](64bit)から読み出した値が書かれる(同)
    • つまり、"orld!\0\0\0"が書き込まれる
  • L1Bメモリの3番地には、PE番号3のPEのGRF0[0~1](64bit)から読み出した値が書かれる(同)
    • つまり、"orld!\0\0\0"が書き込まれる
  • L1Bメモリの4番地には、PE番号0のPEのGRF0[2~3](64bit)から読み出した値が書かれる(同)
    • つまり、"Hello, w"が書き込まれる
  • L1Bメモリの5番地には、PE番号1のPEのGRF0[2~3](64bit)から読み出した値が書かれる(同)
    • つまり、"Hello, w"が書き込まれる
  • L1Bメモリの6番地には、PE番号2のPEのGRF0[2~3](64bit)から読み出した値が書かれる(同)
    • つまり、"Hello, w"が書き込まれる
  • L1Bメモリの7番地には、PE番号3のPEのGRF0[2~3](64bit)から読み出した値が書かれる(同)
    • つまり、"Hello, w"が書き込まれる

となって先頭に"Hello, world!\0"が完成しました(近くのメモリを破壊していますが)。

GRF[2] = "Hell", GRF[3] = "o, w", GRF[0] = "orld", GRF[1] = "!\0\0\0"のようなミドルエンディアン的なレジスタの使い方をしていたのは、こういう理由だったのですね。

また、今度はnop/3が挿入されています。 これは次のL2Bメモリへの転送までに10サイクル空ける必要がある(3.6.3.6 PE → L1BM転送 ⇒ L1BM → L2BM転送/内部マルチキャスト)ためのようです。

l1bmrsbor $llr0v $llb0 # L1BM[0] = PE[0].GRF[0~1], L1BM[1] = PE[1].GRF[0~1], ...
nop/3 # 10サイクル待つ必要がある

L2Bへの転送

l2bmrsbor $lb0 $lc0
nop/3

同じパターンですね。L1Bメモリの内容をL2Bメモリに転送しています。

nopが三つ必要な理由は分かりませんでした。 さすがにnopが一つもないとアセンブルエラーになります(たぶん4サイクルでは浮動小数点数加算での縮約は無理なのでしょう)が、nopが一つあればアセンブルエラーにはなりませんでした。 caddyを用いたジャッジでも問題なく動作したため、一つで十分と思われます。

l2bmrsbor $lb0 $lc0 # L2BM[0] = L1BM[0], L2BM[1] = L1BM[1], ...
nop/3 # 1つは必要だけれど3つはいらない(?)

出力(メモリへの書き込み)

mvp/n64i01 $lc0@.0 $p0
nop; wait i01

(第一オペランド$lcで始まり、第二オペランド$pではじまり@がつかない)mvp命令は、L2Bメモリの内容をPDMに転送する命令です(3.5.8.5 L2BM → PDM並列個別転送命令)。 ここで、PDMとはPCIe Interface Unit Data Memoryのことで、SRAMだそうです(1.1 構成)。 第0番グループのPDMはホストとPCIeインタフェースで接続(1.1 構成)されるらしいので、ここがMN-Core 2的な入出力になるのでしょう。 命令の動作としては、0番のL2Bについて、L2Bメモリの0番地から64bit値を64個、PDMの0番地に書き込む、という動作になるようです。

これだけでは実際にはデータの転送が終わっていない段階でプログラムが終了してしまうのでしょう。 これをふせぐために、wait i01があるようです。 wait命令を使うことで、mvpが終わるのを待つことができるようです(3.6.13 wait - PE命令にMV命令を待機させる)。 どれを待つかも柔軟に指定できるようで、mvp命令の後ろに書かれているi01部分が指定されています。 タグは0番はだめ(3.6.13 wait - PE命令にMV命令を待機させる)らしいので、1番を使っているのでしょう。

mvp/n64i01 $lc0@.0 $p0 # PDM[0] = L2B[0].L2BM[0~63]
nop; wait i01 # データの転送完了を待つ

まとめると、以下のような感じです。

imm ui"0" $n0          # LM1[0]  = (u32)0
imm ui"0x48656c6c" $r2 # GRF0[2] = (u32)"Hell"
imm ui"0x6f2c2077" $r3 # GRF0[3] = (u32)"o, w"
imm ui"0x6f726c64" $r0 # GRF0[0] = (u32)"orld"
imm ui"0x21000000" $r1 # GRF0[1] = (u32)"!\0\0\0"

sor $subpeid $n0 $omr1 # M[1] = PEID == 0 ? 1 : 0
nop # マスクレジスタへの書き込みは6サイクル待つ必要がある(?)ので1NOP

lpassa $lr2 $lr0/$imr1 # GRF0[2~3] -> GRF0[0~1] if M[1] 

nop/2 # 最低6サイクル空ける必要があるためNOPを二つ

l1bmrsbor $llr0v $llb0 # L1BM[0] = PE[0].GRF[0~1], L1BM[1] = PE[1].GRF[0~1], ...
nop/3 # 10サイクル待つ必要がある

l2bmrsbor $lb0 $lc0 # L2BM[0] = L1BM[0], L2BM[1] = L1BM[1], ...
nop/3 # 1つは必要だけれど3つはいらない(?)

mvp/n64i01 $lc0@.0 $p0 # PDM[0~63] = L2B[0].L2BM[0~63]
nop; wait i01 # データの転送完了を待つ

ゴルフしていく

出力の仕組みが分かったので、簡単なゴルフ問題から解いていきます。

anarchy golf - echo

1Byteのvsmファイルで正答できます。 一応endless問題なのでネタバレしないようにしておきます(といっても256通りしかないんですが)。

anarchy golf - Lower ASCII value

Post mortemの中で、出力を並列に計算できる問題です。 入力を読んだ後に計算する必要がありますが、計算の単位は2Byteなので二冪でありMN-Core 2の機械語コードで処理するのが容易そうです。

データレイアウトを考える

並列に計算できるのはいいですが、入力の2Byteが出力の1Byteに対応するので、「分配したデータをそのまま結合するとレイアウトを含めて同一のデータが上位階層に完成するようになっている。」(1.2 機械語命令概観)と相性が悪いのをどうにかしないといけなさそうです。 幸い、「L1BM 側アクセス速度と PE 側アクセス速度がともに同じである転送命令は同一の種別となる。同一種別の転送ではレイアウトが保たれる。」(3.6.8.2 L1BM命令式種別と折り返し動作)と書いてあり、逆に言えばアクセス速度が異なる転送命令を組み合わせればデータレイアウトを変えることができそうです。

表3.8をにらんでもよくわからないので、L1BM→PEの転送モードの動作を書き出してみました。

L1BM→PE(行き)

  • 1長語PE放送:L1BMの64bitを、64個あるPEに配る。並列処理には向かなさそう
  • 2長語PE放送:L1BMの連続する128bitを、64個あるPEに配る。並列処理には向かなさそう
  • 1長語16x1MAB放送:L1BMの連続する256bitを、16個あるMABに配る。各MAB内の4つのPEは、おのおの異なる64bitを受け取る。最悪これを使えば16個あるMABが全部同じ動作をしてHello, world!みたいな感じに出力を作れそう
  • 2長語16x1MAB放送:L1BMの連続する512bitを、16個あるMABに配る。各MAB内の4つのPEは、おのおの異なる128bitを受け取る
  • 1長語4x4MAB放送:L1BMの連続する1024bitを、MAB4個のグループ(4つある)に配る。グループ内のPE(16個ある)はおのおの異なる64bitを受け取る
  • 2長語4x4MAB放送:L1BMの連続する2048bitを、MAB4個のグループ(4つある)に配る。グループ内のPE(16個ある)はおのおの異なる128bitを受け取る
  • 分配:L1BMの連続する4096bitが読み出され、64個あるPEはおのおの異なる64bitを受け取る

データレイアウトを考える(再)

以下の組み合わせを使えばうまくいきそうともくろみました。

  • L1BM→PEでは、「2長語4x4MAB放送」を使う。L1BM側速度32に対してPE側速度2
    • 最高速なのは「分配」(64:1)だが、半分サイズで受け取る命令が存在しないので、入力と出力でデータレイアウトを変えることができない
    • 4倍無駄に同じ計算をすることになるがしょうがない
  • PE→L1BMでは、「1長語4x4MAB縮約」を使う。L1BM側速度16に対してPE側速度1
    • 対応する命令の半分のサイズを受け取る

コードを書く

データ読み込み部分

上記方針だと、L1Bあたり2048bit(256Byte)を処理して128Byteを生成することができます。 入力は500Byteくらいなので、2サイクルで全部の処理が完了します。 PE命令が長さ4のベクトル命令なので、4長語の処理まではノーコストです。 PEの担当する処理がちょうど2長語×2サイクル、なのでいいかんじです。

まず、上記方針であっているかデバッグプリントで調べます。

mvp/n64i01 $p0 $lc0@.0
nop; wait i01

l2bmb $lc0 $lb0
nop/3

l1bmm4 $llb0 $llr0v
nop/3

d get $lr0n0c0b0 2

これを実行した結果、PEでは以下のように受け取れたようです。

  • MAB番号0, 1, 2, 3のPE0のGRF[0~1] = 0x4974206973206120 = "It is a "
  • MAB番号0, 1, 2, 3のPE1のGRF[0~1] = 0x706572696F64206F = "period o"
  • MAB番号0, 1, 2, 3のPE2のGRF[0~1] = 0x6620636976696C20 = "f civil "
  • MAB番号0, 1, 2, 3のPE3のGRF[0~1] = 0x7761722E20526562 = "war. Reb"
  • MAB番号0, 1, 2, 3のPE0のGRF[2~3] = 0x656C207370616365 = "el space"

……

  • MAB番号4, 5, 6, 7のPE0のGRF[0~1] = 0x696464656E206261 = "idden ba"(64Byte目付近)

……

  • MAB番号0, 1, 2, 3のPE0のGRF[4~5] = 0x7265642073706163 = "red spac"(256Byte目付近)

……

  • MAB番号12, 13, 14, 15のPE2のGRF[6~7] = 0x2E2E2E2E00000000 = "....\0\0\0\0"(末尾)

えー、そういう感じで分配されるのかぁ。 連続する2長語がGRF[0~3]に入ることを期待していましたが、そうではなくて、連続する4長語が4つのPEのGRF[0~1]に入り、その次の連続する4長語が4つのPEのGRF[2~3]に入り、その次の連続する4長語は次のMABに送られて……という感じのようです。 つまり、AoSではなくSoAになっています。 まぁ計算機の命令セットとしてはSoAの方が簡単なので、そういう仕組みになっているのでしょう。

PE内ロジック

ここからPE内のロジックを記述していきます。 8bitシフトとマスクを使って整えた後、min命令を16bit×4のSIMD命令として使えばきれいに書けそうです。

即値付きの演算命令なんていう高級(?)なものはないようなので、シフト量とマスクは定数としておいておく必要があります。 ここで、即値を置く場所が重要です。 なぜなら、

複数の命令式が同一のPEオペランドから読み出している場合、それら全てでアクセスする領域が全サイクルで同一である(3.6.4 並列実行条件)

という制約があるからです。 この記述は意味不明ですが、命令をバンドル(VLIW命令)にまとめられる条件として、同じPEメモリ(GRF0とかLM1とか)に対して2read以上ができないということを意味しているようです(MN-Core Architecture Deep Diveの50ページ付近)。 明示的には書かれていないですが、この制約は2オペランド命令の第一オペランドと第二オペランドが同じPEメモリであってはいけないという制約にもなっているようです(実際アセンブルエラーになりますし、SRAM回路的にもそれが自然です)。

さて、smin命令を使うところまで書きました。 同一のPEオペランドが使えないというのがなかなか厳しい制約です(自分でバンク調停をやっていることになるので)。 LMは書き込みから読み込みまで2ステップ空けなければならず、GRFは書き込みから読み込みまで1ステップ空ければよいので、二オペランド命令を使う時は

  • 第一のステップで、LMに書き込む
  • 第二のステップで、GRFに書き込む
  • 第三のステップは、おやすみ(nop
  • 第四のステップで、LMとGRFをオペランドにして命令を実行

とやればポート衝突とデータ競合と1read制約を全部満たせていいかんじです。

imm us"8" $ln0
imm us"0xff" $lm0

mvp/n64i01 $p0 $lc0@.0
nop; wait i01

l2bmb $lc0 $lb0
nop/3

l1bmm4 $llb0 $llr0v
nop/3

# d get $lr0n0c0b0 2

llsr $lr0v $ln0 $ls0v
land $lr0v $lm0 $ln10v
land $ls0v $lm0 $lr10v
nop
smin $ln10v $lr10v $ls10v
nop

……。 "immop can't output to LM0"というエラーが出てしまいました。 そのような記述はソフトウェア開発者マニュアルは見当たりませんでしたが、アセンブルエラーになるのでそうなのでしょう。

即値命令が発行されているならば LM0 がアクセスされていない(3.6.4 並列実行条件)

との記載がありました。並列実行条件じゃないじゃん!わかりづらすぎる……。

マスク値をLM1に移動させます。 その影響で他の部分のレジスタ割り付けも変わってしまいました。

imm us"8" $ln0
imm us"0xff" $ln2

mvp/n64i01 $p0 $lc0@.0
nop; wait i01

l2bmb $lc0 $lb0
nop/3

l1bmm4 $lb0 $llr0v
nop/3

# d get $lr0n0c0b0 2

llsr $lr0v $ln0 $ls0v
land $lr0v $ln2 $lm10v
land $ls0v $ln2 $lr10v
nop
smin $lm10v $lr10v $ls10v
nop

d get $lr0n0c0b0m0p0 1
d get $ls0n0c0b0m0p0 1
d get $lm10n0c0b0m0p0 1
d get $lr10n0c0b0m0p0 1
d get $ls10n0c0b0m0p0 1

うまくいっていそうです。

次に、データをパックしなおします。 今は16bitのデータが4つで64bitになっていますが、これを8bitのデータ4つで32bitの状態に直します。

# 定数定義に追加

imm ui"0xffff" $ln4
imm us"16" $ln6

# さっきのコードのつづき

land $ls10v $ln4 $lm20v
ilsr $ls10v $ln0 $lr20v
nop
lor $lr20v $lm20v $ls20v
nop
llsr $ls20v $ln6 $lr30v
nop
lor $ls20v $lr30v $lm30v
nop/2

d get $ls10n0c0b0m0p0 1
d get $lm20n0c0b0m0p0 1
d get $lr20n0c0b0m0p0 1
d get $ls20n0c0b0m0p0 1
d get $lr30n0c0b0m0p0 1
d get $lm30n0c0b0m0 1

nopだらけになりましたが、とりあえずコードは書けました。 こういう密な依存があるときはnopが避けがたいのかなと思っていましたが、どうやら$alufというオペランドを指定することでフォワーディングを受け取れるようです(3.6.1.16 aluf - ALU演算結果フォワーディング)。 また、フォワーディングに渡すだけで終了する(PEメモリには書き込まない)場合は$nowriteをデスティネーションに指定する必要があるようです(3.6.1.19 nowrite - フォワーディング用ダミー出力)。 これを使って以下のように書き換えることができました。

imm us"8" $ln0
imm us"0xff" $ln2
imm ui"0xffff" $ln4
imm us"16" $ln6


mvp/n64i01 $p0 $lc0@.0
nop; wait i01

l2bmb $lc0 $lb0
nop/3

l1bmm4 $llb0 $llr0v
nop/3

# d get $lr0n0c0b0 4

llsr $lr0v $ln0 $nowrite
land $aluf $ln2 $lr10v
land $lr0v $ln2 $nowrite
smin $aluf $lr10v $ls10v

land $aluf $ln4 $lr20v
ilsr $ls10v $ln0 $nowrite
lor  $aluf $lr20v $ls20v
llsr $aluf $ln6 $nowrite
lor  $aluf $ls20v $lr30v
nop

d get $lr30n0c0b0m0 2

いいかんじです。

MAB内PE間データ交換

以下の交換をやって、結果の順序を合わせます。

  • PE[0].GRF0[31v4] → PE[0].GRF1[30v2](passa
  • PE[1].GRF0[31v4] → PE[0].GRF1[31v2](msr
  • PE[2].GRF0[31v4] → PE[1].GRF1[30v2](msr
  • PE[3].GRF0[31v4] → PE[1].GRF1[31v2](msr×2)
  • PE[0].GRF0[33v4] → PE[2].GRF1[30v2](msl×2)
  • PE[1].GRF0[33v4] → PE[2].GRF1[31v2](msl
  • PE[2].GRF0[33v4] → PE[3].GRF1[30v2](msl
  • PE[3].GRF0[33v4] → PE[3].GRF1[31v2](passa

これには、以下の順序で命令を並べれば最小命令数になりそうです(問題の入力が短いためベクトル長2で十分であることや、msr×2とmsl×2が同じであることを生かせば、同じPEへの転送や2個となりのPEへの転送を統合してもう一命令ずつ減らせるかもしれません)。

ipassa $r31v4 $s30v2 if PE == 0 (don't care = 1,2,3)
imsr $r31v4 $s31v2 if PE == 0 (don't care = 1,2,3)
ipassa $aluf $s30v2 if PE == 1 (don't care = 2,3)
imsr $aluf $s31v2 if PE == 1 (don't care = 2,3)

ipassa $r33v4 $s31v2 if PE == 3 (don't care = 2)
imsl $r33v4 $s30v2 if PE == 3 (don't care = 2)
ipassa $aluf $s31v2 if PE == 2
imsl $aluf $s30v2 if PE == 2

このマスクを作っていきます。 以下のようにすると前準備なしで必要なマスクが作れます。 また、$alufを除いては副作用がありません。

idec $subpeid $omr1 # [1,2,3]
idec $aluf $omr2 # [2,3]
iand $aluf $aluf $omr3 # [2]

# ↓確認用↓

zero $r0v
ipassa $n4 $r1/$imr1
ipassa $n4 $r2/$imr2
ipassa $n4 $r3/$imr3

d geth $r0n0c0b0m0 4

で、実行してみるとNo valid operation in this line.と怒られます。 どうやらmsr/msl命令は長語の転送しかないという罠を踏んだようです。

例えばmslの演算種別はuntypedなので精度指定は付けられず、(3.6.12 ALU命令式)

交換は長語でやって、短語単位での位置調整(インタリーブするように並べ替え。アンパックというらしい)は別の命令でやることにします。

lpassa $lr30v2 $ls30v
msr $lr30v2 $ls40v
lpassa $aluf $ls30v/$imr1
msr $aluf $ls40v/$imr1

lpassa $lr32v2 $ls40v/$imr2
msl $lr32v2 $ls30v/$imr2
lpassa $aluf $ls40v/$imr3
msl $aluf $ls30v/$imr3
nop
ipassa $s31v2 $s40v2

これではだめだった

で、このコードをcaddyで動かしてみると、先頭の方は正しい答えなのですが、行が変わったあたりから出力が変です。 どうも2Byteを取る組み合わせがおかしいっぽい挙動をしていて、改行コードの違い("\012"と"\015\012"の違い)が原因かと思いましたが、どうやらこの問題は行ごとに解かないといけない問題だった ようです。 つまり、input1の一行目It is a period of civil war. Rebel spaceships, striking from a hidden base, have won their first victory against the evil Galactic Empire.\nは改行文字を含めて奇数文字ですが、次の行は仕切り直して"Du/ri/ng/ t/he/..."の組み合わせで解かなければいけません。

よく確認しておくべきでした。

これに対処するためには、自分の担当文字列に'\n'が何個あるかを数え、それを自分より後ろの文字列の担当者に伝えて適切な量だけシフトすればよさそうです。 つまりprefix sumですが、MAB間の直接通信はないっぽく、縮約や結合はできてもprefix sumを計算するモードはないようなので、全対全の通信(実際には基数4での縮約と放送?)が必要そうです。

また、シフトではみ出した分の“袖交換”も必要です。 一気にたくさん転送するタイプのL1BM命令にはアライン制約があり、ちょっとだけずれた部分を転送するというのが無理そうです。 これに対しては、分配(3.6.8.20 l1bmd - 分配)命令についているラウンドロビンずらし機能を行きでは使って帰りでは使わないみたいなことをすると、ずらしたデータをL1Bメモリ内に生成することができそうです。

ただ、これを完全に実装するのはものすごい労力がかかりそうなのであきらめました。 MN-Core 2の機械語を手で書くことの無謀さがわかってきました。 とりあえず、PE内でのprefix sumだけ実装しました。

prefix sumの実装

以下のような感じでprefix sumできるようです(PEが担当する8文字中、奇数文字目4か所のうちには高々一つの改行しか存在しないと仮定しました。これはテストケースでは満たされます*2)。

なんだか便利そうなTレジスタ(Temporary Registerの略だそうです:1.2 機械語命令概観)をこれまで使っていなかったので、たくさん使ってみました。 書き込んでから読めるまでに1ステップ空けなければいけないので、使い勝手はGRFとあまり変わりませんが、使ってないアドレスに捨てないといけないとか、第一オペランドと第二オペランドを違うPEメモリにしないといけないとかの制約がゆるくなるので楽できます。 1ステップ空けないといけない制約は、そもそも直前のステップの結果は$alufでとれるのであまり気になりません。むしろ、「次の命令と次の次の命令で使いたい」みたいなときに二か所目として使えていいかんじです。 逆に、うしろにvっていうのがつかないのでベクトル命令なのかスカラ命令なのかわかりづらいのがやや不満です。

imm us"0x0a00" $ln8
imm us"0xff00" $ls8
imm us"1" $ln10

# d get $lr0n0c0b0 4

land $lr0v $ls8 $nowrite       # 偶数バイト目をマスク
sxor $aluf $ln8 $omr4          # 2Byteごとに0x0a00であればMR4[i] = 1
zero $ls0v                     # ゼロクリア(たぶん初期化時に0になっているけれど、MR4が利用できるまでに時間が余っているので)
spassa $ln10 $ls0v/$imr4 # 0x0aの所だけを1にする
zero $llt                      # ゼロクリア(たぶん初期化時に0になっているけれど、GRF0[0~7]が利用できるまでに時間が余っているので)
lpassa $ls0v $omr5             # 奇数文字目のどこにも0x0aがなければMR5=1(8Byte単位のマスク)
linc $lt $ls50v                # 64bit整数としての1を作る(imm命令では作れない※)
zero $ls50v/$imr5              # マスク付きの0書き込み=奇数文字目のどこかに0x0aがあれば1のまま
nop                            # $ls50vが利用可能になるまで待つ
msl $ls50v $lt                 # 次の8Byteを担当するPEに渡す
ladd $ls50v $aluf $lr50v/$imr1 # PE番号が1,2,3なら足す
msl $lt $lt                    # 受け取ったのをさらに次の8Byteを担当するPEに渡してT regは上書き
ladd $lr50v $aluf $lr50v/$imr2 # PE番号が2,3なら足す
msl $lt $lt                    # 受け取ったのをさらに次の8Byteを担当するPEに渡してT regは上書き
ladd $lr50v $aluf $lr50v/$imr2 # PE番号が2, 3なら足して……
nop                            # $lr50vが利用可能になるまで待つ
lsub $lr50v $lt $lr50v/$imr3   # PE番号が2なら引く(合わせて、PE番号が3なら足したことになる)

# d get $lr50n0c0b0 4 # 通信完了!

※imm命令のuはunsignedではなくupperの意味っぽい。imm命令を短語書き込みで使えばよかったのかも?

マスクレジスタを使っている部分は、LMへのTレジスタ間接アクセス(3.6.1.6 m - LM0 (ベースアドレスレジスタ書き込みを除く))を使うともう少し短くできるのかもしれません。

他のアイデア

浮動小数点数演算は使わないからマニュアルを読まなくていいかなと思っていたんですが、どうも行列演算ユニットがそれなりに便利に使えるようなので、これをうまく使うと命令数が減らせるのかもしれません(整数命令ばかりやる場合は余っているので)。 まず、4つのPEではMAUを共有しているっぽい(明確な記載はないが行列レジスタ書き込み命令や行列ベクトル積和演算命令はそのように動作するみたい:3.6.10 行列レジスタ書き込み命令式・3.6.9.2 dmfma - 倍精度行列ベクトル積和演算の基本動作)ので、適切な行列を事前に入れておけば、これを無理やり使ったシャッフルやprefix sumが可能なのかもしれません。 また、行列レジスタ転置読み出しというのが存在して(3.6.1.14 x, y - 行列レジスタ。命令名は行列レジスタ読み出し命令っぽい:3.6.11 行列レジスタ読み出し命令式。非常にわかりづらいが、どうも書いたデータをそのまま読み出すことはできなくて転置読み出ししかできないみたい)、これを使うと他のPEが書いたデータが手に入るっぽいので、これを使ったデータ交換も可能なのかもしれません。

水平方向命令は存在しないのかなと思っていましたが、ブロックフロート化命令(3.6.12.11 bfe, bfn - ブロックフロート化命令)は横方向のデータと相互作用のあるSIMD命令なので、これを使うと語長を読み替えてマスクレジスタを作って……みたいなことをしなくてもうまくできる可能性があります。 うまい感じに奇数バイト目が重要で半精度浮動小数点数の指数部になっているのも、いいかんじです。 でもこれってもしかしてPE内で閉じてないで4つのPEが相互作用する命令だったりするのでしょうか? 明確な記述がなくてわかりませんが、倍精度用の命令が1長語読み出して1長語出力する命令であるからには何らかの作用があるはず(相互作用がなければ恒等変換なので無意味)、と思うと多分そうなっているのでしょう。 幸い、この問題の入力は32Byte中にも二つ以上の改行文字があることはないので、それでも使えそうです。 PE方向に縮約するL1BM命令は存在しないので、これは便利かもしれません。

上で書いた後気づいたんですが、ゼロフラッシュマスク(3.6.2.2 ゼロフラッシュマスク適用)が存在するのでゼロ初期化はいらなかったようです。

あと気づいたんですが、ゴルフ的には$nowriteよりも$tで捨てたほうが短くなりますね。 また、$ltなどとは書かなくていいようです(3.6.1.12 t - T レジスタ)。

anarchy golf - hello world

くやしいのでHello, world!のゴルフをやります。 一応endlessの問題ですが、shinhさんも公開していますし、そもそもそんなぎりぎりまで詰めた解でもなさそうなので(Preferred Networks社内でMN-Core機械語を日夜(?)書いている方々に勝てるわけがありません)、参加者が増えることを期待して、公開することにします。

とりあえずshinhさんの解をもとにして、以下のような感じにしてみました。

sor $subpeid $r0 $omr1
imm i"0x6f726c64" $r0
imm f"0x1p-61" $r1
imm f"234929.69" $r0/$imr1
imm i"0x6f2c2077" $r1v/$imr1
nop
l1bmd $lr0 $lb0
nop/3
l2bmd $lb0 $lc0
nop
mvd/n64 $lc0 $p0@0

結構怪しいコードができました。 まず一行目では、エミュレータが初期化時にゼロクリアしていると思われる(実機では何が入っているか保証されなさそうな)$r0を読み出すことで、即値0を作る手間を省いています。 これを仮定出来ない場合、zero $tが一命令かつ文字数も少なくてよさそうです。 基本的にはここで作ったPE=0の時だけ書き込まれないマスクを使って出力文字列を作っていくのですが、いったん書いてからlpassaで上書きするのではなく、imm命令の出力をそのままマスク書き込みする感じにしました。 あと、即値のu指定はいらなかったようです。これは何に使うんだろう?

また、完全に無意味なゴルフテクニックですが、0x21000000は縮められます。 0x21000000は10文字ですが553648128とすれば9文字です。 さらに、今回は末尾が全部0なので浮動小数点数としてみなせば仮数部が空っぽということで、浮動小数点数リテラルを使えば短く書くことができます。 絶対値がかなり小さい数なので十進数で書くとかえって長くなりますが、十六進浮動小数点数リテラルがうまくはまり、7文字にすることができます。 234929.69も同様で、10文字が9文字に減らせます。仮数部が24bitなので十進八桁で指定できて、残りの指数部を小数点の位置でエンコードしているので縮まっているようです。

$r1vとしているのは、上で話したハザードをよけるテクニックです。 これにより、命令数が1つ減ります。 また、ゴルフ的にも、vで1文字増えますがnop/2nopになって2文字減るので、合計で1文字縮まります。

PEからL1Bに送り出すのは、ニーモニックが短いl1bmd命令を使うことにしました。 また、語長も短くていいので、$lr0とすることで、lを削減できます。

最後のMV命令は、mvd命令を使うと@.0@0にできて1文字削減できます。

また、少なくともcaddyの環境では、最後のウェイトは特に必要ないようです。 これは、caddyが後ろにつけるnop; wait i02で待ってくれることが原因かもしれません。

このソースコードは、末尾に改行コードがないと正しく動作しません。 これは、caddyが後ろにつけるfと連結されてしまうからです。 caddyが後ろにつけているcefといった命令(dのようなデバッグコマンドの一種?)の説明はソフトウェア開発者マニュアルにはなかったので、これが何なのかはわかりませんでした。

これにより、185Byte(28/127/30)、11命令の解答を作ることができました。

比較的ゴルフ出来そうなPost mortemな問題

基本的には、以下の条件を一つでも満たすと、MN-Core 2の機械語コードで出力を作るのはかなり面倒そうです。

  • (一入力が「複数ケース」的な場合)ケースごとに入力長が異なる
  • 出力の各行の長さが異なる
  • 出力の各行の長さが同じでも、その長さが入力に応じて変化する(二冪で変化する場合は除く)

以下では、そのような条件を満たさない、比較的MN-Core 2の機械語コードに優しそうな問題を集めてみました。 だれか挑戦してみてください。

anarchy golf - Permutations

Post mortemの中で、出力を並列に計算できそうな問題です。 この問題は入力がなく、計算により出力を構築する問題です。 一行12文字(改行文字含む)の出力を720行すれば正答できるので、比較的MN-Core 2の機械語で解くのが容易そうです。 ただ、64と比較的大きな二冪の倍数ですが、全体としては二冪ではないため、それなりに面倒そうです。

anarchy golf - Missing Digit

Post mortemの中で、出力を並列に計算できそうな問題です。 この問題は入力から出力を計算する問題です。 一行あたり10文字(改行文字含む)を読み2文字(改行文字含む)出力する問題なので、比較的MN-Core 2の機械語で解くのが容易そうです。 処理の単位が10Byteで二冪ではないため、それなりに面倒そうです。

anarchy golf - Digit depth

Post mortemの中で、出力を並列に計算できそうな問題です。 この問題は入力長と出力長が同じ、文字列変換の問題です。 括弧のネスト段数に応じて出力文字を変えるという問題なので、prefix sumアルゴリズム(対数段の通信)で解くことができます。 MN-Core 2のアセンブリを直書きして性能を引き出すために習得しないといけない技能の勉強にはもってこいという感じがします。 ただ、ゴルフ的には、問題サイズが小さいので直列にやった方が縮まりそうな気はします……。

もうちょっと面倒そうな問題

先ほど挙げた条件を多少満たしているものの、なんとかなりそうな問題を集めてみました。

anarchy golf - The 9801 debacle

ぶっちゃけ数を0から9999まで出力するみたいな問題です。 leading zeros付きで出力すればよいので出力単位もそろっていて並列計算可能です。 むずかしいのは入力に応じて出力の単位が8Byte/6Byte/4Byteと変化する部分ですが、二冪が2つ含まれる3通りなので比較的なんとかなるかもしれません。

anarchy golf - Morse decode

モールス信号を解読する問題です。 “袖”付きで入力をPEに持ってくれば、出力を計算すること自体は可能です。 むずかしいのは入力の何バイト目と出力の何バイト目が対応しているのか自明でないことですが、これもprefix sumとかでなんとかするのでしょう。

まとめ

参考文献

[1] J. Makino, K. Hiraki, M. Inaba, Proc. of the 2007 ACM/IEEE Conf. on Supercomput. 21(18), 1-11.
[2] K. Namura, et al., 2021 Symposium on VLSI Circuits, 35(89), 1-2.

*1:「MN-CoreはGRAPE-DRの拡張である」とMN-Coreの論文[2]に記載されていました。

*2:競プロerへ:あなごるは高々三つ与えられるテストケースさえあっていれば良いという文化です。他の提出へのチャレンジ(いわゆるHack)はありません(乱数を利用した提出へのチャレンジとしてのリジャッジは有効化されている問題もあります)。そもそも問題文がものすごく曖昧だし……。

60bit整数を高速に素因数分解する

Library CheckerFactorizeという問題でfastestを取ったので、そのコードの解説です(速度ランキング)。 Factorizeは、 10^{18}までの整数が高々100個与えられ、それらを素因数分解せよ、という問題です。

概要

上位陣(13ms~15ms)の解法では、ECM楕円曲線を用いた素因数分解法)が使われているようですが、よくわからないのでこれには手を出しません。 ECMを用いていない提出の中で最速の提出はhqztrue氏の16msであり、20%以上の差があります。 この提出では、ポラードのロー法という手軽に実装できるアルゴリズムが用いられていました。 これを効率化してみます。

ポラードのロー法には、乱数を連鎖的に生成する部分が存在します。 この部分がレイテンシに支配されていて並列性に乏しいという問題点に着目しました。 これに対し、異なる二つの系列の乱数列を同時に生成する手法を導入します。 これにより、プログラムの並列度をスループットに支配されるところまで上げることができます。 並列化のオーバーヘッドがかなりあり、二倍速くなるというわけではないですが、それでも12msと上位陣の解法を上回る速度を達成できました。

基礎知識

純粋にアルゴリズムの基礎知識です。 CPUの特性を生かした高速化に関する記述はほぼないので、高速化テクニックが知りたい場合は「本題:私の提出コードの解説」まで読み飛ばしても多分大丈夫です。 ポラードのロー法を“貼って”いるけれど、何をやっているのかあまりよく知らない、という場合は読んでみると面白いかもしれません。

ポラードのロー法

ポラードのロー法は、合成数 nを入力すると、その非自明な因子( n 1 nでない約数) aを返す乱択アルゴリズム(ラスベガス法)です。 ここで、返ってくる a素数とは限らない点に注意してください。 とはいえ、素数判定と組み合わせ、素数でなければ再帰的にロー法を適用する、とすることで素因数分解を行うことができます。

乱数を使って素因数分解するというのはどういうことでしょうか。 単純な方法として、 [0, n)の範囲の整数を一様ランダムに生成し、それと nの最大公約数(greatest common divisor, GCD)を求める、という方法があります。 これを繰り返していれば、そのうち aが返ってくることがあるでしょう。 しかし、これには平均で aに比例した回数の乱数生成が必要です。

誕生日のパラドックスを考えると、次のような方法で乱数生成回数を減らせます。  [0, n)の範囲の整数をランダムに \sqrt a個程度生成し、それら同士の差を計算します。 ここで計算された「差」一つ一つに対して、 nとの間の最大公約数を求めれば、乱数生成は \sqrt aに比例した回数で済みます。 しかし、乱数を二つ取り出す組み合わせの数を考えると、「差」は a/2個ほどあるので、計算量はやはり aに比例してしまいます。 なんとか無意味な組み合わせの「差」を求めることを回避できないでしょうか。

これに対してはうまい方法があり、それがポラードのロー法です。 ここまでは一様ランダムな真の乱数を考えていましたが、そうではなく、あえて疑似乱数を使うというものです。 例えば、なんらかの整数 Cに依存した疑似乱数系列 x_i \mapsto x_i^2 + C \bmod nがよく用いられます。‘

「差」と nの間の最大公約数を求めたときに aが返ってくるのは、「差」を取った二つの乱数が \bmod aで一致していた場合なので、乱数の aでの剰余に着目します。 いきなり aでの剰余に着目してもいいのですが、それだと当たり前のことを言っているだけのようにも感じてしまうかもしれません。 そこで、あえて任意の Kでの剰余の場合にどうなるかを考え、その後 K = aという特殊な状況ではどのように性質が変わるかを見てみます。

さて、この疑似乱数系列を \mod\! Kで見た y_i = x_i \bmod Kを考えます。 この時、 y_i \mapsto y_{i+1}という変換は一意に決まる(写像になる)でしょうか?  Kがなんでもない数であれば、これは一意に決まるとは限りません。  x_i y_iよりも多くの情報を持っているので、 y_iが決まっても y_{i+1}には複数の可能性があります。

例で考えてみます。  x_i \mapsto x_i^2 + 1\bmod 15という疑似乱数生成式は、  x_1 = 1に対して1→2→5→11→2→5→11→……という疑似乱数系列を生み出します。 ここで、 K = 6で割った余りを考えると、1→2→5→5→2→5→5→……となっています。 つまり、5の次が5だったり2だったりしています。  y_i \mapsto y_{i+1}という変換が一意でないということです。

しかし、 K nの約数の場合は状況が違います。 たしかに x_i y_iより多くの情報を持っていますが、それは y_i y_{i+1}の関係には影響を及ぼしません。 なぜなら、 y_i \mapsto y_i^2 + C \bmod Kが成り立つからです( \bmod rsで成り立つ式は \bmod rでも成り立つ、という性質を使いました)。 よって、  K nの約数の場合は変換が一意です。 変換が一意なので、系列中に同じ値が出現すると、以降は全く同じ系列が繰り返されることになります(上記例では5→5となったからといって、その後も5→5→5→5→……となるわけではないのと対比してみてください)。 以下では、系列が繰り返される周期を Pと表記します。

話を戻すと、乱数 x_iと乱数 x_jの差 x_i - x_jを求め nとの最大公約数を求める手順で、無意味な組み合わせを回避したいという話でした。 無意味な組み合わせとは、 x_i - x_j aの倍数にならないことが明らかな組み合わせのことです。 これは、 K=aで考えれば、 y_i y_jが必ず異なることが明らかな組み合わせのことです。 ポラードのロー法では、疑似乱数系列に周期があることをうまく活用し、無駄な組み合わせを回避します。 例えば、 y_{i+1} y_{i+3}が異なっていれば、 y_i y_{i+2}も異なっています(対偶を考えれば明らかです)。 よって、 x_{i+1}-x_{i+3} nの最大公約数が1だった場合、 x_i-x_{i+2} nの最大公約数を求める必要はありません。

この考え方、つまり添え字の差が同じペアは調べても無駄、という事実を活用すると、調べるべきペアの個数を P個程度に減らすことができます。  x_i x_{i+1}の差と nの最大公約数を求める、  x_i x_{i+2}の差、 x_i x_{i+3}の差、 x_i x_{i+4}の差、……とP回繰り返せばよいからです(※)。 乱数を P個生成した時、ペアは P^2/2個程度ありますが、調べる必要があるのは P個程度でしかない、ということです。

周期 P K=a以下*1であることが明らかですが、 y_{i+1} = x_{i+1} \bmod K = y_i^2 + C \bmod Kという疑似乱数生成式を使った場合、周期 P \sqrt K=\sqrt a程度になることが経験的に知られています。  nの任意の因子 cについて、 \sqrt c個程度のペアを調べると cが見つかるということなので、基本的には小さい因子から見つかっていきます。  nの最大の素因子はこの方法で見つけなくていい( nを見つかった因子で割るということを繰り返して残った数が素数かを判定すればよい)ので、 nの二番目に大きな素因子の大きさの正の平方根に比例するくらいの時間がかかるとみてよいでしょう。  nの二番目に大きな素因子の大きさは高々 \sqrt nなので、 n^\frac14に比例した回数の乱数生成により素因数分解できることが経験的に知られています。

※実際には、 iをどう選ぶかという問題があります。  y_1が循環に含まれるとは限らないからです(例えば、上で紹介した1→2→5→11→2→5→11→……の例では、1は循環に含まれていません)。 これは、適宜 i自体を進めていく、という方法で解決可能です。 ポラードが提案した進め方(二回に一回進める)よりも乱数生成回数が少なくなる進め方( 2^t回に一回、系列を 2^{t-1}だけ進める)をブレントが提案しています。

参考

以下の記事に、アニメーションを用いたわかりやすい解説があります。

ポラード・ロー素因数分解法について #競技プログラミング - Qiita

モンゴメリ乗算

 \mod\!nの計算は、普通にやると除算命令が必要です。 除算命令は一般に遅いのでどうにか回避したいです。

モンゴメリ乗算は、バレット還元と同様に、いくらかの事前計算を行うことで、剰余演算を除算命令無しに行う方法です。 モンゴメリ乗算では、モンゴメリ表現と呼ばれる数値表現で値を保持したまま計算が進みます。 モンゴメリ表現は、 \bmod nの代表元 Xの代わりに、 2^{64}などの切りのいい数を掛けた 2^{64}X \bmod nを使います。

事前計算で求めるものは、モンゴメリ表現から通常の表現に戻すときに使う値 Rです(他にも用途ごとに他の数を事前計算することがあります)。 これは \bmod 2^{64}での nでの逆元であり、  nR \equiv 1 \mod 2^{64}を解けば求まります。 この Rを用いると、乗算を以下のように行うことができます。

元の表現 モンゴメリ表現
 X  A = 2^{64}X \bmod n
 Y  B = 2^{64}Y \bmod n
 XY  C = 2^{64}XY\bmod n
 \;\;\;\:= 2^{-64}AB\bmod n

さて、この 2^{-64}AB\bmod nはどのように求めたらよいでしょうか? ここで \bmod nが発生していたら意味がありません。 これを事前計算された Rを使って回避するのが、モンゴメリ乗算です。  Z \equiv 2^{-64}D \mod nとなる Zは次のように計算できます。

 Z = \frac{D}{2^{64}} - \frac{(DR \bmod 2^{64})n}{2^{64}}

この手順には 2^{64}での除算や 2^{64}での剰余を求める部分がありますが、64bitCPUであれば実質ノーコストである点に注意します。

なぜこれでいいのでしょうか?まず、 Zが整数であることから確認します。

 D - (DR \bmod 2^{64})n \equiv D - DRn \equiv D - D\times1 \equiv 0 \mod 2^{64}

なので、 D - (DR \bmod 2^{64})n 2^{64}の倍数であり、 2^{64}で割っても整数です。 また、  D (DR \bmod 2^{64})n 2^{64}の倍数とは限りませんが、その差が 2^{64}の倍数であるということは、下位64bitは同じになっています。 よって、もし \frac{D}{2^{64}} \frac{(DR \bmod 2^{64})n}{2^{64}}を切り捨て除算で計算しても、結果は同じになります。  Z = \frac{D - (DR \bmod 2^{64})n}{2^{64}}のように計算すると D - (DR \bmod 2^{64})nの下位64bitの計算が一見必要そう(下位64bitは除算で消えるにせよ、ボローが生じるかを計算しないといけなさそう)ですが、実はそれは不要、ということを意味しています。

つぎに、 Z \equiv 2^{-64}D\mod nであることを確認します。

 Z = 2^{-64}D - 2^{-64}(DR \bmod 2^{64})n \equiv 2^{-64}D \mod n

大丈夫そうです。

このように、 2^{64}などの切りのいい数を掛けた表現で持つことで、多少の事前計算と引き換えに剰余演算を除去するのが、モンゴメリ乗算です。

なお、単に 2^{64}を掛けた表現で持っているだけなので、加減算は普通の剰余環でのものと同じようにして実行可能です。

ポラードのロー法では非常に多くの剰余演算が使われるため、高速化のためには事実上必須のテクニックです。 以下では、ポラードのロー法といったとき、暗黙でこのテクニックが適用されているものとします。

素数判定

ポラードのロー法で素因数分解するためには、素数であるかを(素因数分解によらずに)判定する必要があります。 いろいろな高速な手法が提案されていますが、素因数分解にかかる時間が圧倒的なため、素数判定の高速化に躍起になる必要はありません。 今回はミラーラビン素数判定法を愚直に実装したものを使いました。 ミラーラビン素数判定法は確率的判定法ですが、 2^{64}までの整数の素数判定に対しては、特定の7つの底を使えば間違いなく判定できることが知られています。

Binary GCD

Binary GCD(Stein's algorithm)は、最大公約数を求めるアルゴリズムです。 二進法を採用している計算機(今日のほぼすべての計算機)に適しており、高速に動作します。

最大公約数を求めるアルゴリズムとして有名なユークリッドの互除法は、剰余演算が含まれるため、あまり高速ではありません。 アルゴリズムの進行とともに入力依存で法が変化するため、剰余演算を事前計算で高速化することも不可能です。 Binary GCDは除算を含まないため、高速に動作します。

以前の記事(最大公約数をもっと高速に求める番外編(その2) - よーる)では、C言語コードを工夫して書けばclangが最適なコードを吐くようにできることを紹介しました(インラインアセンブリに手を出さなくて済むということです)。 当時はgccだと最適なコードとならない問題がありましたが、現在のgcc(ジャッジサーバーで使われているのはgcc 13.2)ではほぼ最適なコードが出るようになっているようです。

60bit整数同士の最大公約数を求める場合を実験してみると、このコードのループは平均で40回程度実行されるようです。 Binary GCD (GNU MP 6.1.2)によれば1bitあたり0.68回の反復が必要、と書かれているのでこれと合致します。 ループ一周は5cycleで実行可能*2なため、最大公約数を求めるのにかかる時間は200cycle程度です。

ポラードのロー法で、最大公約数をまとめて計算するテクニック

ポラードのロー法を愚直に実装すると、「乱数を生成し、差を計算し、それと nの最大公約数を求める」というのを繰り返すコードになると思います。 しかし、これはかなり遅くなります。 先ほど示したように、最大公約数を求めるのにかかる時間は200cycle程度です。 これは乱数の生成にかかる時間が10cycleや20cycleであるのに対して10倍も大きいです。 これではアルゴリズムの実行時間の大半が最大公約数を求めることに費やされてしまいます。

これには面白い解決策があります。  nとの間の最大公約数を求める数を、1つの差とするのではなく、 M個の差を全部掛け合わせたもの(の nでの剰余)とするのです。 これにより、最大公約数を求めるのにかかる時間を 1/Mにすることができます。 ただし、全部掛け合わせると \bmod n 0になってしまう場合(≒複数の因数が M個の中で同時に発見された場合)の対処が追加で必要になります。

この手法は非常に大きな性能向上につながるため、高速化のためには事実上必須のテクニックです。 以下では、ポラードのロー法といったとき、暗黙でこのテクニックが適用されているものとします。

本題:私の提出コードの解説

ボトルネックの解析

ポラードのロー法のボトルネックは、疑似乱数を順次作っていく部分にあります。  x_iモンゴメリ表現から x_{i+1} = x_i^2 + C \bmod nモンゴメリ表現を求めるには、モンゴメリ乗算が1回必要です。 モンゴメリ乗算には、以下が必要です。

  • 64bit×64bit=128bitの乗算3回
    • フルの128bitが必要なのが1回
    • 下位64bitだけ得られれば良いのが1回
    • 上位64bitだけ得られれば良いのが1回

これらの乗算はいずれも、IntelSandy Bridge以降のプロセッサでは3cycleで実行可能です。 これらの乗算は直列に依存しているため、これだけで9cycleが必要です。 さらに一回の加減算が必須で、あわせて10cycleのレイテンシが避けがたいです。

差を求めて掛け合わせる部分でもモンゴメリ乗算が1回必要なので、合計で一つの乱数生成当たり6乗算必要(ブレントの改良方式の場合)ということになります。 乗算器はパイプライン化されているので、10cycleの間にスループットとして10回の乗算が可能であり、たった6回しか乗算を行わないのは非常にもったいないです。 このように、疑似乱数を順次作っていく部分がクリティカルパスとなっており、並列に実行できる命令が十分にないため、CPUはかなりの時間遊んでいます(とはいっても平均IPCが2を超えるような命令レベル並列性が豊富なプログラムなのですが、最近のCPUの最大IPCは4以上あるので、それには届かないということです)。

ボトルネックの解決策

乗算器が余っていることに着目し、並列に二つの疑似乱数系列を生成していくようにすることで、高速化が期待できます。 このようにすると、ループ中に12乗算含まれる形になります。 IntelのCPUは1cycleにスループットとして1回の乗算が可能であり、このループの逆数スループットは12cycleとなります。 これは疑似乱数生成のレイテンシは10cycleよりも長いので、理想的にはスループットに支配されたプログラムとなります。 ただし、逆数スループット12cycleはクリティカルパスレイテンシ10cycleとあまり差がないため、あまりにも雑な命令スケジューリングが行われるとレイテンシが見えるようになってしまう点に注意する必要があります。 普通のCPUの命令スケジューリングでは古い命令が優先的に選択されることを考えると、掛け合わせる部分の命令列よりも十分前に疑似乱数を生成する部分の命令列があれば良さそうです。 この目的のためには、ソフトウェアパイプライン化を施せばよいです。

実は、以下の条件を満たす場合はソフトウェアパイプライン化を施さなくてもそこそこの命令列が生成されていました。

  • clangでコンパイルする
  • 掛け合わせる部分のコードが実行されるのが、乱数生成直後ではない(例えば、差を計算するのに時間がかかるなど)

gccコンパイルすると遅いのはなんでだろうとずっと悩んでいたのですが、clangを使っていても差を計算する部分を軽量化(後述)したときにかえって遅くなったことから、命令スケジューリングに問題があることに気づくことができました。 ソフトウェアパイプライン化を施したコードは、上記条件を満たすときのclangの出力コードよりも5%程度高速であり、これが単独fastestを取る詰めの一手となったようです。

(あまり効果のない)軽量化

モンゴメリリダクションしない

普通の感覚では、差を掛け合わせた数と nの最大公約数を求めるとき、まずは差を掛け合わせた数のモンゴメリ表現qに対し元の表現Qを求め(この手順をモンゴメリリダクションといいます)、その後Qnの最大公約数を求めると思います。 しかし、qQ 2^{64}倍した数を \bmod nで考えているだけなので、nの奇数の約数aについて、Qaの倍数であればqaの倍数です。 よって、元の表現Qを求めなくても、qnの最大公約数を求めれば目的を達成できます。

正規化しない引き算

普通の感覚では、モンゴメリ表現xモンゴメリ表現zの差を求める際、x < z ? x - z + mod : x - zなどと求めるはずです。 しかし、モンゴメリ表現は [0, n)の範囲に正規化されている必要はないため、x < zかにかかわらず常にx - z + modとしても大丈夫です。 また、xはかなり長い間固定なので、x + modを計算しておけば、差を求めるのはsub命令一つでできるようになります。

なお、これにより掛け合わせる部分のコードが実行されるのが乱数生成直後になってしまいかえって遅くなってしまったことが、命令スケジューリングに問題があることに気づけたきっかけでした。

謎の数のモンゴメリ表現を使う

普通の感覚では、 x_{i+1} = x_i^2 + C \bmod nの計算にモンゴメリ乗算を使う場合、 Cモンゴメリ表現を求める必要があるように思われます。 しかし、 Cの具体的な値は重要ではなく、 Cには任意の整数が使えます。 よって、「モンゴメリ表現が1になるような謎の定数 Cを使った疑似乱数生成器」を使っているんだと考えることで、モンゴメリ表現に変換するコードを記述しなくて済むようになります。 この定数 C nに依存するためアルゴリズムが正しく動作するのか怪しいところがありますが、少なくとも 2^{40}までの合成数の非自明な因数を見つけることは可能なようです(Intel Core i9-12900Kにて24スレッド並列で1日くらいかけて確認)。

この技術は、あらゆるところに使われています。 乱数のシードZ1, Z2モンゴメリ表現である必要がありますが、いかにも普通の表現のような顔をして1などと書かれています。 掛け合わせていくときの初期値q1, q2は、普通は1のモンゴメリ表現を使いますが、ここでもわざと普通の表現と混同して1などと書いています。

これらは、プログラム中の生存変数を減らすことを目的に行っています。 ただ、頑張って減らしても結局gccが発生させるスピルを取り除くことはできませんでした(clangだと生存変数が多くてもスピルも無駄なmov命令もほぼなかったのに……)。 最初はgccの生成するコードが遅い原因をスピルなのではないかと疑っていた(命令スケジューリングに問題があることに気づいていなかった)ため、生存変数を減らす努力をしていたのですが、結構無意味な努力だったようです。

疑似乱数生成器を固定する

一般に一つの疑似乱数生成器だけでは素因数分解ができない合成数があるため、疑似乱数生成器を[n]( uint64_t x ) { return __uint128_t(x) * x % n + 1; }などとハードコーディングするのは正しくない(無限ループしうる)ことを以前に指摘しました(間違ったポラードのロー法の使い方 - よーる)。 しかし、今回は二つの疑似乱数生成器を並列で使っているため、このような問題は多分起こりません(確証なし)。 よって、プログラムの命令数を減らすために、疑似乱数生成器をハードコーディングすることができそうです。 プログラムの命令数を減らすことが高速化につながっているのかはよくわかっていません。

ポラードの方法(フロイドの循環検出法)でもブレントの方法でもない循環検出法

そんな大層なものではなく、ブレントの方法の亜種です。

周期を Pとして、循環を検出する( nとの最大公約数が1でない「差」を見てくる)のに必要な乱数生成回数は、フロイドの循環検出法では 3P回程度、ブレントの改良方式では 2P回程度です(正確には周期だけでなく循環に達するまでの反復回数(ρのしっぽの部分の長さ)を考えなければなりません。これに関するわかりやすい説明が前掲のポラード・ロー素因数分解法について #競技プログラミング - Qiitaに載っています)。 このことからブレントの改良方式が優れているとの考えが主流ですが、これは必ずしも正しくありません。 二つの手法における計算量を確認してみます。

フロイドの循環検出法では、一周当たり20cycle(乱数生成を直列に2回やるレイテンシ。もう一系列の乱数生成のレイテンシは並列実行で隠蔽可能)で乗算が12回(3乗算×乱数生成3回+3乗算×掛け合わせ1回)です。 循環検出までにかかる周回数は P回程度です。

ブレントの改良方式では、一周当たり10cycle(乱数生成を1回やるレイテンシ)で乗算が6回(3乗算×乱数生成1回+3乗算×掛け合わせ1回)です。 巡回検出までにかかる周回数は 2P回程度です。

このようにしてみると、どちらが特に優れているなどとは言えなさそうに見えてくると思います。

ブレントの改良はもう一つあり、差を取る乱数の片側(提出コードにおけるx側)を変更した直後は掛け合わせを省くという最適化です。 前掲のポラード・ロー素因数分解法について #競技プログラミング - Qiitaによれば、

もし r/2 < k < 3r/4 の間に約数が検出されるとすると、残りの 3r/4 ≦ k < r の間でも再度同じ約数が検出されるため

という説明がされています。 これはアルゴリズムの正当性(および反復回数を上から抑える議論)としては正しいものですが、近年のCPUにおいて高速化することの説明としては正しくなさそうです。 ブレントの改良方式では、一周当たり10cycleで乗算が6回なので、ここから掛け合わせを省いて乗算を3回に減らしても特に高速化しないからです(クリティカルパス上に無い命令が減ってもレイテンシに支配されているプログラムは高速化しません)。 逆にこの“スキップ”により、せっかく乱数を生成しているにもかかわらず因数を見つけるチャンスを失っていることになります。 この機会損失と最大公約数計算の手間のどちらが大きいかは簡単には分かりませんが、単純に「仕事量が減っているから高速化するはず」という雑な議論では済まなさそうです。

私の導入した二つの疑似乱数系列の並列生成を導入した場合、以下のようになります。

  • 掛け合わせを省かない場合、一周当たり12乗算なのでスループットで支配されて12cycleかかる。よって一周は12cycle
  • 掛け合わせを省く場合、一周当たり6乗算なのでスループットは十分で、レイテンシの10cycleがみえる。よって一周は10cycle

つまり、掛け合わせを省くことで多少は高速化するようです。 しかし、実験してみるとこのような先ほどよりも良い条件にも関わらず、素因数分解の速度の向上は見られませんでした。 やはり、“スキップ”中の因子発見チャンスを捨てている損失がある程度存在するようです。

フロイドの循環検出法を用いたものも実験してみましたが、ほぼ同じかやや遅いくらいの性能になりました。 原因は正確には調査できていないですが、どうも「全部掛け合わせると \bmod n 0になってしまう場合」という運の悪いパターンの発生率が上がっているようでした。 ソフトウェアパイプライン化がうまくいっていないなどの原因も考えられます。

以上の実験により、提出コードは“スキップ”を省いたブレントの改良方式を使っています。

パラメータの調整

差を512個掛け合わせてから最大公約数を求める

差を掛け合わせる個数 Mを大きくすることには、以下のような得失があります。

  • 最大公約数を求める回数が減るため、高速化する
  • 「全部掛け合わせると \bmod n 0になってしまう場合」の発生率が上がる(発生するとかなり低速になる)

実験を行った結果、 M=512程度を選ぶと実行時間が最小化されるという結論に達しました。

最悪実行時間に最適化

ブレントの改良方式では、差を取る乱数のうち固定側(提出コードにおけるx側)を更新するのは乱数生成の1, 2, 4, 8, ...回目にするのが普通です。 しかしこれは増やしていくという部分が重要なのであって、最初は1回で更新、などとする必要はありません。 逆にこのようにすると、頻繁に最大公約数を求める手順が挟まるため、大きな合成数素因数分解では不利そうです。 よって、提出コードでは、乱数生成の512, 1024, 2048, 4096, ...回目に更新することにしました。 これにより、小さい合成数に対する速度がかなり低下しますが、(少なくともこの問題に対しては)小さい合成数に対する速度が上がっても得にならないので、このような最適化を行いました。 そもそも、小さい合成数素因数分解を速く行いたいのなら、ロー法ではなく試し割り法の方が向いているという事実もあります。

その他

決定的素数判定法のコードに関する疑惑

本題と関係ありませんが、ミラーラビン素数判定法でif( n <= a ) { return true; }とやるのはまずいという話を見つけました(ミラーラビン素数判定法についてメモのコメント)。 これは一般には多分そうなのですが、提出コードで使われている64bit整数に対する決定的素数判定法(競技プログラミングでよく使われる方法)では問題ないはずです。 まず、n == 2の場合は素数、そうでなくn % 2 == 0の場合は素数ないと判定する部分は問題ありません。 その後、まずはa = 2でミラーラビン素数判定が行われます。 底2での強擬素数(ミラーラビン素数判定が誤判定する合成数)は最小でも2047なので、このテストが通っている時点でif( n <= 7 ) { return true; }if( n <= 61 ) { return true; }は正当です。 七つの底を使う方では、そもそもそちらが実行される条件がn >= 4'759'123'141なので、if( n <= 1'795'265'022 ) { return true; }などの条件が成立することは決してありません。

プログラムの実行時間の比率と乗算器利用効率

おおよそ以下のような感じでした。

  • 素数判定1.5%
  • 最大公約数を求める部分4.3%
  • 最大公約数を求め終わった時の分岐予測ミスペナルティ0.4%
  • 乱数生成部分93.1%
  • その他0.7%

乱数生成部分は一周あたり12cycleで512回繰り返されるので6144cycle、最大公約数の計算は200cycle程度なので、比率のオーダー感としてはあっています。 分岐予測ミスペナルティは20cycleくらいなので、最大公約数の計算にかかる時間との比率もあっています。 乱数生成部分で乗算器を使った回数をカウントし、それを実行にかかった時間で割って得た乗算器利用効率は、86.3%でした。 乱数生成部分が占める時間は93.1%なので、92.7%ほどは乗算器を有効利用できているということになりそうです。

なぜ二並列なのに二倍速くならないのか

まず、二並列にした時点でレイテンシ支配ではなくスループット支配に変わっています。 レイテンシは一周当たり10cycle、逆数スループットは12cycle(つまり一系列当たり6cycle)なので、この時点で1.67倍までしか速くならないことになります。

また、ロー法で循環を検出するまでにかかる平均周回数が半分になるわけではないのも原因のようです。 ロー法で循環を検出するまでにかかる周回数が指数分布に従っているとすれば、(仕事を並列化するのではなく)二系列同時にやって早く巡回を検出したほうを採用、という手順で平均周回数が半分になって、二倍高速化するはずです。 しかし、ロー法で巡回を検出するまでにかかる周回数は指数分布から外れているらしく、二系列同時にやっても平均周回数は半分にはならず、2/3程度にしかならないようです。 分布の標準偏差は平均値の半分程度(指数分布なら平均値と同じくらいのはず)であり、短い巡回が見つかる可能性がかなり低いことを示しています。 非巡回部分(ρのしっぽの部分)に原因があるのかもしれません。 また、歪度は1程度です(これは短い巡回には1という自明な下限があるにもかかわらず長い巡回には上限がないことが原因で、普通です。指数分布なら2あるはずなので、やはり短い巡回が少ないことを意味しているようです)。 尖度は1(四次中央モーメントで書くと4)程度であり、裾の重い分布であることを示しています。

これらから計算すると、提出コードは二並列にしているにもかかわらず実行時間は 2/3\times12/10 = 0.8倍程度、あるいは実行速度は1.25倍程度、にしかならないとわかります。

まとめ

  •  10^{18}までの整数(60bit整数)を素因数分解するプログラムを高速化し、Library Checkerにおいて最高速を達成した
  • ECM楕円曲線を用いた素因数分解法)を使うのではなく、ポラードのロー法をチューニングした、非常にわかりやすい実装
  • 乱数を連鎖的に生成する部分がレイテンシ支配になっていてもったいないことに着目
  • 二つの乱数系列を並列に生成することで高速化できた
  • レイテンシが見えてこないよう、ソフトウェアパイプラインを組む必要があった
  • ブレントによる改良型循環検出アルゴリズムは、別にそんなにすぐれているわけではない。また“スキップ”はやらない方がよさそうだった

*1: x_i^2の部分で正負の違いが無視されて同じ系列に合流するということは、巡回には正の数と負の数のどちらもが入ることはないということで、周期の最大は (a+1)/2な気がします。

*2:trailing zero countが3cycle、引き算命令とシフト命令が1cycleと仮定

Unityのgit管理をWSL2でやる

最近久しぶりにUnityでゲーム作りをしています。

共同で開発するためにgitを使うのですが、gitコマンドをWSL2から使用するといまいちな感じになったので、何とかする方法のメモです。

プロジェクトはWSL2の外に作る

UnityプロジェクトをWSL2の中(/home/lpha/以下など)に作るとうまくいかなかったので、Windows領域(/mnt/c/Users/lpha/Desktop/以下など)に置くことにしました。 いちいちcd /mnt/c/Users/lpha/Desktop/などと打ち込むのも面倒なので、cdsというエイリアスを作りました。

パーミッションの管理

WindowsにはUNIX系のようなファイルパーミッションの概念がないようで、この部分で問題が発生します。 新しく作ったファイルやUnityが自動生成したファイル(.metaファイルなど)のパーミッションは、VisualStudioやSourceTreeなどのWindows向けgitクライアントからコミットすると、644になります。 一方で、このリポジトリをWSL2のgitコマンドで見ると、コミット直後にもかかわらず変更があるように表示されます。 これは、Windows領域のファイルが755扱いで見えているからです。 つまり、gitの認識している状態(HEADではファイルのパーミッションは644)とWSL2の認識している状態(ファイルのパーミッションは755である)が異なるので、変更があるように表示されます。

また、git reset --hard HEADなどとしても、変更がある状態のままになってしまいます。 これは、WSL2はWindows領域のファイルのパーミッションを(デフォルトだと)記録せず、常に755であると認識するためです(記録されないので変更することもできません)。 これに対しては対処方法があります。 /etc/wsl.confに以下の記述を行うことで、Windows領域のファイルのパーミッションを記録することができるようになります(metadataの部分が重要)。

[automount]
options="metadata,umask=22,fmask=11"

なお、この記述を加えた後、WSL2を再起動する必要があります。 そのためには、コマンドプロンプトから以下のコマンドを実行します。

wsl --shutdown

パーミッションの管理2

このように設定しても、デフォルトのパーミッションが一貫していない点は変わりがありません。 なので、VisualStudioやSourceTreeなどからコミットした直後は、WSL2のgitコマンドで見ると変更がある状態になっています。 これは、Windows領域で新規作成されたファイルのパーミッションは何も記録されていないので、WSL2としては755であると判断するしかないためです。 しかし、git reset --hard HEADなどとすれば、ファイルのパーミッションが644に変更され、“正しい”状態になります。 これにより変更がある状態は解消され、以降のgitコマンドはいつも通り使える状態になります。

また、WSL2のgitコマンドでgit commitする場合は、ファイルのパーミッションが755になっていないか気を付けたほうが良いでしょう。 そうしないと、パーミッションが変更されるだけのコミットができたりして、コミット履歴が汚くなってしまいます。 これは多分pre-commitフックなどを使って対処するべき問題のような気がしますが、コミット作成はGUIツールの方が楽なのでgit commitコマンドを直に使用することがなく、放置しています。

メモ

どうもパーミッションを無視するという設定もできるようです。 ただ、私がWSL2からgitコマンドを使用したいという理由で共同開発者全員に強制するのも忍びなく、WSL2を使う側の負担のみで可能な方法を紹介しました。

なぜドルコスト平均法しないほうがいいか

新年になりました。

タイトルに関してですが、なんだか話があったので軽くまとめておきます。

要点としては、

  • 投資の期待値が標準偏差に比して大きい(インデクスファンドなどはこの条件を満たす)投資商品がある場合、
  • 余剰資金があるなら、
  • 何も考えずにその資金を即投入したほうが良い

というものです。

これに対して、「ドルコスト平均法は(高値で買ってしまう)リスクを削減しているから有効なはずだ」という反論がありますが、これが間違っていることを説明します。

なお、「投資の期待値が標準偏差に比して大きい」が満たされない場合、ドルコスト平均法は一定の合理性があります*1

直感的な理解

ドルコスト平均法を行うと、現金で持っている期間が存在するため、リターンを得る機会を逃します。 行いたかったことは「投資を行って、リターンを得る」だったはずなのに、「高値で買うリスクを回避」に目的が変わってしまっています。

その上、「高値で買うリスクを回避」と言っているにもかかわらず、値上がりしそうな投資商品を選んでいるため、待っている間に値上がりする蓋然性が高いです。 待てば待つほど、実のところ、高値で買うリスクを引き受けていることにすらなります。

統計的な理解

たしかにドルコスト平均法はリスクを削減します。 しかしその効果のほとんどは、「現金で持っている期間が長い」ことによります。 ドルコスト平均法では投資額面が平均で半分になるため、リスクを半分に下げる一方で、リターンも半分に下がってしまいます。

ドルコスト平均法特有の効果としては、

  • 投資商品の分散が大きくなるほど、高値で買うリスクを回避したことによる、追加のリスク削減が得られる
  • 一方で、投資商品の期待値が大きくなるほど、複利効果を逃したことによる、追加のリターン削減が発生する

のふたつがあげられます。

半額を投資する方法との比較

リスクを半分にしリターンを半分にするなら、余剰資金の半分を期初で一括で投資に回し、残りを現金で持っておく、という方法でも同じ効果が得られます。 これに対しドルコスト平均法は、

  • 分散がやや下がる
  • 期待リターンがやや下がる
  • あわせると、シャープレシオはわずかに上がる
  • リターンの中央値はほぼ同じ
  • リターンの第一四分位・第三四分位は下がる
  • 期初でのリスクが低い
  • 期末でのリスクが倍
    • 期末には全額を投資に回していることになるので

という違いがあります。

*1:正確にはこれは誤謬で、流動性が低くなっている時に取引をしている(安くなっている時には売り手が少なく、高くなっている時には買い手が少ない)ため売買が成立しにくい、という問題があります。

世界一簡単な物理レジスタ方式アウトオブオーダーCPUの作り方

この記事は自作CPU Advent Calendar 2023の24日目の記事です。

みなさんが使っているパソコンのCPUは、そのほとんどがアウトオブオーダー実行をしているはずです。 しかし、アウトオブオーダーCPUの作り方が書かれた教科書はあまりないような気がします。 トマスロのアルゴリズムというものが紹介されることもありますが、これは現代のアウトオブオーダー実行方式とはかなり違います。 この記事では、より洗練された現代的なアウトオブオーダー実行をするプロセッサを作ってみます。

全体的な構成

今回作るプロセッサは、FPGA上で動き、RV32IZmmul命令セットを動かすことができます。 東大のプロセッサ実験と同じルール(命令メモリ32KiB、データメモリ64KiB)で作ります(プロセッサ設計で一番難しいキャッシュを考えなくていいのですごく楽です)。

パイプラインは、以下のように設計しました。

  • 分岐予測ステージ
  • フェッチステージ
  • デコード+リネーム+ディスパッチステージ
  • スケジュールステージ
  • 発行ステージ
  • 実行ステージ
    • 普通の整数命令のレイテンシは1サイクル
    • 乗算命令のレイテンシは2サイクル
    • ロード命令のレイテンシは3サイクル
  • コミットステージ

今回は、作るのが大変なのでスーパースカラではないアウトオブオーダープロセッサを作ります。 この記事の作り方はスーパースカラであっても適用可能ですが、スーパースカラだといろいろ考えることが多くて大変です。

なお、スーパースカラではないのでIPCはどう頑張っても1.0を超えません。 普通にインオーダーのスーパースカラプロセッサを作ったほうが高いIPCを出せます。 性能が高いプロセッサを作ることが目的なら、こんな面倒なことをしてアウトオブオーダー実行を導入する必要はないです。 今回の話は、今後スーパースカラのアウトオブオーダープロセッサを作る土台作りということでやっています。

以下、各ステージを詳しく見ていきます。 なお、プロセッサが同時に取り扱える命令の数(リオーダーバッファサイズ)は16としました。物理レジスタ数は48です。

分岐予測ステージ

プロセッサの中で性能に最も影響を与えるステージです。

方向分岐予測器には、hashed perceptron分岐予測器を採用しました。 Hashed perceptron分岐予測器は、高性能分岐予測器の中では作るのが最も簡単です。 現在最強と考えられているのはTAGE分岐予測器ですが、hashed perceptron分岐予測器もそれとほとんど変わらないような性能を達成します。 TAGE分岐予測器は実装がとても大変なので、今回はhashed perceptron分岐予測器を使うことにします。

分岐先予測器は、分岐先バッファ(BTB)のみとしました。 置換とかのめんどくさいことを考えたくないので、32KiBの命令メモリに入る8192命令すべてに対して4Byteの飛び先を記録できるだけのBTB(32KiB)を用意しました。 こういう無駄遣いをすることで、アウトオブオーダープロセッサを簡単に作っていきます。 消費資源量の最適化は、正しく動いた後でやりましょう。

精度高く分岐予測するためには、方向分岐予測器で使うグローバル履歴をまだ結果がわかっていないにもかかわらず、予測に基づいて=投機的に更新する必要があります。 しかも、分岐予測ミスが起こった時にこの投機的更新を正しく取り消さないと、それ以降の分岐予測がめちゃくちゃになります。 これの対策は結構簡単で、グローバル履歴をリングバッファに持っておき、各命令でリングバッファのheadポインタを控えておけばよいです。 今回採用したhashed perceptron分岐予測器の最大分岐履歴長は90にしましたが、グローバル履歴のリングバッファサイズはきりよく128にしました。 これにより、38個までの投機的更新を行っても、それを取り消すことが可能です。 実際には高々16回しか投機的更新は行われないため、全然余裕です。

分岐予測ステージの結果は、

  • 次のサイクルの分岐予測の“種”として、分岐予測ステージへ
  • 次のサイクルのフェッチするアドレスとして、フェッチステージへ
  • 今フェッチした命令の「予測されたnext pc」(これは分岐を実行した時に分岐予測が正しいか判定するために必要)として、デコードステージへ

それぞれ送られます。

フェッチステージ

32KiBの命令メモリから4Byte命令をフェッチするだけで、特に工夫はありません。 命令キャッシュではないため、失敗はありません(常に命令をフェッチすることができます)。

デコード+リネーム+ディスパッチステージ

デコード

デコードするだけです。なにもむずかしいことはありません。

今回の作り方では、命令のレイテンシ情報、ソースレジスタを何個持つか、デスティネーションレジスタを持つか、だけをデコードしています。 発行ステージでデコードをする余裕があるため、ALUコードなどはそこでデコードをすることにして、4Byteの命令をそのまま後段に送っています。 このように作ることで、デバッグがやりやすくなっています。

リネーム

レジスタリネーミングを行います。つまり、命令に含まれる論理レジスタ番号を物理レジスタ番号に変換します。 アウトオブオーダー実行した結果を格納するために物理レジスタを使うのが、現代的なアウトオブオーダー実行方式です。

そもそも、アウトオブオーダー実行した結果をそのままアーキテクチャレジスタ(論理レジスタのこと)に書き込んでしまうと分岐予測ミスなどが起こった時に取り返しのつかないことになるので、アウトオブオーダー実行した結果はアーキテクチャレジスタとは別の場所に書いておく必要があります。 これの古典的な実装として、「アウトオブオーダー実行した結果はリオーダーバッファに書き込んでおき、分岐予測ミス等で取り消されないことが確定し次第順番にアーキテクチャレジスタに書き戻す」というものがあります。 しかし、この方式はデータの移動が発生してうれしくないです。 これより優れているのは、アーキテクチャレジスタとアウトオブオーダー実行の結果を書き込むバッファを区別しないで、一体の物理レジスタを使う、という方式です。 論理レジスタの値が物理レジスタのどこに入っているかは、そのポインタだけを記録しておけば、データ本体を動かす必要がなくなります。 「indirectは全てを解決する」というやつです。

ソース論理レジスタ番号を物理レジスタ番号に変換するのは簡単です。 その対応関係が記録された表(レジスタマップテーブル、RMT)を参照するだけです。

デスティネーション論理レジスタ番号は、未使用な物理レジスタ番号に変換される必要があります。 なぜなら、使用中の物理レジスタをデスティネーションとして割り当ててしまうと、他の命令が読もうとしていたデータが上書きされて別のデータに変わってしまい、プロセッサの動作がおかしくなる可能性があるからです。 ここで、「未使用」とか「使用中」とかは、「記録されているデータがもう二度と読まれない」「記録されているデータがまだ読まれる可能性がある」の意味です。

未使用な物理レジスタ番号は、フリーリストとよばれるリングバッファに記録しておくことにして、割り当てが必要になった場合はそこから一つ取り出します。 一般には未使用な物理レジスタが存在しない場合がありますが、今回は十分な物理レジスタを用意しました。 よって、この操作も失敗することがなく、簡単に実装できます。

後続の命令のリネームが正しくなるように、この命令に割り当てられた物理レジスタ番号をRMTに書き込んで更新します。 例えば、この命令のデスティネーション論理レジスタ番号が13だとして、それに物理レジスタ番号41が割り当てられたとします。 RMTにこのことを書いておけば、後続の命令でソース論理レジスタ番号が13の命令は物理レジスタ番号41を読みに行けばいい、と知ることができます。

この書き込みを行う前に、デスティネーション論理レジスタ番号に「今」対応する物理レジスタ番号も読み出しておきます。 これは、この命令がコミット(確定)された時に、永久に読み出されなくなる物理レジスタです。 よって、この命令をコミットするときに、その物理レジスタ番号をフリーリストに記録すれば、「未使用な物理レジスタ番号がフリーリストに入っている」状態を維持できます。 このようにして、物理レジスタガベージコレクション(といってもstd::unique_ptrと同様に自明ですが)が行われます。 物理レジスタ番号をフリーリストから一つ取得した命令は、そのうち物理レジスタ番号を必ず一つフリーリストへ返却するため、未使用な物理レジスタが枯渇することはありません。 RMTに書き込んでしまうと何を返却すべきかの情報が失われるので、この時点で読み出して保存しておく必要があります。

さて、アウトオブオーダープロセッサで最も面倒くさいのは、分岐予測ミスが起こった時にRMTをその分岐命令の時点の状況に戻さないといけない点です。 今回は贅沢に、すべての命令の時点でのRMTをまるまるコピーしておく(すべての命令でチェックポイントを取っておく)という方法を採用しました。 これはかなり贅沢で、RMTだけで〈6bit×32エントリ×16バージョン〉も必要になっています。 これは物理レジスタファイルに匹敵するサイズであり、しかも大量のコピーが発生しているので「実際のデータを移動させずポインタを変えるだけでいいから優れている」とか言ってたのを満たせていない気がしますが、気にしないことにします。

メモリ依存予測

本来であれば、リネーム(命令間のレジスタを介した依存の解析に相当)に加えて、メモリ依存予測(命令間のメモリを介した依存の解析)を行ったほうが良いです。 アウトオブオーダー実行をしている以上、プログラム順で後にあるロード命令が、プログラム順で前にあるストア命令よりも先に実行されることがあります。 この時、アドレスが一致していれば実行結果がおかしくなっている可能性が大です(メモリ順序違反)。 これを検出する機構は当然搭載されるため、間違った実行結果にはなりません(そういうことが起こったら分岐予測ミス同様にフラッシュしてやりなおします)。 とはいえ、頻繁にフラッシュが発生すると性能に影響があるので、「このロード命令はこのストア命令の後にやった方がいい」みたいな予測(これをメモリ依存予測と言います)をするわけです。

ただ、メモリ依存予測を実装するのはすごく大変なので、今回は実装しないことにしました。 実際の所、まともに最適化が行われていれば、ストアしてすぐロードすることはほぼありません。 したがって、命令ウィンドウ(命令を並べ替える範囲)が小さければ、メモリ順序違反はほぼ発生しないため、メモリ依存予測を行わなくてもそれほど性能が落ちることはありません。

また、メモリ順序違反を検出する機構であるロードストアキューのエントリの確保もこの辺のステージで行うのが普通ですが、今回は十分な数(16)のエントリを用意しました。 よって、確保が失敗することはなく、簡単です。

ディスパッチ

スケジューラ(どの命令が実行可能であるかを判断し、その中から適切なものを選ぶ機構)に命令を登録します。 また、スケジューリングに関係ない情報(命令のPCや即値など)を脇によけておきます(payload RAMに格納しておきます)。 これは、スケジュールステージは一番複雑なステージであり、なるべく小さく作りたいからです。

普通のスケジューラの作り方をすると、「この命令が実行された(のでこの命令に依存している命令は実行可能になった)」という信号は1サイクルしか流れません。 なので、

  • ディスパッチをする前に「すでにこの命令は実行されている(のでこの命令に依存している命令はとっくに実行可能である)」という情報を手に入れる必要がある
    • しかもこの情報を適切に更新する必要がある
  • ディスパッチをするその1サイクルに流れる「この命令が実行された(のでこの命令に依存している命令は実行可能になった)」という情報を逃さないようにする必要がある

という二点をちゃんと作らないといけないので大変です。

今回は、「どの物理レジスタが読み出せる状態になっているかを管理する」という贅沢な方式でスケジューラを作ることで、ディスパッチを大幅に簡単に作れるようにしました。 また、スケジューラに入れられる命令数を16と十分に大きくとったため、ディスパッチが失敗することもなく、やはり簡単に作ることができます。

このステージの結果(高々二つのソース物理レジスタ番号と高々一つのデスティネーション物理レジスタ番号と命令の実行レイテンシ情報)はスケジューラに書き込みます。 また、デスティネーション物理レジスタについて、「この物理レジスタはまだ読み出せる状態ではない」と記録します(表に0を書き込みます)。

スケジュールステージ

スケジューラ内にある高々16命令のうち、どの命令が実行可能であるか判断し、実行可能な命令の中から最も古い命令を一つ選びます。

スケジューラ内にある16の命令は各々、エントリのvalidビットに加え「自分のソース物理レジスタは読み出せる状態か」を二つ確認し、それらのandを取ることで実行可能状態(ready)であるかを判定します。 スケジューラはリングバッファになっていて、古い命令から順番に並んでいます。 よって16本のready信号を確認して、その中で最も古いものがどれかを調べることはそれほど難しくありません。 これにより、実行可能である最も古い命令を一つ選ぶことができます。

選ばれた命令は、自身のエントリのvalidビットを0にすることで、再び選択されることを防ぎます。 その後、選ばれた命令のレイテンシー1サイクルだけ待ってから(つまり、レイテンシが1サイクルの命令なら選ばれたサイクルの終了時に)、選ばれた命令のデスティネーション物理レジスタについて、「この物理レジスタは読み出せる状態である」と記録します(表に1を書き込みます)。

スケジューラが管理してる「物理レジスタは読み出せる状態である」は実際には2サイクル(1サイクル+発行レイテンシ分)未来の状態を表しています。 普通はこれはレイテンシ予測に基づくもので、未来の状態が本当にそうなるかはわかりません。 未来予測が正しくなくなる典型的なパターンはキャッシュミスした場合です。 キャッシュミスすればデータは得られないので、「物理レジスタは読み出せる状態である」というのは間違いになります(レイテンシ予測ミス)。 間違った未来予測に基づいてスケジュールすると、正しいデータを使って実行できないので、プロセッサの動作がおかしくなります。 したがって、未来予測が間違った場合の対処を考える必要があります。 これは、アウトオブオーダープロセッサを作るうえで、難しい問題です(リザベーションステーションを使った方式はこの問題が起こらないという利点が一応あります。また、レイテンシが不確かなものについては実際にデータが得られるまで表を更新しない(投機的発行しない)という方式もあります)。

しかし、今回はデータキャッシュがなくてレイテンシ予測は完全に当たるので、簡単に作ることができます。

「物理レジスタは読み出せる状態である」の表は48個ある物理レジスタに対して1bitの情報を管理しているだけですが、ポート数がすごく多いので贅沢なつくりになっています(とはいえFPGAを前提にするとこれより小さく作る方法はそう簡単にはなさそうです。この作り方は結構いい方法なのでは?)。 この表は、スケジューラ内の16個の命令が2つのソースについて確認するので、32readになっています。 しかも、ディスパッチの1writeに加えて、レイテンシが1cycle/2cycle/3cycleの三パターンがあるので3writeが必要で、合計で4writeになります。 つまり、この表だけで〈1bit×48エントリ×32read×4write〉が必要になっています。

ここで選ばれた命令(実際にはその命令の情報が格納されているpayload RAMのインデクス番号)は、発行ステージへ送られます。

発行ステージ

Payload RAMを読み出し、PCと命令語などを得ます。 命令語はここでデコードします。 また、物理レジスタからの読み出しを行います(分散RAMで作って非同期読み出しします)。 さらに、必要であればフォワーディングを受け取ります。 そして、それらのうち適切な方をパイプラインレジスタに書き込みます。

実行ステージ

整数命令

命令を実行するだけです。 命令の実行結果は、即フォワーディングネットワークに流します。

一方、レジスタファイルへの書き込みはまだ行いません。 これは、レジスタファイルの書き込みポートを節約するためです(後から考えると、これは「資源を贅沢に使って簡単に作る」に反していました……)。 整数命令の実行レイテンシは1サイクル、最長の実行レイテンシを持つのはロード命令の3サイクル、ということで、差分の2サイクルは待ちます。 これは普通の五段パイプラインと同じ設計です。 レジスタファイルに書き込まれるまでの間は、フォワーディング経由で値を受け取れるので、これが問題になることはありません。

また、レジスタファイルに書き込むのと同時に、リオーダーバッファに実行完了フラグを書き込みます。 リオーダーバッファといっても、実行結果データを書き込むタイプではありません。

もし分岐命令であれば、実行の結果得られたnext pcが予測結果と一致しているかを確認します。 予測結果と一致していてもしていなくても、分岐予測器を学習します。 また、予測結果と一致していなければ、分岐予測ミスです。 例外の発生するサイクルをロードストア命令と合わせるため、第二サイクルで分岐予測ミス信号を送出します。 これは、複数の実行ユニットから出てきた例外の中で最古のものを選ぶ、という面倒事を回避するためです。 さすがに次のサイクルから正しい命令をフェッチするのはレイテンシ的に厳しいため、このサイクルではmispredフリップフロップに1を書き込むのにとどめ、次のサイクルから分岐予測ミスの対処を開始します。

乗算命令

2サイクルかけて乗算を行うだけです。 乗算結果は即フォワーディングネットワークに流しますが、レジスタファイル・リオーダーバッファへの書き込みは、1サイクル待ちます。 ロード命令に合わせて3サイクルにしておけばこの辺やスケジューラを多少簡単にできたので、レイテンシを2サイクルにしたのは失敗でした(乗算の結果をまた乗算に使うということはあまり発生しないので、レイテンシが長くてもIPCへの影響は軽微です)。

ロード命令

第一サイクルではメモリアドレスを計算します。TLBがないので単純です。

第二サイクルではメモリアドレスを使って、データメモリから読み出します。キャッシュがないので単純です。 さらに、メモリアドレスを使って、ロードストアキューを連想検索します。 もし、ロードストアキューに当該メモリアドレスに対するストアが記録されていれば、その中で最新のストア命令の書いた値を読み出します(ストアフォワーディング)。 また、ロードストアキューに、アドレスを記録します。

第三サイクルでは、ロードストアキューを連想検索して値が得られたならそれを、得られなかったならデータメモリから読み出した値を、零拡張または符号拡張します。 その結果は、フォワーディングネットワークに流します。 また、第三サイクルが終わる瞬間にレジスタファイル・リオーダーバッファに書き込みます。 これにより、発行ステージでのレジスタ値の選択が「レジスタファイルから読み出した値」「今完了した整数命令の結果」「今完了した乗算命令の結果/1サイクル前に完了した整数命令の結果」「今完了したロード命令の結果/1サイクル前に完了した乗算命令の結果/2サイクル前に完了した整数命令の結果」の4択になります。

※ストアフォワーディングはパーシャルメモリアクセス(SB命令で1byteストアしたのにLD命令で4byteロードするなど)を考えると面倒です。メモリアドレスの下位2bitごとにロードストアキューを作ることで対処しています。

ストア命令

第一サイクルではメモリアドレスを計算します。TLBがないので単純です。

第二サイクルでは、メモリアドレスとストアデータをロードストアキューに書き込みます。 さらに、メモリアドレスを使って、ロードストアキューを連想検索します。 もし、ロードストアキューに当該メモリアドレスに対するロードが記録されていれば、メモリ順序違反が起こったことを示しています。 その中で最古のロード命令を探し出し、その命令からフェッチしなおします(リフェッチ)。 このサイクルではmispredフリップフロップに1を書き込むにとどめ、次のサイクルからメモリ順序違反の対処を始めます。

第三サイクルは何もする必要がありません。

コミットステージ

アーキテクチャステートを不可逆的に更新するのがコミットステージです。 リオーダーバッファを確認し、最古の命令に対応するエントリに実行完了フラグが立っていればコミットを行います。 ストア命令であれば、ストアキューに書き込まれている値を読み出し、データメモリに書き込みます。 キャッシュがないので、これは確実に成功します。 ストア以外の命令であれば、何もする必要はありません。

その後、リオーダーバッファのheadポインタ*1 を1インクリメントします。 このheadポインタは、以下の情報を得るために使われます。

  • リオーダーバッファ中の最古の命令はどれか
  • スケジューラ内に入っている命令の中で「新しい命令」と「古い命令」の境界はどこか

また、分岐予測するごとに1インクリメントされるtailポインタと組み合わせて、以下の判定に使われます。

  • リオーダーバッファに空きがあるか
    • 分岐予測を進めていいかの判定で必要
  • ロードストアキューの中で有効な範囲はどこか
    • メモリ順序違反が起こったかの判定で必要

本当はここで、記録されている物理レジスタ番号をフリーリストに返却する必要があるはずです。 しかし、実は逆の発想で「リオーダーバッファの(headではなく)tailを進めるときに、物理レジスタ番号をフリーリストに返却」としても大丈夫です。 こうすると物理レジスタ番号がフリーリストに返却されるのが遅くなってうれしくないですが(というか物理レジスタ不足でストールが発生してtailが進められなくなればデッドロックしてしまう)、今回は十分な物理レジスタがあるため問題ありません。 パイプラインフラッシュが起こった時のことを考えるとこうしたほうがやりやすいため、この方式を採用しました。

パイプラインフラッシュの実装

さて、ここまででパイプラインの“正常”動作を書いてきましたが、実際には分岐予測ミスやメモリ順序違反などの投機ミスが起きるので、パイプラインフラッシュが発生します。 この時にプロセッサの状態を正しく戻さないと正しく動かないので、頑張って作ります。 といっても、今まで贅沢に作ってきたのはパイプラインフラッシュが起こった時にほとんど何もしなくていいようにするためなので、実装はそんなに難しくありません。

分岐予測ステージ

  • グローバル履歴のリングバッファのポインタを戻す
    • 一応、これをちゃんとやらなくてもプロセッサは正しく動くけれど、分岐予測がめちゃくちゃになる
  • 次のサイクルの分岐予測の“種”となるPCをセットしなおす

フェッチステージ

  • 正しいPCを受け取る(まだフェッチはしない)

デコード・リネーム・ディスパッチステージ

  • RMTの状態をチェックポイントから復元する(ポインタを正しい値に修正するだけ)
  • 今デコード中の命令は捨てる

スケジュールステージ

  • フラッシュの原因となった命令以降の命令のエントリについて、有効な命令が入っているかを修正する(validビットを0にする)
  • 物理レジスタのreadyの状態は、戻さなくてよい(ディスパッチ時に正しい値に修正されるので)

発行ステージ

  • 何もしなくてよい

実行ステージ

  • リオーダーバッファのheadとtailの範囲に無い命令であれば
    • mispredフリップフロップに書き込まないようにする
      • これは必須
    • レジスタファイルへの書き込みを行わないようにする
    • フォワーディングネットワークに値を流すのもやめる
      • この二つは今回のパイプライン設計だとやらなくても大丈夫な気がする

コミットステージ

  • 何もしなくてよい(もともとリオーダーバッファのtailより先をコミットしないので、tailを戻しておけば大丈夫)

物理レジスタガベージコレクション

分岐予測が行われて命令がリオーダーバッファに格納される(本当はまだ命令がフェッチされていないのでエントリが確保されるだけ)たびに、以下のどちらかを行えばよいです。

  • 「フラッシュされた」ビットが立っていれば、その命令が上書きするリオーダーバッファエントリに記録されている命令のデスティネーション物理レジスタ番号をフリーリストに返却(投機ミスからの回復動作)
  • 「フラッシュされた」ビットが立っていなければ、その命令が上書きするリオーダーバッファエントリに記録されている命令のデスティネーション論理レジスタにもともと対応していた物理レジスタ番号をフリーリストに返却(通常動作)

こうすることで、フリーリストの書き込みポート数を1のままにすることができます。 というかよく考えると、リオーダーバッファに記録されている情報(フラッシュビット+旧物理レジスタ番号+新物理レジスタ番号)がフリーリストの情報を特定するのに十分なものになっているので、フリーリストをわざわざ実装する必要はなかったようです……。

アウトオブオーダープロセッサの実行の様子

スーパースカラでないので、レイテンシが1サイクルの命令が並んでいるところはインオーダープロセッサと同じ感じになります。 CoreMarkというベンチマークを動かしてみたところ、以下の二か所ではアウトオブオーダー実行がうまく働いていそうでした。 パイプラインは、Konataというパイプライン可視化ツールで可視化しました。

図1: 乗算を行っている部分のパイプラインチャート。暇なサイクルでmove命令が実行できている(分岐予測ミスが起こるかもしれないのに先に実行できる!)。でもこの命令よく見ると無意味では……?
図2: ロードの連鎖がある部分のパイプラインチャート。5命令に6サイクルかかる(これは最古の命令を優先するという戦略がいまいちなため)のでリオーダーバッファ不足で性能が下がっているわけではない

まとめ

キャッシュなどの複雑化の要因をなくし、それでも複雑な部分は大量の資源を使って解決することで、簡単にアウトオブオーダープロセッサを作ることができることがわかりました。 来年はちゃんと省資源化(とできればスーパースカラ化)してみたいと思います。

*1:いろいろな流儀があるようですが、一番古い命令が「待ち行列の先頭」なのでhead、一番新しい命令が「待ち行列の末尾」なのでtail、という流儀を採用しています。パイプラインチャート的にも上がheadになっていい感じです。でも、headの方がなんか最新感があるし、動物は普通は頭がある方向に進む(蛇とか魚とか)のでheadがキューの先頭(のびる方向)であり一番新しい命令を指すべき、という意見も否定しがたいです。

0.1+0.2≠0.3について

浮動小数点数を扱う時に注意しなければいけない例として、0.1 + 0.20.3にならない、というものがあります。 具体的には、0.1 + 0.20.30000000000000004になってしまい、0.3にならないプログラミング言語が多数存在します。 この問題について、以下の記事では十進数から二進数に変換するときの誤差が原因であるとしています(この説は他にも多くのウェブサイトで唱えられています)。

qiita.com

たしかに十進数から二進数に変換するときの誤差は一因なのですが、より大きな寄与を持っているのは足し算の丸め誤差です。 そのことを確認していきたいと思います。

(わかっている読者への注:以下では、浮動小数点数と言ったらIEEE754で規定されたbinary64(二進倍精度浮動小数点数)を指すこととします)

浮動小数点数の基礎

基本的には、分子が 2^{52}以上 2^{53}未満の整数で、分母が二の整数乗であれば、浮動小数点数として表せます。 この形式で表せない数は、浮動小数点数として表せる数で近似します。 この近似を行うことを、丸めと言います。 普通は、最も近い浮動小数点数で近似します(最近接丸め。端数が0.5ぴったりの時のタイブレークの方式の違いに四捨五入や偶数丸めがあります)。 特別な用途(例えばお金が虚空から生まれてはいけないとか*1)では、「真の値を超えない範囲で」最も近い浮動小数点数に丸める(下向き丸め)などを使うこともあります。逆に「真の値を下回らない範囲で」とする方法もあります(上向き丸め)。

ULPはunit of last placeの略で、上記のような分数として表したときの、分子の1の違いを単位としたものです。 誤差は0.38ULP、のような形で誤差の解析によく使います。

計算の過程

0.1を浮動小数点数に変換する

0.1に最も近い浮動小数点数は、 \frac{7205759403792794}{72057594037927936} = \frac{\lceil0.1\times2^{56}\rceil}{2^{56}}です(十六進浮動小数点数リテラル表記で0x1.999999999999ap-4)。 この時点で、分子が7205759403792793.6であるべきところが7205759403792794になってしまっているので、+0.4ULPの誤差を生じています。

0.2を浮動小数点数に変換する

0.2に最も近い浮動小数点数は、 \frac{7205759403792794}{36028797018963968} = \frac{\lceil0.2\times2^{55}\rceil}{2^{55}}です(十六進浮動小数点数リテラル表記で0x1.999999999999ap-3)。 この時点で、分子が7205759403792793.6であるべきところが7205759403792794になってしまっているので、+0.4ULPの誤差を生じています。

0.3を浮動小数点数に変換する

0.3に最も近い浮動小数点数は、 \frac{5404319552844595}{18014398509481984} = \frac{\lfloor0.3\times2^{54}\rfloor}{2^{54}}です(十六進浮動小数点数リテラル表記で0x1.3333333333333p-2)。 この時点で、分子が5404319552844595.2であるべきところが5404319552844595になってしまっているので、-0.2ULPの誤差を生じています。

0.1と0.2を足す

丸める前の加算の結果は、 \frac{5404319552844595.5}{36028797018963968}です。 ここで、分子は 2^{52}以上 2^{53}未満ですが、整数になっていないので丸めを行う必要があります。 四捨五入なり偶数丸めなりを適用すると、 \frac{5404319552844596}{36028797018963968} = \frac{\lceil0.3\times2^{54}\rceil}{2^{54}}を得ます。

これは、0.3に最も近い浮動小数点数である \frac{5404319552844595}{18014398509481984} = \frac{\lfloor0.3\times2^{54}\rfloor}{2^{54}}と分子が1だけ異なります。

誤差の解析

0.1の指数部は0.3の指数部に対して2だけ小さいので、変換の誤差+0.4ULPが最終結果に与える影響は1/4となって+0.1ULPです。 0.2の指数部は0.3の指数部に対して1だけ小さいので、変換の誤差+0.4ULPが最終結果に与える影響は1/2となって+0.2ULPです。 加算の丸め誤差が最終結果に与える影響は+0.5ULPです。

よって、合計で+0.8ULPの誤差が生じます。

このように、この誤差の主要因は加算の丸め誤差です。

いつでもfaithfulになるか

Faithful rounding(信頼可能丸め*2)とは、真の値を上向き丸めして得られる値か真の値を下向き丸めして得られる値の、どちらかが返ってくるような丸めのことです。 言い換えると、真の値を挟む二つの浮動小数点数のうちのどちらかが返ってくるということです。 「最近接の浮動小数点数に丸めてね」というのは要求が厳しすぎるので、少し要件を緩和したということです。

先ほどの例では、真の値0.3に対し、最も近い浮動小数点数 \frac{\lfloor0.3\times2^{54}\rfloor}{2^{54}}を計算結果として得ることができませんでした。 しかし、得られた結果 \frac{\lceil0.3\times2^{54}\rceil}{2^{54}}は真の値を上向き丸めして得られる値です。 したがって、先ほどの計算は、最近接丸めではありませんでしたが信頼可能丸めではあったようです。

さて、十進小数同士の足し算結果が必ずしも最近接丸めではないということがわかりましたが、要件を信頼可能丸めに緩めれば満たされるでしょうか?

答えは否です。 まず、桁落ちが生じる場合は明らかにかなりの誤差が生じ得ます。 桁落ちが生じない場合に限っても、信頼可能丸めですらなくなる例は存在します。

0.17 + 0.28がその一例です。

0.17を浮動小数点数に変換すると、 \frac{6124895493223875}{36028797018963968} = \frac{\lceil0.17\times2^{55}\rceil}{2^{55}}となり、+0.44ULPの誤差が生じます。 0.28を浮動小数点数に変換すると、 \frac{5044031582654956}{18014398509481984} = \frac{\lceil0.28\times2^{54}\rceil}{2^{54}}となり、+0.48ULPの誤差が生じます。 真の値である0.45を浮動小数点数に変換すると、 \frac{8106479329266893}{18014398509481984} = \frac{\lceil0.45\times2^{54}\rceil}{2^{54}}となり、+0.2ULPの誤差が生じます。

二つの浮動小数点数が表す値を加算すると、 \frac{8106479329266893.5}{18014398509481984}となり、四捨五入なり偶数丸めなりをすると \frac{8106479329266894}{18014398509481984}となります。 これは真の値 \frac{8106479329266892.8}{18014398509481984} =  \frac{0.45\times2^{54}}{2^{54}}を挟む二つの浮動小数点数ですらありません。 つまり、加算を実行するだけで、信頼可能丸めにすらならないケースがあることがわかりました。

このケースの誤差は+1.2ULPです。 内訳としては、0.17を変換した時の誤差+0.44ULPが最終結果に与える影響は1/2されて+0.22ULP、0.28を変換した時の誤差+0.48ULPが最終結果に与える影響は+0.48ULP、加算の丸め誤差が0.5ULP、で合計+1.2ULPです。


このような場合、最大の誤差は1.25ULPとなります。 まず、変換の誤差は最大で0.5ULPです。 加算のオペランドの変換誤差が最終結果に与える影響は、最大で0.5ULPですが、両オペランドともそうなることはありません。 なぜなら、そのような場合は繰上りが発生して、変換誤差が最終結果に与える影響が1/2になってしまうからです。 よって最悪ケースでは、片方のオペランドの変換誤差が最終誤差に与える影響が0.5ULP、もう片方のオペランドの変換誤差が最終結果に与える影響が0.25ULPとなります。 これに加えて加算の丸め誤差が最大で0.5ULPなので、合計で最終結果に乗る誤差は最大で1.25ULPであるということになります。

最大誤差が1.0ULP未満であれば信頼可能丸めになりますが、1.0ULPを超えているため、信頼可能丸めにはなりません。 ただし、最大誤差が1.5ULP未満なので、「真の値を最近接丸めした値か、またはその前後の浮動小数点数」が計算結果となることは保証されます*3

*1:言うまでもなく、お金を扱う時に浮動小数点数を使うべきではないですが、例として

*2:忠実丸め、隣接浮動小数点数丸め、F丸め、などの訳も存在

*3:指数部の境界付近ではULPを使った議論は怪しいです。反例があったら教えてください。