数学库之二中应用：SSE/SSE2

引自：http://hi.baidu.com/sige_online/blog/item/b48db235ef12171b91ef39ad.html

下面将通过几个简单的运算例子介绍 SSE intrinsic的使用。首先，使用SSE需要一个新的头文件

#include

 

 

 

 

里面定义了一个新的数据类型，__m128 ，这是一个 128位、４个32位单精度浮点数的结构，如果你正在使用VC.net，你会看到它是一个关键字，被当作一种基本数据类型。要是你不打算使用汇编SSE，那么就没必要深究编译器在幕后到底如何处理__m128 类型的数据，你只需要知道里面能存放四个float ，而这四个float 可以进行并行运算。

在定义了__m128 后，文件声明一大堆对__m128 进行运算的函数，如_mm_add_ps、_mm_sub_ps等等，这就是 SSE运算指令的声明。使用SSE优化在这些声明的帮助下变得非常简单，如计算两个向量之和，平时需要每一个元素进行一次加法运算，现在只需要简单地：

__m128 a , b , c ;

 

c = _mm_add_ps ( a , b );

 

这样等价于：

float a [ 4] , b [ 4] , c [ 4];

 

for ( int i = 0 ; i < 4 ; ++ i )

 

    c [ i ] = a [ i ] + b [ i ];

 

但前者的运算是并行的，在一般情况下效率远比后者要高。况且在描述复杂的运算的时候，如：

      a = b * c + d / e ;

则可以直接写成：

__m128 a = _mm_add_ps ( _mm_mul_ps ( b , c ) , _mm_div_ps ( d , e ) );

 

咋看之下，很多效率至上的人马上就会大叫“昂贵的函数调用啊！ Bad smell code!”。其实我正要告诉你，我也是效率至上派的。前面已经说过了，这些看上去貌似函数的调用实际上并非函数，而是所谓intrinsic，它们在编译优化中将被解释为单条或多条SSE指令，而且编译器会自动调节调用顺序以使其最大并行效率。

不过除了直接使用这些 intrinsic以外，我们还可以把它们封装到类里面，重载运算符，这样就可以把运算写成可读性更强的算术式。如果你不愿意自己动手封装，也可以使用Intel封装好了的F32vec4类，它提供了完备的运算符重载，完全使用SSE，非常方便。

虽然 Intel封装好的类已经很完善了，但还有一大堆数学运算需要我们自己动手进行编写，如内积（点积）和外积（叉积）。

首先来看一个比较实用的运算，求倒数。求倒数在很多数学库里都有专门的优化，通常原理都是先求出一个近似值，然后通过 Newton-Raphson逼近法求出较精确值，下面的代码摘自NV的fastmath.cpp:

#define FP_ONE_BITS 0x3F800000

 

// r = 1/p

 

#define FP_INV ( r , p )                                                    /

 

{                                                                           /

 

    int _i = 2 * FP_ONE_BITS - *( int *)&( p );                   /

 

    r = *( float *)& _i ;                                                 /

 

    r = r * ( 2.0f - ( p ) * r );                                        /

 

}

 

而在 SSE里也提供了两条求倒数的指令rcpss/rcpps（对应的intrinsic是_mm_rcp_ss与_mm_rcp_ps），不过这两条指令求的并非是精确值，而是近似值，所以我们需要对它的结果进行逼近处理。

float __rcp < float >( const float & a ) {

 

    register float r ;

 

    __m128 rcp = _mm_load_ss ( & a );

   rcp = _mm_rcp_ss ( rcp );

    _mm_store_ss ( & r , rcp );

    /* [2 * rcpps(x) - (x * rcpps(x) * rcpps(x))] */

    r = 2.0f * r - ( a * r * r );

    return r ;

}

原理一致，只不过我们还可以用_mm_rcp_ps并行求四分量的倒数。如果你还对SSE的威力有所保留，那我建议你设计一个测试单元测试一下使用除法求倒数与使用SSE求倒数，看效率到底是谁更高、高多少。当然，我自己已经测试过很多次了J。

然后我们把注意力放到一条非常特殊的指令shufps（对应intrinsic是_mm_shuffle_ps）上面。这是一条非常有用的指令，它可以把两个操作数的分量以特定的顺序排列并赋予给目标数。比如

__m128 b = _mm_shuffle_ps ( a , a , 0 );

　　

