Wie berechnet C sin () und andere mathematische Funktionen?

248

Ich habe .NET-Disassemblies und den GCC-Quellcode durchgesehen, kann aber anscheinend nirgendwo die tatsächliche Implementierung sin()und andere mathematische Funktionen finden ... sie scheinen immer auf etwas anderes zu verweisen.

Kann mir jemand helfen, sie zu finden? Ich halte es für unwahrscheinlich, dass ALLE Hardware, auf der C ausgeführt wird, Triggerfunktionen in der Hardware unterstützt. Es muss also irgendwo einen Softwarealgorithmus geben , oder?


Ich bin mir dessen bewusst , dass mehrere Möglichkeiten , Funktionen können berechnet werden, und habe meine eigene Routinen geschrieben Rechenfunktionen Taylorreihe für Spaß. Ich bin gespannt, wie real Produktionssprachen das machen, da alle meine Implementierungen immer um mehrere Größenordnungen langsamer sind, obwohl ich meine Algorithmen für ziemlich clever halte (offensichtlich nicht).

Strang
quelle
2
Bitte beachten Sie, dass diese Implementierung abhängig ist. Sie sollten angeben, an welcher Implementierung Sie am meisten interessiert sind.
Jason
3
Ich habe .NET und C markiert, weil ich an beiden Stellen gesucht habe und es auch nicht herausgefunden habe. Wenn man sich die .NET-Disassemblierung ansieht, sieht es so aus, als würde es sich um nicht verwaltetes C handeln, soweit ich weiß, haben sie dieselbe Implementierung.
Hank

Antworten:

213

In GNU libm ist die Implementierung von sinsystemabhängig. Daher finden Sie die Implementierung für jede Plattform irgendwo im entsprechenden Unterverzeichnis von sysdeps .

Ein Verzeichnis enthält eine Implementierung in C, die von IBM bereitgestellt wurde. Seit Oktober 2011 ist dies der Code, der tatsächlich ausgeführt wird, wenn Sie sin()ein typisches x86-64-Linux-System aufrufen . Es ist anscheinend schneller als die fsinMontageanleitung. Quellcode: sysdeps / ieee754 / dbl-64 / s_sin.c , suchen Sie nach __sin (double x).

Dieser Code ist sehr komplex. Kein Software-Algorithmus ist so schnell wie möglich und auch über den gesamten Bereich von x- Werten genau. Daher implementiert die Bibliothek mehrere verschiedene Algorithmen. Ihre erste Aufgabe besteht darin, x zu betrachten und zu entscheiden, welcher Algorithmus verwendet werden soll.

  • Wenn x sehr, sehr nahe bei 0 liegt, sin(x) == xist dies die richtige Antwort.

  • Etwas weiter draußen sin(x)wird die bekannte Taylor-Serie verwendet. Dies ist jedoch nur in der Nähe von 0 genau, also ...

  • Wenn der Winkel mehr als ungefähr 7 ° beträgt, wird ein anderer Algorithmus verwendet, der Taylor-Reihen-Näherungen sowohl für sin (x) als auch für cos (x) berechnet und dann Werte aus einer vorberechneten Tabelle verwendet, um die Näherung zu verfeinern.

  • Wann | x | > 2, keiner der oben genannten Algorithmen würde funktionieren, daher beginnt der Code damit, einen Wert zu berechnen, der näher an 0 liegt und dem sinoder cosstattdessen zugeführt werden kann.

  • Es gibt noch einen weiteren Zweig, der sich damit befasst, dass x ein NaN oder eine Unendlichkeit ist.

Dieser Code verwendet einige numerische Hacks, die ich noch nie gesehen habe, obwohl sie meines Wissens unter Gleitkomma-Experten bekannt sein könnten. Manchmal brauchten einige Codezeilen mehrere Absätze, um zu erklären. Zum Beispiel diese beiden Zeilen

double t = (x * hpinv + toint);
double xn = t - toint;

werden (manchmal) verwendet, um x auf einen Wert nahe 0 zu reduzieren , der sich von x um ein Vielfaches von π / 2 unterscheidet, insbesondere xn× π / 2. Die Art und Weise, wie dies ohne Teilung oder Verzweigung geschieht, ist ziemlich klug. Aber es gibt überhaupt keinen Kommentar!


Ältere 32-Bit-Versionen von GCC / glibc verwendeten den fsinBefehl, was für einige Eingaben überraschend ungenau ist. Es gibt einen faszinierenden Blog-Beitrag, der dies mit nur 2 Codezeilen veranschaulicht .

Die Implementierung von fdlibm sinin reinem C ist viel einfacher als die von glibc und wird gut kommentiert. Quellcode: fdlibm / s_sin.c und fdlibm / k_sin.c

Jason Orendorff
quelle
35
Um zu sehen, dass dies wirklich der Code ist, der auf x86 ausgeführt wird: Kompilieren Sie ein Programm, das aufruft sin(). tippe gdb a.outdann break sin, dann run, dann disassemble.
Jason Orendorff
5
@ Henry: Machen Sie nicht den Fehler zu denken, dass dies ein guter Code ist. Es ist wirklich schrecklich , lerne nicht so zu codieren!
Thomas Bonini
2
@Andreas Hmm, Sie haben Recht, der IBM-Code sieht im Vergleich zu fdlibm ziemlich schrecklich aus. Ich habe die Antwort bearbeitet, um Links zur Sinusroutine von fdlibm hinzuzufügen.
Jason Orendorff
3
@Henry: __kernel_sinist jedoch in k_sin.c definiert und es ist reines C. Klicken Sie erneut darauf - ich habe die URL beim ersten Mal verpfuscht.
Jason Orendorff
3
Der verknüpfte Sysdeps-Code ist besonders interessant, da er korrekt gerundet ist. Das heißt, es gibt anscheinend die bestmögliche Antwort für alle Eingabewerte, was erst vor relativ kurzer Zeit möglich geworden ist. In einigen Fällen kann dies langsam sein, da möglicherweise viele zusätzliche Ziffern berechnet werden müssen, um eine korrekte Rundung sicherzustellen. In anderen Fällen ist es extrem schnell - für ausreichend kleine Zahlen ist die Antwort nur der Winkel.
Bruce Dawson
66

