文档内容简介：
SSE（為 Streaming SIMD Extensions 的縮寫）是由 Intel 公司，在 1999 年推出 Pentium III 處理器時，同時推出的新指令集。如同其名稱所表示的，SSE 是一種 SIMD 指令集。所謂的 SIMD 是指 single instruction, multiple data，也就是一個指令同時對多個資料進行相同的動作。較早的 MMX 和 AMD 的 3DNow! 也都是 SIMD 指令集。因此，SSE 本質上是非常類似一個向量處理器的。SSE 指令包括了四個主要的部份：單精確度浮點數運算指令、整數運算指令（此為 MMX 之延伸，並和 MMX 使用同樣的暫存器）、Cache 控制指令、和狀態控制指令。這裡主要是介紹浮點數運算指令和 Cache 控制指令。
目录：
[Part 1]
[Part 2]
[Part 3]
[Part 4]
[Part 5]
[Part 6]
[Part 7]
[Part 8]
文档内容：
[Part 1]
SSE 新增了八個新的 128 位元暫存器，xmm0 ~ xmm7。這些 128 位元的暫存器，可以用來存放四個 32 位元的單精確度浮點數。SSE 的浮點數運算指令就是使用這些暫存器。和之前的 MMX 或 3DNow! 不同，這些暫存器並不是原來己有的暫存器（MMX 和 3DNow! 均是使用 x87 浮點數暫存器），所以不需要像 MMX 或 3DNow! 一樣，要使用 x87 指令之前，需要利用一個 EMMS 指令來清除暫存器的狀態。因此，不像 MMX 或 3DNow! 指令，SSE 的浮點數運算指令，可以很自由地和 x87 指令，或是 MMX 指令共用。但是，這樣做的主要缺點是，因為多工作業系統在進行 context switch 時，需要儲存所有暫存器的內容。而這些多出來的新暫存器，也是需要儲存的。因此，既存的作業系統需要修改，在 context switch 時，儲存這八個新暫存器的內容，才能正確支援 SSE 浮點運算指令。下圖是 SSE 新增的暫存器的示意圖：
目前 Intel 新推出的 x86 CPU 均已支援 SSE 指令集，包括 Pentium III、Pentium 4、Celeron 533A 及之後的處理器。AMD 方面，在其 Athlon 處理器上，新增的 Enhanced 3DNow! 指令集中，即包含了 SSE 中的整數運算指令和 Cache 控制指令，但仍不包含浮點數指令和狀態控制指令。但在最新的 Athlon XP 處理器上，已將 SSE 中的浮點運算指令和狀態控制指令加入。因此 Athlon XP 已經是完全支援 SSE 指令集了。
[Part 2]
SSE 的浮點數運算指令，大致上可以分成兩種：packed 和 scalar。Packed 指令是一次對 XMM 暫存器中的四個浮點數（即 DATA0 ~ DATA3）均進行計算，而 scalar 則只對 XMM 暫存器中的 DATA0 進行計算。如下圖所示：
SSE 指令和一般的 x86 指令很類似，基本上包括兩種定址方式：reg-reg 和 reg-mem。下面是兩個例子：
addps xmm0, xmm1 ; reg-reg-
addps xmm0, [ebx] ; reg-mem
指令的運算結果會覆蓋到第一個參數中。例如，以上面的第一個例子來說，xmm0 暫存器會存放最後計算的結果。
另外，絕大部份需要存取記憶體的 SSE 指令，都要求位址是 16 的倍數（也就是對齊在 16 bytes 的邊上）。如果不是的話，就會導致 exception。這是非常重要的。因為，一般的 32 位元浮點數只會對齊在 4 bytes 或 8 bytes 的邊上（根據 compiler 的設定而不同）。另外，若是處理陣列中的數字，也需要特別注意這個問題。
[Part 3]
目前要在 C 或 C++ 程式中，利用 SSE 指令，有很多方法。最一般的方式，是利用內嵌式組合語言。現在包括 Intel C++ Compiler 5.0 和 Microsoft Visual C++ 6.0 Processor Pack 的內嵌式組合語言都支援 SSE 指令集。另外，它們也都支援 SSE 指令集的 intrinsics。Intrinsics 是一種外表類似一般的函數，但是實際上會被 compiler 直接編譯成組合語言的東西。以下是一些例子：
內嵌式組合語言：
__asm addps xmm0, xmm1
__asm movaps [ebx], xmm0
...
__m128 data;
...
__asm
{
lea ebx, data
addps xmm0, xmm1
movaps [ebx], xmm0
}
使用 intrinsics：
__m128 data1, data2;
...
__m128 out = _mm_add_ps(data1, data2);
...
使用 intrinsics 可以增加程式的可讀性，也比較容易使用。不過，在某些情形下，compiler 可能沒辦法產生最好的程式碼，而且，其產生的程式碼的效率，也會隨著 compiler 的不同而有改變。但是，對於大部份的應用來說，使用 intrinsics 的好處通常是很明顯的。因此，在這裡我們也以 intrinsics 為主，而不討論內嵌式組合語言。
SSE 的 intrinsics 的名稱看起來雖然怪異，但是它是有一定的規則的。基本上一個 SSE 的 intrinsics 的型式為：
_mm_<opcode>_<suffix>
其中，<opcode> 是指令的類別，像是 add、sub 等等。而 <suffix> 則是資料的種類。在 SSE 浮點運算指令中，只有兩個種類：ps 和 ss。ps 是指 Packed Single-precision，也就是這個指令對暫存器中的四個單精度浮點數進行運算。而 ss 則是指 Scalar Single-precision，也就是這個指令只對暫存器中的 DATA0 進行運算。所以，像上面的 _mm_add_ps 指令，就是把兩個四維向量相加的指令。
SSE 也需要一個新的資料型態，也就是上面所看到的 __m128。__m128 是一個 16 bytes（128 bits）的資料型態，對應 SSE 的 128 位元暫存器。幾乎所有的 SSE 浮點運算的 intrinsics 都是使用這個資料型態。比如說，_mm_add_ps 這個 intrinsic 的函數宣告為：
__m128 _mm_add_ps(__m128 a, __m128 b);
可以看到，它的參數和傳回值的型態都是 __m128。基本上，這個 intrinsic 的動作就是把兩個參數相加，並把結果以傳回值的方式傳回。
另外，要使用 SSE 的 intrinsics 之前，要記得先包含 xmmintrin.h 這個 header 檔。如下：
#include <xmmintrin.h>
這樣才能使用這些 instrinsics。另外，這個檔案也包含了一些方便的巨集定義，用來配合一些 SSE 的 intrinsics 使用。
[Part 4]
現在我們就先以一個很簡單的例子，來說明 SSE 浮點運算指令的使用。下面是一個簡單的程式片斷：
- float input1[4] = { 1.2f, 3.5f, 1.7f, 2.8f };
-
- float input2[4] = { -0.7f, 2.6f, 3.3f, -0.8f };
- float output[4];
-
- for(int i = 0; i < 4; i++) {
-
- output[i] = input1[i] + input2[i];
-
-
- }
-
這個程式片斷的動作非常簡單，只是把 input1 和 input2 中的四組數字兩兩相加，並把結果寫到 output 中。在執行完後，output 的內容應該是 { 0.5f, 6.1f, 5.0f, 2.0f }。
上面的程式很適合用 SSE 浮點運算指令來做。它所進行的四個加法運算，剛好可以用一個 SSE 浮點運算指令，也就是 addps 來完成。以下是修改後的程式：
- float input1[4] = { 1.2f, 3.5f, 1.7f, 2.8f };
-
- float input2[4] = { -0.7f, 2.6f, 3.3f, -0.8f };
- float output[4];
-
- __m128 a = _mm_load_ps(input1);
-
- __m128 b = _mm_load_ps(input2);
-
- __m128 t = _mm_add_ps(a, b);
-
- _mm_store_ps(output, t);
-
上面的程式中，先用兩個 _mm_load_ps 把浮點數資料載入 __m128 的變數中。要注意的是，這裡是假設這兩個浮點數陣列都是對齊在 16 bytes 的邊上。如果不是的話，就不能用 _mm_load_ps 這個 intrinsic 來載入，而要改用 _mm_loadu_ps 這個 intrinsic。它是專門用來處理沒有對齊在 16 bytes 邊上的資料的。但是，它的速度會比較慢。
另外，因為 x86 的 little endian 特性，位址較低的 byte 會放在暫存器的右邊。也就是說，若以上面的 input1 為例，在載入到 XMM 暫存器後，暫存器中的 DATA0 會是 1.2，而 DATA1 是 3.5，DATA2 是 1.7，DATA3 是 2.8。如果需要以相反的順序載入的話，可以用 _mm_loadr_ps 這個 intrinsic。當然，在這個例子中，順序並不影響結果，所以不需要利用這個 intrinsic。
一般來說，宣告一個 float 的陣列，並不會對齊在 16 bytes 的邊上。如果希望它會對齊在 16 bytes 的邊上，以便利用 SSE 指令的話，Visual C++ 6.0 Processor Pack 和 Intel C++ compiler 都可以指定對齊方式。指定的方式如下：
-
- __declspec(align(16)) float input1[4];
在上例中，input1 會被對齊在 16 bytes 的邊上。這樣就可以直接用較快的 _mm_load_ps 來載入資料了。因為 SSE 浮點數指令常常需要資料對齊在 16 bytes 的邊上，所以在 xmmintrin.h 也定義了一個巨集 _MM_ALIGN16，是同樣的意義。因此，上面的程式也可以寫成：
-
- _MM_ALIGN16 float input1[4];
利用 SSE 浮點數指令會帶來多大的差別呢？對 1,024 個浮點數進行測試的結果，結果如下（使用 Intel C++ Compiler 5.0 for Windows 編譯，在 Intel Pentium III 1.0B Ghz 上執行）：
可以看到利用 SSE 浮點運算可以得到大約兩倍於 x87 浮點運算的效率。事實上，因為 Pentium III 架構上的因素，所以，如果同時使用加法和乘法運算，最多可以得到四倍的效率。在後面會有一個例子。
[Part 5]
看了簡單的例子之後，我們現在列出一些 SSE 浮點運算指令所支援的一些基本的運算：
| 指令 | Intrinsic | 功能 |
addps/addss | _mm_add_ps/_mm_add_ss | 加法 |
subps/subss | _mm_sub_ps/_mm_sub_ss | 減法 |
mulps/mulss | _mm_mul_ps/_mm_mul_ss | 乘法 |
divps/divss | _mm_div_ps/_mm_div_ss | 除法 |
sqrtps/sqrtss | _mm_sqrt_ps/_mm_sqrt_ss | 平方根 |
maxps/maxss | _mm_max_ps/_mm_max_ss | 逐項取最大值 |
minps/minss | _mm_min_ps/_mm_min_ss | 逐項取最小值 |
這些運算基本上都符合 IEEE 754 中的規範，根據其計算結果，設定適當的旗標（divide-by-zero、invalid 等等），或是產生 exception。這可以由一個 MXCSR 暫存器來設定。MXCSR 暫存器是一個 32 位元的旗標暫存器，可以設定是否要產生各種 exception，並會記錄上次的計算中，發生了哪些情況。下面是 MXCSR 暫存器的說明圖：
以下是 MXCSR 各欄位的說明：
| 欄位 | 名稱 | Intrinsic 巨集 | 說明 |
| IE | Invalid Operation Flag | _MM_EXCEPT_INVALID | 運算中發生 invalid operation |
| DE | Denormal Flag | _MM_EXCEPT_DENORMA | 試圖載入 denormal value 或對 denormal value 進行運算 |
| ZE | Divide-by-Zero Flag | _MM_EXCEPT_DIV_ZERO | 運算中發生 divide by zero |
| OE | Overflow Flag | _MM_EXCEPT_OVERFLOW | 運算中發生 overflow |
| UE | Underflow Flag | _MM_EXCEPT_UNDERFLOW | 運算中發生 underflow |
| PE | Precision Flag | _MM_EXCEPT_INEXACT | 運算結果不精確（inexact） |
| DAZ | Denormals Are Zeros | N/A | 強制將 denormals 視為 0 |
| IM | Invalid Operation Mask | _MM_MASK_INVALID | 關閉 invalid operation exception |
| DM | Denormal Operation Mask | _MM_MASK_DENORM | 關閉 denormal operation exception |
| ZM | Divide-by-Zero Mask | _MM_MASK_DIV_ZERO | 關閉 divide by zero exception |
| OM | Overflow Mask | _MM_MASK_OVERFLOW | 關閉 overflow exception |
| UM | Underflow Mask | _MM_MASK_UNDERFLOW | 關閉 underflow exception |
| PM | Precision Mask | _MM_MASK_INEXACT | 關閉 inexact exception |
| RC | Rounding Control | N/A | 設定捨去方式 |
| FZ | Flush to Zero | N/A | 打開 Flush-to-zero 模式 |
上表中有幾個地方是需要特別注意的：
-
- DAZ 若設為 1，則會開啟 Denormals Are Zeros 模式。當此模式開啟時，若來源參數中有任何 denormals，都會被視為 0，而且不會設定 DE 或產生 denormal operation exception。這個模式和 IEEE 754 規範不相容，但是通常會得到更好的效率。另外 DAZ 目前只在 Intel Pentium 4 和 Intel Xeon 處理器上有支援。有關判斷 DAZ 是否支援的方式，請參考 IA-32 Intel Architecture Software Developer's Manual, Volume 1: Basic Architecture, Section 11.6.3。
- RC 可以設定成四種方式：00 為 Round to nearest (even)，即將運算結果捨去至最接近數值。如果兩個數值同樣接近，則捨去至偶數。01 為 Round down，即往負無限大方向捨去。10 為 Round up，即往正無限大方向捨去。11 則是 Round toward zero (truncate)，即往零的方向捨去。
- FZ 若設為 1，則會開啟 Flush-to-zero 模式。當此模式開啟，且 underflow exception 被關閉時，若在運算時發生 underflow，則結果為 0（符號同真正結果），且會設定 PE 和 UE。若 underflow exception 未關閉，則此模式無任何影響。注意此模式和 IEEE 754 規範不相容（規範要求 underflow 需產生 denormal 結果），但通常會得到更好的效率。
設定和讀取 MXCSR 同樣有 intrinsic 可以用，即 _mm_getcsr 和 _mm_setcsr。不過，因為 MXCSR 基本上是由 bit field 組成，因此在 xmmintrin.h 裡面也定義了一些方便的巨集。這些巨集包括：
| 巨集 | 功用 |
_MM_SET_EXCEPTION_STATE | 設定 exception 狀態（IE ~ PE） |
_MM_GET_EXCEPTION_STATE | 取得 exception 狀態（IE ~ PE） |
_MM_SET_EXCEPTION_MASK | 設定 exception mask（IM ~ PM） |
_MM_GET_EXCEPTION_MASK | 取得 exception mask（IM ~ PM） |
_MM_SET_ROUNDING_MODE | 設定捨去方式 |
_MM_GET_ROUNDING_MODE | 取得目前的捨去方式 |
_MM_SET_FLUSH_ZERO_MODE | 設定 Flush-to-zero 模式（_MM_FLUSH_ZERO_ON 和 _MM_FLUSH_ZERO_OFF） |
_MM_GET_FLUSH_ZERO_MODE | 取得 Flush-to-zero 模式設定 |
捨去方式的設定包括：
_MM_ROUND_NEAREST | Round to nearest (even) |
_MM_ROUND_DOWN | Round down |
_MM_ROUND_UP | Round up |
_MM_ROUND_TOWARD_ZERO | Round toward zero (truncate) |
在開機時，所有的 exception mask 都是設定為 1，也就是所有的 exception 都是關閉的。如果有需要的話，可以將其開啟。下面的程式片斷是一個例子：
- _MM_ALIGN16 float test1[4] = { 0, 0, 0, 1 };
-
- _MM_ALIGN16 float test2[4] = { 1, 2, 3, 0 };
- _MM_ALIGN16 float out[4];
-
- _MM_SET_EXCEPTION_MASK(0); // 打開所有的 exception
-
- __try {
-
- __m128 a = _mm_load_ps(test1);
-
- __m128 b = _mm_load_ps(test2);
-
- a = _mm_div_ps(a, b);
-
- _mm_store_ps(out, a);
-
-
- }
-
- __except(EXCEPTION_EXECUTE_HANDLER) { // exception handler
-
- if(_mm_getcsr() & _MM_EXCEPT_DIV_ZERO) {
-
- cout << "Divide by zero" << endl;
-
-
- }
-
- return;
- if(_mm_getcsr() & _MM_EXCEPT_DIV_ZERO) {
-
-
- }
-
- cout << out[3] << endl;
- __except(EXCEPTION_EXECUTE_HANDLER) { // exception handler
上面的程式在執行時，會產生 divide by zero 的 exception，並進入 exception handler。在 exception handler 中，因為 ZE 會被設定，因此會顯示出 "Divide by zero" 的錯誤訊息。如果把 _MM_SET_EXCEPTION_MASK(0); 這一行去掉，則不會產生 exception，而會直接輸出正無限大（1.#INF）的結果。
[Page 6]
現在我們來看看一個實際的例子，是一個 3D 繪圖運算中很常需要進行的動作，也就是三維向量的內積。
兩個三維向量 <x1, y1, z1> 和 <x2, y2, z2> 的定義為：
x1x2 + y1y2 + z1z2
但是 SSE 浮點運算指令是以四維向量為單位的，而且也沒有計算內積的指令。要直接對一個三維向量計算出內積，會需要相當多的資料重排動作，反而會浪費很多時間。
而且，如果三維向量在記憶體中是以 x1, y1, z1, x2, y2, z2, ... 的方式排列的話，還會有對齊的問題。因為一個三維向量只需要 12 bytes，所以無法讓每個向量都對齊在 16 bytes 的邊上。這也會降低效率。
一個簡單的方法是改變資料的排列方式。如果我們把資料改成以下面的方式排列：
| x1, x2, x3, x4, ... |
| y1, y2, y3, y4, ... |
| z1, z2, z3, z4, ... |
那我們就同時解決了兩個問題。首先，因為一次讀入四個 x 值（及 y 值和 z 值），所以不會有對齊的問題。另外，運算也變得容易。下面的程式片斷一次可以算出四個向量內積：
-
- __m128 x1 = _mm_load_ps(vec1_x);
- __m128 y1 = _mm_load_ps(vec1_y);
- __m128 z1 = _mm_load_ps(vec1_z);
- __m128 x2 = _mm_load_ps(vec2_x);
- __m128 y2 = _mm_load_ps(vec2_y);
- __m128 z2 = _mm_load_ps(vec2_z);
- __m128 t1 = _mm_mul_ps(x1, x2);
- __m128 t2 = _mm_mul_ps(y1, y2);
- t1 = _mm_add_ps(t1, t2);
- t2 = _mm_mul_ps(z1, z2);
- t1 = _mm_add_ps(t1, t2);
- _mm_store_ps(output, t1);
由於 SIMD 指令的特性，很多時候，不要把工作當成向量來考慮。而應該是把多個工作，結合到一個向量中來進行。以上面的例子來說，一個三維向量的內積，是很難有效率地用 SSE 指令來完成的。事實上，配合 SSE 的資料重排指令（後面會介紹），一個四維向量的內積，需要八個 SSE 浮點運算指令才能完成。這並不是很有效率。但是，如果我們不要以向量為單位，而改成考慮一次對四個向量進行運算的話（如同上面的例子），四個三維向量的內積，由五個指令就可以完成了。四個四維向量的內積也只需要七個指令。
以下是以 1,024 個三維向量，以上述的資料排列方式，在 Intel Pentium III 1.0B Ghz 上測試的結果（使用 Intel C++ Compiler 5.0 for Windows 編譯）：
可以看到使用 SSE 浮點運算指令，得到的效率比使用 x87 浮點數要高出甚多，幾乎達到三倍。
[Page 7]
前面的向量內積的測試，是以 1,024 個三維向量，對一個固定的三維向量進行內積得到的結果。它的資料量算是相當小的，甚至可以整個放到 L1 cache 裡面。但是，實際情形常常不是這樣。我們常常需要對一大堆資料進行運算，因此資料是無法放在 cache 裡面的，而需要從主記憶體中取得。這樣一來，速度會變得相當的慢。下面是對 102,400 個三維向量進行同樣動作的結果：
可以看到效率差了很多，這是因為資料都需要從主記憶體中讀取，因而讀取的時間成為瓶頸，而不是計算的時間。在這種情形下，SSE 浮點運算指令速度雖快，也沒有什麼幫助。而且，因為內積的計算是相當簡單的，所以這也讓 SSE 和 x87 的結果差距變得很小。
不過，現在我們會討論到 SSE 的 "streaming" 部份。SSE 除了運算的指令之外，還支援了一些 cache 控制指令。我們在這裡先介紹兩個較簡單的指令：prefetch 和 movntps。prefetch 指令實際上有四個不同的指令，包括 prefetch0、prefetch1、prefetch2、和 prefetchnta。不過，它們都是用同一個 intrinsic 表示的，也就是 _mm_prefetch。
prefetch 指令的主要目的，是提前讓 CPU 載入稍後運算所需要的資料。通常是在對目前的資料進行運算之前，告訴 CPU 載入下一筆資料。這樣就可以讓目前的運算，和載入下一筆資料的動作，可以同時進行。如果運算的複雜度夠高的話，這樣可以完全消除讀取主記憶體的 latency。不同的 prefetch 指令則是告訴 CPU 將資料載入不同層次的 cache。不過，最常用的還是 prefetchnta，這個指令會把資料載入到離 CPU 最近的 cache 中（通常是 L1 cache 或 L2 cache），適用於資料在短時間內就會用到的情形。
另外 prefetch 指令不會產生任何 exception。它本質上只是一個 hint，CPU 並不一定會真的進行載入的動作。所以，即使 prefetch 一個不合法的記憶體位址，也不會產生錯誤。這讓程式可以不用處理討厭的邊界問題。
除了 prefetch 之外，另一個指令是 movntps，它的 intrinsics 是 _mm_stream_ps。這個指令的用途，是要求 CPU 在寫入資料的時候，不要把資料寫到 cache 中，而是直接將資料寫到主記憶體中。實際上它以 write combining 的方式寫入的。為什麼要這樣做呢？這是因為，很多時候計算的結果並不是立刻需要用到的，通常是很久以後才會用到。所以，這些資料如果被放在 cache 中，完全是浪費空間。而且，更糟的是，它們可能會把 cache 中有用的資料擠掉，而使得這些資料常常需要重新從主記憶體中載入。因此，如果讓這些資料不要被放在 cache 中，就可以避免這種問題。
對上面的程式，加上適當的 prefetch，並利用 movntps 指令，可以修改成類似下面的程式片斷：
-
- __m128 x1 = _mm_load_ps(vec1_x);
- __m128 y1 = _mm_load_ps(vec1_y);
- __m128 z1 = _mm_load_ps(vec1_z);
- __m128 x2 = _mm_load_ps(vec2_x);
- __m128 y2 = _mm_load_ps(vec2_y);
- __m128 z2 = _mm_load_ps(vec2_z);
- _mm_prefetch((const char*)(vec1_x + next), _MM_HINT_NTA);
- _mm_prefetch((const char*)(vec1_y + next), _MM_HINT_NTA);
- _mm_prefetch((const char*)(vec1_z + next), _MM_HINT_NTA);
- _mm_prefetch((const char*)(vec2_x + next), _MM_HINT_NTA);
- _mm_prefetch((const char*)(vec2_y + next), _MM_HINT_NTA);
- _mm_prefetch((const char*)(vec2_z + next), _MM_HINT_NTA);
- __m128 t1 = _mm_mul_ps(x1, x2);
- __m128 t2 = _mm_mul_ps(y1, y2);
- t1 = _mm_add_ps(t1, t2);
- t2 = _mm_mul_ps(z1, z2);
- t1 = _mm_add_ps(t1, t2);
- _mm_stream_ps(output, t1);
同樣對 102,400 個三維向量進行測試，結果為：
由於 x87 無法利用 movntps 指令，所以效率較差。可以看到，加上 prefetch 和 movntps 指令後，SSE 浮點運算的效率提高了約 50%。事實上，它所使用的記憶體頻寬達到約 500MB/s，已經是理論頻寬（1,066MB/s）的一半左右了。
[Page 8]
前面介紹的計算都是蠻簡單的計算。現在我們來看一個比較複雜的例子。假設現在有一堆數字，都是單精確度的浮點數（float），而且範圍都在正負 0.5 之間。然後現在要計算它們的自然指數（exp）。當然我們知道 SSE 並沒有這個指令。事實上除了簡單的四則運算（加、減、乘、除）之外，SSE 和計算有關的指令，只有平方根而已。其它複雜的東西像是各種超越函數都是沒有的。
不過，這並不表示要算超越函數就不能用 SSE 了。這只是表示我們得自己算而已。以自然指數來說，最簡單的算法當然是用 Taylor series。事實上有更快的算法，不過我們為了舉例方便，所以就直接用 Taylor series 了。因為我們的數字範圍是在正負 0.5 之間，所以只需要取到第九項，也就是 x 的八次方項，就已經很接近單精確度浮點數所能表示的精確度了。下面是利用 Taylor series 算 exp(x) 的程式片斷：
-
- int i;
- float result = coeff[8] * x;
- for(i = 7; i >= 2; i--) {
- result += coeff[i];
- result *= x;
- }
- return (result + 1) * x + 1;
上面的程式片斷中，coeff 陣列需要事先設定成 1/n! 的值。這可以事先算好，直接放在程式裡面，所以不需要花時間。這個程式片斷是利用 Horner 法來計算 Taylor series，以減少計算量。在上例中，每個 X 需要 8 個乘法和 8 個加法，共 16 個計算。
這個程式要改成用 SSE 來寫，並不會很困難。不過，coeff 陣列需要修改成 __m128 的陣列，而且每項存放四個相同的數字（也就是 1/n!）。這樣可以改成下面的程式：
-
- int i;
- __m128 X = _mm_load_ps(data);
- __m128 result = _mm_mul_ps(coeff_sse[8], X);
- for(i = 7; i >=2; i--) {
- result = _mm_add_ps(result, coeff_sse[i]);
- result = _mm_mul_ps(result, X);
- }
- result = _mm_add_ps(result, sse_one);
- result = _mm_mul_ps(result, X);
- result = _mm_add_ps(result, sse_one);
- _mm_store_ps(out, result);
上面的程式片斷和 x87 的版本長得幾乎一樣，當然它是一次對四個浮點數做運算，所以應該會快四倍。下圖是測試結果：
這些測試結果是對 1,024 個隨機產生的數字進行的。第一個結果是直接用 exp 函式，它裡面是用 x87 的指令像 F2XM1 和 FYL2X 等指令來完成計算的。當然這個函式並沒有範圍的限制，而且它得到的結果也是最精確的。不過，它的速度非常的慢，平均每個數字需要 128 個 cycles 才能算完。第二個結果是用 x87 版的 Taylor series 來計算，也就是第一個程式片斷。它很明顯要快得多，平均每個數字只需要 66 個 cycles，幾乎是快了一倍。不過，其實這裡就產生一個問題：為什麼 16 個計算會需要 66 個 cycles 才能完成？
第三個結果是用 SSE 計算 Taylor series，也就是第二個程式片斷。它比 x87 版本快了近四倍，平均每個數字只需要 17 cycles。可是，因為 SSE 指令一次對四個數字進行運算，所以其實每四個數字共花了 68 cycles。這裡也是同樣的問題：為什麼 16 個計算會需要 68 cycles？這表示平均一個 SSE 乘法指令或加法指令，花了四個 cycles 執行。但是，理論上 Pentium III 不是可以每 cycle 執行一個 SSE 指令嗎？
也許問題是出在上面的程式片斷裡面的迴圈（for(i = 7; i >=2; i--)）。所以，也許可以把它展開試試。它會得到上面的第四個結果，實際上也是快了不少，平均每個四個數字花了約 43 cycles。也就是說，平均每個運算花了 2.7 cycles。這樣已經算是相當不錯了。可是，它和「每 cycle 執行一個 SSE 指令」還有一段距離。
事實上，這是因為這些指令都有相依性。在上面的程式片斷中，result 這個變數不但從頭用到尾，而且每個指令都用到這個變數，也會修改這個變數。因此，第二個指令要等到第一個指令執行完後，才能執行。在實際執行時，因為 Pentium III 的 instruction reorder buffer 夠大，所以可以達到一部份程度的平行執行，但是效率還不是非常理想。
要儘可能達到平行執行的效率，就要用另一個方法來展開迴圈。也就是說，不展開內部的迴圈，而是一次對多組資料進行計算。因為多組資料的計算之間並沒有相依性，要平行執行就容易多了。上圖中的第五個結果，就是一次對四組數字進行計算（也就是共 16 個數字），而內部的迴圈並沒有展開。這樣平均每四個數字只需要 34 cycles，也就是平均每個運算花了約 2.1 cycles。比起第四個結果，又要再快了一些。
那如果一次對八組數字進行計算，會再變快嗎？答案是不會，它反而會變慢。為什麼呢？因為 SSE 只有八個浮點數暫存器，所以如果一次對八組數字進行運算，那一定會有些東西無法放到暫存器中，而需要放到記憶體裡面。雖然這些資料幾乎是一定會放在快速的 L1 cache 中，但是它還是會需要一些額外的 load/store 指令。所以反而不會增加效率。不過，我們還是可以把內部的迴圈也展開，因為這個迴圈並不大，所以展開並不會太誇張。把內部的迴圈也展開後，得到的就是第六個結果，也就是平均每四個數字只需要 22.6 cycles，平均每個運算只需要 1.4 cycles！這已經離「每 cycle 執行一個 SSE 指令」不遠了！
還有辦法再加快嗎？恐怕不容易了。因為它已經達到約 2,837MFLOPS（在 Pentium III 1.0B Ghz 上執行）。事實上，這個結果比前面計算向量內積的速度還要快。所以，要再加快應該是很難了。
所以，這裡要討論的重點就是，當計算相當冗長，而且相依性很高的時候，光是展開內部的迴圈是不夠的，應該要考慮也展開外部的迴圈（也就是一次對多組資料進 行同樣的計算）。這樣通常比展開內部的迴圈要更有效。當然，如果情況許可的話，也可以連內部的迴圈也展開，通常可以達到更高的效率。
另外，可能會有人問：這樣計算自然指數，精確度有多高呢？和用 exp 函數比較又是如何？在這個例子中，1,024 個隨機數字中，誤差都小於 1.2E-7，也就是單精確度浮點數的 epsilon。這表示有誤差的情形，都只有最後一個位元不同。因此，這應該已經是非常高的精確度了。
扫一扫
1141
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
