Raytrace world

Chaos 的光线跟踪世界

  C++博客 :: 首页 :: 新随笔 :: 联系 :: 聚合  :: 管理 ::
  4 随笔 :: 0 文章 :: 63 评论 :: 0 Trackbacks

Chaos on Graphics

About 数学库之一 SSE/SSE2

 

Chaos Chiao

 

 

毫无疑问,数学库是图形程序的基石,是图形程序运行效率的关键之一。一个优秀的数学库可以让图形程序运行得更流畅,甚至要快上几十倍上百倍。有时候替换一条除法运算会带来成倍的效率增长,比如用乘以 1/op 替换 vector 里的 operator /。当然,更高级的优化是使用 SIMD 优化海量运算,这就是本文的中心——SSE/SSE2 优化。

在描述 SSE/SSE2 优化前,我先介绍一般的 vector/matrix 库构造。当然,在 OpenEXR 里已经有一个非常优秀的 Imath 实现了,数学库的实现细节可以参照它。

在图形程序里我们经常会遇到向量运算,这是标准C++编译器所不能直接支持的,如三维空间向量。传统的C图形程序会使用“数组+宏”的实现方式:

 

typedef float vector[3];

  

到了C++时代,一般会封装成:

 

class Vector

{

private:

    float x , y , z;

};

 

然后加入通常的各种method,如

 

const float& X( void ) const {

    return x;

}

 

标准的向量算法如内积、外积、单位化、长度、运算等等,都可以封装为成员函数。

 

Vector operator + ( const Vector& a , const float & b ) {

    return Vector( a.x + b , a.y + b , a.z + b );

}

 

类似的数学库可以在Aqsis等一些开源的图形程序里找到。不过这些结构并不适合接下来我们要讨论的SSE/SSE2优化。

 

SSE – Streaming SIMD Extension,是IntelPIII开始加入的一种x86扩展指令集。在SSE以前,x86的浮点运算都是以栈式FPU完成的,有一定x86汇编经验的人应该不会对那些复杂的fldfst指令陌生吧。而SSE一方面让浮点运算可以像整数运算的模式、如 add eax , ebx 那样通过直接访问寄存器完成,绕开了讨厌的栈,另一方面引入了SIMD这个概念。SIMD – Single Instruction Multiply Data,顾名思义,它可以同时让一条指令在多个数据上执行,这种体系结构在一度在大型机上非常流行,需要经常进行海量运算的大型机器通常会通过一个数学SIMD虚拟机加快处理速度,比如同时让一组数据执行一个变换,数据的规模有上百万之巨,而SIMD则可以优化数据的存储与运算,减免某些切换Context的开销。

在硬件层面上面支持SIMD,某程度是因游戏需要的驱使,因为越来越多的3D游戏涉及大量的向量操作,一般的浮点运算优化已经不能再适应这种并行运算的需要了,而直接从指令上支持SIMD操作则可以进一步简化向量运算的优化,提高指令执行效率。像 addps 这样的SSE指令,可以并行执行四个32位浮点数的加法运算,而延迟只有4 cycle;相比之下,原来的fadd指令光执行一个32位单精度浮点数加法的延迟已经达到了3 cycle了,还没计算fst等存储指令的延迟。(具体见后面的指令执行单元表)

显然,SSE能给图形程序带来极大的优化,其提高远胜于基于整数的MMX与双单元单精度浮点数的3DNow!。但SSE对数据组织的要求是苛刻的,若要发挥SSE的最大威力,我们还需要进行对齐向量数据,把向量对齐到16字节。如果我们正在使用一般的三分量向量,那么就意味着有要浪费四分之一的存储空间来换取速度。当然,这4字节还可以有很多用途,只是你必须处理得非常小心,因为任何运算都将同时应用到四个分量上。

要使用SSE,必须先确认你的编译器是否支持新的指令集。VC6 sp6VC.net.net 2003ICLGCC nasm 都支持SSE指令集。我推荐使用ICL,它的优化做得最棒,生成的指令最紧凑、效率最高。使用SSE有两种途径,一是直接编写汇编代码,但难度较大,需要有一定的汇编经验;二是使用SSE intrinsic,一种直接在C/C++里使用SSE指令的伪函数调用。在图形运算的核心环节上、如raytrace核心,我建议使用汇编,这样才能极大地体现出SSE的优势、与x86指令混合使用,并充分使用它的并行性。而在大多数场合下则推荐使用intrinsic,它的可读性高,而且编译器会在最后把函数调用替换成SSE指令,这样既不需要写内嵌汇编代码,又可以保证代码的执行效率。

 

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

 

#include

 

里面定义了一个新的数据类型,__m128,这是一个128位、4个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逼近法求出较精确值,下面的代码摘自NVfastmath.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( __m128 a , __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( __m128 a )

{

    / 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 __m128 a ) {

    // get length square

    __m128 l = _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( __m128 a )

{

    static const union { int i[4]; __m128 m; } __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有两兄弟,分别是locuthicut

 

__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对数学运算进行优化,并尽可能地不分拆向量,这样可以保证8128位的xmm寄存器可以满足大部分运算。不过SSE intrinsic始终受编译器生成代码的质量好坏影响,没能真正发挥出SSE的全部威力,接下来我们将讨论SSE汇编的用法与优化。

 


to be continued...

posted on 2005-10-22 02:39 Chaos Chiao 阅读(4554) 评论(7)  编辑 收藏 引用 所属分类: Chaos on Graphics

评论

# re: Chaos on Graphics ~ About 数学库之一 SSE/SSE2 2006-02-05 16:21 yyy
看了5遍,越看越好看!  回复  更多评论
  

# re: Chaos on Graphics ~ About 数学库之一 SSE/SSE2 2006-04-25 03:22 信仔
请告之 现在做什么??

佩服!  回复  更多评论
  

# re: Chaos on Graphics ~ About 数学库之一 SSE/SSE2 2006-04-25 19:39 丁勇
希望与博客主人及各位交流,科技创新,http://dycomputer.superbuckler.com/,或者http://blog.sina.com.cn/u/1224675400  回复  更多评论
  

# re: Chaos on Graphics ~ About 数学库之一 SSE/SSE2 2006-04-25 19:44 晴天
2030就是欧美软件外包专才培养计划,http://www.i2030.org/  回复  更多评论
  

# re: Chaos on Graphics ~ About 数学库之一 SSE/SSE2 2007-02-11 13:18 pb
写的很好。谢谢分享。
不知道这个系列还有没:)  回复  更多评论
  

# re: Chaos on Graphics ~ About 数学库之一 SSE/SSE2 2008-04-04 20:19 Suzuki
i think there is practical significances only in 3D program. as for general program, vector is seldom used. nonetherless it is sound for boast. i wonder if there is hardware multiplier.  回复  更多评论
  

# re: Chaos on Graphics ~ About 数学库之一 SSE/SSE2 2008-09-28 17:14 luguo
那么使用AMD的是否也有类似的指令来优化
使用这些代码在AMD上应该不能运行才对..  回复  更多评论
  


只有注册用户登录后才能发表评论。
【推荐】超50万行VC++源码: 大型组态工控、电力仿真CAD与GIS源码库
网站导航: 博客园   IT新闻   BlogJava   知识库   博问   管理