浮動小数点の限界を知るためにビットパターンで演算を確認する

この記事について

浮動小数点はほぼ全てのコンピューター上でサポートされており数値計算をする上で非常に便利ですが、演算精度の予測が困難であり、移植性に課題があります。

その課題を理解するため、浮動小数点で表現される数値について単精度浮動小数点を例にビットパターンレベルで確認します。

その上で、浮動小数点の演算で誤差が生じる原因である端数処理(丸め処理)について実験コードを用いて記載し、浮動小数点演算のイメージをつかみます。

浮動小数点の問題

• $x, y$が単精度浮動小数点の場合、$ x+y-y \neq x $となる場合がある
• 0.01を100回加算すると、結果は0.999999となる

普通の計算においては全く成立しないような結果が、浮動小数点の計算では発生します。これは浮動小数点が実際は実数ではなく、有限個の値を表現した数値表現であることが原因です。

浮動小数点の定義

単精度浮動小数点の定義を定めるIEEE754-2008では単精度浮動小数点(binary32 / float)を以下のように定義しています。

• 符号部 (sign : b31)
• 指数部 (exponent : b30 - b23)
• 仮数部 (fraction : b22 - b0)

実験

binary32が数値保持や演算において、どのように振る舞うか実験コードを用意して確認します。この例では、binary32のみ対象にしていますが、倍精度浮動小数点(binary64 / double)でも桁数が異なるだけで同じ振る舞いとなります。

実験コード

符号部、指数部、仮数部を指定して単精度浮動小数点(binary32)のビットパターンを生成して値を確認します。

// 単精度浮動小数点(binary32)の数値とビットパターン(16進数)を表示
//
// 小数点以下46桁、有効桁数7桁の指数表示と
// ビットパターンを16進数表示します
static void print_float(const char* name, float v){
    printf("%s = %.46f (%e)\n", name, v, v);
    printf("bit pattern : %#x\n", *(unsigned int*)&v);
}

// ビットパターンをfloat型の値に変換
static float convert_float_from(unsigned int ptn){
    float val;  
    val = *((float*)&ptn); // メモリ上のパターンをbinary32として読み直します
    return val;
}

// binary32の定義式を表示
static void print_float_debug(unsigned int _sign, int _expo, unsigned int _frac){

    printf("[sign=%d, expo=%d(%d), frac=%d] := ", _sign, _expo, _expo-127, _frac);
    if( 1 <= _expo && _expo <= 254 ) // 指数部が0x01～0x0FE
    {
        printf("(-1)^(%d) x (1.0 + %d x 2^(-23)) x 2^(%d)\n", _sign, _frac, _expo-127);
    }
    else if( _expo == 0 ) // 指数部が0x00
    {
        if( _frac != 0 )
        {
            printf("(-1)^(%d) x (%d x 2^(-23-126))\n", _sign, _frac);
        }
        else
        {
            // Zero
        }
    }
    else if( _expo == 255 ) // 指数部が0xFF
    {
        if( _frac == 0 )
        {
            // Inf
        }
        else
        {
            // NaN
        }
    }
}

// 符号、指数、仮数の数値を与えて浮動小数点のビットパターンを作成
static float make_float_from(unsigned int _sign, int _expo, unsigned int _frac){
    unsigned int val = 0;
    unsigned int s, e, f;

    print_float_debug(_sign, _expo, _frac);

    s = _sign & 0x1;
    e = (_expo + 127) & 0xFF;
    f = _frac & 0x7FFFFF;

    val |= s << 31;
    val |= e << 23;
    val |= f;

    return convert_float_from(val);
}

実験プログラム

実験コードを用いて、binary32で表現される色々な数値を確認します。

最大値

binary32の最大値は、符号=0(+)、指数=127(254)、仮数=8388607 (=$2^{23}-1$) で構成された値となり、約$2^{128}$ = 約$3.4028234 \times 10^{38}$となります。(正確には、$(1+(2^{23}-1)/2^{23}) \times 2^{127}$となります。)

この数値は、float.hでFLT_MAXとして定められています。これよりも大きな値は正の無限大として扱われます。

