[C++] Schnelleres sin/cos/tan

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

[Krishty]

Ich würde mich ja kaputtlachen, wenn es nicht so traurig wäre:

179,99999991618096828460693359375° (0x7FFFFFFF)
calc:   -1,4629180792671596820949490891972e-9
VC-CRT: -1.4629183534635463e-009
ich:    -1.4629180e-009

180° (0x80000000)
calc:   0
VC-CRT: -1.2246467991473532e-016
ich:    0.00000000

Es ist immer das gleiche Muster. Der Fehler tritt ausschließlich bei Winkeln >179 ° auf. Das bedeutet: Die VC-CRT benutzt für [0, 180[ ° eine ähnliche Annäherungsfunktion wie ich; darüber allerdings ausschließlich mit positiven Winkeln nach Modulo – und vertut sich dadurch (wie es wohl auch meine Implementierung bis heute abend getan hatte).

Sie ist an diesen Stellen nur auf 3 ULPs genau; und die liegen auch noch so ungünstig, dass selbst nach dem Runden zu float noch eine ULP Abweichung bestehen bleibt.

Hier ist mein Quelltext; die Funktion ist in dieser Form unter x64 doppelt so schnell wie das VC-CRT-tan():

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
	//		• 0x60000000: 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²²
	//
	// Sector 1 is computed by a similar term, but with numerator and denominator flipped after assigning
	//
	//		x = pi/2 - x
	//
	//				1 - 0.429135022 x² + 0.00971659383 x²²
	//		tan(x) = ————————————————————————————————————
	//					0.999999986 x - 0.0958010197 x³
	//
	// When converting the angle to a floating-point number, one must multiply by (2 pi) / 2³². 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²²
	//
	// 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
		Float8B padding[3];
	} 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(int(angle.unorm));
	} else if(0x20000000u == sectorID || 0x40000000u == sectorID) {
		x = Float8B(int(0x40000000u - angle.unorm));
	} else { assert("binary logic broke", 0x60000000u == sectorID);
		x = Float8B(int(angle.unorm - 0x80000000));
	}

	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:
	if(0x00000000u == sectorID || 0x60000000u == sectorID) {
		return Float4B(termA / termB);
	} else { assert("binary logic broke", 0x20000000u == sectorID || 0x40000000u == sectorID);
		return Float4B(termB / termA);
	}
}

Kann die Abweichungen jemand reproduzieren? Nicht, dass ich mich hier die ganze Zeit lächerlich mache … das wäre nämlich schon ein Zufall, wenn die VC-CRT die gleiche Annäherung benützte wie ich.

[Krishty]

HAH, geht doch!
Das Problem ist die Multiplikation 0x7FFFFFFF * 1.4629180792671596810513378043098e-9. Das Ergebnis davon liegt 1 ULP neben dem mathematisch richtigen Ergebnis. Gibt man den Winkel direkt und ohne Multiplikation in die VC-CRT-tan() ein, ist das Ergebnis auch richtig.

Jetzt muss ich nur noch herausfinden, warum die Multiplikation bei 3/2•pi danebenliegt, bei 1/2•pi aber nicht. Wenn die Genauigkeit nicht ausreichen würde, müsste ja jeder Winkel >= 180 ° falsch sein – sind sie aber nicht, sondern nur 7 Sonderfälle.

Nachtrag: 0xBFFFFFFE als Multiplikand ist daneben, 0xBFFFFFFF aber nicht. Wenn 0xBFFFFFFE falsch ist, müsste auch 0x3FFFFFFE danebenliegen – tut es aber nicht.

[BeRsErKeR]

Jetzt muss ich nur noch herausfinden, warum die Multiplikation bei 3/2•pi danebenliegt, bei 1/2•pi aber nicht. Wenn die Genauigkeit nicht ausreichen würde, müsste ja jeder Winkel >= 180 ° falsch sein – sind sie aber nicht, sondern nur 7 Sonderfälle.

Wieso? Die Größe der Zahl hat doch nichts mit der Genauigkeit zu tun. 1/3 ist ja auch nicht genau darstellbar, 1/2 und 1/10 hingegen schon.

[Krishty]

Jetzt muss ich nur noch herausfinden, warum die Multiplikation bei 3/2•pi danebenliegt, bei 1/2•pi aber nicht. Wenn die Genauigkeit nicht ausreichen würde, müsste ja jeder Winkel >= 180 ° falsch sein – sind sie aber nicht, sondern nur 7 Sonderfälle.

Wieso? Die Größe der Zahl hat doch nichts mit der Genauigkeit zu tun. 1/3 ist ja auch nicht genau darstellbar, 1/2 und 1/10 hingegen schon.

