Каков наиболее эффективный способ для плавающего и двойного сравнения?

Сравнение с эпсилонным значением – это то, что делают большинство людей (даже в игровом программировании).

Вы должны немного изменить свою реализацию:

 bool AreSame(double a, double b) { return fabs(a - b) < EPSILON; } 

Изменить: Christer добавила стопку большой информации по этой теме в недавнем сообщении в блоге . Наслаждаться.

Я обнаружил, что Google C ++ Testing Framework содержит отличную кросс-платформенную реализацию на основе шаблонов AlmostEqual2sComplement, которая работает как по удвоению, так и по плаванию. Учитывая, что он выпущен под лицензией BSD, использование его в вашем собственном коде не должно быть проблемой, если вы сохраняете лицензию. Я извлек приведенный ниже код из http://code.google.com/p/googletest/source/browse/trunk/include/gtest/internal/gtest-internal.h https://github.com/google/googletest/blob /master/googletest/include/gtest/internal/gtest-internal.h и добавлена ​​лицензия сверху.

Обязательно #define GTEST_OS_WINDOWS к некоторому значению (или изменить код, в котором он используется для чего-то, что соответствует вашей кодовой базе), это лицензия BSD в конце концов).

Пример использования:

 double left = // something double right = // something const FloatingPoint lhs(left), rhs(right); if (lhs.AlmostEquals(rhs)) { //they're equal! } 

