Schnellster Algorithmus zur Identifizierung des kleinsten und größten x, der die Doppelpräzisionsgleichung x + a == b true

Question

Im Rahmen der statischen Analyse bin ich daran interessiert, die Werte von x im Dann-Zweig der Bedingung zu bestimmen unten:

double x; 
x = …; 
if (x + a == b) 
{ 
    … 

a und b kann mit doppelter Genauigkeit Konstanten sein (verallgemeinernd auf beliebige Ausdrücke ist der einfachste Teil des Problems), und die Compiler IEEE 754 streng folgen kann angenommen werden, davon ausgegangen werden (FLT_EVAL_METHOD 0). Der Rundungsmodus zur Laufzeit kann als nächster-gerade angenommen werden.

Wenn die Berechnung mit Rationals billig wäre, wäre es einfach: Die Werte für x wären die Zahlen mit doppelter Genauigkeit im rationalen Intervall (b - a - 0,5 * ulp1 (b) ... b - a + 0,5 * ulp2 (b)). Die Grenzen sollten enthalten sein, wenn b gerade ist, ausgeschlossen, wenn b ungerade ist, und ulp1 und ulp2 sind zwei leicht unterschiedliche Definitionen von "ULP", die identisch genommen werden können, wenn es einem nichts ausmacht, ein wenig Genauigkeit bei Zweierpotenzen zu verlieren.

Leider kann die Berechnung mit Rationals teuer sein. Betrachten Sie, dass eine andere Möglichkeit darin besteht, jede der Grenzen durch Dichotomie in 64 Additionen mit doppelter Genauigkeit zu erhalten (jede Operation entscheidet über ein Bit des Ergebnisses). 128 Gleitkomma-Additionen, um die untere und obere Grenze zu erhalten, können durchaus schneller sein als jede mathematische Lösung.

Ich frage mich, ob es eine Möglichkeit gibt, über die "128 Gleitkomma-Ergänzungen" Idee zu verbessern. Eigentlich habe ich meine eigene Lösung mit Änderungen des Rundungsmodus und nextafter Anrufe, aber ich würde niemanden Stil zu krampfen und ihnen eine elegantere Lösung verpassen wollen als die, die ich derzeit habe. Ich bin mir auch nicht sicher, ob das zweimalige Ändern des Rundungsmodus tatsächlich billiger ist als 64 Gleitkomma-Additionen.

Answer 1

Du hast bereits eine schöne und elegante Lösung in Frage:

Wenn Computing mit rationals billig war, wäre es einfach sein: die Werte für x die Zahlen mit doppelter Genauigkeit in den rationalen enthalten sein würden Intervall (b - a - 0,5 * ulp1 (b) ... b - a + 0,5 * ulp2 (b)). Die Grenzen sollten eingeschlossen werden, wenn b gerade ist, ausgeschlossen, wenn b ungerade ist, und ulp1 und ulp2 sind zwei leicht unterschiedliche Definitionen von "ULP", die identisch genommen werden können, wenn es etwas dagegen macht, ein wenig Genauigkeit auf Potenzen von zu verlieren zwei.

Das Folgende ist eine halbwegs begründete Skizze einer Teillösung für das auf diesem Absatz basierende Problem. Hoffentlich bekomme ich eine Chance, es bald auszufüllen. Um eine echte Lösung zu finden, musst du mit Subnormalen, Nullen, NaNs und all den anderen lustigen Dingen umgehen. Ich gehe davon aus, dass a und b sind sagen, dass 1e-300 < |a| < 1e300 und 1e-300 < |b| < 1e300, so dass keine Verrücktheit an irgendeinem Punkt auftritt.

Kein Überlauf und Unterlauf, Sie können ulp1(b) von b - nextafter(b, -1.0/0.0) erhalten. Sie können ulp2(b) von nextafter(b, 1.0/0.0) - b erhalten.

Wenn b/2 <= a <= 2b, dann sagt Sterbenz Theorem, dass b - a ist genau. So wird die nächste double der unteren Grenze und (b - a) + ulp2/2 wird die nächste double der oberen Grenze sein.Probieren Sie diese Werte und die Werte unmittelbar davor und danach aus und wählen Sie das breiteste Intervall aus, das funktioniert.

Wenn b > 2a, b - a > b/2. Der berechnete Wert von b - a ist um höchstens einen halben ul. Eine ulp1 ist höchstens zwei ulp, wie auch eine ulp2, so dass das rationale Intervall, das Sie gaben, höchstens zwei ul breit ist. Finde heraus, welche der fünf am nächsten liegenden Werte zu b-a funktionieren.

Wenn a > 2b, ist ein ul b-a mindestens so groß wie ein ulp b; Wenn etwas funktioniert, wette ich, dass es unter den drei nächsten Werten zu b-a sein muss. Ich stelle mir den Fall vor, wo a und b unterschiedliche Zeichen haben, funktioniert ähnlich.

Ich schrieb einen kleinen Stapel von C++ - Code, der diese Idee implementiert. Es fehlte nicht an zufälligen Fuzz-Tests (in ein paar verschiedenen Bereichen), bevor mir das Warten langweilig wurde. Hier ist es:

void addeq_range(double a, double b, double &xlo, double &xhi) { 
    if (a != a) return; // empty interval 
    if (b != b) { 
    if (a-a != 0) { xlo = xhi = -a; return; } 
    else return; // empty interval 
    } 
    if (b-b != 0) { 
    // TODO: handle me. 
    } 

    // b is now guaranteed to be finite. 
    if (a-a != 0) return; // empty interval 

    if (b < 0) { 
    addeq_range(-a, -b, xlo, xhi); 
    xlo = -xlo; 
    xhi = -xhi; 
    return; 
    } 

    // b is now guaranteed to be zero or positive finite and a is finite. 
    if (a >= b/2 && a <= 2*b) { 
    double upulp = nextafter(b, 1.0/0.0) - b; 
    double downulp = b - nextafter(b, -1.0/0.0); 
    xlo = (b-a) - downulp/2; 
    xhi = (b-a) + upulp/2; 
    if (xlo + a == b) { 
     xlo = nextafter(xlo, -1.0/0.0); 
     if (xlo + a != b) xlo = nextafter(xlo, 1.0/0.0); 
    } else xlo = nextafter(xlo, 1.0/0.0); 
    if (xhi + a == b) { 
     xhi = nextafter(xhi, 1.0/0.0); 
     if (xhi + a != b) xhi = nextafter(xhi, -1.0/0.0); 
    } else xhi = nextafter(xhi, -1.0/0.0); 
    } else { 
    double xmid = b-a; 
    if (xmid + a < b) { 
     xhi = xlo = nextafter(xmid, 1.0/0.0); 
     if (xhi + a != b) xhi = xmid; 
    } else if (xmid + a == b) { 
     xlo = nextafter(xmid, -1.0/0.0); 
     xhi = nextafter(xmid, 1.0/0.0); 
     if (xlo + a != b) xlo = xmid; 
     if (xhi + a != b) xhi = xmid; 
    } else { 
     xlo = xhi = nextafter(xmid, -1.0/0.0); 
     if (xlo + a != b) xlo = xmid; 
    } 
    } 
}