program

print_float("MAX", make_float_from(0,127,8388607) );

output

[sign=0, expo=254(127), frac=8388607] := (-1)^(0) x (1.0 + 8388607 x 2^(-23)) x 2^(127)
MAX = 340282346638528859811704183484516925440.0000000000000000000000000000000000000000000000 (3.402823e+38)
bit pattern : 0x7f7fffff

正の最小値

正規化数の場合

正規化された状態での正の最小値は$2^{-126}$ = 約$1.175494 \times 10^{−38}$の値まで表現できます。

program

print_float("Normalization MIN", make_float_from(0,-126,0) );

output

[sign=0, expo=1(-126), frac=0] := (-1)^(0) x (1.0 + 0 x 2^(-23)) x 2^(-126)
Normalization MIN = 0.0000000000000000000000000000000000000117549435 (1.175494e-38)
bit pattern : 0x800000

非正規化数の場合

正規化数の最小値以下の数値を表現する非正規化数では、指数部は0となり$2^{-126}$であることを表します。仮数部の値がそのまま正規化されていない2進数の小数を表すことになり、最小で$2^{-126} \times 2^{-23}$ = $2^{-149}$ = 約$1.4 \times 10^{-45}$の値まで表現できます。

この数値は、float.hでFLT_TRUE_MINとして定められています。これより小さな正の数値はゼロとして扱われます。

program

print_float("Denormalization MIN", make_float_from(0,-127,1) );

output

[sign=0, expo=0(-127), frac=1] := (-1)^(0) x (1 x 2^(-23-126))
Denormalization MIN = 0.0000000000000000000000000000000000000000000014 (1.401298e-45)
bit pattern : 0x1

非正規化数はより小さな値を扱えるメリットの半面、ハードウェアサポートを得られず非常に低速になる可能性があります。また、非正規化数を含んだ演算を行うと有効桁数が極端に低下してしまうデメリットもあります。(参考： FLP05-C 非正規化数を使用しない)。そのため、必要に応じて、コンパイラの指定等で非正規化数を無効化することができるようになっています。

数値間の間隔

binary32では正の数値として$10^{-45}～10^{38}$オーダーという広範囲の数値を表現できることがわかりました。ただし、この範囲の全ての数値を表現できるわけではありません。

binary32では、仮数部が23ビットあることから、各$2^{n}$～$2^{n+1}$の区間を$2^{23}$個ずつ分割した数値を表すことができます。$2^{n}～2^{n+1}$の区間において、仮数部が1つ増えるごとに、$2^{-23}$と$2^{n}$を乗算した数値分だけ値が増えていく形になります。

例えば、指数部が-1(126)の場合、$2^{-23} \times 2^{-1}$ = $2^{-24}$ = 約0.000000059(5.9e-8)が仮数部1つ分が表す間隔の値となります。

program

float f_1, f_2;
f_1 = make_float_from(0, -1, 0); // 0.5(0x3f000000)
f_2 = make_float_from(0, -1, 1); // 0.5(0x3f000000)から1間隔ずつずらした値
print_float( "f_1", f_1 );
print_float( "f_2", f_2 );
print_float( "sub", f_2 - f_1 ); // 引き算した結果は仮数部1つ分の差

output

[sign=0, expo=126(-1), frac=0] := (-1)^(0) x (1.0 + 0 x 2^(-23)) x 2^(-1)
[sign=0, expo=126(-1), frac=1] := (-1)^(0) x (1.0 + 1 x 2^(-23)) x 2^(-1)

f_1 = 0.5000000000000000000000000000000000000000000000 (5.000000e-01)
bit pattern : 0x3f000000

f_2 = 0.5000000596046447753906250000000000000000000000 (5.000001e-01)
bit pattern : 0x3f000001

sub = 0.0000000596046447753906250000000000000000000000 (5.960464e-08)
bit pattern : 0x33800000

指数部を1増やして0とした場合、仮数部が1つ増えるごとに、$2^{-23}$ = 約1.19e-7ずつ増えていく形になります。先程の例と比べて、間隔は2倍となります。(指数部が1増えるごとに間隔は2倍となります)。