Вот код:

 // Copyright 2005, Google Inc. // All rights reserved. // // Redistribution and use in source and binary forms, with or without // modification, are permitted provided that the following conditions are // met: // // * Redistributions of source code must retain the above copyright // notice, this list of conditions and the following disclaimer. // * Redistributions in binary form must reproduce the above // copyright notice, this list of conditions and the following disclaimer // in the documentation and/or other materials provided with the // distribution. // * Neither the name of Google Inc. nor the names of its // contributors may be used to endorse or promote products derived from // this software without specific prior written permission. // // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT // OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT // LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, // DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY // THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT // (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. // // Authors: [email protected] (Zhanyong Wan), [email protected] (Sean Mcafee) // // The Google C++ Testing Framework (Google Test) // This template class serves as a compile-time function from size to // type. It maps a size in bytes to a primitive type with that // size. eg // // TypeWithSize<4>::UInt // // is typedef-ed to be unsigned int (unsigned integer made up of 4 // bytes). // // Such functionality should belong to STL, but I cannot find it // there. // // Google Test uses this class in the implementation of floating-point // comparison. // // For now it only handles UInt (unsigned int) as that's all Google Test // needs. Other types can be easily added in the future if need // arises. template  class TypeWithSize { public: // This prevents the user from using TypeWithSize with incorrect // values of N. typedef void UInt; }; // The specialization for size 4. template <> class TypeWithSize<4> { public: // unsigned int has size 4 in both gcc and MSVC. // // As base/basictypes.h doesn't compile on Windows, we cannot use // uint32, uint64, and etc here. typedef int Int; typedef unsigned int UInt; }; // The specialization for size 8. template <> class TypeWithSize<8> { public: #if GTEST_OS_WINDOWS typedef __int64 Int; typedef unsigned __int64 UInt; #else typedef long long Int; // NOLINT typedef unsigned long long UInt; // NOLINT #endif // GTEST_OS_WINDOWS }; // This template class represents an IEEE floating-point number // (either single-precision or double-precision, depending on the // template parameters). // // The purpose of this class is to do more sophisticated number // comparison. (Due to round-off error, etc, it's very unlikely that // two floating-points will be equal exactly. Hence a naive // comparison by the == operation often doesn't work.) // // Format of IEEE floating-point: // // The most-significant bit being the leftmost, an IEEE // floating-point looks like // // sign_bit exponent_bits fraction_bits // // Here, sign_bit is a single bit that designates the sign of the // number. // // For float, there are 8 exponent bits and 23 fraction bits. // // For double, there are 11 exponent bits and 52 fraction bits. // // More details can be found at // http://en.wikipedia.org/wiki/IEEE_floating-point_standard. // // Template parameter: // // RawType: the raw floating-point type (either float or double) template  class FloatingPoint { public: // Defines the unsigned integer type that has the same size as the // floating point number. typedef typename TypeWithSize::UInt Bits; // Constants. // # of bits in a number. static const size_t kBitCount = 8*sizeof(RawType); // # of fraction bits in a number. static const size_t kFractionBitCount = std::numeric_limits::digits - 1; // # of exponent bits in a number. static const size_t kExponentBitCount = kBitCount - 1 - kFractionBitCount; // The mask for the sign bit. static const Bits kSignBitMask = static_cast(1) << (kBitCount - 1); // The mask for the fraction bits. static const Bits kFractionBitMask = ~static_cast(0) >> (kExponentBitCount + 1); // The mask for the exponent bits. static const Bits kExponentBitMask = ~(kSignBitMask | kFractionBitMask); // How many ULP's (Units in the Last Place) we want to tolerate when // comparing two numbers. The larger the value, the more error we // allow. A 0 value means that two numbers must be exactly the same // to be considered equal. // // The maximum error of a single floating-point operation is 0.5 // units in the last place. On Intel CPU's, all floating-point // calculations are done with 80-bit precision, while double has 64 // bits. Therefore, 4 should be enough for ordinary use. // // See the following article for more details on ULP: // http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm. static const size_t kMaxUlps = 4; // Constructs a FloatingPoint from a raw floating-point number. // // On an Intel CPU, passing a non-normalized NAN (Not a Number) // around may change its bits, although the new value is guaranteed // to be also a NAN. Therefore, don't expect this constructor to // preserve the bits in x when x is a NAN. explicit FloatingPoint(const RawType& x) { u_.value_ = x; } // Static methods // Reinterprets a bit pattern as a floating-point number. // // This function is needed to test the AlmostEquals() method. static RawType ReinterpretBits(const Bits bits) { FloatingPoint fp(0); fp.u_.bits_ = bits; return fp.u_.value_; } // Returns the floating-point number that represent positive infinity. static RawType Infinity() { return ReinterpretBits(kExponentBitMask); } // Non-static methods // Returns the bits that represents this number. const Bits &bits() const { return u_.bits_; } // Returns the exponent bits of this number. Bits exponent_bits() const { return kExponentBitMask & u_.bits_; } // Returns the fraction bits of this number. Bits fraction_bits() const { return kFractionBitMask & u_.bits_; } // Returns the sign bit of this number. Bits sign_bit() const { return kSignBitMask & u_.bits_; } // Returns true iff this is NAN (not a number). bool is_nan() const { // It's a NAN if the exponent bits are all ones and the fraction // bits are not entirely zeros. return (exponent_bits() == kExponentBitMask) && (fraction_bits() != 0); } // Returns true iff this number is at most kMaxUlps ULP's away from // rhs. In particular, this function: // // - returns false if either number is (or both are) NAN. // - treats really large numbers as almost equal to infinity. // - thinks +0.0 and -0.0 are 0 DLP's apart. bool AlmostEquals(const FloatingPoint& rhs) const { // The IEEE standard says that any comparison operation involving // a NAN must return false. if (is_nan() || rhs.is_nan()) return false; return DistanceBetweenSignAndMagnitudeNumbers(u_.bits_, rhs.u_.bits_) <= kMaxUlps; } private: // The data type used to store the actual floating-point number. union FloatingPointUnion { RawType value_; // The raw floating-point number. Bits bits_; // The bits that represent the number. }; // Converts an integer from the sign-and-magnitude representation to // the biased representation. More precisely, let N be 2 to the // power of (kBitCount - 1), an integer x is represented by the // unsigned number x + N. // // For instance, // // -N + 1 (the most negative number representable using // sign-and-magnitude) is represented by 1; // 0 is represented by N; and // N - 1 (the biggest number representable using // sign-and-magnitude) is represented by 2N - 1. // // Read http://en.wikipedia.org/wiki/Signed_number_representations // for more details on signed number representations. static Bits SignAndMagnitudeToBiased(const Bits &sam) { if (kSignBitMask & sam) { // sam represents a negative number. return ~sam + 1; } else { // sam represents a positive number. return kSignBitMask | sam; } } // Given two numbers in the sign-and-magnitude representation, // returns the distance between them as an unsigned number. static Bits DistanceBetweenSignAndMagnitudeNumbers(const Bits &sam1, const Bits &sam2) { const Bits biased1 = SignAndMagnitudeToBiased(sam1); const Bits biased2 = SignAndMagnitudeToBiased(sam2); return (biased1 >= biased2) ? (biased1 - biased2) : (biased2 - biased1); } FloatingPointUnion u_; }; 

