Warum ist SSE skalar sqrt (x) langsamer als rsqrt (x) * x?

106

Ich habe einige unserer Kernmathematiken auf einem Intel Core Duo profiliert und bei der Betrachtung verschiedener Ansätze zur Quadratwurzel etwas Seltsames festgestellt: Mit den skalaren SSE-Operationen ist es schneller, eine reziproke Quadratwurzel zu nehmen und zu multiplizieren um das sqrt zu erhalten, muss der native sqrt-opcode verwendet werden!

Ich teste es mit einer Schleife wie:

inline float TestSqrtFunction( float in );

void TestFunc()
{
  #define ARRAYSIZE 4096
  #define NUMITERS 16386
  float flIn[ ARRAYSIZE ]; // filled with random numbers ( 0 .. 2^22 )
  float flOut [ ARRAYSIZE ]; // filled with 0 to force fetch into L1 cache

  cyclecounter.Start();
  for ( int i = 0 ; i < NUMITERS ; ++i )
    for ( int j = 0 ; j < ARRAYSIZE ; ++j )
    {
       flOut[j] = TestSqrtFunction( flIn[j] );
       // unrolling this loop makes no difference -- I tested it.
    }
  cyclecounter.Stop();
  printf( "%d loops over %d floats took %.3f milliseconds",
          NUMITERS, ARRAYSIZE, cyclecounter.Milliseconds() );
}

Ich habe dies mit ein paar verschiedenen Körpern für die TestSqrtFunction versucht, und ich habe einige Timings, die mir wirklich am Kopf kratzen. Das mit Abstand Schlimmste war, die native Funktion sqrt () zu verwenden und den "intelligenten" Compiler "optimieren" zu lassen. Bei 24 ns / float war dies mit der x87-FPU erbärmlich schlecht:

inline float TestSqrtFunction( float in )
{  return sqrt(in); }

Das nächste, was ich versuchte, war die Verwendung eines Intrinsic, um den Compiler zu zwingen, den skalaren Sqrt-Opcode von SSE zu verwenden:

inline void SSESqrt( float * restrict pOut, float * restrict pIn )
{
   _mm_store_ss( pOut, _mm_sqrt_ss( _mm_load_ss( pIn ) ) );
   // compiles to movss, sqrtss, movss
}

Dies war besser bei 11,9 ns / float. Ich habe auch Carmacks verrückte Newton-Raphson-Approximationstechnik ausprobiert , die mit 4,3 ns / float sogar noch besser lief als die Hardware, allerdings mit einem Fehler von 1 zu 2 10 (was für meine Zwecke zu viel ist).

Der Trottel war, als ich die SSE-Operation für die reziproke Quadratwurzel ausprobierte und dann eine Multiplikation verwendete, um die Quadratwurzel zu erhalten (x * 1 / √x = √x). Auch wenn diese beiden abhängigen Operationen nimmt, war es die schnellste Lösung bei weitem, bei 1.24ns / Schwimmer und genau 2 -14 :

inline void SSESqrt_Recip_Times_X( float * restrict pOut, float * restrict pIn )
{
   __m128 in = _mm_load_ss( pIn );
   _mm_store_ss( pOut, _mm_mul_ss( in, _mm_rsqrt_ss( in ) ) );
   // compiles to movss, movaps, rsqrtss, mulss, movss
}

Meine Frage ist im Grunde, was gibt ? Warum ist der in die Hardware integrierte Quadratwurzel-Opcode von SSE langsamer als die Synthese aus zwei anderen mathematischen Operationen?

Ich bin mir sicher, dass dies wirklich die Kosten für die Operation selbst sind, da ich Folgendes überprüft habe:

  • Alle Daten passen in den Cache und die Zugriffe erfolgen nacheinander
  • Die Funktionen sind inline
  • Das Abrollen der Schleife macht keinen Unterschied
  • Compiler-Flags sind auf vollständige Optimierung gesetzt (und die Assembly ist gut, habe ich überprüft)

( edit : stephentyrone weist korrekt darauf hin, dass Operationen an langen Zahlenfolgen die vektorisierenden SIMD-gepackten Operationen verwenden sollten, wie rsqrtps- aber die Array-Datenstruktur hier dient nur zu Testzwecken: Was ich wirklich zu messen versuche, ist die skalare Leistung für die Verwendung in Code das kann nicht vektorisiert werden.)

