[C++] Schnelleres sin/cos/tan

[Krishty]

Pah, ihr seid doch alles Kunstbanausen.

Die richtige Lösung™ wäre ein Fixkomma-Datentyp der automatisch von 360° auf 0° umwrappt, zusammen mit einem Satz optimierter trigonometrischer AVX-Vektorfunktionen (sin(·) bzw. cos(·) usw. geben dann für · als dieser Datentyp einen half/float/double zurück). Oder man macht stattdessen was sinnvolles mit seiner Zeit.

Eben angepackt (allerdings mit unsigned (normalized) int). Wie gut sich das anfühlt … den Quadranten eines Winkels für die eigenen trigonometrischen Funktionen bestimmen? Bloß die oberen Bits abfragen!

Doof nur, dass Quaternionen in einen Winkel bis zu 720° umgewandelt werden können.

Wo genau, außer bei Quaternion-Interpolation, ist das von Nutzen?

Bogenmaß ist dermaßen ätzend zu lesen und zu bearbeiten, dass ich darüber nachdenke, auf Grad oder gleich Gon (360° == 400 gon) umzustellen.

Das hier zu Common7\Packages\Debugger\autoexp.dat unter [Visualizers] hinzugefügt:

    Math::Angle{
        preview (
            #(
                $e.unorm * 1.4629180792671596810513378043098e-9, " (", $e.unorm * 0.00000008381903171539306640625, " rad)"
            )
        )
    }

ang.png (5.88 KiB) 16354 mal betrachtet

[Jörg]

Da ist wohl ein G verschwunden ...

[Krishty]

Ohja. Das war mal andersrum … toll, wie blind man für eigene Fehler ist.

[eXile]

Eben angepackt (allerdings mit unsigned (normalized) int). Wie gut sich das anfühlt … den Quadranten eines Winkels für die eigenen trigonometrischen Funktionen bestimmen? Bloß die oberen Bits abfragen!

Das ist schön zu hören. :) Hast du dir neue Minimax-Polynome gebastelt, um die höhere Genauigkeit (in Vergleich zu float) zu nutzen?

Wo genau, außer bei Quaternion-Interpolation, ist das von Nutzen?

Naja, genau da. Und zwar sowohl beim SLERP mit Quaternionen wie auch der Konvexkombination von Quaternionen. Ich wollte eher davor warnen, im Bereich der Quaternionen auch nur in irgendeiner Weise mit Winkeln zu rechnen, genauso wenig wie man im Bereich der Rotationsmatrizen niemals plötzlich anfangen sollte, mit Winkeln zu rechnen.

AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAH. :twisted:

Du hast noch nie Code von Mathematikern gesehen. Ich sage nur: Jeder Variablenname ist maximal zwei Buchstaben lang. AAAAAAAAA!

[Krishty]

Eben angepackt (allerdings mit unsigned (normalized) int). Wie gut sich das anfühlt … den Quadranten eines Winkels für die eigenen trigonometrischen Funktionen bestimmen? Bloß die oberen Bits abfragen!

Das ist schön zu hören. :) Hast du dir neue Minimax-Polynome gebastelt, um die höhere Genauigkeit (in Vergleich zu float) zu nutzen?

Nein. Zum einen sind meine Mathematikfähigkeiten bestenfalls mit legasthenisch zusammenzufassen, so dass ich nach fünf Papern dazu immernoch keinen Plan habe, wie man auf die Dinger kommt außer durch Brute Force; zum anderen weiß ich nicht, warum Integer-Winkel andere Polynome gebrauchen können sollten als Gleitkommawinkel, und zu guter Letzt fällt es mir sehr schwer, eine anständige Referenz für die tatsächlichen Werte zu finden.

