Motivation
In der Spieleprogrammierung werden häufig trigonometrische Funktionen, insbesondere Kosinus und Sinus verwendet. Vor allem im 3D Bereich bauen viele Berechnungen auf trigonometrischen Funktionen auf. Demnach könnten diese Funktionen nicht schnell genug sein. Um mehr Leistung aus dem Programm herauszuholen greifen viele auf Lookup Tabellen zurück. Das ist keine schlechte Lösung, wenn das Ergebnis nicht sehr genau sein soll. Möchte man jedoch mindestens eine Genauigkeit von einer Nachkommastelle haben, so braucht man einen 3600 Werte großen float-Array und je größer das Array, desto kleiner wird der Geschwindigkeitsgewinn gegenüber der Standardimplementierung sein.
Es gibt einen anderen Weg: die Approximation(deutsch: Annäherung) einer Funktion. Durch diese mathematische Methode kann man eine Funktion gewinnen, die eine bestimmte Funktion in einem vordefinierten Intervall annähernd beschreibt.
In unserem Fall würden wir eine Funktion gewinnen, die den Sinus im Intervall von –PI bis + PI gut beschreibt. Ein größeres Intervall wird nicht gebraucht, da der Sinus 2Pi-periodisch ist.
Sinus
Wie genau man auf die Näherungsfunktion kommt, werde ich hier nicht behandeln.
Wer sich dafür interessiert, sollte nach Begriffen Modellierung, Approximation, Potenzreihenentwicklung, Taylor’sche Formel oder Maclaurin’sche Formel suchen. Das Thema Modellierung wird normalerweise auf einem Gymnasium in der zwölften Klasse behandelt.
Stattdessen nehmen wir gleich die fertige Potenzreihe und setzen sie in C++ um.
[math] sin(x) = \sum_{n=1}^{\infty} (-1)^n \frac{x^{2n+1}}{(2n+1)!} [/math]
Diese Funktion modelliert Sinus um den Nullpunkt herum. Je mehr Potenzreihenglieder aufsummiert werden, desto besser beschreibt diese Funktion den Sinus. Dies kann man sehr schön an folgenden Funktionsplots sehen.
Legende:
Orange = sin(x)
Schwarz = [math]x – \frac{x^3}{3!} + \frac{x^5}{5!} [/math]
Wie man sieht, wird der Sinus im Koordinatenursprung gut beschrieben, nicht aber an den PI-Grenzen. Wir addieren einen weiteren Term dazu.
Legende:
Orange = sin(x)
Schwarz = [math]x – \frac{x^3}{3!} + \frac{x^5}{5!} – \frac{x^7}{7!}[/math]
Schon besser, aber die Ungenauigkeit doch noch zu groß. Wir erweitern weiter.
Legende:
Orange = sin(x)
Schwarz = [math]x – \frac{x^3}{3!} + \frac{x^5}{5!} – \frac{x^7}{7!} + \frac{x^9}{9!}[/math]
Diese Näherungsfunktion beschreibt den Sinus im Intervall von –PI bis PI schon ziemlich genau, was man daran erkennt, dass in diesem Intervall der orangenfarbene Graph komplett überdeckt wird.
Hier noch mal ein größerer Ausschnitt.
Wie gut diese Modellierung wirklich ist, schauen wir uns später an konkreten Zahlen an.
Demnach lautet unsere Funktion [math]sin(x) = x – \frac{x^3}{3!} + \frac{x^5}{5!} – \frac{x^7}{7!} + \frac{x^9}{9!}[/math]. Bevor wir sie aber in C++ umsetzen, müssen wir bedenken, dass sie nur im Intervall von –PI bis PI gilt. Wird ein Argument übergeben, der nicht in diesem Intervall liegt, so liefert die Funktion falsche Werte, was man an dem letzten Graphen sehr gut zu erkennen ist.
Die Sinusfunktion ist genauso wie der Kosinus [math]2\pi[/math] periodisch. Das heißt, wenn man zum Funktionsargument z-mal [math]2\pi[/math] dazu addiert oder von dem Funktionsargument z-mal [math]2\pi[/math] abzieht, wird sich der Wert der Funktion sin(x) nicht verändert (z ist eine Ganzzahl).
Man könnte auch schreiben:
[math] sin(x) = sin(x + z * 2\pi) [/math] oder [math] sin(x) = sin(x -z * 2\pi) [/math]
Für uns bedeutet das, dass ein Funktionsargument, der nicht in dem Intervall von –PI bis PI liegt, auf dieses Intervall zurückgerechnet werden kann ohne dass das Ergebnis der Funktion verfälscht wird.
Dazu einige Beispiele (Argument => Umrechnung => neuer Wert ):
85,222 => 85,222-14*2*PI => ~ -2,743
12,982 => 12,982-2*PI => ~ 0,416
-29, 314 => -29, 314+5*2*PI => ~ 2.102
Wie kommt man auf die z für die Umrechnung?
Es wird zwischen zwei Fällen unterschieden:
- x ist größer als PI
- int T = x + PI
- T wird durch 2*PI geteilt, das Ergebnis ohne Nachkommastelle(!) ist z
- x ist kleiner als PI
- int T = x – PI
- T wird durch 2*PI geteilt, das Ergebnis ohne Nachkommastelle(!) ist z
Als nächstes wird das Argument durch die Formel [math]x = x-z*2\pi[/math] auf das Intervall –PI bis PI zurückgerechnet. Ob das z positiv oder negativ ist, spielt keine Rolle, den wenn z negativ ist wird der Term [math]x = x-(-z)*2\pi[/math] zu [math]x = x+z*2\pi[/math].
Nachdem das Argument im richtigen Intervall liegt, können wir unsere Taylor-Reihe benutzen. Doch zuerst schreiben wir die Taylor-Reihe [math]x – \frac{x^3}{3!} + \frac{x^5}{5!} – \frac{x^7}{7!} + \frac{x^9}{9!}[/math] so um, dass die Fakultät- und Division-Operationen verschwinden, da diese sehr langsamer sind.
Wir bekommen
[math] sin(x) = x – 0.1666667*x*x*x + 0.0083333*x*x*x*x*x – 0.0001984*x*x*x*x*x*x*x + 0.0000027*x*x*x*x*x*x*x*x*x[/math]
und durch Umformung
[math]sin(x) = x*(1 – x*x*(0.1666667 – x*x*(0.0083333 – x*x*(0.0001984 – 0.0000027*x*x)))) [/math]
An dieser Stelle kommt unsere erste Sinus-Funktion in C++.
const float PI = 3.14159265f;
const float PI2 = PI*2;
// schnelle Sinus Implementierung
inline float _fastcall fastSin(float x)
{
if(PI < x)
{
x = x-static_cast<int>((x+PI)/(PI2))*PI2;
}
else if(x < -PI)
{
x = x-static_cast<int>((x-PI)/(PI2))*PI2;
}
return x*(1 - x*x*(0.16666667f - x*x*(0.00833333f - x*x*(0.0001984f - x*x*0.0000027f))));
}
Ich habe einen kleinen Geschwindigkeitstest durchgeführt, bei dem std::sin(x) und unsere Implementierung 10-Millionen Mal aufgerufen wurde. Das Ergebnis ist ein Mittelwert aus 5 Durchläufen: std::sin(): 1174ms, fastSin(): 656ms. Selbstverständlich können diese Werte nur qualitativ betrachtet werden.
Unsere Implementierung weist eine Genauigkeit von 2-6 Nachkommastellen auf, die für Spieleprogrammierung akzeptabel sein sollte. Wie man erkennt, ist die Ungenauigkeit an den [math]\pi[/math]-Grenzen am höchsten.
| Argument | std::sin | fastSin | Abweichung |
|---|---|---|---|
| 0 | 0.000000 | 0.000000 | 0 |
| ¼ PI | 0.707107 | 0.707107 | 0 |
| ½ PI | 1.000000 | 1.000000 | 0 |
| ¾ PI | 0.707107 | 0.707287 | 0.000180 |
| PI | 0.000000 | 0.005301 | 0.005301 |
Um eine höhere Genauigkeit zu erzielen, muss man das Polynom um weitere Glieder erweitern: [math]x – \frac{x^3}{3!} + \frac{x^5}{5!} – \frac{x^7}{7!} + \frac{x^9}{9!} – \frac{x^11}{11!}[/math].
Als nächsten Schritt werden wir die Methode in Inline Assembler schreiben, wobei dabei Intels SSE Technik verwendet wird. SSE wird ab Pentium 3 und Athon XP unterstützt. Ich werde hier nicht erklären, wie SSE funktioniert, weil es den Rahmen dieses Artikels sprengen würde. Es handelt sich um 1:1 Konvertierung der C++ Version der Funktion.
// schnelle Sinus Implementierung mit Inline Assembler
inline float _fastcall fastSinSSE(float x)
{
static const float _0_16666 = -1/6.0f;
static const float _0_00833 = 1/120.0f;
static const float _0_00019 = -1/5040.0f;
static const float _0_0000027 = 1/362880.0f;
_asm
{
movss xmm0, x
comiss xmm0, PI // if (x > PI)
jbe SHORT ln1
addss xmm0, PI
jmp SHORT ln2
ln1:
movss xmm1, _PI
comiss xmm1, xmm0 // if ( x < -PI)
jbe SHORT ln3
subss xmm0, PI
ln2:
divss xmm0, PI2
cvttss2si ecx, xmm0 // = z ohne rundung!
cvtsi2ss xmm1, ecx // z wieder in float umwandeln
movss xmm0, x
mulss xmm1, PI2
subss xmm0, xmm1
ln3:
// xmm0 beinhaltet immer das aktuelle ergebnis
// xmm1 beinhaltet zu jeder zeit x^2
// xmm2 beinhaltet zu jeder zeit x
// x + x*x*x*(_0_16666 + x*x*(_0_00833 + x*x*(_0_00019f + _0_0000027*x*x))));
movss xmm2, xmm0 // x in xmm0 legen
movss xmm1, xmm0 // x in xmm1 legen
mulss xmm1, xmm1 // xmm1 = x*x
movss xmm0, _0_0000027
mulss xmm0, xmm1 // xmm0 = _0_0000027*x*x
addss xmm0, _0_00019
mulss xmm0, xmm1 // xmm0 = x*x*(_0_00019f + _0_0000027*x*x)
addss xmm0, _0_00833
mulss xmm0, xmm1 // xmm0 = x*x*(_0_00833 + x*x*(_0_00019f + _0_0000027*x*x))
addss xmm0, _0_16666
mulss xmm0, xmm1 // xmm0 = x*x*(_0_16666 + x*x*(_0_00833 + x*x*(_0_00019 + _0_0000027*x*x)))
mulss xmm0, xmm2
addss xmm0, xmm2
movss x, xmm0 // das Ergebnis in x rüberschieben
}
return x;
}
Durch Assembler wurde unsere Funktion noch ein wenig schneller: 10-Millionen Berechnungen in 493ms. Also rund zwei Mal schneller als die Standardimplementierung. Der Preis dafür ist die Ungenauigkeit.
Als nächstes schauen wir uns im Schnelldurchlauf den Kosinus an.
Kosinus
Das Schreiben der Kosinus-Funktion funktioniert analog zur Sinus-Funktion. Man kann aber auch die Gesetzmäßigkeit [math]cos(x)=sin( \frac{\pi}{2} – x)[/math] ausnutzen. Man schreibt einfach eine neue Funktion, die unsere Sinus-Funktion mit dem Argument [math]\frac{\pi}{2} – x[/math] aufruft:
float fastCos(float x)
{
return fastSin(PI/2 – x);
}
Ja, so schnell kann es gehen. Der Nachteil von dieser Funktion ist, dass das Aufrufen der Funktion länger dauert als die Berechnung selbst, also machen wir es analog zum Sinus – Nährungsfunktion in Inline Assembler implementieren.
Die Taylor-Reihe für Kosinus:
[math] sin(x) = \sum_{n=1}^{\infty} (-1)^n \frac{x^{2n}}{(2n)!} [/math]
Wir nehmen die ersten sechs Glieder.
Legende:
Orange = cos(x)
Schwarz = [math]1 – \frac{x^2}{2!} + \frac{x^4}{4!} – \frac{x^6}{6!}+ \frac{x^8}{8!}+ \frac{x^10}{10!}[/math]
Schon besser, aber die Ungenauigkeit doch noch zu groß. Wir erweitern weiter.
Wir schreiben das Polynom in leistungsstärkere Variante um:
[math]cos(x) = 1-x*x*(0.5-x*x*(0.04166667-x*x*(0.00138889-x*x*(0.00002480-x*x*0.000000275))))[/math]
// schnelle Kosinus Implementierung
inline float _fastcall fastCos(float x)
{
if(PI < x)
{
x = x-static_cast<int>((x+PI)/(PI2))*PI2;
}
else if(x < -PI)
{
x = x-static_cast<int>((x-PI)/(PI2))*PI2;
}
return 1-x*x*(0.5f-x*x*(0.04166667f-x*x*(0.00138889f-x*x*(0.00002480f-x*x*0.000000275f))));
}
Die Genauigkeit von Kosinus ist besser als die von Sinus und liegt bei 3-6 Nachkommastellen. Das liegt daran, dass wir mehr Reihenglieder benutzt haben.
| Argument | std::cos | fastCos | Abweichung |
|---|---|---|---|
| 0 | 1.000000 | 1.000000 | 0 |
| ¼ PI | 0.707107 | 0.707107 | 0 |
| ½ PI | 1.000000 | 1.000000 | 0 |
| ¾ PI | -0.707107 | -0.707287 | 0.000058 |
| PI | -1.000000 | -1.001790 | 0.001790 |
Beim Geschwindigkeitstest (10 Millionen Mal aufgerufen) braucht die Standardmplementierung von Kosinus 2305ms und fastCos 1337ms.
Die Inline Assembler Version schafft 10 Millionen Aufrufe in 996ms.
// schnelle Kosinus Implementierung mit Inline Assembler
inline float _fastcall fastCosSSE(float x)
{
static const float _0_5 = -1/2.0f;
static const float _0_0416 = 1/24.0f;
static const float _0_001387 = -1/720.0f;
static const float _0_0000248 = 1/40320.0f;
static const float _0_000000275 = -1/3629000.0f;
static const float _1 = 1.0f;
// Wenn x groeßer oder kleiner als PI,
// dann wird x auf das Intervall -PI bis PI zurückgerechnet
_asm
{
movss xmm0, x
comiss xmm0, PI // if (x > PI)
jbe SHORT ln1
addss xmm0, PI
jmp SHORT ln2
ln1:
movss xmm1, _PI
comiss xmm1, xmm0 // if ( x < -PI)
jbe SHORT ln3
subss xmm0, PI
ln2:
divss xmm0, PI2
cvttss2si ecx, xmm0 // = z ohne rundung!
cvtsi2ss xmm1, ecx // z wieder in float umwandeln
movss xmm0, x
mulss xmm1, PI2
subss xmm0, xmm1
ln3:
// xmm0 beinhaltet immer das aktuelle ergebnis
// xmm1 beinhaltet zu jeder zeit x^2
// 1-x*x*(0.5f-x*x*(0.04166667f-x*x*(0.00138889f-x*x*(0.00002480f-x*x*0.000000275f))));
movss xmm1, xmm0 // x in xmm1 legen
mulss xmm1, xmm0 // xmm1 = x*x
movss xmm0, _0_000000275
mulss xmm0, xmm1 // xmm0 = _0_000000275*x*x
addss xmm0, _0_0000248
mulss xmm0, xmm1 // xmm0 = x*x*(_0_0000248 + _0_000000275*x*x)
addss xmm0, _0_001387
mulss xmm0, xmm1 // xmm0 = x*x*(_0_001387 + x*x*(_0_0000248 + _0_000000275*x*x))
addss xmm0, _0_0416
mulss xmm0, xmm1 // xmm0 = x*x*(_0_0416 + x*x*(_0_001387 + x*x*(_0_0000248 + _0_000000275*x*x)))
addss xmm0, _0_5
mulss xmm0, xmm1 // xmm0 = x*x*(_0_5+x*x*(_0_0416 + x*x*(_0_001387 + x*x*(_0_0000248 + _0_000000275*x*x))))
addss xmm0, _1 // xmm0++
movss x, xmm0 // das ergebnis in x rüberschieben
}
return x;
}
Auch hier gilt, dass weitere Reihenglieder die Genauigkeit erhöhen. Aber auch schon jetzt setzt der Datentyp float eine Genauigkeitsgrenze fest, da dieser nur 7-8 Nachkommastellen und sogar nur ungenau speichern kann.
Als Abschluss kann man hier die Funktionen als C++ Code mit den dazugehörigen Tests downloaden.
Schnelle trigonometrische Funktionen
Hi
Ich bin hier so reingestolpert. Darf ich dir mal einen Tip geben? Du meinst wahrscheinlich bei einer „Genauigkeit“ von einer Nachkommastelle und 3.600 Werte eine Auflösung der Stützstellen von 0,1°? Wenn das der Fall ist, dann brauchst du keine 3.600 Werte sondern nur 1/4 davon. Denn für ein
0 ≤ x ≤ p/2 ⤇ y = Sin(x)
pi/2 ≤ x ≤ pi ⤇ y = Sin(pi-x)
pi ≤ x ≤ 3pi/2 ⤇ y = -Sin(x-pi)
3pi&2 ≤ x ≤ 2pi ⤇ y = -Sin(2pi/x)
Damit lassen sich alle Werte ermitteln, wenn man die Werte zwischen 0 und pi/2 kennt.
Und das sind nur 900 Punkte, die man abspeichern muss.
Und bei einer Reihenentwicklung brauchst du eigentlich nur um pi/4 zu entwickeln und hier dürfte eine Reihe mit einem Element reichen.
Danke, beides ist mir bewusst =)
Es ist nur so, dass man dann bei Lockup Tabellen noch mehr Fälle zu unterscheiden müsste, was wahrscheinlich länger dauern würde, als wenn man etwas mehr Werte abspeichert.
Bei der Reihenentwicklung ist es ähnlich. Das Umrechnen auf das gültige Intervall dauert schon in dieser Version länger als die Berechnung selbst.
Man könnte aber trotzdem etwas experimentieren und schauen ob es eine leistungsfähige Lösung gibt.
Wenn man diese Regel nutzt, die ich genannt habe, dann genügt es das Näherungspolynom in den ersten beiden Termen zu berechnen – also x – x3/3!. Denn dieses Polynom liefert von -π/2 bis +π/2 eine gute Näherung. Du hast bis x9 ausgewertet, weil du eine Auswertung brauchtest von -π bis π.
Du hast aber die Formel „sin(x + 2zπ)“ genannt und damit erklärt, dass ein beliebiges x abbildest auf das Intervall 0 bis 2π. Mit einer ähnlichen Argumentation kannst du auch auf -π bis π abbilden.
Wenn dann dein x die Bedingung |x| ≤ π/2 erfüllt, dann gilt die deine erste Formel. Wir sparen also mit einem Vergleich drei Terme des Näherungspolynoms.
Das offene Problem sind die Punkte oberhalb von π/2 und unterhalb von -π/2. Der erste Vergleich oben ist gemacht, wir haben einmal gesagt, was zu tun ist, wenn |x| ≤ π/2 gilt. Gilt das nicht, dann sind wir im zeiten Teil, wir sind also immer noch im ersten Vergleich.
Nun brauche ich nur den x-Wert auf x‘ = sign(x) * (|x| – π/2) abbilden. Etwas „beschleunigt“ könnte man sagen:
If x < 0 then
x' = π/2 – x
Else
x' = x – π/2
End If
Also ein zweiter Vergleich und eine weitere Additon.
In der Summe ersetzen zwei Vergleiche und eine Addition die Terme x5/5!, x7/7! und x9/9! sowie drei Additionen. Ok, du hast selbst die Berechnung dieser Terme auch schon optimiert. Also die Einsparung ist nicht sooo groß wie hier dargestellt. Aber trotzdem, dein Polynom wird kürzer.
Noch ein Tip: Du solltest im WordPress <sup> und <sub> freischalten!
Schön, dass einer mitdenkt ^^
Ich werde es mir auf jeden Fall näher anschauen ;)
Habe momentan noch Klausurstress…
Zu WordPress: Ja, muss sowieso die erlaubten Tags anzeigen lassen.
Ich muss noch eine Korrektur anfügen. Ich habe auf die Cosinusfunktion gesehen und behauptet ich hätte vom Sinus geschrieben:
In der Schleife lautet es:
If x < 0 then (meint x = π/2)
x‘ = π – x
End If
Aber kümmere dich erst mal um deine Klausuren und viel Erfolg. …
Sorry das WP hat was geschluckt aber ich denke mal du findest es raus ..
So, ich habe mir noch mal alle Kommentare durchgelesen und meine Antwort, die im zweiten Kommentar steht, hat sich nicht verändert.
Wenn ich noch eine Fallunterscheidung einbaue, dann spare ich mir eine Konstantendeklaration und vier SSE-Anweisungen, dafür kommen aber bereits durch den zweiten Vergleich mindestens vier weitere Anweisungen dazu. Ich habe es auch ausprobiert, es bringt überhaupt keinen Geschwindigkeitsvorteil.
Noch mehr Reihenglieder zu entfernen geht auch nicht, denn es müssen mindestens drei sein, also x – x3/3! + x5/5!. Ansonsten ist die Ungenauigkeit einfach zu groß.
Ich habe die Assembler-Version weiter optimiert und konnte so einige Anweisungen sparen.
Sehr interessanter Artikel und extrem nützlich, da in Spielen ja enorm viele solcher trigonometrischen Funktionen aufgerufen werden.
Könntest du vllt ein paar Stichwörter bzw. Links liefern, wie sich das Ganze mit dem Tangens und dem Sin^-1 und Cos^-1 verhält? Gibt es dort ebenfalls solche Approximationsfunktionen oder andere Optimierungstricks?
Außerdem schreibst du an einer Stelle „Es handelt sich um 1:1 Konvertierung der C++ Version der Funktion.“. Wieso ist dann die Assembler Implementierung trotzdem schneller? Wandelt der Compiler den C++ Code nicht in exakt diesen Code um?
Ich fände es auch mal interessant, ob es Möglichkeiten gibt den Sinus und Kosinus in „einem Atemzug“ zu berechnen, denn für einige Dinge braucht man sowieso beides (z.B. Rotationsmatrizen). Zwar lässt sich aus dem Sinus der cos errechnen: cos(x) = sqrt(1.0f – sin(x) * sin(x)), doch dies ist bedauerlicherweise langsamer als die Berechnung des cos direkt…
Sorry, ich muss meine letzte Aussage zurücknehmen, sqrt(1.0f – sin(x) * sin(x)) ist doch schneller als cos(x) direkt zuberechnen (auch als deine optimierte Methode), ich hatte einen Fehler beim testen gemacht…
Mit Cos^-1 meinst du sicherlich Arcus Cosinus. Arccos(x) ist nicht(!) Cos(x)^-1 = 1/Cos(x), wenn es auch so auf dem Taschenrechner steht. Nur mal so als Hinweis.
Man man kann Arccos und Arcsin genauso leicht durch Reihenentwicklung http://de.wikipedia.org/wiki/Arcussinus#Reihenentwicklungen implementieren.
Ich beziehe mich auf meine C++-Funktion, nicht die von STL. Compiler erzeugen nicht immer den optimalen Code (wobei ich nicht behaupten kann, dass mein Code optimal ist).
Wenn könntest deine Gleichung zuerst analytisch als Reihenentwicklung aufschreiben und dann Zusammenfassen und schauen, ob sich daraus ein Vorteil ziehen lässt.
Ist dein Programm zu langsam?
Dass cos(x)^-1 != 1/cos(x) ist ist mir klar :D
Danke für den Link, das werde ich mir dann mal anschauen und es versuchen zu implementieren. Ich bezog mich schon auf deine Funktion, ich meine wenn der Assembler code eine 1:1 Umsetzung deiner C++ Implementierung ist dann gehe ich davon aus, dass der Compiler den gleichen Code erzeugen würde => keine Performanceunterschiede.
Und nein mein Programm ist nicht zu langsam, ich beschäftige mich nur zur Zeit mit Optimierungsmöglichkeiten in C++ im Rahmen der Spieleprogrammierung, da ich die Thematik interessant finde. ;)
Jetzt bist du aber ins Fettnäpfchen getreten, denn cos(x)^-1 == 1/cos(x) != arccos(x) :P
An welche Gleichung mit trigonometrischen Funktionen im Bezug auf die Rotationsmatrix hast du gedacht?
erwischt -.-
Naja, Rotationsmatrizen im 3D Raum (egal um welche Achsen) beinhalten allgeimen sin(x) und cos(x) innerhalb von einer Matrix. D.h. zu beginn der Erzeugung dieser Matrix ist es nützlich sowohl den Sinus als auch den Cosinus bereits einmal zu berechnen. Und wenn man dann mit der in deinem Artikel vorgestellten Methode den Sinus berechnet, so lässt sich über sqrt(1.0f – sin(x) * sin(x)), indem man eine Implementierung der sqrt in asm verwendet (http://www.indiedev.de/wiki/Fast_Sqrt) der Kosinus berechnen, etwas schnelleres fiele mir jetzt nicht ein… :D
Mit Matrizen hatte ich mich vor einigen Tagen intensiver beschäftigt, und das Resultat war dieser Artikel (http://www.indiedev.de/wiki/Matrix#Rotation) in dem du auch noch einmal den Aufbau der Rotationsmatrizen findest, falls du sie nicht im Kopf hast, dann wirst du sofort sehen, dass in jeder dieser Matrizen sowohl der Sinus als auch der Cosinus vorkommt und deshalb dachte ich gleich an Optimierungsbedarf ;)
Ja, das habe ich mir gedacht, dass du es meinst. Solange Kosinus und Sinus die gleichen Argumente haben, dann kann man es schön zusammenfassen und in einer Reihe ausrechnen. Bei der Rotationsmatrix mit mehreren Winkeln hat man viele Terme und man muss abwiegen ob es sich lohnt es als Reihe zu berechnen. Vor allem bei einer Rotationsmatrix kann es durch Ungenauigkeiten zu unschönen „Zuckungen“ in der Rotation kommen.
So, inzwischen habe ich auch noch etwas in asm rumexperimentiert und bin auf folgendes gekommen:
inline float __declspec(naked) __fastcall Sin(float f) {
_asm {
fld dword ptr[esp+4]
fsin
ret 4
}
}
Diese Methode ist bisher mit Abstand das schnellste und schlägt in der Geschwindigkeit sogar deine Approximation (und ist exakt). Oder habe ich irgendwelche Nachteile übersehen?
Zur schnellen Berechnung von SinCos gibts hier etwas: http://www.gamedev.net/topic/389251-sse2-replacements-for-some-maths-functions/ (letzter post)
Na das ist nichts anderes als std::sin.
Nach meinem Testcode ist die Näherung ist immer noch etwa 2-3 mal schneller (wobei komischerweise die SSE Version langsamer als die C++ Version ist, hat wohl etwas mit dem Prozessor zu tun). Zeige mir deinen Testcode (der Testcode im Beitrag ist nicht ganz korrekt, wie ich sehe).
Was stimmt denn an dem code im beitrag nicht (Ich kenne mich in Assembler noch nicht so gut aus)?
Also bei mir ist diese Methode fast doppelt so schnell wie die des std::sin. Getestet wurde mit 1.000.000 sinus Berechnungen unterschiedlicher Winkel… 35ms:65ms
Also ich verwende VS2010 mit normalen Release-Einstellungen+SSE2 Satz auf nem Prozessor Intel i520 und die Näherung ist um einiges schneller.
Man muss sicherstellen, dass die berechneten Werte nicht wegoptimiert werden. Füge doch oberhalb der Schleife eine float Variable dazu und summiere die berechneten Werte immer auf.
Zeige doch einfach den Code, denn meine Kristallkugel ist kaputt ;)
Ja, das mit dem wegoptimieren habe ich schon beachtet. Das Problem könnte evtl sein, dass ich im Debug Mode getestet habe, nur im Release gibts grad noch ein paar Linker errors, da muss ich erst nochmal näher reinschauen.
Naja, es hatte schon einen Grund wieso ich den Code nicht zeigen wollte…er nutzt einige meiner Klassen und Makros. Außerdem lassen sich Kommentare hier nicht sehr schön formatieren :D
inline float __declspec(naked) __fastcall Sin(float angle) {
_asm {
fld dword ptr[esp+4]
fmul dword ptr[ToRadians]
fsin
ret 4
}
}
//…
float f = 0.0f;
float angle = 0.0f;
BEGIN_PROFILE(„StdSin“);
for(int i = 0; i < tests; i++) {
f *= std::sin(angle * ToRadians);
angle ++;
}
END_PROFILE("StdSin");
angle = 0.0f;
BEGIN_PROFILE("MySin");
for(int j = 0; j FindProfilerInfo(„StdSin“).AvgTime; //65
double t2 = Profiling->FindProfilerInfo(„MySin“).AvgTime; //35
Winkel sind in Grad, statt Bogenmaß. Ich werde jetzt aber erstmal schauen ob ichs im Release Build zum laufen bringe (da das Projekt noch einige andere Libs (DirectX, meine Engine, etc.) einbindet)…
Seltsam…er hat
j < tests; j++) {
f *= Sin(angle);
angle ++;
}
END_PROFILE("MySin");
angle = 0.0f;
double t1 =
weg geschnitten…
Hier auf Pastebin: http://pastebin.com/tDPJpSFQ … hoffentlich lesbarer -.-
@Kevin:
Und hast du es noch mal getestet?
Anmerkung zu deinen Artikeln auf http://www.indiedev.de: Die Quellenangabe (Referenz) muss direkt an der Stelle sein (spätestens am Satzende), an der du dich auf sie beziehst. Ein Link am Ende des Artikels reicht nicht aus (siehe auch z.B. die Vorgehensweise bei Wikipedia). Beziehst du dich in mehreren Abschnitten (die nicht direkt nebeneinander liegen bzw. anderes Thema oder anderen Gedanken betreffen) auf die gleiche Quelle, so muss die Referenz wiederholt gesetzt werden.
Wenn du mal später eine wissenschaftliche Arbeit schreibst, wird darauf sehr viel Wert gelegt ;)
Seit 1-2 Tagen habe ich das Internet durchforstet nach einer Sinus-Implementation die schnell ist, nicht unbedingt präzise sein muss und deren Theorie ich verstehe. Ich bin recht schnell auf den Ansatz der Taylorreihe gestoßen, doch war der immer recht kompliziert erklärt worden. Selbst wenn du nicht den Code hier gezeigt hättest könnte ich nun die Sinus-Funktion selber implementieren (Auch wenn man mit 1/4 der Werte auskommt, wie in den Kommentaren besprochen).
Vielen Dank für diesen tollen Artikel!
Kann mir mal bitte jemand erklären, wieso für den Wert -2.968630 (knapp vor -170 Grad) zwar der Cosinus-Ersatz perfekt dasteht, der für den Sinus aber nicht? Statt korrekter -0.172 liefert er -0.077. Das ist weit weg von 2-6 Stellen hinter dem Komma…
Bei mir funktioniert das obige Programm recht genau (wie der Test beweist)
function fsin(s)
local x=s
if x>180.0 then
x=(x+180)%360-180
elseif x90.0 then
x=180.0-x
elseif x<-90.0 then
x=-180-x
end
-- print (x)
x=math.pi*x/180.0
local p=x
local f=1
local r=p
-- print (x)
-- 6 Stellen hinter dem Komma fuer 1/4-Kreis
for j=2,10,2 do
p=-p * x * x
f=f*j*(j+1)
r=r+p/f
end
return(r)
end
function fcos(c)
return fsin(90-c)
end
-- TEST:
for w=-360,360,1 do
print (1-fsin(w)*fsin(w)-fcos(w)*fcos(w))
end