sse 指令最佳化叉乘

sse referece : http://docs.sun.com/app/docs/doc/817-5477/6mkuavhrm?a=view

 

下文轉自：http://goutie.blog.sohu.com/82010206.html

 

4.2.3 叉乘
　　另外一個向量計算中經常要實現的就是叉乘（譯者註：呵呵，忘的差不多了，與向量內積不同，叉乘的結果為一個向量，詳細請見空間解析幾何相關知識。）。如果在三維空間有兩個任意向量向量，通過叉乘操作就可以得到第三個向量，有意思的是，如果這兩個向量不平行，那麼這第三條向量就會與另外兩條向量垂直。
　　需要該操作的一個樣本是求一個平面或多邊形上的法線向量，以便計算光影效果。再一次，SSE將在C++面前展現效率的魔力。
　　inline void ZFXVector::Cross(const ZFXVector &u, const ZFXVector &v)
　　{
　　  if(!g_bSSE)
　　  {
　　    x = u.y * v.z - u.z * v.y;
　　    y = u.z * v.x - u.x * v.z;
　　    z = u.x * v.y - u.y * v.x;
　　    w = 1.0f
　　    // 譯者註：以上對照叉乘公式可得。另pdf中，u為v1，v為v2，懷疑作者為著作權之故是故意為之，請多加註意！！！
　　    //         同時，也請自學的朋友好自為之。
　　  }
　　  else
　　  {
　　    __asm {
　　      mov esi, u
　　      mov edi, v
　　     
　　      movups xmm0, [esi]
　　      movups xmm1, [edi]
　　      movaps xmm2, xmm0
　　      xovaps xmm3, xmm1
　　     
　　      shufps xmm0, xmm0, 0xc9
　　      shufps xmm1, xmm1, 0xd2
　　      mulps xmm0, xmm1
　　     
　　      shufps xmm2, xmm2, 0xd2
　　      shufps xmm3, xmm3, 0xc9
　　      mulps xmm2, xmm3
　　     
　　      subps xmm0, xmm2
　　      mov esi, this
　　      movups [esi], xmm0
　　    }
　　    w = 1.0f
　　  }
　　}
　　
　　以上代碼中，首先要把兩個需要做叉乘的向量，每個都要分別放入兩個XMM寄存器中，因為我們必須對其中一個寄存器中值進行混排，然後再乘到另一個寄存器中。因此，最終的乘積結果位於XMM0及XMM2中。如果對叉乘沒啥概念，還好前面提供了對應的C++代碼，希望你能看的明白:)，哦，還有最後一步，用SUBS實現兩個向量的減法運算，此時XMM0中的才是最終的叉乘結果。當然，記得把結果拷貝回記憶體中儲存，以備不時之需。
　　注意！注意！這裡的混排控制值是以0x打頭的，而不是原來的h尾碼。如果不是那樣的話，編譯器識別十六進位時，就有大麻煩了:)（譯者註：注意啊，我記得前面有一段代碼中就用了h尾碼，我想這對於真正的程式員來說，應該不是什麼問題。稍帶說一句，任何問題的解決都要靠自己的努力，而不是隨便一個問題，自己沒怎麼剖析就問別人，那樣一個人永遠成長不了。尋求協助有時是必要的，但不是最關鍵的。）
　　
　　4.2.4 向量矩陣乘法
　　
　　我們經常會需要用一個4x4矩陣乘以一個向量，因此，我們也需要一個快些的操作實現。在這裡，我們先用一下ZFXMatrix庫中所實現的4x4矩陣類，以便能簡要說明這一操作實現。
　　還有，我假定你懂得向量與矩陣的乘法知識，否則，還是通過以下函數中的C++代碼學習吧。
　　ZFXVector ZFXVector::operator * (const ZFXMatrix &m) const
　　{
　　  ZFXVector vcResult;
　　 
　　  if(!g_bSSE)
　　  {
　　    vcResult.x = x*m._11 + y*m._21 + z*m._31 + m._41;
　　    vcResult.y = x*m._12 + y*m._22 + z*m._32 + m._42;
　　    vcResult.z = x*m._13 + y*m._23 + z*m._33 + m._43;
　　    vcResult.w = x*m._14 + y*m._24 + z*m._34 + m._44;
　　   
　　    vcResult.x = vcResult.x/vcResult.w;
　　    vcResult.y = vcResult.y/vcResult.w;
　　    vcResult.z = vcResult.z/vcResult.w;
　　    vcResult.w = 1.0f
　　  }
　　  else
　　  {
　　    float *ptrRet = (float *)&vcResult;
　　    _asm {
　　      mov ecx, this     ; 向量
　　      mov edx, m        ; 矩陣
　　      movss xmm0, [ecx]
　　      mov eax, ptrRet   ; 結果向量
　　      shufps xmm0, xmm0, 0
　　      movss xmm1, [ecx+4]
　　      mulps xmm0, [edx]
　　      shufps xmm1, xmm1, 0
　　      movss xmm2, [ecx+8]
　　      mulps xmm1, [edx+16]
　　      shufps xmm2, xmm2, 0
　　      movss xmm3, [ecx+12]
　　      mulps xmm2, [edx+32]
　　      shufps xmm3, xmm3, 0
　　      addps xmm0, xmm1
　　      mulps xmm3, [edx+48]
　　      addps xmm2, xmm3
　　      addps xmm0, xmm2
　　      movups [eax], xmm0  ;儲存結果
　　      mov [eax+3], 1      ; w= 1
　　    }
　　  }
　　  return vcResult;
　　}
　　
　　我們前面所學到的SSE知識，足夠能使你實現一個向量和矩陣的乘法。看上去有很多混排指令，其實只是廣播而已。
　　呵呵，到目前可以說，你已經開闢了自己的前進之路，圍繞這向量類一個個函數的建立，漸漸精通了SSE指令及其實現。需要牢記的是，使用過於頻繁的部分值得我們去最佳化。“過於頻繁”意味著也許一幀之中你會調用到上百次。因此，不應為向量的SSE實現而驚訝，雖然有點複雜，但是的確管用。好的，關於向量就介紹到這裡，現在我們可以開始關注矩陣了吧？！

 

