株式会社ヘキサドライブ | HEXADRIVE

やっと春らしくなってきつつある今日この頃、いかがお過ごしでしょうか？コンドウです。
最近は家で暖房器具に電源入れる機会もぐっと減りエコしていますが、
まだまだ片づけるには早いですね
さて第４回となりました浮動小数点シリーズですが、今回は乗算について解説させていただきます。
最初から読まれていない方はまず、こちらからご覧いただきますとより理解が深まります。
浮動小数点のカラクリ
浮動小数点のカラクリ２
浮動小数点のカラクリ３
以前にも書いたかと思うのですが、実は乗算って加減算よりもずっと簡単です。
例えば、1234 × 528 について考えてみます。
まずは浮動小数点と同じように指数部と仮数部に分離します。
(分かりやすくする為に10進数表記で説明します)

　　1234 = 1.234 × 10³
　　528 = 5.28 × 10²
　　1234 × 528 = (1.234×5.28) × 10⁽³⁺²⁾ = 6.51552 × 10⁵ = 651552

という感じで計算します。ここで仮数部と指数部に分離しましたが、加減算のように指数部合わせの
フローが不要となり単純に処理できます。以上をまとめますと、

① 指数部同士を足す
② 仮数部同士を掛ける
③ ②の結果を正規化し指数部に反映
④ 符号の処理

順番に説明します。上と同じ、1234 × 528 を使って計算します。

　　1234 = (1+0.205078125) x 2¹⁰ = 指数部10(+127)、仮数部1.205078125
　　528 = (1+0.03125) x 2⁹ = 指数部9(+127)、仮数部1.03125

①は+127のバイアスが入っているので実際には137(10+127)と136(9+127)になります。
指数部同士を単純に足すと137+136=273ととても大きくなるので、本来なら

　　137-127=10 (バイアス引く)
　　136-127=9 (バイアス引く)
　　10+9=19 (指数部の足し算)
　　19+127=146 (結果にバイアスかける)

という手続きになりますが、

　　137+136=273 (バイアス引く前の段階で足す)
　　273-127=146 (結果からバイアス分引く)

でも同じ結果になりますのでこちらの計算を採用します。
②は仮数部を取り出した状態の時には+1が省略されているので加算の時と同じように+1します。

1234の仮数部 1.205078125 = 9A4000(16進数) = 100110100100000000000000(２進数)
528の仮数部 1.03125 = 840000(16進数) = 100001000000000000000000(２進数)

２進数で表記すると桁が多くなるので16進数で計算します。

　　0x9A4000 × 0x840000 = 0x4F8900000000

で、23bit目に小数点があるので小数点位置を合わせる為に結果を23bit右シフトします。
この時に前回説明しましたGRSビットへの反映も内部では忘れず行います。

　　0x4F8900000000 >> 23 = 0x9F1200 = 100111110001001000000000(２進数)

③ですが、結果は23bitで収まっているので正規化の必要はありません。23bitをオーバー
している場合は適宜指数部への反映を行います。
最後に④の符号処理ですが乗算は同じ符号であればプラスに、異なる符号ならマイナスに
なりますので排他的論理和とれば簡単に再現できます。上記の場合はどちらもプラスなので
結果もプラスになります。
結果は

　　符号 0
　　指数 146
　　仮数 0x1F1200(小数点では0.24273681640625)

となり、(1+0.24273681640625) x 2^(146-127) = 651552
と晴れて答えが求まりました。
いつものようにプログラム作ってみました。前回作った加算でGRSビットからの切り上げが
発生した時に正規化する処理を忘れていましたので今回は正規化する処理を関数化しました。

