SSE优化在数学库中的应用之四
引自：http://hi.baidu.com/sige_online/blog/item/d8fdfffc8f0033f7fd037fac.html
我和这个家伙打交道已经有差不多七年时间了，所以脾性非常熟悉。首先求分量的平方和，判断是否为0（问我为什么不直接用 if( len == 0 )？好样的，请先去复习一下浮点数的基本知识），然后再求倒数，最后反映到分量上。在把它写成SSE intrinsic格式前，我先引入另外一个能极大提升运算效率的函数，求平方根的倒数。有数值运算编成经验的人都知道，如果说除法是恶魔的话，那么平方根就是撒旦了，而平方根的倒数简直就是撒旦他妈。虽然上面提供了倒数的逼近方法，但仅仅使用它还是绕不开最主要的开销、平方根运算。幸好，SSE提供了一个直接计算平方根倒数近似值的指令，rsqrtss/rsqrtps（即_mm_rsqrt_ss和_mm_rsqrt_ps）。照搬倒数求法，可以轻松得出：/* r = 1 / sqrt(a) */
/* 0.5 * rsqrtss * (3 - x * rsqrtss(x) * rsqrtss(x)) */
__m128 _mm_rsqrt( __m128a )
{
/／ divisor
static const __m128 _05 = _mm_set1_ps( 0.5f );
static const __m128 _3 = _mm_set1_ps( 3.f );
__m128 rsqrt = _mm_rsqrt_ss( a );
rsqrt =
_mm_mul_ss(
_mm_mul_ss( _05 , rsqrt ) ,
_mm_sub_ss( _3 , _mm_mul_ss( a , _mm_mul_ss( rsqrt , rsqrt ) ) ) );
return rsqrt;
}
那么就可以轻松得出单位化向量的函数了：
// normalize & return value
__m128 _mm_normalize( const __m128a ) {
// get length square
__m128l = _mm_square_ps( a );
// test if length is zero
if( _mm_iszero_ss( l ) )
return z;
// length inverse
l = _mm_rsqrt( l );
// shuffle
l = _mm_shuffle_ps( l , l , 0 );
// multiply to vector
return _mm_mul_ps( a , l );
}
SSE除了以上这些数学运算操作外，还提供了位运算。位运算？想到什么了吗？对！比较与选择。首先来看一个最简单的，求绝对值。通常我们会把 abs 写成非常简洁的形式：
float abs( float a ) {
a >= 0 ? a : -a;
}
但当我们已经Pack了一个向量到__m128结构里，而又不希望Unpack他们进行浮点数的比较，那么就可以使用SSE的位操作。
/* abs */
__m128 _mm_abs_ps( __m128a )
{
static const union { int i[4]; __m128m; } __mm_abs_mask_cheat_ps
= {0x7fffffff, 0x7fffffff, 0x7fffffff, 0x7fffffff};
return _mm_and_ps( a, __mm_abs_mask_cheat_ps.m );
}
还记得单精度浮点数的符号存放在什么位上面吗？我们只需把它除掉，然后就可以很轻松地得到了正值了。
图形程序很多时候会用到32位浮点色彩，其值域通常为[0,1]，所以clamp函数出现的频率也十分频繁。要将rgba的值同时clamp到值域内，毫无疑问，SSE的并行特性又得到了发挥的机会。先来看cut函数，它负责把超出值域的值干掉，但为了更灵活，我们一次只cut一边的区间，所以cut有两兄弟，分别是locut和hicut。
__m128 _mm_locut_ps( __m128 val , __m128 bound )
{
__m128 mask = _mm_cmplt_ps( val , bound );
return _mm_or_ps( _mm_and_ps( mask , bound ) , _mm_andnot_ps( mask , val ) );
}
__m128 _mm_hicut_ps( __m128 val , __m128 bound )
{
__m128 mask = _mm_cmpgt_ps( val , bound );
return _mm_or_ps( _mm_and_ps( mask , bound ) , _mm_andnot_ps( mask , val ) );
}
_mm_cmp**_ps是一系列的比较函数，非常丰富，也很好用，如果替换成相应的if-else，并行性将被大大削弱。不过_mm_cmp**_ps的最大缺点就是灵活度不够，返回值是一系列位标记，其具体的用法将在SSE汇编中讨论。有了这俩哥们，clamp变得非常简单：
__m128 _mm_clamp_ps( __m128 val , __m128 min , __m128 max )
{
return _mm_locut_ps( _mm_hicut_ps( val , max ) , min );
}
以上只是一些很简单的实现，利用了SSE intrinsic对数学运算进行优化，并尽可能地不分拆向量，这样可以保证8个128位的xmm寄存器可以满足大部分运算。
