浮動小数点数の誤差を考慮した比較【double/float型の正しい比較方法】

double型やfloat型の浮動小数点数同士の比較は、演算時に発生する誤差の影響によって、比較演算子による正確な比較が行なえなくなることがあります。

double a = 1, b = 1;

b += 0.001;
b -= 0.001; // 誤差が発生(b: 0.99999999999999988)

if (a == b) {
  printf("same"); // 実行されない
  // bに誤差が発生しているためaと一致しない
}

浮動小数点数の誤差を考慮した等価比較

次の方法で、浮動小数点数の誤差を考慮した等価比較を実現することができます。

// #include <float.h> // DBL_EPSILON
// #include <math.h>  // fabs, fmax
double a = 1, b = 1 + 0.001 - 0.001;

if (fabs(a - b) <= DBL_EPSILON * fmax(1, fmax(fabs(a), fabs(b)))) {
   printf("equal"); // 実行される
}

二つの値の差が誤差の範囲内に収まっているかどうかを判定しています。

より具体的には、両値の差の絶対値が極めて小さな値かどうかを判定しています。上記の例では、計算機イプシロン(DBL_EPSILON)という限りなくゼロに近い値を許容誤差として用いていますが、場合によっては0.000011e-10等のより精度の低い値を利用することも可能です。

上記の式に関するより具体的な解説や計算機イプシロン、fmaxの乗算に関する説明については、次項以降で行います。

double/float/long double型別の比較

double型やfloat型、long double型で、それぞれ異なる関数や定数を用いる必要があります。

浮動小数点数型絶対値関数計算機イプシロン
floatfabsfFLT_EPSILON
doublefabsDBL_EPSILON
long doublefabslLDBL_EPSILON
// #include <float.h> // FLT_EPSILON, DBL_EPSILON, LDBL_EPSILON
// #include <math.h>  // fabsf, fabs, fabsl, fmaxf, fmax, fmaxl

/* float型(倍精度浮動小数点数) */
if (fabsf(a - b) <= FLT_EPSILON * fmaxf(1.f, fmaxf(fabsf(a), fabsf(b)))) {}
/* double型(単精度浮動小数点数) */
if (fabs(a - b) <= DBL_EPSILON * fmax(1, fmax(fabs(a), fabs(b)))) {}
/* long double型(四倍精度浮動小数点数) */
if (fabsl(a - b) <= LDBL_EPSILON * fmaxl(1.L, fmaxl(fabsl(a), fabsl(b)))) {}

計算機イプシロンによる等価比較について

実数型の値に対して、丸め誤差や桁落ち、情報落ちを考慮した等価比較を行うためには、浮動小数点数の差が「計算機イプシロン(machine epsilon, マシンイプシロン、機械エプシロン)」という特殊な値を下回っているかどうかを判断する必要があります

以下のコードはdouble型に対する等価比較とその基本的な例です。

// #include <float.h> // DBL_EPSILON
// #include <math.h>  // fabs

double x = 0.1, y = 1 - 0.9; // {x: 0.10000000000000001, y: 0.099999999999999978}

if (fabs(x - y) < DBL_EPSILON) {
  printf("equal");
}

//    x - y    ≒ 0.00000000000000002
// DBL_EPSILON ≒ 0.00000000000000022

上記のfabs(x - y) < DBL_EPSILONが誤差を考慮した比較処理の例です。DBL_EPSILONが実際の計算機イプシロンを表す定数です。fabs関数は絶対値を求めるための関数です。DBL_EPSILONとfabsはそれぞれfloat.h, math.hヘッダーに定義されています。

#include <float.h>
#define DBL_EPSILON 2.2204460492503131E-16 // (0.00000000000000022204460492503131)

#include <math.h>
double fabs(double x); // 引数`x`の符号を取り除いた正の数を返す

実際の開発ではDBL_EPSILONによるイプシロンを誤差の厳格な許容範囲として捉え、fabs(x - y) <= DBL_EPSILONによる式を用いることが合理的です。また相対的な誤差を考慮する場合には、イプシロンに対して絶対値を乗算させる必要があります。

double x = 100 * (  0.1  ); // x == 10
double y = 100 * (1 - 0.9); // y == 9.9999999999999982

if (fabs(x - y) <= DBL_EPSILON * fmax(1, fmax(fabs(x), fabs(y)))) {
   printf("equal");
}

また複数回の演算によって、計算機イプシロンの範囲を超えた誤差が蓄積されるケースも考えられます。その場合は計算機イプシロンではなく、誤差の許容範囲を独自に決めて等価比較する必要があります。次項を参考にしてください。

許容範囲の独自指定

実際の開発では、扱う浮動小数点数の桁数に応じて、適した限界値を指定することが有効です。

たとえば、浮動小数点数型変数の値が0.11.0の範囲で増減するような処理では、限界値を0.010.001に指定すること等が考えられます。

double a = 1, b = 1 + 0.9 - 0.9 + 0.8 - 0.8 + 0.6 - 0.6;
assert( fabs(a - b) <= 0.01        ); // ok
assert( fabs(a - b) <= DBL_EPSILON ); // failed

許容誤差を取る比較関数の定義

以下は、浮動小数点数の誤差を考慮した比較関数の実装例です。

int dblcmpt(double x, double y, double tolerance) {
   return (x > y + tolerance) ? 1 : (y > x + tolerance) ? -1 : 0;
}

dblcmpt関数一つで「以上、以下、超過、未満、等価、不等」全ての判定を行うことができます。

assert( dblcmpt(1, 1.02, 0.01) <  0 ); // 未満
assert( dblcmpt(1, 1.01, 0.01) <= 0 ); // 以下
assert( dblcmpt(1, 0.99, 0.01) == 0 ); // 等価
assert( dblcmpt(1, 0.98, 0.01) != 0 ); // 不等
assert( dblcmpt(1, 0.99, 0.01) >= 0 ); // 以上
assert( dblcmpt(1, 0.98, 0.01) >  0 ); // 超過

double a = 100 * (  0.1  );
double b = 100 * (1 - 0.9);
assert( dblcmpt(a, b, DBL_EPSILON * fmax(1, fmax(fabs(a), fabs(b)))) == 0 );
広告