Crashworks
quelle
13
x / sqrt (x) = sqrt (x). Oder anders ausgedrückt: x ^ 1 * x ^ (- 1/2) = x ^ (1 - 1/2) = x ^ (1/2) = sqrt (x)
Crashworks
6
natürlich , inline float SSESqrt( float restrict fIn ) { float fOut; _mm_store_ss( &fOut, _mm_sqrt_ss( _mm_load_ss( &fIn ) ) ); return fOut; }. Dies ist jedoch eine schlechte Idee, da es leicht zu einem Load-Hit-Store-Stall kommen kann, wenn die CPU die Floats in den Stack schreibt und sie dann sofort zurückliest - insbesondere für den Rückgabewert vom Vektorregister in ein Float-Register jonglieren ist eine schlechte Nachricht. Außerdem nehmen die zugrunde liegenden Maschinen-Opcodes, die die SSE-Intrinsics darstellen, ohnehin Adressoperanden an.
Crashworks
4
Wie wichtig LHS ist, hängt von der jeweiligen Generation und dem Schritt eines bestimmten x86 ab: Ich habe die Erfahrung gemacht, dass das Verschieben von Daten zwischen Registersätzen (z. B. FPU zu SSE zu eax) bei i7 bis zu einem sehr schlechten Wert ist, während ein Roundtrip zwischen xmm0 und Stack erfolgt und zurück nicht, wegen Intels Store-Weiterleitung. Sie können es selbst zeitlich festlegen, um sicher zu sehen. Im Allgemeinen ist es am einfachsten, potenzielle LHS zu erkennen, indem Sie sich die emittierte Baugruppe ansehen und feststellen, wo Daten zwischen Registersätzen jongliert werden. Ihr Compiler macht möglicherweise das Schlaue oder nicht. Was die Normalisierung von Vektoren betrifft
Crashworks
2
Für den PowerPC ja: IBM verfügt über einen CPU-Simulator, der LHS und viele andere Pipeline-Blasen durch statische Analyse vorhersagen kann. Einige PPCs verfügen auch über einen Hardware-Zähler für LHS, den Sie abrufen können. Es ist schwieriger für das x86; Gute Profiling-Tools sind seltener (VTune ist heutzutage etwas kaputt) und die neu geordneten Pipelines sind weniger deterministisch. Sie können versuchen, es empirisch zu messen, indem Sie Anweisungen pro Zyklus messen, was genau mit den Hardware-Leistungsindikatoren erfolgen kann. Die Register "Anweisungen zurückgezogen" und "Gesamtzyklen" können beispielsweise mit PAPI oder PerfSuite ( bit.ly/an6cMt ) gelesen werden .
Crashworks
2
Sie können auch einfach einige Permutationen auf eine Funktion schreiben und sie zeitlich festlegen, um festzustellen, ob sie besonders unter Ständen leiden. Intel veröffentlicht nicht viele Details über die Funktionsweise ihrer Pipelines (dass LHS überhaupt ein schmutziges Geheimnis ist). Ich habe also viel gelernt, indem ich mir ein Szenario angesehen habe, das andere Bögen (z. B. PPC) zum Stillstand bringt ) und dann ein kontrolliertes Experiment erstellen, um zu sehen, ob das x86 es auch hat.
Crashworks

Antworten:

216

sqrtssergibt ein korrekt gerundetes Ergebnis. rsqrtssgibt eine Annäherung an das Reziproke, genau auf ungefähr 11 Bits.

sqrtssgeneriert ein weitaus genaueres Ergebnis, wenn Genauigkeit erforderlich ist. rsqrtssexistiert für die Fälle, in denen eine Annäherung ausreicht, aber Geschwindigkeit erforderlich ist. Wenn Sie die Dokumentation von Intel lesen, finden Sie auch eine Befehlssequenz (reziproke Quadratwurzel-Approximation, gefolgt von einem einzelnen Newton-Raphson-Schritt), die nahezu vollständige Genauigkeit liefert (~ 23 Bit Genauigkeit, wenn ich mich richtig erinnere) und immer noch etwas ist schneller als sqrtss.

Bearbeiten: Wenn die Geschwindigkeit kritisch ist und Sie dies für viele Werte wirklich in einer Schleife aufrufen, sollten Sie die vektorisierten Versionen dieser Anweisungen verwenden rsqrtpsoder sqrtpsbeide vier Floats pro Anweisung verarbeiten.

