Jaki jest najskuteczniejszy sposób porównywania zmiennoprzecinkowego i podwójnego?


524

Jaki byłby najbardziej efektywny sposób porównania dwóch doublelub dwóch floatwartości?

Po prostu robienie tego jest nieprawidłowe:

bool CompareDoubles1 (double A, double B)
{
   return A == B;
}

Ale coś takiego:

bool CompareDoubles2 (double A, double B) 
{
   diff = A - B;
   return (diff < EPSILON) && (-diff < EPSILON);
}

Wydaje się do przetwarzania odpadów.

Czy ktoś zna mądrzejszego porównywania pływaka?


2
> czy bardziej efektywne byłoby dodanie ... na początku funkcji? <invoke Knuth>Przedwczesna optymalizacja jest źródłem wszelkiego zła. </invoke Knuth>Po prostu skorzystaj z abs (ab) <EPS, jak wspomniano powyżej, jest to jasne i łatwe do zrozumienia.
Andrew Coleson,


2
Jedyną nieoptymalną implementacją oryginalnego plakatu jest to, że zawiera on dodatkowy oddział w &&. Odpowiedź Dz.U. jest optymalna. fabs jest wewnętrzną funkcją, która jest pojedynczą instrukcją na x87 i przypuszczam, że na prawie wszystko inne. Zaakceptuj już odpowiedź OJ!
3lat

3
Jeśli możesz, upuść zmiennoprzecinkowy i użyj punktów stałych. Przykład: użyj {stałego punktu} milimetrów zamiast {zmiennoprzecinkowego} metra.
Thomas Matthews,

33
„Po prostu robienie tego nie jest poprawne” - To zwykłe śmieci, oczywiście używanie ==może być całkowicie poprawne, ale to całkowicie zależy od kontekstu nie podanego w pytaniu. Do czasu poznania tego kontekstu ==nadal pozostaje „najbardziej wydajnym sposobem” .
Christian Rau,

Odpowiedzi:


459

Zachowaj szczególną ostrożność, korzystając z innych sugestii. Wszystko zależy od kontekstu.

Spędziłem dużo czasu na śledzeniu błędów w systemie, który zakładał, że a==bjeśli |a-b|<epsilon. Podstawowymi problemami były:

  1. Domniemane domniemanie w algorytmie, że jeśli a==bi b==cwtedy a==c.

  2. Korzystanie z tego samego epsilonu dla linii mierzonych w calach i linii mierzonych w milsach (0,001 cala). To jest a==bale 1000a!=1000b. (Właśnie dlatego AlmostEqual2sComplement prosi o epsilon lub max ULPS).

  3. Zastosowanie tego samego epsilonu zarówno dla cosinusa kątów, jak i długości linii!

  4. Korzystanie z takiej funkcji porównywania do sortowania elementów w kolekcji. (W tym przypadku użycie wbudowanego operatora C ++ == dla podwójnych dało poprawne wyniki.)

Tak jak powiedziałem: wszystko zależy od kontekstu i oczekiwanego rozmiaru ai b.

BTW, std::numeric_limits<double>::epsilon() to „maszyna epsilon”. Jest to różnica między 1,0 a następną wartością reprezentowaną przez podwójny. Myślę, że można go użyć w funkcji porównania, ale tylko wtedy, gdy oczekiwane wartości są mniejsze niż 1. (Jest to odpowiedź na odpowiedź @ cdv ...)

Ponadto, jeśli w zasadzie masz intarytmetykę doubles(tutaj używamy podwójnych wartości do utrzymania wartości int w niektórych przypadkach), twoja arytmetyka będzie poprawna. Na przykład 4.0 / 2.0 będzie taki sam jak 1.0 + 1.0. Tak długo, jak nie robisz rzeczy, które powodują frakcje (4.0 / 3.0) lub nie wychodzą poza rozmiar int.


10
+1 za wskazanie rzeczy oczywistych (które często są ignorowane). W przypadku metody ogólnej można ustawić epsilon względem, fabs(a)+fabs(b)ale z kompensacją NaN, sumy 0 i przelewu, staje się to dość skomplikowane.
peterchen

4
Musi być coś, czego nie rozumiem. Typowy float/ doubleto MANTISSA x 2 ^ EXP . epsilonbędzie zależeć od wykładnika. Na przykład jeśli mantysa ma 24 bity, a wykładnik ma 8 bitów, to dla niektórych wartości jest 1/(2^24)*2^127lub ; czy to odnosi się do minimalnego epsilon ? ~2^103epsilon
artless noise

3
Poczekaj sekundę. Czy to, co powiedziałem, masz na myśli? Mówisz dlaczego |a-b|<epsilon, to nie jest poprawne. Dodaj ten link do swojej odpowiedzi; jeśli zgadzasz się z cygnus-software.com/papers/comparingfloats/comparingfloats.htm i mogę usunąć moje głupie komentarze.
artless noise

3
To bardzo długi komentarz, a nie odpowiedź sama w sobie. Czy istnieje (zbiór) kanonicznych odpowiedzi na wszystkie konteksty?
Merlyn Morgan-Graham

2
Stary link wydaje się przestarzały, nowa strona jest tutaj randomascii.wordpress.com/2012/02/25/…
Marson Mao

174

Porównanie z wartością epsilon jest tym, co robi większość ludzi (nawet w programowaniu gier).

Powinieneś jednak nieco zmienić swoją implementację:

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

Edycja: Christer dodał stos świetnych informacji na ten temat w ostatnim poście na blogu . Cieszyć się.


@OJ: czy coś jest nie tak z pierwszą próbką kodu? Myślałem, że jedynym problemem jest sytuacja taka jak ta: float a = 3.4; if(a == 3.4){...}tzn. Kiedy porównujesz zapisany zmiennoprzecinkowy z literałem | W takim przypadku obie liczby są przechowywane, więc będą miały taką samą reprezentację, jeśli będą równe, więc jaka szkoda to zrobić a == b?
Lazer

