Im Zusammenhang mit der statischen Analyse bin ich daran interessiert, die Werte von x im then-Zweig der folgenden Bedingung zu bestimmen:
a und b können als Konstanten mit doppelter Genauigkeit angenommen werden (die Verallgemeinerung zu beliebigen Ausdrücken ist der einfachste Teil des Problems), und der Compiler kann streng nach IEEE 754 folgen ( FLT_EVAL_METHOD ist 0) ). Der Rundungsmodus zur Laufzeit kann als nächster-gerade angenommen werden.
Wenn das Berechnen mit Rationals billig wäre, wäre es einfach: Die Werte für x wären die Zahlen mit doppelter Genauigkeit, die im rationalen Intervall enthalten sind (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 einem nichts ausmacht, ein wenig Genauigkeit bei Potenzen zu verlieren zwei.
Leider kann das Berechnen 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 Fließkomma-Additionen, um die untere und obere Grenze zu erhalten, sind möglicherweise schneller als jede mathematische Lösung.
Ich frage mich, ob es eine Möglichkeit gibt, die Idee "128 Gleitkomma-Additionen" zu verbessern. Eigentlich habe ich meine eigene Lösung mit Änderungen im Rundungsmodus und in nextafter Calls, aber ich würde nicht den Style von irgendjemandem verkrampfen wollen und sie dazu bringen, eine elegantere Lösung zu verpassen als die, die ich momentan habe. Ich bin mir auch nicht sicher, ob das zweimalige Ändern des Rundungsmodus tatsächlich billiger ist als 64 Gleitkomma-Additionen.
Sie haben in Ihrer Frage schon eine schöne und elegante Lösung gefunden:
Wenn Computing mit Rationals billig wäre, wäre es einfach: die Werte für x wären die Zahlen mit doppelter Genauigkeit, die im rationalen enthalten sind Intervall (b - a - 0,5 * ulp1 (b) ... b - a + 0,5 * ulp2 (b)). Die Grenzen sollte enthalten sein, wenn b gerade ist, ausgeschlossen, wenn b ungerade ist, und ulp1 und ulp2 sind zwei leicht unterschiedliche Definitionen von "ULP", die genommen werden können identisch, wenn es einem nichts ausmacht, ein wenig Präzision bei den Potenzen zu verlieren zwei.
Was folgt, ist eine halbwegs begründete Skizze einer Teillösung des Problems, die auf diesem Absatz basiert. 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 sagen, dass 1e-300 < |a| < 1e300 und 1e-300 < |b| < 1e300 sind, so dass an keiner Stelle Verrücktheit 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 der Satz von Sterbenz, dass b - a genau ist. So wird (b - a) - ulp1 / 2 der unteren Grenze am nächsten double und (b - a) + ulp2 / 2 der oberen Grenze am nächsten double 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 höchstens um eine halbe ul. Entfernt. Ein ulp1 ist höchstens zwei ulp, wie auch ein ulp2 , also ist das von Ihnen angegebene rationale Intervall höchstens zwei ul breit. Finde heraus, welche der fünf am nächsten kommenden Werte für b-a funktionieren.
Wenn a > 2b , ein ulp von b-a ist mindestens so groß wie ein ulp von b ; Wenn etwas funktioniert, wette ich, dass es zu den drei engsten Werten von b-a gehören muss. Ich stelle mir vor, dass der Fall, in dem a und b unterschiedliche Vorzeichen haben, ähnlich funktioniert.
Ich habe einen kleinen Stapel C ++ - Code geschrieben, 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:
%Vor%Tags und Links c floating-point ieee-754