Zuerst: 1/3 und 1/10 sind nicht genau darstellbar (weil sie binär keine endlichen Brüche sind), 1/2 schon.

Dann: Doch! doubles haben eine 53-Bit-Mantisse, d.h., sie können jede Ganzzahl bis 2^53 Bits verlustfrei darstellen. Da die Darstellung verlustfrei ist, erwarte ich auch von der Multiplikation nicht mehr als ein halbes ULP Ungenauigkeit – ich kriege aber mehr.


Außerdem habe ich die Berechnung nun ins Horner-Format umgewandelt. Das ist erstmal sinnlos, weil die Hardware kein Multiply-Add anbietet – aber der geneigte Optimizer (nicht VC++) kann nun erkennen, wie das Ganze zu vektorisieren ist:

Code: Alles auswählen

auto const x² = x * x;
auto const termA = (constants.b * x² + constants.a) * x;
auto const termB = (constants.d * x² + constants.c) * x² + constants.one;

if(0x00000000u == sectorID || 0x60000000u == sectorID) {
	return Float4B(termA / termB);
} else {
	return Float4B(termB / termA);
}

[BeRsErKeR]

Hab mal ne Verständnisfrage zu deinem Code. Die 4 Werte für die Sektoren sind ja 0x0..., 0x2..., 0x4... und 0x8... (laut Kommentar). Wenn du aber die Maske 0x6... nutzt kann ja niemals der Wert 0x8... für die SektorID bei rauskommen. Wäre das oberste Bit gesetzt so würde die Maske dieses ja bei Verundung eliminieren und du hättest wieder 0x0..., also den 1. Sektor. Und weiter unten prüfst du dann auf == 0x6... Versteh nicht so ganz was dahinter steckt. Objektiv sieht das aber merkwürdig aus. Wie stellst du dann fest ob der 4. Sektor gemeint ist?

[Krishty]

Scheiße! Das ist ein Fehler im Kommentar, den ich schon zwei Mal korrigiert habe. Durch meine Unart, nach nicht funktionierenden Optimierungen einfach ein Revert zu machen, ist er aber dringeblieben. Dort sollte 0x60000000 für den vierten Sektor stehen. Entschuldigung dafür.

[Krishty]

Ich habe von hier die Parameter meiner sin()-/cos()-Funktion zusammenkopiert. Sie sollten eigentlich auf 13 Stellen genau sein; trotzdem fallen mir zwei Anomalien auf:

1. Natürlich treten an denselben Stellen wie bei tan() Ungenauigkeiten bei der Umrechnung von Ganz- zu Gleitkommazahlen auf.
2. Etwa zwanzig Mal pro Grad weicht ein Wert eine ULP vom sin() / cos() der Standardbibliothek ab.

[Krishty]

Lustig: Mein sincos() ist bloß 15 % schneller als je ein Aufruf an das CRT-sin() und -cos() zusammen.

Aber jetzt kommt’s: Daraufhin habe ich mit Intrinsics alles in astreinen SSE2-Text verwandelt, der Sinus und Kosinus als Vektor zweier doubles schön vektorisiert berechnet –

– und jetzt dauert ein Aufruf so lange wie vier CRT-sin() zusammen. Jedoch: Ich initialisiere noch nicht ganz sauber; das läuft über den Stapel.


Seriell:

Code: Alles auswählen