Einsatztauglich ist bisher nur die float tan()-Funktion, und auch die wartet noch ihre finale Prüfung ab. Obwohl die Werte absolut konsistent zu einem von Greens Papern sind, habe ich dauernd 1 ULPs Abweichung gegenüber dem double tan() der VC-CRT gemessen, und ich weiß nicht, wer nun recht hat. Weil ich keine Arbitrary-Precision-Number-Bibliothek gefunden habe, die sich via Download & Include direkt unter Visual Studio kompilieren ließ, hätte ich auch nicht die geringste Referenz zum Testen etwaiger double-Implementierungen (auf die ich im Augenblick sowieso verzichten will, weil ich sonst nie mehr vorankomme).

Für heute abend habe ich mir http://www.ganssle.com/approx.htm vorgenommen. Sieht so aus, als ob man dort gepflegt abkupfern könnte statt sich Minimax selber zu bestimmen oder das bereits außer Druck befindliche Standardwerk von Hart zum Nachschlagen zu bemühen, das jemand meiner mathematischen Bildung eh nicht bedienen könnte.

Für die Optimierung auf SSE / AVX habe ich eh keine Zeit, zumal ich den Wechsel zu Clang im Auge habe und mich Intrinsics dabei absolut abschrecken. Von nicht vektorisiertem Text, der Sinus und Cosinus desselben Winkels simultan berechnet und sich durch die Festkommadarstellung der Winkel das leidige Modulo spart, verspreche ich mir bereits mehr Leistung als von normalem seriellen float sin(), cos().

Letztendlich geht es auch weniger um Geschwindigkeit als viel mehr um Abnabelung von der C++-Laufzeitbibliothek und darum, es richtig zu machen statt so, wie es alle tun.

[eXile]