11
@DonReba: Tylko jeśli EPSILONjest zdefiniowany jako DBL_EPSILON. Zwykle będzie to konkretna wartość wybierana w zależności od wymaganej dokładności porównania.
Nemo157,

7
EPSILONporównanie nie działa, gdy zmiennoprzecinkowe są duże, ponieważ różnica między kolejnymi zmiennoprzecinkowymi również staje się duża. Zobacz ten artykuł .
kevintodisco

22
Nic dziwnego, że w niektórych grach dochodzi do walki z Z, gdy tekstury / obiekty migają daleko, jak w Battlefield 4. Porównywanie różnic EPSILONjest prawie bezużyteczne. Musisz porównać z progiem, który ma sens dla dostępnych jednostek. Użyj również, std::absponieważ jest przeciążony dla różnych typów zmiennoprzecinkowych.
Maxim Egorushkin

11
Przegłosowałem, ponieważ przykładowy kod pokazuje typowy błąd, który jest powtarzany przez większość programistów. Punkt zmiennoprzecinkowy zawsze dotyczy błędów względnych, ponieważ jest to zmiennoprzecinkowy (nie stały punkt). Tak więc nigdy nie będzie działać poprawnie ze stałym błędem (epsilon).
user2261015

115

Odkryłem, że platforma testowa Google C ++ zawiera ładną wieloplatformową implementację AlmostEqual2sComplement, która działa zarówno na podwójnych, jak i zmiennoprzecinkowych. Biorąc pod uwagę, że jest on wydany na licencji BSD, użycie go we własnym kodzie nie powinno stanowić problemu, pod warunkiem zachowania licencji. Wyodrębniłem poniższy kod z 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 i dodał licencję na górze.

Pamiętaj, aby # zdefiniować GTEST_OS_WINDOWS do pewnej wartości (lub zmienić kod, w którym jest on przyzwyczajony, na coś, co pasuje do twojej bazy kodowej - w końcu jest to licencja BSD).

Przykład użycia:

double left  = // something
double right = // something
const FloatingPoint<double> lhs(left), rhs(right);

if (lhs.AlmostEquals(rhs)) {
  //they're equal!
}

Oto kod:

// 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: wan@google.com (Zhanyong Wan), eefacm@gmail.com (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. e.g.
//
//   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 <size_t size>
class TypeWithSize {
 public:
  // This prevents the user from using TypeWithSize<N> 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 <typename RawType>
class FloatingPoint {
 public:
  // Defines the unsigned integer type that has the same size as the
  // floating point number.
  typedef typename TypeWithSize<sizeof(RawType)>::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<RawType>::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<Bits>(1) << (kBitCount - 1);

  // The mask for the fraction bits.
  static const Bits kFractionBitMask =
    ~static_cast<Bits>(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_;
};

EDYCJA: Ten post ma 4 lata. Prawdopodobnie jest nadal poprawny, a kod jest fajny, ale niektórzy znaleźli ulepszenia. Najlepiej skorzystaj z najnowszej wersji AlmostEqualskodu źródłowego Google Test, a nie tej, którą tutaj wkleiłem.


3
+1: Zgadzam się, że to jest poprawne. Jednak nie wyjaśnia dlaczego. Zobacz tutaj: cygnus-software.com/papers/comparingfloats/comparingfloats.htm Przeczytałem ten post na blogu po tym, jak napisałem tutaj komentarz do najlepszego wyniku; Wierzę, że mówi to samo i zapewnia racjonalne / rozwiązanie, które jest zaimplementowane powyżej. Ponieważ jest tak dużo kodu, ludzie przegapią odpowiedź.
artless noise

Jest kilka paskudnych rzeczy, które mogą się zdarzyć, gdy dochodzi do niejawnych rzutów, powiedzmy FloatPoint <double> fp (0.03f). Wprowadziłem kilka modyfikacji, aby temu zapobiec. szablon <typename U> wyraźny FloatingPoint (const U & x) {if (typeid (U) .name ()! = = typid (RawType) .name ()) {std :: cerr << "Dokonujesz niejawnej konwersji za pomocą FloatingPoint, Don't "<< std :: endl; assert (typeid (U) .name () == typeid (RawType) .name ()); } u_.value_ = x; }
JeffCharter

2
Dobre znalezisko! Wydaje mi się jednak, że najlepiej byłoby przekazać je do testu Google, skąd ten kod został skradziony. Zaktualizuję post, aby odzwierciedlić, że prawdopodobnie jest nowsza wersja. Jeśli faceci Google zachowują się swędząco, czy możesz umieścić to np. W GitHub? Link do tego też też.
skrebbel

3
Najnowszy fragment kodu znajduje się tutaj i tutaj .
Jaege

1
Wyodrębniłem niezbędne wiersze do pliku gist. Każdy może stąd dotrzeć .
Yusuf Tarık Günaydın

111

Porównanie liczb zmiennoprzecinkowych dla zależy od kontekstu. Ponieważ nawet zmiana kolejności operacji może dawać różne wyniki, ważne jest, aby wiedzieć, jak „równe” mają być liczby.

Porównywanie liczb zmiennoprzecinkowych przez Bruce'a Dawsona jest dobrym miejscem do rozpoczęcia, patrząc na porównanie liczb zmiennoprzecinkowych.

Poniższe definicje pochodzą ze sztuki programowania komputerowego Knutha :

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);
}

Oczywiście wybór epsilon zależy od kontekstu i określa, jak równe mają być liczby.

Inną metodą porównywania liczb zmiennoprzecinkowych jest spojrzenie na ULP (jednostki na ostatnim miejscu) liczb. Nie zajmując się konkretnie porównaniami, artykuł To, co każdy informatyk powinien wiedzieć o liczbach zmiennoprzecinkowych, jest dobrym źródłem informacji na temat tego, jak działa zmiennoprzecinkowy i jakie są pułapki, w tym czym jest ULP.


