Algorithmus von boost :: math :: erf

Question

Gibt es irgendwelche Details für den Algorithmus hinter der ERF-Funktion von Boost? Die Dokumentation des Moduls ist nicht sehr präzise. Alles, was ich herausgefunden habe, ist, dass mehrere Methoden gemischt sind. Für mich sieht es aus wie Variationen von Abramowitz und Stegun.

• Welche Methoden sind gemischt?
• Wie sind die Methoden gemischt?
• Was ist die Komplexität der Erf-Funktion (konstante Zeit)?

Sebastian

Answer 1

Die Dokumentation für Boost Math Toolkit hat eine lange Liste von references, unter denen Abramowitz und Stegun. Die Schnittstelle erf-function enthält einen Vorlagenparameter policy, mit dem die numerische Genauigkeit (und damit die Laufzeitkomplexität) gesteuert werden kann.

#include <boost/math/special_functions/erf.hpp> 
namespace boost{ namespace math{ 

template <class T> 
calculated-result-type erf(T z); 

template <class T, class Policy> 
calculated-result-type erf(T z, const Policy&); 

template <class T> 
calculated-result-type erfc(T z); 

template <class T, class Policy> 
calculated-result-type erfc(T z, const Policy&); 

}} // namespaces 

UPDATE:

Unterhalb einer wörtlich Kopie des Abschnitts "Umsetzung" der früheren vorgesehen Bezug auf die erf-Funktion:

Implementierung

Alle Versionen von Diese Funktionen verwenden zuerst die üblichen Reflexionsformeln, um ihre Argumente positiv zu machen:

erf(-z) = 1 - erf(z); 

erfc(-z) = 2 - erfc(z); // preferred when -z < -0.5 

erfc(-z) = 1 + erf(z); // preferred when -0.5 <= -z < 0 

Die generischen Versionen dieser Funktionen sind im Hinblick auf die unvollständige Gammafunktion implementiert.

Wenn die Signifikanz (Mantisse) Größe erkannt wird (derzeit für 53, 64 und 113-Bit Reelle, plus Single-Precision 24-Bit behandelt über die Förderung zu verdoppeln), dann eine Reihe von rationalen Näherungen von JM entwickelt.

Für z < = 0,5 dann eine rationale Annäherung an erf verwendet wird, basiert auf der Beobachtung, daß erf eine ungerade Funktion ist und daher erf berechnet unter Verwendung von:

erf(z) = z * (C + R(z*z)); 

wo die rationale Approximation R (z * z) ist auf absoluten Fehler optimiert: Solange sein absoluter Fehler klein genug ist im Vergleich zu der Konstante C, wird jeder Rundungsfehler, der während der Berechnung von R (z * z) auftritt, effektiv aus dem Ergebnis verschwinden. Infolgedessen ist der Fehler für erf und erfc in diesem Bereich sehr niedrig: das letzte Bit ist nur in einer sehr kleinen Anzahl von Fällen falsch.

Für z> 0,5 wir beobachten, daß über einen kleinen Intervall [a, b) dann:

erfc(z) * exp(z*z) * z ~ c 

für eine Konstante c.

Daher wird für z> 0,5 berechnen wir erfc Verwendung:

erfc(z) = exp(-z*z) * (C + R(z - B))/z; 

Wieder R (z - B) ist für absoluten Fehler optimiert, und die Konstante C ist der Durchschnitt der erfc (z) * exp (z * z) * z an den Endpunkten des Bereichs.Noch einmal, solange der absolute Fehler in R (z - B) im Vergleich zu c klein ist, dann wird c + R (z - B) korrekt gerundet, und der Fehler im Ergebnis hängt nur von der Genauigkeit des exp ab Funktion. In der Praxis beschränkt sich der Fehler in allen außer einer sehr kleinen Anzahl von Fällen auf das letzte Bit des Ergebnisses. Die Konstante B wird so gewählt, daß das linke Ende des Bereichs der vernünftigen Näherung ist 0.

Für große z über einen Bereich [a, + ∞] Die obige Annäherung wird modifiziert, um:

erfc(z) = exp(-z*z) * (C + R(1/z))/z; 

Die rationalen Näherungen sind in excruciating detail erläutert. Wenn Sie weitere Informationen benötigen, können Sie sich immer die source code ansehen.