このトピックは StackOverflow で何度も取り上げられてきましたが、これは新しい見方だと思います。はい、私はBruce Dawson の記事とWhat Every Computer Scientist Should Know About Floating-Point Arithmeticとthis nice answerを読みました。
私が理解しているように、典型的なシステムでは、浮動小数点数が等しいかどうかを比較するときに 4 つの基本的な問題があります。
- 浮動小数点の計算は正確ではありません
- が「小さい」かどうかは、と
a-bのスケールに依存します。ab a-bが「小さい」かどうかは、aandの型b(例: float、double、long double)によって異なります。- 通常、浮動小数点には +-無限大、NaN、および非正規化表現があり、これらはいずれも単純な定式化を妨げる可能性があります
この答え- 別名。「Google アプローチ」 -- 人気があるようです。それはすべてのトリッキーなケースを処理します。また、比較を非常に正確にスケーリングし、2 つの値が互いに一定数のULP内にあるかどうかをチェックします。したがって、たとえば、非常に大きな数は、「ほぼ等しい」と無限大を比較します。
でも：
- 私の意見では、それは非常に面倒です。
- 内部表現に大きく依存し、ユニオンを使用して float からビットを読み取るなど、特に移植性はありません。
- 単精度と倍精度の IEEE 754 のみを処理します (特に、x86 long double は処理しません)。
似たようなものが欲しいのですが、標準の C++ を使用し、長い double を処理します。「標準」とは、可能であれば C++03、必要であれば C++11 を意味します。
これが私の試みです。
#include <cmath>
#include <limits>
#include <algorithm>
namespace {
// Local version of frexp() that handles infinities specially.
template<typename T>
T my_frexp(const T num, int *exp)
{
typedef std::numeric_limits<T> limits;
// Treat +-infinity as +-(2^max_exponent).
if (std::abs(num) > limits::max())
{
*exp = limits::max_exponent + 1;
return std::copysign(0.5, num);
}
else return std::frexp(num, exp);
}
}
template<typename T>
bool almostEqual(const T a, const T b, const unsigned ulps=4)
{
// Handle NaN.
if (std::isnan(a) || std::isnan(b))
return false;
typedef std::numeric_limits<T> limits;
// Handle very small and exactly equal values.
if (std::abs(a-b) <= ulps * limits::denorm_min())
return true;
// frexp() does the wrong thing for zero. But if we get this far
// and either number is zero, then the other is too big, so just
// handle that now.
if (a == 0 || b == 0)
return false;
// Break the numbers into significand and exponent, sorting them by
// exponent.
int min_exp, max_exp;
T min_frac = my_frexp(a, &min_exp);
T max_frac = my_frexp(b, &max_exp);
if (min_exp > max_exp)
{
std::swap(min_frac, max_frac);
std::swap(min_exp, max_exp);
}
// Convert the smaller to the scale of the larger by adjusting its
// significand.
const T scaled_min_frac = std::ldexp(min_frac, min_exp-max_exp);
// Since the significands are now in the same scale, and the larger
// is in the range [0.5, 1), 1 ulp is just epsilon/2.
return std::abs(max_frac-scaled_min_frac) <= ulps * limits::epsilon() / 2;
}
このコードは、(a) 関連するすべてのケースを処理し、(b) IEEE-754 単精度および倍精度の Google 実装と同じことを行い、(c) 完全に標準の C++ であると主張します。
これらの主張の 1 つまたは複数は、ほぼ確実に間違っています。できれば修正を加えて、そのようなことを示す回答を受け入れます。適切な回答には、次の 1 つ以上が含まれている必要があります。
- 最後の場所の単位よりも
ulps大きく異なるが、この関数が true を返す特定の入力 (差が大きいほど良い) ulps最終位の単位までは異なるが、この関数が false を返す特定の入力(差が小さいほど良い)- 私が見逃したすべてのケース
- このコードが未定義の動作に依存する方法、または実装定義の動作に応じて壊れる方法。（可能であれば、関連する仕様を引用してください。）
- 特定した問題の修正
- コードを壊さずに単純化する方法
私は、この問題に重大な報奨金を与えるつもりです。