Mmmh ja, ich stelle auch gerade fest, dass meine intuitive Idee ja schwere Schwierigkeiten mitbringen würde. Meine Idee war, man nehme sein Minimax-Polynom (das man übrigens mit dem Remez-Algorithmus unter bestimmten Voraussetzungen wie Stetigkeit usw. bestimmen kann) mal in Standardform (welche ja in der Regel nicht die effizienteste Form ist):
\($$P(x) = a_0 x + a_1 x^3 + a_2 x^5 + a_3 x^7 + a_4 x^9 + a_5 x^{11} + a_6 x^{13} + a_7 x^{15}$$\)Ich dachte jetzt, oh ja, man kann ja den Konvertierungsfaktor von Math::Angle einfach in die \($a_i$\) reinbaken. Tja, falsch gedacht. Wenn man die Konvertierung von unserem Math::Angle \($\alpha \in [0, \mathrm{0xffffffff}]$\) in einen float \($x \in [0, 2\pi[$\) mittels
\($$x = \alpha \cdot c \quad \text{mit } c \approx 1.4629180792671596810513378043098 \cdot 10^{-9}$$\) durchführt, kriegt man
\($$P(x) = (a_0c)\alpha + (a_1c^3) \alpha^3 + (a_2c^5) \alpha^5 + (a_3c^7) \alpha^7 + (a_4c^9) \alpha^9 + (a_5c^{11}) \alpha^{11} + (a_6c^{13}) \alpha^{13} + (a_7c^{15}) \alpha^{15}$$\)Spätestens der gebakte Koeffizient \($a_7c^{15}$\) ist jedoch nicht mehr sinnvoll speicherbar (liegt bei ca. \($1.46 \cdot 10^{-135}$\)).

DIe intuitive Lösung, d.h. einfach \($$P(x) = a_0 (\alpha c) + a_1 {(\alpha c)}^3 + a_2 {(\alpha c)}^5 + a_3 {(\alpha c)}^7 + a_4 {(\alpha c)}^9 + a_5 {(\alpha c)}^{11} + a_6 {(\alpha c)}^{13} + a_7 {(\alpha c)}^{15}$$\)zu berechnen, ist somit die richtige Lösung. Keine Probleme mit irrwitzigen Zahlen. Hier ist Reinbaken in die Koeffizienten eine Scheißidee. (Und ja, noch einmal: Zum Berechnen nicht die Standardform nehmen. Aber das Problem ist hier sehr gut darzustellen und tritt auch bei anderen Polynomformen auf. Kein Reinbaken!)

Super, nun habe ich so einen langen Beitrag geschrieben, nur um dazulegen, warum ich vor ein paar Stunden falsch lag. ;)

[Krishty]

Hmm. Inwiefern ist e-135 nicht sinnvoll speicherbar? Für mehr als fünf Stellen Genauigkeit muss man eh mit double-Konstanten rechnen; und dabei reicht der Exponent bis e-308.

[eXile]

Hmm. Inwiefern ist e-135 nicht sinnvoll speicherbar? Für mehr als fünf Stellen Genauigkeit muss man eh mit double-Konstanten rechnen; und dabei reicht der Exponent bis e-308.

Meines Wissens nach hat double ein Maschinenepsilon von \($2^{-53} \approx 1.11 \cdot 10^{-16}$\) (daher reichen 17 Nachkommastellen in Dezimalschreibweise). Die größe bzw. kleinste mittels double darstellbare Zahl beträgt ca. \($\pm 1.79\cdot10^{308}$\). Für ein Maschinenepsilon in Größenordnungen um die \($10^{-135}$\) braucht man nach meinen Berechnungen ca. 449 Mantisse-Bits, d.h. ca. 512-bittige Gleitkommazahlen.

[B.G.Michi]

das "Maschinenepsilon" das du meinst ist der Abstand von 1.0 * 10^0 zum nächst größeren, darstellbaren Wert (1.0 + 2^-53) = 1.000..011 * 10^0. Also immer mit Exponent 0. Desswegen ist 1.46 * 10^-135 trotzdem speicherbar.

[eXile]

Tatsache. Mein Wissen über Gleitkommazahlen ist schon wieder zu sehr eingerostet. Richtig, die 17 Dezimalstellen beziehen sich auf signifikante Dezimalstellen, d.h. keine vorläufigen Nullen.

Auf der Sonnenseite: Meine Optimierung funktioniert damit wohl. Fröhliches baketieren.

[Matthias Gubisch]

Ich find die Diskussin hier echt interessant, aber könnte man die nicht in einen anderer Thread auslagern?
Hier geht mir das zu schnell verloren

[Krishty]

Hier ist mein 32-Bit-tan(). Ich bin noch nicht dazu gekommen, alle 2³² Werte zu testen oder mir den Maschinentext anzusehen, weil ich im Augenblick vor einem x86-only-Eee PC sitze. 0° ergibt 0; 90° ergibt -INF; dazwischen scheint immer mal wieder ein ULPs Abweichung im Vergleich zur Visual C++ CRT-Routine zu sein. Ich weiß nicht, ob das an mir, an Herr Green, oder an der CRT liegt.

Code: Alles auswählen

struct Angle { unsigned int unorm; };
typedef float Float4B;
typedef double Float8B;

Float4B tangentOf(
	Angle angle
) {
	// Tangent can be divided into four sectors:
	//		• 0: [  0°,  45°[
	//		• 1: [ 45°,  90°[
	//		• 2: [ 90°, 135°[ (inversion of sector 1)
	//		• 3: [135°, 180°[ (inversion of sector 0)
	//	At 180°, the function repeats.
	// Since the angle is given as an unsigned normalized integer, and each sector is exactly the 8th of the overall range, the
	//	sector can be determined from bits 30 and 29 (mask 0x60000000):
	//		• 0x00000000: 0
	//		• 0x20000000: 1
	//		• 0x40000000: 2
	//		• 0x80000000: 3
	auto const sectorID = angle.unorm & 0x60000000;

	// Sector 0 is computed via Green, Robin — "Faster Math Functions", pg. 15 (http://www.research.scea.com/research/pdfs/RGREENfastermath_GDC02.pdf):
	//
	//					0.999999986 x - 0.0958010197 x³
	//		tan(x) = ————————————————————————————————————
	//				1 - 0.429135022 x² + 0.00971659383 x²²
	//
	// When converting the angle to a floating-point number, one must multiply by (2 pi) / 360. This multiplication can be saved
	//	by applying it to the constant factors:
	//
	//					1.4629180587863065713111022695911e-9 x - 2.9993707578799427148695747771533e-28 x³
	//		tan(x) = ——————————————————————————————————————————————————————————————————————————————————————
	//				1 - 9.1840443709068308636681822066321e-19 x² + 4.4503490744640484968825016740227e-38 x²²
	//
	// Sector 1 is computed by a similar term, but with numerator and denominator flipped after assigning
	//
	//		x = 0x40000000 - x
	//
	//	where 0x40000000 is pi/4:
	//
	//				1 - 9.1840443709068308636681822066321e-19 x² + 4.4503490744640484968825016740227e-38 x²²
	//		tan(x) = ——————————————————————————————————————————————————————————————————————————————————————
	//					1.4629180587863065713111022695911e-9 x - 2.9993707578799427148695747771533e-28 x³
	//
	// Gather all constants in the same memory block, which is aligned on a cache line boundary and does not exceed cache line
	//	size. Since constants will always be read from read-only memory (even trivial ones like 1.0), this will make them
	//	consume just a single cache line. If they were scattered all over the read-only section, the CPU would in worst case be
	//	touching one cache line per constant.
	__declspec(align(64)) static struct {
		Float8B one; //  1.0
		Float8B a;	 //  0.999999986   * 1.4629180792671596810513378043098e-9,
		Float8B b;	 // -0.0958010197  * 3.1308338546629715204060346522108e-27,
		Float8B c;	 // -0.429135022   * 2.1401293066467156958511258973014e-18,
		Float8B d;	 //  0.00971659383 * 4.5801534491681520631005954830806e-36
	} const constants = {
		 1.0,
		 1.4629180587863065713111022695911e-9,
		-2.9993707578799427148695747771533e-28,
		-9.1840443709068308636681822066321e-19,
		 4.4503490744640484968825016740227e-38
	};

	// Modulate the range to 180° only AFTER the sector is determined. This will hopefully break any dependency of branching on
	//	modified angle, thus suppress mispredictions.
	angle.unorm = angle.unorm & 0x7FFFFFFF;

	// Adjust x for the computation according to the angle's sector.
	Float8B x;
	if(0x00000000u == sectorID) {
		x = Float8B(angle.unorm);
	} else if(0x20000000u == sectorID) {
		x = Float8B(0x40000000u - angle.unorm);
	} else if(0x40000000u == sectorID) {
		x = Float8B(angle.unorm - 0x40000000u);
	} else { assert("binary logic broke", 0x60000000u == sectorID);
		x = Float8B(0x80000000 - angle.unorm);
	}

	// These are always required:
	auto const x² = x * x;
	auto const x³ = x² * x;
	auto const x²² = x² * x²;
	auto const termA = constants.a * x + constants.b * x³;
	auto const termB = constants.one + constants.c * x² + constants.d * x²²;

	// Perform final division and adjust sign according to the sector:
	// TODO: bake negation into constants; operator - will touch global memory
	if(0x00000000u == sectorID) {
		return Float4B(termA / termB);
	} else if(0x20000000u == sectorID) {
		return Float4B(termB / termA);
	} else if(0x40000000u == sectorID) {
		return -Float4B(termB / termA);
	} else { assert("binary logic broke", 0x60000000u == sectorID);
		return -Float4B(termA / termB);
	}
}

[Krishty]

Soo, die x86-Version ist getestet. Sie weicht überall 0 oder 1 ULP von der VC-CRT-Version ab, außer bei den Winkeln 0x3FFFFFFF und 0x40000000 (89,99999991618096828460693359375 und 90 °), wo der Unterschied weitaus höher ausfällt. Die Funktion ist auf einem Intel Atom ungefähr 3,5× so schnell wie tan(double) der VC-CRT, wobei das aber nichts über die Leistung auf PCs aussagt.

[dot]

Was ich noch interessant finden würde: Ist sie auch schneller als der Intrinsic?

[Krishty]

Als wer?

[dot]

Naja, eben als die fptan Instruction der FPU...

EDIT: Ok, wie es aussieht lässt sich MSVC nicht dazu bewegen, fptan zu generieren...

[Krishty]

Keine Ahnung. Interessiert mich auch nicht, denn 32-Bit-x86 ist tot. (Ja – bigott, das von einem 32-Bit-Atom aus zu schreiben.)

Ich meine aber mal gelesen zu haben, dass die Leistung von fsin / fcos / ftan im Vergleich zu den SSE-optimierten VC-CRT-Funktionen nicht sehr berauschend sei.

[dot]

Ah ok, macht Sinn, ich dachte, dass die entsprechenden FPU Instructions evtl. relativ schnell wären, aber fptan scheint je nach CPU eine Latency in der Gegend von 100-200 Cycles zu haben, das ist natürlich elendslahm...

[Krishty]

Ja; wobei es aber perfekt lokal und kompakt ist, während die Funktion erstmal einen Satz Text, Daten und Register verschlingt. In Ausnahmefällen kann es vielleicht doch schneller sein, aber im Normalfall ist es tatsächlich elendslahm.

[Krishty]

Auf einem Core i7 x86 ist die neue Funktion 3× so schnell wie tan(double); als x64 zwar nochmal 30 % schneller als ihr x86-Gegenstück, aber nur noch doppelt so schnell wie das VC-CRT-tan(double).

[BeRsErKeR]

Ich find die Diskussin hier echt interessant, aber könnte man die nicht in einen anderer Thread auslagern?
Hier geht mir das zu schnell verloren

Das Abspalten ist ja schön und gut, aber irgendwie weiß man jetzt nicht mehr so ganz um was es hier geht. An was bastelt Krishty denn? Mir fehlt ein wenig der Zusammenhang und ich bin gerade zu faul den Jammer-Thread zu durchforsten.

[Artificial Mind]

Wenn ich das richtig verstande habe, dann geht es um Winkelberechnungen (und Konsorten) auf einem Datentyp, der von selbst nach 360° "wrappt". Und dafür ist ein unsigned Integer geeignet (0 -> 0°, 0xfff... -> 360° (-eps?)).

[BeRsErKeR]

Naja aber das klang irgendwie so, als ob Krishty irgendeine Lib baut oder ähnliches.

[Matthias Gubisch]

Krishty versucht eben genau diesen Datentyp zu bauen.

Soweit ich das verstanden habe weil das 1. schneller ist als die von der Runtime zur verfügung gestellten Funktionen und er 2. die Runtime loswerden will ;)

