Ich war überrascht, dass es schwierig ist, Informationen zu diesem Problem auf Russisch zu finden, als ob es nur wenige Menschen interessiert, dass mathematische Bibliotheken, die in modernen Compilern verwendet werden, manchmal kein korrekt gerundetes Ergebnis liefern. Ich mache mir Sorgen um diese Situation, da ich gerade an der Entwicklung solcher mathematischer Bibliotheken arbeite. In der ausländischen Literatur ist dieses Problem gut abgedeckt, daher habe ich beschlossen, es in populärwissenschaftlicher Form auf Russisch zu präsentieren, wobei ich mich auf westliche Quellen und noch ein wenig persönliche Erfahrung stütze.

---

Freunde, für Ihre Bequemlichkeit ist der Artikel auch im Videopräsentationsformat (fast 34 Minuten) erhältlich. Dieses Format eignet sich eher für Leser, die Schwierigkeiten haben, die erforderlichen mathematischen Bilder in ihren Köpfen zu erstellen, da die Präsentation viel illustratives Material enthält. Die Informationen im Video sind vollständig identisch mit dem Inhalt des Artikels. Bitte handeln Sie nach Belieben.

---

Ich wiederhole, dass dies kein wissenschaftlicher, sondern ein populärwissenschaftlicher Artikel ist. Nach dem Lesen werden Sie sich kurz damit vertraut machen.

• Transzendentale Elementarfunktionen (exp, sin, log, cosh und andere), die mit Gleitkomma-Arithmetik arbeiten, werden falsch gerundet, manchmal machen sie im letzten Bit einen Fehler.
• Der Grund für Fehler liegt nicht immer in der Faulheit oder geringen Qualifikation der Entwickler, sondern in einem grundlegenden Umstand, den die moderne Wissenschaft noch nicht überwinden konnte.
• «», - .
• , , , , exp2(x) pow(2.0, x).

Um diesen Artikel zu verstehen, müssen Sie mit dem Gleitkommaformat IEEE-754 vertraut sein. Es reicht aus, wenn Sie nur verstehen, dass dies beispielsweise Folgendes ist: 0x400921FB54442D18 - Nummer pi im Format mit doppelter Genauigkeit (binär64 oder doppelt), das heißt, Sie verstehen nur, was ich mit diesem Datensatz meine; Ich muss nicht in der Lage sein, solche Transformationen im laufenden Betrieb durchzuführen. Und ich werde Sie an die Rundungsmodi in diesem Artikel erinnern, dies ist ein wichtiger Teil der Geschichte. Es ist auch wünschenswert, "Programmierer" Englisch zu sprechen, da es Begriffe und Zitate aus der westlichen Literatur gibt, aber Sie können mit einem Online-Übersetzer auskommen.



Beispiele zuerst, damit Sie sofort verstehen, worum es im Gespräch geht. Jetzt werde ich den Code in C ++ geben, aber wenn dies nicht Ihre Sprache ist, dann werden Sie sicher immer noch leicht verstehen, was geschrieben steht. Bitte schauen Sie sich diesen Code an:

#include <stdio.h>
#include <cmath>

int main() {
  float x = 0.00296957581304013729095458984375f;  // ,  .
  float z;
  z = exp2f(x);  // z = 2**x  .
  printf ("%.8f\n", z);  //      8   .
  z = powf(2.0f, x);  // z = 2**x  
  printf ("%.8f\n", z);  //   .
  return 0;
}

Die Zahl x wird absichtlich mit einer solchen Anzahl von signifikanten Stellen geschrieben, dass sie im Float-Typ genau darstellbar ist, dh dass der Compiler sie ohne Rundung in einen Binärcode konvertiert. Schließlich wissen Sie sehr gut, dass einige Compiler nicht fehlerfrei runden können (wenn Sie es nicht wissen, geben Sie in den Kommentaren an, dass ich einen separaten Artikel mit Beispielen schreiben werde). Als nächstes müssen wir im Programm 2 ^x berechnen , aber lassen Sie uns dies auf zwei Arten tun: die Funktion exp2f (x) und die explizite Potenzierung von zwei powf (2.0f, x). Das Ergebnis wird natürlich anders sein, da ich oben gesagt habe, dass Elementarfunktionen nicht in allen Fällen korrekt funktionieren können, und ich habe speziell ein Beispiel ausgewählt, um dies zu zeigen. Hier ist die Ausgabe:

1.00206053
1.00206041

Vier Compiler gaben mir diese Werte: Microsoft C ++ (19.00.23026), Intel C ++ 15.0, GCC (6.3.0) und Clang (3.7.0). Sie unterscheiden sich in einem niedrigstwertigen Bit. Hier ist der Hexadezimalcode für diese Zahlen:

0x3F804385  // 
0x3F804384  // 

Bitte denken Sie an dieses Beispiel. Darauf werden wir uns etwas später mit der Essenz des Problems befassen. Um jedoch einen klareren Eindruck zu erhalten, lesen Sie zunächst die Beispiele für den Datentyp mit doppelter Genauigkeit (double, binary64) mit einigen anderen Elementarfunktionen. Ich präsentiere die Ergebnisse in der Tabelle. Richtige Antworten (falls verfügbar) haben * am Ende.