1
dzięki za opublikowanie, jak ustalić, która liczba jest mniejsza / większa!
Pomidor,

1
fabs(a - b) <= ( (fabs(a) < fabs(b) ? fabs(b) : fabs(a)) * epsilon);uratował mi życie. LOL Zauważ, że ta wersja (nie sprawdziłem, czy dotyczy również innych) również uwzględnia zmianę, która może wystąpić w integralnej części liczby zmiennoprzecinkowej (przykład: 2147352577.9999997616 == 2147352576.0000000000gdzie wyraźnie widać, że istnieje prawie różnica między 2między dwiema liczbami), co jest całkiem miłe! Dzieje się tak, gdy skumulowany błąd zaokrąglenia przepełnia dziesiętną część liczby.
rbaleksandar

Bardzo miły i pomocny artykuł Bruce'a Dawsona, dzięki!
BobMorane,

2
Biorąc pod uwagę, że to pytanie jest oznaczone jako C ++, twoje czeki byłyby łatwiejsze do odczytania w formie std::max(std::abs(a), std::abs(b))(lub z std::min()); std::absw C ++ jest przeciążony typami float i double, więc działa dobrze (zawsze możesz zachować fabsczytelność).
Razakhel

1
Okazało się, że problem tkwił w moim kodzie, różnica między pierwotną oczekiwaną wartością a analizowanym ciągiem.
mwpowellhtx

47

Aby uzyskać bardziej szczegółowe podejście, przeczytaj Porównanie liczb zmiennoprzecinkowych . Oto fragment kodu z tego linku:

// 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;    
}

14
Jaka jest sugerowana wartość maxUlps?
unj2

6
Czy „ *(int*)&A;” naruszy ścisłą zasadę aliasingu?
osgx

3
Według gtest (poszukiwanie ULP) 4 jest liczbą dopuszczalną.
Maj Oakes

4
A oto kilka aktualizacji artykułu Bruce'a Dawsona (jeden z nich znajduje się we wstępie do artykułu): randomascii.wordpress.com/2012/02/25/... i randomascii.wordpress.com/2012/06/26/...
Michael Burr,

3
Zajęło mi trochę czasu, aby dowiedzieć się, co na ULP było: Jednostki na ostatnim miejscu
JeffCharter

27

Uświadomienie sobie, że jest to stary wątek, ale ten artykuł jest jednym z najprostszych, jakie znalazłem na temat porównywania liczb zmiennoprzecinkowych, a jeśli chcesz dowiedzieć się więcej, ma również bardziej szczegółowe odniesienia, a strona główna zawiera pełny zakres zagadnień radzenie sobie z liczbami zmiennoprzecinkowymi Przewodnik zmiennoprzecinkowy: porównanie .

Możemy znaleźć nieco bardziej praktyczny artykuł w zmienionym zmiennym punkcie tolerancji i zauważa, że ​​istnieje test absolutnej tolerancji , który sprowadza się do tego w C ++:

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

i test tolerancji względnej :

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

Notatki artykule, że bezwzględna test nie powiedzie się, gdy xi ysą duże i nie we względnym przypadku, gdy są one małe. Zakładając, że tolerancja bezwzględna i względna jest taka sama, test łączony wyglądałby tak:

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<double>::epsilon()*maxXYOne ;
}

25

Przenośnym sposobem na uzyskanie epsilon w C ++ jest

#include <limits>
std::numeric_limits<double>::epsilon()

Wtedy staje się funkcja porównania

#include <cmath>
#include <limits>

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

34
Będziesz chciał wielokrotności tego epsilonu najprawdopodobniej.
user7116

11
Nie możesz po prostu użyć std :: abs? AFAIK, std :: abs jest przeciążony również dla gier podwójnych. Ostrzegaj mnie, jeśli się mylę.
kolistivra

3
@kolistivra, mylisz się. Litera „f” w „fabs” nie oznacza typu float. Prawdopodobnie myślisz o funkcjach C fabsf () i fabsl ().
jcoffland

9
W rzeczywistości z powodów opisanych w artykule Bruce'a epsilon zmienia się wraz ze wzrostem wartości zmiennoprzecinkowej. Zobacz część, w której mówi: „W przypadku liczb większych niż 2,0 różnica między zmiennoprzecinkowymi powiększa się, a jeśli porównasz zmiennoprzecinkowe za pomocą FLT_EPSILON, wtedy przeprowadzasz tylko droższą i mniej oczywistą kontrolę równości”.
bobobobo,

5
Wiem, że to stare, ale std :: abs jest przeciążone dla typów zmiennoprzecinkowych w cmath.
mholzmann

18

Skończyłem spędzać sporo czasu na przeglądaniu materiału w tym wielkim wątku. Wątpię, czy wszyscy chcą spędzać tak dużo czasu, dlatego chciałbym podkreślić podsumowanie tego, czego się nauczyłem i rozwiązania, które wdrożyłem.

Szybkie podsumowanie

  1. Czy 1e-8 jest w przybliżeniu taki sam jak 1e-16? Jeśli patrzysz na głośne dane z czujnika, prawdopodobnie tak, ale jeśli wykonujesz symulację molekularną, być może nie! Konkluzja: Zawsze musisz myśleć o wartości tolerancji w kontekście konkretnego wywołania funkcji, a nie tylko uczynić ją ogólną stałą zakodowaną na całej aplikacji.
  2. W przypadku ogólnych funkcji bibliotecznych nadal dobrze jest mieć parametr z domyślną tolerancją . Typowy wybór to numeric_limits::epsilon()taki sam, jak FLT_EPSILON w float.h. Jest to jednak problematyczne, ponieważ epsilon do porównywania wartości takich jak 1.0 nie jest taki sam jak epsilon dla wartości takich jak 1E9. FLT_EPSILON jest zdefiniowany dla 1.0.
  3. Oczywista implementacja sprawdzania, czy liczba mieści się w tolerancji, fabs(a-b) <= epsilonnie działa jednak, ponieważ domyślna epsilon jest zdefiniowana dla 1.0. Musimy skalować epsilon w górę lub w dół w kategoriach a i b.
  4. Istnieją dwa rozwiązania tego problemu: albo ustawisz proporcjonalność epsilon, max(a,b)albo możesz uzyskać kolejne reprezentatywne liczby wokół a, a następnie sprawdź, czy b mieści się w tym zakresie. Ten pierwszy nazywa się metodą „względną”, a później nazywa się metodą ULP.
  5. Obie metody i tak naprawdę zawodzą w porównaniu z 0. W takim przypadku aplikacja musi zapewnić prawidłową tolerancję.

