Mengapa SSE scalar sqrt (x) lebih lambat dari rsqrt (x) * x?

106

Saya telah membuat profil beberapa matematika inti kami pada Intel Core Duo, dan ketika melihat berbagai pendekatan untuk akar kuadrat, saya telah melihat sesuatu yang aneh: menggunakan operasi skalar SSE, lebih cepat mengambil akar kuadrat timbal balik dan mengalikannya untuk mendapatkan sqrt, daripada menggunakan opcode sqrt asli!

Saya mengujinya dengan loop seperti:

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() );
}

Saya sudah mencoba ini dengan beberapa badan berbeda untuk TestSqrtFunction, dan saya punya beberapa pengaturan waktu yang benar-benar menggaruk kepala saya. Yang terburuk dari semuanya sejauh ini adalah menggunakan fungsi sqrt () asli dan membiarkan kompiler "pintar" "mengoptimalkan". Pada 24ns / float, menggunakan x87 FPU, ini sangat buruk:

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

Hal berikutnya yang saya coba adalah menggunakan intrinsik untuk memaksa kompiler menggunakan opcode sqrt skalar SSE:

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
}

Ini lebih baik, pada 11.9ns / float. Saya juga mencoba teknik perkiraan Newton-Raphson yang aneh dari Carmack , yang berjalan bahkan lebih baik daripada perangkat kerasnya, pada 4.3ns / float, meskipun dengan kesalahan 1 dalam 2 10 (yang terlalu berlebihan untuk tujuan saya).

Doozy adalah ketika saya mencoba operasi SSE untuk akar kuadrat timbal balik , dan kemudian menggunakan perkalian untuk mendapatkan akar kuadrat (x * 1 / √x = √x). Meskipun ini membutuhkan dua operasi yang bergantung, ini adalah solusi tercepat sejauh ini, pada 1,24ns / float dan akurat hingga 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
}

Pertanyaan saya pada dasarnya adalah apa yang memberi ? Mengapa opcode akar kuadrat bawaan ke perangkat keras SSE lebih lambat daripada menyintesisnya dari dua operasi matematika lainnya?

Saya yakin ini benar-benar biaya operasi itu sendiri, karena saya telah memverifikasi:

  • Semua data cocok dengan cache, dan aksesnya berurutan
  • fungsinya sebaris
  • membuka gulungan loop tidak ada bedanya
  • bendera kompiler disetel ke optimasi penuh (dan perakitannya bagus, saya centang)

( edit : stephentyrone dengan benar menunjukkan bahwa operasi pada string angka yang panjang harus menggunakan operasi yang dikemas SIMD vektor, seperti rsqrtps- tetapi struktur data array di sini hanya untuk tujuan pengujian: apa yang sebenarnya saya coba ukur adalah kinerja skalar untuk digunakan dalam kode yang tidak dapat divektorisasi.)

Crashworks
sumber
13
x / akar persegi (x) = akar persegi (x). Atau, dengan kata lain: x ^ 1 * x ^ (- 1/2) = x ^ (1 - 1/2) = x ^ (1/2) = sqrt (x)
Crashworks
6
tentu saja inline float SSESqrt( float restrict fIn ) { float fOut; _mm_store_ss( &fOut, _mm_sqrt_ss( _mm_load_ss( &fIn ) ) ); return fOut; }. Tetapi ini adalah ide yang buruk karena dapat dengan mudah menyebabkan kemacetan pemuatan-hit-store jika CPU menulis float ke stack dan kemudian segera membacanya kembali - beralih dari register vektor ke register float untuk nilai kembalian khususnya adalah berita buruk. Selain itu, opcode mesin yang mendasari yang diwakili oleh SSE intrinsik take address operand.
Crashworks
4
Seberapa penting LHS bergantung pada gen tertentu dan langkah dari x86 tertentu: pengalaman saya adalah bahwa pada apa pun hingga i7, memindahkan data antar set register (misalnya FPU ke SSE ke eax) sangat buruk, sementara perjalanan bolak-balik antara xmm0 dan stack dan tidak kembali, karena penerusan toko Intel. Anda bisa mengatur waktunya sendiri untuk memastikannya. Umumnya cara termudah untuk melihat potensi LHS adalah dengan melihat rakitan yang dipancarkan dan melihat di mana data disulap di antara set register; kompiler Anda mungkin melakukan hal yang cerdas, atau mungkin juga tidak. Untuk menormalkan vektor, saya menulis hasil saya di sini: bit.ly/9W5zoU
Crashworks
2
Untuk PowerPC, ya: IBM memiliki simulator CPU yang dapat memprediksi LHS dan banyak gelembung pipa lainnya melalui analisis statis. Beberapa PPC juga memiliki penghitung perangkat keras untuk LHS yang dapat Anda jajak pendapat. Lebih sulit untuk x86; alat pembuatan profil yang baik lebih langka (VTune agak rusak akhir-akhir ini) dan saluran pipa yang disusun ulang kurang deterministik. Anda dapat mencoba mengukurnya secara empiris dengan mengukur instruksi per siklus, yang dapat dilakukan secara tepat dengan penghitung kinerja perangkat keras. Register "instruksi pensiun" dan "siklus total" dapat dibaca dengan misalnya PAPI atau PerfSuite ( bit.ly/an6cMt ).
Crashworks
2
Anda juga dapat dengan mudah menulis beberapa permutasi pada suatu fungsi dan mengatur waktunya untuk melihat apakah ada yang menderita terutama dari kios. Intel tidak mempublikasikan banyak detail tentang cara kerja pipeline mereka (bahwa mereka LHS sama sekali adalah semacam rahasia kotor), jadi banyak yang saya pelajari adalah dengan melihat skenario yang menyebabkan macet di arch lain (misalnya PPC ), lalu membuat eksperimen terkontrol untuk melihat apakah x86 juga memilikinya.
Crashworks