CosineAndSine<Float4B> const cosineAndSineOf(
	Angle angle // pass by value: 4 % faster
) {
	CosineAndSine<Float4B> result;

	UnsignedInt4B sinUNorm;
	if(0x40000000 <= angle.unorm && 0xC0000000 > angle.unorm) {
		sinUNorm = 0x80000000u - angle.unorm;
	} else {
		sinUNorm = angle.unorm;
	}

	UnsignedInt4B cosUNorm;
	if(0x80000000 <= angle.unorm) {
		cosUNorm = 0x40000000u + angle.unorm;
	} else {
		cosUNorm = 0x40000000u - angle.unorm;
	}

	static double const
		a0 =  1.4629180792671596810513378043098e-9,   // ^1
		a1 = -5.218056424355326540518689143791e-28,   // -1.666666666640169148537065260055e-1  * 3.1308338546629715204060346522109e-27  // ^3
		a2 =  5.5836577275526616254281114471759e-47,  //  8.333333316490113523036717102793e-3  * 6.7003892866059294987769581698842e-45  // ^5
		a3 = -2.8451779180169008034005793100212e-66,  // -1.984126600659171392655484413285e-4  * 1.4339699478207029913663242067659e-62  // ^7
		a4 =  8.4568853391845024057023404189569e-86,  //  2.755690114917374804474016589137e-6  * 3.0688811101817481779794040453891e-80  // ^9
		a5 = -1.6438192896934818124147706229392e-105, // -2.502845227292692953118686710787e-8  * 6.5678024025144678446613914998034e-98  // ^11
		a6 =  2.16283153455215654880174036432e-125;   //  1.538730635926417598443354215485e-10 * 1.4055946401885921624309340128322e-115 // ^13


	__m128i unsignedNormalized;
	unsignedNormalized.m128i_i32[0] = cosUNorm;
	unsignedNormalized.m128i_i32[1] = sinUNorm;

	{
		auto sinAngle = double(int(sinUNorm));
		auto sinAngle² = sinAngle * sinAngle;
		auto x3 = sinAngle² * sinAngle;
		auto x4 = sinAngle² * sinAngle²;
		auto x8 = x4 * x4;
		auto x9 = x8 * sinAngle;
		auto A = x3 * (a1 + sinAngle² * (a2 + sinAngle² * a3));
		auto B = a4 + sinAngle² * (a5 + sinAngle² * a6);
		auto C = a0 * sinAngle;
		result.sine = Float4B(A + C + x9 * B);
	}

	{
		auto cosAngle = double(int(cosUNorm));
		double cosAngle² = cosAngle * cosAngle;
		auto x3 = cosAngle² * cosAngle;
		auto x4 = cosAngle² * cosAngle²;
		auto x8 = x4 * x4;
		auto x9 = x8 * cosAngle;
		auto A = x3 * (a1 + cosAngle² * (a2 + cosAngle² * a3));
		auto B = a4 + cosAngle² * (a5 + cosAngle² * a6);
		auto C = a0 * cosAngle;
		result.cosine = Float4B(A + C + x9 * B);
	}

	return result;
}

Und hier die vektorisierte Variante:

Code: Alles auswählen

	__m128i unsignedNormalized;
	unsignedNormalized.m128i_i32[0] = cosUNorm;
	unsignedNormalized.m128i_i32[1] = sinUNorm;

	auto const vx = _mm_cvtepi32_pd(unsignedNormalized);
	auto const vx² = _mm_mul_pd(vx, vx);
	auto const va0 = _mm_set1_pd(a0);
	auto const va1 = _mm_set1_pd(a1);
	auto const va2 = _mm_set1_pd(a2);
	auto const va3 = _mm_set1_pd(a3);
	auto const va4 = _mm_set1_pd(a4);
	auto const va5 = _mm_set1_pd(a5);
	auto const va6 = _mm_set1_pd(a6);
	auto const vx3 = _mm_mul_pd(vx², vx);
	auto const vx4 = _mm_mul_pd(vx², vx²);
	auto const vx8 = _mm_mul_pd(vx4, vx4);
	auto const vx9 = _mm_mul_pd(vx8, vx);
	auto const vA = _mm_mul_pd(vx3, _mm_add_pd(va1, _mm_mul_pd(vx², _mm_add_pd(va2, _mm_mul_pd(vx², va3)))));
	auto const vB = _mm_add_pd(va4, _mm_mul_pd(vx², _mm_add_pd(va5, _mm_mul_pd(vx², va6))));
	auto const vC = _mm_mul_pd(va0, vx);
	auto const resDouble = _mm_add_pd(vA, _mm_add_pd(vC, _mm_mul_pd(vx9, vB)));
	auto const resFloat = _mm_cvtpd_ps(resDouble);
	result.cosine = resFloat.m128_f32[0];
	result.sine = resFloat.m128_f32[1];

[Krishty]

Ich brauche mal eure Hilfe, weil ich allein nicht mehr wirklich weiterkomme. Das hier ist der Maschinentext der seriellen Version:

Code: Alles auswählen