Funktionen wie Sinus und Cosinus sind im Mikrocode in Mikroprozessoren implementiert. Intel-Chips haben zum Beispiel Montageanleitungen für diese. Der AC-Compiler generiert Code, der diese Assembly-Anweisungen aufruft. (Im Gegensatz dazu wird ein Java-Compiler dies nicht tun. Java wertet Triggerfunktionen eher in Software als in Hardware aus und läuft daher viel langsamer.)

Chips verwenden keine Taylor-Reihen, um Triggerfunktionen zu berechnen, zumindest nicht vollständig. Zunächst verwenden sie CORDIC , aber sie können auch eine kurze Taylor-Reihe verwenden, um das Ergebnis von CORDIC zu verbessern, oder für spezielle Fälle wie die Berechnung des Sinus mit hoher relativer Genauigkeit für sehr kleine Winkel. Weitere Erläuterungen finden Sie in dieser StackOverflow-Antwort .

John D. Cook
quelle
10
Transzendentale mathematische Funktionen wie Sinus und Cosinus können in Mikrocode oder als Hardwareanweisungen in aktuellen 32-Bit-Desktop- und Serverprozessoren implementiert werden. Dies war nicht immer der Fall, bis beim i486 (DX) alle Gleitkommaberechnungen in Software ("Soft-Float") für die x86-Serie ohne separaten Coprozessor durchgeführt wurden. Nicht alle (FPUs) enthielten transzendentale Funktionen (z. B. Weitek 3167).
mctylr
1
Kannst du genauer sein? Wie "poliert" man eine Annäherung mit einer Taylor-Reihe?
Hank
4
Angenommen, Sie berechnen sowohl Sinus als auch Cosinus, um eine Antwort zu "aufpolieren". Angenommen, Sie kennen den genauen Wert von beiden an einem Punkt (z. B. von CORDIC), möchten aber den Wert an einem nahe gelegenen Punkt. Dann können Sie für eine kleine Differenz h die Taylor-Näherungen f (x + h) = f (x) + h f '(x) oder f (x + h) = f (x) + h f' (x) anwenden. + h ^ 2 f '' (x) / 2.
John D. Cook
6
x86 / x64-Chips haben eine Montageanleitung zur Berechnung des Sinus (fsin), aber diese Anweisung ist manchmal ziemlich ungenau und wird daher selten mehr verwendet. Weitere Informationen finden Sie unter randomascii.wordpress.com/2014/10/09/… . Die meisten anderen Prozessoren tun dies keine Anweisungen für Sinus und Cosinus, da die Berechnung in Software mehr Flexibilität bietet und möglicherweise sogar schneller ist.
Bruce Dawson
3
Das Cordic-Material in den Intel-Chips wird im Allgemeinen NICHT verwendet. Erstens ist die Genauigkeit und Auflösung des Vorgangs für viele Anwendungen äußerst wichtig. Cordic ist notorisch ungenau, wenn Sie die 7. Stelle erreichen, und unvorhersehbar. Zweitens habe ich gehört, dass es einen Fehler in ihrer Implementierung gibt, der noch mehr Probleme verursacht. Ich habe mir die Sin-Funktion für Linux GCC angesehen und sie verwendet Chebyshev. Das eingebaute Zeug wird nicht verwendet. Außerdem ist der Cordic-Algorithmus im Chip langsamer als die Softwarelösung.
Donald Murray
63

OK Kinder, Zeit für die Profis ... Dies ist eine meiner größten Beschwerden bei unerfahrenen Software-Ingenieuren. Sie berechnen transzendentale Funktionen von Grund auf neu (unter Verwendung von Taylors Reihen), als hätte noch niemand in ihrem Leben diese Berechnungen durchgeführt. Nicht wahr. Dies ist ein genau definiertes Problem, das von sehr cleveren Software- und Hardware-Ingenieuren tausende Male angegangen wurde und eine genau definierte Lösung hat. Grundsätzlich verwenden die meisten transzendentalen Funktionen Chebyshev-Polynome, um sie zu berechnen. Welche Polynome verwendet werden, hängt von den Umständen ab. Erstens ist die Bibel zu diesem Thema ein Buch mit dem Titel "Computer Approximations" von Hart und Cheney. In diesem Buch können Sie entscheiden, ob Sie einen Hardware-Addierer, Multiplikator, Teiler usw. haben, und entscheiden, welche Operationen am schnellsten sind. zB Wenn Sie einen wirklich schnellen Teiler hatten, Der schnellste Weg zur Berechnung des Sinus könnte P1 (x) / P2 (x) sein, wobei P1, P2 Chebyshev-Polynome sind. Ohne den schnellen Teiler könnte es nur P (x) sein, wobei P viel mehr Terme als P1 oder P2 hat ... also wäre es langsamer. Der erste Schritt besteht also darin, Ihre Hardware zu bestimmen und festzustellen, was sie kann. Dann wählen Sie die geeignete Kombination von Chebyshev-Polynomen (hat normalerweise die Form cos (ax) = aP (x) für Cosinus, wobei P wiederum ein Chebyshev-Polynom ist). Dann entscheiden Sie, welche Dezimalgenauigkeit Sie möchten. Wenn Sie beispielsweise eine Genauigkeit von 7 Stellen wünschen, schlagen Sie dies in der entsprechenden Tabelle in dem von mir erwähnten Buch nach und erhalten (für Genauigkeit = 7,33) eine Zahl N = 4 und eine Polynomzahl 3502. N ist die Reihenfolge der Polynom (also ist es p4.x ^ 4 + p3.x ^ 3 + p2.x ^ 2 + p1.x + p0), weil N = 4. Dann suchen Sie den tatsächlichen Wert von p4, p3, p2, p1, p0-Werte im hinteren Teil des Buches unter 3502 (sie werden in Gleitkommazahlen angegeben). Dann implementieren Sie Ihren Algorithmus in Software in der Form: (((p4.x + p3) .x + p2) .x + p1) .x + p0 .... und so würden Sie den Kosinus auf 7 Dezimalstellen berechnen Orte auf dieser Hardware.