Jawaban:

216

sqrtssmemberikan hasil yang benar. rsqrtssmemberikan perkiraan ke timbal balik, akurat sekitar 11 bit.

sqrtssmenghasilkan hasil yang jauh lebih akurat, karena diperlukan akurasi. rsqrtssada untuk kasus-kasus ketika perkiraan sudah cukup, tetapi kecepatan diperlukan. Jika Anda membaca dokumentasi Intel, Anda juga akan menemukan urutan instruksi (pendekatan akar kuadrat timbal balik diikuti oleh langkah Newton-Raphson tunggal) yang memberikan presisi hampir penuh (~ 23 bit akurasi, jika saya ingat dengan benar), dan masih agak lebih cepat dari sqrtss.

sunting: Jika kecepatan sangat penting, dan Anda benar-benar memanggil ini dalam satu lingkaran untuk banyak nilai, Anda harus menggunakan versi vektor dari instruksi ini, rsqrtpsatau sqrtps, keduanya memproses empat pelampung per instruksi.

Stephen Canon
sumber
3
Langkah n / r memberi Anda akurasi 22-bit (menggandakannya); 23-bit akan menjadi akurasi penuh.
Jasper Bekkers
7
@Jasper Bekkers: Tidak, tidak akan. Pertama, float memiliki presisi 24 bit. Kedua, sqrtssadalah benar bulat , yang membutuhkan ~ 50 bit sebelum pembulatan, dan tidak dapat dicapai dengan menggunakan sederhana N / R iterasi dalam presisi tunggal.
Stephen Canon
1
Ini pasti alasannya. Untuk memperluas hasil ini: Proyek Embree Intel ( software.intel.com/en-us/articles/… ), menggunakan vektorisasi untuk matematikanya. Anda dapat mengunduh sumbernya di tautan itu dan melihat bagaimana mereka melakukan Vektor 3/4 D mereka. Normalisasi vektornya menggunakan rsqrt diikuti dengan iterasi newton-raphson, yang kemudian sangat akurat dan masih lebih cepat dari 1 / ssqrt!
Brandon Pelfrey
7
Peringatan kecil: x rsqrt (x) menghasilkan NaN jika x adalah nol atau tak terbatas. 0 * rsqrt (0) = 0 * INF = NaN. INF rsqrt (INF) = INF * 0 = NaN. Untuk alasan ini, CUDA pada GPU NVIDIA menghitung perkiraan akar kuadrat presisi tunggal sebagai timbal (rsqrt (x)), dengan perangkat keras memberikan perkiraan cepat ke akar kuadrat timbal balik dan timbal balik. Jelas, pemeriksaan eksplisit yang menangani dua kasus khusus juga dimungkinkan (tetapi akan lebih lambat di GPU).
njuffa
@BrandonPelfrey Di file manakah Anda menemukan langkah Newton Rhapson?
fredoverflow
7

Ini juga berlaku untuk divisi. MULSS (a, RCPSS (b)) jauh lebih cepat dari DIVSS (a, b). Nyatanya masih lebih cepat bahkan saat Anda meningkatkan presisi dengan iterasi Newton-Raphson.

Intel dan AMD sama-sama merekomendasikan teknik ini dalam manual pengoptimalan mereka. Dalam aplikasi yang tidak memerlukan kepatuhan IEEE-754, satu-satunya alasan untuk menggunakan div / sqrt adalah keterbacaan kode.