Ich hoffe, Sie haben nicht den Eindruck, dass ich absichtlich einige völlig einzigartige Tests gemacht habe, die Sie kaum finden können? Wenn ja, lassen Sie uns auf unseren Knien eine vollständige Aufzählung aller möglichen Bruchargumente für die 2 ^x -Funktion für den Float-Datentyp kochen . Es ist klar, dass wir nur an x-Werten zwischen 0 und 1 interessiert sind, da andere Argumente ein Ergebnis erzeugen, das sich nur im Wert im Exponentenfeld unterscheidet und nicht von Interesse ist. Sie selbst verstehen:

$$

$$2^x = 2^{[x]}\cdot2^{\{x\}}.$$

Nachdem ich ein solches Programm geschrieben hatte (der versteckte Text wird unten sein), überprüfte ich die exp2f-Funktion und wie viele fehlerhafte Werte sie im Intervall x von 0 bis 1 erzeugt.







Aus dem folgenden Programm geht hervor, dass die Anzahl der getesteten Argumente x 197612997 betrug. Es stellt sich heraus, dass beispielsweise Microsoft C ++ die 2 ^x -Funktion für fast ein Prozent von ihnen falsch berechnet . Freut euch nicht, liebe Fans von GCC und Clang, es ist nur so, dass diese Funktion in diesen Compilern korrekt implementiert ist, aber in anderen voller Fehler.





Ich werde den Leser mit diesen Beispielen nicht langweilen. Die Hauptsache hier war zu zeigen, dass moderne Implementierungen von transzendentalen Funktionen das letzte Bit falsch abrunden können und verschiedene Compiler an verschiedenen Stellen Fehler machen, aber keiner von ihnen wird richtig funktionieren. Übrigens erlaubt der IEEE-754-Standard diesen Fehler im letzten Bit (worüber ich weiter unten sprechen werde), aber es scheint mir immer noch seltsam: ok double, dies ist ein großer Datentyp, aber float kann mit roher Gewalt überprüft werden! War es so schwer zu tun? Gar nicht so schwer, und ich habe bereits ein Beispiel gezeigt.



Unser Aufzählungscode enthält eine "selbstgeschriebene" Funktion zur korrekten Berechnung 2 ^xunter Verwendung eines Approximationspolynoms 10. Grades, und es wurde in wenigen Minuten geschrieben, da solche Polynome automatisch abgeleitet werden, beispielsweise im Maple-Computeralgebrasystem. Es reicht aus, eine Bedingung für das Polynom festzulegen, um eine Genauigkeit von 54 Bit bereitzustellen (für diese Funktion 2 ^x ). Warum 54? Aber Sie werden es bald herausfinden, gleich nachdem ich Ihnen die Essenz des Problems erklärt habe und warum es im Prinzip jetzt unmöglich ist, schnelle und korrekte transzendentale Funktionen für den Datentyp der Vierfachgenauigkeit (binär128) zu erstellen, obwohl es bereits Versuche gibt, dieses Problem theoretisch anzugreifen.

Standardrundung und das Problem damit

Wenn Sie nicht in die Entwicklung mathematischer Bibliotheken vertieft sind, ist es nichts Falsches, die Standardrundungsregel für Gleitkommazahlen gemäß dem IEEE-754-Standard zu vergessen. Deshalb werde ich Sie daran erinnern. Wenn Sie sich gut an alles erinnern, schauen Sie sich trotzdem mindestens das Ende dieses Abschnitts an, Sie werden überrascht sein: Ich werde Ihnen eine Situation zeigen, in der das Aufrunden einer Zahl sehr schwierig sein kann.



Sie können sich leicht daran erinnern, was unter dem Namen "aufrunden" (auf plus unendlich), "abrunden" (auf minus unendlich) oder "auf null runden" ist (wenn überhaupt, gibt es Wikipedia)). Die Hauptschwierigkeiten für Programmierer ergeben sich beim Runden "auf den nächsten, aber bei gleichem Abstand vom nächsten - auf den mit der letzten geraden Ziffer". Ja, so wird dieser Rundungsmodus übersetzt, den die westliche Literatur kurz nennt: "Runde am nächsten an gerade".



Dieser Rundungsmodus wird standardmäßig verwendet und funktioniert wie folgt. Wenn sich als Ergebnis von Berechnungen herausstellt, dass die Länge der Mantisse größer ist, als der resultierende Datentyp aufnehmen kann, wird auf den nächsten von zwei möglichen Werten gerundet. Es kann jedoch vorkommen, dass sich herausstellt, dass die ursprüngliche Zahl genau in der Mitte zwischen den beiden nächsten liegt. Dann wird das Ergebnis ausgewählt, für das sich herausstellt, dass das letzte Bit (nach dem Runden) gerade ist, dh gleich Null. Betrachten Sie vier Beispiele, bei denen Sie nach dem binären Dezimalpunkt auf zwei Bits runden müssen:

1. Runde 1.00 1 001. Das dritte Bit nach dem Dezimalpunkt ist 1, aber dann gibt es ein weiteres 6. Bit, das 1 ist, was bedeutet, dass die Rundung erhöht wird, da die ursprüngliche Zahl näher an 1.01 als an 1.00 liegt.
2. 1,001000. , 1,00 1,01, .
3. 1,011000. 1,01 1,10. , .
4. 1,010111. , 1,01, 1,10.

Aus diesen Beispielen geht hervor, dass alles einfach ist, aber nicht. Tatsache ist, dass wir manchmal nicht sicher sagen können, ob wir wirklich in der Mitte zwischen zwei Werten liegen. Siehe ein Beispiel. Lassen Sie uns noch einmal auf zwei Bits nach dem Dezimalpunkt runden:



1.00 1 000000000000000000000000000000000000001



Es ist Ihnen jetzt klar, dass die Rundung auf die Zahl 1.01 erfolgen sollte. Sie sehen jedoch eine Zahl mit 40 Bit nach dem Dezimalpunkt. Was wäre, wenn Ihr Algorithmus keine Genauigkeit von 40 Bit liefern könnte und nur 30 Bit erreicht? Dann wird eine andere Nummer ausgegeben :



1.00 1 000000000000000000000000000



Da Sie nicht wissen, dass es an der 40. Position (die der Algorithmus nicht berechnen kann) eine geschätzte gibt, runden Sie diese Zahl ab und erhalten 1,00, was falsch ist. Sie haben das letzte Stück falsch aufgerundet - das ist das Thema unserer Diskussion. Aus dem Obigen geht hervor, dass Sie die Funktion bis zu 40 Bit berechnen müssen, um nur das 2. Bit korrekt zu erhalten! Wow! Und wenn sich herausstellt, dass die "Lokomotive" der Nullen noch länger ist? Darüber werden wir im nächsten Abschnitt sprechen.



Dies ist übrigens der Fehler, den viele Compiler machen, wenn sie die Dezimalschreibweise einer Gleitkommazahl in das resultierende Binärformat konvertieren. Wenn die ursprüngliche Dezimalzahl im Programmcode zu nahe an der Mitte zwischen zwei genau darstellbaren Binärwerten liegt, wird sie nicht richtig gerundet. Dies ist jedoch nicht das Thema dieses Artikels, sondern ein Grund für eine separate Geschichte.

Die Essenz des Problems der Rundung des letzten signifikanten Bits

Das Problem tritt aus zwei Gründen auf. Das erste ist die bewusste Ablehnung zeitaufwändiger Berechnungen zugunsten der Geschwindigkeit. In diesem Fall ist es eine Nebensache, solange die angegebene Genauigkeit eingehalten wird und welche Bits in der Antwort enthalten sind. Der zweite Grund ist das Dilemma des Tischmachers, das das Hauptthema unseres Gesprächs ist. Lassen Sie uns beide Gründe genauer betrachten.

Erster Grund

Sie verstehen natürlich, dass die Berechnung transzendentaler Funktionen durch einige Näherungsmethoden implementiert wird, beispielsweise durch die Methode der Approximation von Polynomen oder sogar (selten) durch Reihenexpansion. Um die Berechnungen so schnell wie möglich durchzuführen, verpflichten sich die Entwickler, so wenige Iterationen der numerischen Methode wie möglich durchzuführen (oder ein Polynom mit dem geringstmöglichen Grad zu verwenden), sofern der Algorithmus einen Fehler zulässt, der den halben Wert des letzten Bits der Mantisse nicht überschreitet. In der Literatur wird dies als 0,5 ulp geschrieben (ulp = Einheit an letzter Stelle ).



Wenn es sich beispielsweise um eine Zahl x vom Typ float im Intervall (0,5; 1) handelt, ist der Wert ulp = 2 ^-23 . Im Intervall (1; 2) ist ulp = 2 ^-22 . Mit anderen Worten, wenn x im Intervall (0; 1) liegt, dann 2 ^xwird auf dem Intervall (1,2), und eine Genauigkeit von 0.5ulp zu gewährleisten, müssen Sie, grob gesagt, auf EPS = 2 wählen ^-23 (so werden wir die Konstante „Epsilon“, in dem gemeinen Volk als „Fehler“ oder „Genauigkeit“ bezeichnen, die wie Sie möchten, finden Sie bitte keinen Fehler).



Für angewandte Berechnungen ist dies ausreichend, aber die Tatsache, dass die letzten Bits möglicherweise nicht mit dem absoluten Ergebnis übereinstimmen, ist für fast 100% der Programmierer nicht wichtig, da es für sie nicht wichtig ist, wie die Bits sein werden, sondern wie genau sie sein werden.