下面的轉自：http://www.gameres.com/document.asp?TopicID=84488

高效率的3D圖形數學庫(1) ----Vector概覽
潛水了很長時間，該是做點貢獻的時候了，最近寫的，發上來給各位拍磚：

 

    最近研究彙編比較多，看自己C++代碼的彙編源碼簡直是一種折磨，這迫使我將所有數學庫重新用彙編指令實現，當然，包括對CPUID的檢測和使用擴充指令集。測試結果是與D3DX9的數學函數比較的，效果另人滿意，除了矩陣相乘的演算法總是與D3DXMatrixMultiply函數有7%的差距外，其餘都是持平甚至遙遙領先（也許是我瘋了，有新的看官可以自己測一下）。 由於本人技術淺薄，測試效率的方法又比較簡陋，所以還請高手指正！
第一步是介紹我的Vector類，以下是聲明：

struct __declspec(dllexport) Vector
{

/******************變數********************/

float x, y, z, w;

/******************構造*******************/

// 建構函式
Vector() {}
// 建構函式
Vector(const float* v);
// 建構函式
Vector(float _x, float _y, float _z, float _w);

/******************方法*******************/

// 設定向量
void SetVector(const float* v);
// 設定向量
void SetVector(float _x, float _y, float _z, float _w);
// 減法
void Difference(const Vector* pSrc, const Vector* pDest);
// 反向量
void Inverse();
// 單位化向量
void Normalize();
// 是否單位向量
bool IsNormalized();
// 向量長度(慢)
float GetLength();
// 向量長度的平方(快)
float GetLengthSq();
// 通過兩向量求叉乘,結果儲存在該向量中
void Cross(const Vector* pU, const Vector* pV);
// 求兩向量夾角
float AngleWith(Vector& v);

/*************運算子多載*****************/

// 運算子多載
void operator += (Vector& v);
// 運算子多載
void operator -= (Vector& v);
// 運算子多載
void operator *= (float v);
// 運算子多載
void operator /= (float v);
// 運算子多載
Vector operator + (Vector& v) const;
// 運算子多載
Vector operator - (Vector& v) const;
// 運算子多載
float operator * (Vector& v) const;
// 運算子多載
void operator *= (GaiaMatrix& m);
// 運算子多載
Vector operator * (float f) const;
// 運算子多載
bool operator ==(Vector& v);
// 運算子多載
bool operator !=(Vector& v);
// 運算子多載
//void operator = (Vector& v);
};