Bertengkar
sumber
1
Broadwell dan yang lebih baru memiliki kinerja pembagian FP yang lebih baik, sehingga penyusun seperti clang memilih untuk tidak menggunakan timbal balik + Newton untuk skalar pada CPU terbaru, karena biasanya tidak lebih cepat. Dalam kebanyakan loop, divbukan satu-satunya operasi, jadi total throughput uop sering kali menjadi bottleneck bahkan ketika ada divpsatau divss. Lihat pembagian floating point vs perkalian floating point , di mana jawaban saya memiliki bagian tentang mengapa rcppsthroughput tidak lagi menang. (Atau kemenangan latensi), dan angka pada pembagian throughput / latensi.
Peter Cordes
Jika persyaratan akurasi Anda sangat rendah sehingga Anda dapat melewati iterasi Newton, ya a * rcpss(b)bisa lebih cepat, tetapi masih lebih uops daripada a/b!
Peter Cordes
5

Alih-alih memberikan jawaban, itu sebenarnya mungkin salah (saya juga tidak akan memeriksa atau berdebat tentang cache dan hal lainnya, katakanlah keduanya identik) Saya akan mencoba mengarahkan Anda ke sumber yang dapat menjawab pertanyaan Anda.
Perbedaannya mungkin terletak pada bagaimana akar dan akar dihitung. Anda dapat membaca lebih lanjut di sini http://www.intel.com/products/processor/manuals/ . Saya sarankan untuk memulai dari membaca tentang fungsi prosesor yang Anda gunakan, ada beberapa info, terutama tentang rsqrt (cpu menggunakan tabel pencarian internal dengan perkiraan besar, yang membuatnya lebih mudah untuk mendapatkan hasilnya). Tampaknya, rsqrt jauh lebih cepat daripada sqrt, sehingga 1 operasi mul tambahan (yang tidak terlalu mahal) mungkin tidak mengubah situasi di sini.

Sunting: Beberapa fakta yang mungkin layak untuk disebutkan:
1. Setelah saya melakukan beberapa optimalisasi mikro untuk perpustakaan grafik saya dan saya telah menggunakan rsqrt untuk menghitung panjang vektor. (alih-alih sqrt, saya telah mengalikan jumlah kuadrat saya dengan rsqrt, yang persis seperti yang telah Anda lakukan dalam pengujian), dan hasilnya lebih baik.
2. Menghitung rsqrt menggunakan tabel lookup sederhana mungkin lebih mudah, seperti untuk rsqrt, ketika x pergi ke tak terhingga, 1 / sqrt (x) pergi ke 0, jadi untuk x kecil nilai fungsinya tidak berubah (banyak), sedangkan untuk sqrt - itu menuju tak terbatas, jadi itu kasus sederhana itu;).

Juga, klarifikasi: Saya tidak yakin di mana saya menemukannya di buku yang saya tautkan, tetapi saya cukup yakin saya telah membaca bahwa rsqrt menggunakan beberapa tabel pencarian, dan itu harus digunakan hanya, ketika hasilnya tidak perlu persisnya, meskipun - saya mungkin juga salah, seperti beberapa waktu lalu :).

Marcin Deptuła
sumber
4

Newton-Raphson konvergen ke nol f(x)menggunakan kenaikan sama dengan di -f/f' mana f'turunannya.

Karena x=sqrt(y), Anda dapat mencoba memecahkan f(x) = 0untuk xmenggunakan f(x) = x^2 - y;

Kemudian kenaikannya adalah: dx = -f/f' = 1/2 (x - y/x) = 1/2 (x^2 - y) / x yang memiliki pembagian lambat di dalamnya.

Anda dapat mencoba fungsi lain (seperti f(x) = 1/y - 1/x^2) tetapi akan sama rumitnya.

Mari kita simak 1/sqrt(y)sekarang. Anda dapat mencoba f(x) = x^2 - 1/y, tetapi akan sama rumitnya: dx = 2xy / (y*x^2 - 1)misalnya. Satu pilihan alternatif yang tidak jelas f(x)adalah:f(x) = y - 1/x^2

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

Ah! Ini bukan ekspresi yang sepele, tetapi Anda hanya memiliki perkalian di dalamnya, tidak ada pembagian. => Lebih cepat!

Dan: langkah pembaruan penuh new_x = x + dxkemudian berbunyi:

x *= 3/2 - y/2 * x * x yang juga mudah.

skal
sumber
2

Ada sejumlah jawaban lain untuk ini dari beberapa tahun yang lalu. Inilah yang konsensusnya benar:

  • Instruksi rsqrt * menghitung pendekatan ke akar kuadrat timbal balik, hingga sekitar 11-12 bit.
  • Ini diimplementasikan dengan tabel pencarian (yaitu ROM) yang diindeks oleh mantissa. (Faktanya, ini adalah tabel pencarian terkompresi, mirip dengan tabel matematika lama, menggunakan penyesuaian pada bit orde rendah untuk menghemat transistor.)
  • Alasan mengapa ini tersedia adalah karena ini adalah perkiraan awal yang digunakan oleh FPU untuk algoritme akar kuadrat "sebenarnya".
  • Ada juga instruksi timbal balik perkiraan, rcp. Kedua instruksi ini adalah petunjuk tentang bagaimana FPU mengimplementasikan akar dan pembagian kuadrat.

