高效率的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, xmm0movaps xmm3, xmm1shufps xmm0, xmm0, 0xc9shufps xmm1, xmm1, 0xd2mulps xmm0, xmm1shufps xmm2, xmm2, 0xd2shufps xmm3, xmm3, 0xc9mulps xmm2, xmm3subps xmm0, xmm2mov eax, thismovups [eax], xmm0mov [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 _DEBUGassert(w!=1.0f && w!=0.0f);#endifif (g_bUseSSE2){_asm{mov ecx, this;mov edx, m;movss xmm0, [ecx];//lea eax, vr;shufps xmm0, xmm0, 0; // xmm0 = x,x,x,xmovss xmm1, [ecx+4];mulps xmm0, [edx];shufps xmm1, xmm1, 0; // xmm1 = y,y,y,ymovss xmm2, [ecx+8];mulps xmm1, [edx+16];shufps xmm2, xmm2, 0; // xmm2 = z,z,z,zmovss xmm3, [ecx+12];mulps xmm2, [edx+32];shufps xmm3, xmm3, 0; // xmm3 = w,w,w,waddps xmm0, xmm1;mulps xmm3, [edx+48];addps xmm0, xmm2;addps xmm0, xmm3; // xmm0 = resultmovups [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，*= 和AngleWithGetLengthSq有潜在危险，因为我是根据.Net2003的编辑器来写的代码，我知道ecx==this，知道float的返回值是直接从浮点栈寄存器fstp到外面参数的，所以，我会用这种方法来写，甚至没有写返回值！而看此文的您可能不会使用与我一样的编辑器，所以，在理解了实质之后，运用合理的算法来实现你的数学库。后面的函数都使用了编辑器无关的方法写的。*= 的运算符重载的潜在危险在于，Vector是4D的，可以表示3D的向量或者3D空间点坐标。如果是向量，则w==0，这样就只会受到旋转和缩放的影响。而如果是表示空间点，w==1，就会受到所有类型的变动，如平移、旋转和缩放。由于向量是不能平移的，处于对运算效率的考虑，这时候就需要数学库的调用者自己注意了。AngleWith函数之所以不对其进行内联化，是因为在以后的文章中，我会去进一步优化这里的代码。GetLength和acosf都不是内联函数，我必须要将其展开，以汇编实现，并重新组织编码。这个函数好像在D3DX9的数学库中是没有的~~没办法比较了。以上几个函数的效率与D3DX库比较结果大致是这样的：GetLengthSq微高于D3DXGetLength是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和浮点指令的理解，以及最有用的矩阵相乘算法。To be Continued.....