EDIT: Это сообщение 4 года. Это, вероятно, все еще актуально, и код хорош, но некоторые люди нашли улучшения. Лучше всего получить последнюю версию AlmostEquals прямо из исходного кода Google Test, а не тот, который я вставил здесь.

Сравнение чисел с плавающей запятой зависит от контекста. Поскольку даже изменение порядка операций может привести к различным результатам, важно знать, как «равно» вы хотите, чтобы цифры были.

Сравнение чисел с плавающей точкой от Bruce Dawson – хорошее место для начала при сравнении с плавающей точкой.

Следующие определения относятся к искусству компьютерного программирования Кнутом :

 bool approximatelyEqual(float a, float b, float epsilon) { return fabs(a - b) <= ( (fabs(a) < fabs(b) ? fabs(b) : fabs(a)) * epsilon); } bool essentiallyEqual(float a, float b, float epsilon) { return fabs(a - b) <= ( (fabs(a) > fabs(b) ? fabs(b) : fabs(a)) * epsilon); } bool definitelyGreaterThan(float a, float b, float epsilon) { return (a - b) > ( (fabs(a) < fabs(b) ? fabs(b) : fabs(a)) * epsilon); } bool definitelyLessThan(float a, float b, float epsilon) { return (b - a) > ( (fabs(a) < fabs(b) ? fabs(b) : fabs(a)) * epsilon); } 

Конечно, выбор epsilon зависит от контекста и определяет, насколько равным вы хотите, чтобы номера были.

Другим методом сравнения чисел с плавающей запятой является просмотр ULP (единиц в последнем месте) чисел. Хотя речь идет не только о сравнении, то статья, которую каждый компьютерный ученый должна знать о числах с плавающей запятой, является хорошим ресурсом для понимания того, как работает точка с плавающей запятой и каковы подводные камни, в том числе и то, что представляет собой ULP.

Для более глубокого подхода прочитайте Сравнение чисел с плавающей запятой . Вот fragment кода из этой ссылки:

 // Usable AlmostEqual function bool AlmostEqual2sComplement(float A, float B, int maxUlps) { // Make sure maxUlps is non-negative and small enough that the // default NAN won't compare as equal to anything. assert(maxUlps > 0 && maxUlps < 4 * 1024 * 1024); int aInt = *(int*)&A; // Make aInt lexicographically ordered as a twos-complement int if (aInt < 0) aInt = 0x80000000 - aInt; // Make bInt lexicographically ordered as a twos-complement int int bInt = *(int*)&B; if (bInt < 0) bInt = 0x80000000 - bInt; int intDiff = abs(aInt - bInt); if (intDiff <= maxUlps) return true; return false; } 

Портативный способ получить epsilon в C ++ – это

 #include  std::numeric_limits::epsilon() 

Тогда функция сравнения становится

 #include  #include  bool AreSame(double a, double b) { return std::fabs(a - b) < std::numeric_limits::epsilon(); } 

Понимая, что это старый stream, но эта статья является одной из самых прямых, которые я нашел при сравнении чисел с плавающей запятой, и если вы хотите изучить больше, у нее есть более подробные ссылки, а также основной сайт охватывает весь спектр проблем с числами с плавающей запятой Руководство по плавающей запятой: Сравнение .

