Chaos on Graphics
About 數學庫之一 SSE/SSE2
Chaos Chiao
毫無疑問，數學庫是圖形程序的基石，是圖形程序運行效率的關鍵之一。一個優秀的數學庫可以讓圖形程序運行得更流暢，甚至要快上幾十倍上百倍。有時候替換一條除法運算會帶來成倍的效率增長，比如用乘以 1/op 替換 vector 里的 operator /。當然，更高級的優化是使用 SIMD 優化海量運算，這就是本文的中心——SSE/SSE2 優化。
在描述 SSE/SSE2 優化前，我先介紹一般的 vector/matrix 庫構造。當然，在 OpenEXR 里已經有一個非常優秀的 Imath 實現了，數學庫的實現細節可以參照它。
在圖形程序里我們經常會遇到向量運算，這是標準C++編譯器所不能直接支持的，如三維空間向量。傳統的C圖形程序會使用“數組+宏”的實現方式：
typedef float vector[3];
到了C++時代，一般會封裝成：
class Vector
{
private:
float x , y , z;
};
然后加入通常的各種method，如
const float& X( void ) const {
return x;
}
標準的向量算法如內積、外積、單位化、長度、運算等等，都可以封裝為成員函數。
Vector operator + ( const Vector& a , const float & b ) {
return Vector( a.x + b , a.y + b , a.z + b );
}
類似的數學庫可以在Aqsis等一些開源的圖形程序里找到。不過這些結構并不適合接下來我們要討論的SSE/SSE2優化。
SSE – Streaming SIMD Extension，是Intel從PIII開始加入的一種x86擴展指令集。在SSE以前，x86的浮點運算都是以棧式FPU完成的，有一定x86匯編經驗的人應該不會對那些復雜的fld、fst指令陌生吧。而SSE一方面讓浮點運算可以像整數運算的模式、如 add eax , ebx 那樣通過直接訪問寄存器完成，繞開了討厭的棧，另一方面引入了SIMD這個概念。SIMD – Single Instruction Multiply Data，顧名思義，它可以同時讓一條指令在多個數據上執行，這種體系結構在一度在大型機上非常流行，需要經常進行海量運算的大型機器通常會通過一個數學SIMD虛擬機加快處理速度，比如同時讓一組數據執行一個變換，數據的規模有上百萬之巨，而SIMD則可以優化數據的存儲與運算，減免某些切換Context的開銷。
在硬件層面上面支持SIMD，某程度是因游戲需要的驅使，因為越來越多的3D游戲涉及大量的向量操作，一般的浮點運算優化已經不能再適應這種并行運算的需要了，而直接從指令上支持SIMD操作則可以進一步簡化向量運算的優化，提高指令執行效率。像 addps 這樣的SSE指令，可以并行執行四個32位浮點數的加法運算，而延遲只有4 cycle；相比之下，原來的fadd指令光執行一個32位單精度浮點數加法的延遲已經達到了3 cycle了，還沒計算fst等存儲指令的延遲。（具體見后面的指令執行單元表）
顯然，SSE能給圖形程序帶來極大的優化，其提高遠勝于基于整數的MMX與雙單元單精度浮點數的3DNow!。但SSE對數據組織的要求是苛刻的，若要發揮SSE的最大威力，我們還需要進行對齊向量數據，把向量對齊到16字節。如果我們正在使用一般的三分量向量，那么就意味著有要浪費四分之一的存儲空間來換取速度。當然，這４字節還可以有很多用途，只是你必須處理得非常小心，因為任何運算都將同時應用到四個分量上。
要使用SSE，必須先確認你的編譯器是否支持新的指令集。VC6 sp6、VC.net、.net 2003、ICL、GCC 、nasm 都支持SSE指令集。我推薦使用ICL，它的優化做得最棒，生成的指令最緊湊、效率最高。使用SSE有兩種途徑，一是直接編寫匯編代碼，但難度較大，需要有一定的匯編經驗；二是使用SSE intrinsic，一種直接在C/C++里使用SSE指令的偽函數調用。在圖形運算的核心環節上、如raytrace核心，我建議使用匯編，這樣才能極大地體現出SSE的優勢、與x86指令混合使用，并充分使用它的并行性。而在大多數場合下則推薦使用intrinsic，它的可讀性高，而且編譯器會在最后把函數調用替換成SSE指令，這樣既不需要寫內嵌匯編代碼，又可以保證代碼的執行效率。
下面將通過幾個簡單的運算例子介紹SSE intrinsic的使用。首先，使用SSE需要一個新的頭文件
#include
里面定義了一個新的數據類型，__m128，這是一個128位、４個32位單精度浮點數的結構，如果你正在使用VC.net，你會看到它是一個關鍵字，被當作一種基本數據類型。要是你不打算使用匯編SSE，那么就沒必要深究編譯器在幕后到底如何處理__m128類型的數據，你只需要知道里面能存放四個float，而這四個float可以進行并行運算。
在定義了__m128后，文件聲明一大堆對__m128進行運算的函數，如_mm_add_ps、_mm_sub_ps等等，這就是SSE運算指令的聲明。使用SSE優化在這些聲明的幫助下變得非常簡單，如計算兩個向量之和，平時需要每一個元素進行一次加法運算，現在只需要簡單地：
__m128 a , b , c;
c = _mm_add_ps( a , b );
這樣等價于：
float a[4] , b[4] , c[4];
for( int i = 0 ; i < 4 ; ++ i )
c[i] = a[i] + b[i];
但前者的運算是并行的，在一般情況下效率遠比后者要高。況且在描述復雜的運算的時候，如：
a = b * c + d / e;
則可以直接寫成：
__m128 a = _mm_add_ps( _mm_mul_ps( b , c ) , _mm_div_ps( d , e ) );
咋看之下，很多效率至上的人馬上就會大叫“昂貴的函數調用啊！Bad smell code!”。其實我正要告訴你，我也是效率至上派的。前面已經說過了，這些看上去貌似函數的調用實際上并非函數，而是所謂intrinsic，它們在編譯優化中將被解釋為單條或多條SSE指令，而且編譯器會自動調節調用順序以使其最大并行效率。
不過除了直接使用這些intrinsic以外，我們還可以把它們封裝到類里面，重載運算符，這樣就可以把運算寫成可讀性更強的算術式。如果你不愿意自己動手封裝，也可以使用Intel封裝好了的F32vec4類，它提供了完備的運算符重載，完全使用SSE，非常方便。
雖然Intel封裝好的類已經很完善了，但還有一大堆數學運算需要我們自己動手進行編寫，如內積（點積）和外積（叉積）。
首先來看一個比較實用的運算，求倒數。求倒數在很多數學庫里都有專門的優化，通常原理都是先求出一個近似值，然后通過Newton-Raphson逼近法求出較精確值，下面的代碼摘自NV的fastmath.cpp:
#define FP_ONE_BITS 0x3F800000
// r = 1/p
#define FP_INV(r,p) \
{ \
int _i = 2 * FP_ONE_BITS - *(int *)&(p); \
r = *(float *)&_i; \
r = r * (2.0f - (p) * r); \
}
而在SSE里也提供了兩條求倒數的指令rcpss/rcpps（對應的intrinsic是_mm_rcp_ss與_mm_rcp_ps），不過這兩條指令求的并非是精確值，而是近似值，所以我們需要對它的結果進行逼近處理。
float __rcp<float>( const float& a ) {
register float r;
__m128 rcp = _mm_load_ss( &a );
rcp = _mm_rcp_ss( rcp );
_mm_store_ss( &r , rcp );
/* [2 * rcpps(x) - (x * rcpps(x) * rcpps(x))] */
r = 2.0f * r - ( a * r * r );
return r;
}
原理一致，只不過我們還可以用_mm_rcp_ps并行求四分量的倒數。如果你還對SSE的威力有所保留，那我建議你設計一個測試單元測試一下使用除法求倒數與使用SSE求倒數，看效率到底是誰更高、高多少。當然，我自己已經測試過很多次了J。
然后我們把注意力放到一條非常特殊的指令shufps（對應intrinsic是_mm_shuffle_ps）上面。這是一條非常有用的指令，它可以把兩個操作數的分量以特定的順序排列并賦予給目標數。比如
__m128 b = _mm_shuffle_ps( a , a , 0 );
則 b 的所有分量都是 a 中下標為0的分量。第三個參數控制分量分配，是一個8bit的常量，這個常量的1~8位分別控制了從兩個操作數中選擇分量的情況，具體怎么控制將在后面討論SSE匯編中一并說明，而在使用intrinsic的時候，最好使用_MM_SHUFFLE宏，它可以定義分配情況。下面我們來復習一下叉積的求法。
c = a x b
可以寫成：
Vector cross(const Vector& a , const Vector& b ) {
return Vector(
( a[1] * b[2] - a[2] * b[1] ) ,
( a[2] * b[0] - a[0] * b[2] ) ,
( a[0] * b[1] - a[1] * b[0] ) );
}
那么寫成SSE intrinsic形式則是：
/* cross */
__m128 _mm_cross_ps( __m128 a , __m128 b ) {
__m128 ea , eb;
// set to a[1][2][0][3] , b[2][0][1][3]
ea = _mm_shuffle_ps( a , a , _MM_SHUFFLE( 3 , 0 , 2 , 1 ) );
eb = _mm_shuffle_ps( b , b , _MM_SHUFFLE( 3 , 1 , 0 , 2 ) );
// multiply
__m128 xa = _mm_mul_ps( ea , eb );
// set to a[2][0][1][3] , b[1][2][0][3]
a = _mm_shuffle_ps( a , a , _MM_SHUFFLE( 3 , 1 , 0 , 2 ) );
b = _mm_shuffle_ps( b , b , _MM_SHUFFLE( 3 , 0 , 2 , 1 ) );
// multiply
__m128 xb = _mm_mul_ps( a , b );
// subtract
return _mm_sub_ps( xa , xb );
}
這就是shuffle強大的地方，它可以直接在寄存器里直接調整分量的順序。而且配合_mm_movehl_ps，我們可以輕松解決點積的運算。_mm_movehl_ps把操作數高位兩個分量賦予目標數的低位兩分量，而目標數的高位兩分量值不變，相當于：
a[0] = b[2];
a[1] = b[3];
三分量的向量求點積，可以寫成：
float dot( const float& a , const float& b ) const {
return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
}
則用SSE intrinsic可以寫成：
/* x[0] * x[1] + y[0] * y[1] + z[0] * z[1] */
__m128 _mm_dot_ps( __m128 x , __m128 y ) {
__m128 s , r;
s = _mm_mul_ps( x , y );
r = _mm_add_ss( s , _mm_movehl_ps( s , s ) );
r = _mm_add_ss( r , _mm_shuffle_ps( r , r , 1 ) );
return r;
}
通過這兩個例子，可以留意到向量內元素的垂直相加一般形式，即：
/* x[0] + x[1] + x[2] + x[3] */
__m128 _mm_sum_ps( __m128 x ) {
__m128 r;
r = _mm_add_ps( x , _mm_movehl_ps( x , x ) );
r = _mm_add_ss( r , _mm_shuffle_ps( r , r , 1 ) );
return r;
}
那么通過擴展，可以得到求向量長度的函數，首先是求分量平方和函數：
/* x[0] * x[0] + y[0] * y[0] + z[0] * z[0] */
__m128 _mm_square_ps( __m128 x ) {
__m128 s , r;
s = _mm_mul_ps( x , x );
r = _mm_add_ss( s , _mm_movehl_ps( s , s ) );
r = _mm_add_ss( r , _mm_shuffle_ps( r , r , 1 ) );
return r;
}
然后就可以直接把結果求平方根，可得長度。解決了長度，接下來則是很重要的單位化了。可以說單位化是最重要的一個函數，它經常被調用到，而函數內的陷阱卻又最多。求單位化其實并不難，就是分量除以向量長度，可以寫成：
void normalize( const Vector& a ) {
float len = a[0] * a[0] + a[1] * a[1] + a[2] * a[2];
if( is_zero( len ) )
return;
len = 1 / len;
a[0] *= len;
a[1] *= len;
a[2] *= len;
}
我和這個家伙打交道已經有差不多七年時間了，所以脾性非常熟悉。首先求分量的平方和，判斷是否為0（問我為什么不直接用 if( len == 0 )？好樣的，請先去復習一下浮點數的基本知識），然后再求倒數，最后反映到分量上。在把它寫成SSE intrinsic格式前，我先引入另外一個能極大提升運算效率的函數，求平方根的倒數。有數值運算編成經驗的人都知道，如果說除法是惡魔的話，那么平方根就是撒旦了，而平方根的倒數簡直就是撒旦他媽。雖然上面提供了倒數的逼近方法，但僅僅使用它還是繞不開最主要的開銷、平方根運算。幸好，SSE提供了一個直接計算平方根倒數近似值的指令，rsqrtss/rsqrtps（即_mm_rsqrt_ss和_mm_rsqrt_ps）。照搬倒數求法，可以輕松得出：
/* r = 1 / sqrt(a) */
/* 0.5 * rsqrtss * (3 - x * rsqrtss(x) * rsqrtss(x)) */
__m128 _mm_rsqrt( __m128 a )
{
/／ divisor
static const __m128 _05 = _mm_set1_ps( 0.5f );
static const __m128 _3 = _mm_set1_ps( 3.f );
__m128 rsqrt = _mm_rsqrt_ss( a );
rsqrt =
_mm_mul_ss(
_mm_mul_ss( _05 , rsqrt ) ,
_mm_sub_ss( _3 , _mm_mul_ss( a , _mm_mul_ss( rsqrt , rsqrt ) ) ) );
return rsqrt;
}
那么就可以輕松得出單位化向量的函數了：
// normalize & return value
__m128 _mm_normalize( const __m128 a ) {
// get length square
__m128 l = _mm_square_ps( a );
// test if length is zero
if( _mm_iszero_ss( l ) )
return z;
// length inverse
l = _mm_rsqrt( l );
// shuffle
l = _mm_shuffle_ps( l , l , 0 );
// multiply to vector
return _mm_mul_ps( a , l );
}
SSE除了以上這些數學運算操作外，還提供了位運算。位運算？想到什么了嗎？對！比較與選擇。首先來看一個最簡單的，求絕對值。通常我們會把 abs 寫成非常簡潔的形式：
float abs( float a ) {
a >= 0 ? a : -a;
}
但當我們已經Pack了一個向量到__m128結構里，而又不希望Unpack他們進行浮點數的比較，那么就可以使用SSE的位操作。
/* abs */
__m128 _mm_abs_ps( __m128 a )
{
static const union { int i[4]; __m128 m; } __mm_abs_mask_cheat_ps
= {0x7fffffff, 0x7fffffff, 0x7fffffff, 0x7fffffff};
return _mm_and_ps( a, __mm_abs_mask_cheat_ps.m );
}
還記得單精度浮點數的符號存放在什么位上面嗎？我們只需把它除掉，然后就可以很輕松地得到了正值了。
圖形程序很多時候會用到32位浮點色彩，其值域通常為[0,1]，所以clamp函數出現的頻率也十分頻繁。要將rgba的值同時clamp到值域內，毫無疑問，SSE的并行特性又得到了發揮的機會。先來看cut函數，它負責把超出值域的值干掉，但為了更靈活，我們一次只cut一邊的區間，所以cut有兩兄弟，分別是locut和hicut。
__m128 _mm_locut_ps( __m128 val , __m128 bound )
{
__m128 mask = _mm_cmplt_ps( val , bound );
return _mm_or_ps( _mm_and_ps( mask , bound ) , _mm_andnot_ps( mask , val ) );
}
__m128 _mm_hicut_ps( __m128 val , __m128 bound )
{
__m128 mask = _mm_cmpgt_ps( val , bound );
return _mm_or_ps( _mm_and_ps( mask , bound ) , _mm_andnot_ps( mask , val ) );
}
_mm_cmp**_ps是一系列的比較函數，非常豐富，也很好用，如果替換成相應的if-else，并行性將被大大削弱。不過_mm_cmp**_ps的最大缺點就是靈活度不夠，返回值是一系列位標記，其具體的用法將在SSE匯編中討論。有了這倆哥們，clamp變得非常簡單：
__m128 _mm_clamp_ps( __m128 val , __m128 min , __m128 max )
{
return _mm_locut_ps( _mm_hicut_ps( val , max ) , min );
}
以上只是一些很簡單的實現，利用了SSE intrinsic對數學運算進行優化，并盡可能地不分拆向量，這樣可以保證8個128位的xmm寄存器可以滿足大部分運算。不過SSE intrinsic始終受編譯器生成代碼的質量好壞影響，沒能真正發揮出SSE的全部威力，接下來我們將討論SSE匯編的用法與優化。
to be continued...