# この記事は、2015年12月11日にYahoo!ブログへ投稿した記事を、Yahoo!ブログの終了に伴って2019年8月に移行したものです。

CPU実験 Advent Calendar 2015の記事です。

読者層として3年の夏休み頃のIS 16erを想定して書かれています。

この記事では、CPU実験で動かすことになるraytracer/min-rt.mlの中で必要な浮動小数点数周りの演算についての軽い解説を行います。解説しなくてもその内分かるかとは思いますが、一応記録のためです。

以下の記事では2015年度のCPU実験のルールを元にした記述がたくさんあります。ルールが変わっているかもしれないので注意してください。

浮動小数点数演算について見る前に、IEEE 754 単精度浮動小数点数について復習します。

IEEE 754 単精度浮動小数点数は、32bitで実数を近似的に表現するための規格で、符号部1bit, 指数部8bit, 仮数部23bitからなります。

指数部は127を0とするオフセット付きの表現になっており、また、仮数部はけち表現(1.xxx…のxxx…しか保存しない表現)をすることによって1bit節約されています。

単精度で表現できる通常の実数(正規化数)の他に、±inf, NaN, 非正規化数があります。また0には2種類(+0と-0)あります。

ここらへんについて怪しい方は是非今の内に“Wikipedia - IEEE 754”や生の規格(東大内だとpdfで落とせます)を読んでおいてください。以下では上の知識は既知とします。

「浮動小数点数演算」と言いましたがこれは曖昧で、実際に実装するときには「ハードウェア実装する浮動小数点数演算」「コンパイラのライブラリ関数としてソフトウェア実装する浮動小数点数演算」「min-caml言語として実装する浮動小数点数演算」の3種類があります。これらの演算をFPU係やコンパイラ係が実装することになります。

どこまでの演算をどの分類にするかは、班のISA(命令セットアーキテクチャ)やコンパイラ・ライブラリの仕様によります。どのような設計にするとより高速にレイトレーシングが実行できるか、などを目標にこれらを定めることになります。

この記事の目標は、必要な演算について適当な解説をして、1st アーキテクチャの設計前になるべくクリアなイメージを持ってもらうことにあります。

では、どのような演算が必要になるのか列挙していきます。

加算命令faddと乗算命令fmulは、全ての班が1st アーキテクチャでハードウェア実装することになるでしょう。おそらくハードウェア実験でfaddの実装はしたことになっているでしょう。まだしていない人はFPU係として実装することになります。faddの実装にはハードウェア構成法の方の小林先生の著書である『ディジタル・ハードウェア設計の基礎と実践』が多いに参考になるでしょう。地下に数冊あるので是非読んでみて下さい。faddのC実装がちゃんとできているかどうかは綿密にチェックしましょう。14erのb-inaryさんが作られたverify.cが参考になるかもしれません(14er wiki/ハードウェア実験にあります)。ただし安全のためにはもっと標本数を増やすべきです。faddが実装できればfmulの実装は容易にできるでしょう。2004年に書かれた神FPU資料のfmul実装は古いです。16erの皆さんが扱う基盤がリニューアルされたものなら、他の部分の記述も古いものになってしまっているかもしれないので注意が必要です。

faddのC実装はある程度真面目に作っておく必要があります。というのも、CPU実験の狭義完動条件が「シミュレータの出力と実機の出力のdiffが0」なので、シミュレータに組み込まれたfaddの実装が雑だとVHDL実装とのdiffが出てしまいます。

NaNや非正規化数などにどこまで対応するかどうかは、CPU実験のレギュレーションやraytracer/min-rt.mlの中で行われる演算と相談しながら行うことになります。レイトレーシングに満足せずもっと豊富なことをしたいのであれば実装しても良いでしょう。はじめは面倒くさいし実装していなくていいんじゃないかあと思います。実際にレイトレでNaNや非正規化数が現れるのかどうかは、コンパイラ、アセンブラ、シミュレータ完成後にシミュレータで確認することはできます。本当にCPU実験のレギュレーションかちかちで作るのであれば、丸めの計算をroundingをnearest evenなどといったものでなくもっとサボることもできるらしいのですが、僕は実際に実装していないので分かりません。速さを求めるのであれば神FPU資料なども参考に研究することになります。ただしfaddの誤差をサボると三角関数などのソフトウェア実装した関数で誤差が大きくなって困ることになるかもしれないので加減が必要です。

VHDL実装をパイプライン方式で行ったとして、66MHzでfaddが3 clock、fmulが2 clockで終わればとりあえずは大丈夫です。ただしその内クロックを上げることになるので、それに対応することになるでしょう。また、うまく丸め計算をサボるともう少し高速化することもできるようです。13er grafiさんの記事には2パス加算器の話が載っています(ヘネパタにアイディアが載っているらしい。こちらも参考になる)。

