高效率的3D图形数学库

http://dev.gameres.com/

Vector概览
潜水了很长时间，该是做点贡献的时候了，最近写的，发上来给各位拍砖：  SSE与矩阵相乘
这次将介绍SSE扩展指令集，以及矩阵乘法的优化，不喜汇编者请发送WM_CLOSE消息！！
闲话就不说了，SSE指令的历史到处都是，主要说说我对指令集的原理、作用和用法的理解。

 

    最近研究汇编比较多，看自己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试一下，我在这儿接受任何人的拍砖！

    下一次我将详细说说我对SSE和浮点指令的理解，以及最有用的矩阵相乘算法。

 

SSE指令集的最大有点就是能够 4个float并行运算，与其说这恰好符合图形算法，到不如说就是为图形编程设计的。比如说，两个向量相加，x1+x2, y1+y2, z1+z2, w1+w2。用了4条加法指令，不如来一条省事。这就是SSE能为我们做的。这里要介绍8个寄存器：xmm0, xmm1, xmm2...xmm6, xmm7，这些寄存器都是128位的，每个寄存器都可以存放4个float，所以，x1, y1, z1, w1这4个float型可以放在xmm0中，而x2, y2, z2, w2可以放在xmm1中，于是 xmm0 += xmm1，结果就保存在xmm0寄存器了！下面是代码:

struct Vector
{
       float x, y, z, w;

       void Add(const Vector* pIn)
       {
          _asm
          {
             mov  eax, pIn;   // 这里其实应该是 mov eax, dword ptr[pIn];
             mov  ecx, this;

             movups  xmm0, [ecx];   // 把this的xyzw放入xmm0
             movups  xmm1, [eax];   // 把pIn的xyzw放入xmm1
             addps   xmm0, xmm1;    // xmm0 += xmm1
             movups  [ecx], xmm0;     // 把xmm0的值给this
             mov     [ecx+12], 3F800000h    // 别忘了给w = 1.0f
          }
       }
}

其实也就这么简单，movups是把内存中的向量放入寄存器，或者把寄存器中的向量放回内存，总之是一次移动128个位的一条指令,这个指令比普通的MOV指令要慢70%左右，比浮点乘法要慢的多了，大量的时间花在了movups上。所以说，简单的算法不值得去用SSE指令。

这个向量加法指令很可能会比下面的算法更慢：
void Add(const Vector* pIn)
{
    x += pIn->x;
    y += pIn->y;
    z += pIn->z;
    w += pIn->w;
}

那么，如果我要让一个向量的xyzw这4个float同时去加上同一个float值怎么办呢？比方说，我要让在xmm0中的xyzw同时加上float型变量h，就会像下面这样去做：

float h = ...;
_asm
{
movss  xmm1, h
shufps xmm1, xmm1, 0
addps  xmm0, xmm1
}

好了，出现2个新的指令： movss和shufps。

movss是将一个32位的值移动到一个128位寄存器的低32位中。也就是说，这时xmm1中的4个32位区块中只有第0个区块是存放了h的值。

这时，我想让其他3个区块也存放同样的值，以便对xmm0进行并行加法处理，于是用到了下面的指令————

shufps是用来将128位寄存器中4个区块的数值进行相互调换、拷贝或覆盖，就像洗牌一样，所以叫做洗牌指令。我们现看该指令的第3个操作数是 0，该操作数是8位的，0 = 00 00 00 00，你看，我把8位数分成了4份，每一份都表示了一个寄存器区间。
00就是第0个区块
01就是第1个区块
10就是第2个区块
11就是第3个区块

shufps  dest, src, 00 00 00 00
我们将根据第3个操作数，从src中选取区间，并把该区块中的数给拖到dest相应的区块内。这里的操作数4个都是0，所以，dest的4个区块中存放的都是src中第0个区间的值。而dest和src是同一个寄存器的时候，就会是自我洗牌。于是xmm1中的4个区块存放的都是xmm1第0个区块的值。下面的示例会帮助你更好的理解这个问题：

dest =  w1, z1, y1, x1
src  =  w2, z2, y2, x2

经过该指令： shufps  dest, src, 01 11 00 10
得到该结果： dest =  y2, w2, x2, z2
看下面的会更清晰：

dest第3区块 <<------ src中01区块(1)
dest第2区块 <<------ src中11区块(3)
dest第1区块 <<------ src中00区块(0)
dest第0区块 <<------ src中10区块(2)

注意：这里的第3操作数所代表的dest寄存器区块顺序是3210，而不是0123，x通常存放在0区块，内存中地址越小的就放在编号越小的区块~~总之别给搞反了就好，这里的确很容易混淆。

而这里其实有一个万恶的限制!!!：如果dest和src是不同的寄存器，那么3、4区块的值会取dest自己区块中的，而不是src区块中的值。所以，上面执行过的dest实际是这样的：
dest =  y1, w1, x2, z2

dest第3区块 <<------ dest中01区块(1)
dest第2区块 <<------ dest中11区块(3)
dest第1区块 <<------ src中00区块(0)
dest第0区块 <<------ src中10区块(2)

如果还没明白的话.....  可以发论坛中的短消息给我。

SSE基本指令的计算有“加减乘除”四则运算，分别是： addps, subps, mulps, divps。很好记的，都是x86基本指令后面加个ps后缀。如果是ss后缀，则代表只有第0个区间做运算，速度会比ps的快20%左右，在AMD的CPU上会比浮点指令要慢，所以，传说中的性能提升400%就是纯属口胡，因为那只是理论值。但是300%还是可以做到的，特别是复杂的算法就更容易做到，比如说矩阵相乘。