Beachten Sie, dass die meisten Hardware-Implementierungen von transzendentalen Operationen in einer FPU normalerweise einen Mikrocode und Operationen wie diese beinhalten (abhängig von der Hardware). Chebyshev-Polynome werden für die meisten Transzendentalen verwendet, aber nicht für alle. Beispiel Quadratwurzel ist schneller, wenn zuerst eine doppelte Iteration der Newton-Raphson-Methode unter Verwendung einer Nachschlagetabelle verwendet wird. Auch dieses Buch "Computer Approximations" wird Ihnen das sagen.

Wenn Sie diese Funktionen implementieren möchten, würde ich jedem empfehlen, eine Kopie dieses Buches zu erhalten. Es ist wirklich die Bibel für diese Art von Algorithmen. Beachten Sie, dass es eine Reihe alternativer Methoden zur Berechnung dieser Werte gibt, z. B. Cordics usw., diese eignen sich jedoch am besten für bestimmte Algorithmen, bei denen Sie nur eine geringe Genauigkeit benötigen. Um die Präzision jedes Mal zu gewährleisten, sind die Chebyshev-Polynome der richtige Weg. Wie gesagt, gut definiertes Problem. Wurde seit 50 Jahren gelöst ..... und so wird es gemacht.

Nun gibt es jedoch Techniken, mit denen die Chebyshev-Polynome verwendet werden können, um ein Ergebnis mit einfacher Genauigkeit mit einem Polynom niedrigen Grades zu erhalten (wie im obigen Beispiel für Cosinus). Dann gibt es andere Techniken zum Interpolieren zwischen Werten, um die Genauigkeit zu erhöhen, ohne zu einem viel größeren Polynom gehen zu müssen, wie zum Beispiel "Gal's Accurate Tables Method". Auf diese letztere Technik bezieht sich der Beitrag, der sich auf die ACM-Literatur bezieht. Aber letztendlich werden die Chebyshev-Polynome verwendet, um 90% des Weges dorthin zu erreichen.

Genießen.

Donald Murray
quelle
6
Ich konnte den ersten Sätzen nicht mehr zustimmen. Es sei auch daran erinnert, dass die Berechnung spezieller Funktionen mit garantierter Präzision ein schwieriges Problem ist . Die klugen Leute, die Sie erwähnen, verbringen den größten Teil ihres Lebens damit. Technisch gesehen sind Min-Max-Polynome das begehrte Graal, und Chebyshev-Polynome sind einfachere Proxys für sie.
Alexandre C.
161
-1 für den unprofessionellen und weitläufigen (und leicht unhöflichen) Ton und für die Tatsache, dass der tatsächliche nicht redundante Inhalt dieser Antwort, der von Streifzug und Herablassung befreit ist, im Wesentlichen auf "Sie verwenden häufig Chebyshev-Polynome; siehe dieses Buch" hinausläuft Für mehr Details ist es wirklich gut! " Was, wissen Sie, durchaus richtig sein mag, aber es ist nicht wirklich die Art von in sich geschlossener Antwort, die wir hier auf SO wollen. Auf diese Weise zusammengefasst hätte es jedoch einen anständigen Kommentar zu der Frage abgegeben.
Ilmari Karonen
2
In den frühen Jahren der Spieleentwicklung wurde dies normalerweise mit Nachschlagetabellen durchgeführt (kritisches Bedürfnis nach Geschwindigkeit). Wir haben normalerweise nicht die Standard-lib-Funktionen für diese Dinge verwendet.
Topspin
4
Ich verwende in eingebetteten Systemen häufig Nachschlagetabellen und Bittians (anstelle von Radians), aber dies ist für eine spezielle Anwendung (wie Ihre Spiele). Ich denke, der Typ interessiert sich dafür, wie der c-Compiler die Sünde für Gleitkommazahlen berechnet ...
Donald Murray
1
Ah, vor 50 Jahren. Ich fing an, mit solchen auf den Burroughs B220 mit McLaren-Serie zu spielen. Später CDC-Hardware und dann Motorola 68000. Arcsin war chaotisch - ich wählte den Quotienten aus zwei Polynomen und entwickelte Code, um die optimalen Koeffizienten zu finden.
Rick James
15

Für sinInsbesondere würde Taylorentwicklung mit Ihnen:

sin (x): = x - x ^ 3/3! + x ^ 5/5! - x ^ 7/7! + ... (1)

Sie würden so lange Begriffe hinzufügen, bis entweder der Unterschied zwischen ihnen geringer als ein akzeptiertes Toleranzniveau ist oder nur für eine begrenzte Anzahl von Schritten (schneller, aber weniger genau). Ein Beispiel wäre so etwas wie:

float sin(float x)
{
  float res=0, pow=x, fact=1;
  for(int i=0; i<5; ++i)
  {
    res+=pow/fact;
    pow*=-1*x*x;
    fact*=(2*(i+1))*(2*(i+1)+1);
  }

  return res;
}

Anmerkung: (1) funktioniert aufgrund der Annäherung sin (x) = x für kleine Winkel. Für größere Winkel müssen Sie immer mehr Terme berechnen, um akzeptable Ergebnisse zu erzielen. Sie können ein while-Argument verwenden und für eine bestimmte Genauigkeit fortfahren:

double sin (double x){
    int i = 1;
    double cur = x;
    double acc = 1;
    double fact= 1;
    double pow = x;
    while (fabs(acc) > .00000001 &&   i < 100){
        fact *= ((2*i)*(2*i+1));
        pow *= -1 * x*x; 
        acc =  pow / fact;
        cur += acc;
        i++;
    }
    return cur;

}
Blind
quelle
1
Wenn Sie die Koeffizienten ein wenig anpassen (und sie fest in ein Polynom codieren), können Sie etwa 2 Iterationen früher stoppen.
Rick James
14

Ja, es gibt auch Softwarealgorithmen zur Berechnung sin. Grundsätzlich erfolgt die Berechnung solcher Daten mit einem digitalen Computer normalerweise mit numerischen Methoden wie der Approximation der Taylor-Reihe, die die Funktion darstellt.

Numerische Methoden können Funktionen auf ein beliebiges Maß an Genauigkeit approximieren. Da das Maß an Genauigkeit, das Sie in einer schwebenden Zahl haben, endlich ist, eignen sie sich ziemlich gut für diese Aufgaben.

Mehrdad Afshari
quelle
12
Eine echte Implementierung wird wahrscheinlich keine Taylor-Serie verwenden, da es effizientere Möglichkeiten gibt. Sie müssen nur in der Domäne [0 ... pi / 2] korrekt approximieren, und es gibt Funktionen, die eine gute Approximation effizienter liefern als eine Taylor-Reihe.
David Thornley
2
@ David: Ich stimme zu. Ich war vorsichtig genug, um das Wort "wie" in meiner Antwort zu erwähnen. Die Taylor-Erweiterung ist jedoch einfach, um die Idee hinter Methoden zu erklären, die Funktionen approximieren. Trotzdem habe ich Software-Implementierungen gesehen (nicht sicher, ob sie optimiert wurden), die Taylor-Serien verwendeten.
Mehrdad Afshari
1
Tatsächlich sind Polynomnäherungen eine der effizientesten Methoden zur Berechnung trigonometrischer Funktionen.
Jeremy Salwen
13

Verwenden Sie die Taylor-Reihe und versuchen Sie, eine Beziehung zwischen den Begriffen der Reihe zu finden, damit Sie die Dinge nicht immer wieder berechnen

Hier ist ein Beispiel für Cosinus:

double cosinus(double x, double prec)
{
    double t, s ;
    int p;
    p = 0;
    s = 1.0;
    t = 1.0;
    while(fabs(t/s) > prec)
    {
        p++;
        t = (-t * x * x) / ((2 * p - 1) * (2 * p));
        s += t;
    }
    return s;
}

damit können wir den neuen Term der Summe mit dem bereits verwendeten erhalten (wir vermeiden die Fakultät und x 2p )

Erläuterung

Hannoun Yassir
quelle
2
Wussten Sie, dass Sie die Google Chart-API verwenden können, um solche Formeln mit TeX zu erstellen? code.google.com/apis/chart/docs/gallery/formulas.html
Gab Royer
11

Es ist eine komplexe Frage. Intel-ähnliche CPUs der x86-Familie haben eine Hardware-Implementierung der sin()Funktion, sie ist jedoch Teil der x87-FPU und wird im 64-Bit-Modus nicht mehr verwendet (stattdessen werden SSE2-Register verwendet). In diesem Modus wird eine Softwareimplementierung verwendet.

Es gibt mehrere solcher Implementierungen. Eine ist in fdlibm und wird in Java verwendet. Soweit ich weiß, enthält die glibc-Implementierung Teile von fdlibm und andere Teile, die von IBM bereitgestellt wurden.

Software-Implementierungen von transzendentalen Funktionen, wie sie sin()typischerweise Annäherungen durch Polynome verwenden, werden häufig aus Taylor-Reihen erhalten.

Thomas Pornin
quelle
3
SSE2 - Register sind nicht zu berechnen sin (), weder in x86 noch in x64 - Modus und natürlich verwendet, sin in Hardware berechnet wird , unabhängig von der Betriebsart. Hey, es ist 2010, in dem wir leben :)
Igor Korkhov
7
@Igor: Das hängt davon ab, welche Mathematikbibliothek Sie sich ansehen. Es stellt sich heraus, dass die am besten optimierten Mathematikbibliotheken auf x86 SSE-Softwareimplementierungen verwenden sinund cosdiese schneller sind als die Hardwareanweisungen auf der FPU. Einfachere, naivere Bibliotheken verwenden in der Regel die Anweisungen fsinund fcos.
Stephen Canon
@ Stephen Canon: Haben diese schnellen Bibliotheken eine 80-Bit-Genauigkeit wie FPU-Register? Ich habe den sehr hinterhältigen Verdacht, dass sie Geschwindigkeit gegenüber Präzision bevorzugen, was natürlich in vielen Szenarien vernünftig ist, zum Beispiel in Spielen. Und ich glaube, dass die Berechnung des Sinus mit 32-Bit-Genauigkeit unter Verwendung von SSE und vorberechneten Zwischentabellen schneller sein kann als unter Verwendung FSINmit voller Genauigkeit. Ich wäre Ihnen sehr dankbar, wenn Sie mir die Namen dieser schnellen Bibliotheken nennen würden. Es ist interessant, einen Blick darauf zu werfen.
Igor Korkhov
@Igor: Auf x86 im 64-Bit-Modus, zumindest auf allen mir bekannten Unix-ähnlichen Systemen, ist die Genauigkeit auf 64 Bit beschränkt, nicht auf die 79 Bit der x87-FPU. Die Software-Implementierung von sin()ist ungefähr doppelt so schnell wie die fsinBerechnung (gerade weil sie mit weniger Präzision durchgeführt wird). Beachten Sie, dass der x87 bekanntermaßen eine etwas geringere tatsächliche Genauigkeit aufweist als die angekündigten 79 Bit.
Thomas Pornin
1
In der Tat verwenden sowohl 32-Bit- als auch 64-Bit-Implementierungen von sin () in den msvc-Laufzeitbibliotheken nicht den FSIN-Befehl. Tatsächlich liefern sie unterschiedliche Ergebnisse, zum Beispiel sin (0.70444454416678126). Dies führt in einem 32-Bit-Programm zu 0,64761068800896837 (rechts mit einer Toleranz von 0,5 * (eps / 2)) und in einem 64-Bit-Programm zu 0,64761068800896848 (falsch).
e.tadeu
9