Für diejenigen, die nicht verstehen, werde ich ein Beispiel im Dezimalzahlensystem geben. Hier sind zwei Zahlen: 1.999999 und 2.0. Nehmen wir an, das erste ist das, was der Programmierer erhalten hat, und das zweite ist der Standard dessen, was hätte passieren sollen, wenn wir unbegrenzte Möglichkeiten hätten. Der Unterschied zwischen ihnen beträgt nur ein Millionstel, dh die Antwort wird mit einem Fehler von EPS = 10 ^-6 berechnet . Diese Antwort enthält jedoch keine einzige korrekte Nummer. Ist es schlimm? Nein, aus Sicht des Anwendungsprogramms ist dies lila, der Programmierer rundet die Antwort beispielsweise auf zwei Dezimalstellen und erhält 2,00 (zum Beispiel ging es um Währung, 2,00 USD), er braucht nicht mehr, aber die Tatsache, dass er Setzen Sie EPS = 10 ^-6 in mein Programm , dann gut gemacht, haben Sie eine Marge für den Fehler der Zwischenberechnung genommen und das Problem richtig gelöst.



Mit anderen Worten, seien Sie nicht verwirrt: Die Genauigkeit und die Anzahl der richtigen Bits (oder Ziffern) sind zwei verschiedene Dinge. Diejenigen, die Genauigkeit benötigen (dies sind fast 100% der Programmierer), betrifft das diskutierte Problem überhaupt nicht. Jeder, der eine Bitfolge benötigt, um mit einer korrekt gerundeten Referenz übereinzustimmen, ist über dieses Problem sehr besorgt, beispielsweise Entwickler von Bibliotheken elementarer Funktionen. Trotzdem ist es für alle nützlich, dies für die allgemeine Entwicklung zu wissen.



Ich möchte Sie daran erinnern, dass dies die erste Richtung des Problems war: Die letzten Teile der Antwort können falsch sein, da dies eine absichtliche Lösung ist. Die Hauptsache ist, die Genauigkeit von 0,5 ulp (oder höher) beizubehalten. Daher wird der numerische Algorithmus nur dann aus dieser Bedingung ausgewählt, wenn er nur extrem schnell funktioniert. In diesem Fall erlaubt der Standard die Implementierung von Elementarfunktionen ohne korrekte Rundung des letzten Bits. Ich zitiere [1, Abschnitt 12.1] (Englisch):

In der Version 1985 des IEEE 754-Standards für Gleitkomma-Arithmetik wurde nichts bezüglich der Elementarfunktion angegeben. Dies lag daran, dass seit Jahren angenommen wird, dass korrekt gerundete Funktionen zumindest für einige Eingabeargumente viel zu langsam sind. Die Situation hat sich seitdem geändert und die Version 2008 des Standards empfiehlt (erfordert jedoch nicht), dass einige Funktionen korrekt gerundet werden.

Die folgenden Funktionen werden empfohlen, müssen jedoch nicht richtig gerundet werden:

Der zweite Grund

Schließlich kamen wir zum Gesprächsthema: Table Maker's Dilemma (abgekürzt als TMD). Ich konnte seinen Namen nicht angemessen ins Russische übersetzen, er wurde von William Kahan (Gründungsvater von IEEE-754) in Artikel [2] eingeführt. Wenn Sie den Artikel lesen, werden Sie vielleicht verstehen, warum der Name genau so ist. Kurz gesagt, das Wesentliche des Dilemmas ist, dass wir eine absolut genaue Rundung der Funktion z = f (x) erhalten müssen, als ob uns eine unendliche Bitaufzeichnung des perfekt berechneten Ergebnisses z zur Verfügung stünde. Aber es ist jedem klar, dass wir keine unendliche Folge bekommen können. Wie viele Bits müssen dann genommen werden? Oben habe ich ein Beispiel gezeigt, in dem 40 Bits des Ergebnisses angezeigt werden müssen, um nach dem Runden mindestens 2 korrekte Bits zu erhalten. Und das Wesentliche des TMD-Problems ist, dass wir es nicht im Voraus wissenbis zu wie viele Bits, um den Wert von z zu berechnen, damit nach dem Runden so viele Bits korrekt sind, wie wir benötigen. Was ist, wenn es hundert oder tausend gibt? Wir wissen es nicht im Voraus!



Wie ich bereits sagte, müssen wir für die Funktion 2 ^x für den Datentyp float, bei dem der Bruchteil der Mantisse nur 23 Bit hat, die Berechnung mit einer Genauigkeit von 2 bis ^{54 durchführen,} damit die Rundung für alle möglichen x-Argumente ausnahmslos korrekt erfolgt. Es ist nicht schwierig, diese Schätzung durch umfassende Suche zu erhalten, aber für die meisten anderen Funktionen, insbesondere für Typen doppelt oder lang doppelt (setzen Sie "Klasse", wenn Sie wissen, was es ist), sind solche Schätzungen unbekannt .



Lassen Sie uns bereits verstehen, warum dies geschieht. Ich habe absichtlich das allererste Beispiel in diesem Artikel mit dem Datentyp float angegeben und Sie gebeten, sich daran zu erinnern, da es in diesem Typ nur 32 Bit gibt und es einfacher ist, es zu betrachten, in anderen Datentypen ist die Situation ähnlich.



Wir haben mit der Zahl x = 0,00296957581304013729095458984375 begonnen. Dies ist eine genau darstellbare Zahl im Float-Datentyp, dh sie ist so geschrieben, dass sie ohne Rundung in das binäre Float-System konvertiert werden kann. Wir berechnen 2 ^x , und wenn wir einen Taschenrechner mit unendlicher Genauigkeit hatten, sollten wir erhalten (Sie können mich überprüfen, die Berechnung erfolgt im Online-System WolframAlpha ):