然後是簡單的內嵌函式:

// 建構函式
inline Vector::Vector(const float* v)
: x(v[0])
, y(v[1])
, z(v[2])
, w(v[3])
{
}

// 建構函式
inline Vector::Vector(float _x, float _y, float _z, float _w)
: x(_x)
, y(_y)
, z(_z)
, w(_w)
{
}

// 設定向量
inline void Vector::SetVector(const float* v)
{
x = v[0]; y = v[1]; z = v[2];
}

// 設定向量
inline void Vector::SetVector(float _x, float _y, float _z, float _w)
{
x = _x; y = _y; z = _z; w = _w;
}

// 減法
inline void Vector::Difference(const Vector* pSrc, const Vector* pDest)
{
x = pDest->x - pSrc->x;
y = pDest->y - pSrc->y;
x = pDest->z - pSrc->z;
}

// 反向量
inline void Vector::Inverse()
{
x = -x; y = -y; z = -z;
}

// 是否單位向量
inline bool Vector::IsNormalized()
{
return CmpFloatSame(x*x+y*y+z*z, 1.0f);
}

// 運算子多載
inline void Vector::operator += (Vector& v)
{
x += v.x; y += v.y; z += v.z;
}
// 運算子多載
inline void Vector::operator -= (Vector& v)
{
x -= v.x; y -= v.y; z -= v.z;
}
// 運算子多載
inline void Vector::operator *= (float f)
{
x *= f; y *= f; z *= f;
}
// 運算子多載
inline void Vector::operator /= (float f)
{
f = 1.0f/f;
x *= f; y *= f; z *= f;
}
// 運算子多載
inline Vector Vector::operator + (Vector& v) const
{
return Vector(x+v.x, y+v.y, z+v.z, w);
}
// 運算子多載
inline Vector Vector::operator - (Vector& v) const
{
return Vector(x-v.x, y-v.y, z-v.z, w);
}
// 運算子多載
inline float Vector::operator * (Vector& v) const
{
return (x*v.x + y*v.y + z*v.z);
}
// 運算子多載
inline Vector Vector::operator * (float f) const
{
return Vector(x*f, y*f, z*f, w);
}
// 運算子多載
inline bool Vector::operator ==(Vector& v)
{
return ((((x-v.x)<FLOAT_EPS && (x-v.x)>-FLOAT_EPS) || ((y-v.y)<FLOAT_EPS && (y-v.y)>-FLOAT_EPS) || ((z-v.z)<FLOAT_EPS && (z-v.z)>-FLOAT_EPS))? false:true);
}
// 運算子多載
inline bool Vector::operator !=(Vector& v)
{
return ((((x-v.x)<FLOAT_EPS && (x-v.x)>-FLOAT_EPS) || ((y-v.y)<FLOAT_EPS && (y-v.y)>-FLOAT_EPS) || ((z-v.z)<FLOAT_EPS && (z-v.z)>-FLOAT_EPS))? true:false);
}

這裡比較重要的最佳化有幾點，也可以作為寫代碼的原則，非常非常重要：

1、可以用const的地方一定要用！編輯器會拿這個來最佳化的。
2、return返回一個值的時候，如果可以的話，就一定要以建構函式的形式傳回值。如：
return Vector(x+v.x, y+v.y, z+v.z, w);
3、多個數除以同一個數時，一定要按照如Vector::operator /= (float f)中的形式寫。
4、這樣的小函數一定是要inline的！