Chebyshev-Polynome sind, wie in einer anderen Antwort erwähnt, die Polynome, bei denen der größte Unterschied zwischen der Funktion und dem Polynom so gering wie möglich ist. Das ist ein ausgezeichneter Start.

In einigen Fällen ist der maximale Fehler nicht das, woran Sie interessiert sind, sondern der maximale relative Fehler. Zum Beispiel sollte für die Sinusfunktion der Fehler in der Nähe von x = 0 viel kleiner sein als für größere Werte; Sie möchten einen kleinen relativen Fehler. Sie würden also das Chebyshev-Polynom für sin x / x berechnen und dieses Polynom mit x multiplizieren.

Als nächstes müssen Sie herausfinden, wie das Polynom ausgewertet wird. Sie möchten es so auswerten, dass die Zwischenwerte klein sind und daher Rundungsfehler klein sind. Andernfalls können die Rundungsfehler viel größer werden als die Fehler im Polynom. Und bei Funktionen wie der Sinusfunktion ist es möglich, dass das Ergebnis, das Sie für sin x berechnen, größer ist als das Ergebnis für sin y, selbst wenn x <y ist, wenn Sie nachlässig sind. Daher ist eine sorgfältige Auswahl der Berechnungsreihenfolge und die Berechnung der Obergrenzen für den Rundungsfehler erforderlich.

Zum Beispiel sin x = x - x ^ 3/6 + x ^ 5/120 - x ^ 7/5040 ... Wenn Sie naiv sin x = x * (1 - x ^ 2/6 + x ^ 4 / berechnen) 120 - x ^ 6/5040 ...), dann nimmt diese Funktion in Klammern ab, und es kommt vor , dass wenn y die nächstgrößere Zahl zu x ist, sin y manchmal kleiner als sin x ist. Berechnen Sie stattdessen sin x = x - x ^ 3 * (1/6 - x ^ 2/120 + x ^ 4/5040 ...), wo dies nicht passieren kann.

Bei der Berechnung von Chebyshev-Polynomen müssen Sie beispielsweise die Koeffizienten normalerweise auf doppelte Genauigkeit runden. Während ein Chebyshev-Polynom optimal ist, ist das Chebyshev-Polynom mit auf doppelte Genauigkeit gerundeten Koeffizienten nicht das optimale Polynom mit Koeffizienten mit doppelter Genauigkeit!

Zum Beispiel für sin (x), wo Sie Koeffizienten für x, x ^ 3, x ^ 5, x ^ 7 usw. benötigen, gehen Sie wie folgt vor: Berechnen Sie die beste Näherung von sin x mit einem Polynom (ax + bx ^ 3 +) cx ^ 5 + dx ^ 7) mit höherer als doppelter Genauigkeit, dann runde a auf doppelte Genauigkeit, was A ergibt. Der Unterschied zwischen a und A wäre ziemlich groß. Berechnen Sie nun die beste Näherung von (sin x - Ax) mit einem Polynom (bx ^ 3 + cx ^ 5 + dx ^ 7). Sie erhalten unterschiedliche Koeffizienten, weil sie sich an die Differenz zwischen a und A anpassen. Runde b mit doppelter Genauigkeit B. Dann approximieren Sie (sin x - Ax - Bx ^ 3) mit einem Polynom cx ^ 5 + dx ^ 7 und so weiter. Sie erhalten ein Polynom, das fast so gut ist wie das ursprüngliche Chebyshev-Polynom, aber viel besser als Chebyshev, das auf doppelte Genauigkeit gerundet ist.

Als nächstes sollten Sie die Rundungsfehler bei der Wahl des Polynoms berücksichtigen. Sie haben ein Polynom mit minimalem Fehler im Polynom gefunden, das Rundungsfehler ignoriert, aber Sie möchten Polynom plus Rundungsfehler optimieren. Sobald Sie das Chebyshev-Polynom haben, können Sie die Grenzen für den Rundungsfehler berechnen. Angenommen, f (x) ist Ihre Funktion, P (x) ist das Polynom und E (x) ist der Rundungsfehler. Sie möchten | nicht optimieren f (x) - P (x) | möchten Sie optimieren | f (x) - P (x) +/- E (x) |. Sie erhalten ein etwas anderes Polynom, das versucht, die Polynomfehler niedrig zu halten, wenn der Rundungsfehler groß ist, und die Polynomfehler ein wenig lockert, wenn der Rundungsfehler klein ist.

All dies führt leicht zu Rundungsfehlern von höchstens dem 0,55-fachen des letzten Bits, wobei +, -, *, / Rundungsfehler von höchstens dem 0,50-fachen des letzten Bits aufweisen.

gnasher729
quelle
1
Dies ist eine schöne Erklärung, wie man kann sin berechnen (x) effizient, aber es scheint nicht wirklich die OP Frage zu beantworten, die speziell darüber, wie häufig C - Bibliotheken / Compiler kann es berechnet wird .
Ilmari Karonen
Chebyshev-Polynome minimieren den maximalen Absolutwert über ein Intervall, minimieren jedoch nicht den größten Unterschied zwischen einer Zielfunktion und dem Polynom. Minimax-Polynome machen das.
Eric Postpischil
9

In Bezug auf trigonometrische Funktionen wie sin(), cos(), tan()nach 5 Jahren hat es nicht erwähnt worden, nachdem einen wichtigen Aspekt von hohen Qualität trigonometrischen Funktionen: Reichweite Reduzierung .

Ein früher Schritt bei einer dieser Funktionen besteht darin, den Winkel im Bogenmaß auf einen Bereich von 2 * π zu reduzieren. Aber π ist irrational, so dass einfache Reduktionen wie das x = remainder(x, 2*M_PI)Einführen eines Fehlers M_PIoder die Maschine pi eine Annäherung an π sind. Also, wie geht das?x = remainder(x, 2*π) ?