Stephen Canon
quelle
3
Der n / r-Schritt gibt Ihnen eine Genauigkeit von 22 Bit (er verdoppelt ihn). 23 Bit wären genau die volle Genauigkeit.
Jasper Bekkers
7
@ Jasper Bekkers: Nein, das würde es nicht. Erstens hat float eine Genauigkeit von 24 Bit. Zweitens sqrtsswird korrekt gerundet , was vor dem Runden ~ 50 Bit erfordert, und kann nicht mit einer einfachen N / R-Iteration mit einfacher Genauigkeit erreicht werden.
Stephen Canon
1
Dies ist definitiv der Grund. Um dieses Ergebnis zu erweitern: Intels Embree-Projekt ( software.intel.com/en-us/articles/… ) verwendet für seine Mathematik die Vektorisierung. Sie können die Quelle unter diesem Link herunterladen und sehen, wie sie ihre 3/4 D-Vektoren machen. Ihre Vektornormalisierung verwendet rsqrt, gefolgt von einer Iteration von Newton-Raphson, die dann sehr genau und immer noch schneller als 1 / ssqrt ist!
Brandon Pelfrey
7
Eine kleine Einschränkung: x rsqrt (x) führt zu NaN, wenn x entweder Null oder unendlich ist. 0 * rsqrt (0) = 0 * INF = NaN. INF rsqrt (INF) = INF * 0 = NaN. Aus diesem Grund berechnet CUDA auf NVIDIA-GPUs ungefähre Quadratwurzeln mit einfacher Genauigkeit als Rezept (rsqrt (x)), wobei die Hardware sowohl eine schnelle Annäherung an die reziproke als auch an die reziproke Quadratwurzel bietet. Natürlich sind auch explizite Überprüfungen der beiden Sonderfälle möglich (auf der GPU wären sie jedoch langsamer).
Njuffa
@BrandonPelfrey In welcher Datei haben Sie den Newton Rhapson-Schritt gefunden?
Fredoverflow
7

Dies gilt auch für die Teilung. MULSS (a, RCPSS (b)) ist viel schneller als DIVSS (a, b). Tatsächlich ist es immer noch schneller, selbst wenn Sie die Präzision mit einer Newton-Raphson-Iteration erhöhen.

Intel und AMD empfehlen diese Technik in ihren Optimierungshandbüchern. In Anwendungen, für die keine IEEE-754-Konformität erforderlich ist, ist der einzige Grund für die Verwendung von div / sqrt die Lesbarkeit des Codes.

Spat
quelle
1
Broadwell und später haben eine bessere FP-Divide-Leistung, daher verwenden Compiler wie clang bei neueren CPUs kein reziprokes + Newton für Skalar, da es normalerweise nicht schneller ist. In den meisten Schleifen divist dies nicht die einzige Operation, sodass der gesamte UOP-Durchsatz häufig der Engpass ist, selbst wenn ein divpsoder vorhanden ist divss. Siehe Gleitkommadivision vs. Gleitkommamultiplikation , wo meine Antwort einen Abschnitt darüber enthält, warum rcppskein Durchsatzgewinn mehr ist. (Oder ein Latenzgewinn) und Zahlen zum Teilungsdurchsatz / zur Latenz.
Peter Cordes
Wenn Ihre Genauigkeitsanforderungen so niedrig sind, dass Sie eine Newton-Iteration überspringen können, a * rcpss(b)kann dies zwar schneller sein, aber es sind immer noch mehr Ups als a/b!
Peter Cordes
5

Anstatt eine Antwort zu geben, die möglicherweise falsch ist (ich werde auch nicht nach Cache und anderen Dingen suchen oder darüber streiten, sagen wir, sie sind identisch), werde ich versuchen, Sie auf die Quelle zu verweisen, die Ihre Frage beantworten kann.
Der Unterschied könnte darin liegen, wie sqrt und rsqrt berechnet werden. Weitere Informationen finden Sie hier http://www.intel.com/products/processor/manuals/ . Ich würde vorschlagen, zunächst die von Ihnen verwendeten Prozessorfunktionen zu lesen. Es gibt einige Informationen, insbesondere zu rsqrt (die CPU verwendet eine interne Nachschlagetabelle mit großer Annäherung, wodurch es viel einfacher wird, das Ergebnis zu erhalten). Es scheint, dass rsqrt so viel schneller als sqrt ist, dass 1 zusätzliche Mehrfachoperation (die nicht zu kostspielig ist) die Situation hier möglicherweise nicht ändert.

