浮動小数点数演算器のつくりかた(その1)

最近浮動小数点数演算器を作っているので、作るときのポイントを忘れないようにメモしておきたいと思います。

以下の全4回で書いていく予定です。

  1. 最近接丸めの浮動小数点数加算器の作り方(今回)
  2. 最近接丸めの浮動小数点数乗算器の作り方
  3. RISC-Vの規格に従った(つまり、IEEE754規格にも従った)浮動小数点数積和演算器の作り方
  4. RISC-Vの規格に従った浮動小数点数除算・平方根両対応演算器の作り方

コード

以下の74行のVerilog HDLコードで、x86の行う最近接丸めの浮動小数点数加算と同じ結果が得られるはずです。

メモ:これを書いたのは2022/10/06で、3cycleパイプライン化したものがRSDで採用されています。

module fp_adder(lhs, rhs, result);
    input  wire[31:0] lhs;
    input  wire[31:0] rhs;
    output wire[31:0] result;

    wire       lhs_sign = lhs[31];
    wire       rhs_sign = rhs[31];
    wire [7:0] lhs_expo = lhs[30:23];
    wire [7:0] rhs_expo = rhs[30:23];
    wire[22:0] lhs_mant = lhs[22:0];
    wire[22:0] rhs_mant = rhs[22:0];

    // NaN handling
    wire lhs_is_nan = lhs_expo == 8'hff & lhs_mant != 0;
    wire rhs_is_nan = rhs_expo == 8'hff & rhs_mant != 0;
    wire lhs_is_inf = lhs_expo == 8'hff & lhs_mant == 0;
    wire rhs_is_inf = rhs_expo == 8'hff & rhs_mant == 0;
    wire res_is_nan = lhs_is_nan | rhs_is_nan | (lhs_sign != rhs_sign & lhs_is_inf & rhs_is_inf);
    wire[31:0]  nan = lhs_is_nan ? lhs | 32'h00400000 : rhs_is_nan ? rhs | 32'h00400000 : 32'hffc00000; // qNan

    // Preparation
    wire swap_is_needed   = lhs[30:0] < rhs[30:0];
    wire is_subtraction   = lhs_sign != rhs_sign;
    wire[31:0] s_lhs      = swap_is_needed ? rhs : lhs;
    wire[31:0] s_rhs      = swap_is_needed ? lhs : rhs;
    wire [7:0] s_lhs_expo = s_lhs[30:23];
    wire [7:0] s_rhs_expo = s_rhs[30:23];
    wire [7:0] s_lhs_offs = s_lhs_expo == 0 ? 1 : s_lhs_expo;
    wire [7:0] s_rhs_offs = s_rhs_expo == 0 ? 1 : s_rhs_expo;
    wire[23:0] s_lhs_mant = { s_lhs_expo != 8'h00, s_lhs[22:0] }; // s_lhs_expo != 8'h00 is the hidden bit of a normalized number
    wire[23:0] s_rhs_mant = { s_rhs_expo != 8'h00, s_rhs[22:0] }; // s_rhs_expo != 8'h00 is the hidden bit of a normalized number
    wire       prec_loss  = is_subtraction & s_lhs_expo - s_rhs_expo <= 1 & s_lhs_expo != 8'hff;

    // When precision loss does not occur
    wire[26:0] adder_lhs   = is_subtraction ? { s_lhs_mant, 1'b0, 2'b0 } : { 1'b0, s_lhs_mant, 2'b0 }; // round logic commonizing trick: (3 bits for guard+round+sticky) or (1 bit for carry and 2 bits for guard+sticky)
    wire       large_diff  = s_lhs_expo - s_rhs_expo > 31;
    wire [4:0] offs_diff   = large_diff ? 31 : s_lhs_offs[4:0] - s_rhs_offs[4:0];
    wire[48:0] shifted_rhs = { s_rhs_mant, 25'h0 } >> offs_diff >> !is_subtraction; // when offs_diff = 24 addition, the hidden bit should be shifted to the guard bit position. Example: 0x1.p+24 + 0x1.p+0 => GS=10, 0x1.p+24 + 0x1.000002p+0 => GS=11
    wire[26:0] adder_rhs   = { shifted_rhs[48:24], shifted_rhs[23], shifted_rhs[22:0] != 0 }; // Last 2 bits are the guard bit and the sticky bit. Note: the sticky bit may be incorrect but is not used when |rhs| < 0.5ULP(|lhs|).

    wire[26:0] adder_result  = is_subtraction ? adder_lhs - adder_rhs : adder_lhs + adder_rhs;
    wire       round_away_a  = adder_result[26] ? adder_result[2] & (adder_result[3] | adder_result[1] | adder_result[0])
                                                : adder_result[1] & (adder_result[2] |                   adder_result[0]); // round to nearest, ties to even
    wire       exp_plus_one  = (s_lhs[30:23] == 8'h00 & adder_result[25]) | adder_result >= 27'h3fffffe; // when the sum of two subnormal numbers is a normal number or a carry is generated with rounding taken into account

    wire[22:0] final_mant_a  = (adder_result[26] ? adder_result[25:3] : adder_result[24:2]) + { 22'h0, round_away_a }; // No special treatment is required even if an overflow occurs since the answer will be 0 and it will be correct.
    wire [7:0] final_expo_a  = s_lhs_expo + { 7'h0, exp_plus_one } - { 7'h0, is_subtraction }; // No overflow occurs because (is_subtraction ? 2 : 0) <= s_lhs_expo <= 254 holds.
    wire       res_is_inf    = s_lhs_expo == 8'hff | final_expo_a == 8'hff;
    wire[31:0] inf           = { s_lhs[31], 8'hff, 23'h0 };

    wire[31:0] final_result_a = res_is_inf ? inf : { s_lhs[31], final_expo_a, final_mant_a };

    // When precision loss occurs
    wire[24:0] suber_lhs    = { s_lhs_mant, 1'b0 };
    wire[24:0] suber_rhs    = { s_rhs_mant, 1'b0 } >> (s_lhs_offs[0] != s_rhs_offs[0]); // At most, 1-bit shift. s_lhs_offs[0] != s_rhs_offs[0] is equal to s_lhs_expo - s_rhs_expo because s_lhs_expo - s_rhs_expo <= 1.
    wire[24:0] suber_result = suber_lhs - suber_rhs;

    reg  [4:0] lzc; always@* for(lzc = 0; lzc <= 24; lzc = lzc+1) if( suber_result[24-lzc] ) break; // leading zeros count (however, if suber_result == 0, i is 24, not 25.)
    wire       round_away_s = suber_result[1] & suber_result[0]; // round to nearest, ties to even
    wire       subnormal    = { 3'b0, lzc } >= s_lhs_offs;

    wire[22:0] final_mant_s = s_lhs_offs == 1 ? suber_result[23:1] :
                              subnormal       ? suber_result[22:0] << (s_lhs_offs[4:0]-5'd2) :
                              lzc == 0        ? suber_result[23:1] + { 22'h0, round_away_s } // never overflows
                                              : suber_result[22:0] << (lzc-5'd1);
    wire [7:0] final_expo_s = subnormal ? 0 : s_lhs_offs - { 3'b0, lzc };
    wire       res_is_zero  = suber_result == 0;
    wire[31:0] zero         = 32'h00000000;

    wire[31:0] final_result_s = res_is_zero ? zero : { s_lhs[31], final_expo_s, final_mant_s };

    assign result = res_is_nan ? nan :
                    prec_loss  ? final_result_s : final_result_a;
endmodule

コードの方針

Two path加算器です。 つまり、「シフトしてから加減算する回路」と「減算してからシフトする回路」の両方があるような構成です。 前者は指数部が遠いときに選ばれるのでfar pathと呼ばれます。 後者は指数部が近いときに選ばれるのでclose pathと呼ばれます。

大まかなコードはの構成は以下のようになっています。

  1. NaNを返すべき例外ケースを判別する。
  2. 絶対値の大きい方をs_lhs、絶対値の小さい方をs_rhsに取り直す。
  3. (far path) 実質的な加算または指数部の差が2以上の実質的な減算の場合、シフトしてから加減算する。シフトの時に追い出されたビットはsticky bitとしてまとめる。加算の時と減算の時で1 bitずらすことで、round bitを不要とする。
  4. (close path) 指数部の差が1以内の実質的な引き算の場合、引き算してからleading zeros countを行い、その分だけシフトする。丸めはほとんど自明である。
  5. 終結果を、NaN、Inf、far pathの結果、+0.0、close pathの結果、のいずれかから選ぶ。

コードの詳細

コードのスタイルについて

N_FLOAT_MANTISSA_BITSみたいな読みづらい上に何を意味しているかわかりづらい定数を使わず、マジックナンバーを直書きしているのがこだわりポイントです。 仮数部のビット数+αみたいなビット数のwireを定義する場所が多いため、定数を使ってしまうとものすごく読みにくくなります。 また、「仮数部のビット数」と言われてもhidden bitが含まれているのかいないのかは明確でないと思います。

浮動小数点数演算器のHDLを読めるプログラマであれば[22:0]が何を意味しているかは一目瞭然です。 また逆に、どんなに浮動小数点数演算器のHDLに詳しかったとしても、wire [N_FLOAT_MANTISSA_BITS + 3 - 1 : 0]などと書かれているコードを見たとき、なぜ3でよいのかをすぐに理解するのは困難でしょう。 よってその近くにはなぜ余分に3 bit必要なのかがコメントが書かれているはずです。 その場合、wire [25:0]と書かれていても、それが何を意味しているかは明白です。 それなら、ぱっと見で把握しやすくてわかりやすい表記にすべきです。 何も考えずにマジックナンバーを排除するのではなく、何が読みやすいのかを考えてコードを書きましょう。

また、[25:3][3+:23]のように書くことはしませんでした。 [25:3]などと書いたほうが、どのビットを使用済かが明確にわかるためです(近くに[25]が出てきていたら1 bitずれている可能性が高い、と判断できます)。 [3+:23]の方が23 bitの信号線であることが明白という考え方もありそうで、このあたりは好みが分かれるところかもしれません。

準備

説明不要だと思います。

module fp_adder(lhs, rhs, result);
    input  wire[31:0] lhs;
    input  wire[31:0] rhs;
    output wire[31:0] result;

    wire       lhs_sign = lhs[31];
    wire       rhs_sign = rhs[31];
    wire [7:0] lhs_expo = lhs[30:23];
    wire [7:0] rhs_expo = rhs[30:23];
    wire[22:0] lhs_mant = lhs[22:0];
    wire[22:0] rhs_mant = rhs[22:0];

NaNを返すべき例外ケースの判別

これもほとんど明らかだと思います。

返すNaNは、x86の挙動に合わせて左オペランド優先で伝搬させています。 両オペランドがNaNではないのにNaNを返すべきケース、つまり無限大同士の引き算の場合、x86の挙動にあわせて0xffc00000(符号部が負で指数部が最大値で仮数部の先頭に1が立っている(=quietな)NaN)を返しています。

RISC-Vの仕様では、常に0x7fc00000canonical NaN。符号部が正で指数部が最大値で仮数部の先頭に1が立っている(=quietな)NaN)を返せばよいようです。

    // NaN handling
    wire lhs_is_nan = lhs_expo == 8'hff & lhs_mant != 0;
    wire rhs_is_nan = rhs_expo == 8'hff & rhs_mant != 0;
    wire lhs_is_inf = lhs_expo == 8'hff & lhs_mant == 0;
    wire rhs_is_inf = rhs_expo == 8'hff & rhs_mant == 0;
    wire res_is_nan = lhs_is_nan | rhs_is_nan | (lhs_sign != rhs_sign & lhs_is_inf & rhs_is_inf);
    wire[31:0]  nan = lhs_is_nan ? lhs | 32'h00400000 : rhs_is_nan ? rhs | 32'h00400000 : 32'hffc00000; // qNan

計算の準備

用語の定義

  • ビット表現に含まれる指数部:浮動小数点数のビット表現の「指数部」とされている部分のビット表現。下駄履きしていない値。単精度浮動小数点数だと8 bit分。
  • ビット表現に含まれる仮数部:浮動小数点数のビット表現の「仮数部」とされている部分のビット表現。hidden bitを含まない。単精度浮動小数点数だと23 bit分。
  • 実質的な指数部:ビット表現に含まれる指数部が非零(正規化数)であればその値、零(非正規化数)であれば1。
  • 実質的な仮数部:ビット表現に含まれる仮数部と、hidden bit(0または1)を連結した値。
  • 実質的な減算:true subtraction。同符号数同士の引き算または異符号数同士の足し算のこと。
  • 実質的な加算:true addition。同符号数同士の足し算または異符号数同士の引き算のこと。

コードの説明

引き算の結果が負になったりするとややこしいので、わかりやすくオペランドを交換します。 どちらの絶対値が大きいかを調べるためには、符号部は無視して、指数部と仮数部を辞書式順序で比較すればよいです。 31 bit整数の比較を行うのはいかにも遅そうでためらわれますが、いずれにせよ指数部の差の計算(8 bit整数の引き算)は行わないといけないので、それと比べるとそれほど遅いというわけではありません。 ※David Lutz, "ARM Floating Point 2019: Latency, Area, Power", 2019 IEEE 26th Symposium on Computer Arithmetic (ARITH), 26, pp. 97–98, 2019. によれば、絶対値の比較をせず、しかも指数部の引き算もせず、a-b, a-(b>>1), b-a, b-(a>>1)の4通りの計算をいきなり始めてしまって、後で選択するという速度に特化した実装がARMプロセッサでは採用されているようです……。

次に、非正規化数への対処を行います。 非正規化数は、実質的な指数部(offs)が1で、hidden bit込みの仮数部(mant)の先頭(つまりhidden bit)が0であるとして取り扱えばよいです。 逆に、正規化数の場合は、実質的な指数部(offs)はビット表現に含まれる指数部(expo)そのままで、hidden bit込みの仮数部(mant)は1とビット表現に含まれる仮数部を連結したものとして取り扱います。

最後に、close pathを使うかfar pathを使うかを選択します。 実質的な引き算で、かつ指数部の差が1以下であればclose pathを使います。 それ以外の場合はfar pathを使います(指数部が近くても実質的な加算であれば使うのはfar pathであるという点が紛らわしいので注意してください)。

なお、交換後の左オペランドが無限大の場合は、指数部の差が1以下の実質的な引き算であってもfar pathを使います。 これは、far path側に無限大判定回路と無限大出力回路があるため、それを流用したほうが楽だからです。

    // Preparation
    wire swap_is_needed   = lhs[30:0] < rhs[30:0];
    wire is_subtraction   = lhs_sign != rhs_sign;
    wire[31:0] s_lhs      = swap_is_needed ? rhs : lhs;
    wire[31:0] s_rhs      = swap_is_needed ? lhs : rhs;
    wire [7:0] s_lhs_expo = s_lhs[30:23];
    wire [7:0] s_rhs_expo = s_rhs[30:23];
    wire [7:0] s_lhs_offs = s_lhs_expo == 0 ? 1 : s_lhs_expo;
    wire [7:0] s_rhs_offs = s_rhs_expo == 0 ? 1 : s_rhs_expo;
    wire[23:0] s_lhs_mant = { s_lhs_expo != 8'h00, s_lhs[22:0] }; // s_lhs_expo != 8'h00 is the hidden bit of a normalized number
    wire[23:0] s_rhs_mant = { s_rhs_expo != 8'h00, s_rhs[22:0] }; // s_rhs_expo != 8'h00 is the hidden bit of a normalized number
    wire       prec_loss  = is_subtraction & s_lhs_expo - s_rhs_expo <= 1 & s_lhs_expo != 8'hff;

Far path

まず、実質的な指数部(offs)の差を調べます。 ただし、ビット表現に含まれる指数部(expo)の差が31超の場合は、31にクランプします。 ここは一貫していないように見えますが、以下の二点によるものです。

  • s_lhs_expo - s_rhs_expoを計算する回路はすでに存在している(close pathを使う条件の一つである「指数部の差が1以下」を判定するために必要)
  • クランプする値が31である必要は特になく、十分多く右シフトされればいいだけ。±1の誤差にこだわる必要はない。シフト回路が大きくならないように5 bitに収めるのが重要

その後、実質的な指数部の差に従って、右オペランドをシフトします。 この時、is_subtractionによってシフト量が違いますが、これはadder_lhsの位置に合わせるためです。 adder_lhsは、引き算の時はs_lhs_mantを左端にそろえていますが、足し算の時はs_lhs_mantを左端から1 bit空けた位置に配置しています。 これは、以下の理由によります(浮動小数点数演算における保護桁について - よーるで説明した方法です)。

  • 実質的な加算の時は繰り上がりがありえるものの、繰り下がりはあり得ない。よって上位に1 bitの余分な空間が必要な一方、下位にはguard bitとsticky bitだけで十分で、round bitは不要
  • 実質的な減算の時は繰り下がりがありえるものの、繰り上がりはあり得ない。よって上位はギリギリに詰めていい一方、下位にはguard bitとsticky bitに加えてround bitが必要

実質的な加算の場合と実質的な減算の場合で1 bitずらして配置することで、丸めの回路をこの二つで共通化することができます。

Sticky bitは、シフトしてはみ出た部分全ての論理和として作成します。 たくさんシフトするとはみ出た部分がsticky bitに回収されずに虚空へ消えますが、そもそもそんなに大きくシフトされる場合は0.5 ULPに満たないので、最近接丸めを実装する上では問題ありません。 Guard bitが立ちうる限界は、hidden bit込みの仮数部の最上位ビット(つまりhidden bitの位置)がguard bitの位置までシフトされた時です。 よって、sticky bitを計算するために必要なのはhidden bit抜きの仮数部のビット数、つまり23 bitです。

次に、足し算または引き算を実行します。 繰り上がりや繰り下がりの発生は、adder_result[26](最上位ビット)を見ればわかります。 最上位ビットが1であればそれは繰り上がりが発生した・繰り下がりが発生しなかったことを意味します。 最上位ビットが0であればそれは繰り上がりが発生しなかった・繰り下がりが発生したことを意味します。 最上位ビットが1の場合、adder_result[3]仮数部の最下位のビットになります(丸めたことによる指数部の繰り上がりが発生しない場合)。その下のadder_result[2]がguard bit、|adder_result[1:0]がsticky bit相当です。 最上位ビットが0の場合、far pathを使っている仮定(※)から、adder_result[2]仮数部の最下位のビットになります(同)。その下のadder_result[1]がguard bit、adder_result[0]がsticky bitです。

※実質的な加算であれば、adder_lhs仮数部の最下位の位置より低くなることはありえません。実質的な減算の場合、far pathを使うということは実質的な指数部の差が2以上あるということであり、最大で1桁の繰り下がりしか生じません。よって、最上位ビットが0だった場合、adder_lhs仮数部の最下位の位置より1つ下が計算結果の最下位の位置です。ビット表現に含まれる指数部の差が2以上あるということは左オペランドの実質的な指数部は2以上であり、1桁の繰り下がりを生じても非正規化数になることはないため、最下位の位置の計算が狂うことはない点にも注意してください。

これを使って、round to nearest, ties to evenで丸めた場合、切り上げか切り捨てかを判定します。 より正確には、「絶対値が大きくなる方に丸める」か「絶対値が小さくなる方に丸める」のどちらであるかを判定します(round_away_a)。 絶対値が小さくなる方に丸めるのは簡単で、単に下位ビットを削ればいいだけです。 絶対値が大きくなる方に丸めるためには、適切な位置に1を足す必要があります。

また、指数部に繰り上がりが生じるかも判定します。 adder_result >= 27'h3fffffeという式で判定できます(adder_resultの最下位はsticky bitですが、繰り上がる側が“偶数”になっているので、sticky bitが立っていなくても繰り上がります)。 これは丸めを考慮した上で指数部が繰り上がるかを判定する式になっているので、round_away_aを追加で考慮に入れる必要はありません。 なお、2桁の繰上りが生じることはありえません。

続いて、指数部と仮数部を求めます。 仮数部のビット表現は、(adder_result[26] ? adder_result[25:3] : adder_result[24:2]) + { 22'h0, round_to_away }で求められます。 指数部のビット表現は、s_lhs_expo + { 7'h0, exp_plus_one } - { 7'h0, is_subtraction }で求められます。 仮数部を求める式ではオーバーフローが発生し得ますが、オーバーフローが発生するのは指数部が繰り上がって仮数部が23'h0になるべき時なので、結果的にOKです。 指数部を求める式で- { 7'h0, is_subtraction }としているのは、実質的な加算の時の桁をそろえるために、繰り下がりが発生することがデフォルトであると取り扱っているためです。

最後に、指数部が8'hffになっていれば、それはオーバーフロー(計算結果の絶対値が0x1.ffffffp+127以上=0x1.fffffep+127よりも0x1.p+128の方が近い)ということなので、左オペランドと同じ符号を持つ無限大を返すようにしておきます。 ついでに、この回路を流用して、左オペランドが元々無限大であったときも、左オペランドと同じ符号を持つ無限大を返すようにしておきました。

なお、このままでは非正規化数の取り扱いが微妙に正しくありません。 「指数部が繰り上がる」の判定に「左オペランドが非正規化数かつadder_result[25]1になっている」を加えれば正しくなります。 非正規化数同士の加算結果が0x1.p-126を超えない時、丸めが生じないため(※)、特に丸めの結果繰り上がることもなく、adder_result[25]を見るだけで判定できます。

※非正規化数は(というか単精度浮動小数点数全てがそうですが)全て0x1.p-149の倍数である点に注意します。

    // When precision loss does not occur
    wire[26:0] adder_lhs   = is_subtraction ? { s_lhs_mant, 1'b0, 2'b0 } : { 1'b0, s_lhs_mant, 2'b0 }; // round logic commonizing trick: (3 bits for guard+round+sticky) or (1 bit for carry and 2 bits for guard+sticky)
    wire       large_diff  = s_lhs_expo - s_rhs_expo > 31;
    wire [4:0] offs_diff   = large_diff ? 31 : s_lhs_offs[4:0] - s_rhs_offs[4:0];
    wire[48:0] shifted_rhs = { s_rhs_mant, 25'h0 } >> offs_diff >> !is_subtraction; // when offs_diff = 24 addition, the hidden bit should be shifted to the guard bit position. Example: 0x1.p+24 + 0x1.p+0 => GS=10, 0x1.p+24 + 0x1.000002p+0 => GS=11
    wire[26:0] adder_rhs   = { shifted_rhs[48:24], shifted_rhs[23], shifted_rhs[22:0] != 0 }; // Last 2 bits are the guard bit and the sticky bit. Note: the sticky bit may be incorrect but is not used when |rhs| < 0.5ULP(|lhs|).

    wire[26:0] adder_result  = is_subtraction ? adder_lhs - adder_rhs : adder_lhs + adder_rhs;
    wire       round_away_a  = adder_result[26] ? adder_result[2] & (adder_result[3] | adder_result[1] | adder_result[0])
                                                : adder_result[1] & (adder_result[2] |                   adder_result[0]); // round to nearest, ties to even
    wire       exp_plus_one  = (s_lhs[30:23] == 8'h00 & adder_result[25]) | adder_result >= 27'h3fffffe; // when the sum of two subnormal numbers is a normal number or a carry is generated with rounding taken into account

    wire[22:0] final_mant_a  = (adder_result[26] ? adder_result[25:3] : adder_result[24:2]) + { 22'h0, round_away_a }; // No special treatment is required even if an overflow occurs since the answer will be 0 and it will be correct.
    wire [7:0] final_expo_a  = s_lhs_expo + { 7'h0, exp_plus_one } - { 7'h0, is_subtraction }; // No overflow occurs because (is_subtraction ? 2 : 0) <= s_lhs_expo <= 254 holds.
    wire       res_is_inf    = s_lhs_expo == 8'hff | final_expo_a == 8'hff;
    wire[31:0] inf           = { s_lhs[31], 8'hff, 23'h0 };

    wire[31:0] final_result_a = res_is_inf ? inf : { s_lhs[31], final_expo_a, final_mant_a };

Close path

指数部の差は、s_lhs_offs[0] != s_rhs_offs[0]で調べることができます。 なぜなら、close pathが使われるのは、指数部の差が1以下の時だからです。 これを使って右オペランドをシフトし、左オペランドとの差を計算します。

次に、引き算の結果のleading zeros countを求めます。 引き算の結果が0の場合は特別扱いされるため、0の場合のleading zeros countはドントケアであり、何が返ってきてもよいです。

指数部の差が1以下なので、計算結果には0.5 ULP単位の端数しか出ません。 よって、丸めの判定は単純にround_away_s = suber_result[1] & suber_result[0]で行えます(suber_result[1]仮数部の最下位、suber_result[0]がguard bit相当です)。

終結果の仮数部は、以下の5通りの中から選ぶことになります。

  1. 引き算の結果が0になった場合
    • IEEE754規格で、符号部が正のゼロを返すことに決まっている。
  2. オペランドの実質的な指数部が1の時
    • 結果の実質的な指数部も1になるはずなので、suber_result[23:1]がビット表現に含まれる仮数部となる。suber_result[24]1であれば正規化数、0であれば非正規化数である。
    • 繰り下がりが起こらなかったパターンの亜種。
  3. オペランドの実質的な指数部は2以上だったのに結果が非正規化数になる場合
    • 少なくとも2桁は繰り下がりが発生している。leading zeros countではなく左オペランドの実質的な指数部に応じた数だけ左シフトする必要があり、suber_result[22:0] << (s_lhs_offs-2)がビット表現に含まれる仮数部となる。
  4. それ以外で、繰り下がりが起こらなかった場合
    • suber_result[23:1] + { 22'h0, round_away_s }がビット表現に含まれる仮数部となる。
    • suber_result[24]1のときにsuber_result[23:0]24'hffffffになることはないため、ここでオーバーフローすることはない。
  5. それ以外の場合(繰り下がりが起きて、結果が正規化数になる場合)
    • suber_result[22:0] << (lzc-1)がビット表現に含まれる仮数部となる。

終結果の指数部は、非正規化数なら8'h00とし、正規化数なら左オペランドの実質的な指数部からlzcだけ引いた数とすればOKです。 もちろん、ゼロを返す場合の指数部は8'h00です。

    // When precision loss occurs
    wire[24:0] suber_lhs    = { s_lhs_mant, 1'b0 };
    wire[24:0] suber_rhs    = { s_rhs_mant, 1'b0 } >> (s_lhs_offs[0] != s_rhs_offs[0]); // At most, 1-bit shift. s_lhs_offs[0] != s_rhs_offs[0] is equal to s_lhs_expo - s_rhs_expo because s_lhs_expo - s_rhs_expo <= 1.
    wire[24:0] suber_result = suber_lhs - suber_rhs;

    reg  [4:0] lzc; always@* for(lzc = 0; lzc <= 24; lzc = lzc+1) if( suber_result[24-lzc] ) break; // leading zeros count (however, if suber_result == 0, lzc is 24, not 25.)
    wire       round_away_s = suber_result[1] & suber_result[0]; // round to nearest, ties to even
    wire       subnormal    = { 3'b0, lzc } >= s_lhs_offs;

    wire[22:0] final_mant_s = s_lhs_offs == 1 ? suber_result[23:1] :
                              subnormal       ? suber_result[22:0] << (s_lhs_offs[4:0]-5'd2) :
                              lzc == 0        ? suber_result[23:1] + { 22'h0, round_away_s } // never overflows
                                              : suber_result[22:0] << (lzc-5'd1);
    wire [7:0] final_expo_s = subnormal ? 0 : s_lhs_offs - { 3'b0, lzc };
    wire       res_is_zero  = suber_result == 0;
    wire[31:0] zero         = 32'h00000000;

    wire[31:0] final_result_s = res_is_zero ? zero : { s_lhs[31], final_expo_s, final_mant_s };

終結果の選択

NaNを出力すべき場合はnanを、close pathの結果を使う場合はfinal_result_sを、far pathの結果を使う場合はfinal_result_aを返すだけです。

    assign result = res_is_nan ? nan :
                    prec_loss  ? final_result_s : final_result_a;
endmodule

確認

入力が共に「ビット表現に含まれる指数部が0,1,2,3,21,22,23,24,25,26,27,128,129,130,250,251,252,253,254,255、ビット表現に含まれる仮数部が0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,32,64,98,99,100,101,102,128,256,512,1024,2048,4096,8192,16384,32768,65536,13072,262144,0x7fff000,0x7fff800,0x7fffc00,0x7ffe00,0x7fff00,0x7fff80,0x7fffc0,0x7fffe0,0x7ffff0,0x7ffff1,0x7ffff2,0x7ffff3,0x7ffff4,0x7ffff5,0x7ffff6,0x7ffff7,0x7ffff8,0x7ffff9,0x7ffffa,0x7ffffb,0x7ffffc,0x7ffffd,0x7ffffe,0x7fffff、符号部が正か負」であるような計算の実行結果がx86とビットレベルで一致することを確認しました。 この中には、正負のゼロや無限大、NaN、非正規化数などが含まれており、多くのコーナーケースをカバーしているはずです。

また、疑似乱数生成器std::mt19937で生成した 6.4×10^{12}個の入力に対する実行結果もx86とビットレベルで一致することを確認しました。

タイミング

RSDの3cycleパイプライン方式を踏まえてどのくらいの遅延なのかを書き出してみます。

計算の準備

31 bit比較の後、その結果に応じて選択した結果を使った8 bit減算・比較を行い、その後多少の論理演算を行う部分がクリティカルパスです。 XilinxFPGAなら、クリティカルパスはCARRY4x8+(LUT+CARRY4x2)+LUTになると思います。

なお、31bit比較の後、その結果に応じて選択する部分のみが真のクリティカルパスに乗っており、その後の8 bit減算・比較・論理演算はクリティカルパスには乗っていません。 prec_lossの計算は、フォワーディングなどの関係で第一サイクルに計算するのが苦しければ、後ろ倒しすることが可能です。 その場合、クリティカルパスはCARRY4x8+LUTになりますが、lhs_expo == 0rhs_expo == 0を明示的に事前に計算しておく必要があるかもしれません。

Far path

前半

最大31 bit分のシフトを行った後、sticky bitを作るために23 bitの論理和リダクションを行う部分がクリティカルパスです。 XilinxFPGAなら、クリティカルパスはLUTが5段分だと思います。 16 bitシフト、8/4 bitシフト、2/1 bitシフト、の順でやればLUT 4段分にできそうです。

後半

27 bitの加減算を行い、繰り上がりに応じて結果の選択を行った後、丸め方向に応じて23 bit整数に0か1を足す部分がクリティカルパスです。 一応その後にInfを選ぶかの選択が一段ありますが、これは最終段階での選択のLUTとマージできそうです。 XilinxFPGAなら、クリティカルパスはCARRY4x7+LUT+CARRY4x6になると思います。

Close path

前半

最大1 bitのシフトを行った後、25 bitの加算を行い、leading zeros countする部分がクリティカルパスです。 XilinxFPGAなら、クリティカルパスは(LUT+CARRY4x6)+LUT 4段になると思います(多分)。

後半

シフト量に定数を加算した後、最大31 bitのシフトを行った後、5つから選択する部分がクリティカルパスです。 XilinxFPGAなら、クリティカルパスはLUT 6段になると思います(多分)。 シフト量への定数加算はパディングを入れればノーコストにできます。 5つから選択する部分は早く計算される3つを先に選択しておき、遅く計算されるシフト系の2つと合わせて3つの中から選択、のような構成にすればクリティカルパス上のLUTは1段で済ませることができます。 よって、合わせてLUT 4段で作ることは可能だと思います。

ちゃんとシフタを共有する構成にすれば、最後のLUTを最終段階の選択のLUTとマージすることも可能で、LUT 3段まで縮められると思います。 シフト量を選択する部分が増えているように見えますが、シフタの初段の2:1セレクタ部分が3入力LUTから5入力LUTになるだけなので段数は変わりません。 パディングの量が異なるのは問題になりますが、1 bitシフトを最初にやるなら5入力LUTが6入力LUTになるだけですし、そもそもlzc1だけずらした定義にするという方法もあります。

シフトが必要な場合は丸めが不要なので、シフトと丸めの両方がクリティカルパスに乗ることはありません。 また、一般にシフトは加算より重い処理です。 よって、close pathでは丸めとそれに伴う加算はクリティカルパスにはならず、シフトを行う部分がクリティカルパスになります。 なお、その気になれば+1加算は前倒しすることもできる(round_away_s ? suber_result[23:1] : suber_result[23:1] + 1のようにできる)ので、leading zeros countと正規化シフトを行っている裏で並列に計算する時間が十分にあります。 よって丸めとそれに伴う加算がクリティカルパスとなって困ることは決してありません。

終結果の選択

3つの計算結果から選択する部分がクリティカルパスです。 XilinxFPGAなら、クリティカルパスはLUT 1段分です。

場合分け一覧

テストベクトルを作るときの参考になるかと思って、代表的な場合分けを列挙してみました。 このほかにも丸めに関するコーナーケースのテストが必要です。

  • NaNを返す
    • 入力のどちらか一方がNaN
    • 実質的な減算で両オペランドが無限大
  • 無限大を返す
    • 正規化数を含む実質的な加算で繰り上がり指数部が255になる(実質的な指数部が254である数を含む場合のみ発生)
    • 片方のオペランドが無限大の実質的な加算(無限大を0x1.p+128とみなして上の経路で無限大を返すようにも設計できるが、(+Inf)+(+Inf)で指数部がオーバーフローするためやめたほうがよい。下の経路と共有するべき)
    • 片方のオペランドが無限大の実質的な減算(常に無限大を返すべきである。多くの場合には無限大を0x1.p+128とみなして計算しても同じ結果になるので、チェック漏れが発生しがち)
  • Far pathの結果を返す
    • 正規化数を含む実質的な加算で繰り上がらない(正規化数になる)
    • 正規化数を含む実質的な加算で繰り上がるが指数部が254以下になる(正規化数になる)
    • 非正規化数同士の実質的な加算で繰り上がる(正規化数になる。なお、ここでの「繰り上がる」は上の「繰り上がる」とは1 bitずれていることに注意)
    • 非正規化数同士の実質的な加算で繰り上がらない(非正規化数になる。(±0.0)+(±0.0)(複合同順)は(本当は非正規化数ではないけれど)この経路で±0.0(複合同順)が返ってくる)
    • “遠い”実質的な減算(必ず正規化数になる)
  • 符号が正のゼロを返す
    • オペランドの絶対値が一致している実質的な減算(“近い”実質的な減算に該当する。符号が正のゼロを返すのが正しい。(±0.0)+(∓0.0)(複合同順)もこれに該当する)
  • Close pathの結果を返す
    • 実質的な指数部が1である数同士の実質的な減算(“近い”実質的な減算に該当する。実質的に固定小数点数になっている領域であり、シフトと丸めの両方が不要となる。hidden bit相当のビット位置が0か1かで非正規化数か正規化数かがわかる。下の繰り下がらない経路と共有することも可能)
    • それ以外の“近い”実質的な減算で繰り下がらない(正規化数になる。シフトが不要となる)
    • “近い”実質的な減算で繰り下がりが発生し、繰り下がった桁数が左オペランドの実質的な指数部未満(正規化数になる。lzcに応じたシフトが必要。丸めは不要)
    • “近い”実質的な減算で繰り下がりが発生し、繰り下がった桁数が左オペランドの実質的な指数部以上(非正規化数になる。s_lhs_offsに応じたシフトが必要。丸めは不要)

まとめ

(答えを示されれれば)特に難しいところはなかったと思いますが、最近接丸めを行う浮動小数点数加算器の作り方でした。