公開日：2021/06/05

16th Season

#16-5　半精度浮動小数点演算ユニット

　FPGAボード(Sipeed Lichee Tang Primer)用に作成したオリジナル16ビットCPUコアMoscoviumに関しては、乗除算コプロセッサを追加したことで、CPU命令で整数の四則演算を実施することが可能となっています。
　次に強化するのは、浮動小数点演算です。

　PC等では、通常、32ビット単精度浮動小数点演算(float型)や、64ビット倍精度浮動小数点演算(double型)が使用されています。しかし、16ビットCPUコアであるMoscoviumにとっては、32ビットのfloat型でさえ、手に余るのです。
　「さすがにMoscoviumに浮動小数点演算は厳しいな……」、そう思っていたところ、面白いものを見つけました。16ビット半精度浮動小数点演算(half型)と言う物がある様です。

IEEE 754 binary16

　PC等でfloat型やdouble型の表現に用いられていて、一般的となっているのは、IEEE 754と言う規格です。この中に、binary16と言うフォーマットが規定されていて、これが今回対応したいと考えている浮動小数点表現です。

　　IEEE 754 binary16を検索する(別ウインドウが開きます)

　half、half-float、FP16等の名前で呼ばれているようですので、本HPではhalf型と表記したいと思います。

　実は、IEEE 754の中では、実際に演算に使用するfloat型やdouble型のフォーマットと異なり、half型のフォーマットは「交換形式」として定義されていて、データ転送やストレージに使用される形式とされています。

half型バイナリ表現

　half型のバイナリデータは以下図のような形式となっています。



　符号(Sx)1ビット、指数(Ex’)5ビット、仮数(Fx’)10ビットの計16ビットです(↓)。

　　half型X＝(−1)^Sx×0x1.<Fx’>×2^Ex’−15

• 符号Sx：
• X≧0の場合、Sx＝0。
• X≦0の場合、Sx＝1。
• 不等号の両方に＝が入っているが、正しい。half型には、(＋)0と−0の両方が存在する。
• 指数Ex’：
• 指数Exは、−14〜15の値を取ることが可能で、＋15のオフセット(下駄)を履いたEx’＝Ex＋15が、5ビット符号無整数フィールドとして記録されている。
• Ex’＝1〜30であり、5ビットの符号無整数では、0と31が抜けている。この二値は以下の値として意味が定義されている。
• Ex’＝0でFx’＝0の場合は、ゼロ(0)。
• Ex’＝0でFx’≠0の場合は、非正規化数。
• 絶対値が2⁻¹⁴よりも小さな数値の表現に割り当てられている。
• Ex’＝31でFx’＝0の場合は、INF(無限大、∞)。
• Ex’＝31でFx’≠0の場合は、NaN(非数、Not a Number)。
• 仮数Fx’：
• 仮数Fx’は、上で図示したFxの数式の通り、最上位ビットの上に1と小数点が省略されている。
• 従って、実際の数値Fx＝0x1.<Fx’>である(<Fx’>は左詰め)。
• Fx’が10ビット幅なので、省略された小数点上の1と合わせて、Fxは11ビット符号無固定小数点数となる。
• ただし、省略された1は常に1なので、ダイナミックレンジとしては紛れもなく、Fx’の10ビット分しかない。
• 10進数換算の有効数字は、log₁₀(2¹⁰)≒3.010と、3桁程度。

　10進数換算の有効数字としては、たったの3桁程度と言う事で、円周率だとπ＝3.14が関の山です。ことによると近似値22/7≒3.14285714286と変わりません（１＾＾；）。22/7なら、整数演算できますので、「こんな精度で浮動小数点演算要るのか？」ってのが第一印象です。

※単精度浮動小数点数(float型)の解説は、こちらのページにあります。half型とfloat型は、指数と仮数のビット数が異なるだけで基本的には同じ構造をしています。

用途

　half型で表現可能な数値は、以下の様になります。

