[C++] Schnelleres sin/cos/tan

Design Patterns, Erklärungen zu Algorithmen, Optimierung, Softwarearchitektur
Forumsregeln
Wenn das Problem mit einer Programmiersprache direkt zusammenhängt, bitte HIER posten.
Benutzeravatar
Krishty
Establishment
Beiträge: 8229
Registriert: 26.02.2009, 11:18
Benutzertext: state is the enemy
Kontaktdaten:

Re: [C++] Schnelleres sin/cos/tan

Beitrag von 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.
Zuletzt geändert von Krishty am 11.09.2012, 17:07, insgesamt 1-mal geändert.
seziert Ace Combat, Driver, und S.T.A.L.K.E.R.   —   rendert Sterne
Benutzeravatar
Krishty
Establishment
Beiträge: 8229
Registriert: 26.02.2009, 11:18
Benutzertext: state is the enemy
Kontaktdaten:

Re: [C++] Schnelleres sin/cos/tan

Beitrag von 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.
seziert Ace Combat, Driver, und S.T.A.L.K.E.R.   —   rendert Sterne
Benutzeravatar
BeRsErKeR
Establishment
Beiträge: 689
Registriert: 27.04.2002, 22:01

Re: [C++] Schnelleres sin/cos/tan

Beitrag von BeRsErKeR »

Krishty hat geschrieben: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.
Ohne Input kein Output.
Benutzeravatar
Krishty
Establishment
Beiträge: 8229
Registriert: 26.02.2009, 11:18
Benutzertext: state is the enemy
Kontaktdaten:

Re: [C++] Schnelleres sin/cos/tan

Beitrag von Krishty »

BeRsErKeR hat geschrieben:
Krishty hat geschrieben: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);
}
seziert Ace Combat, Driver, und S.T.A.L.K.E.R.   —   rendert Sterne
Benutzeravatar
BeRsErKeR
Establishment
Beiträge: 689
Registriert: 27.04.2002, 22:01

Re: [C++] Schnelleres sin/cos/tan

Beitrag von 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?
Ohne Input kein Output.
Benutzeravatar
Krishty
Establishment
Beiträge: 8229
Registriert: 26.02.2009, 11:18
Benutzertext: state is the enemy
Kontaktdaten:

Re: [C++] Schnelleres sin/cos/tan

Beitrag von 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.
seziert Ace Combat, Driver, und S.T.A.L.K.E.R.   —   rendert Sterne
Benutzeravatar
Krishty
Establishment
Beiträge: 8229
Registriert: 26.02.2009, 11:18
Benutzertext: state is the enemy
Kontaktdaten:

Re: [C++] Schnelleres sin/cos/tan

Beitrag von 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.
seziert Ace Combat, Driver, und S.T.A.L.K.E.R.   —   rendert Sterne
Benutzeravatar
Krishty
Establishment
Beiträge: 8229
Registriert: 26.02.2009, 11:18
Benutzertext: state is the enemy
Kontaktdaten:

Re: [C++] Schnelleres sin/cos/tan

Beitrag von 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];
seziert Ace Combat, Driver, und S.T.A.L.K.E.R.   —   rendert Sterne
Benutzeravatar
Krishty
Establishment
Beiträge: 8229
Registriert: 26.02.2009, 11:18
Benutzertext: state is the enemy
Kontaktdaten:

Re: [C++] Schnelleres sin/cos/tan

Beitrag von 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 …
seziert Ace Combat, Driver, und S.T.A.L.K.E.R.   —   rendert Sterne
Benutzeravatar
dot
Establishment
Beiträge: 1734
Registriert: 06.03.2004, 18:10
Echter Name: Michael Kenzel
Kontaktdaten:

Re: [C++] Schnelleres sin/cos/tan

Beitrag von 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!?
Zuletzt geändert von dot am 12.09.2012, 17:50, insgesamt 1-mal geändert.
Benutzeravatar
Krishty
Establishment
Beiträge: 8229
Registriert: 26.02.2009, 11:18
Benutzertext: state is the enemy
Kontaktdaten:

Re: [C++] Schnelleres sin/cos/tan

Beitrag von Krishty »

Core i7.
seziert Ace Combat, Driver, und S.T.A.L.K.E.R.   —   rendert Sterne
Benutzeravatar
Krishty
Establishment
Beiträge: 8229
Registriert: 26.02.2009, 11:18
Benutzertext: state is the enemy
Kontaktdaten:

Re: [C++] Schnelleres sin/cos/tan

Beitrag von 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 …
seziert Ace Combat, Driver, und S.T.A.L.K.E.R.   —   rendert Sterne
Benutzeravatar
Krishty
Establishment
Beiträge: 8229
Registriert: 26.02.2009, 11:18
Benutzertext: state is the enemy
Kontaktdaten:

Re: [C++] Schnelleres sin/cos/tan

Beitrag von 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.
seziert Ace Combat, Driver, und S.T.A.L.K.E.R.   —   rendert Sterne
Benutzeravatar
eXile
Establishment
Beiträge: 1136
Registriert: 28.02.2009, 13:27

Re: [C++] Schnelleres sin/cos/tan

Beitrag von eXile »

Krishty hat geschrieben: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.
Benutzeravatar
Krishty
Establishment
Beiträge: 8229
Registriert: 26.02.2009, 11:18
Benutzertext: state is the enemy
Kontaktdaten:

Re: [C++] Schnelleres sin/cos/tan

Beitrag von 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.
seziert Ace Combat, Driver, und S.T.A.L.K.E.R.   —   rendert Sterne
Benutzeravatar
Krishty
Establishment
Beiträge: 8229
Registriert: 26.02.2009, 11:18
Benutzertext: state is the enemy
Kontaktdaten:

Re: [C++] Schnelleres sin/cos/tan

Beitrag von Krishty »

Krishty hat geschrieben: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.
Weil ich in letzter Zeit viel mit SSE zu tun habe, habe ich mich dem Problem nochmal angenommen. Die SSE-optimierte Variante, die sincos() vektorisiert ausrechnet, ist tatsächlich gut 30 % schneller als die skalare Variante. Und 2,3× so schnell wie cos() und sin() der Visual C++ 2013-CRT. Ich habe es diesmal wieder auf einem Core i7 getestet, jedoch auf was aus der Ivy Bridge-Ära (der Core i7 der früheren Tests hat mittlerweile ein halbes Jahrzehnt auf dem Buckel). Also sowohl eigenes sincos() als auch SSE-vektorisierung lohnen sich. Der Quelltext, den ich früher gepostet habe, sieht auch völlig in Ordnung aus.

Mit 40 % Zeitanteil ist der langsamste Teil der Funktion jetzt übrigens das Kopieren des Rückgabewerts – wegen diesem „keine NRVO mit float“-Bug in Visual C++. Bin gespannt, ob sie den im 2015er Compiler behoben haben.
seziert Ace Combat, Driver, und S.T.A.L.K.E.R.   —   rendert Sterne
Antworten