[dot]

Nun, es ist vor allem auch recht elegant. Ein potentielles Problem, das ich dabei sehe, ist allerdings, dass man damit z.B. keine Winkelgeschwindigkeiten größer als 360°/s darstellen kann...

[Krishty]

Der Kram fließt einfach so in mein Framework ein. Es geht mir in erster Linie darum,

es richtig zu machen statt so, wie es alle tun

, dann um

Abnabelung von der C++-Laufzeitbibliothek

und schließlich um

Geschwindigkeit

.

Das Argument mit den Winkelgeschwindigkeiten ist schwach. Zum einen sind Geschwindigkeiten, Kräfte, usw. nicht wirklich für Festkommadarstellung geeignet (lies: das würde ich nicht machen) und zum anderen weiß ich nicht, warum man den Sinus/Kosinus/Tangens einer Winkelgeschwindigkeit ermitteln sollte.

[dot]

Darum ja auch nur potentielle Probleme ;)
Klar macht der Cosinus einer Winkelgeschwindigkeit an sich jetzt nicht unbedingt Sinn. Die von mir implizierte Frage ausgesprochen: Winkel existieren ja normalerweise nicht in einem Vakuum, bevor ich jetzt den Cosinus eines Winkels ausrechne, will in in der Regel z.B. erstmal den Winkel an sich ausrechnen, wie effizient lassen sich denn also solche Fixkommawinkel mit normalen float Berechnungen mischen?