float f_1, f_2;
f_1 = make_float_from(0, 0, 0);  // 1.0(0x3f800000)
f_2 = make_float_from(0, 0, 1);  // 1.0(0x3f800000)から1間隔ずつずらした値
print_float( "f_1", f_1 );
print_float( "f_2", f_2 );
print_float( "sub", f_2 - f_1 ); // 先ほどの結果より2倍となるはず

[sign=0, expo=127(0), frac=0] := (-1)^(0) x (1.0 + 0 x 2^(-23)) x 2^(0)
[sign=0, expo=127(0), frac=1] := (-1)^(0) x (1.0 + 1 x 2^(-23)) x 2^(0)

f_1 = 1.0000000000000000000000000000000000000000000000 (1.000000e+00)
bit pattern : 0x3f800000

f_2 = 1.0000001192092895507812500000000000000000000000 (1.000000e+00)
bit pattern : 0x3f800001

sub = 0.0000001192092895507812500000000000000000000000 (1.192093e-07)
bit pattern : 0x34000000

値が大きくなってくると間隔も無視できない程度に大きくなります。例えば、指数部が22の場合、間隔は$2^{-23} \times 2^{22}$ = 1/2 = 0.5となります。$2^{22}$ = 4194304以上の数字では、0.5よりも小さな間隔の数値を表現することができなくなります。

program

float f_1, f_2;
f_1 = make_float_from(0, 22, 0); // 2^22(=4194304)
f_2 = make_float_from(0, 22, 1); // 2^22から1間隔ずつずらした値
print_float( "f_1", f_1 );
print_float( "f_2", f_2 );
print_float( "sub", f_2 - f_1 ); // 0.5となるはず

output

[sign=0, expo=149(22), frac=0] := (-1)^(0) x (1.0 + 0 x 2^(-23)) x 2^(22)
[sign=0, expo=149(22), frac=1] := (-1)^(0) x (1.0 + 1 x 2^(-23)) x 2^(22)

f_1 = 4194304.0000000000000000000000000000000000000000000000 (4.194304e+06)
bit pattern : 0x4a800000

f_2 = 4194304.5000000000000000000000000000000000000000000000 (4.194304e+06)
bit pattern : 0x4a800001

sub = 0.5000000000000000000000000000000000000000000000 (5.000000e-01)
bit pattern : 0x3f000000

誤差の原因

丸め誤差

先の例で、binary32では、指数部が22以上($2^{22}$以上)の領域では、0.5より小さな間隔は表現できないことがわかりました。この領域での演算過程で0.5より絶対値が小さい数値の加減算が入ってきた場合、演算結果が一定のルールによって丸められます。これにより、例えば、binary32の演算では、以下のようなことが発生します。

program

float x, y;
int e = 22;                            // 指数部の値です
x = make_float_from(0, e, 0);          // 2^22
y = make_float_from(0, (-23+e)-1, 0);  // 0.25f (2^(-1)より小さな値です)
print_float( "x        ", x );
print_float( "x + y    ", x + y );     // 2^22 + 0.25f -> 表現できないため丸められます
print_float( "x + y - y", x + y - y ); // 本来なら、xと同じ値になるはずです

output

[sign=0, expo=149(22), frac=0] := (-1)^(0) x (1.0 + 0 x 2^(-23)) x 2^(22)
[sign=0, expo=125(-2), frac=0] := (-1)^(0) x (1.0 + 0 x 2^(-23)) x 2^(-2)

x         = 4194304.0000000000000000000000000000000000000000000000 (4.194304e+06)
bit pattern : 0x4a800000

x + y     = 4194304.0000000000000000000000000000000000000000000000 (4.194304e+06)
bit pattern : 0x4a800000

x + y - y = 4194303.7500000000000000000000000000000000000000000000 (4.194304e+06)
bit pattern : 0x4a7fffff

このプログラムは以下の事を示しています。