Frühe Bibliotheken verwendeten erweiterte Präzision oder gestaltete Programmierung, um qualitativ hochwertige Ergebnisse zu erzielen, jedoch immer noch über einen begrenzten Bereich von double. Wenn ein großer Wert wie angefordert wurde sin(pow(2,30)), waren die Ergebnisse bedeutungslos oder 0.0und möglicherweise mit einem Fehlerflag, das auf einen TLOSSvollständigen Genauigkeitsverlust oder einen PLOSSteilweisen Genauigkeitsverlust gesetzt war.

Eine gute Entfernungsreduzierung großer Werte auf ein Intervall wie -π bis π ist ein herausforderndes Problem, das den Herausforderungen der grundlegenden Triggerfunktion wie sin()sich selbst Konkurrenz macht .

Ein guter Bericht ist die Argumentreduktion für große Argumente: Gut bis zum letzten Bit (1992). Es behandelt das Problem gut: Erörtert die Notwendigkeit und den Stand der Dinge auf verschiedenen Plattformen (SPARC, PC, HP, 30+ andere) und bietet einen Lösungsalgorithmus, der Qualitätsergebnisse für alle double von -DBL_MAXbis liefert DBL_MAX.


Wenn die ursprünglichen Argumente in Grad angegeben sind und dennoch einen großen Wert haben können, verwenden Sie sie fmod()zuerst, um die Genauigkeit zu verbessern. Ein Gut fmod()führt zu keinem Fehler und bietet somit eine hervorragende Reichweitenreduzierung.

// sin(degrees2radians(x))
sin(degrees2radians(fmod(x, 360.0))); // -360.0 < fmod(x,360) < +360.0

Verschiedene Triggeridentitäten und remquo()bieten noch mehr Verbesserungen. Beispiel: sind ()

chux - Monica wieder einsetzen
quelle
6

Die eigentliche Implementierung der Bibliotheksfunktionen hängt vom jeweiligen Compiler und / oder Bibliotheksanbieter ab. Ob es sich um Hardware oder Software handelt, ob es sich um eine Taylor-Erweiterung handelt oder nicht usw., ist unterschiedlich.

Mir ist klar, dass das absolut keine Hilfe ist.

John Bode
quelle
5

Sie werden normalerweise in Software implementiert und verwenden in den meisten Fällen nicht die entsprechenden Hardware-Aufrufe (dh Assemblierungsaufrufe). Wie Jason jedoch betonte, sind diese implementierungsspezifisch.

Beachten Sie, dass diese Softwareroutinen nicht Teil der Compilerquellen sind, sondern sich in der entsprechenden Bibliothek wie clib oder glibc für den GNU-Compiler befinden. Sehen http://www.gnu.org/software/libc/manual/html_mono/libc.html#Trig-Functions

Wenn Sie mehr Kontrolle wünschen, sollten Sie sorgfältig abwägen, was Sie genau benötigen. Einige der typischen Methoden sind die Interpolation von Nachschlagetabellen, der Assemblyaufruf (der häufig langsam ist) oder andere Approximationsschemata wie Newton-Raphson für Quadratwurzeln.

Mnemosyn
quelle
5

Wenn Sie eine Implementierung in Software und nicht in Hardware wünschen, suchen Sie in Kapitel 5 der Numerischen Rezepte nach einer endgültigen Antwort auf diese Frage . Meine Kopie befindet sich in einer Box, daher kann ich keine Details angeben, aber die Kurzversion (wenn ich mich an dieses Recht erinnere) ist, dass Sie tan(theta/2)als Ihre primitive Operation nehmen und die anderen von dort aus berechnen. Die Berechnung erfolgt mit einer Seriennäherung, die jedoch viel schneller konvergiert als eine Taylor-Serie.

Entschuldigung, ich kann mich nicht mehr erinnern, ohne das Buch in die Hand zu nehmen.

Norman Ramsey
quelle
5

Es gibt nichts Schöneres, als die Quelle zu treffen und zu sehen, wie jemand es tatsächlich in einer allgemein verwendeten Bibliothek getan hat. Schauen wir uns insbesondere eine Implementierung der C-Bibliothek an. Ich habe mich für uLibC entschieden.

Hier ist die Sündenfunktion:

http://git.uclibc.org/uClibc/tree/libm/s_sin.c

Dies sieht so aus, als würde es einige Sonderfälle behandeln und dann eine Argumentreduktion durchführen, um die Eingabe dem Bereich [-pi / 4, pi / 4] zuzuordnen (Aufteilung des Arguments in zwei Teile, einen großen Teil und einen Schwanz). vor dem Anruf

http://git.uclibc.org/uClibc/tree/libm/k_sin.c

die dann auf diese beiden Teile arbeitet. Wenn es keinen Schwanz gibt, wird eine ungefähre Antwort unter Verwendung eines Polynoms vom Grad 13 erzeugt. Wenn es einen Schwanz gibt, erhalten Sie eine kleine Korrekturaddition, die auf dem Prinzip basiert, dasssin(x+y) = sin(x) + sin'(x')y

Moschops
quelle
4

Wann immer eine solche Funktion bewertet wird, gibt es auf einer bestimmten Ebene höchstwahrscheinlich entweder:

  • Eine Wertetabelle, die interpoliert wird (für schnelle, ungenaue Anwendungen - z. B. Computergrafik)
  • Die Bewertung einer Reihe, die zum gewünschten Wert konvergiert - wahrscheinlich keine Taylor-Reihe, eher etwas, das auf einer ausgefallenen Quadratur wie Clenshaw-Curtis basiert.

Wenn keine Hardwareunterstützung vorhanden ist, verwendet der Compiler wahrscheinlich die letztere Methode und gibt nur Assembler-Code (ohne Debug-Symbole) aus, anstatt eine AC-Bibliothek zu verwenden. Dies macht es für Sie schwierig, den tatsächlichen Code in Ihrem Debugger aufzuspüren.