#define FRAC_LEN        23
#define GRS_LEN         3
#define BIAS            127
typedef unsigned __int64    uint64;
typedef signed __int64      sint64;
union FLOATING_VALUE{
float   FloatValue;
struct {
unsigned int Fraction:FRAC_LEN;
unsigned int Exp:31-FRAC_LEN;
unsigned int Sign:1;
}s;
};
// 正規化
void normalize( int *exp, unsigned int *frac )
{
unsigned int value;
int Exp = *exp, Frac = *frac;
// 1.0～1.99999…に収まるよう正規化
if(Frac == 0){                      // ０？
Frac = 0;
Exp = 0;
}else{
if( Frac >= (1<<(FRAC_LEN+GRS_LEN+1))){ // 桁あふれ (2.0以上)
do{
value = Frac & 1;   // シフトによって捨てられる値
Frac >>= 1;
Frac |= value;  // スティッキーbitの補充
Exp++;
}while(Frac >= (1<<(FRAC_LEN+GRS_LEN+1)) && Exp < 255);
}else if(Frac < (1<<(FRAC_LEN+GRS_LEN))){   // 桁借り(1.0未満)
do{
value = Frac & 1;
Frac = (Frac - value) << 1;
Frac |= value;  // スティッキーbitの補充
Exp--;
}while(Frac < (1<<(FRAC_LEN+GRS_LEN)) && Exp > 0);
}
if(Exp>=255){   // オーバーフロー？
// 無限大とする
Exp=255;
Frac=0;
}else if(Exp<0){    // アンダーフロー？
// 非正規値とする
Exp=0;
}
}
*exp = Exp;
*frac = Frac;
}
float FMul(float a, float b)
{
FLOATING_VALUE va,vb,vc;
int value;
va.FloatValue = a;
vb.FloatValue = b;
// どちらかが０なら計算せず０を返す
va.s.Fraction==0 && va.s.Exp==0) || (vb.s.Fraction==0 && vb.s.Exp==0)){
return 0.0f;
}
unsigned int kasuuA,kasuuB;
kasuuA = va.s.Fraction + (1<<FRAC_LEN); // 省略されている1.0(1<<23)を加える
kasuuB = vb.s.Fraction + (1<<FRAC_LEN); // 省略されている1.0(1<<23)を加える
int exp = va.s.Exp + vb.s.Exp - BIAS;                   // 指数部の計算
// 仮数部の計算(32bitで計算するとオーバーフローするので64bit型で計算)
uint64 kasuuL = static_cast<uint64>(kasuuA) * static_cast<uint64>(kasuuB);
// 右シフトによって捨てられるbit分
unsigned int stub = kasuuL & (( 1 << (FRAC_LEN - GRS_LEN)) - 1 );
kasuuL >>= (FRAC_LEN - GRS_LEN);        // 小数点位置に合わせる(GRSビットは残しておく)
kasuuA = static_cast<unsigned int>(kasuuL | (stub != 0));   // スティッキーbitの処理
normalize( &exp, &kasuuA );                             // 正規化
value = kasuuA & ((1<<GRS_LEN)-1);      // GRSビット抽出
if(value >= 4){                         // 0b100以上なら切り上げ、0b011以下なら切り下げ
kasuuA += 1<<GRS_LEN;               // +1切り上げ
normalize( &exp, &kasuuA );         // 再度正規化
}
kasuuA >>= GRS_LEN;                     // GRSビット捨てる
// 結果を格納
vc.s.Fraction = kasuuA;
vc.s.Exp = exp;
vc.s.Sign = va.s.Sign ^ vb.s.Sign;      // 符号の処理
return vc.FloatValue;
}

どちらかが０の場合は最初の方で０をリターンする処理が入っています。これは０を０として処理しなかった
場合は1.0×2^-127 の限りなく０に近いけど０で無い値として処理されてしまうからです。

---

CPUに付いている浮動小数点コプロセッサの中には、積和計算(乗算＋加算)や積差計算(乗算＋減算)を
１つの命令で行うものがあります。しかも大抵は加算や乗算と同じ所要クロックで行える為、ちょっと
お得な命令です。(乗算の結果に対して加算を行うので同じクロックで終えるのは謎です)
この積和命令って、内部では多bit状態(GRSビットなどを保持)のまま、計算し続けているので乗算と
加算を別々で行った場合に比べて計算精度が高かったりするコプロセッサもあります。２つの計算を
同時にしてくれてかつ精度が高いとはとてもお得感がありますね。但しIEEE754の仕様と異なる
(IEEE754には積和計算など無い)のでCのコンパイルオプションでは選択できるようになっているものも
あります。(浮動小数点の高速計算ON/OFFとか)
コプロセッサには便利な命令(例えば逆数を計算したり、平方根計算等)が積んでいたりするので使用の
コンパイラで組み込み関数が無いかもしくはアセンブラで直接書くすると色々と発見できます。
次回は四則演算最後の除算の予定です。では