• $x$は数値を表現しただけなので、正しく保持されています
• $x + y$の演算過程において、$2^{22}+0.25$という値を表現することがbinary32ではできずに$2^{22}$に丸められてしまっています(この段階で丸め誤差が発生しています)
• その後の0.25の減算は$2^{22}-0.25$として計算されています
• 減算の演算結果は指数部が21である$2^{21}$の領域に入ります
• $2^{21}$の領域の範囲では隣り合う数値の間隔は0.25となります
• そのため、$2^{22}-0.25$は正しく保持されます。

指数部がnの時、binary32の隣り合う数値の差は$2^{-23} \times 2^{n} = 2^{-23+n}$です。これよりも小さな値の加減算は、計算した結果をbinary32で表現することができず、丸められた結果になります。

表現できない値

浮動小数点は有限個の値を表現した値となっています。そのため、解像度に限界があり、実数のほとんどの値が表現できないことになります。

例えば、0.01は$2^{-7}～2^{-6}$の間にある数値です。この領域の値は、$2^{-7}$から$2^{-30}$の倍数ごと離れた値が表現できることになります。0.01はこの倍数の値とならず、丸められてしまいます。

program

float f_1, f_2;
f_1 = convert_float_from(0x3c23d70a); // 0.01fの一歩手前
f_2 = convert_float_from(0x3c23d70b); // 0.01fの一歩後
print_float( "f_1", f_1 );
print_float( "f_2", f_2 );
print_float( "0.1", 0.01f );          // f_1かf_2に丸められる
print_float( "sub", f_2 - f_1 );      // 2^(-30) (約9.313226e-10となるはず)

output

f_1  = 0.0099999997764825820922851562500000000000000000 (1.000000e-02)
bit pattern : 0x3c23d70a

f_2  = 0.0100000007078051567077636718750000000000000000 (1.000000e-02)
bit pattern : 0x3c23d70b

0.01 = 0.0099999997764825820922851562500000000000000000 (1.000000e-02)
bit pattern : 0x3c23d70a

sub = 0.0000000009313225746154785156250000000000000000 (9.313226e-10)
bit pattern : 0x30800000

丸められて0.01より僅かに小さな値となっています。このため、このまま繰り返し加算すると、その誤差が蓄積してしまいます。

丸めモード

浮動小数点の丸めの方法には以下のモードがあります。

• 正の無限大方向へ丸める（切り上げ）
• 負の無限大方向へ丸める（切り下げ）
• 0の方向へ丸める（切り捨て）
• 最も近い値の方向へ丸める（ISO丸め：中間値は偶数方向へ丸める）

現在の丸めモードはfenv.hに定義されているfegetround関数で取得できます(※C99以降)。

print_round_mode.c

#include <fenv.h>

void print_round_mode(){
    int mode;
    mode = fegetround();
    switch (mode) {
        case FE_UPWARD:
            printf("切り上げ\n");
            break;
        case FE_DOWNWARD:
            printf("切り下げ\n");
            break;
        case FE_TOWARDZERO:
            printf("切り捨て\n");
            break;
        case FE_TONEAREST:
            printf("ISO丸めｆ\n");
            break;
        default:
            break;
    }
}

移植性の問題の原因

これまでの結果より、IEEE754形式に準拠している浮動小数点をサポートしている環境間においても、さまざまな環境依存の変数があることを記載しました。これらの要因により、移植した際に同じ結果が得られるプログラムとならならい可能性があります。浮動小数点のプログラムを移植する際には、以下の点に注意が必要となります。

• 非正規化数の有効化 / 無効化
• 丸めモード
• コンパイラの最適化設定
  • 計算過程において精度を拡張した状態で内部演算を実行する / しない

演算精度が必要な場合の対策

浮動小数点の演算過程での丸め処理を正確に把握して、誤差を正確に予測することは困難なものとなります。そのため、精度が重要な場合は、誤差の発生予測が比較的容易な固定小数点の体系を採用することを検討すべきとされています。

• ダイナミックレンジの自由度が重要ではなく、演算結果の数値範囲や桁数があらかじめわかっている場合は、小数点以下がなくなるように倍数オフセットをかけて整数として計算する
• BCD (2進化10進コード)形式を採用する