在做图像处理的SSE优化时，也会经常遇到一些小的过程、数值优化等代码，本文分享一些个人收藏或实现的代码片段给大家。
一、快速求对数运算
对数运算在图像处理中也是个经常会遇到的过程，特备是在一些数据压缩和空间转换时常常会用到，而且是个比较耗时的函数，标准的SSE库里并没有提供该函数的实现，如果需要高精度的SSE版本，网络上已经有了，参考：https://github.com/to-miz/sse_mathfun_extension/blob/master/sse_mathfun.h，这个的精度和标准库的精度基本一致了，稍作整理后的代码如下：
// 对数函数的SSE实现，高精度版
inline __m128 _mm_log_ps(__m128 x)
{
)) ] = { 0x00800000, 0x00800000, 0x00800000, 0x00800000 };
)) ] = { ~0x7f800000, ~0x7f800000, ~0x7f800000, ~0x7f800000 };
)) ] = { 0x7f, 0x7f, 0x7f, 0x7f };
)) ] = { 1.0f, 1.0f, 1.0f, 1.0f };
)) ] = { 0.5f, 0.5f, 0.5f, 0.5f };
)) ] = { 0.707106781186547524f, 0.707106781186547524f, 0.707106781186547524f, 0.707106781186547524f };
)) ] = { 7.0376836292E-2f, 7.0376836292E-2f, 7.0376836292E-2f, 7.0376836292E-2f };
)) ] = { -1.1514610310E-1f, -1.1514610310E-1f, -1.1514610310E-1f, -1.1514610310E-1f };
)) ] = { 1.1676998740E-1f, 1.1676998740E-1f, 1.1676998740E-1f, 1.1676998740E-1f };
)) ] = { -1.2420140846E-1f, -1.2420140846E-1f, -1.2420140846E-1f, -1.2420140846E-1f };
)) ] = { 1.4249322787E-1f, 1.4249322787E-1f, 1.4249322787E-1f, 1.4249322787E-1f };
)) ] = { -1.6668057665E-1f, -1.6668057665E-1f, -1.6668057665E-1f, -1.6668057665E-1f };
)) ] = { 2.0000714765E-1f, 2.0000714765E-1f, 2.0000714765E-1f, 2.0000714765E-1f };
)) ] = { -2.4999993993E-1f, -2.4999993993E-1f, -2.4999993993E-1f, -2.4999993993E-1f };
)) ] = { 3.3333331174E-1f, 3.3333331174E-1f, 3.3333331174E-1f, 3.3333331174E-1f };
)) ] = { -2.12194440e-4f, -2.12194440e-4f, -2.12194440e-4f, -2.12194440e-4f };
)) ] = { 0.693359375f, 0.693359375f, 0.693359375f, 0.693359375f };
__m128 one = *(__m128*)_ps_1;
__m128 invalid_mask = _mm_cmple_ps(x, _mm_setzero_ps());
/* cut off denormalized stuff */
x = _mm_max_ps(x, *(__m128*)_ps_min_norm_pos);
__m128i emm0 = _mm_srli_epi32(_mm_castps_si128(x), );
/* keep only the fractional part */
x = _mm_and_ps(x, *(__m128*)_ps_inv_mant_mask);
x = _mm_or_ps(x, _mm_set1_ps(0.5f));
emm0 = _mm_sub_epi32(emm0, *(__m128i *)_pi32_0x7f);
__m128 e = _mm_cvtepi32_ps(emm0);
e = _mm_add_ps(e, one);
__m128 mask = _mm_cmplt_ps(x, *(__m128*)_ps_sqrthf);
__m128 tmp = _mm_and_ps(x, mask);
x = _mm_sub_ps(x, one);
e = _mm_sub_ps(e, _mm_and_ps(one, mask));
x = _mm_add_ps(x, tmp);
__m128 z = _mm_mul_ps(x, x);
__m128 y = *(__m128*)_ps_log_p0;
y = _mm_mul_ps(y, x);
y = _mm_add_ps(y, *(__m128*)_ps_log_p1);
y = _mm_mul_ps(y, x);
y = _mm_add_ps(y, *(__m128*)_ps_log_p2);
y = _mm_mul_ps(y, x);
y = _mm_add_ps(y, *(__m128*)_ps_log_p3);
y = _mm_mul_ps(y, x);
y = _mm_add_ps(y, *(__m128*)_ps_log_p4);
y = _mm_mul_ps(y, x);
y = _mm_add_ps(y, *(__m128*)_ps_log_p5);
y = _mm_mul_ps(y, x);
y = _mm_add_ps(y, *(__m128*)_ps_log_p6);
y = _mm_mul_ps(y, x);
y = _mm_add_ps(y, *(__m128*)_ps_log_p7);
y = _mm_mul_ps(y, x);
y = _mm_add_ps(y, *(__m128*)_ps_log_p8);
y = _mm_mul_ps(y, x);
y = _mm_mul_ps(y, z);
tmp = _mm_mul_ps(e, *(__m128*)_ps_log_q1);
y = _mm_add_ps(y, tmp);
tmp = _mm_mul_ps(z, *(__m128*)_ps_0p5);
y = _mm_sub_ps(y, tmp);
tmp = _mm_mul_ps(e, *(__m128*)_ps_log_q2);
x = _mm_add_ps(x, y);
x = _mm_add_ps(x, tmp);
x = _mm_or_ps(x, invalid_mask); // negative arg will be NAN
return x;
}
看上去有一大堆代码，不过实测这个的速度越是标准库（本文是指启动增强指令集选项设置为：未设置，设计上编译器在此种情况下会自动设置为SSE2增强，这可以从反编译logf函数看到，因此，这里的速度比较还不是和纯Fpu实现的比较）的2倍，如果稍微降低点精度，比如_ps_log_p5到_ps_log_p8之间的代码，还能提高点速度。
另外，在很多场合我们还可以使用另外一种低精度的log函数，其C代码如下所示：
//https://*.com/questions/9411823/fast-log2float-x-implementation-c
inline float IM_Flog(float val)
{
union
{
float val;
int x;
} u = { val };
) & ) - );
u.x &= ~( << );
u.x += ( << );
log_2 += ((-0.34484843f) * u.val + 2.02466578f) * u.val - 0.67487759f;
return log_2 * 0.69314718f;
}
这个函数大概有小数点后2位精度。
上述代码大约也是标准函数的2倍速度左右。但是上述函数是可以向量化的，我们来尝试实现。
我们首先来看联合体，其实这个东西就是两个东西占同一个内存空间，然后外部用不同的规则去读取他，在SSE里，有着丰富的cast函数，他也是干这个事情的，比如这里的联合体就可以用_mm_castps_si128来转换，而实际上这个Intrinsic并不会产生任何的汇编语句。
那么后面的那些移位、或运算、非运算、加减乘除之类的就是直接翻译了，毫无难处，完整的代码如下所示：
inline __m128 _mm_flog_ps(__m128 x)
{
__m128i I = _mm_castps_si128(x);
__m128 log_2 = _mm_cvtepi32_ps(_mm_sub_epi32(_mm_and_si128(_mm_srli_epi32(I, ), _mm_set1_epi32()), _mm_set1_epi32()));
I = _mm_and_si128(I, _mm_set1_epi32(-)); // 255 << 23
I = _mm_add_epi32(I, _mm_set1_epi32()); // 127 << 23
__m128 F = _mm_castsi128_ps(I);
__m128 T = _mm_add_ps(_mm_mul_ps(_mm_set1_ps(-0.34484843f), F), _mm_set1_ps(2.02466578f));
T = _mm_sub_ps(_mm_mul_ps(T, F), _mm_set1_ps(0.67487759f));
return _mm_mul_ps(_mm_add_ps(log_2, T), _mm_set1_ps(0.69314718f));
}
经过实测，这个速度可以达到标准库的7到8倍的优势。
二、快速求幂运算
一般图像编程中有log出现的地方就会有exp出现，因此exp的优化也尤为重要，同样在sse_mathfun.h中也有exp的优化（还有sin,cos的SSE优化语句呢），我这里就不贴那个的代码了，我们同样关注下用联合体实现的近似快速算法，其C代码如下所示：
inline float IM_Fexp(float Y)
{
union
{
double Value;
];
} V;
V.X[] = ( + + 0.5F);
V.X[] = ;
return (float)V.Value;
}
测试这个和标准的exp库函数速度居然差不多，不晓得为啥，但我们来试下他的SSE优化版本了。
V.X[1] = (int)(Y * 1512775 + 1072632447 + 0.5F);这句话没啥难度，直接翻译就可以了，注意几个强制类型转化就可以了，如下所示：
__m128i T = _mm_cvtps_epi32(_mm_add_ps(_mm_mul_ps(Y, _mm_set1_ps()), _mm_set1_ps()));
由于我们想一次性处理4个float类型的数据，因此也就需要4个union的空间，这样就需要2个__m128i变量来保存数据，每个XMM寄存器的数据应该分别为：
T1 0 T0 0 + T3 0 T2 0 （高位----》低位）
这个可以使用unpack来实现，具体如下：
__m128i TL = _mm_unpacklo_epi32(_mm_setzero_si128(), T);
__m128i TH = _mm_unpackhi_epi32(_mm_setzero_si128(), T);
最后我们认为__m128i里的数据是double数据，直接一个cast就可以了，然后因为我们只需要单精度的数据，再使用_mm_cvtpd_ps将double转换为float类型，注意这个时候还需要将他们连接再一起形成一个完整的__m128变量，最终的代码如下：
inline __m128 _mm_fexp_ps(__m128 Y)
{
__m128i T = _mm_cvtps_epi32(_mm_add_ps(_mm_mul_ps(Y, _mm_set1_ps()), _mm_set1_ps()));
__m128i TL = _mm_unpacklo_epi32(_mm_setzero_si128(), T);
__m128i TH = _mm_unpackhi_epi32(_mm_setzero_si128(), T);
return _mm_movelh_ps(_mm_cvtpd_ps(_mm_castsi128_pd(TL)), _mm_cvtpd_ps(_mm_castsi128_pd(TH)));
}
实测这个的提速大概有10倍。
如果要求double的exp，其SSE代码你会了吗？
三、pow函数的优化。
一种常用的近似算法如下所示：
inline float IM_Fpow(float a, float b)
{
union
{
double Value;
];
} V;
V.X[] = (] - ) + );
V.X[] = ;
return (float)V.Value;
}
和exp很类似，留给有兴趣的人自己实现。
四：两个求倒数函数的优化误区
SSE提供了连个快速求倒数的函数，_mm_rcp_ps，_mm_rsqrt_ps，他们都是近似值，只有12bit的精度，如果想通过他们得到精确的倒数值，需要牛顿 - 拉弗森方法，比如利用_mm_rcp_ps求精确倒数的代码如下：
__forceinline __m128 _mm_prcp_ps(__m128 a)
{
__m128 rcp = _mm_rcp_ps(a); // 此函数只有12bit的精度.
return _mm_sub_ps(_mm_add_ps(rcp, rcp), _mm_mul_ps(a, _mm_mul_ps(rcp, rcp))); // x1 = x0 * (2 - d * x0) = 2 * x0 - d * x0 * x0，使用牛顿 - 拉弗森方法这种方法可以提高精度到23bit
}
但是实测这个还不如直接用_mm_div_ps的速度，即使是下面的函数：
__forceinline __m128 _mm_fdiv_ps(__m128 a, __m128 b)
{
return _mm_mul_ps(a, _mm_rcp_ps(b));
}
似乎速度也不够好，而且精度还低了。
特别低，如果使用_mm_rcp_ps和_mm_rsqrt_ps联合求近似sqrt，即如下代码，速度好像还慢了，真搞不明白为什么。
__forceinline __m128 _mm_fsqrt_ps(__m128 a)
{
return _mm_rcp_ps(_mm_rsqrt_ps(a));
}
五、其他参考
在http://www.alfredklomp.com/programming/sse-intrinsics/以及 http://www.itkeyword.com/doc/0326039046115117x827/c++-sse2-intrinsics-comparing-unsigned-integers等网站上还有很多参考的资料，希望大家自己去学习下。