Bearbeiten:
Einige Fakten, die erwähnenswert sein könnten: 1. Nachdem ich einige Mikrooptimierungen für meine Grafikbibliothek durchgeführt und rsqrt zur Berechnung der Länge von Vektoren verwendet habe. (Anstelle von sqrt habe ich meine Quadratsumme mit rsqrt multipliziert, was genau das ist, was Sie in Ihren Tests getan haben), und es hat eine bessere Leistung erbracht.
2. Das Berechnen von rsqrt mithilfe einer einfachen Nachschlagetabelle ist möglicherweise einfacher, da bei rsqrt, wenn x gegen unendlich geht, 1 / sqrt (x) auf 0 geht, sodass sich bei kleinen x die Funktionswerte nicht ändern (viel), während bei sqrt - es geht bis ins Unendliche, also ist es dieser einfache Fall;).

Klarstellung: Ich bin nicht sicher, wo ich es in Büchern gefunden habe, die ich verlinkt habe, aber ich bin mir ziemlich sicher, dass ich gelesen habe, dass rsqrt eine Nachschlagetabelle verwendet, und es sollte nur verwendet werden, wenn das Ergebnis vorliegt muss aber nicht genau sein - ich könnte mich auch irren, wie es vor einiger Zeit war :).

Marcin Deptuła
quelle
4

Newton-Raphson konvergiert gegen die Null f(x), -f/f' wenn Inkremente verwendet werden, f'die der Ableitung entsprechen.

Für x=sqrt(y), können Sie versuchen zu lösen , f(x) = 0für xVerwendung f(x) = x^2 - y;

Dann ist das Inkrement: dx = -f/f' = 1/2 (x - y/x) = 1/2 (x^2 - y) / x das eine langsame Teilung hat.

Sie können andere Funktionen (wie f(x) = 1/y - 1/x^2) ausprobieren , diese sind jedoch ebenso kompliziert.

Schauen wir uns 1/sqrt(y)jetzt an. Sie können es versuchen f(x) = x^2 - 1/y, aber es wird genauso kompliziert sein: dx = 2xy / (y*x^2 - 1)zum Beispiel. Eine nicht offensichtliche alternative Wahl für f(x)ist:f(x) = y - 1/x^2

Dann: dx = -f/f' = (y - 1/x^2) / (2/x^3) = 1/2 * x * (1 - y * x^2)

Ah! Es ist kein trivialer Ausdruck, aber Sie haben nur Multiplikationen, keine Teilung. => Schneller!

Und: Der vollständige Aktualisierungsschritt new_x = x + dxlautet dann:

x *= 3/2 - y/2 * x * x das ist auch einfach.

skal
quelle
2

Es gibt bereits vor einigen Jahren eine Reihe anderer Antworten darauf. Hier ist, was der Konsens richtig gemacht hat:

  • Die rsqrt * -Anweisungen berechnen eine Annäherung an die reziproke Quadratwurzel, die gut zu ungefähr 11-12 Bits ist.
  • Es wird mit einer Nachschlagetabelle (dh einem ROM) implementiert, die von der Mantisse indiziert wird. (Tatsächlich handelt es sich um eine komprimierte Nachschlagetabelle, die den alten mathematischen Tabellen ähnelt und Anpassungen an den niederwertigen Bits verwendet, um Transistoren zu sparen.)
  • Der Grund, warum es verfügbar ist, ist, dass es die anfängliche Schätzung ist, die von der FPU für den "echten" Quadratwurzel-Algorithmus verwendet wird.
  • Es gibt auch eine ungefähre gegenseitige Anweisung, rcp. Beide Anweisungen sind ein Hinweis darauf, wie die FPU Quadratwurzel und Division implementiert.

Hier ist, was der Konsens falsch gemacht hat:

  • FPUs aus der SSE-Ära verwenden Newton-Raphson nicht zur Berechnung von Quadratwurzeln. Es ist eine großartige Methode in Software, aber es wäre ein Fehler, sie auf diese Weise in Hardware zu implementieren.

Der NR-Algorithmus zum Berechnen der reziproken Quadratwurzel hat diesen Aktualisierungsschritt, wie andere angemerkt haben:

x' = 0.5 * x * (3 - n*x*x);

Das sind viele datenabhängige Multiplikationen und eine Subtraktion.

Was folgt, ist der Algorithmus, den moderne FPUs tatsächlich verwenden.

In Anbetracht b[0] = n, nehmen wir eine Reihe von Zahlen finden können , Y[i]so dass b[n] = b[0] * Y[0]^2 * Y[1]^2 * ... * Y[n]^2Ansätze 1. Dann betrachten:

x[n] = b[0] * Y[0] * Y[1] * ... * Y[n]
y[n] = Y[0] * Y[1] * ... * Y[n]