1.0020604729652405753669743044108123031635398201893943954577320057 ...



Lassen Sie uns diese Zahl in eine Binärzahl übersetzen.



Angenommen , 64 Bit reichen aus: 1.00000000100001110000100 1 000000000000000000000000000001101111101



Das Rundungsbit (24. Bit nach dem Dezimalpunkt) ist unterstrichen. Frage: Wo runden? Hoch oder runter? Natürlich wissen Sie das, weil Sie genug Teile sehen und eine Entscheidung treffen können. Aber schau genau hin ...



Nach dem Rundungsbit haben wir 29 Nullen. Dies bedeutet, dass wir uns sehr, sehr nahe an der Mitte zwischen den beiden nächsten Zahlen befinden und es ausreicht, sich nur ein wenig nach unten zu bewegen, da sich die Rundungsrichtung ändert. Aber die Frage ist: Wo wird diese Verschiebung sein? Der numerische Algorithmus kann sich nacheinander Schritt für Schritt dem genauen Wert von verschiedenen Seiten nähern, und bis wir alle diese 29 Nullen passieren und eine Genauigkeit erreichen, die den Wert der allerletzten Null in dieser "Lokomotive" überschreitet, kennen wir die Rundungsrichtung nicht ... Was ist, wenn in Wirklichkeit die richtige Antwort lauten sollte:



1.00000000100001110000100 0 11111111111111111111111111111?



Dann wird die Rundung nach unten sein.



Wir wissen das erst, wenn unsere Genauigkeit das 54. Bit nach dem Dezimalpunkt erreicht. Wenn das 54. Bit genau bekannt ist, wissen wir genau, welcher der beiden nächsten Zahlen wir tatsächlich näher sind. Solche Zahlen werden als am härtesten zu-rund-Punkte [1, Abschnitt 12.3] (kritische Punkte für die Rundung), und die Zahl 54 bezeichnet Härte-zu-rund, und wird durch den Buchstaben M bezeichnet in dem Buch zitiert.



Die Komplexität der Rundung (m) ist die Anzahl der Bits, die das Minimum ist, um sicherzustellen, dass für alle Argumente einer bestimmten Funktion f (x) und für einen vorgewählten Bereich die Funktion f (x) korrekt auf das letzte Bit gerundet wird (für verschiedene Rundungsmodi kann es unterschiedliche geben Wert m). Mit anderen Worten, für den Datentyp float und für das x-Argument aus dem Bereich (0; 1) für den Rundungsmodus "nächste gerade" beträgt die Rundungskomplexität m = 54. Dies bedeutet, dass wir für absolut alle x aus dem Intervall (0; 1) die gleiche Genauigkeit ESP = 2 ^-54 in den Algorithmus ^einfügen können und alle Ergebnisse nach dem binären Dezimalpunkt korrekt auf 23 Bit gerundet werden.



In der Tat sind einige Algorithmen in der Lage, ein genaues Ergebnis zu liefern, und basierend auf 53 und sogar 52 Bit zeigt Brute Force dies, aber theoretisch benötigen Sie genau 54. Wenn es nicht die Möglichkeit gäbe, die Brute Force herauszudrehen, könnten wir nicht "schummeln". und speichere ein paar Bits, wie ich es im obigen Brute-Force-Programm getan habe. Ich habe ein Polynom mit einem Grad genommen, der niedriger ist als es sollte, aber es funktioniert immer noch, nur weil ich Glück hatte.



Unabhängig vom Rundungsmodus gibt es also zwei Situationen: Entweder gibt es im Rundungsbereich eine "Lokomotive" mit Nullen oder eine "Lokomotive" mit Einsen. Die Aufgabe des richtigen Algorithmus zur Berechnung der transzendentalen Funktion f (x) besteht darin, den Wert dieser Funktion zu verfeinern, bis die Genauigkeit den Wert des letzten Bits dieser "Lokomotive" überschreitet, und bis genau klar wird, dass infolge nachfolgender Schwankungen des numerischen Algorithmus zur Berechnung von f (x) Nullen werden nicht zu Einsen oder umgekehrt. Sobald sich alles stabilisiert hat und der Algorithmus eine Genauigkeit erreicht hat, die über die Grenzen der "Dampflokomotive" hinausgeht, können wir runden, als hätten wir eine unendliche Anzahl von Bits. Und diese Rundung erfolgt mit dem richtigen letzten Bit. Aber wie kann das erreicht werden?

"Krücken"

Wie erwähnt, besteht das Hauptproblem darin, den Algorithmus dazu zu bringen, die Lokomotive von Nullen oder Einsen zu überwinden, die unmittelbar nach dem Rundungsbit kommt. Wenn die Lokomotive überwunden ist und wir sie als Ganzes sehen, entspricht dies der Tatsache, dass diese Nullen oder Einsen bereits genau berechnet wurden und wir bereits genau wissen , in welche Richtung die Rundung jetzt erfolgen wird. Aber wenn wir die Länge der Lokomotive nicht kennen, wie können wir dann einen Algorithmus entwerfen?