加算命令と乗算命令があるなら減算命令fsubと除算命令fdivも欲しくなる所ですが、もしISAをなるべく小さくしたいのであればこれらはとりあえず必要ありません。代わりに正負反転命令fnegと逆数命令finvが必要です。

fnegの実装は簡単です。finvの実装はNewton-Raphson法を用いるか、線形近似の値を保存しておくテーブルを用意しておいて、テーブル引きを使って値を求める方法の2種類に分かれます(他の方法もあるのかもしれませんが知りません)。

Newton-Raphson法で実装する場合“Wikipedia:en - Division algorithm”が参考になります。テーブル引きで実装する場合は神FPU資料が参考になるでしょう。finvには大体2 clockかかります。ただし分散RAMを使って1 clockで終わらせる実装も発見されていて、こちらはハードウェア資源との相談になります。

テーブル引きの誤差を減らす方法はよく知りません。FPU係の誰かが記事にしてくれるかも?

コンパイラとしてmin-camlをそのまま使うなら、浮動小数点関係の中間表現としてあらかじめ存在するものはFMov, FNeg, FAdd, FSub, FMul, FDiv, LdF, StF, IfFEq, IfFLEの10個です。楽をするなら1stではこの10個さえできればOKです。FSubはfaddとfnegで表現できます。fdivは実はx/yをx * (1/y)と最適化するのは正しくない(たとえばx=y=0)のですが、CPU実験のレギュレーション的にはOKなので大丈夫です。

僕の班の場合、1stが完動した後にシミュレータを使ってraytracer/min-rt.mlでの色んな命令の出現頻度を測った所fsubが結構多かったので2ndではfsubを命令セットに加えることにしました。fdivは今の所入れていません。ここらへんは各班で色々画策して決めることになるかと思います。

LdF, StFはGPR/FPR(整数レジスタと浮動小数点数レジスタ)を分けている班は用意する必要のある命令です。分けていない班は中間表現ごと消して大丈夫(なはず)。面白命令入れたい班の人はx87の定数ロード命令などいかがでしょうか(速くなるかどうかは知りません)。役に立たなさそうな面白命令として、上位/下位n bitだけが即値オペランドで指定できる浮動小数点数ロード命令とか(無いところは0埋め)。

fmovはGPR/FPR分けてる班なら有ると思うのでOKかと思います。ちなみに凄い細かい話をすると”fmov x to y”を”y = fadd 0.0 x”で実装するのは間違いです。xが-0.0のときを考えて下さい。

浮動小数点数によるブランチ命令は最低限equalとless than (もしくはless than or equal to)があれば大丈夫です。いや待てよ、本当に最小の命令セットにするならequalも不要ですね(ニッコリ)。コンパイラに優しいのはeq, neq, lt, leq全部あるパターンかな?

平方根命令fsqrtの計算はfinvと同じくNewton-Raphson法かテーブル引きを使うでしょう。ここらへんからハードウェア実装するかソフトウェア実装するかが分かれてくると思います。

面白命令として、高速に平方根の逆数を求めるfast InvSqrtというのが知られています。特にfsqrtをソフトウェア実装している場合「平方根の逆数を求めてから乗除算に使う」という方法で高速化の可能性があります。計算には有名なマジックナンバーがあるので参考になるかも。ただし誤差がちゃんとレギュレーションを満たすように注意してください。1 iterationでは満たさないと思います。

三角関数のsin/cosは入力xの値によって誤差が大きくなるので、適当に誤差の小さい区間に帰着させてあとはマクローリン展開(をちょっと弄ったもの)によって計算するという方法が神FPU資料に載っています。また、15erには自力で色々試してみている人もいたようです。

sin/cosは1stではソフトウェア実装の人が多かった印象です。レイトレーシングではそんなに呼ばれない[要出典]ということもありハードウェア実装のメリットがそんなに無い班も多いです。自分たちの班のISAと相談することになります。

sin/cosが出来上がったら誤差をulpでちゃんと表示してあげることをオススメします(参考)。もしかすると、たとえばsinでx=nπの周辺において誤差が大きくなっているかもしれませんが、実はそれはCPU実験のレギュレーション内です。しっかりレギュレーションを確認すると、定数cというとても1に近い数が設定されていることに気づきます。このcの意味に気づけば、cを含んだ意味でcos/sinがレギュレーションを満たしていることを全数探索できます。

あ、そうそう、連続系アルゴリズムの講義で習うかとは思いますが、級数の計算はHorner法でやると良いです。

tanもレイトレの中で使われていますが、raytracer/min-rt.mlの中身でtan = sin/cosと定義されてしまっています。ただしCPU実験のレギュレーションにより、勝手に別のものに置き換えても良いことになっています。