James
quelle
4

Wie viele Leute betonten, ist es implementierungsabhängig. Soweit ich Ihre Frage verstehe, waren Sie an einer echten Softwareimplementierung von mathematischen Funktionen interessiert , haben es aber einfach nicht geschafft, eine zu finden. Wenn dies der Fall ist, sind Sie hier:

  • Laden Sie den Glibc-Quellcode von http://ftp.gnu.org/gnu/glibc/ herunter.
  • Sehen Sie sich die Datei dosincos.cim entpackten Ordner glibc root \ sysdeps \ ieee754 \ dbl-64 an
  • Ebenso finden Sie Implementierungen des Restes der Mathematikbibliothek. Suchen Sie einfach nach der Datei mit dem entsprechenden Namen

Sie können sich auch die Dateien mit der .tblErweiterung ansehen. Ihr Inhalt besteht lediglich aus riesigen Tabellen mit vorberechneten Werten verschiedener Funktionen in binärer Form. Aus diesem Grund ist die Implementierung so schnell: Anstatt alle Koeffizienten der von ihnen verwendeten Serien zu berechnen, führen sie einfach eine schnelle Suche durch, die viel schneller ist. Übrigens verwenden sie Tailor-Reihen, um Sinus und Cosinus zu berechnen.

Ich hoffe das hilft.

Igor Korkhov
quelle
4

Ich werde versuchen, den Fall sin()in einem C-Programm zu beantworten , das mit dem C-Compiler von GCC auf einem aktuellen x86-Prozessor kompiliert wurde (sagen wir ein Intel Core 2 Duo).

In der C - Sprache umfasst die Standard - C - Bibliothek gemeinsame mathematische Funktionen, nicht in der Sprache enthielt selbst (zB pow, sinund cosfür Strom, Sinus, Kosinus und jeweils). Die Überschriften davon sind in math.h enthalten .

Auf einem GNU / Linux-System werden diese Bibliotheksfunktionen von glibc (GNU libc oder GNU C Library) bereitgestellt. Der GCC-Compiler möchte jedoch, dass Sie mithilfe des Compiler-Flags eine Verknüpfung zur Mathematikbibliothek ( libm.so) herstellen -lm, um die Verwendung dieser Mathematikfunktionen zu ermöglichen. Ich bin mir nicht sicher, warum es nicht Teil der Standard-C-Bibliothek ist. Dies wäre eine Softwareversion der Gleitkommafunktionen oder "Soft-Float".

Nebenbei: Der Grund für die Trennung der mathematischen Funktionen ist historisch und sollte meines Wissens lediglich die Größe ausführbarer Programme in sehr alten Unix-Systemen reduzieren , möglicherweise bevor gemeinsam genutzte Bibliotheken verfügbar waren.

Jetzt kann der Compiler die Standardfunktion der C-Bibliothek sin()(bereitgestellt von libm.so) optimieren, die durch einen Aufruf einer nativen Anweisung für die integrierte sin () -Funktion Ihrer CPU / FPU ersetzt wird, die als FPU-Anweisung ( FSINfür x86 / x87) vorhanden ist Neuere Prozessoren wie die Core 2-Serie (dies ist bereits im i486DX richtig). Dies würde von Optimierungsflags abhängen, die an den gcc-Compiler übergeben werden. Wenn der Compiler angewiesen würde, Code zu schreiben, der auf einem i386 oder einem neueren Prozessor ausgeführt wird, würde er eine solche Optimierung nicht vornehmen. Das -mcpu=486Flag würde den Compiler darüber informieren, dass eine solche Optimierung sicher ist.