• 絶対値最大値：
• Ex＝15、Fx＝0x1.ffcの場合が絶対値最大で、この時のXmax＝65,504。
• 符号ビットが独立しているので、表現可能な数値の範囲は、−65,504≦X≦65,504となる。
• 16ビット符号無し整数の最大値が65,535なので、「ほぼ同程度」というべきか、「uint16にも満たない」と言うべきか……。
• 演算結果がXmaxを上回る(オーバフロウする)と、INFとなる。
• 絶対値最小値：
• Ex＝−14、Fx＝0x1.000の場合が絶対値最小で、この時のXmin＝6.10352×10⁻⁵。
• 演算結果がXminを下回る(アンダフロウする)と、0となる。
• ※実は、Xminよりさらに小さい数を扱うための「非正規化数」という表現があるのだが、今回は非対応なので、説明省略。

　絶対値で、約6.10×10⁻⁵〜65,504と言う、非常に広いダイナミックレンジを扱えることから、0〜65,535までの整数しか扱う事の出来ないuint16型で不足する用途での交換形式として使用されているようです。何しろhalf型は、uint16型と同じ16ビット幅ですからね。
　具体的には、GPU(Graphics Processing Unit)などで使われているのだとか。へぇ〜。

　次節でもう少し具体的な設計を進めますが、仮数Fxが11ビットしかないので、演算回路はfloat型やdouble型と比べてかなり小さくすることが出来ます。このためhalf型の演算回路を多並列に構成して、並列処理のスループット向上などにも応用が進んでいるようです。

※最近のPC用CPUが採用している64ビットのアーキテクチャ相当であれば、その中に16ビットhalf型のデータを64/16＝4個も格納することが出来て、4数の並列演算が可能となります。128ビットなら、更にその倍。

FPUを作ろう

　然して、今回はMoscovium周辺演算回路として、FPU(Floating Point Unit)を設計しようと思っています。

　16ビットCPUであるMoscoviumには、16ビット汎用レジスタが8個しかありません。32ビットfloat型を汎用レジスタに保持しようとすると、2個のレジスタをペアで使うしか無い上、データバスも16ビットしかないので、データの転送にも2バスサイクルを要してしまいます。
　これでは、折角のFPUが泣きます。
　ところが、16ビットhalf型であれば、16ビット汎用レジスタにピタリ収まるし、16ビットバスでの1バスサイクル転送も可能で、とてもHappyなのです。

　以下、Moscovium用Half precision FPU(HFPU)の目標仕様とします。

　　目標仕様：

• IEEE 754 binary16をベースとする。
• ただし、非正規化数には対応しない。
• もし、非正規化数が入力された場合には、ゼロに縮退して演算することとする。
• INF(無限大)、NaN(非数)をサポートする。
• NaNは、half型バイナリ表現で0xfe00のみを発生させる。
• それ以外のNaN表現が入力された場合には、0xfe00に縮退して出力する。
• 加減乗除算を実行可能なFPUとする。
• Moscoviumのコプロセッサインターフェイスに接続するコプロセッサの一種としてHFPUを設計する。
• 浮動小数点レジスタは、Moscoviumの汎用レジスタと共用とする。

浮動小数点演算の算数

　以下、HFPUの設計に必要となる、浮動小数点演算の算数について、簡単にまとめておきます。以下、A，Bをhalf型の浮動小数点数値とします。

