使用AVX的四元数平方

Question

是否有人知道如何使用AVX对这个功能进行向量化

void cuadradoYSumaNormal(quaternion* a, quaternion* b, quaternion* c) {
          c->w = a->w*a->w - a->x*a->x - a->y*a->y - a->z*a->z + b->w;
          c->x = 2.*a->w*a->x + b->x;
          c->y = 2.*a->w*a->y + b->y;
          c->z = 2.*a->w*a->z + b->z;
    }

我可以假设a，b和c的单位长度

quaternion是以下结构：

struct quaternion{
  double w;
  double x;
  double y;
  double z;
};

该函数必须做的是将四元数*a平方（使用四元数乘法规则），然后添加四元数*b并将结果存储在*c中。

Answer 1

a具有单位长度，即aw^2+ax^2+ay^2+az^2 == 1时，此解决方案有效>

在那种情况下，c->w的计算等同于计算2*a->w*a->w - 1.0 + b->w，这使得向量化变得更加容易。可以通过将a（或a->w）加到自身上来实现2的乘法。为了减少延迟链，应将-1.0添加到b->w。可能的实现：

inline __m256d unit(double value = 1.0)
{
    return _mm256_set_pd(0,0,0,value);
}

void cuadradoYSumaNormal_avx(quaternion* a, quaternion* b, quaternion* c) {

    __m256d aw = _mm256_broadcast_sd(&a->w);
    __m256d a_ = _mm256_loadu_pd(&a->w);
    __m256d b_ = _mm256_loadu_pd(&b->w);

    __m256d a_squared_plus_one = _mm256_mul_pd(aw, _mm256_add_pd(a_,a_));
    __m256d c_ = _mm256_add_pd(a_squared_plus_one, _mm256_add_pd(b_, unit(-1.0)));

    _mm256_storeu_pd(&c->w, c_);
}
如果除了AVX之外，您还可以使用FMA，则可以加入一些加法和乘法运算到

(aw * a + [-0.5,0,0,0]) * 2.0 + b
仅产生两个FMA（一个广播和一些负载）。可能的实现：

void cuadradoYSumaNormal_fma(quaternion* a, quaternion* b, quaternion* c) {

    __m256d aw = _mm256_broadcast_sd(&a->w);
    __m256d a_ = _mm256_loadu_pd(&a->w);
    __m256d b_ = _mm256_loadu_pd(&b->w);

    __m256d a_squared_half = _mm256_fmadd_pd(aw, a_, unit(-0.5));
    __m256d c_ = _mm256_fmadd_pd(a_squared_half, _mm256_set1_pd(2.0), b_);

    _mm256_storeu_pd(&c->w, c_);
}