高效率的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) || ((y-v.y)-FLOAT_EPS) || ((z-v.z)-FLOAT_EPS))? false:true); } // 运算符重载 inline bool Vector::operator !=(Vector& v) { return ((((x-v.x)-FLOAT_EPS) || ((y-v.y)-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和浮点指令的理解，以及最有用的矩阵相乘算法。