• 加算：A＋Bまたは、−(A＋B)
• half型は符号付数なので、上述の様に、結果として加算となるA＋Bと、−(A＋B)を加算と定義する。
• 加減算では、A≧0，B≧0として、符号部分だけ外に掃き出して表現している。
• 従って、−(A＋B)＝(−A)＋(−B)と言う意味で、Aの符号Sa＝1(マイナス)、Bの符号Sb＝1(マイナス)の加算と、理解すること。
• こうして整理すると、AとBが同符号の場合は加算、異符号の場合が減算と考えられる。
• 演算方法に関しては、10進数の演算を思い浮かべると理解がし易い。
• 1.23×10²＋2.34×10¹＝3.57×10³
• は、正しくない。
• 10の指数が異なるのに、仮数をそのまま足してしまった。：　1.23＋2.34＝3.57
• 正しくは、1.23×10²＋0.234×10²＝1.464×10²
• 桁(指数)を合わせてから加算する。
• half型も同じ。１）まず指数を合わせて、２）それから仮数を加算する、必要がある。
• X＝A＋B＝Fa×2^Ea＋Fb×2^Eb
• Ea≧Eb：　X＝(Fa＋Fb×2^Eb−Ea)×2^Ea
• Ea＜Eb：　X＝(Fa×2^Ea−Eb＋Fb)×2^Eb
• 2進数で2のべき乗を掛け算する事はシフト演算に置き換えることが出来るので、互いの指数(EaとEb)が等しくなるまで、仮数FaまたはFbをシフトして桁を合わせる。後は普通に足してしまって良い。
• 最後に符号を決める。A＋Bの場合も、−(A＋B)の場合も、結局、加算結果Xの符号Sx＝Saである。
• 減算：A−Bまたは、−(A−B)
• half型は符号付数なので、上述の様に、結果として減算となるA−Bと、−(A−B)を減算と定義する。
• 桁を合わせるところまでは加算と同じ。
• X＝A−B＝Fa×2^Ea−Fb×2^Eb
• Ea≧Eb：　X＝(Fa−Fb×2^Eb−Ea)×2^Ea
• Ea＜Eb：　X＝(Fa×2^Ea−Eb−Fb)×2^Eb
• この時、減算では、減算結果のFxが負数(2の補数表現)になってしまう場合がある。
• この時は、Fxの符号を反転して正の数に調整すると共に、符号Sxを反転する。
• 符号：
• Fx＜0だった場合：　Sx＝！Sa，Fx＝−Fx
• Fx≧0だった場合：　Sx＝Sa
• 乗算：X＝A×B
• これも、10進数の乗算を思い返すと、容易に解る。
• X＝A×B＝Fa×2^Ea×Fb×2^Eb
• X＝(Fa×Fb)×2^Ea＋Eb
• 仮数(Fa，Fb)同士は、11ビット符号無固定小数点乗算。
• 乗算結果は22ビットだが、小数点位置を挟む有効な11ビットを取り出して演算結果とする。
• 指数(Ea，Eb)は加算となる。
• Ex＝Ea＋Ebの値の範囲は−28≦Ex≦30だが、half型で表現可能な範囲である−14≦Ex≦15を超えた場合は、0またはINFとなる。
• 符号Sx＝Sa＾Sb(EOR)である。
• 除算：X＝A÷B
• 除算も同様。
• X＝A÷B＝Fa×2^Ea÷(Fb×2^Eb)
• X＝(Fa÷Fb)×2^Ea−Eb
• 仮数(Fa，Fb)は、符号無固定小数点除算。商として11ビットの精度を確保する必要があるので、この点を踏まえて除算を実行する。
• 指数(Ea，Eb)は減算となる。
• Ex＝Ea−Ebの値の範囲は−29≦Ex≦29だが、half型で表現可能な範囲である−14≦Ex≦15を超えた場合は、0またはINFとなる。
• 符号Sx＝Sa＾Sb(EOR)である。

　実は桁合わせ処理が必要になる分、浮動小数点演算では乗除算よりも加減算の方が少し処理が複雑となって、その分実行サイクル数もかかります。面白いですね。除算はステップ除算を使うと、そこそこサイクル数がかさみますから、恐らく一番速いのは乗算でしょう(加減算を抑えて)。
　しかも加減算には、「桁落ち」という演算精度を低下させる現象が起きるので、浮動小数点演算の中で加減算はとても曲者です。

非有限値の演算

　half型には、整数型には存在しなかったINFやNaN等の非有限値が存在するので、これらの演算にも対応する必要があります。