现在发一段矩阵相乘的代码，出于是重点介绍SSE指令用法的缘故，主要注意了代码的可读性，所以没有对指令进行重排，所以会比我的最优化版慢25%左右的速度，但即便如此，已然是一般算法的3倍速了，看官们可以考回去实验一下，重排后的代码会在以后放出。

void MultMatrix(const Matrix* pOut, const Matrix* pIn1, const Matrix* pIn2)
{
if (!g_bUseSSE2)
{
// [edx]   =   xmm0   *   xmm4    +    xmm1   *   xmm5    +    xmm2  *   xmm6   +    xmm3   *   xmm7
pOut->_11 = pIn1->_11*pIn2._11 + pIn1->_12*pIn2._21 + pIn1->_13*pIn2._31 + pIn1->_14*pIn2._41;
pOut->_12 = pIn1->_11*pIn2._12 + pIn1->_12*pIn2._22 + pIn1->_13*pIn2._32 + pIn1->_14*pIn2._42;
pOut->_13 = pIn1->_11*pIn2._13 + pIn1->_12*pIn2._23 + pIn1->_13*pIn2._33 + pIn1->_14*pIn2._43;
pOut->_14 = pIn1->_11*pIn2._14 + pIn1->_12*pIn2._24 + pIn1->_13*pIn2._34 + pIn1->_14*pIn2._44;

pOut->_21 = pIn1->_21*pIn2._11 + pIn1->_22*pIn2._21 + pIn1->_23*pIn2._31 + pIn1->_24*pIn2._41;
pOut->_22 = pIn1->_21*pIn2._12 + pIn1->_22*pIn2._22 + pIn1->_23*pIn2._32 + pIn1->_24*pIn2._42;
pOut->_23 = pIn1->_21*pIn2._13 + pIn1->_22*pIn2._23 + pIn1->_23*pIn2._33 + pIn1->_24*pIn2._43;
pOut->_24 = pIn1->_21*pIn2._14 + pIn1->_22*pIn2._24 + pIn1->_23*pIn2._34 + pIn1->_24*pIn2._44;

pOut->_31 = pIn1->_31*pIn2._11 + pIn1->_32*pIn2._21 + pIn1->_33*pIn2._31 + pIn1->_34*pIn2._41;
pOut->_32 = pIn1->_31*pIn2._12 + pIn1->_32*pIn2._22 + pIn1->_33*pIn2._32 + pIn1->_34*pIn2._42;
pOut->_33 = pIn1->_31*pIn2._13 + pIn1->_32*pIn2._23 + pIn1->_33*pIn2._33 + pIn1->_34*pIn2._43;
pOut->_34 = pIn1->_31*pIn2._14 + pIn1->_32*pIn2._24 + pIn1->_33*pIn2._34 + pIn1->_34*pIn2._44;

pOut->_41 = pIn1->_41*pIn2._11 + pIn1->_42*pIn2._21 + pIn1->_43*pIn2._31 + pIn1->_44*pIn2._41;
pOut->_42 = pIn1->_41*pIn2._12 + pIn1->_42*pIn2._22 + pIn1->_43*pIn2._32 + pIn1->_44*pIn2._42;
pOut->_43 = pIn1->_41*pIn2._13 + pIn1->_42*pIn2._23 + pIn1->_43*pIn2._33 + pIn1->_44*pIn2._43;
pOut->_44 = pIn1->_41*pIn2._14 + pIn1->_42*pIn2._24 + pIn1->_43*pIn2._34 + pIn1->_44*pIn2._44;
} 
else
{
_asm
{
mov edx, pIn2; // 这时保存的是pIn2
movups xmm4, [edx]; //pIn2的第1行
movups xmm5, [edx+16]; //pIn2的第2行
movups xmm6, [edx+32]; //pIn2的第3行
movups xmm7, [edx+48]; //pIn2的第4行

mov eax, pIn1; // 这时保存的是pIn1
mov edx, pOut;

mov ecx, 4; // 循环4次

LOOPIT: // 开始循环
movss xmm0, [eax]; xmm0 = pIn1->x
shufps xmm0, xmm0, 0; 洗牌xmm0 = pIn1->x, pIn1->x, pIn1->x, pIn1->x
mulps xmm0, xmm4;

movss xmm1, [eax+4]; xmm1 = pIn1->y
shufps xmm1, xmm1, 0; 洗牌xmm1 = pIn1->y, pIn1->y, pIn1->y, pIn1->y
mulps xmm1, xmm5;

movss xmm2, [eax+8]; xmm2 = pIn1->z
shufps xmm2, xmm2, 0; 洗牌xmm2 = pIn1->z, pIn1->z, pIn1->z, pIn1->z
mulps xmm2, xmm6;

movss xmm3, [eax+12]; xmm3 = pIn1->w
shufps xmm3, xmm3, 0; 洗牌xmm3 = pIn1->w, pIn1->w, pIn1->w, pIn1->w
mulps xmm3, xmm7;

addps xmm0, xmm1;
addps xmm2, xmm3;
addps xmm0, xmm2; 最终结果行保存在xmm0

movups [edx], xmm0; 将结果保存到pOut中
add edx, 16;
add eax, 16; 作为变址用

loop LOOPIT;
}
}
}

这里的汇编对照上面的一般算法就很容易理解。下回我会说说矩阵类的关键算法。