以上4點一定要遵守，否則做出的彙編代碼慘不忍睹！效率自然也是一落千丈,切記切記。

接下來是Vector的進階函數部分：

// 向量長度的平方(快)
float Vector::GetLengthSq() // 潛在危險
{
_asm
{
fld         dword ptr [ecx];
fmul        dword ptr [ecx];
fld         dword ptr [ecx+4];
fmul        dword ptr [ecx+4];
faddp       st(1),st;
fld         dword ptr [ecx+8] ;
fmul        dword ptr [ecx+8] ;
faddp       st(1),st ;
}
//return x*x + y*y + z*z;
}

// 向量長度(慢)
float Vector::GetLength()
{
float f;
if (g_bUseSSE2)
{
_asm
{
lea ecx, f;
mov eax, this;
mov dword ptr [eax+12], 0; // w = 0.0f;

movups xmm0, [eax];
mulps xmm0, xmm0;
movaps xmm1, xmm0;
shufps xmm1, xmm1, 4Eh; 洗牌
addps xmm0, xmm1;
movaps xmm1, xmm0;
shufps xmm1, xmm1, 11h; 洗牌
addss xmm0, xmm1;

sqrtss xmm0, xmm0; 第一個單元求開方
movss dword ptr [ecx], xmm0; 第一個單元的值給ecx指向的記憶體空間

mov dword ptr [eax+12], 3F800000h; // 3F800000h == 1.0f
}
}
else
{
f = (float)sqrt(x*x+y*y+z*z);
}
return f;
}

// 單位化向量
void Vector::Normalize()
{
if (g_bUseSSE2)
{
_asm
{
mov eax, this;
mov dword ptr[eax+12], 0;

movups xmm0, [eax];
movaps xmm2, xmm0;
mulps xmm0, xmm0;
movaps xmm1, xmm0;
shufps xmm1, xmm1, 4Eh;
addps xmm0, xmm1;
movaps xmm1, xmm0;
shufps xmm1, xmm1, 11h;
addps xmm0, xmm1;

rsqrtps xmm0, xmm0;
mulps xmm2, xmm0;
movups [eax], xmm2;

mov dword ptr [eax+12], 3F800000h;
}
}
else
{
float f = (float)sqrt(x*x+y*y+z*z);
if (f != 0.0f)
{
f = 1.0f/f;
x*=f; y*=f; z*=f;
}
}
}

// 通過兩向量求叉乘,結果儲存在該向量中
void Vector::Cross(const Vector* pU, const Vector* pV)
{
if (g_bUseSSE2)
{
_asm
{
mov eax, pU;
mov edx, pV;

movups xmm0, [eax]
movups xmm1, [edx]
movaps xmm2, xmm0
movaps xmm3, xmm1

shufps xmm0, xmm0, 0xc9
shufps xmm1, xmm1, 0xd2
mulps xmm0, xmm1

shufps xmm2, xmm2, 0xd2
shufps xmm3, xmm3, 0xc9
mulps xmm2, xmm3

subps xmm0, xmm2

mov eax, this
movups [eax], xmm0

mov [eax+12], 3F800000h;
}
}
else
{
x = pU->y * pV->z - pU->z * pV->y;
y = pU->z * pV->x - pU->x * pV->z;
z = pU->x * pV->y - pU->y * pV->x;
w = 1.0f;
}
}