• NaNとの演算は全てNaN(0xfe00)。
• 以下、有限値A(≧0)とする。明示的にゼロ(0)との演算が記載されない場合、Aはゼロを含む。
• 加算：
• half型は符号付数なので、結果として加算となる、A＋Bおよび、−(A＋B)を加算とする(A≧0，B≧0)。
• INF＋INF＝INF
• INF＋A＝INF
• (−0)＋(−0)＝−0
• 0＋0＝0
• 加算は可換なので、被加数と加数を交換しても演算結果は同じ。
• 減算：
• half型は符号付数なので、結果として減算となる、A−Bおよび、−(A−B)を減算とする(A≧0，B≧0)。
• INF−INF＝NaN
• INF−A＝INF
• A−INF＝−INF
• A−A＝0
• (−A)−(−A)＝0
• 乗算：
• INF×INF＝INF
• INF×0＝NaN
• INF×A＝INF
• NaN以外の演算では、別途符号が計算される。
• 乗算は可換なので、被乗数と乗数を交換しても演算結果は同じ。
• 除算：
• INF÷INF＝NaN
• INF÷0＝INF
• INF÷A＝INF
• 0÷INF＝0
• 0÷0＝NaN
• 0÷A＝0
• A÷INF＝0
• A÷0＝INF
• NaN以外の演算では、別途符号が計算される。
• 全演算で、演算結果が有限値Aの絶対値最大値＝65,504を超えた場合はINFとなる。

　以上で全て網羅できていると思うし……、合っているとは思うけど……。
　数学の極限値計算(lim等)を踏襲した演算ルールみたいですが、なかなか難しくて正解が解りませんでした。CygwinのGCCでプログラムを組んで、実際に計算させてみた結果に基づいた、つ・も・り……(↓)。

　　#include　<stdio.h>
　　int　main (void)
　　{
　　　　float　inf=1./0.;
　　　　float　nan=0./0.;
　　　　printf ("%g %g %g\n",inf,nan, 123. - inf);
　　　　return (0);
　　}

　上コードは一例です。加減乗除算、±0、±A、±INF、NaNの組み合わせをすべて計算させて表示してみました。
　まあ、ここで設計したHFPUの演算結果が万一間違っていたとしても、as isで使ってください(無保証ですから〜)ｍ（１＿＿）ｍ。

浮動小数点演算回路の基本設計

　以下、前節の算数に基づいて、浮動小数点演算回路の設計を進めて行きます。