Мы можем найти несколько более практичную статью в допуске с плавающей точкой, и есть замечание о абсолютном допуске , которое сводится к этому в C ++:

 bool absoluteToleranceCompare(double x, double y) { return std::fabs(x - y) <= std::numeric_limits::epsilon() ; } 

и относительной толерантности :

 bool relativeToleranceCompare(double x, double y) { double maxXY = std::max( std::fabs(x) , std::fabs(y) ) ; return std::fabs(x - y) <= std::numeric_limits::epsilon()*maxXY ; } 

В статье отмечается, что абсолютный тест терпит неудачу, когда x и y являются большими и терпят неудачу в относительном случае, когда они малы. Предполагая, что абсолютная и относительная толерантность одинакова, комбинированный тест будет выглядеть следующим образом:

 bool combinedToleranceCompare(double x, double y) { double maxXYOne = std::max( { 1.0, std::fabs(x) , std::fabs(y) } ) ; return std::fabs(x - y) <= std::numeric_limits::epsilon()*maxXYOne ; } 

Код, который вы написали, прослушивается:

 return (diff < EPSILON) && (-diff > EPSILON); 

Правильный код:

 return (diff < EPSILON) && (diff > -EPSILON); 

(… и да это другое)

Я удивляюсь, если фабрики не заставят вас потерять ленивую оценку в каком-то случае. Я бы сказал, это зависит от компилятора. Возможно, вы захотите попробовать оба. Если они в среднем эквивалентны, возьмите реализацию с фабриками.

Если у вас есть информация о том, какой из двух поплавков, скорее всего, будет больше, чем другой, вы можете играть по порядку сравнения, чтобы лучше использовать ленивую оценку.

Наконец, вы можете получить лучший результат, добавив эту функцию. Скорее всего, не улучшится, хотя …

Edit: OJ, спасибо за исправление кода. Я удалил свой комментарий соответственно

