読者です 読者をやめる 読者になる 読者になる

WGC1 第4章 浮動小数点表現 4.7浮動小数点演算

C

4.7 浮動小数点演算

浮動少数点の加算と減算はC/C++などの高級言語
浮動小数点の乗算と除算は高級言語よりもアセンブラ言語の方が簡単

浮動小数点の加算と減算のコード

#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つ注意すべきなのは、両者の符合が異なる場合は、
     * 実際は一方の値から他方の値を減算する必要があることである
     * (浮動少数点形式は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

これだけのMakefileでコンパイルできたんだ。。。知らんかった。 (..;

実行結果
[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

/Tp ファイルを .cpp としてコンパイルする オプション。
VC++C言語でinline使えないんだっけ??

実行結果
[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宛 で、お願いしまます。











Randall Hyde、鵜飼 文敏、まつもと ゆきひろ、後藤 正徳、トップスタジオ






広告を非表示にする