• 加減算前の桁合わせ回路：
• 加減算を実行する時には、最初にAとBの桁を合わせる必要がある。
• while (Eb＜Ea)｛Eb＋＋；Fb>>＝1；｝
• while (Ea＜Eb)｛Ea＋＋；Fa>>＝1；｝
• 上の処理は、指数EaとEbの内、より大きい方の指数により小さい方の指数を合わせて、仮数をシフトしている。このため、仮数のシフトは右シフトだけとなる。
• 一回、Ea−Ebを計算してしまうと、必要なシフト量はおのずと明らかなので、1サイクルでFxのシフトとExのインクリメントを終了することも可能。
• なんとなく、費用対効果(速度と論理規模のバランス)が合わない気がするため、順序回路で実装することに決める。
• つまり、今回のHFPUでは加減算前の桁合わせ(シフト)には複数サイクル掛かる。
• 加算回路：
• 最終的にAとBの符号を比較して、同符号(Sa＝Sb)の場合に加算を実行する。
• 以下例では、A≧0，B≧0とする。
• 例）AからＢを引く場合、A−Ｂ＝Ａ＋(−B)と見て、異符号→減算。
• 例）−AからBを引く場合、−Ａ−Ｂ＝(−A)＋(−B)と見て、同符号→加算。
• 例）Aと−Bを足す場合、Ａ＋(−B)なので、異符号→減算。
• 例）−Aと−Bを足す場合、(−Ａ)＋(−B)なので、同符号→加算。
• 以下、桁合わせ後の指数と仮数を使用する。
• Sx＝Sa；Ex＝Ea；Fx＝Fa＋Fb
• Norm(Ex，Fx)
• Norm()は、正規化処理。後で回路が登場する。
• 減算回路：
• 最終的にAとBの符号を比較して、異符号(Sa≠Sb)の場合に減算を実行する。
• 以下、桁合わせ後の指数と仮数を使用する。
• Sx＝Sa；Ex＝Ea；Fx＝Fa−Fb
• 仮数Fxの減算は、2の補数を使った符号付減算となる。……普通の符号付固定小数点減算。
• Fa＜Fbの場合、Fx＜0となる。この時は、Fx＝−Fx；Sx＝！Sx、とする。
• 更に減算時にはFx＝0となる場合がある。この場合はX＝0とする。
• Norm(Ex，Fx)
• 乗算回路：
• Sx＝Sa＾Sb；Ex＝Ea＋Eb；Fx＝Fa×Fb
• Norm(Ex，Fx)
• 0x1.000≦(Fa，Fb)≦0x1.ffcなので、0x1.00000≦Fa×Fb≦0x3.ff001が保証される。
• 仮数は11×11＝22ビット符号無固定小数点乗算となる。
• 演算結果22ビットから、固定小数点位置を挟んだ11ビット幅を積として出力する。
• 乗算回路はTD IDEのDSPマクロを使用(↓)。
  
  
  
  
• (整数型)乗算コプロセッサを作った時、16×16ビット符号付／符号無乗算を実行するのに使った乗算器マクロ(フラッシュ乗算器)は、17×17＝34ビット乗算器だったので、ずっと小さい規模で済む。
• 除算回路：
• Sx＝Sa＾Sb；Ex＝Ea−Eb；Fx＝Fa÷Fb
• Norm(Ex，Fx)
• 0x1.000≦(Fa，Fb)≦0x1.ffcなので、0x0.800≦Fa÷Fb≦0x1.ffcが保証される。
• 仮数は11ビットだが、除算コプロセッサと同様なハイブリッド除算器(4ビットフラッシュ除算器)を使うため、4の倍数に切り上げて24÷12ビット除算を実行し、12ビットの商を求める。
• (整数型)除算コプロセッサを作った時の16÷16ビット除算器は、実は内部では32÷16ビット除算を実行している。
• 被除数のビット数が除数のビット数の×2倍になるのはステップ除算器の特徴。
• こちらの図(ステップ除算器)で、Rem：Denレジスタ(16×2＝32ビット)をDsoレジスタ(16ビット)で割り算している状況が解る。
• 32÷16ビット除算器と比べて、若干だが小さい24÷12ビット除算器で浮動小数点除算を実行できる。論理規模→小と、演算速度→向上が期待できる。
• 商が12ビットなので、ハイブリッド除算器のステップ数は3ステップとなり、32ビット除算器の32/4＝8ステップと比べて半分以下のステップ数で終了できる。
• しかも、仮数は符号無固定小数点数なので、符号付除算の時の様な被除数・除数・商の符号反転などの処理も不要で、高速に演算可能。
• 正規化回路：
• 上の各演算回路でNorm()として参照してきた回路。
• half型の仮数Fxは、必ず、0x1.xxxの形になっている必要がある。
• 演算の結果、0x3.xxxとなって桁上げが有ったり、逆に桁借りによって0x0.xxxとなった場合、小数点上の1の位に"1"が立つ様に仮数Fxをシフトして、それに合わせて指数Exを増減する必要がある。
• この処理を正規化と呼ぶ。
• while (Fx≧0x2.000)｛Fx>>＝1；Ex＋＋；｝
• Ex＞15となったら、オーバフロウ。X＝±INF。
• while (Fx＜0x1.000))｛Fx<<＝1；Ex−−；｝
• Ex＜−14となったら、アンダフロウ。X＝±0。
• こちらも、加減算前のシフト処理と同様、今回は順序回路で組んでおく事とする。
• 演算性能を追求する場合は、場合分けして1サイクルで正規化処理を終了することも可能だが、これもまた費用対効果が合わないと判断した。