则 b 的所有分量都是 a 中下标为0的分量。第三个参数控制分量分配，是一个8bit的常量，这个常量的1~8位分别控制了从两个操作数中选择分量的情况，具体怎么控制将在后面讨论SSE汇编中一并说明，而在使用intrinsic的时候，最好使用_MM_SHUFFLE 宏，它可以定义分配情况。下面我们来复习一下叉积的求法。

c = a x b

可以写成：

Vector cross ( const Vector & a , const Vector & b ) {

    return Vector (

        ( a [ 1] * b [ 2] - a [ 2] * b [ 1] ) ,

        ( a [ 2] * b [ 0] - a [ 0] * b [ 2] ) ,

        ( a [ 0] * b [ 1] - a [ 1] * b [ 0] ) );

}

那么写成SSE intrinsic形式则是：

/* cross */

__m128 _mm_cross_ps ( __m128a , __m128 b ) {

    __m128 ea , eb ;

    // set to a[1][2][0][3] , b[2][0][1][3]

    ea = _mm_shuffle_ps ( a , a , _MM_SHUFFLE ( 3 , 0 , 2 , 1 ) );

    eb = _mm_shuffle_ps ( b , b , _MM_SHUFFLE ( 3 , 1 , 0 , 2 ) );

    // multiply

    __m128 xa = _mm_mul_ps ( ea , eb );

    // set to a[2][0][1][3] , b[1][2][0][3]

    a = _mm_shuffle_ps ( a , a , _MM_SHUFFLE ( 3 , 1 , 0 , 2 ) );

    b = _mm_shuffle_ps ( b , b , _MM_SHUFFLE ( 3 , 0 , 2 , 1 ) );

    // multiply

    __m128 xb = _mm_mul_ps ( a , b );

    // subtract

    return _mm_sub_ps ( xa , xb );

}

这就是shuffle强大的地方，它可以直接在寄存器里直接调整分量的顺序。而且配合_mm_movehl_ps，我们可以轻松解决点积的运算。_mm_movehl_ps把操作数高位两个分量赋予目标数的低位两分量，而目标数的高位两分量值不变，相当于：

a [ 0] = b [ 2];

a [ 1] = b [ 3];

三分量的向量求点积，可以写成：

float dot ( const float & a , const float & b ) const {

    return a [ 0] * b [ 0] + a [ 1] * b [ 1] + a [ 2] * b [ 2];

}

则用SSE intrinsic可以写成：

/* x[0] * x[1] + y[0] * y[1] + z[0] * z[1] */

__m128 _mm_dot_ps ( __m128 x , __m128 y ) {

    __m128 s , r ;

    s = _mm_mul_ps ( x , y );

    r = _mm_add_ss ( s , _mm_movehl_ps ( s , s ) );

    r = _mm_add_ss ( r , _mm_shuffle_ps ( r , r , 1 ) );

    return r ;

}

通过这两个例子，可以留意到向量内元素的垂直相加一般形式，即：

/* x[0] + x[1] + x[2] + x[3] */

__m128 _mm_sum_ps ( __m128 x ) {

    __m128 r ;

    r = _mm_add_ps ( x , _mm_movehl_ps ( x , x ) );

    r = _mm_add_ss ( r , _mm_shuffle_ps ( r , r , 1 ) );

    return r ;

}

那么通过扩展，可以得到求向量长度的函数，首先是求分量平方和函数：

/* x[0] * x[0] + y[0] * y[0] + z[0] * z[0] */

__m128 _mm_square_ps ( __m128 x ) {

    __m128 s , r ;

    s = _mm_mul_ps ( x , x );

    r = _mm_add_ss ( s , _mm_movehl_ps ( s , s ) );

    r = _mm_add_ss ( r , _mm_shuffle_ps ( r , r , 1 ) );

    return r ;

}

然后就可以直接把结果求平方根，可得长度。解决了长度，接下来则是很重要的单位化了。可以说单位化是最重要的一个函数，它经常被调用到，而函数内的陷阱却又最多。求单位化其实并不难，就是分量除以向量长度，可以写成：

void normalize ( const Vector & a ) {

    float len = a [ 0] * a [ 0] + a [ 1] * a [ 1] + a [ 2] * a [ 2];

    if ( is_zero ( len ) )

        return ;

    len = 1 / len ;

    a [ 0] *= len ;

    a [ 1] *= len ;

    a [ 2] *= len ;

}

我和这个家伙打交道已经有差不多七年时间了，所以脾性非常熟悉。首先求分量的平方和，判断是否为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寄存器可以满足大部分运算。不过SSE intrinsic始终受编译器生成代码的质量好坏影响，没能真正发挥出SSE的全部威力，接下来我们将讨论SSE汇编的用法与优化。