[Krishty]

Guter Punkt: Wo kommen denn deine Winkel her?

Rotationen, Winkelgeschwindigkeiten und Drehmomente sind ja vektorielle Größen. Die stellt man als Quaternions o.ä. dar, aber (hoffentlich!) nicht mit Euler-Winkeln. Hier also erstmal nicht viele Spuren von sin(), cos() und tan(). Wo kommen die Quaternions her, die dort aufaddiert werden? Bei mir gibt es nur zwei Stellen, wo ich tatsächlich mit Winkeln rechne:

Zum einen die Benutzeroberfläche („wird die Maus nach rechts bewegt, dreht sich der Betrachter um x Grad um die Y-Achse“ oder „Dieses HUD-Element zeigt die Kompassrichtung an“) – Benutzereingaben kommen aber eh meist als Festkommazahlen an (WinAPI), Kompassrichtungen sind selten.

Zum anderen das Laden fremder Dateiformate. Die geben Rotationen meist in float an; ich habe auch eins, das nativ mit 16-Bit-Festkommazahlen arbeitet.

Die Konvertierung float->int geht dank Hardware-Unterstützung (seit SSE2) recht zügig. Detailliert schreibe ich darüber aber im Status-Update, das ich gerade vorbereite (da habe ich eine böse Überraschung erlebt).