ブロック図と、各ブロックの機能設計

　ここまでの検討に基づいて作成したブロック図は以下の通りです。



　以下、各ブロックの機能設計を進めます。

• hlfa，hlfb：　浮動小数点レジスタA，B
• CPUとのコプロセッサI/F(Aバス、Bバス)からデータを受け取るレジスタ。
• 内部表現である、Ex＝Ex’−15，Fx＝0x1.<Fx’>に変換して、hadd／hmul／hdiv演算器への入力バスへ出力する。
• 加減算実行時の桁合わせ(シフト)回路を含む。
• シフト回路を使って、half型→uint16型変換も実施できるようにしておく(hlfaブロック)。
• 減算コマンド(HSUB)だった場合、hlfbブロックの出力データにおいて、符号(Sb)を反転して出力する。
• hadd：　加減算回路
• 加減算と符号処理を行う。
• A，B入力の符号を判定して、同符号(Sa＝Sb)なら加算、異符号(Sa≠Sb)なら減算を実施する。
• Fa，Fbは符号無固定小数点数だが、減算があるので、符号ビットを付加して2の補数による符号付加減算(Fx＝Fa±Fb)を実施する。
• 減算した結果、Fx＜0(2の補数)となった場合、Fx＝−Fx、Sx＝！Sxとして出力する。
• 非有限数(INFやNaN等)との加減算も、このブロックで実行する。
• hmul：　乗算回路
• 符号(Sa，Sb)と、指数(Ea，Eb)、仮数(Fa，Fb)の乗算を実行する。
• 符号の乗算Sx＝Sa＾Sb。
• 指数の乗算Ex＝Ea＋Eb。
• 仮数の乗算Fx＝Fa×Fbは、11×11＝22ビット符号無固定小数点乗算となる。
• 演算結果22ビットから、固定小数点位置を挟んだ11ビット幅を積として出力する。
• 乗算自体は、TD IDEのDSPマクロを使用する。
• 非有限数(INFやNaN等)との乗算も、このブロックで実行する。
• hdiv：　除算回路
• 符号(Sa，Sb)と、指数(Ea，Eb)、仮数(Fa，Fb)の除算を実行する。
• 符号の除算Sx＝Sa＾Sb。
• 指数の除算Ex＝Ea−Eb。
• 仮数の除算Fx＝Fa÷Fbは、24÷12ビット符号無除算で、ハイブリッド除算器を使用する。
• 非有限数(INFやNaN等)との除算も、このブロックで実行する。
• norm：　正規化回路。
• 各演算器(hadd／hmul／hdiv)からの演算結果出力を正規化してCバスへ出力する。
• オーバフロウでINF、アンダフロウで0をはじめ、各演算器での非有限値演算結果を受けて、非有限値の出力も担う。
• half型→uint16型変換等の変換機能(処理はhlfaブロック)も、このブロックを通ってCバス出力となる。
• uint16型→half型変換機能に関しては、正規化処理そのもの、なので、本ブロックが実行する。
• hctl：　制御ステートマシン
• HFPU全体を制御するステートマシン。
• CPUからのコプロセッサコマンドccmd[4:0]を受けて、各演算を制御実行する。
• 加減算の桁合わせ、除算のステップ、正規化のシフトなど、全てを制御する。
• CPUとのコプロセッサインターフェイスに対して、レディ信号crdyを出力し、HFPU内部処理が終了するまでCPUにウェイトをかける。

HFPU用コプロセッサコマンド

　CPUからのコマンド入力(ccmd[4:0])は、他のコプロセッサ(MULC16、DIVC32)と共用なので、こちらのコプロセッサコマンドと競合しないよう、以下表の通り、コマンドを割り当てます。

　四則演算に留まらず、整数⇔half型浮動小数点数変換など、最低限必要な機能を割り当ててあります。


　今回はここまでで一旦閉じて、次回動作検証してHFPUとそれを搭載したMoscovium CPUを完成させます💥。