[C++] Schnelleres sin/cos/tan
Forumsregeln
Wenn das Problem mit einer Programmiersprache direkt zusammenhängt, bitte HIER posten.
Wenn das Problem mit einer Programmiersprache direkt zusammenhängt, bitte HIER posten.
- dot
- Establishment
- Beiträge: 1734
- Registriert: 06.03.2004, 18:10
- Echter Name: Michael Kenzel
- Kontaktdaten:
Re: [C++] Schnelleres sin/cos/tan
Naja, eben als die fptan Instruction der FPU...
EDIT: Ok, wie es aussieht lässt sich MSVC nicht dazu bewegen, fptan zu generieren...
EDIT: Ok, wie es aussieht lässt sich MSVC nicht dazu bewegen, fptan zu generieren...
- Krishty
- Establishment
- Beiträge: 8238
- Registriert: 26.02.2009, 11:18
- Benutzertext: state is the enemy
- Kontaktdaten:
Re: [C++] Schnelleres sin/cos/tan
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.
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
- Establishment
- Beiträge: 1734
- Registriert: 06.03.2004, 18:10
- Echter Name: Michael Kenzel
- Kontaktdaten:
Re: [C++] Schnelleres sin/cos/tan
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
- Establishment
- Beiträge: 8238
- Registriert: 26.02.2009, 11:18
- Benutzertext: state is the enemy
- Kontaktdaten:
Re: [C++] Schnelleres sin/cos/tan
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
- Establishment
- Beiträge: 8238
- Registriert: 26.02.2009, 11:18
- Benutzertext: state is the enemy
- Kontaktdaten:
Re: [C++] Schnelleres sin/cos/tan
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).
Re: Jammer-Thread
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.Matthias Gubisch hat geschrieben: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
Ohne Input kein Output.
- Artificial Mind
- Establishment
- Beiträge: 802
- Registriert: 17.12.2007, 17:51
- Wohnort: Aachen
Re: [C++] Schnelleres sin/cos/tan
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?)).
Re: [C++] Schnelleres sin/cos/tan
Naja aber das klang irgendwie so, als ob Krishty irgendeine Lib baut oder ähnliches.
Ohne Input kein Output.
-
- Establishment
- Beiträge: 470
- Registriert: 01.03.2009, 19:09
Re: [C++] Schnelleres sin/cos/tan
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 ;)
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 ;)
Bevor man den Kopf schüttelt, sollte man sich vergewissern einen zu haben
- dot
- Establishment
- Beiträge: 1734
- Registriert: 06.03.2004, 18:10
- Echter Name: Michael Kenzel
- Kontaktdaten:
Re: [C++] Schnelleres sin/cos/tan
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
- Establishment
- Beiträge: 8238
- Registriert: 26.02.2009, 11:18
- Benutzertext: state is the enemy
- Kontaktdaten:
Re: [C++] Schnelleres sin/cos/tan
Der Kram fließt einfach so in mein Framework ein. Es geht mir in erster Linie darum,
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.
, dann umKrishty hat geschrieben:es richtig zu machen statt so, wie es alle tun
und schließlich umAbnabelung von der C++-Laufzeitbibliothek
.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
- Establishment
- Beiträge: 1734
- Registriert: 06.03.2004, 18:10
- Echter Name: Michael Kenzel
- Kontaktdaten:
Re: [C++] Schnelleres sin/cos/tan
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?
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
- Establishment
- Beiträge: 8238
- Registriert: 26.02.2009, 11:18
- Benutzertext: state is the enemy
- Kontaktdaten:
Re: [C++] Schnelleres sin/cos/tan
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).
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
- Establishment
- Beiträge: 8238
- Registriert: 26.02.2009, 11:18
- Benutzertext: state is the enemy
- Kontaktdaten:
Re: [C++] Schnelleres sin/cos/tan
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.
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
- Establishment
- Beiträge: 8238
- Registriert: 26.02.2009, 11:18
- Benutzertext: state is the enemy
- Kontaktdaten:
Re: [C++] Schnelleres sin/cos/tan
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.
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
- Establishment
- Beiträge: 8238
- Registriert: 26.02.2009, 11:18
- Benutzertext: state is the enemy
- Kontaktdaten:
Re: [C++] Schnelleres sin/cos/tan
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);
}
}
Zuletzt geändert von Krishty am 11.09.2012, 17:07, insgesamt 1-mal geändert.
- Krishty
- Establishment
- Beiträge: 8238
- Registriert: 26.02.2009, 11:18
- Benutzertext: state is the enemy
- Kontaktdaten:
Re: [C++] Schnelleres sin/cos/tan
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.
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.
Re: [C++] Schnelleres sin/cos/tan
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 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.
Ohne Input kein Output.
- Krishty
- Establishment
- Beiträge: 8238
- Registriert: 26.02.2009, 11:18
- Benutzertext: state is the enemy
- Kontaktdaten:
Re: [C++] Schnelleres sin/cos/tan
Zuerst: 1/3 und 1/10 sind nicht genau darstellbar (weil sie binär keine endlichen Brüche sind), 1/2 schon.BeRsErKeR hat geschrieben: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 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.
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);
}
Re: [C++] Schnelleres sin/cos/tan
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.
- Krishty
- Establishment
- Beiträge: 8238
- Registriert: 26.02.2009, 11:18
- Benutzertext: state is the enemy
- Kontaktdaten:
Re: [C++] Schnelleres sin/cos/tan
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
- Establishment
- Beiträge: 8238
- Registriert: 26.02.2009, 11:18
- Benutzertext: state is the enemy
- Kontaktdaten:
Re: [C++] Schnelleres sin/cos/tan
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:
- Natürlich treten an denselben Stellen wie bei tan() Ungenauigkeiten bei der Umrechnung von Ganz- zu Gleitkommazahlen auf.
- Etwa zwanzig Mal pro Grad weicht ein Wert eine ULP vom sin() / cos() der Standardbibliothek ab.
- Krishty
- Establishment
- Beiträge: 8238
- Registriert: 26.02.2009, 11:18
- Benutzertext: state is the enemy
- Kontaktdaten:
Re: [C++] Schnelleres sin/cos/tan
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:Und hier die vektorisierte Variante:
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;
}
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
- Establishment
- Beiträge: 8238
- Registriert: 26.02.2009, 11:18
- Benutzertext: state is the enemy
- Kontaktdaten:
Re: [C++] Schnelleres sin/cos/tan
Ich brauche mal eure Hilfe, weil ich allein nicht mehr wirklich weiterkomme. Das hier ist der Maschinentext der seriellen Version:
… und dies jener der vektorisierten Variante:
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 …
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
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
- dot
- Establishment
- Beiträge: 1734
- Registriert: 06.03.2004, 18:10
- Echter Name: Michael Kenzel
- Kontaktdaten:
Re: [C++] Schnelleres sin/cos/tan
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.
- Krishty
- Establishment
- Beiträge: 8238
- Registriert: 26.02.2009, 11:18
- Benutzertext: state is the enemy
- Kontaktdaten:
Re: [C++] Schnelleres sin/cos/tan
Core i7.
- Krishty
- Establishment
- Beiträge: 8238
- Registriert: 26.02.2009, 11:18
- Benutzertext: state is the enemy
- Kontaktdaten:
Re: [C++] Schnelleres sin/cos/tan
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 …
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
- Establishment
- Beiträge: 8238
- Registriert: 26.02.2009, 11:18
- Benutzertext: state is the enemy
- Kontaktdaten:
Re: [C++] Schnelleres sin/cos/tan
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.
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.
Re: [C++] Schnelleres sin/cos/tan
Richtig, deswegen sagen sowohl die Präsentationen von Green wie auch der zitierte Blog-Eintrag von dir: Einfach haufenweise Permutationen ausprobieren.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?!
- Krishty
- Establishment
- Beiträge: 8238
- Registriert: 26.02.2009, 11:18
- Benutzertext: state is the enemy
- Kontaktdaten:
Re: [C++] Schnelleres sin/cos/tan
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.
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.