[Krishty]

Aaaalso:

Ich habe mit Entzücken festgestellt, dass Herr Greens Näherungsfunktion nicht für [0, pi/4] definiert ist, sondern für [-pi/4, pi/4]. Das bedeutet: Die Negierung am Ende der Funktion entfällt, wenn man eine negative Gleitkommazahl als Eingabe herstellen kann. Wie schafft man das?

Hier kommt eine weitere glückliche Entdeckung ins Spiel: unsigned int zu double ist deutlich langsamer als int zu double. Der Grund ist, dass SSE2 mit CVTSI2SD eine Funktion anbietet, die nur für int funktioniert. Eine unsigned int-Konvertierung behandelt der Compiler, indem er erst als int konvertiert; dann testet, ob das höchstwertige Bit gesetzt ist; und dementsprechend die double-Konstante 4294967296.0 aufaddiert oder es sein lässt.

Wird der Eingabewinkel also zu int konvertiert, und die Konstanten für Sektor 3 und 4 derart verändert, dass eine negative Zahl herauskommt, beschleunigt das einerseits die Konvertierung zu double und schmeißt andererseits die Negierung am Schluss raus.

Interessanterweise wird es dadurch aber keinen Deut schneller, obwohl der Ausführungspfad kürzer geworden ist und weniger Sprünge beinhaltet. Ich forsche weiter.

[Krishty]

Kurze Frage: VS-CRT-x64-tan() gibt mir für 0x00000005 * 1.4629180792671596810513378043098e-9 einen Wert von

     7.3145903963357981e-009

aus. Für den gleichen Input negiert allerdings

    -7.3145909558344820e-009

(die roten Stellen weichen vom Windows-Taschenrechner ab). Bei dem Winkel handelt es sich um (-)0,00000041909515857696533203125 °, oder 0x00000005 bzw. 0xFFFFFFFB in meinen Winkeln. Das Ergebnis des positiven Winkels ist richtig, das des negativen Winkels jedoch ab der achten Stelle Unsinn. Fehler in der VS-CRT?

Komischerweise ist mir der Fehler nicht aufgefallen, als ich noch mit ausschließlich positiven unsigned ints (also die Version, die auf Seite 1 gepostet ist) gearbeitet hatte …

Ich habe insgesamt vier bis sechs Stellen über die 2³² Winkel, an denen mir Tests Abweichungen gegenüber der VS-CRT von mehr als einem ULP melden. Das ist, wohlgemerkt, der Vergleich meiner 32-Bit-genauen Tangensfunktion mit der 64-Bit-genauen VC-CRT-Implementierung. Da sind die nicht definierten Winkel wie 90 ° bei; trotzdem beginne ich an der Qualität der VS-Mathefunktionen zu zweifeln.