Math::cosineAndSineOf:
 sub         rsp,38h  
 movaps      xmmword ptr [rsp+20h],xmm9  
 lea         eax,[rcx-40000000h]  
 movaps      xmmword ptr [rsp+10h],xmm10  
 movaps      xmmword ptr [rsp],xmm11  
 cmp         eax,7FFFFFFFh  
 ja          Math::cosineAndSineOf+2Bh (13F337A4Bh)  
 mov         edx,80000000h  
 sub         edx,ecx  
 jmp         Math::cosineAndSineOf+2Dh (13F337A4Dh)  
 mov         edx,ecx  
 cmp         ecx,80000000h  
 jb          Math::cosineAndSineOf+3Dh (13F337A5Dh)  
 lea         eax,[rcx+40000000h]  
 jmp         Math::cosineAndSineOf+44h (13F337A64h)  
 mov         eax,40000000h  
 sub         eax,ecx  
 movd        xmm9,edx  
 movd        xmm11,eax  
 cvtdq2pd    xmm9,xmm9  
 cvtdq2pd    xmm11,xmm11  
 movapd      xmm0,xmm9  
 movapd      xmm10,xmm11  
 mulsd       xmm0,xmm9  
 mulsd       xmm10,xmm11  
 movapd      xmm1,xmm0  
 mulsd       xmm1,mmword ptr [a6 (13F3553B8h)]  
 subsd       xmm1,mmword ptr [__real@2a2e2928c15c82bc (13F355740h)]  
 mulsd       xmm1,xmm0  
 addsd       xmm1,mmword ptr [a4 (13F3553B0h)]  
 mulsd       xmm1,xmm0  
 subsd       xmm1,mmword ptr [__real@32532d2c9034c37d (13F355738h)]  
 mulsd       xmm1,xmm0  
 addsd       xmm1,mmword ptr [a2 (13F3553A8h)]  
 mulsd       xmm1,xmm0  
 subsd       xmm1,mmword ptr [__real@3a44abbce62454fc (13F355730h)]  
 mulsd       xmm1,xmm0  
 addsd       xmm1,mmword ptr [a0 (13F3553A0h)]  
 mulsd       xmm1,xmm9  
 movaps      xmm9,xmmword ptr [rsp+20h]  
 cvtsd2ss    xmm0,xmm1  
 movapd      xmm1,xmm10  
 mulsd       xmm1,mmword ptr [a6 (13F3553B8h)]  
 subsd       xmm1,mmword ptr [__real@2a2e2928c15c82bc (13F355740h)]  
 movss       dword ptr [rsp+4Ch],xmm0  
 mulsd       xmm1,xmm10  
 addsd       xmm1,mmword ptr [a4 (13F3553B0h)]  
 mulsd       xmm1,xmm10  
 subsd       xmm1,mmword ptr [__real@32532d2c9034c37d (13F355738h)]  
 mulsd       xmm1,xmm10  
 addsd       xmm1,mmword ptr [a2 (13F3553A8h)]  
 mulsd       xmm1,xmm10  
 subsd       xmm1,mmword ptr [__real@3a44abbce62454fc (13F355730h)]  
 mulsd       xmm1,xmm10  
 movaps      xmm10,xmmword ptr [rsp+10h]  
 addsd       xmm1,mmword ptr [a0 (13F3553A0h)]  
 mulsd       xmm1,xmm11  
 movaps      xmm11,xmmword ptr [rsp]  
 cvtsd2ss    xmm0,xmm1  
 movss       dword ptr [rsp+48h],xmm0  
 mov         rax,qword ptr [rsp+48h]  
 add         rsp,38h  
 ret  

… und dies jener der vektorisierten Variante:

Code: Alles auswählen

Math::cosineAndSineOf:
 sub         rsp,18h  
 lea         eax,[rcx-40000000h]  
 cmp         eax,7FFFFFFFh  
 ja          Math::cosineAndSineOf+1Ah (13F627A3Ah)  
 mov         edx,80000000h  
 sub         edx,ecx  
 jmp         Math::cosineAndSineOf+1Ch (13F627A3Ch)  
 mov         edx,ecx  
 cmp         ecx,80000000h  
 jb          Math::cosineAndSineOf+2Ch (13F627A4Ch)  
 lea         eax,[rcx+40000000h]  
 jmp         Math::cosineAndSineOf+33h (13F627A53h)  
 mov         eax,40000000h  
 sub         eax,ecx  
 movapd      xmm2,xmmword ptr [va6 (13F645400h)]  
 movapd      xmm1,xmmword ptr [va3 (13F6453D0h)]  
 mov         dword ptr [rsp],eax  
 mov         dword ptr [rsp+4],edx  
 xor         ecx,ecx  
 mov         qword ptr [rsp+8],rcx  
 cvtdq2pd    xmm3,mmword ptr [rsp]  
 movapd      xmm4,xmm3  
 mulpd       xmm4,xmm3  
 mulpd       xmm2,xmm4  
 mulpd       xmm1,xmm4  
 movapd      xmm0,xmm4  
 addpd       xmm2,xmmword ptr [va5 (13F6453F0h)]  
 addpd       xmm1,xmmword ptr [va2 (13F6453C0h)]  
 mulpd       xmm0,xmm4  
 mulpd       xmm2,xmm4  
 mulpd       xmm0,xmm0  
 mulpd       xmm1,xmm4  
 addpd       xmm2,xmmword ptr [va4 (13F6453E0h)]  
 addpd       xmm1,xmmword ptr [va1 (13F6453B0h)]  
 mulpd       xmm0,xmm3  
 mulpd       xmm4,xmm3  
 mulpd       xmm2,xmm0  
 mulpd       xmm1,xmm4  
 movapd      xmm0,xmmword ptr [va0 (13F6453A0h)]  
 mulpd       xmm0,xmm3  
 addpd       xmm2,xmm0  
 addpd       xmm2,xmm1  
 cvtpd2ps    xmm0,xmm2  
 movaps      xmmword ptr [rsp],xmm0  
 mov         rax,qword ptr [rsp]  
 add         rsp,18h  
 ret  

