Schnelle trigonometrische Funktionen

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

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.

$sin(x) = \sum_{n=1}^{\infty} (-1)^n \frac{x^{2n+1}}{(2n+1)!}$

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.

sinus modellierung plot 1

Legende:
Orange = sin(x)
Schwarz = $x - \frac{x^3}{3!} + \frac{x^5}{5!}$

Wie man sieht, wird der Sinus im Koordinatenursprung gut beschrieben, nicht aber an den PI-Grenzen. Wir addieren einen weiteren Term dazu.

sinus modellierung plot 2

Legende:
Orange = sin(x)
Schwarz = $x - \frac{x^3}{3!} + \frac{x^5}{5!} - \frac{x^7}{7!}$
Schon besser, aber die Ungenauigkeit doch noch zu groß. Wir erweitern weiter.

sinus modellierung plot 3

Legende:
Orange = sin(x)
Schwarz = $x - \frac{x^3}{3!} + \frac{x^5}{5!} - \frac{x^7}{7!} + \frac{x^9}{9!}$
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.

sinus modellierung plot 4

Wie gut diese Modellierung wirklich ist, schauen wir uns später an konkreten Zahlen an.

Demnach lautet unsere Funktion $sin(x) = x - \frac{x^3}{3!} + \frac{x^5}{5!} - \frac{x^7}{7!} + \frac{x^9}{9!}$ . 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 $2\pi$ periodisch. Das heißt, wenn man zum Funktionsargument z-mal $2\pi$ dazu addiert oder von dem Funktionsargument z-mal $2\pi$ abzieht, wird sich der Wert der Funktion sin(x) nicht verändert (z ist eine Ganzzahl).

Man könnte auch schreiben:

$sin(x) = sin(x + z * 2\pi)$ oder $sin(x) = sin(x -z * 2\pi)$

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 $x = x-z*2\pi$ 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 $x = x-(-z)*2\pi$ zu $x = x+z*2\pi$ .

Nachdem das Argument im richtigen Intervall liegt, können wir unsere Taylor-Reihe benutzen. Doch zuerst schreiben wir die Taylor-Reihe $x - \frac{x^3}{3!} + \frac{x^5}{5!} - \frac{x^7}{7!} + \frac{x^9}{9!}$ so um, dass die Fakultät- und Division-Operationen verschwinden, da diese sehr langsamer sind.

Wir bekommen
$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$

und durch Umformung

$sin(x) = x*(1 - x*x*(0.1666667 - x*x*(0.0083333 - x*x*(0.0001984 - 0.0000027*x*x))))$

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 $\pi$ -Grenzen am höchsten.

fastSin Genauigkeit
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: $x - \frac{x^3}{3!} + \frac{x^5}{5!} - \frac{x^7}{7!} + \frac{x^9}{9!} - \frac{x^11}{11!}$ .

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

kosinus

Das Schreiben der Kosinus-Funktion funktioniert analog zur Sinus-Funktion. Man kann aber auch die Gesetzmäßigkeit $cos(x)=sin( \frac{\pi}{2} - x)$ ausnutzen. Man schreibt einfach eine neue Funktion, die unsere Sinus-Funktion mit dem Argument $\frac{\pi}{2} - x$ 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:

$sin(x) = \sum_{n=1}^{\infty} (-1)^n \frac{x^{2n}}{(2n)!}$

Wir nehmen die ersten sechs Glieder.
kosinus modellierung

Legende:
Orange = cos(x)
Schwarz = $1 - \frac{x^2}{2!} + \frac{x^4}{4!} - \frac{x^6}{6!}+ \frac{x^8}{8!}+ \frac{x^10}{10!}$
Schon besser, aber die Ungenauigkeit doch noch zu groß. Wir erweitern weiter.

Wir schreiben das Polynom in leistungsstärkere Variante um:
$cos(x) = 1-x*x*(0.5-x*x*(0.04166667-x*x*(0.00138889-x*x*(0.00002480-x*x*0.000000275))))$

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

fastCos Genauigkeit
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 2.06 kB -

Besucher, die diesen Beitrag fanden, suchten nach:

negativer sinus negativer cosinus
◊
sinuskurven modellieren
◊
- sinus
◊
sinuskurve intervall 20
◊
c++ sse sinus
◊
sse trigonometrische funktionen
◊
modellieren mit dem sinus formeln
◊
funktionen
◊
sinus wert bei pi/6
◊
sinusfunktion intervall (-2,2) darstellung

Tags: assembler, C++, kosinus, optimierung, sinus, SSE, taylor
Kategorie: C++
Letzte Aktualisierung am: 18. August 2010
8 Kommentare

8 Kommentare zu “Schnelle trigonometrische Funktionen”

# Wolfgangam 19. August 2010 um 22:45 Uhr

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.
# Maximam 19. August 2010 um 23:04 Uhr

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.
# Wolfgangam 20. August 2010 um 10:55 Uhr

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.
# Wolfgangam 20. August 2010 um 10:56 Uhr

Noch ein Tip: Du solltest im WordPress <sup> und <sub> freischalten!
# Maximam 20. August 2010 um 15:01 Uhr

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.
# Wolfgangam 21. August 2010 um 08:49 Uhr

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. …
# Wolfgangam 21. August 2010 um 08:51 Uhr

Sorry das WP hat was geschluckt aber ich denke mal du findest es raus ..
# Maximam 26. August 2010 um 16:58 Uhr

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 – x³/3! + x⁵/5!. Ansonsten ist die Ungenauigkeit einfach zu groß.

Ich habe die Assembler-Version weiter optimiert und konnte so einige Anweisungen sparen.

Trackback URI | Kommentare als RSS

virtual | Maxim

Schnelle trigonometrische Funktionen

8 Kommentare zu “Schnelle trigonometrische Funktionen”

Einen Kommentar schreiben

Ein Thema vorschlagen

Letzte Kommentare

Ähnliche Beiträge:

Kategorien

Zufalls Blogroll

Tag Cloud