WGC1 第4章 浮動小数点表現 4.7浮動小数点演算
浮動小数点の加算と減算のコード
#include <stdio.h> #include <assert.h> typedef long unsigned real; #define asreal(x) (* ((float *)&x)) void fpadd(real left, real right, real *dest); void fpsub(real left, real right, real *dest); void fpsub(real left, real right, real *dest) { /* 右オペランドの符号ビットを反転させる */ right = right ^ 0x80000000; fpadd(left, right, dest); } inline int extractSign(real from) { return(from >>31); } inline int extractExponent(real from){ return((from>>23) & 0xff) - 127; } inline int extractMantissa(real from) { if ((from &0x7fffffff) == 0 ) return 0; return ((from &0x7FFFFF) | 0x800000); } void shiftAndRound(unsigned long *valToShift, int bitsToShift) { /* マスクを使ってビットをマスクアウトし、スティキービットを調べる */ static unsigned masks[24] = { 0, 1, 3, 7, 0xf, 0x1f, 0x3f, 0x7f, 0xff, 0x1ff, 0x3ff, 0x7ff, 0xfff, 0x1fff, 0x3fff, 0x7fff, 0xffff, 0x1ffff, 0x3ffff, 0x7ffff, 0xfffff, 0x1fffff, 0x3fffff, 0x7fffff }; /* HOmasks:masksエントリによってマスクされた値のHOビットをマスクアウトする */ static unsigned HOmasks[24] = { 0, 1, 2, 4, 0x8, 0x10, 0x20, 0x40, 0x80, 0x100, 0x200, 0x400, 0x800, 0x1000, 0x2000, 0x4000, 0x8000, 0x10000, 0x20000, 0x40000, 0x80000, 0x100000, 0x200000, 0x400000 }; /* shiftedOut:非正規化操作時に仮数からシフトアウトされる * 値を保持する (非正規化された値を丸めるために使用) */ int shiftedOut; assert(bitsToShift <= 23); /* 最初に、シフトアウトするビットを取得する（それにより、 * シフト後にこの値を丸める方法を判断できるようにする) */ shiftedOut = *valToShift & masks[bitsToShift]; /* 指定されたビット数だけ値を右シフトする * 注:ビッ31は常に0なので、Cコンパイラが * 右に論理シフトを行うか算術シフトを行うかは問題ではない */ *valToShift = *valToShift >> bitsToShift; /* 必要に応じて、値を丸める*/ if (shiftedOut > HOmasks[bitsToShift]) { /* シフトアウトしたビットがLOビットの半分より大きい場合は、 * 1を加えて値を大きい値に丸める */ *valToShift = *valToShift +1; } else if (shiftedOut == HOmasks[bitsToShift]) { /* シフトアウトしたビットがLOビットの値のちょうど半分の場合は、 * LOビットが0になる最近値に丸める */ *valToShift = *valToShift + (*valToShift & 1); } else { /* それ以外の場合は値の前の前の値に丸める。 * 現在の値はすでに切り捨てられて * (小さい値に丸められて)いるので * 何もする必要はない */ } } /* 符号、指数、仮数の各フィールドを32ビットの * 「実数」値にパックする。正規化値、非正規化値、および0では * 有効に機能するが、NaNおよび無限大では機能しない */ inline real packFP(int sign, int exponent, int mantissa) { return (real) ( (sign <<31) |((exponent + 127) << 23) |(mantissa & 0x7fffff) ); } /* 次の計算を行う: * dest = left + right * 3つのオペランドはすべて「実数」値(32ビット浮動小数点数値) */ void fpadd(real left, real right, real *dest) { /* 次の各変数には、左オペランドに関連するフィールドが格納される: */ int Lexponent; long unsigned Lmantissa; int Lsign; /* 次の各変数には、右オペランドに関連するフィールドが格納される: */ int Rexponent; long unsigned Rmantissa; int Rsign; /* 次の各変数には、結果の個別のフィールドが格納される: */ int Dexponent; long unsigned Dmantissa; int Dsign; /* 操作しやすいようにフィールドを抽出する: */ Lexponent = extractExponent(left); Lmantissa = extractMantissa(left); Lsign = extractSign(left); Rexponent = extractExponent(right); Rmantissa = extractMantissa(right); Rsign = extractSign(right); /* 特殊なオペランド(無限大とNaN)を処理するコード */ if ( Lexponent == 128 ) { if ( Lmantissa == 0x800000 ) { /* 右オペランドが無限大の場合、結果は * 右オペランドの値によって異なる */ if ( Rexponent == 128 ) { /* 指数のビットがすべて1(非バイアス後に128)の場合は * この値が無限大か(仮数部が0x800000)、QNaN(仮数=0xc00000)か、 * SNaN(仮数部が0以外で0x800000と等しくない)かが * 仮数によって決まる */ if ( Rmantissa == 0x800000 ) { /* 値は無限大か? */ /* 無限大 + 無限大 = 無限大 * -無限大 - 無限大 = -無限大 * -無限大 + 無限大 = NaN * 無限大 - 無限大 = NaN */ if ( Lsign == Rsign ) { *dest = right; } else { *dest = 0x7fC00000; /* +QNaN */ } } else { /* Rmantissaは0x800000以外なので、NaNである */ /* 右はNaNなので、それを伝播する */ *dest = right; } } else { *dest = left; } } else { /* Lmantissaは0以外で、Lexponentはすべて1 */ /* 左オペランドがNaNの場合は、結果も同じNaN */ *dest = left; } /* 結果はすでに計算したので、このまま戻る */ return; } else if (Rexponent == 128 ) { /* 2つの場合がある:右がNaNの場合(この場合は * 左の値にかかわらずNaNを伝播する必要がある)と * 左は「通常の」数値なので * * 無限大も伝播することになる。通常の数値に * 無限大を加えると無限大になるからである。 */ /* 右はNaNなので、それを伝播する */ *dest = right; return; } /* これで2つの浮動少数点数値の用意ができたので、値を加算する * 最初に、2つのオペランドの指数が異なる場合は、一方の * オペランドを「非正規化」する必要がある(値を加算または減算する * ときは指数が同じでなければならない) * * アルゴリズム:指数が小さいほうの値を選択する。 * その値の仮数を、2つの指数の差で指定されるビット数だけ * 右にシフトする */ Dexponent = Rexponent; if ( Rexponent > Lexponent ) { shiftAndRound(&Lmantissa, (Rexponent - Lexponent)); } else if ( Rexponent < Lexponent) { shiftAndRound(&Rmantissa, (Lexponent - Rexponent)); Dexponent = Lexponent; } /* 仮数を加算する。１つ注意すべきなのは、両者の符合が異なる場合は、 * 実際は一方の値から他方の値を減算する必要があることである * (浮動少数点形式は1の補数なので、小さいほうの仮数から大きいほうの * 仮数を減算し、元の符号の値と大きい方の仮数の組み合わせに * 応じて結果の符合をセットする */ if ( Rsign ^ Lsign ) { /* 符号が異なる場合は、一方の値から他方の値を減算する必要がある */ if ( Lmantissa > Rmantissa ) { /* 左の値の方が大きいので、結果は左オペランドの * 符号を継承する */ Dmantissa = Lmantissa - Rmantissa; Dsign = Lsign; } else { /* 右の値のほうが大きいので、結果は右オペランドの * 符号を継承する */ Dmantissa = Rmantissa - Lmantissa; Dsign = Rsign; } } else { /* 符号は同じなので、値を加算する */ Dsign = Lsign; Dmantissa = Lmantissa + Rmantissa; } /*ここで結果を正規化する * * 加算/減算時に、1ビットのオーバーフローが起きる可能性がある * ここでそれに対処する(オーバーフローが発生した場合は、 * 仮数を右の位置に1つシフトし、その分だけ指数をインクリメント * して調整する)。指数をインクリメントしたときにオーバーフローが * 発生した場合、このコードは無限大を返すことに注意 * (無限大は指数の値が$FFになった状態) */ if ( Dmantissa >= 0x1000000 ) { /* 加算/減算を行った時に余分なビットが2つ以上出ることはない * ここで使用している浮動小数点形式において、 * 加算または減算によって返される仮数の値は * 0x1fffffeが上限である。したがって、この値を丸めるときに * 25番目のビットにオーバーフローすることはない */ /* 結果を24ビットに移動する * シフト操作で2による除算が行われた * これでシフトの影響を相殺する * (指数をインクリメントして値を2掛ける) */ shiftAndRound(&Dmantissa, 1); ++Dexponent; }else { /* HOビットがクリアの場合は、ビットを左側に * シフトし、同時に指数をデクリメントして * 結果を正規化する。ゼロは非常によくある * 結果なので特殊なケースとして扱う */ if ( Dmantissa != 0 ) { /* whileループ(左シフトによって)仮数に * 2を掛けた後、数値全体を(指数をデクリメントして) * 2で割る。DmantissaのHOビットがセットされるか、 * 指数が-127(バイアス127形式の0)になるまで * これが継続される。 * Dexponentが-128に達した場合は値が正規化されているので * ループを停止する */ while ( (Dmantissa < 0x800000) && ( Dexponent > -127 )) { Dmantissa = Dmantissa << 1; --Dexponent; } } else { /* 仮数が0になった場合は、ほかのビットもすべてクリアする */ Dsign = 0; Dexponent = -127; } } /* 結果を再構築して格納する */ *dest = packFP(Dsign , Dexponent, Dmantissa); } int main(int argc, char **argv) { real l, r, d; asreal(l) = 1.0; asreal(r) = 2.0; fpadd(l, r, &d); printf("dest = %x\n", d); printf("dest = %12E\n", asreal(d)); l = d; asreal(r) = 4.0; fpsub(l, r, &d); printf("dest2 = %x\n", d); printf("dest2 = %12E\n", asreal(d)); return 0; }
gcc
コンパイル
Floating-Point_Addition_and_Subtraction:
[oc@centos5 4.7.2_Floating-Point_Addition_and_Subtraction]$ make cc -g Floating-Point_Addition_and_Subtraction.c -o Floating-Point_Addition_and_Subtraction
実行結果
[oc@centos5 4.7.2_Floating-Point_Addition_and_Subtraction]$ ./Floating-Point_Addition_and_Subtraction dest = 40400000 dest = 3.000000E+00 dest2 = bf800000 dest2 = -1.000000E+00
途中、結果が間違ってたのに気づかなかった。。。orz
1.0+2.0=3.0
3.0-4.0=-1.0
0、fが足りなかったり、配列の要素が抜けてたり
写経だとこういうミスが置きやすいな。
デバッグしても、どう動くべきかがわからないと原因を見つけるのが大変だ。
まぁ、そのおかげでemacsでのデバッグになれたかな
Visucla C++
環境
Visual C++ 9 (Microsoft(R) 32-bit C/C++ Optimizing Compiler Version 15.00.21022.08 for 80x86)
Windows XP
コンパイル
[D:\workspace\C_SandBox\10BOOKs\WGC1\4.7.2_Floating-Point_Addition_and_Subtraction]cl /Tp Floating-Point_Addition_and_Subtraction.c Microsoft(R) 32-bit C/C++ Optimizing Compiler Version 15.00.21022.08 for 80x86 Copyright (C) Microsoft Corporation. All rights reserved. Floating-Point_Addition_and_Subtraction.c Floating-Point_Addition_and_Subtraction.c : warning C4819: ファイルは、現在のコード ページ (932) で表示できない文字を含んでいます。データの損失を防ぐために、ファイルを Unicode 形式で保存してください。 Microsoft (R) Incremental Linker Version 9.00.21022.08 Copyright (C) Microsoft Corporation. All rights reserved. /out:Floating-Point_Addition_and_Subtraction.exe Floating-Point_Addition_and_Subtraction.obj
実行結果
[D:\workspace\C_SandBox\10BOOKs\WGC1\4.7.2_Floating-Point_Addition_and_Subtraction]./Floating-Point_Addition_and_Subtraction.exe dest = 40400000 dest = 3.000000E+000 dest2 = bf800000 dest2 = -1.000000E+000
書籍だと intになってるが、コンパイルエラーになるので、longに変えた。
void shiftAndRound(unsigned ing *valToShift, int bitsToShift) ↓ void shiftAndRound(unsigned long *valToShift, int bitsToShift) }}
[D:\workspace\C_SandBox\10BOOKs\WGC1\4.7.2_Floating-Point_Addition_and_Subtraction]cl /Tp Floating-Point_Addition_and_Subtraction.c
Microsoft(R) 32-bit C/C++ Optimizing Compiler Version 15.00.21022.08 for 80x86
Copyright (C) Microsoft Corporation. All rights reserved.
Floating-Point_Addition_and_Subtraction.c
Floating-Point_Addition_and_Subtraction.c(194) : error C2664: 'shiftAndRound' : 1 番目の引数を 'unsigned long *' から 'unsigned int *' に変換できません。(新しい機能 ; ヘルプを参照)
指示された型は関連がありません。変換には reinterpret_cast、C スタイル キャストまたは関数スタイルのキャストが必要です。
Floating-Point_Addition_and_Subtraction.c(196) : error C2664: 'shiftAndRound' : 1 番目の引数を 'unsigned long *' から 'unsigned int *' に変換できません。(新しい機能 ; ヘルプを参照)
指示された型は関連がありません。変換には reinterpret_cast、C スタイル キャストまたは関数スタイルのキャストが必要です。
Floating-Point_Addition_and_Subtraction.c(248) : error C2664: 'shiftAndRound' : 1 番目の引数を 'unsigned long *' から 'unsigned int *' に変換できません。(新しい機能 ; ヘルプを参照)
指示された型は関連がありません。変換には reinterpret_cast、C スタイル キャストまたは関数スタイルのキャストが必要です。
添削歓迎
ここ間違ってるよ
こうした方がよくないか？
こういうことなじゃないかな。
この環境だとこうなるよ
などなど
方法は、コメント、トラックバック、はてブ、Twitter @orange_clover宛 で、お願いしまます。
