我读了一些关于使用 SSE 内在函数的文章，并尝试用双打实现四元数旋转。下面是我写的普通函数和 SSE 函数，
void quat_rot(quat_t a, REAL* restrict b){
///////////////////////////////////////////
// Multiply vector b by quaternion a //
///////////////////////////////////////////
REAL cross_temp[3],result[3];
cross_temp[0]=a.el[2]*b[2]-a.el[3]*b[1]+a.el[0]*b[0];
cross_temp[1]=a.el[3]*b[0]-a.el[1]*b[2]+a.el[0]*b[1];
cross_temp[2]=a.el[1]*b[1]-a.el[2]*b[0]+a.el[0]*b[2];
result[0]=b[0]+2.0*(a.el[2]*cross_temp[2]-a.el[3]*cross_temp[1]);
result[1]=b[1]+2.0*(a.el[3]*cross_temp[0]-a.el[1]*cross_temp[2]);
result[2]=b[2]+2.0*(a.el[1]*cross_temp[1]-a.el[2]*cross_temp[0]);
b[0]=result[0];
b[1]=result[1];
b[2]=result[2];
}
与上交所
inline void cross_p(__m128d *a, __m128d *b, __m128d *c){
const __m128d SIGN_NP = _mm_set_pd(0.0, -0.0);
__m128d l1 = _mm_mul_pd( _mm_unpacklo_pd(a[1], a[1]), b[0] );
__m128d l2 = _mm_mul_pd( _mm_unpacklo_pd(b[1], b[1]), a[0] );
__m128d m1 = _mm_sub_pd(l1, l2);
m1 = _mm_shuffle_pd(m1, m1, 1);
m1 = _mm_xor_pd(m1, SIGN_NP);
l1 = _mm_mul_pd( a[0], _mm_shuffle_pd(b[0], b[0], 1) );
__m128d m2 = _mm_sub_sd(l1, _mm_unpackhi_pd(l1, l1));
c[0] = m1;
c[1] = m2;
}
void quat_rotSSE(quat_t a, REAL* restrict b){
///////////////////////////////////////////
// Multiply vector b by quaternion a //
///////////////////////////////////////////
__m128d axb[2];
__m128d aa[2];
aa[0] = _mm_load_pd(a.el+1);
aa[1] = _mm_load_sd(a.el+3);
__m128d bb[2];
bb[0] = _mm_load_pd(b);
bb[1] = _mm_load_sd(b+2);
cross_p(aa, bb, axb);
__m128d w = _mm_set1_pd(a.el[0]);
axb[0] = _mm_add_pd(axb[0], _mm_mul_pd(w, bb[0]));
axb[1] = _mm_add_sd(axb[1], _mm_mul_sd(w, bb[1]));
cross_p(aa, axb, axb);
_mm_store_pd(b, _mm_add_pd(bb[0], _mm_add_pd(axb[0], axb[0])));
_mm_store_sd(b+2, _mm_add_pd(bb[1], _mm_add_sd(axb[1], axb[1])));
}
旋转基本上是使用函数完成的，
然后我运行以下测试来检查每个函数执行一组旋转所需的时间，
int main(int argc, char *argv[]){
REAL a[] __attribute__ ((aligned(16))) = {0.2, 1.3, 2.6};
quat_t q = {{0.1, 0.7, -0.3, -3.2}};
REAL sum = 0.0;
for(int i = 0; i < 4; i++) sum += q.el[i] * q.el[i];
sum = sqrt(sum);
for(int i = 0; i < 4; i++) q.el[i] /= sum;
int N = 1000000000;
for(int i = 0; i < N; i++){
quat_rotSSE(q, a);
}
printf("rot = ");
for(int i = 0; i < 3; i++) printf("%f, ", a[i]);
printf("\n");
return 0;
}
我使用 gcc 4.6.3 和 -O3 -std=c99 -msse3 编译。
使用 unix 的正常功能的时间time是 18.841s 和 SSE 的 21.689s。
我是否遗漏了什么，为什么 SSE 实施比正常实施慢 15%？在哪些情况下，双精度的 SSE 实现会更快？
编辑：从评论中得到建议，我尝试了几件事，
- -O1 选项给出非常相似的结果。
- 尝试使用
restrict函数cross_p并添加一个 __m128d 来保存第二个叉积。这在生产的组件中没有区别。 - 据我了解，为正常功能生成的程序集仅包含标量指令，除了一些
movapd.
为 SSE 函数生成的汇编代码仅比普通代码少 4 行。
编辑：添加到生成的程序集的链接，