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

ゴールデンウィークも終わり、今年入社した新人はそろそろ五月病を
発症している頃ではないでしょうか。私といえば季節の変わり目で
風邪を引き、病気に対して抵抗力が落ち込んでいるコンドウです
さて第回となりました浮動小数点シリーズですが、今回は四則演算最後の
除算について解説させていただきます。
最初から読まれていない方はまず、こちらからご覧いただきますとより理解が深まります。
浮動小数点のカラクリ
浮動小数点のカラクリ２
浮動小数点のカラクリ３
浮動小数点のカラクリ４
今回は除算です。除算は精度さえ気にしなければ除数の逆数での乗算に置き換える事ができます。
つまり、Ａ÷Ｂ　⇒　Ａ×（１÷Ｂ）に置き換えて計算すればいい訳です。
(逆数はニュートン法を用いると除算使わないで計算できます。)
しかし精度が必要な場合は除算がどうしても必要になります。
それでは除算方法について考えてみます。
除算は乗算に反対の処理を行う感じで進めます。
具体的な手順は以下の通りです。

① 指数部同士を引く
② 仮数部同士を割る
③ ②の結果を正規化し指数部に反映
④ 符号の処理

順番に説明します。1540 ÷ 14 を使って計算します。

1540 = (1+0.50390625) x 2¹⁰ = 指数部10、仮数部1.50390625
14 = (1+0.75) x 2³ = 指数部3、仮数部1.75

①は+127のバイアスが入っているので実際には137(10+127)と130(3+127)になります。
指数部同士を単純に引くと137-130=7ととても小さくなるので、本来なら

　　137-127=10 (バイアス引く)
　　130-127=3 (バイアス引く)
　　10-3=7 (指数部の引き算)
　　7+127=134 (結果にバイアス足す)

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

　　137-130=7 (バイアス引く前の段階で引く)
　　7+127=134 (結果にバイアス分足す)

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

1540の仮数部 1.50390625 = C08000(16進数) = 110000001000000000000000(２進数)
14の仮数部 1.75 = E00000(16進数) = 111000000000000000000000(２進数)

このまま普通に除算してしまうと0か1にしかならなくなるので
本来は計算結果に対して23bit左シフト(小数点の位置が23bit目にあるので)するところを
計算前の被除数に対して23bit左シフトして精度を保ちます。
２進数で表記すると桁が多くなるので16進数で計算します。

(0xC08000<<23) / 0xE00000 = 0x604000000000 / 0xE00000 = 0x6E0000
③ですが、0x6E0000は0.859375で１未満となるので正規化の必要があります。
　　仮数部　0x6E0000 << 1 = 0xDC0000 　　指数部　134 - 1 = 133 最後に④の符号処理は乗算と同じく同じ符号であればプラスに、異なる符号ならマイナスになりますので 排他的論理和とれば簡単に再現できます。上記の場合はどちらもプラスなので結果もプラスになります。 結果は
　　符号 0
　　指数 133
　　仮数 0x5C0000(1.0が加算されているので引いて、0xDC0000 – 0x800000 = 0x5C0000)

となり、(1+0.71875) x 2^(133-127) = 110
と晴れて答えが求まりました。
いつものようにプログラム作ってみました。プログラムでは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 FDiv(float a, float b)
{
FLOATING_VALUE va,vb,vc;
int value;
va.FloatValue = a;
vb.FloatValue = b;
unsigned int kasuuA,kasuuB;
// aが０なら計算せず０を返す
if(va.s.Fraction==0 && va.s.Exp==0){
vc.s.Fraction = 0;
vc.s.Exp = 0;
vc.s.Sign = va.s.Sign ^ vb.s.Sign;      // 符号の処理
return vc.FloatValue;
}
// bが０なら計算せずエラーを返す(divied by 0)
if(vb.s.Fraction==0 && vb.s.Exp==0){
// 無限大を返す
vc.s.Fraction = 0;
vc.s.Exp = 255;
vc.s.Sign = va.s.Sign ^ vb.s.Sign;      // 符号の処理
return vc.FloatValue;
}
kasuuA = va.s.Fraction + (1<<FRAC_LEN); // 省略されている1.0(1<<23)を加える
kasuuB = vb.s.Fraction + (1<<FRAC_LEN); // 省略されている1.0(1<<23)を加える
kasuuA <<= GRS_LEN;                     // GRSビット分桁上げ
kasuuB <<= GRS_LEN;                     // GRSビット分桁上げ
int exp = va.s.Exp - vb.s.Exp + BIAS;   // 指数部の計算
// 仮数部の計算(GRSビット分底上げ)
uint64 kasuuL = (static_cast<uint64>(kasuuA)<<(FRAC_LEN+GRS_LEN)) / static_cast<uint64>(kasuuB);
kasuuA = static_cast<unsigned int>(kasuuL);
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;
}

除数が０の場合、計算不能になるので事前に初期段階ではじいています。仮数部同士の除算では０チェック
していませんが、両方とも24bit目が必ず1になるので０チェックは省いています。
被除数が０の場合も除数の値に関わらず必ず０になるので、事前に０を返しています。

---

除算は加算や乗算のようにパイプライン化できない演算(同じ計算回路を何度も再利用しないといけない為)
なので非常に時間のかかる処理になってしまいます。しかもパイプライン構造でないので１回の除算が
完了するまで次の除算が計算できず待たされてしまうなどデメリットが非常に大きい計算です。
その為、可能であれば冒頭でも書きました逆数の乗算に置き換える等工夫する必要があります。
除算がどの位遅いかというとプロセッサにもよりますが乗算に比べると数倍～十数倍(実計算時間)の差があります。
例えば透視変換で、XとYをZで割る処理の場合
　　SX = X / Z;
　　SY = Y / Z;
となりますが、これを
　　INVZ = 1 / Z;
　　SX = X * INVZ;
　　SY = Y * INVZ;
に置き換えてもまだこっちの方が早かったりします。２個以上の除算は検討する必要がありますね。
最適化かけたCコンパイラなら定数の除数の場合は自動的に逆数の乗算などに置き換えてくれたりします。
しかし除数が変数の場合はそのまま除算命令になるので手動で最適化する必要があります。
さて５回に渡って浮動小数点の仕組みについて説明してきましたが、これらの計算部分はほとんど最適化
が施されていませんので実用的ではありません。あくまで浮動小数点の仕組みを理解する事が趣旨であり
分かりやすくした結果、速度が犠牲になっています。もしちゃんと実用的なロジックでシステムに組み込
みたいって思っている(よっぽど古いプロセッサ以外は必要ありませんが)方は是非gcc(gnuのcコンパイラ)
のソースを入手して解析してみて下さい。このシリーズでは説明していない浮動小数点⇒10進数表示など
も調べる事ができるかと思います。製作されるならビット長128bitなど超精度浮動小数点など面白いかも
しれませんね。
今回の浮動小数点ネタは私自身も深く理解する事が出来(資料作るにあたって色々調べました)とても楽しく
進める事ができました。またこのようなネタがありましたら私自身の勉強も兼ねて公表したいと思います。
次回も何かプログラムネタを探したいと思います。では