毫无疑问,数学库是图形程序的基石,是图形程序运行效率的关键之一。一个优秀的数学库可以让图形程序运行得更流畅,甚至要快上几十倍上百倍。有时候替换一条除法运算会带来成倍的效率增长,比如用乘以 1/op 替换 vector 里的 operator /。当然,更高级的优化是使用 SIMD 优化海量运算,这就是本文的中心——SSE/SSE2 优化。
在描述 SSE/SSE2 优化前,我先介绍一般的 vector/matrix 库构造。当然,在 OpenEXR 里已经有一个非常优秀的 Imath 实现了,数学库的实现细节可以参照它。
SSE – Streaming SIMD Extension,是Intel从PIII开始加入的一种x86扩展指令集。在SSE以前,x86的浮点运算都是以栈式FPU完成的,有一定x86汇编经验的人应该不会对那些复杂的fld、fst指令陌生吧。而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,必须先确认你的编译器是否支持新的指令集。VC6 sp6、VC.net、.net 2003、ICL、GCC 、nasm 都支持SSE指令集。我推荐使用ICL,它的优化做得最棒,生成的指令最紧凑、效率最高。使用SSE有两种途径,一是直接编写汇编代码,但难度较大,需要有一定的汇编经验;二是使用SSE intrinsic,一种直接在C/C++里使用SSE指令的伪函数调用。在图形运算的核心环节上、如raytrace核心,我建议使用汇编,这样才能极大地体现出SSE的优势、与x86指令混合使用,并充分使用它的并行性。而在大多数场合下则推荐使用intrinsic,它的可读性高,而且编译器会在最后把函数调用替换成SSE指令,这样既不需要写内嵌汇编代码,又可以保证代码的执行效率。
#include <xmmintrin.h></xmmintrin.h><xmmintrin.h></xmmintrin.h>
里面定义了一个新的数据类型,__m128,这是一个128位、4个32位单精度浮点数的结构,如果你正在使用VC.net,你会看到它是一个关键字,被当作一种基本数据类型。要是你不打算使用汇编SSE,那么就没必要深究编译器在幕后到底如何处理__m128类型的数据,你只需要知道里面能存放四个float,而这四个float可以进行并行运算。
在定义了__m128后,文件声明一大堆对__m128进行运算的函数,如_mm_add_ps、_mm_sub_ps等等,这就是SSE运算指令的声明。使用SSE优化在这些声明的帮助下变得非常简单,如计算两个向量之和,平时需要每一个元素进行一次加法运算,现在只需要简单地:
__m128a , 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;
则可以直接写成:
__m128a = _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)上面。这是一条非常有用的指令,它可以把两个操作数的分量以特定的顺序排列并赋予给目标数。比如
__m128b = _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汇编的用法与优化。
感兴趣者可以到Chaos on Graphics阅读作者其他文章
相关推荐
Culture On The Edge Of Chaos,Cultural Algorithms and the Foundations of Social Intelligence, 2019.
Chaos.NaCl, Chaos.NaCl 加密库 Chaos.NaClChaos.NaCl 是用 C# 编写的加密库。 它是基于djb的NaCl 。目前它支持:Ed25519签名使用 Curve25519 ( 蒙哥马利形式) 或者 Ed25519 public 密钥的密钥交换
混沌测试chaosblade 2022年new version
chaos1.7.1混沌测试工具
chaosblade-operator镜像
Chaos是一个基于Linux平台, reactor模式的网络事件库, 目前仅支持TCP传输协议, 仅在x86_64下编译, 并遵循3-clause BSD开源协议. 在使用上, 可以说它很像boost asio, 可能是由于我对boost asio的接口设计很有爱...
chaosblade帮助文档
standish group chaos report 2020 2020年standish group chaos简报
资源分类:Python库 所属语言:Python 使用前提:需要解压 资源全名:chaos_cases-0.0.2-py3-none-any.whl 资源来源:官方 安装方法:https://lanzao.blog.csdn.net/article/details/101784059
分3个part,这个是part 2. [from amazon] Topics covered: Overview of fractals and chaos theory, feedback and multiple reduction copy machines (MRCMs), the Cantor Set, the Sierpinski Gasket and Carpet, ...
A color image cryptosystem based on dynamic DNA encryption and chaos
这是一个matlab的程序包,可以用于计算非线性问题时计算李雅普诺夫指数,用来判断运动的状态。
UE5 Chaos物理破碎使用教程
该项目是一个基于的chaosblade执行器,用于通过增强类在Java应用程序上进行混沌实验。 演习可以通过cli刀片进行,有关详细信息,请参见项目。 编译中 在项目根目录中,执行以下命令进行编译 make 编译结果将存储在...
chaosblade-1.6.1-linux-amd64.tar.gz
阿里云ChaosBlade项目分享,了解ChaosBlade运行原理以及阿里这方面想要打造的生态圈,ChaosBlade解决了混沌工程中故障注入部分的实现。
image encryption using chaos
Topics covered: Overview of fractals and chaos theory, feedback and multiple reduction copy machines (MRCMs), the Cantor Set, the Sierpinski Gasket and Carpet, the Pascal Triangle, the Koch Curve, ...
George M. Zaslavsky的《The physics of chaos in Hamiltonian systems》第二版
Wireless Communication with Chaos, 供无线通信相关人员参考