Warum ist die vektorisierte Variante so deutlich langsamer (immernoch fast doppelt so lange Ausführungszeit)? Der einzige Grund, den ich erkennen könnte, ist, dass dort doppelt so viel über den Stapel wandert …

[dot]

Bist du immer noch auf deinem Atom unterwegs oder hast du das auch schonmal auf einem Core i7 getestet? iirc muss der Atom viele Vectorinstruktionen mangels genug ALU auf mehrere Zyklen aufteilen (könnte mir vorstellen, dass es bei double noch besonders schlimm is). Wohl kaum der einzige Grund, aber möglicherweise mitunter ein Grund!?

[Krishty]

Core i7.

[Krishty]

Seit wann benutzt DirectX Math denn ein Minimax-Polynom? Als ich mich Anfang des Jahres von XNA Math verabschiedet habe, wurden noch die Standard-sin()- und -cos()-Funktionen der CRT aufgerufen …

Interessant, dass sie einen 7. Grad als Annäherung und einen 11. als tatsächlichen Wert benutzen. Ich war glatt davon ausgegangen, dass der 7. Grad (je nach Quelle 12–13 dezimale Stellen) für float ausreichend wäre …

[Krishty]

Den Quelltext von Hand umzusortieren bringt Vorteile. WTF?! Bin ich im Jahr 1996 gelandet oder auf einen ATI-Shader-Compiler getreten oder was ist das für eine Scheiße?!

Ich hatte mich so gefreut, dass VC++ mit Gleitkommaarithmetik für x64 endlich mal stabilen Maschinentext erzeugt, und jetzt passiert sowas. Nur, damit ihr eine Vorstellung bekommt: Ob die Funktion 9 oder 14 Register belegt ergibt sich daraus, in welcher Reihenfolge ich die voneinander völlig unabhängigen Sinus- und Kosinusberechnungen anordne. Ich wette, wenn ich alle Permutationen der Zeilen durchginge, würde ich eine noch viel effizientere Version finden. Bah.

[eXile]

Den Quelltext von Hand umzusortieren bringt Vorteile. WTF?! Bin ich im Jahr 1996 gelandet oder auf einen ATI-Shader-Compiler getreten oder was ist das für eine Scheiße?!

Richtig, deswegen sagen sowohl die Präsentationen von Green wie auch der zitierte Blog-Eintrag von dir: Einfach haufenweise Permutationen ausprobieren.

[Krishty]

Sooo … ich habe nun auch ein auf 0.34 ULPs genaues atan(). Geklaut ist es von http://www.ganssle.com/approx/approx.pdf (der Web-Version des Artikels fehlt der Range Reduction-Text) und so angepasst, dass es meine Ganzzahl-Winkel aussppuckt. Die Genauigkeit ist also 0,34 ÷ 2³² Vollkreise – oder 0,0000000285°.

Das Branching dadrin ist wirklich widerlich (kein Wunder, wo Herr Ganssle doch so ein x86-Gegner ist). Herr Green führt eine Range Reduction auf 8 statt 12 Sektoren durch, und nähert darin mit viel einfacheren Polynomen an (keine Divisionen!); darum gehe ich davon aus, dass man die Leistung noch deutlich steigern kann. Weil ich aber 11 Stellen Genauigkeit brauche statt den 6, die Herr Green produziert, und der Vergleich meiner Ganzzahlwinkel mit den Gleitkommawinkeln der Standardbibliothek quälend schwierig ist, fahre ich erstmal die sichere Schiene.

Ich werde nun zusehen, dass ich ein atan2() schreibe, das die atan2()-Range Reduction mit der von atan() kombiniert. Wenn ich das geschickt mache, kann ich noch mehr Genauigkeit rausholen; und dann gibt es auch wieder Quelltext.