Die komplementäre Fehlerfunktion erfc ist eine spezielle Funktion, die eng mit der Standardnormalverteilung verknüpft ist. Es wird häufig in der Statistik und den Naturwissenschaften (z. B. Diffusionsproblemen) verwendet, wo die "Schwänze" dieser Verteilung berücksichtigt werden müssen, und die Verwendung der Fehlerfunktion erf ist daher nicht geeignet.
Die komplementäre Fehlerfunktion wurde in der Standard-Mathematikbibliothek ISO C99 als Funktionen erfcf , erfc und erfcl ; Diese wurden anschließend in ISO C ++ übernommen. Somit kann Quellcode leicht in Open-Source-Implementierungen dieser Bibliothek gefunden werden, zum Beispiel in glibc .
Viele bestehende Implementierungen sind jedoch skalarer Natur, während moderne Prozessorhardware SIMD-orientiert ist (entweder explizit, wie in x86-CPU oder implizit wie in GPUs). Aus Gründen der Performance ist daher eine vektorisierbare Implementierung sehr wünschenswert. Dies bedeutet, dass Verzweigungen vermieden werden müssen, außer im Rahmen einer ausgewählten Zuweisung. Ebenso ist eine umfangreiche Verwendung von Tabellen nicht angezeigt, da eine parallelisierte Suche oft ineffizient ist.
Wie würde man eine effiziente vektorisierbare Implementierung der Single-Precision-Funktion erfcf() erstellen? Die Genauigkeit, gemessen in ulp , sollte ungefähr der glibc skalaren Implementierung entsprechen, die einen maximalen Fehler von 3.12575 ulps aufweist (bestimmt durch erschöpfende Tests). Die Verfügbarkeit von Fused Multiply-Add (FMA) kann angenommen werden, da dies zu diesem Zeitpunkt von allen gängigen Prozessorarchitekturen (CPUs und GPUs) angeboten wird. Während die Verarbeitung von Floating-Point-Statusflags und errno ignoriert werden kann, sollten Denormals, Unfinities und NaNs gemäß den IEEE 754-Bindungen für ISO C behandelt werden.
Nachdem wir verschiedene Ansätze untersucht haben, ist der Algorithmus, der am besten geeignet erscheint, der folgende Algorithmus:
M. M. Shepherd und JG Laframboise, "Chebyshev Approximation von (1 + 2 x ) exp ( x 2 ) erfc x in 0 ≤ x & lt; ∞. " Mathematik der Berechnung , Band 36, Nr. 153, Januar 1981, S. 249-253 ( Online-Kopie )
Der Grundgedanke des Artikels besteht darin, eine Annäherung an (1 + 2 x ) exp ( x 2 ) erfc (< em> x ), von dem wir erfcx ( x ) berechnen können, indem wir einfach durch (1 + 2 x ) und dann erfc (x) dividieren Multiplikation mit exp (- x 2 ). Der eng begrenzte Bereich der Funktion mit Funktionswerten grob in [1, 1.3] und seine allgemeine "Flachheit" eignen sich gut für die polynome Approximation. Numerische Eigenschaften dieses Ansatzes werden durch die Eingrenzung des Approximationsintervalls weiter verbessert: Das ursprüngliche Argument x wird durch q = ( x - K) / ( x + K), wobei K eine geeignet gewählte Konstante ist, gefolgt von der Berechnung p ( q ), wobei p ist ein Polynom.
Da erfc - x = 2 - erfc x müssen wir nur das Intervall [0, ∞] berücksichtigen, das dem Intervall [-1, 1 zugeordnet ist ] durch diese Transformation. Bei IEEE-754-Einfachgenauigkeit verschwindet erfcf() (wird zu Null) für x & gt; 10.0546875, so muss man nur x ∈ [0, 10.0546875) berücksichtigen. Was ist der "optimale" Wert von K für diesen Bereich? Ich kenne keine mathematische Analyse, die die Antwort liefern würde, das Papier schlägt K = 3,75 basierend auf Experimenten vor.
Man kann leicht feststellen, dass für Berechnungen mit einfacher Genauigkeit eine Minimax-Polynom-Approximation vom Grad 9 für verschiedene Werte von K in dieser allgemeinen Umgebung ausreichend ist. Durch systematisches Erzeugen solcher Näherungen mit dem Remez-Algorithmus, wobei K zwischen 1,5 und 4 in Schritten von 1/16 variiert, wird der niedrigste Näherungsfehler für K = {2, 2,625, 3,3125} beobachtet. Von diesen ist K = 2 die vorteilhafteste Wahl, da es sich für eine sehr genaue Berechnung von ( x - K) / ( x + K) eignet, wie in diese Frage .
Der Wert K = 2 und die Eingabedomäne für x legen nahe, dass Variante 4 von meine Antwort , jedoch kann man einmal experimentell nachweisen, dass die günstigere Variante 5 hier die gleiche Genauigkeit erreicht, was wahrscheinlich auf die sehr flache Steigung der approximierten Funktion für q & gt; -0.5, wodurch sich der Fehler im Argument q um etwa den Faktor zehn verringert.
Da die Berechnung von erfc() zusätzlich zur anfänglichen Approximation Nachverarbeitungsschritte erfordert, ist es klar, dass die Genauigkeit dieser beiden Berechnungen hoch sein muss, um ein ausreichend genaues Endergebnis zu erhalten. Fehlerkorrekturtechniken müssen verwendet werden.
Man beobachtet, dass der bedeutendste Koeffizient in der Polynomapproximation von (1 + 2 x ) exp ( x 2 ) erfc (< em> x ) hat die Form (1 + s), wobei s & lt; 0.5. Dies bedeutet, dass wir den führenden Koeffizienten genauer darstellen können, indem wir 1 abspalten und nur s im Polynom verwenden. Anstatt also ein Polynom p (q) zu berechnen, dann multipliziert man es mit dem reziproken r = 1 / (1 + 2 x ), ist es mathematisch äquivalent, aber numerisch vorteilhaft zu berechnen die Kernannäherung als p ( q ) + 1, und verwenden Sie p , um fma (p, r, r) zu berechnen.
Die Genauigkeit der Division kann verbessert werden, indem ein Anfangsquotient q aus dem reziproken r berechnet wird, um den Rest e = p +1 - q * (1 + 2 x ) mit Hilfe einer FMA, dann e anwenden die Korrektur q = q + ( e * r ), wiederum unter Verwendung einer FMA.
Die Potenzierung hat Fehlervergrößerungseigenschaften, daher muss die Berechnung von e - x 2 sorgfältig durchgeführt werden. Die Verfügbarkeit von FMA erlaubt trivialerweise die Berechnung von - x 2 als ein doppeltes float s hohes : s niedriges . e x ist seine eigene Ableitung, so dass man e s high berechnen kann: s low als e s hoch + e s hoch * s niedrig . Diese Berechnung kann mit der Multiplikation des vorherigen Zwischenergebnisses r kombiniert werden, um r = r * e s hoch zu ergeben + r * e s hoch * s niedrig .Durch Verwendung von FMA wird sichergestellt, dass der höchstwertige Term r * e s high so genau wie möglich berechnet wird.
Kombiniert man die obigen Schritte mit ein paar einfachen Auswahlen, um Ausnahmefälle und negative Argumente zu behandeln, kommt man zum folgenden C-Code:
%Vor% Ich habe meine eigene Implementierung von expf() im obigen Code verwendet, um meine Arbeit von Unterschieden in den expf() Implementierungen auf verschiedenen Computerplattformen zu isolieren. Aber jede Implementierung von expf() , deren maximaler Fehler in der Nähe von 0,5 ul liegt, sollte gut funktionieren. Wie oben gezeigt, hat my_expf() bei Verwendung von my_erfcf() einen maximalen Fehler von 2.65712 ulps. Vorausgesetzt, dass ein vektorisierbarer expf() verfügbar ist, sollte der obige Code problemlos vektorisieren. Ich habe eine kurze Überprüfung mit dem Intel-Compiler 13.1.3.198 gemacht. Ich habe einen Aufruf von my_erfcf() in einer Schleife aufgerufen, #include <mathimf.h> hinzugefügt, den Aufruf von my_expf() durch einen Aufruf von expf() ersetzt und dann mit diesen Befehlszeilenoptionen kompiliert:
Der Intel-Compiler berichtete, dass die Schleife vektorisiert wurde, was ich durch Überprüfung des disassemblierten Binärcodes zweimal überprüft habe.
Da my_erfcf() nur reziproke statt vollständige Divisionen verwendet, ist es der Verwendung schneller reziproker Implementierungen zugänglich, vorausgesetzt, sie liefern nahezu korrekt gerundete Ergebnisse. Für Prozessoren, die eine schnelle reziproke Approximation mit einfacher Genauigkeit in Hardware bereitstellen, kann dies leicht erreicht werden, indem dies mit einer Halley-Iteration mit kubischer Konvergenz gekoppelt wird. Ein (skalares) Beispiel für diesen Ansatz für x86-Prozessoren ist:
Tags und Links algorithm c math floating-point