Ich hatte vor einiger Zeit einen Code geschrieben, der versuchte, ohne Verwendung von Bibliotheksfunktionen zu berechnen . Gestern habe ich den alten Code überprüft und versucht, ihn so schnell wie möglich (und korrekt) zu machen. Hier ist mein bisheriger Versuch:
const double ee = exp(1);
double series_ln_taylor(double n){ /* n = e^a * b, where a is an non-negative integer */
double lgVal = 0, term, now;
int i, flag = 1;
if ( n <= 0 ) return 1e-300;
if ( n * ee < 1 )
n = 1.0 / n, flag = -1; /* for extremely small n, use e^-x = 1/n */
for ( term = 1; term < n ; term *= ee, lgVal++ );
n /= term;
/* log(1 - x) = -x - x**2/2 - x**3/3... */
n = 1 - n;
now = term = n;
for ( i = 1 ; ; ){
lgVal -= now;
term *= n;
now = term / ++i;
if ( now < 1e-17 ) break;
}
if ( flag == -1 ) lgVal = -lgVal;
return lgVal;
}
Hier versuche ich, zu finden , so dass etwas mehr als n, und dann füge ich den Logarithmus Wert von , das kleiner als 1 ist. Zu diesem Zeitpunkt kann die Taylor-Expansion vonlog(1-x)ohne Bedenken verwendet werden.
Ich habe in letzter Zeit ein Interesse an numerischen Analysen entwickelt, und deshalb kann ich mir die Frage stellen, wie viel schneller dieses Codesegment in der Praxis ausgeführt werden kann, während es korrekt genug ist. Muss ich auf einige andere Methoden wechseln, beispielsweise Kettenbruch verwenden, wie dies ?
Die mit der C-Standardbibliothek bereitgestellte Funktion ist fast 5,1-mal schneller als diese Implementierung.
UPDATE 1 : Unter Verwendung der in Wikipedia erwähnten hyperbolischen Arctan-Reihe scheint die Berechnung fast 2,2-mal langsamer zu sein als die Protokollfunktion der C-Standardbibliothek. Obwohl ich die Leistung nicht ausführlich überprüft habe, scheint meine aktuelle Implementierung bei größeren Zahlen WIRKLICH langsam zu sein. Ich möchte sowohl meine Implementierung auf fehlergebundene als auch die durchschnittliche Zeit für eine Vielzahl von Zahlen überprüfen, wenn ich dies verwalten kann. Hier ist meine zweite Anstrengung.
double series_ln_arctanh(double n){ /* n = e^a * b, where a is an non-negative integer */
double lgVal = 0, term, now, sm;
int i, flag = 1;
if ( n <= 0 ) return 1e-300;
if ( n * ee < 1 ) n = 1.0 / n, flag = -1; /* for extremely small n, use e^-x = 1/n */
for ( term = 1; term < n ; term *= ee, lgVal++ );
n /= term;
/* log(x) = 2 arctanh((x-1)/(x+1)) */
n = (1 - n)/(n + 1);
now = term = n;
n *= n;
sm = 0;
for ( i = 3 ; ; i += 2 ){
sm += now;
term *= n;
now = term / i;
if ( now < 1e-17 ) break;
}
lgVal -= 2*sm;
if ( flag == -1 ) lgVal = -lgVal;
return lgVal;
}
Jeder Vorschlag oder jede Kritik wird geschätzt.
double series_ln_better(double n){ /* n = e^a * b, where a is an non-negative integer */
double lgVal = 0, term, now, sm;
int i, flag = 1;
if ( n == 0 ) return -1./0.; /* -inf */
if ( n < 0 ) return 0./0.; /* NaN*/
if ( n < 1 ) n = 1.0 / n, flag = -1; /* for extremely small n, use e^-x = 1/n */
/* the cutoff iteration is 650, as over e**650, term multiplication would
overflow. For larger numbers, the loop dominates the arctanh approximation
loop (with having 13-15 iterations on average for tested numbers so far */
for ( term = 1; term < n && lgVal < 650 ; term *= ee, lgVal++ );
if ( lgVal == 650 ){
n /= term;
for ( term = 1 ; term < n ; term *= ee, lgVal++ );
}
n /= term;
/* log(x) = 2 arctanh((x-1)/(x+1)) */
n = (1 - n)/(n + 1);
now = term = n;
n *= n;
sm = 0;
/* limiting the iteration for worst case scenario, maximum 24 iteration */
for ( i = 3 ; i < 50 ; i += 2 ){
sm += now;
term *= n;
now = term / i;
if ( now < 1e-17 ) break;
}
lgVal -= 2*sm;
if ( flag == -1 ) lgVal = -lgVal;
return lgVal;
}
quelle
frexpKirills Antwort berührte bereits eine Vielzahl relevanter Themen. Ich möchte einige davon auf der Grundlage praktischer Erfahrungen im Entwurf von Mathematikbibliotheken erweitern. Ein Hinweis im Vorfeld: Entwickler von Mathematikbibliotheken verwenden in der Regel jede veröffentlichte algorithmische Optimierung sowie viele maschinenspezifische Optimierungen, von denen nicht alle veröffentlicht werden. Der Code wird häufig in Assemblersprache geschrieben, anstatt kompilierten Code zu verwenden. Es ist daher unwahrscheinlich, dass eine einfache und kompilierte Implementierung mehr als 75% der Leistung einer vorhandenen hochwertigen Mathematikbibliotheksimplementierung erzielt, wenn identische Funktionssätze vorausgesetzt werden (Genauigkeit, Behandlung von Sonderfällen, Fehlerberichterstattung, Unterstützung im Rundungsmodus).
Die Genauigkeit wird in der Regel durch Vergleich mit einer Referenz von (Drittanbietern) mit höherer Genauigkeit bewertet. Funktionen mit einfacher Genauigkeit und einfachem Argument können leicht ausführlich getestet werden, andere Funktionen erfordern das Testen mit (gerichteten) zufälligen Testvektoren. Natürlich kann man nicht unendlich genaue Referenzergebnisse berechnen, aber die Untersuchung des Table-Maker-Dilemmas legt nahe, dass es für viele einfache Funktionen ausreicht, eine Referenz mit einer Genauigkeit von etwa dem Dreifachen der Zielgenauigkeit zu berechnen. Siehe zum Beispiel:
Vincent Lefèvre, Jean-Michel Müller, "Schlimmste Fälle für eine korrekte Rundung der Elementarfunktionen in doppelter Präzision". In Proceedings 15. IEEE Symposium on Computer Arithmetic , 2001, 111-118). (Preprint online)
In Bezug auf die Leistung muss zwischen der Optimierung der Latenz (wichtig, wenn man die Ausführungszeit abhängiger Operationen betrachtet) und der Optimierung des Durchsatzes (relevant, wenn die Ausführungszeit unabhängiger Operationen berücksichtigt wird) unterschieden werden. In den letzten zwanzig Jahren hat die Verbreitung von Hardware-Parallelisierungstechniken wie Parallelität auf Befehlsebene (z. B. superskalar, Prozessoren außerhalb der Reihenfolge), Parallelität auf Datenebene (z. B. SIMD-Befehle) und Parallelität auf Thread-Ebene (z. B. Hyper-Threading) zugenommen. Multi-Core-Prozessoren) hat dazu geführt, dass der Rechendurchsatz als relevantere Metrik im Vordergrund steht.
Die Fused Multiply-Add-Operation ( FMA ), die vor 25 Jahren von IBM eingeführt wurde und jetzt auf allen wichtigen Prozessorarchitekturen verfügbar ist, ist ein entscheidender Baustein für die Implementierung moderner Mathematikbibliotheken. Es bietet eine Reduzierung von Rundungsfehlern, einen begrenzten Schutz gegen subtraktive Löschung und vereinfacht die Doppel-Doppel-Arithmetik erheblich .
C99log()C99fma()quelle