Ich passe momentan Fließkommazahlen an, um einen Wert abzuschätzen. (Es ist: p(k,t) für diejenigen, die interessiert sind.) Im Grunde kann das Dienstprogramm nie ergibt eine Unterschätzung dieses Wertes: Die Sicherheit der wahrscheinlichen Primegeneration hängt von einer numerisch robusten Implementierung ab. Während die Ausgabeergebnisse mit den veröffentlichten Werten übereinstimmen, habe ich den DBL_EPSILON -Wert verwendet, um sicherzustellen, dass insbesondere die Division ein Ergebnis liefert, das nie kleiner als der wahre Wert ist:
Betrachte: double x, y; /* assigned some values... */
Die Auswertung: r = x / y; tritt häufig auf, aber diese (endliche Genauigkeit) Ergebnisse können signifikante Stellen vom wahren Ergebnis abschneiden - eine möglicherweise unendliche rationale Erweiterung. Ich versuche derzeit, dies zu mildern, indem ich eine Vorspannung auf den Zähler anlege, d. H.,
Wenn Sie etwas über dieses Thema wissen, ist p(k,t) normalerweise viel kleiner als die meisten Schätzungen - aber es ist einfach nicht gut genug, um das Problem mit dieser "Beobachtung" zu verwerfen. Ich kann natürlich sagen:
Natürlich muss ich sicherstellen, dass das "voreingenommene" Ergebnis größer oder gleich dem "exakten" Wert ist. Obwohl ich sicher bin, dass es mit der Manipulation oder Skalierung von DBL_EPSILON zu tun hat, möchte ich natürlich, dass das "verzerrte" Ergebnis das "genaue" Ergebnis durch ein Minimum übersteigt - nachweisbar unter IEEE-754 arithmetischen Annahmen.
Ja, ich habe nach Goldbergs Papier gesucht und nach einer robusten Lösung gesucht. Bitte schlagen Sie nicht die Manipulation von Rundungsmodi vor. Idealerweise bin ich auf der Suche nach jemandem mit einem sehr guten Verständnis von Fließkomma-Theoremen oder kenne ein sehr gut illustriertes Beispiel.
>
BEARBEITEN: Zur Klärung, (((1.0 + DBL_EPSILON) * x) / y) oder ein Formular (((1.0 + c) * x) / y) , ist keine Voraussetzung. Das war einfach ein Ansatz, den ich als "wahrscheinlich gut genug" benutzte, ohne dafür eine solide Grundlage geschaffen zu haben. I kann angeben, dass der Zähler und der Nenner keine speziellen Werte sind: NaNs, Infs usw., und der Nenner nicht Null sein wird.
Erstens: Ich weiß, dass Sie den Rundungsmodus nicht einstellen wollen, aber das sollte wirklich gesagt werden
In Bezug auf die Präzision, wie andere bemerkt haben, wird die Einstellung des Rundungsmodus so gut wie möglich zu einer Antwort führen. Konkret vorausgesetzt, dass x und y beide positiv sind (was der Fall zu sein scheint, aber in Ihrer Frage nicht explizit angegeben wurde), ist das Folgende ein Standard-C-Snippet mit dem gewünschten Effekt [1]:
Nun, abgesehen davon, gibt es legitime Gründe, den Rundungsmodus nicht ändern zu wollen (einige Plattformen unterstützen kein Round-to-Plus-Infinity, auf manchen Plattformen ändert sich der Rundungsmodus, führt zu einem großen Serialisierungs-Stall usw.) ), und Ihr Wunsch, dies nicht zu tun, sollte nicht so beiläufig beiseite geschoben werden. Also, respektieren Sie Ihre Frage, was können wir noch tun?
Wenn Ihre Plattform Fused Multiply-Add unterstützt, steht Ihnen eine sehr elegante Lösung zur Verfügung:
%Vor%Auf Plattformen mit Hardware-FMA-Unterstützung ist dies sehr effizient. Auch wenn fma () in Software implementiert ist, kann es akzeptabel sein. Dieser Ansatz hat den Vorteil, dass er das gleiche Ergebnis liefert wie der Rundungsmodus. das heißt, die engste mögliche Grenze.
Wenn die C-Bibliothek Ihrer Plattform antediluvian ist und fma nicht zur Verfügung stellt, besteht immer noch Hoffnung. Ihre Behauptung ist korrekt (unter der Annahme, dass es keine Denormalwerte gibt - ich müsste mehr darüber nachdenken, was bei Denormalen passiert); (1.0+DBL_EPSILON)*x/y ist wirklich immer größer oder gleich dem unendlich genauen x / y. Es ist manchmal ein ul größer als der kleinste Wert mit dieser Eigenschaft, aber das ist eine sehr kleine und wahrscheinlich akzeptable Marge. Der Beweis für diese Behauptungen ist ziemlich pingelig und wahrscheinlich nicht für StackOverflow geeignet, aber ich gebe eine kurze Skizze:
(1,0 + eps) * x & gt; = x + eps & gt; x. Um dies zu sehen, beachten Sie:
%Vor%Sei P das mathematisch genaue x / y. Wir haben:
%Vor%Nun, y ist oben durch 2 begrenzt, also gibt uns das:
%Vor% was ausreichend ist, um zu garantieren, dass das Ergebnis auf einen Wert & gt; = P rundet. Dies zeigt uns auch den Weg zu einer engeren Grenze. Wir könnten stattdessen nextafter(x,INFINITY)/y verwenden, um in vielen Fällen den gewünschten Effekt mit einer engeren Grenze zu erhalten. ( nextafter(x,INFINITY) ist immer x + ulp, während (1.0 + eps)*x die Hälfte der Zeit x + 2ulp ist. Wenn Sie die nextafter Bibliotheksfunktion nicht aufrufen möchten, können Sie stattdessen (x + (0.75*DBL_EPSILON)*x) verwenden, um das gleiche Ergebnis zu erhalten. unter der Arbeitsannahme positiver Normalwerte).
Um wirklich pedantisch korrekt zu sein, würde dies erheblich komplizierter werden. Niemand schreibt wirklich Code wie diesen, aber es würde in diese Richtung gehen:
%Vor%Tags und Links c floating-point ieee-754 floating-accuracy