2012年02月08日

浮動小数点のカラクリ3

しばらくの間、休暇を頂いていましたのでご無沙汰になります。コンドウです。

健康の為、普段は極力歩くよう心掛けているのですが
寒すぎてついつい車車(セダン)に乗ってしまう怠け者ですもうやだ〜(悲しい顔)

さて前回から浮動小数点関係のネタで攻めていますが、
ここまでいったら四則演算までいっちゃおうかと
決めましたのでこのネタ終わるまでお付き合いの程宜しくお願いしますグッド(上向き矢印)

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

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


実は前回、前々回共にさらーっと書いてた計算による桁落ち処理なんですが、詳しく調べたところ
もうちょっと複雑な仕組みになっているようですのでその辺を含めてまずは説明します。
これらははみ出た端数をどう扱うかの丸め処理にも大きく関係してきます。


・ガードビット、ラウンドビット、スティッキービット

IEEE754では計算の精度を高める為に仮数部の下位に3bitを追加して計算するようです。
それぞれのビットに名前があり、上位からガードビット、ラウンドビット、
スティッキービットと命名されています。

仮数部
1.00000000000000000000000GRS

"1."は隠れている+1.0です。
"00000000000000000000000"は仮数部本体です。
"G"はガードビット
"R"はラウンドビット
"S"はスティッキービット
全部で1+23+1+1+1=27bitになります。


ガードビット(Guard bit)・ラウンドビット(Round bit)

ガードビットとラウンドビットは仮数の最下位ビットから加算や減算などの計算の為に
こぼれてしまった捨てられてしまうビットを一時的に貯蔵する為のビットです。


スティッキービット(Sticky bit)

スティッキービットは先の2bitとはちょっと性質が異なり、ラウンドビットから右シフト
されて落ちた1を保持し続けるビットです。一度1が入ると0に戻ることはありません。
最初は0で初期化されます。

1.00000000000000000110110 を右シフトしていくと以下のようなイメージになります。

1.00000000000000000110110 0 0 0
0.10000000000000000011011 0 0 0
0.01000000000000000001101 1 0 0
0.00100000000000000000110 1 1 0
0.00010000000000000000011 0 1 1
0.00001000000000000000001 1 0 1
0.00000100000000000000000 1 1 1


今回は前回の加算処理にこの3ビットを加味し最近値への丸め(*1)で処理したプログラムに
改造してみたいと思います。
折角改造するなら前回さらーっと端折った負の数や0との加算や無限大などのエラー処理も
付け加えたいと思います。


#define FRAC_LEN 23 // 仮数部のbit数
#define GRS_LEN 3 // 付加される3bit

union FLOATING_VALUE{
float FloatValue;
struct {
unsigned int Fraction:FRAC_LEN;
unsigned int Exp:31-FRAC_LEN;
unsigned int Sign:1;
}s;
};

float FAdd(float a, float b)
{
FLOATING_VALUE va,vb,vc;
int kasuuA,kasuuB;
int shift,maxExp,value;

va.FloatValue = a;
vb.FloatValue = b;
vc.s.Sign = 0;
kasuuA = (va.s.Fraction + (1<<FRAC_LEN)) << GRS_LEN; // 省略されている1.0(1<<23)を加える
kasuuB = (vb.s.Fraction + (1<<FRAC_LEN)) << GRS_LEN; // 省略されている1.0(1<<23)を加える
if(va.s.Sign) kasuuA = -kasuuA; // 負の数なら仮数部を符号反転
if(vb.s.Sign) kasuuB = -kasuuB; // 負の数なら仮数部を符号反転
if(va.s.Exp > vb.s.Exp){ // aの方が大きい
maxExp = va.s.Exp;
shift = va.s.Exp - vb.s.Exp;
if(shift >= (1+FRAC_LEN+GRS_LEN)){ // 仮数部が無くなる程のシフトなら
kasuuB = (kasuuB != 0); // 1のbitが1つでもあればスティッキービットを1に
}else{
value = kasuuB & ((1<<shift)-1);// シフトによって捨てられる値
kasuuB >>= shift; // 仮数部の指数部合せ
kasuuB |= (value!=0); // スティッキーbitの補充
}
}else{ // bの方が大きい
maxExp = vb.s.Exp;
shift = vb.s.Exp - va.s.Exp;
if(shift >= (1+FRAC_LEN+GRS_LEN)){ // 仮数部が無くなる程のシフトなら
kasuuA = (kasuuB != 0); // 1のbitが1つでもあればスティッキービットを1に
}else{
value = kasuuA & ((1<<shift)-1);// シフトによって捨てられる値
kasuuA >>= shift; // 仮数部の指数部合せ
kasuuA |= (value!=0); // スティッキーbitの補充
}
}
kasuuA += kasuuB; // 仮数部の加算
if(kasuuA < 0){ // 負の数になっていたら
kasuuA = -kasuuA; // 正の数化
vc.s.Sign=1; // 符号ビット立てる
}

// 1.0〜1.99999…に収まるよう正規化
if(kasuuA == 0){ // 0?
vc.s.Sign=0;
kasuuA = 0;
maxExp = 0;
}else{
if( kasuuA >= (1<<(FRAC_LEN+GRS_LEN+1))){ // 桁あふれ (2.0以上)
do{
value = kasuuA & 1; // シフトによって捨てられる値
kasuuA >>= 1;
kasuuA |= value; // スティッキーbitの補充
maxExp++;
}while(kasuuA >= (1<<(FRAC_LEN+GRS_LEN+1)) && maxExp < 255);

}else if(kasuuA < (1<<(FRAC_LEN+GRS_LEN))){ // 桁借り(1.0未満)
do{
value = kasuuA & 1;
kasuuA = (kasuuA - value) << 1;
kasuuA |= value; // スティッキーbitの補充
maxExp--;
}while(kasuuA < (1<<(FRAC_LEN+GRS_LEN)) && maxExp > 0);
}
if(maxExp>=255){ // オーバーフロー?
// 無限大とする
maxExp=255;
kasuuA=0;
}else if(maxExp<0){ // アンダーフロー?
// 非正規値とする
maxExp=0;
}
}
value = kasuuA & ((1<<GRS_LEN)-1); // GRSビット抽出
kasuuA >>= GRS_LEN; // GRSビット捨てる
if(value >= 4){ // 0b100以上なら切り上げ、0b011以下なら切り下げ
kasuuA++;
}
// 結果を格納
vc.s.Fraction = kasuuA;
vc.s.Exp = maxExp;
return vc.FloatValue;
}

float FSub(float a, float b)
{
FLOATING_VALUE vb;

vb.FloatValue = b;
vb.s.Sign ^=1;
return FAdd(a,vb.FloatValue);
}


今回は減算(FSub)も用意しました。減算は加算処理の仮数部同士の計算のところを
+から-にするだけで対応できますが、上記のプログラムでは引く値を符号反転して
加算に置き換えています。
今回丸め処理をもうちょっとまじめにやってみましたが、まだWindowsの浮動小数点とは
最大1bitの誤差が出てしまうようです。(丸め処理のところかな?)
その辺も含めて自身への宿題にしたいと思います。
次回こそは乗算したいと思います。乗算計算でもGRSビットが活躍しますexclamation
乞うご期待下さいexclamation×2

*1 丸め処理の1つで正しい値に近い方へ丸める方法です。他にもゼロ方向への丸め、
正の無限大方向への丸め、負の無限大方向への丸めなどがあります。



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