atanもsin/cos同様にマクローリン展開で計算できます。僕はほとんど神FPU資料に従ってしまったのであまり詳しいことは分かりません……。

浮動小数点数と整数の変換命令としてint_of_float(ftoi, f2i)とfloat_of_int(itof, i2f)があります。神FPU資料には8388608を使った実装法が載っていますが、他の実装として、GPR/FPRが分かれていないなら指数部を読んで適当にシフトしてあげるというのもあります。ちなみにCPU実験のレギュレーションに沿ったint_of_floatとOCamlのPervasives.int_of_floatは微妙に違うので注意してください。

int_of_floatには、爆速で終わる実装が発見されています(ただしISA的な注意が必要)。詳しくはこちらを: (1) / (2) / (3) / (4)

似た命令にffloorがあります。ffloorは神FPU資料を見て何も考えずに作ると間違えるので注意してください。「床関数をミスるとレイトレの床がバグる」という超有用な金言があるので、床がバグったらいの一番にffloorを疑って下さい。

int_of_float, float_of_int, ffloorはソフトウェア実装するかハードウェア実装するかの選択で悩むことになるかもしれません。呼び出し回数や1回の呼び出しでの消費クロック数などを考えながら選択することになるでしょう。

raytracer/min-rt.ml特有の命令として、fless, fispos, fisneg, fiszero, fhalf, fsqr, fnegがあります。これらはminiMLRuntime.mlに定義されています。

これらはCPU実験のレギュレーションによりある程度弄くることができます。min-caml言語として実装してインライン展開しても良いですし、たとえばfhalfはGPR/FPRの区別が無いアーキテクチャの場合はfmulしなくても単に指数部を弄るだけで実装できます。ここらへんは小細工です。

FPRを使うライブラリ命令は実はもう1つあって、それがread_floatです。これはRS232Cから送られてきたASCIIの小数を32bit floatに変換してレジスタに格納する命令です。別にread_intも作ると思うのでそれをそのまま拡張すると作ることができます。小数点が無かったりすることもあるので注意してください。ASCIIの1文字は8bitなので、8bit単位で読み書きする命令があると便利です。レギュレーションが変わっていなければ審判サーバーとの通信過程で1byteだけのデータを送信する箇所があるので、I/Oがいつも32bit単位だと思っていると失敗すると思います。

また、審判サーバーが変わっていないのなら、更に32bitのバイナリを読んでその内容をそのままFPRに格納する命令が必要なはずです(そしてこれはすぐ実装できます)。

以下は瑣末な点です。

誤差を測る: 浮動小数点数まわりの演算を実装するときは、誤差がちゃんとレギュレーションを満たすことを確認した方が良いでしょう。そこまで細かくなくても、グラフの概形が合っていることを見るだけでも何かバグってないか見ることができます。1引数の関数は全数検査が現実的な時間でできるのでやってみるといいかもしれません。15erはシミュレータ経由でデータを出させてグラフプロットする人が多かったようですが、14erには実機で直接検査していた方もいたようです。

また、x87との誤差を見るときは注意が必要です。C言語でやるなら、全ての表現がfloatになっていないと、途中でdoubleの計算になってからfloatにキャストされるのは結果が違うかもしれません。math.hの関数を使うときも、たとえばsinの代わりにsinfを使うべきでしょう。また、もう無いと思いますがx64ではなくてx86でやるときは、うまくコンパイルしないと浮動小数点数演算の際に一度80bitに拡張されてから計算されてしまうので誤差がでます(参考1、参考2)。また、OCamlは倍精度です。ついでに言うとOCamlの整数型は1bit無いです(GCに使われている。ビット幅をフルに使うにはnativeintを用いる)。

シミュレータで数える: どの命令が要りそうか判断するために、どの命令が何回実行されたかカウントするのは便利でした。関数呼び出しの回数を測るのはシミュレータの組み方によっては難しいかもしれませんが、うちはたまたま標準だとnopが1回も実行されないので、回数を数えたいところにnopを手で挿入してやりました。

ISAについて考える: どこまでハードウェア実装してどこからソフトウェア実装するかは結構悩みどころかもしれません。でもたとえば“Computer Organization and Design”には「ぶっちゃけISAはそこまで関係ない」と書かれていたりもします。ISAに神経質になりすぎるよりコアの設計やコンパイラの最適化を考えた方がよっぽど速度に寄与するのかもしれません。また、数年前までは神FPU資料に結構頼っていたと思われるFPU実装も、そこまで単純なものには見えなくなってきました。

いくつかの実装については、14er wiki/CPU実験のSUPERHACKERのCPU実験/Yebiや神FPU資料が参考になります。浮動小数点数演算周りの話は、班の全係が関係して色々考えられる面白い部分だと個人的に思っているので、是非色々試して面白実装してみてください(面白いのができたら、是非教えてください!)。