Implementacja funkcji narzędziowych (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<typename TReal>
static bool isApproximatelyEqual(TReal a, TReal b, TReal tolerance = std::numeric_limits<TReal>::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<typename TReal>
static bool isApproximatelyZero(TReal a, TReal tolerance = std::numeric_limits<TReal>::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<typename TReal>
static bool isDefinitelyLessThan(TReal a, TReal b, TReal tolerance = std::numeric_limits<TReal>::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<typename TReal>
static bool isDefinitelyGreaterThan(TReal a, TReal b, TReal tolerance = std::numeric_limits<TReal>::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<typename TReal>
static bool isWithinPrecisionInterval(TReal a, TReal b, unsigned int interval_size = 1)
{
    TReal min_a = a - (a - std::nextafter(a, std::numeric_limits<TReal>::lowest())) * interval_size;
    TReal max_a = a + (std::nextafter(a, std::numeric_limits<TReal>::max()) - a) * interval_size;

    return min_a <= b && max_a >= b;
}

isDefinitelyLessThankontrole diff < tolerance, co oznacza, że ​​aib są prawie równe (a więc a nie jest zdecydowanie mniejsze niż b). Czy nie ma większego sensu sprawdzanie różnicy tolerancji w obu przypadkach? Lub może dodać orEqualToargument, który kontroluje, czy przybliżona kontrola równości powinna zwrócić wartość true, czy nie.
Matt Chambers

14

Napisany kod jest nieprawidłowy:

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

Poprawny kod to:

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

(... i tak, to jest inne)

Zastanawiam się, czy fabs nie sprawi, że w niektórych przypadkach stracisz leniwą ocenę. Powiedziałbym, że to zależy od kompilatora. Możesz spróbować obu. Jeśli są one średnio równoważne, weź implementację z fabs.

Jeśli masz jakieś informacje na temat tego, który z dwóch pływaków jest bardziej prawdopodobny niż inny, możesz grać w kolejności porównania, aby lepiej wykorzystać leniwą ocenę.

Wreszcie możesz uzyskać lepszy wynik, wprowadzając tę ​​funkcję. Prawdopodobnie nie poprawi się znacznie ...

Edycja: OJ, dziękuję za poprawienie kodu. Odpowiednio usunąłem swój komentarz


13

`return fabs (a - b) <EPSILON;

Jest to w porządku, jeśli:

  • rząd wielkości twoich danych wejściowych niewiele się zmienia
  • bardzo małą liczbę przeciwnych znaków można traktować jako równe

Ale w przeciwnym razie popadniesz w kłopoty. Liczby o podwójnej precyzji mają rozdzielczość około 16 miejsc po przecinku. Jeśli dwie liczby, które porównujesz, są większe niż EPSILON * 1.0E16, możesz równie dobrze powiedzieć:

return a==b;

Zbadam inne podejście, które zakłada, że ​​musisz się martwić o pierwszy problem i zakładam, że drugi jest w porządku. Rozwiązaniem byłoby coś takiego:

#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;
}

Jest to drogie obliczeniowo, ale czasem wymaga tego. To właśnie musimy zrobić w mojej firmie, ponieważ mamy do czynienia z biblioteką inżynieryjną, a nakłady mogą się różnić o kilkadziesiąt rzędów wielkości.

W każdym razie chodzi o to (i dotyczy praktycznie każdego problemu programistycznego): oceń swoje potrzeby, a następnie wymyśl rozwiązanie odpowiadające twoim potrzebom - nie zakładaj, że łatwa odpowiedź zaspokoi twoje potrzeby. Jeśli po dokonaniu oceny okaże się, że fabs(a-b) < EPSILONto wystarczy, doskonale - skorzystaj z niego! Pamiętaj jednak o jego wadach i innych możliwych rozwiązaniach.


3
Oprócz literówek (s / - /, / brakuje przecinka w fmax ()), ta implementacja ma błąd dla liczb bliskich zeru, które są w EPSILON, ale jeszcze nie VERYSMALL. Np. AreSame (1.0E-10, 1.0E-9) zgłasza fałsz, ponieważ błąd względny jest ogromny. Możesz zostać bohaterem w swojej firmie.
brlcad

1
@brlcad Nie dostałeś punktu zmiennoprzecinkowego . 1.0E-10 i 1.0E-9 różnią się wielkością 10. Więc to prawda, że ​​nie są takie same. zmiennoprzecinkowa zawsze dotyczy błędów względnych . Jeśli masz system, który uważa 1.0E-10 i 1.0E-9 za prawie równe, ponieważ oba są „całkiem bliskie zeru” (co brzmi rozsądnie dla ludzi, ale nie jest niczym matematycznym), to EPSILON należy odpowiednio dostosować dla takiego systemu.
user2261015,

8

Jak zauważyli inni, użycie epsilonu o stałym wykładniku (np. 0,0000001) będzie bezużyteczne dla wartości odbiegających od wartości epsilon. Na przykład, jeśli dwie wartości to 10000.000977 i 10000, to NIE 32-bitowych wartości zmiennoprzecinkowych między tymi dwiema liczbami - 10000 i 10000.000977 są tak bliskie, jak to tylko możliwe, bez bycia identycznym dla każdego bitu. Tutaj epsilon mniejszy niż 0,0009 jest bez znaczenia; równie dobrze możesz użyć prostego operatora równości.

Podobnie, gdy obie wartości zbliżają się do wielkości epsilon, błąd względny rośnie do 100%.

Zatem próba zmieszania stałej liczby punktów, takiej jak 0,00001 z wartościami zmiennoprzecinkowymi (gdzie wykładnik jest arbitralny), jest bezcelowym ćwiczeniem. To zadziała tylko wtedy, gdy będziesz mieć pewność, że wartości argumentu znajdują się w wąskiej domenie (to znaczy blisko jakiegoś konkretnego wykładnika) i jeśli właściwie wybierzesz wartość epsilon dla tego konkretnego testu. Jeśli wyciągniesz liczbę z powietrza („Hej! 0,00001 jest mały, więc to musi być dobry!”), Jesteś skazany na błędy numeryczne. Spędziłem dużo czasu na debugowaniu złego kodu numerycznego, w którym niektóre kiepskie tępaki rzucają losowymi wartościami epsilon, aby kolejny test zadziałał.

Jeśli wykonujesz programowanie numeryczne i uważasz, że musisz sięgnąć po epsilony o stałym punkcie, PRZECZYTAJ ARTYKUŁ BRUCE NA TEMAT PORÓWNANIA LICZB PUNKTOWYCH .

Porównywanie liczb zmiennoprzecinkowych


5

Qt implementuje dwie funkcje, być może możesz się z nich nauczyć:

static inline bool qFuzzyCompare(double p1, double p2)
{
    return (qAbs(p1 - p2) <= 0.000000000001 * qMin(qAbs(p1), qAbs(p2)));
}

static inline bool qFuzzyCompare(float p1, float p2)
{
    return (qAbs(p1 - p2) <= 0.00001f * qMin(qAbs(p1), qAbs(p2)));
}

Od tego czasu możesz potrzebować następujących funkcji

Zauważ, że porównywanie wartości, w których p1 lub p2 wynosi 0,0, nie będzie działać, podobnie jak porównywanie wartości, w których jedną z wartości jest NaN lub nieskończoność. Jeśli jedną z wartości jest zawsze 0,0, użyj zamiast tego qFuzzyIsNull. Jeśli jedna z wartości prawdopodobnie będzie wynosić 0,0, jednym rozwiązaniem jest dodanie 1,0 do obu wartości.

static inline bool qFuzzyIsNull(double d)
{
    return qAbs(d) <= 0.000000000001;
}

static inline bool qFuzzyIsNull(float f)
{
    return qAbs(f) <= 0.00001f;
}

3

Ogólne porównanie liczb zmiennoprzecinkowych jest na ogół bez znaczenia. Sposób porównania naprawdę zależy od konkretnego problemu. W wielu problemach liczby są wystarczająco dyskrecjonalne, aby umożliwić porównanie ich w ramach danej tolerancji. Niestety jest tyle problemów, w których taka sztuczka tak naprawdę nie działa. Na przykład rozważ rozważenie funkcji Heaviside (krok) dla danej liczby (przychodzą na myśl cyfrowe opcje akcji), gdy obserwacje są bardzo blisko bariery. Przeprowadzenie porównania opartego na tolerancji nie przyniosłoby wiele dobrego, ponieważ skutecznie przeniósłoby problem z pierwotnej bariery na dwie nowe. Ponownie, nie ma ogólnego rozwiązania takich problemów, a konkretne rozwiązanie może wymagać posunięcia aż do zmiany metody numerycznej w celu osiągnięcia stabilności.


3

Niestety, nawet twój „marnotrawny” kod jest niepoprawny. EPSILON to najmniejsza wartość, którą można dodać do 1,0 i zmienić jej wartość. Wartość 1,0 jest bardzo ważna - większe liczby nie zmieniają się po dodaniu do EPSILON. Teraz możesz przeskalować tę wartość do liczb, które porównujesz, aby stwierdzić, czy są różne, czy nie. Prawidłowe wyrażenie do porównania dwóch podwójnych to:

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

To jest minimum. Ogólnie rzecz biorąc, chciałbyś uwzględnić szum w swoich obliczeniach i zignorować kilka najmniej znaczących bitów, więc bardziej realistyczne porównanie wyglądałoby następująco:

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

Jeśli wydajność porównania jest dla Ciebie bardzo ważna i znasz zakres swoich wartości, powinieneś zamiast tego użyć liczb stałych.


2
„EPSILON to najmniejsza wartość, którą można dodać do 1,0 i zmienić jej wartość”: W rzeczywistości ten zaszczyt przypada następcy 0,5 * EPSILON (w domyślnym trybie zaokrąglania do najbliższego). blog.frama-c.com/index.php?post/2013/05/09/FLT_EPSILON
Pascal Cuoq

Jak myślisz, dlaczego EPSILONw pytaniu jest DBL_EPSILONlub FLT_EPSILON? Problem leży w twojej wyobraźni, gdzie podstawiłeś DBL_EPSILON(co byłoby niewłaściwym wyborem) w kodzie, który go nie używał.
Ben Voigt

@BenVoigt, masz rację, w tym momencie coś mi przyszło do głowy i w tym świetle zinterpretowałem to pytanie.
Don Reba,

2

Moja klasa na podstawie wcześniej opublikowanych odpowiedzi. Bardzo podobny do kodu Google, ale używam błędu, który popycha wszystkie wartości NaN powyżej 0xFF000000. To pozwala na szybsze sprawdzenie NaN.

Ten kod ma na celu zademonstrowanie koncepcji, a nie ogólne rozwiązanie. Kod Google już pokazuje, jak obliczyć wszystkie wartości specyficzne dla platformy i nie chciałem tego wszystkiego kopiować. Przeprowadziłem ograniczone testy tego kodu.

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<U32>(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 );
                  }
};

2

Oto dowód, że użycie std::numeric_limits::epsilon()nie jest odpowiedzią - nie powiedzie się dla wartości większych niż jeden:

Dowód mojego komentarza powyżej:

#include <stdio.h>
#include <limits>

double ItoD (__int64 x) {
    // Return double from 64-bit hexadecimal representation.
    return *(reinterpret_cast<double*>(&x));
}

void test (__int64 ai, __int64 bi) {
    double a = ItoD(ai), b = ItoD(bi);
    bool close = std::fabs(a-b) < std::numeric_limits<double>::epsilon();
    printf ("%.16f and %.16f %s close.\n", a, b, close ? "are " : "are not");
}

int main()
{
    test (0x3fe0000000000000L,
          0x3fe0000000000001L);

    test (0x3ff0000000000000L,
          0x3ff0000000000001L);
}

Uruchomienie daje wynik:

0.5000000000000000 and 0.5000000000000001 are  close.
1.0000000000000000 and 1.0000000000000002 are not close.

Zauważ, że w drugim przypadku (jeden i tylko większy niż jeden), dwie wartości wejściowe są tak bliskie, jak to tylko możliwe, i nadal są porównywane jako niezamknięte. Tak więc, dla wartości większych niż 1.0, równie dobrze możesz po prostu użyć testu równości. Naprawione epsilony nie uratują cię podczas porównywania wartości zmiennoprzecinkowych.


Uważam, że return *(reinterpret_cast<double*>(&x));chociaż zwykle działa, jest w rzeczywistości nieokreślonym zachowaniem.
Jaap Versteegh,

Dobry punkt, chociaż ten kod ma charakter poglądowy, więc wystarczający do zademonstrowania problemu numeric_limits<>::epsiloni punktu podłogowego IEEE 754.
Steve Hollasch,

Jest to również słuszna kwestia, ale nie jest mądrym imho publikować informacje o przepełnieniu stosu, oczekując tego rodzaju wglądu. Kod zostanie skopiowany na ślepo, co utrudni wyeliminowanie tego bardzo powszechnego wzorca - razem z trikiem związkowym - którego należy po prostu unikać, tak jak wszystkie UD.
Jaap Versteegh,

1

Znalazłem inną ciekawą implementację na: https://en.cppreference.com/w/cpp/types/numeric_limits/epsilon

#include <cmath>
#include <limits>
#include <iomanip>
#include <iostream>
#include <type_traits>
#include <algorithm>



template<class T>
typename std::enable_if<!std::numeric_limits<T>::is_integer, bool>::type
    almost_equal(T x, T y, int ulp)
{
    // the machine epsilon has to be scaled to the magnitude of the values used
    // and multiplied by the desired precision in ULPs (units in the last place)
    return std::fabs(x-y) <= std::numeric_limits<T>::epsilon() * std::fabs(x+y) * ulp
        // unless the result is subnormal
        || std::fabs(x-y) < std::numeric_limits<T>::min();
}

int main()
{
    double d1 = 0.2;
    double d2 = 1 / std::sqrt(5) / std::sqrt(5);
    std::cout << std::fixed << std::setprecision(20) 
        << "d1=" << d1 << "\nd2=" << d2 << '\n';

    if(d1 == d2)
        std::cout << "d1 == d2\n";
    else
        std::cout << "d1 != d2\n";

    if(almost_equal(d1, d2, 2))
        std::cout << "d1 almost equals d2\n";
    else
        std::cout << "d1 does not almost equal d2\n";
}

0

Byłbym bardzo ostrożny wobec którejkolwiek z tych odpowiedzi, która obejmuje odejmowanie zmiennoprzecinkowe (np. Fabs (ab) <epsilon). Po pierwsze, liczby zmiennoprzecinkowe stają się bardziej rzadkie przy większych wielkościach i przy wystarczająco dużych wielkościach, w których odstępy są większe niż epsilon, równie dobrze możesz po prostu wykonać a == b. Po drugie, odejmowanie dwóch bardzo bliskich liczb zmiennoprzecinkowych (jak to zwykle bywa, biorąc pod uwagę, że szukasz prawie równości) jest dokładnie tym, jak dostajesz katastrofalne anulowanie .

Chociaż nie jest przenośny, myślę, że odpowiedź Grom najlepiej radzi sobie z unikaniem tych problemów.


1
+1 za dobrą informację. Nie widzę jednak, jak można zepsuć porównanie równości, zwiększając błąd względny; IMHO błąd staje się znaczący tylko w wyniku odejmowania, jednak jego rząd wielkości w stosunku do dwóch odejmowanych argumentów powinien nadal być wystarczająco wiarygodny, aby ocenić równość. Chyba że rozdzielczość musi być ogólnie wyższa, ale w takim przypadku jedynym rozwiązaniem jest przejście do reprezentacji zmiennoprzecinkowej z bardziej znaczącymi bitami w mantysie.
sehe

Odjęcie dwóch prawie równych liczb NIE prowadzi do katastrofalnego anulowania - w rzeczywistości nie wprowadza żadnego błędu (twierdzenie qv Sterbzena). Katastrofalny anulowanie nastąpi wcześniej, podczas obliczania ai bsiebie. Nie ma absolutnie żadnego problemu z użyciem odejmowania zmiennoprzecinkowego jako części rozmytego porównania (chociaż jak powiedzieli inni, bezwzględna wartość epsilon może, ale nie musi być odpowiednia dla danego przypadku użycia)
Sneftel

0

W oprogramowaniu numerycznym istnieją przypadki, w których chcesz sprawdzić, czy dwie liczby zmiennoprzecinkowe są dokładnie takie same. Wysłałem to na podobne pytanie

https://stackoverflow.com/a/10973098/1447411

Nie można więc powiedzieć, że „CompareDoubles1” jest ogólnie błędne.


Właściwie to bardzo solidne odniesienie do dobrej odpowiedzi, chociaż jest bardzo wyspecjalizowane, aby ograniczyć kogoś bez wiedzy naukowej lub analizy numerycznej (tj. LAPACK, BLAS) do niezrozumienia kompletności. Innymi słowy, zakłada się, że przeczytałeś coś takiego jak wprowadzenie do przepisów numerycznych lub Analiza numeryczna Burden & Faires.
mctylr

0

To zależy od tego, jak dokładne ma być porównanie. Jeśli chcesz porównać dokładnie ten sam numer, po prostu idź z ==. (Prawie nigdy nie chcesz tego robić, chyba że faktycznie chcesz dokładnie tego samego numeru.) Na dowolnej przyzwoitej platformie możesz również wykonać następujące czynności:

diff= a - b; return fabs(diff)<EPSILON;

ponieważ fabswydaje się być dość szybki. Przez dość szybko mam na myśli, że jest to trochę bitowe ORAZ, więc lepiej być szybkim.

A sztuczki na liczbach całkowitych do porównywania podwójnych i zmiennoprzecinkowych są fajne, ale zwykle utrudniają efektywną obsługę różnych potoków procesora. I na pewno nie jest to szybsze na niektórych architekturach w kolejności ze względu na użycie stosu jako tymczasowego obszaru przechowywania dla często używanych wartości. (Load-hit-store dla tych, którym zależy.)


0

Pod względem skali ilości:

Jeśli epsilonniewielki ułamek wielkości (tj. Wartość względna) w pewnym sensie fizycznym Ai Btypy są porównywalne w tym samym sensie, niż myślę, że następujące jest całkiem poprawne:

#include <limits>
#include <iomanip>
#include <iostream>

#include <cmath>
#include <cstdlib>
#include <cassert>

template< typename A, typename B >
inline
bool close_enough(A const & a, B const & b,
                  typename std::common_type< A, B >::type const & epsilon)
{
    using std::isless;
    assert(isless(0, epsilon)); // epsilon is a part of the whole quantity
    assert(isless(epsilon, 1));
    using std::abs;
    auto const delta = abs(a - b);
    auto const x = abs(a);
    auto const y = abs(b);
    // comparable generally and |a - b| < eps * (|a| + |b|) / 2
    return isless(epsilon * y, x) && isless(epsilon * x, y) && isless((delta + delta) / (x + y), epsilon);
}

int main()
{
    std::cout << std::boolalpha << close_enough(0.9, 1.0, 0.1) << std::endl;
    std::cout << std::boolalpha << close_enough(1.0, 1.1, 0.1) << std::endl;
    std::cout << std::boolalpha << close_enough(1.1,    1.2,    0.01) << std::endl;
    std::cout << std::boolalpha << close_enough(1.0001, 1.0002, 0.01) << std::endl;
    std::cout << std::boolalpha << close_enough(1.0, 0.01, 0.1) << std::endl;
    return EXIT_SUCCESS;
}

0

Używam tego kodu:

bool AlmostEqual(double v1, double v2)
    {
        return (std::fabs(v1 - v2) < std::fabs(std::min(v1, v2)) * std::numeric_limits<double>::epsilon());
    }

2
Nie po to epsilonjest.
Sneftel

1
Dlaczego nie? Możesz to wyjaśnić?
debiut

2
@debuti epsilonjest jedynie odległość między 1 a następnego numeru przedstawialnym po 1. w najlepszym razie, że kod jest po prostu staramy się, by sprawdzić, czy oba numery są dokładnie równe względem siebie, ale ponieważ nie-potęgi 2 są mnożone przez epsilon, to nawet nie robi tego poprawnie.
Sneftel

2
Aha, i std::fabs(std::min(v1, v2))jest niepoprawny - dla wejść ujemnych wybiera ten o większej wielkości.
Sneftel

0

Piszę to dla java, ale może uznasz to za przydatne. Używa długich zamiast podwójnych, ale dba o NaN, podnormalne itp.

public static boolean equal(double a, double b) {
    final long fm = 0xFFFFFFFFFFFFFL;       // fraction mask
    final long sm = 0x8000000000000000L;    // sign mask
    final long cm = 0x8000000000000L;       // most significant decimal bit mask
    long c = Double.doubleToLongBits(a), d = Double.doubleToLongBits(b);        
    int ea = (int) (c >> 52 & 2047), eb = (int) (d >> 52 & 2047);
    if (ea == 2047 && (c & fm) != 0 || eb == 2047 && (d & fm) != 0) return false;   // NaN 
    if (c == d) return true;                            // identical - fast check
    if (ea == 0 && eb == 0) return true;                // ±0 or subnormals
    if ((c & sm) != (d & sm)) return false;             // different signs
    if (abs(ea - eb) > 1) return false;                 // b > 2*a or a > 2*b
    d <<= 12; c <<= 12;
    if (ea < eb) c = c >> 1 | sm;
    else if (ea > eb) d = d >> 1 | sm;
    c -= d;
    return c < 65536 && c > -65536;     // don't use abs(), because:
    // There is a posibility c=0x8000000000000000 which cannot be converted to positive
}
public static boolean zero(double a) { return (Double.doubleToLongBits(a) >> 52 & 2047) < 3; }

Należy pamiętać, że po wielu operacjach zmiennoprzecinkowych liczba może się znacznie różnić od oczekiwanych. Nie ma kodu, aby to naprawić.


0

Co powiesz na to?

template<typename T>
bool FloatingPointEqual( T a, T b ) { return !(a < b) && !(b < a); }

Widziałem różne podejścia - ale nigdy tego nie widziałem, więc jestem ciekawy, czy mogę o nich komentować!


To nie działa dla 1.99999999 i 1.99999998
Mehdi

-1
/// testing whether two doubles are almost equal. We consider two doubles
/// equal if the difference is within the range [0, epsilon).
///
/// epsilon: a positive number (supposed to be small)
///
/// if either x or y is 0, then we are comparing the absolute difference to
/// epsilon.
/// if both x and y are non-zero, then we are comparing the relative difference
/// to epsilon.
bool almost_equal(double x, double y, double epsilon)
{
    double diff = x - y;
    if (x != 0 && y != 0){
        diff = diff/y; 
    }

    if (diff < epsilon && -1.0*diff < epsilon){
        return true;
    }
    return false;
}

Użyłem tej funkcji do mojego małego projektu i działa, ale zwróć uwagę na następujące:

Błąd podwójnej precyzji może cię zaskoczyć. Powiedzmy, że epsilon = 1.0e-6, wtedy 1.0 i 1.000001 NIE powinny być uważane za równe zgodnie z powyższym kodem, ale na moim komputerze funkcja uważa je za równe, ponieważ 1,000001 nie może być dokładnie przetłumaczone na format binarny, jest to prawdopodobnie 1.0000009xxx. Testuję to z 1.0 i 1.0000011 i tym razem otrzymuję oczekiwany wynik.


-1

To kolejne rozwiązanie z lambda:

#include <cmath>
#include <limits>

auto Compare = [](float a, float b, float epsilon = std::numeric_limits<float>::epsilon()){ return (std::fabs(a - b) <= epsilon); };

Jest to dokładnie to samo, co wiele innych odpowiedzi, z wyjątkiem tego, że jest to lambda i nie ma żadnego wyjaśnienia, więc nie ma to dużej wartości jako odpowiedź.
stijn

-2

Moja droga może być niepoprawna, ale przydatna

Konwertuj oba zmiennoprzecinkowe na ciągi, a następnie porównuj ciągi

bool IsFlaotEqual(float a, float b, int decimal)
{
    TCHAR form[50] = _T("");
    _stprintf(form, _T("%%.%df"), decimal);


    TCHAR a1[30] = _T(""), a2[30] = _T("");
    _stprintf(a1, form, a);
    _stprintf(a2, form, b);

    if( _tcscmp(a1, a2) == 0 )
        return true;

    return false;

}

można również wykonać nadpisywanie operatora


+1: hej, nie zamierzam robić z tym programowania, ale pomysł floatów w obie strony pojawił się kilka razy na blogu Bruce'a Dawsona (traktat?: D) na ten temat i jeśli jesteś uwięziony pokoju, a ktoś przykłada broń do twojej głowy i mówi: „hej, musisz porównać dwa pływaki z X znaczącymi liczbami, masz 5 minut, GO!” prawdopodobnie należy to rozważyć. ;)
shelleybutterfly

@shelleybutterfly Znowu pojawiło się pytanie o najbardziej efektywny sposób porównywania dwóch liczb zmiennoprzecinkowych.
Tommy Andersen

@TommyA lol, być może, ale założę się, że podróż w obie strony została odrzucona z powodów niezwiązanych z wydajnością. Chociaż moja intuicja jest taka, że ​​byłaby raczej nieefektywna w porównaniu z matematyką HW fp, ale mówi również, że algorytmy w oprogramowaniu fp raczej nie będą miały co najmniej dużej różnicy. Z niecierpliwością oczekuję, że przeprowadzona analiza wykazała, że ​​obawy dotyczące wydajności są w tym przypadku znaczące. Poza tym czasami nie zawsze optymalna może być cenną odpowiedzią, a ponieważ została oceniona - pomimo tego, że jest to ważna technika, o której wspominał nawet blog Dawsona na ten temat, więc pomyślałem, że zasługuje na pozytywną opinię.
shelleybutterfly

-2

Nie można porównać dwóch doublez ustalonym EPSILON. W zależności od wartości double,EPSILON jest różna.

Lepszym podwójnym porównaniem byłoby:

bool same(double a, double b)
{
  return std::nextafter(a, std::numeric_limits<double>::lowest()) <= b
    && std::nextafter(a, std::numeric_limits<double>::max()) >= b;
}

-2

W bardziej ogólny sposób:

template <typename T>
bool compareNumber(const T& a, const T& b) {
    return std::abs(a - b) < std::numeric_limits<T>::epsilon();
}

4
Ta metoda ma wiele słabych stron, na przykład jeśli liczby ai bsą już mniejsze niż epsilon()tam, różnica może być nadal znacząca. I odwrotnie, jeśli liczby są bardzo duże, to nawet kilka bitów błędu spowoduje, że porównanie się nie powiedzie, nawet jeśli chciałbyś, aby liczby były równe. Ta odpowiedź to dokładnie ten rodzaj „ogólnego” algorytmu porównania, którego chcesz uniknąć.
SirGuy,

-3

Dlaczego nie wykonać bitowego XOR? Dwie liczby zmiennoprzecinkowe są równe, jeśli odpowiadające im bity są równe. Myślę, że decyzja o umieszczeniu bitów wykładniczych przed mantysą została podjęta w celu przyspieszenia porównania dwóch pływaków. Myślę, że wiele odpowiedzi tutaj nie ma sensu w porównaniu epsilon. Wartość Epsilon zależy tylko od tego, z jaką dokładnością porównywane są liczby zmiennoprzecinkowe. Na przykład po wykonaniu arytmetyki z liczbami zmiennoprzecinkowymi otrzymujesz dwie liczby: 2.5642943554342 i 2.5642943554345. Nie są one równe, ale dla rozwiązania liczą się tylko 3 cyfry dziesiętne, więc są one równe: 2,564 i 2,564. W takim przypadku wybierasz epsilon równy 0,001. Porównanie Epsilon jest również możliwe w bitowym XOR. Popraw mnie, jeśli się mylę.


Nie dodawaj tej samej odpowiedzi do wielu pytań. Odpowiedz na najlepszy i oznacz resztę jako duplikaty. Zobacz meta.stackexchange.com/questions/104227/…
Bhargav Rao

Nie sądzę, aby „porównanie epsilon” było możliwe przy użyciu tylko ExOr (i jednego lub dwóch porównań), nawet ograniczone do znormalizowanych reprezentacji w tym samym formacie.
greybeard
Korzystając z naszej strony potwierdzasz, że przeczytałeś(-aś) i rozumiesz nasze zasady używania plików cookie i zasady ochrony prywatności.
Licensed under cc by-sa 3.0 with attribution required.