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试一下，我在这儿接受任何人的拍砖