2012年03月23日

浮動小数点のカラクリ4

やっと晴れ春らしくなってきつつある今日この頃、いかがお過ごしでしょうか?コンドウです。

最近は家で暖房器具に電源入れる機会もぐっと減りエコひらめきしていますが、
まだまだ片づけるには早いですねあせあせ(飛び散る汗)

さて第4回となりました浮動小数点シリーズですが、今回は乗算について解説させていただきます。

最初から読まれていない方はまず、こちらからご覧いただきますとより理解が深まります。

浮動小数点のカラクリ
浮動小数点のカラクリ2
浮動小数点のカラクリ3


以前にも書いたかと思うのですが、実は乗算って加減算よりもずっと簡単です。

例えば、1234 × 528 について考えてみます。

まずは浮動小数点と同じように指数部と仮数部に分離します。
(分かりやすくする為に10進数表記で説明します)


  1234 = 1.234 × 103
  528 = 5.28 × 102

  1234 × 528 = (1.234×5.28) × 10(3+2) = 6.51552 × 105 = 651552


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


@ 指数部同士を足す
A 仮数部同士を掛ける
B Aの結果を正規化し指数部に反映
C 符号の処理


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


  1234 = (1+0.205078125) x 210 = 指数部10(+127)、仮数部1.205078125
  528 = (1+0.03125) x 29 = 指数部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 (結果からバイアス分引く)


でも同じ結果になりますのでこちらの計算を採用します。


Aは仮数部を取り出した状態の時には+1が省略されているので加算の時と同じように+1します。


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


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


  0x9A4000 × 0x840000 = 0x4F8900000000


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


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



Bですが、結果は23bitで収まっているので正規化の必要はありません。23bitをオーバー
している場合は適宜指数部への反映を行います。

最後にCの符号処理ですが乗算は同じ符号であればプラスに、異なる符号ならマイナスに
なりますので排他的論理和とれば簡単に再現できます。上記の場合はどちらもプラスなので
結果もプラスになります。
結果は


  符号 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){ // 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;
// どちらかが0なら計算せず0を返す
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;
}


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






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

次回は四則演算最後の除算の予定です。ではexclamation

posted by 管理人 at 14:13 | 研究・開発