Die erste "Krücke"

Dem Leser mag es so erscheinen, als ob die Antwort offensichtlich ist: Nehmen Sie die Arithmetik mit unendlicher Präzision und geben Sie eine absichtlich übermäßige Anzahl von Bits ein. Wenn dies nicht ausreicht, geben Sie eine andere ein und berechnen Sie sie neu. Im Allgemeinen ist es richtig. Dies geschieht, wenn die Geschwindigkeit und die Ressourcen des Computers keine besondere Rolle spielen. Dieser Ansatz hat einen Namen: Zivs Mehrebenenstrategie [1, Abschnitt 12.3]. Sein Wesen ist äußerst einfach. Der Algorithmus sollte Berechnungen auf mehreren Ebenen unterstützen: eine schnelle vorläufige Berechnung (in den meisten Fällen stellt sie sich als endgültig heraus), eine langsamere, aber genauere Berechnung (spart in den meisten kritischen Fällen), eine noch langsamere, aber noch genauere Berechnung (wenn sie absolut „schlecht“ ist "Musste) und so weiter.



In der überwiegenden Mehrheit der Fälle reicht es aus, die Genauigkeit etwas höher als 0,5 ulp zu halten, aber wenn eine "Lokomotive" erscheint, erhöhen wir sie. Solange die "Dampflokomotive" erhalten bleibt, erhöhen wir die Genauigkeit, bis klar ist, dass weitere Schwankungen der numerischen Methode diese "Dampflokomotive" nicht beeinflussen. Wenn wir zum Beispiel in unserem Fall ESP = 2 ^{-54 erreicht haben} , erscheint an der 54. Position eine Einheit, die die Lok sozusagen vor Nullen "schützt" und garantiert, dass es keine Subtraktion eines Wertes größer oder gleich 2 ^{-53 mehr gibt} und Nullen werden nicht zu Einsen, wodurch das Rundungsbit auf Null gezogen wird.



Es war eine populärwissenschaftliche Präsentation, egal wie der Rundungstest von Ziv, bei der gezeigt wird, wie schnell in einem Schritt überprüft werden kann, ob wir die gewünschte Genauigkeit erreicht haben. Lesen Sie in [1, Kapitel 12] oder in [3, Abschnitt 10.5].



Das Problem bei diesem Ansatz liegt auf der Hand. Es ist notwendig, einen Algorithmus zur Berechnung jeder transzendentalen Funktion f (x) zu entwerfen, damit im Verlauf des Stücks die Genauigkeit der Berechnungen erhöht werden kann. Für die Software-Implementierung ist dies immer noch nicht so beängstigend. Beispielsweise ermöglicht die Newton-Methode grob gesagt, die Anzahl der exakten Bits nach dem Dezimalpunkt bei jeder Iteration zu verdoppeln. Sie können verdoppeln, bis es "genug" wird, obwohl dies ein ziemlich zeitaufwändiger Prozess ist, muss ich zugeben, dass Newtons Methode nicht immer gerechtfertigt ist, da sie die Berechnung der Umkehrfunktion f ^-1 erfordert(x), was in einigen Fällen nicht einfacher sein kann als die Berechnung von f (x) selbst. Für die Hardware-Implementierung ist die "Ziva-Strategie" völlig ungeeignet. Der im Prozessor fest verdrahtete Algorithmus muss eine Reihe von Aktionen mit der bereits voreingestellten Anzahl von Bits ausführen. Dies ist recht problematisch zu implementieren, wenn wir diese Anzahl nicht im Voraus kennen. Bestandsaufnahme? Und wie viel?



Der probabilistische Ansatz zur Lösung des Problems [1, Abschnitt 12.6] ermöglicht es uns, den Wert von m zu schätzen (denken Sie daran, dies ist die Anzahl der Bits, die für eine korrekte Rundung ausreicht). Es stellt sich heraus, dass die Länge der "Lokomotive" im probabilistischen Sinne etwas größer ist als die Länge der Mantisse der Zahl. In den meisten Fällen reicht es daher aus, m etwas mehr als das Doppelte des Wertes der Mantisse zu nehmen, und nur in sehr seltenen Fällen ist es notwendig, noch mehr zu nehmen. Ich zitiere die Autoren dieser Arbeit: "Wir schließen daraus, dass m in der Praxis etwas größer als 2p sein muss" (sie haben p - die Länge der Mantisse zusammen mit dem ganzzahligen Teil, dh p = 24 für float). Weiter im Text zeigen sie, dass die Fehlerwahrscheinlichkeit bei einer solchen Strategie nahe Null ist, aber immer noch positiv, und dies wird durch Experimente bestätigt.



Dennoch gibt es immer noch Fälle, in denen der Wert von m noch weiter genommen werden muss und der schlimmste Fall nicht im Voraus bekannt ist. Es gibt theoretische Schätzungen für den schlimmsten Fall [1, Abschnitt 12.7.2], aber sie liefern undenkbare Millionen von Bits, was nicht gut ist. Hier ist eine Tabelle aus der zitierten Arbeit (dies ist für die Funktion exp (x) im Intervall von -ln (2) bis ln (2)):

Zweite "Krücke"

In der Praxis wird m nicht so schrecklich groß sein. Und um den schlimmsten Fall zu bestimmen, wird eine zweite "Krücke" angewendet, die als "erschöpfende Vorberechnung" bezeichnet wird. Wenn für den Datentyp float (32 Bit) die Funktion f ein Argument (x) hat, können wir alle möglichen Werte von x leicht "ausführen". Das Problem tritt nur bei Funktionen auf, die mehr als ein Argument haben (darunter pow (x, y)), für die wir uns so etwas nicht vorstellen konnten. Nachdem wir alle möglichen Werte von x überprüft haben, berechnen wir unsere Konstante m für jede Funktion f (x) und für jeden Rundungsmodus. Dann sind die Berechnungsalgorithmen, die in Hardware implementiert werden müssen, so ausgelegt, dass sie eine Genauigkeit von 2 ^m liefern . Dann ist die Rundung von f (x) in allen Fällen garantiert korrekt.



Bei einem Doppeltyp (64 Bit) ist eine einfache Aufzählung fast unmöglich. Sie sortieren jedoch! Aber wie? Die Antwort ist in [4] gegeben. Ich werde es Ihnen ganz kurz erzählen.



Die Domäne der Funktion f (x) ist in sehr kleine Segmente unterteilt, so dass innerhalb jedes Segments f (x) durch eine lineare Funktion der Form b-ax ersetzt werden kann (die Koeffizienten a und b sind natürlich für verschiedene Segmente unterschiedlich). Die Größe dieser Segmente wird analytisch berechnet, so dass eine solche lineare Funktion in jedem Segment tatsächlich kaum vom Original zu unterscheiden wäre.



Nach einigen Skalierungs- und Verschiebungsoperationen kommen wir dann zu folgendem Problem: Kann eine gerade Linie b-ax "nahe genug" an einen ganzzahligen Punkt heranreichen?



Es stellt sich heraus, dass es relativ einfach ist, eine Ja- oder Nein-Antwort zu geben. Das heißt "Ja" - wenn potenziell gefährliche Punkte nahe einer geraden Linie liegen, und "Nein" - wenn kein solcher Punkt im Prinzip nahe an die Linie kommen kann. Das Schöne an dieser Methode ist, dass die Antwort "Nein" in der Praxis in den allermeisten Fällen erhalten wird und die Antwort "Ja", die selten erhalten wird, Sie dazu zwingt, das Segment mit einer umfassenden Suche zu durchlaufen, um festzustellen, welche spezifischen Punkte sich als kritisch herausstellten.



Das Durchlaufen der Argumente zu f (x) wird jedoch um ein Vielfaches reduziert und ermöglicht das Erkennen von Bruchstellen für Zahlen wie double (binary64) und long double (80 Bit!). Dies geschieht auf Supercomputern und natürlich auf Grafikkarten ... in Ihrer Freizeit vom Bergbau. Bisher weiß jedoch noch niemand, was mit dem Datentyp binary128 zu tun ist. Ich möchte Sie daran erinnern, dass der Bruchteil der Mantisse solcher Zahlen 112 Bit beträgt . Daher kann man in der bisherigen ausländischen Literatur zu diesem Thema nur halbphilosophische Argumente finden, die mit „wir hoffen ...“ („wir hoffen ...“) beginnen.



Die Details der Methode, mit der Sie den Durchgang einer Linie in der Nähe von ganzzahligen Punkten schnell bestimmen können, sind hier unangemessen. Für diejenigen, die den Prozess genauer lernen möchten, empfehle ich , sich beispielsweise in Artikel [5] mit dem Problem zu befassen, den Abstand zwischen einer geraden Linie und Z ^{2 zu ermitteln} . Es beschreibt einen verbesserten Algorithmus, der im Laufe der Konstruktion dem berühmten Euklid-Algorithmus zum Finden des größten gemeinsamen Teilers ähnelt. Ich werde das gleiche Bild aus [4] und [5] geben, das die weitere Transformation des Problems zeigt:







Es gibt umfangreiche Tabellen, die die schlimmsten Fälle von Rundungen in unterschiedlichen Intervallen für jede transzendentale Funktion enthalten. Sie befinden sich in [1 Abschnitt 12.8.4] und in [3, Abschnitt 10.5.3.2] sowie in separaten Artikeln, beispielsweise in [6].



Ich werde einige Beispiele geben, indem ich zufällige Zeilen aus solchen Tabellen nehme. Ich betone, dass dies nicht die schlimmsten Fälle für alle x sind, sondern nur für einige kleine Intervalle, siehe die Quelle, wenn Sie interessiert sind.







Wie lese ich die Tabelle? Der Wert x wird in hexadezimaler Gleitkomma-Doppelschreibweise angegeben. Zuerst gibt es erwartungsgemäß eine führende, dann 52 Bits des Bruchteils der Mantisse und des Buchstabens P. Dieser Buchstabe bedeutet "multiplizieren mit zwei zu einer Potenz", gefolgt von einem Grad. Zum Beispiel P-23 Mittel die angegebenen Mantisse muss von 2 multipliziert werden , ^-23 .



Stellen Sie sich außerdem vor, dass die Funktion f (x) mit unendlicher Genauigkeit berechnet wird und die ersten 53 Bits davon abgeschnitten werden (ohne Rundung!). Es sind diese 53 Bits (eines davon bis zum Komma), die in der Spalte f (x) angegeben sind. Nachfolgende Bits werden in der letzten Spalte angezeigt. Das "Grad" -Zeichen der Bitfolge in der letzten Spalte bedeutet die Anzahl der Bitwiederholungen, dh beispielsweise 10 ⁵³1011 bedeutet, dass zuerst das Bit gleich 1 kommt, dann 53 Nullen und dann 1011. Dann die Ellipse, was bedeutet, dass wir im Allgemeinen den Rest der Bits überhaupt nicht benötigen.



Außerdem ist es eine Frage der Technologie - wir kennen die schlimmsten Fälle für jedes Intervall einer separat genommenen Funktion und können für dieses Intervall eine solche Annäherung wählen, damit sie den schlimmsten Fall mit ihrer Genauigkeit abdeckt. Mit nur Jahren Supercomputer-Computing ist es möglich, schnelle und genaue Hardware-Implementierungen von Elementarfunktionen zu erstellen. Die Sache ist klein: Es bleibt zumindest den Compiler-Entwicklern beizubringen, diese Tabellen zu verwenden.

Warum wird das benötigt?

Gute Frage! Immerhin habe ich oben wiederholt gesprochen, dass fast 100% der Programmierer keine Elementarfunktion mit einer Genauigkeit auf das korrekt gerundete letzte Bit kennen müssen (oft benötigen sie nicht einmal die Hälfte der Bits). Warum fahren Wissenschaftler Supercomputer und kompilieren Tabellen, um ein „nutzloses“ Problem zu lösen?



Erstens ist die Herausforderung von grundlegender Bedeutung. Es ist ziemlich interessant, keine exakte Rundung zu erhalten, um eine genaue Rundung zu erreichen, sondern im Prinzip zu verstehen, wie dieses interessante Problem gelöst werden kann. Welche Geheimnisse der Computermathematik wird uns ihre Lösung enthüllen? Wie könnten diese Geheimnisse für andere Aufgaben verwendet werden? Grundlagenwissenschaften - sie sind so, man kann jahrzehntelang eine Art "Unsinn" machen, und dann, hundert Jahre später, dank dieses "Unsinns", findet in einem anderen Bereich ein wissenschaftlicher Durchbruch statt.



Zweitens das Problem der Code-Portabilität. Wenn es sich eine Funktion leisten kann, die letzten Bits des Ergebnisses so zu behandeln, wie sie es möchte, bedeutet dies, dass auf verschiedenen Plattformen und auf verschiedenen Compilern leicht unterschiedliche Ergebnisse erzielt werden können (selbst wenn sie innerhalb des angegebenen Fehlers liegen). In einigen Fällen ist dies nicht wichtig, in einigen Fällen kann es jedoch von Bedeutung sein, insbesondere wenn das Programm einen Fehler aufweist, der auf einer Plattform angezeigt wird, auf einer anderen Plattform jedoch nicht genau aufgrund der unterschiedlichen Bits des Ergebnisses. Aber warum beschreibe ich Ihnen die bekannten Kopfschmerzen, die mit unterschiedlichem Programmverhalten verbunden sind? Du weißt das alles ohne mich. Es wäre großartig, ein mathematisches System zu haben, das auf allen Plattformen genau gleich funktioniert, egal wie kompiliert es ist. Das müssen Sie richtig machen runden Sie das letzte Stück ab.

Liste der Quellen

[1] Jean-Michel Muller, „Elementare Funktionen: Algorithmen und Implementierung“, 2016



[2] William Kahan, „ Ein um die Hälfte zu kluger Logarithmus “, 2004



[3] Jean-Michel Muller, „Handbuch der Gleitkomma-Arithmetik“ , 2018



[4] Vincent Lefèvre, Jean-Michel Müller, "Auf dem Weg zu korrekt gerundeten Transzendentalen", IEEE TRANSACTIONS ON COMPUTERS, VOL. 47, NO. 11, NOVEMBER 1998. pp. 1235-1243



[5] Vincent Lefèvre. "Neue Ergebnisse zum Abstand zwischen einem Segment und Z ² ". Anwendung auf die exakte Rundung. 17. IEEE-Symposium für Computerarithmetik - Arith'17, Juni 2005, Cape Cod, MA,

USA. S. 68-75



[6] Vincent Lefèvre, Jean-Michel Müller, „Schlimmste Fälle für eine korrekte Rundung der Elementarfunktionen in doppelter Präzision“, Bericht (INSTITUT NA TIONAL DE RECHERCHE EN INFORMA TIQUE ET EN AUTOMA TIQUE) Nr. 4044 - November 2000 - 19 Seiten.