`return fabs (a – b)

Это нормально, если:

• порядок величин ваших входов существенно не меняется
• очень малые числа противоположных знаков можно рассматривать как равные

Но в противном случае это приведет вас к неприятностям. Номера с двойной точностью имеют разрешение около 16 знаков после запятой. Если два числа, которые вы сравниваете, больше по величине, чем EPSILON * 1.0E16, тогда вы также можете сказать:

 return a==b; 

Я рассмотрю другой подход, предполагающий, что вам нужно беспокоиться о первом выпуске и предположить, что второе отлично подходит для вашего приложения. Решение будет выглядеть примерно так:

 #define VERYSMALL (1.0E-150) #define EPSILON (1.0E-8) bool AreSame(double a, double b) { double absDiff = fabs(a - b); if (absDiff < VERYSMALL) { return true; } double maxAbs = max(fabs(a) - fabs(b)); return (absDiff/maxAbs) < EPSILON; } 

Это дорого стоит вычислительно, но иногда это требует. Это то, что мы должны делать в своей компании, потому что мы имеем дело с инженерной библиотекой, и входы могут варьироваться на несколько десятков порядков.

Во всяком случае, дело в этом (и применимо практически к любой проблеме программирования): Оцените свои потребности, затем придумайте решение для удовлетворения ваших потребностей - не считайте, что простой ответ будет отвечать вашим потребностям. Если после вашей оценки вы обнаружите, что fabs(ab) < EPSILON будет достаточным, идеальным - используйте его! Но имейте в виду его недостатки и другие возможные решения.

Как указывали другие, использование эпсилона фиксированной экспоненты (например, 0,0000001) будет бесполезным для значений от значения эпсилона. Например, если ваши два значения равны 10000.000977 и 10000, то между этими двумя номерами нет 32-разрядных значений с плавающей запятой – 10000 и 10000.000977 находятся как можно ближе, чем вы можете быть без бит-бит-бит. Здесь эпсилон менее 0,0009 не имеет смысла; вы можете использовать оператор прямого равенства.

Аналогично, поскольку эти два значения приближаются к размеру epsilon, относительная погрешность возрастает до 100%.

Таким образом, попытка смешивания числа фиксированной точки, такого как 0,00001 с плавающей точкой (где показатель произвольности), является бессмысленным упражнением. Это будет работать, только если вы можете быть уверены, что значения операнда лежат в узком домене (то есть, близком к определенному показателю), и если вы правильно выбираете значение epsilon для этого конкретного теста. Если вы вытащите номер из эфира («Эй, 0,00001 мало, так что это должно быть хорошо!»), Вы обречены на числовые ошибки. Я потратил много времени на отладку плохого числового кода, где некоторые бедные шмаки бросают в случайных значениях эпсилона, чтобы сделать еще один тестовый пример.

Если вы выполняете численное программирование любого типа и считаете, что вам нужно достичь эпсилон с фиксированной точкой, СТАТЬЯ ЧИТАТЬ БРЮССА О СРАВНЕНИИ ПЛОСКОГО ЧИСЛА .

Сравнение чисел с плавающей запятой

Я закончил тем, что тратил довольно много времени на материал в этой большой теме. Я сомневаюсь, что каждый хочет провести так много времени, поэтому я хотел бы выделить резюме того, что я узнал, и решение, которое я внедрил.

Краткое резюме

1. Есть две проблемы, с которыми сравниваются float: у вас ограниченная точность, а значение «приблизительно нуля» зависит от контекста (см. Следующую точку).
2. Является ли 1E-8 примерно таким же, как 1E-16? Если вы смотрите на шумные данные датчиков, то, вероятно, да, но если вы делаете молекулярное моделирование, то может и не быть! Итог: вам всегда нужно думать о допустимой величине в контексте конкретного вызова функции, а не просто делать ее универсальной жестко-кодированной константой приложения.
3. Для общих функций библиотеки все же приятно иметь параметр с допустимым значением по умолчанию . Типичным выбором является numeric_limits::epsilon() который аналогичен FLT_EPSILON в float.h. Это, однако, проблематично, потому что epsilon для сравнения значений типа 1.0, если не такой же, как epsilon, для значений, подобных 1E9. FLT_EPSILON определен для 1.0.
4. Очевидная реализация для проверки того, находится ли число в пределах допуска, – fabs(ab) <= epsilon однако это не работает, потому что по умолчанию epsilon определен для 1.0. Нам нужно масштабировать epsilon вверх или вниз в терминах a и b.
5. Есть два решения этой проблемы: либо вы устанавливаете epsilon пропорционально max(a,b) либо можете получить следующие отображаемые числа вокруг a, а затем посмотреть, попадает ли b в этот диапазон. Первый называется «относительным» методом, а позже называется методом ULP.
6. Оба метода фактически терпят неудачу в любом случае при сравнении с 0. В этом случае приложение должно предоставить правильный допуск.

Внедрение служебных функций (C ++ 11)

 //implements relative method - do not use for comparing with zero //use this most of the time, tolerance needs to be meaningful in your context template static bool isApproximatelyEqual(TReal a, TReal b, TReal tolerance = std::numeric_limits::epsilon()) { TReal diff = std::fabs(a - b); if (diff <= tolerance) return true; if (diff < std::fmax(std::fabs(a), std::fabs(b)) * tolerance) return true; return false; } //supply tolerance that is meaningful in your context //for example, default tolerance may not work if you are comparing double with float template static bool isApproximatelyZero(TReal a, TReal tolerance = std::numeric_limits::epsilon()) { if (std::fabs(a) <= tolerance) return true; return false; } //use this when you want to be on safe side //for example, don't start rover unless signal is above 1 template static bool isDefinitelyLessThan(TReal a, TReal b, TReal tolerance = std::numeric_limits::epsilon()) { TReal diff = a - b; if (diff < tolerance) return true; if (diff < std::fmax(std::fabs(a), std::fabs(b)) * tolerance) return true; return false; } template static bool isDefinitelyGreaterThan(TReal a, TReal b, TReal tolerance = std::numeric_limits::epsilon()) { TReal diff = a - b; if (diff > tolerance) return true; if (diff > std::fmax(std::fabs(a), std::fabs(b)) * tolerance) return true; return false; } //implements ULP method //use this when you are only concerned about floating point precision issue //for example, if you want to see if a is 1.0 by checking if its within //10 closest representable floating point numbers around 1.0. template static bool isWithinPrecisionInterval(TReal a, TReal b, unsigned int interval_size = 1) { TReal min_a = a - (a - std::nextafter(a, std::numeric_limits::lowest())) * interval_size; TReal max_a = a + (std::nextafter(a, std::numeric_limits::max()) - a) * interval_size; return min_a <= b && max_a >= b; } 

General-purpose comparison of floating-point numbers is generally meaningless. How to compare really depends on a problem at hand. In many problems, numbers are sufficiently discretized to allow comparing them within a given tolerance. Unfortunately, there are just as many problems, where such trick doesn’t really work. For one example, consider working with a Heaviside (step) function of a number in question (digital stock options come to mind) when your observations are very close to the barrier. Performing tolerance-based comparison wouldn’t do much good, as it would effectively shift the issue from the original barrier to two new ones. Again, there is no general-purpose solution for such problems and the particular solution might require going as far as changing the numerical method in order to achieve stability.

Unfortunately, even your “wasteful” code is incorrect. EPSILON is the smallest value that could be added to 1.0 and change its value. The value 1.0 is very important — larger numbers do not change when added to EPSILON. Now, you can scale this value to the numbers you are comparing to tell whether they are different or not. The correct expression for comparing two doubles is:

 if (fabs(a - b) <= DBL_EPSILON * fmax(fabs(a), fabs(b))) { // ... } 

This is at a minimum. In general, though, you would want to account for noise in your calculations and ignore a few of the least significant bits, so a more realistic comparison would look like:

 if (fabs(a - b) <= 16 * DBL_EPSILON * fmax(fabs(a), fabs(b))) { // ... } 

If comparison performance is very important to you and you know the range of your values, then you should use fixed-point numbers instead.

My class based on previously posted answers. Very similar to Google’s code but I use a bias which pushes all NaN values above 0xFF000000. That allows a faster check for NaN.

This code is meant to demonstrate the concept, not be a general solution. Google’s code already shows how to compute all the platform specific values and I didn’t want to duplicate all that. I’ve done limited testing on this code.

 typedef unsigned int U32; // Float Memory Bias (unsigned) // ----- ------ --------------- // NaN 0xFFFFFFFF 0xFF800001 // NaN 0xFF800001 0xFFFFFFFF // -Infinity 0xFF800000 0x00000000 --- // -3.40282e+038 0xFF7FFFFF 0x00000001 | // -1.40130e-045 0x80000001 0x7F7FFFFF | // -0.0 0x80000000 0x7F800000 |--- Valid <= 0xFF000000. // 0.0 0x00000000 0x7F800000 | NaN > 0xFF000000 // 1.40130e-045 0x00000001 0x7F800001 | // 3.40282e+038 0x7F7FFFFF 0xFEFFFFFF | // Infinity 0x7F800000 0xFF000000 --- // NaN 0x7F800001 0xFF000001 // NaN 0x7FFFFFFF 0xFF7FFFFF // // Either value of NaN returns false. // -Infinity and +Infinity are not "close". // -0 and +0 are equal. // class CompareFloat{ public: union{ float m_f32; U32 m_u32; }; static bool CompareFloat::IsClose( float A, float B, U32 unitsDelta = 4 ) { U32 a = CompareFloat::GetBiased( A ); U32 b = CompareFloat::GetBiased( B ); if ( (a > 0xFF000000) || (b > 0xFF000000) ) { return( false ); } return( (static_cast(abs( a - b ))) < unitsDelta ); } protected: static U32 CompareFloat::GetBiased( float f ) { U32 r = ((CompareFloat*)&f)->m_u32; if ( r & 0x80000000 ) { return( ~r - 0x007FFFFF ); } return( r + 0x7F800000 ); } }; 

It depends on how precise you want the comparison to be. If you want to compare for exactly the same number, then just go with ==. (You almost never want to do this unless you actually want exactly the same number.) On any decent platform you can also do the following:

 diff= a - b; return fabs(diff)