Nun , wenn das Programm die Software - Version der sin () Funktion ausgeführt wird , würde sie es so auf der Grundlage eines CORDIC (Coordinate Rotation Digital Computer) oder BKM - Algorithmus oder mehr wahrscheinlich eine Tabelle oder Potenzreihenberechnung , die jetzt häufig verwendet wird , zu berechnen solche transzendentalen Funktionen. [Src: http://en.wikipedia.org/wiki/Cordic#Application]

Jede neuere (seit ca. 2,9x) Version von gcc bietet auch eine integrierte Version von sin, __builtin_sin()mit der der Standardaufruf der C-Bibliotheksversion als Optimierung ersetzt wird.

Ich bin mir sicher, dass das so klar wie Schlamm ist, aber hoffentlich gibt es Ihnen mehr Informationen als erwartet und viele Absprungpunkte, um selbst mehr zu lernen.

mctylr
quelle
3

Wenn Sie sich die tatsächliche GNU-Implementierung dieser Funktionen in C ansehen möchten, sehen Sie sich die neueste Version von glibc an. Siehe die GNU C-Bibliothek .

Chris Tonkinson
quelle
3

Verwenden Sie keine Taylor-Serien. Chebyshev-Polynome sind sowohl schneller als auch genauer, wie einige Leute oben hervorgehoben haben. Hier ist eine Implementierung (ursprünglich aus dem ZX Spectrum ROM): https://albertveli.wordpress.com/2015/01/10/zx-sine/

Albert Veli
quelle
2
Dies scheint die gestellte Frage nicht wirklich zu beantworten. Die OP ist zu fragen , wie trigonometrische Funktionen werden berechnet durch gemeinsame C - Compiler / Bibliotheken (und ich bin ziemlich sicher , dass ZX Spectrum nicht qualifiziert), nicht , wie sie sollten berechnet werden. Dies könnte jedoch ein nützlicher Kommentar zu einigen der früheren Antworten gewesen sein.
Ilmari Karonen
1
Ah, du hast recht. Es hätte ein Kommentar und keine Antwort sein sollen. Ich habe SO schon eine Weile nicht mehr benutzt und vergessen, wie das System funktioniert. Wie auch immer, ich denke, die Spectrum-Implementierung ist relevant, da sie eine sehr langsame CPU hatte und die Geschwindigkeit von entscheidender Bedeutung war. Der beste Algorithmus ist dann sicherlich noch ziemlich gut, daher wäre es für C-Bibliotheken eine gute Idee, Triggerfunktionen unter Verwendung von Chebyshev-Polynomen zu implementieren.
Albert Veli
2

Das Berechnen von Sinus / Cosinus / Tangens ist mit Code mithilfe der Taylor-Reihe sehr einfach. Das Schreiben selbst dauert etwa 5 Sekunden.

Der gesamte Prozess kann hier mit folgender Gleichung zusammengefasst werden:

Sünde und Kostenerweiterung

Hier sind einige Routinen, die ich für C geschrieben habe:

double _pow(double a, double b) {
    double c = 1;
    for (int i=0; i<b; i++)
        c *= a;
    return c;
}

double _fact(double x) {
    double ret = 1;
    for (int i=1; i<=x; i++) 
        ret *= i;
    return ret;
}

double _sin(double x) {
    double y = x;
    double s = -1;
    for (int i=3; i<=100; i+=2) {
        y+=s*(_pow(x,i)/_fact(i));
        s *= -1;
    }  
    return y;
}
double _cos(double x) {
    double y = 1;
    double s = -1;
    for (int i=2; i<=100; i+=2) {
        y+=s*(_pow(x,i)/_fact(i));
        s *= -1;
    }  
    return y;
}
double _tan(double x) {
     return (_sin(x)/_cos(x));  
}
user1432532
quelle
4
Dies ist eine ziemlich schlechte Implementierung, da nicht verwendet wird, dass die aufeinanderfolgenden Terme der Sinus- und Cosinus-Reihe sehr einfache Quotienten haben. Dies bedeutet, dass man die Anzahl der Multiplikationen und Divisionen von O (n ^ 2) hier auf O (n) reduzieren kann. Weitere Reduzierungen werden durch Halbieren und Quadrieren erreicht, wie dies beispielsweise in der Mathematikbibliothek bc (POSIX Multiprecision Calculator) der Fall ist.
Lutz Lehmann
2
Es scheint auch nicht die gestellte Frage zu beantworten; Das OP fragt, wie Triggerfunktionen von gängigen C-Compilern / Bibliotheken berechnet werden, nicht für benutzerdefinierte Neuimplementierungen.
Ilmari Karonen
2
Ich denke, es ist eine gute Antwort, da es den Geist der Frage beantwortet, die (und ich kann natürlich nur raten) die Neugier auf eine ansonsten "Black Box" -Funktion wie sin () weckt. Dies ist die einzige Antwort, die es einem ermöglicht, schnell zu verstehen, was passiert, indem man es in wenigen Sekunden beschönigt, anstatt einen optimierten C-Quellcode zu lesen.
Mike M
Tatsächlich verwenden Bibliotheken die viel optimierte Version, indem sie erkennen, dass Sie, sobald Sie einen Begriff haben, den nächsten Begriff erhalten können, indem Sie einige Werte multiplizieren. Siehe ein Beispiel in Blindys Antwort . Sie berechnen immer wieder die Leistung und die Fakultäten, was weitaus langsamer ist
phuclv
0

Verbesserte Version des Codes aus Blindys Antwort

#define EPSILON .0000000000001
// this is smallest effective threshold, at least on my OS (WSL ubuntu 18)
// possibly because factorial part turns 0 at some point
// and it happens faster then series element turns 0;
// validation was made against sin() from <math.h>
double ft_sin(double x)
{
    int k = 2;
    double r = x;
    double acc = 1;
    double den = 1;
    double num = x;

//  precision drops rapidly when x is not close to 0
//  so move x to 0 as close as possible
    while (x > PI)
        x -= PI;
    while (x < -PI)
        x += PI;
    if (x > PI / 2)
        return (ft_sin(PI - x));
    if (x < -PI / 2)
        return (ft_sin(-PI - x));
//  not using fabs for performance reasons
    while (acc > EPSILON || acc < -EPSILON)
    {
        num *= -x * x;
        den *= k * (k + 1);
        acc = num / den;
        r += acc;
        k += 2;
    }
    return (r);
}
Rugnar
quelle
0

Das Wesentliche dafür liegt in diesem Auszug aus der Angewandten Numerischen Analyse von Gerald Wheatley:

Wenn Ihr Softwareprogramm den Computer auffordert, einen Wert von Geben Sie hier die Bildbeschreibung einoder zu erhalten Geben Sie hier die Bildbeschreibung ein, haben Sie sich gefragt, wie er die Werte erhalten kann, wenn die leistungsstärksten Funktionen, die er berechnen kann, Polynome sind? Diese werden nicht in Tabellen nachgeschlagen und interpoliert! Vielmehr approximiert der Computer jede andere Funktion als Polynome von einem Polynom, das so zugeschnitten ist, dass die Werte sehr genau angegeben werden.

Einige Punkte, die oben erwähnt werden sollten, sind, dass einige Algorithmen tatsächlich aus einer Tabelle interpolieren, wenn auch nur für die ersten paar Iterationen. Beachten Sie auch, wie erwähnt wird, dass Computer Approximationspolynome verwenden, ohne anzugeben, welcher Typ von Approximationspolynomen verwendet wird. Wie andere im Thread hervorgehoben haben, sind Chebyshev-Polynome in diesem Fall effizienter als Taylor-Polynome.

Dean P.
quelle
-1

wenn du willst sindann

 __asm__ __volatile__("fsin" : "=t"(vsin) : "0"(xrads));

wenn du willst cosdann

 __asm__ __volatile__("fcos" : "=t"(vcos) : "0"(xrads));

wenn du willst sqrtdann

 __asm__ __volatile__("fsqrt" : "=t"(vsqrt) : "0"(value));

Warum also ungenauen Code verwenden, wenn die Maschinenanweisungen ausreichen?

user80998
quelle