// 運算子多載
void Vector::operator *= (Matrix& m) // 潛在危險
{
#ifdef _DEBUG
assert(w!=1.0f && w!=0.0f);
#endif

if (g_bUseSSE2)
{
_asm
{
mov ecx, this;
mov edx, m;
movss xmm0, [ecx];
//lea eax, vr;
shufps xmm0, xmm0, 0; // xmm0 = x,x,x,x

movss xmm1, [ecx+4];
mulps xmm0, [edx];
shufps xmm1, xmm1, 0; // xmm1 = y,y,y,y

movss xmm2, [ecx+8];
mulps xmm1, [edx+16];
shufps xmm2, xmm2, 0; // xmm2 = z,z,z,z

movss xmm3, [ecx+12];
mulps xmm2, [edx+32];
shufps xmm3, xmm3, 0; // xmm3 = w,w,w,w

addps xmm0, xmm1;
mulps xmm3, [edx+48];

addps xmm0, xmm2;
addps xmm0, xmm3; // xmm0 = result
movups [ecx], xmm0;
mov [ecx+12], 3F800000h;
}

}
else
{
Vector vr;
vr.x = x*m._11 + y*m._21 + z*m._31 + w*m._41;
vr.y = x*m._12 + y*m._22 + z*m._32 + w*m._42;
vr.z = x*m._13 + y*m._23 + z*m._33 + w*m._43;
vr.w = x*m._14 + y*m._24 + z*m._34 + w*m._44;

x = vr.x;
y = vr.y;
z = vr.z;
w = 1.0f;
}
}

// 求兩向量夾角
float Vector::AngleWith(Vector& v)
{
return (float)acosf((*this * v)/(this->GetLength()*v.GetLength()*2.0f));
}

這裡要說明3個函數：GetLengthSq，*= 和AngleWith
    GetLengthSq有潛在危險，因為我是根據.Net2003的編輯器來寫的代碼，我知道ecx==this，知道float的傳回值是直接從浮點棧寄存器fstp到外面參數的，所以，我會用這種方法來寫，甚至沒有寫傳回值！而看此文的您可能不會使用與我一樣的編輯器，所以，在理解了實質之後，運用合理的演算法來實現你的數學庫。後面的函數都使用了編輯器無關的方法寫的。

    *= 的運算子多載的潛在危險在於，Vector是4D的，可以表示3D的向量或者3D空間點座標。如果是向量，則w==0，這樣就只會受到旋轉和縮放的影響。而如果是表示空間點，w==1，就會受到所有類型的變動，如平移、旋轉和縮放。由於向量是不能平移的，處於對運算效率的考慮，這時候就需要數學庫的調用者自己注意了。

    AngleWith函數之所以不對其進行內聯化，是因為在以後的文章中，我會去進一步最佳化這裡的代碼。GetLength和acosf都不是內嵌函式，我必須要將其展開，以彙編實現，並重新組織編碼。這個函數好像在D3DX9的數學庫中是沒有的~~沒辦法比較了。

以上幾個函數的效率與D3DX庫比較結果大致是這樣的：
    GetLengthSq微高於D3DX
    GetLength是D3DX速度的2倍多，因為D3D庫沒有用SSE指令。
    Normalize和Cross的速度比D3DX的高的太多，有些離譜。同樣是因為D3D庫沒有用SSE指令。
    *=的效率低於D3DXVec3Transform約7%，有進一步提高的可能！高手來看看。D3DX庫用的是3DNow!運算的，居然比SSE快！大概是因為我的AMD3000+的緣故吧...，換在Inter上應該速度差不多了。
    AngleWith沒有辦法評測，因為沒有比照對象。

    很多演算法都經過手工的指令重排，發現指令的順序對效率的影響是非常大的！在改變指令順序時一定要謹慎！最好拷貝一份原來的，否則在排比較長的彙編代碼時會把自己玩暈的~o~
    順便提幾個很多人疑惑的問題：
    1、那個C++庫裡的_mm_mov_ps()類似的代碼，簡直就是垃圾！想要效率就千萬別用那個，好好的學習彙編，然後親手寫代碼。那些庫裡的函數搞出的代碼簡直就是慘不忍睹！
    2、movups和movaps的效率差距幾乎可以忽略不計的！別為了快那麼百分之一的速度就聲明一個_m128的Vector或者Matrix，以後建立數組的時候可有你受的了！
    3、本人的測試方法太菜了，就是迴圈1000萬遍，用timeGetTime()看個大概。多運行幾遍找個平均而已。所以，一旦Release模式的內聯就測不出效率了~有時間的高人們可以去測試一下，估計能內聯的函數都是快接近效率極限的，不太值得最佳化。
如果對我的測試有什麼疑惑，看官們可以考回去自己測試效率，換多種CPU試一下，我在這兒接受任何人的拍磚