Gleitkommavergleich überarbeitet

Dieses Thema ist bei StackOverflow schon oft aufgetaucht, aber ich glaube, das ist eine neue Einstellung. Ja ich habe gelesenBruce Dawsons Artikel undWas jeder Informatiker über Gleitkomma-Arithmetik wissen sollte undDiese schöne Antwort.

Soweit ich weiß, gibt es in einem typischen System vier grundlegende Probleme beim Vergleich von Gleitkommazahlen auf Gleichheit:

Gleitkommaberechnungen sind nicht genauOba-b ist "klein" hängt vom Maßstab aba undbOba-b ist "klein" hängt von der Art aba undb (z. B. float, double, long double)Fließkommazahlen haben typischerweise + -Infinity-, NaN- und denormalisierte Darstellungen, von denen jede eine naive Formulierung stören kann

Diese Antwort - aka. "Der Google-Ansatz" - scheint beliebt zu sein. Es behandelt alle kniffligen Fälle. Und es skaliert den Vergleich sehr genau und prüft, ob zwei Werte innerhalb einer festen Zahl von liegenULPs von einander. So vergleicht zum Beispiel eine sehr große Zahl "fast gleich" mit unendlich.

Jedoch:

Es ist meiner Meinung nach sehr chaotisch.Es ist nicht besonders portabel, stützt sich stark auf interne Darstellungen, verwendet eine Vereinigung, um die Bits von einem Float zu lesen usw.Es verarbeitet nur IEEE 754 mit einfacher und doppelter Genauigkeit (insbesondere kein x86-Long-Double).

Ich möchte etwas ähnliches, aber mit Standard-C ++ und lange Doppelte. Mit "Standard" meine ich wenn möglich C ++ 03 und wenn nötig C ++ 11.

Hier ist mein Versuch.

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

Ich behaupte, dass dieser Code (a) alle relevanten Fälle behandelt, (b) das Gleiche tut wie die Google-Implementierung für IEEE-754 mit einfacher und doppelter Genauigkeit und (c) vollkommen standardmäßiges C ++ ist.

Eine oder mehrere dieser Behauptungen sind mit ziemlicher Sicherheit falsch. Ich werde jede Antwort akzeptieren, die dies demonstriert, vorzugsweise mit einem Fix. Eine gute Antwort sollte eines oder mehrere der folgenden Elemente enthalten:

Spezifische Eingaben unterscheiden sich um mehr alsulps Einheiten an letzter Stelle, für die diese Funktion jedoch true zurückgibt (je größer der Unterschied, desto besser)Spezifische Eingaben unterscheiden sich um bis zuulps Einheiten an letzter Stelle, für die diese Funktion jedoch false zurückgibt (je kleiner der Unterschied, desto besser)Alle Fälle, die ich verpasst habeJede Art und Weise, in der dieser Code auf undefiniertem Verhalten beruht oder je nach implementierungsdefiniertem Verhalten abbricht. (Bitte geben Sie nach Möglichkeit eine entsprechende Spezifikation an.)Behebt alle Probleme, die Sie identifizierenEine Möglichkeit, den Code zu vereinfachen, ohne ihn zu beschädigen

Ich beabsichtige, dieser Frage eine nicht triviale Belohnung zu geben.

Antworten auf die Frage(3)

Ihre Antwort auf die Frage