Inilah yang salah konsensus:

  • FPU era SSE tidak menggunakan Newton-Raphson untuk menghitung akar kuadrat. Ini adalah metode yang bagus dalam perangkat lunak, tetapi akan menjadi kesalahan untuk menerapkannya seperti itu di perangkat keras.

Algoritme NR untuk menghitung akar kuadrat timbal balik memiliki langkah pembaruan ini, seperti yang telah dicatat orang lain:

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

Itu banyak perkalian bergantung data dan satu pengurangan.

Berikut ini adalah algoritme yang sebenarnya digunakan oleh FPU modern.

Diberikan b[0] = n, misalkan kita dapat menemukan serangkaian angka Y[i]sedemikian rupa sehingga b[n] = b[0] * Y[0]^2 * Y[1]^2 * ... * Y[n]^2mendekati 1. Kemudian pertimbangkan:

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

Jelas x[n]pendekatan sqrt(n)dan y[n]pendekatan 1/sqrt(n).

Kita bisa menggunakan langkah pembaruan Newton-Raphson untuk akar kuadrat timbal balik untuk mendapatkan barang Y[i]:

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

Kemudian:

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

dan:

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

Pengamatan kunci berikutnya adalah itu b[i] = x[i-1] * y[i-1]. Begitu:

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

Kemudian:

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]))

Artinya, diberi inisial x dan y, kita dapat menggunakan langkah pembaruan berikut:

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

Atau, lebih bagus lagi, kita bisa atur h = 0.5 * y. Ini adalah inisialisasi:

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

Dan ini adalah langkah pembaruannya:

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

Ini adalah algoritme Goldschmidt, dan memiliki keuntungan besar jika Anda menerapkannya di perangkat keras: "loop dalam" terdiri dari tiga penambahan-perkalian dan tidak ada yang lain, dan dua di antaranya independen dan dapat disambungkan.

Pada tahun 1999, FPU sudah membutuhkan sirkuit add / substract pipelined dan sirkuit multiply pipelined, jika tidak SSE tidak akan terlalu "streaming". Hanya satu dari setiap sirkuit yang dibutuhkan pada tahun 1999 untuk mengimplementasikan loop dalam ini dengan cara yang sepenuhnya terhubung tanpa membuang banyak perangkat keras hanya pada akar kuadrat.

Hari ini, tentu saja, kami telah menyatukan multiply-add yang diekspos ke programmer. Sekali lagi, loop dalam adalah tiga FMA pipelined, yang (sekali lagi) umumnya berguna bahkan jika Anda tidak menghitung akar kuadrat.

Nama samaran
sumber
1
Terkait: Bagaimana sqrt () dari GCC bekerja setelah dikompilasi? Metode root mana yang digunakan? Newton-Raphson? memiliki beberapa tautan ke desain unit eksekusi div / sqrt perangkat keras. RSqrt vektorisasi cepat dan timbal balik dengan SSE / AVX bergantung pada presisi - satu iterasi Newton dalam perangkat lunak, dengan atau tanpa FMA, untuk digunakan dengan _mm256_rsqrt_ps, dengan analisis kinerja Haswell. Biasanya hanya merupakan ide yang bagus jika Anda tidak memiliki pekerjaan lain dalam loop dan akan menghambat throughput pembagi. HW sqrt adalah single uop jadi tidak apa-apa dicampur dengan pekerjaan lain.
Peter Cordes
-2

Ini lebih cepat karena instruksi ini mengabaikan mode pembulatan, dan tidak menangani pengecualian titik floatin atau angka yang dinormalisasi. Untuk alasan ini, jauh lebih mudah untuk melakukan pipeline, berspekulasi, dan mengeksekusi instruksi fp lainnya Rusak.

Witek
sumber
Jelas salah. FMA bergantung pada mode pembulatan saat ini, tetapi memiliki throughput dua per jam di Haswell dan yang lebih baru. Dengan dua unit FMA dengan pipeline penuh, Haswell dapat memiliki hingga 10 FMA dalam penerbangan sekaligus. Jawaban yang benar adalah rsqrt's banyak akurasi yang lebih rendah, yang berarti lebih sedikit pekerjaan yang harus dilakukan (atau tidak sama sekali?) Setelah meja-lookup untuk mendapatkan menebak mulai.
Peter Cordes