Klar x[n]Ansätze sqrt(n)und y[n]Ansätze1/sqrt(n) .

Wir können den Newton-Raphson-Aktualisierungsschritt für die reziproke Quadratwurzel verwenden, um ein gutes Ergebnis zu erhalten Y[i]:

b[i] = b[i-1] * Y[i-1]^2
Y[i] = 0.5 * (3 - b[i])

Dann:

x[0] = n Y[0]
x[i] = x[i-1] * Y[i]

und:

y[0] = Y[0]
y[i] = y[i-1] * Y[i]

Die nächste wichtige Beobachtung ist die folgende b[i] = x[i-1] * y[i-1]. So:

Y[i] = 0.5 * (3 - x[i-1] * y[i-1])
     = 1 + 0.5 * (1 - x[i-1] * y[i-1])

Dann:

x[i] = x[i-1] * (1 + 0.5 * (1 - x[i-1] * y[i-1]))
     = x[i-1] + x[i-1] * 0.5 * (1 - x[i-1] * y[i-1]))
y[i] = y[i-1] * (1 + 0.5 * (1 - x[i-1] * y[i-1]))
     = y[i-1] + y[i-1] * 0.5 * (1 - x[i-1] * y[i-1]))

Das heißt, bei anfänglichem x und y können wir den folgenden Aktualisierungsschritt verwenden:

r = 0.5 * (1 - x * y)
x' = x + x * r
y' = y + y * r

Oder, noch schicker, wir können setzen h = 0.5 * y. Dies ist die Initialisierung:

Y = approx_rsqrt(n)
x = Y * n
h = Y * 0.5

Und das ist der Update-Schritt:

r = 0.5 - x * h
x' = x + x * r
h' = h + h * r

Dies ist der Goldschmidt-Algorithmus, und er hat einen großen Vorteil, wenn Sie ihn in Hardware implementieren: Die "innere Schleife" besteht aus drei Multiplikationsadditionen und nichts anderem, und zwei davon sind unabhängig und können per Pipeline übertragen werden.

Im Jahr 1999 benötigten FPUs bereits eine Pipeline-Additions- / Subtraktionsschaltung und eine Pipeline-Multiplikationsschaltung, da SSE sonst nicht sehr "Streaming" wäre. 1999 wurde nur eine von jeder Schaltung benötigt, um diese innere Schleife vollständig zu implementieren, ohne viel Hardware nur an der Quadratwurzel zu verschwenden.

Heute haben wir natürlich Multiply-Add zusammengeführt, das dem Programmierer ausgesetzt ist. Wiederum besteht die innere Schleife aus drei FMAs mit Pipeline, die (wieder) im Allgemeinen nützlich sind, selbst wenn Sie keine Quadratwurzeln berechnen.

Pseudonym
quelle
1
Verwandte: Wie funktioniert sqrt () von GCC nach dem Kompilieren? Welche Wurzelmethode wird verwendet? Newton-Raphson? hat einige Links zu Hardware-Div / SQL-Ausführungseinheiten-Designs. Schnelle vektorisierte rsqrt und wechselseitig mit SSE / AVX je nach Präzision - eine Newton-Iteration in Software mit oder ohne FMA zur Verwendung _mm256_rsqrt_psmit Haswell-Perf-Analyse. Normalerweise nur eine gute Idee, wenn Sie keine andere Arbeit in der Schleife haben und den Teilerdurchsatz stark beeinträchtigen würden. HW sqrt ist Single Uop, ist also in Ordnung mit anderen Arbeiten gemischt.
Peter Cordes
-2

Dies ist schneller, da diese Anweisungen Rundungsmodi ignorieren und keine Gleitkomma-Ausnahmen oder dernormalisierten Zahlen verarbeiten. Aus diesen Gründen ist es viel einfacher, andere fp-Anweisungen außerhalb der Reihenfolge zu leiten, zu spekulieren und auszuführen.

Witek
quelle
Offensichtlich falsch. FMA hängt vom aktuellen Rundungsmodus ab, hat jedoch bei Haswell und höher einen Durchsatz von zwei pro Takt. Mit zwei FMA-Einheiten mit vollständiger Pipeline kann Haswell bis zu 10 FMAs gleichzeitig im Flug haben. Die richtige Antwort ist rsqrtdie viel geringere Genauigkeit, was bedeutet, dass nach einer Tabellensuche viel weniger Arbeit (oder gar keine?) Zu erledigen ist, um eine erste Vermutung zu erhalten.
Peter Cordes