本記事の目的は Berkeley SoftFloat による浮動小数点数の足し算の実装 (C言語) を理解することです.
- https://github.com/ucb-bar/berkeley-softfloat-3/blob/master/source/f64_add.c
- https://github.com/ucb-bar/berkeley-softfloat-3/blob/master/source/s_addMagsF64.c
話を簡単にするために倍精度 (f64) に限定して議論します. オリジナルのコードは if のネストがわかりにくい上に goto まで登場して, これぞスパゲッティコードという感じです. 私が書き直した Rust コードもある方がわかりやすいと思うので, 本記事に関する関数のみ gist に切り出しておきました.
なお Rust プロジェクト全体はまだ作成中でプライベートです. なのでこの gist では丸めや例外処理等が未実装です. その点に目をつむれば一応動くはずです (丸めによる誤差はあります).
準備
まずは f64 の仕様の確認です. これは 1bit の符号部, 11bit の指数部, 52bit (実質 53bit) の仮数部からなります. 従って, 次の関数により f64 (を f64::to_bits により u64 に変換したもの) を符号ビット, 指数部 i16, 仮数部 u64 に分解することができます.
fn sign_exp_sig(x: u64) -> (u64, i16, u64) {
let sign = u & 0x800_0000000000000; // アンダースコアは (符号部, 指数部) と仮数部の境界
let exp = ((u>>52) & 0x7ff) as i16;
let sig = u & Self::SIGNIFICAND;
(sign, exp, sig)
}
符号 sign と仮数部 sig (significand の略) が紛らわしいですが, 符号 sign はほとんど登場しないので許してください. 細かいことですが, 符号部はそのビットが最下位に来るようにシフトすれば 0 または 1 になってすっきりしますが, 最上位ビット以外の場所の処理に直接関与することはありませんから, わざわざ動かさないことにします.
指数部 exp は 0 から 2047 の値を取り, 0 のとき非正規化数またはゼロ, 2047 のとき無限大または NaN, 1 から 2046 のとき正規化数と言うのでした. 一方の仮数部は 0 以上 $2^{52}$ 未満の整数値を取ります. 正規化数に関しては, もとの浮動小数点数 x は, 次のように書けます.
$$x = (-1)^\mathrm{sign} \times \frac{ 2^{52} + \mathrm{sig} }{ 2^{52} } \times 2^{\mathrm{exp} - 1023}$$
正の正規化数の足し算
まずは話を簡単にするために, 正の正規化数 a: f64, b: f64 の足し算を考えましょう. この二つの数を指数部と仮数部に分解します (符号部は 0 なので忘れていいですね).
let (a, b) = (a.to_bits(), b.to_bits());
let (_, exp_a, sig_a) = sig_exp_sig(a);
let (_, exp_b, sig_b) = sig_exp_sig(b);
a と b が同程度の大きさ (つまり指数部が同じ値) であればさくっと足せそうですが, 違う大きさ (指数部が違う値) であれば, 足し算ができるように桁を揃える必要があります. なので最初のステップは指数部の差を計算し, それがゼロかどうか判定することです.
let exp_diff: i16 = exp_a - exp_b;
if exp_diff == 0 {
甲(exp_a, sig_a, exp_b, sig_b)
} else {
乙(exp_a, sig_a, exp_b, sig_b)
}
甲, 乙 という処理は, その最初の部分で入力が正規化数であるかどうかを判定し, そうでなければ別の処理を呼ぶことになりますが, いまは正規化数ということがわかっているので, それは後回しにしましょう.
甲: 指数部が共通だった場合の処理
まずは指数部が共通だった場合を考えます. このとき, いまの仮定 (正の正規化数) のもとでは, 求める足し算
$$a + b = \frac{ 2^{52} + \mathrm{sig}(a) }{ 2^{52} } \times 2^{\mathrm{exp} - 1023} + \frac{ 2^{52} + \mathrm{sig}(b) }{ 2^{52} } \times 2^{\mathrm{exp} - 1023}$$
(共通の指数部を exp と書きました) は次のように計算できます.
$$a + b = \frac{ 2^{53} + \mathrm{sig}(a) + \mathrm{sig}(b) }{ 2^{52} } \times 2^{\mathrm{exp} - 1023}$$
従って, 求める値は指数部が exp, 仮数部が $2^{52} + \mathrm{sig}(a) + \mathrm{sig}(b)$ ということになります. と思いきや, 仮数部は $2^{52}$ 未満でなければならなかったのですから, これではダメです. 必ず指数部を動かして帳尻を合わせる必要があります. その方法を論じる前に, 指数部が異なる場合の処理について見てみましょう.
乙: 指数部が異なる場合の処理
指数部が異なるふたつの浮動小数点数の足し算は, 特に二つの数の絶対値に大きな隔たりがある場合, 情報落ちと呼ばれる現象が起きることはよく知られています (知らない方は「数値計算の常識」を読みましょう). しかしながら, いま私たちは 53bit の仮数部を表現するのに u64 を使用しているため, 11bit の余白があります. そこで, 演算の間だけ, 仮数部のビットを 9bit 増やして小数点以下 60 桁 (2進数で) の数として扱うことにしましょう.
9bit シフトしたということをわかりやすくするために, 「shift された仮数部」の意味でこれを ssig と書くことにします.
let (ssig_a, ssig_b) = (sig_a << 9, sig_b << 9);
さて, いま指数部が異なる (exp_a != exp_b) 状況を考えていますが, exp_a > exp_b と仮定して一般性を失いません (そうでなければ a と b の役割を入れ替えればよいのです). このとき, a と b はともに正規化数ですから
$$0 < \mathrm{exp} ( b ) < \mathrm{exp} ( a ) < 2047$$
という不等式が成立しています. いま, 小さい方の数 b は
$$b = \frac{ 2^{61} + \mathrm{ssig}(b) }{ 2^{61} } \times 2^{\mathrm{exp} (b) - 1023} =: \frac{ \beta }{ 2^{61} } \times 2^{\mathrm{exp}(b) - 1023}$$
という形に表示されています (beta = ssig_b + 0x2000000000000000 とおきました) が, 指数部を動かしてこれを exp_a に調整してみましょう. この処理は Berkeley SoftFloat では softfloat_shiftRightJam64 という関数に対応します. いま exp_diff = exp_a - exp_b は正の値ですから
$$b = \frac{ \beta }{ 2^{61} } \times 2^{- \mathrm{diff}} \times 2^{\mathrm{exp} (b) + \mathrm{diff} - 1023}$$
つまり
$$b = \frac{ \beta \rhd \mathrm{diff} }{ 2^{61} } \times 2^{\mathrm{exp}(a) - 1023}$$
ということになります( $\rhd$ は右シフト演算子です). ただし注意点があります.
- シフトしたときに消えてしまったビットが 0 でないならば, その分四捨五入 (2進法でも四捨五入って言うんでしょうか) するべきです.
-
exp_diffが 63 より大きいとき, $\beta \rhd \mathrm{diff}$ はゼロになってしまいます. 最初から $\beta = 0$ ならそれで正しいのですが, $\beta$ がゼロでないのにゼロになってしまうというのは, 完全な情報落ちが起こっているということです. BSF にならい, その場合は $\beta \rhd \mathrm{diff}$ をゼロより大きい最小の数, すなわち 1 に置き換えます.
以上のことを考慮すると, $\beta \rhd \mathrm{diff}$ と書いているものは正確には次の計算をしていることになります.
fn shift_right_jam64(beta: u64, diff: i16) -> u64 {
if diff < 63 {
beta >> diff | ((beta << (-diff & 63)) != 0) as u64
} else { (beta != 0) as u64 }
}
これで, ようやく $a + b$ を計算することができました.
$$a + b = \frac{ 2^{61} + \mathrm{ssig}(a) + (\beta \rhd \mathrm{diff}) }{ 2^{61} } \times 2^{\mathrm{exp}(a) - 1023}$$
この分子を以下では $\alpha$ と書きます. いまの状況では $0 \leq \mathrm{ssig}(a) < 2^{61}$, $0 \leq ( \beta \rhd \mathrm{diff} ) < 2^{61}$ が成立しています ($\beta$ の定義に $2^{61}$ が含まれることに注意してください) から, $\alpha$ は
$$2^{61} \leq \alpha < 1.5 \times 2^{62}$$
を満足します. ということは, 入力によって最上位ビットの位置が違うことになります. それでは不便ですから, 乙の処理を抜ける前に, 必要に応じて調整しておくこととしましょう.
if alpha < 0x4000000000000000 {
(exp - 1, alpha << 1))
} else {
(exp, alpha)
}
仮数部の帳尻合わせ
甲乙どちらの処理も完成したので, 最後の帳尻合わせです. この処理は Berkeley SoftFloat では次のファイルに実装されています.
この部分は甲, 乙どちらの分岐を通ってきても共通です. なお甲の場合は, 乙で採用していた 9bit の追加の有効桁数に対応するために 9bit シフトしておきます. その結果, いま私たちの計算結果は
$$a + b = \frac{ \alpha }{ 2^{61} } \times 2^{\mathrm{exp}(a) - 1023}$$
という形に表示できています. $\alpha$ は $2^{62} \leq \alpha < 2^{63}$ を満足します (乙ではそうなるように調整しました. 甲では調整の必要はなく, 常にそうなっていることが簡単に確認できます). 従って, 丸めとして切り捨てを採用するならば, $\alpha$ を 10bit シフトして最上位ビットを落としたものが求める仮数部になります.
同時に, 指数部を 1 増やします. ただしこれにより指数部が正規化数の上限 2046 を超えてしまう可能性があり, その場合 Overflow 処理をしたうえで無限大を返します.
ただし, これらの挙動は丸め方式に依存します. それについても書くつもりだったのですが, 既にこの記事は十分に長くなってしまったので, また別記事にしたいと思います.
非正規化数, 正の無限大, NaN
さて, 正規化数については足し算に成功しました. 非正規化数, 無限大, そして NaN についても足し算ができなければなりません.
非正規化数とゼロ
非正規化数とゼロ (以下まとめて非正規化数と呼びます) は, 正規化数ではケチ表現のために省略されていた最上位ビットの 1 が仮数部に含まれるという点を除けば正規化数と同じです. そのため, 正規化数と非正規化数の足し算においては, $\beta$ の定義が修正されるだけで, あとはまったく同じです.
一方, 非正規化数と非正規化数の足し算は, 仮数部 sig_a, sig_b を足したものを (符号を除き) 足し算の結果の f64 のビット表示と解釈することができます. これは, 足し算の結果が非正規化数ならばすぐにわかりますが, 正規化数になる場合はどうでしょう？ その場合, 仮数部どうしの足し算により繰り上がりが起きますが, その繰り上がった数字は指数部のビットを形成し, まさに正しく欲しかった浮動小数点数のバイナリ表示になります. うまくできていますね.
無限大
有限数と無限大の足し算は無限大になります. 正の無限大と正の無限大の足し算も無限大です.
NaN
入力の少なくとも一方が NaN である場合, 出力も NaN です. 詳細はそのうち例外機構に関する別記事に書くんじゃないかなと思います.
負の数の足し算
お疲れ様でした. これで, 正の数と正の数の足し算はできるようになりました. それでは次に負の数と負の数の足し算について考えましょう.
これはとても簡単です. というのも, $a$ と $b$ が負の数のとき, $a + b = -((-a) + (-b))$ によりこの計算は正の数の足し算 $(-a) + (-b)$ に帰着できるからです. なので, 符号ビットを無視して正の数同士の足し算の処理を呼び出し, 最後に出力の符号ビットを反転させれば出来上がりです.
正と負の数の足し算
最後に, これまで論じてこなかった正と負の数の足し算をフォローする必要があります. けれど, それってよく考えるまでもなく, 引き算ですよね. ですから, この場合の足し算の処理は引き算の処理を呼び出すことによって実現されます. 浮動小数点数の引き算アルゴリズムについてはまた別の記事で扱いたいと思います.
終わりに
以上, 足し算についてお話しました. 長々と書きましたが, 個々の要素自体はさほど難しくはないはずで, それで浮動小数点数の足し算という難しそうなことが実現できるのは素晴らしいですね. 最後までご清聴ありがとうございました.
ライセンス表記
本記事に含まれる Rust コードは修正 BSD ライセンスのもとで配布される Berkeley SoftFloat 3 のC言語コードを Rust に移植したものです.
ライセンス原文
License for Berkeley SoftFloat Release 3e
John R. Hauser
2018 January 20
The following applies to the whole of SoftFloat Release 3e as well as to
each source file individually.
Copyright 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018 The Regents of the
University of California. All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
Redistributions of source code must retain the above copyright notice,
this list of conditions, and the following disclaimer.
Redistributions in binary form must reproduce the above copyright
notice, this list of conditions, and the following disclaimer in the
documentation and/or other materials provided with the distribution.
Neither the name of the University nor the names of its contributors
may be used to endorse or promote products derived from this software
without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS "AS IS", AND ANY
EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, ARE
DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE FOR ANY
DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.