Wie berechnet C sin () und andere mathematische Funktionen?

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).

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

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 .

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.

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;

}

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.

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 )

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.

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.

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 ()

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.

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.

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.

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

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.

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.

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.

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 .

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/

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:

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));  
}

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);
}

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 oder zu erhalten , 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.

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?