如何用最少的指令将两个四元数相乘？

Posted 2023-04-12

【中文标题】如何用最少的指令将两个四元数相乘？

【英文标题】：How to multiply two quaternions with minimal instructions?

【发布时间】：2013-09-03 18:32:47

【问题描述】：

经过一番思考，我想出了以下代码，用于使用 SSE 将两个四元数相乘：

#include <pmmintrin.h> /* SSE3 intrinsics */

/* multiplication of two quaternions (x, y, z, w) x (a, b, c, d) */

__m128 _mm_cross4_ps(__m128 xyzw, __m128 abcd)

    /* The product of two quaternions is:                                 */
    /* (X,Y,Z,W) = (xd+yc-zb+wa, -xc+yd+za+wb, xb-ya+zd+wc, -xa-yb-zc+wd) */

    __m128 wzyx = _mm_shuffle_ps(xyzw, xyzw, _MM_SHUFFLE(0,1,2,3));
    __m128 baba = _mm_shuffle_ps(abcd, abcd, _MM_SHUFFLE(0,1,0,1));
    __m128 dcdc = _mm_shuffle_ps(abcd, abcd, _MM_SHUFFLE(2,3,2,3));

    /* variable names below are for parts of componens of result (X,Y,Z,W) */
    /* nX stands for -X and similarly for the other components             */

    /* znxwy  = (xb - ya, zb - wa, wd - zc, yd - xc) */
    __m128 ZnXWY = _mm_hsub_ps(_mm_mul_ps(xyzw, baba), _mm_mul_ps(wzyx, dcdc));

    /* xzynw  = (xd + yc, zd + wc, wb + za, yb + xa) */
    __m128 XZYnW = _mm_hadd_ps(_mm_mul_ps(xyzw, dcdc), _mm_mul_ps(wzyx, baba));

    /* _mm_shuffle_ps(XZYnW, ZnXWY, _MM_SHUFFLE(3,2,1,0)) */
    /*      = (xd + yc, zd + wc, wd - zc, yd - xc)        */
    /* _mm_shuffle_ps(ZnXWY, XZYnW, _MM_SHUFFLE(2,3,0,1)) */
    /*      = (zb - wa, xb - ya, yb + xa, wb + za)        */

    /* _mm_addsub_ps adds elements 1 and 3 and subtracts elements 0 and 2, so we get: */
    /* _mm_addsub_ps(*, *) = (xd+yc-zb+wa, xb-ya+zd+wc, wd-zc+yb+xa, yd-xc+wb+za)     */

    __m128 XZWY = _mm_addsub_ps(_mm_shuffle_ps(XZYnW, ZnXWY, _MM_SHUFFLE(3,2,1,0)),
                                _mm_shuffle_ps(ZnXWY, XZYnW, _MM_SHUFFLE(2,3,0,1)));

    /* now we only need to shuffle the components in place and return the result      */
    return _mm_shuffle_ps(XZWY, XZWY, _MM_SHUFFLE(2,1,3,0));

    /* operations: 6 shuffles, 4 multiplications, 3 compound additions/subtractions   */

我希望程序集包含最少的指令。但是，当我使用 gcc -msse3 -S 将其编译为程序集时，生成的函数有 67 条指令：

    .text   
    .globl __Z13_mm_cross4_psU8__vectorfS_
__Z13_mm_cross4_psU8__vectorfS_:
LFB594:
    pushq   %rbp
LCFI0:
    movq    %rsp, %rbp
LCFI1:
    subq    $232, %rsp
    movaps  %xmm0, -336(%rbp)
    movaps  %xmm1, -352(%rbp)
    movaps  -336(%rbp), %xmm0
    movaps  -336(%rbp), %xmm1
    shufps  $27, %xmm1, %xmm0
    movaps  %xmm0, -16(%rbp)
    movaps  -352(%rbp), %xmm0
    movaps  -352(%rbp), %xmm1
    shufps  $17, %xmm1, %xmm0
    movaps  %xmm0, -32(%rbp)
    movaps  -352(%rbp), %xmm0
    movaps  -352(%rbp), %xmm1
    shufps  $187, %xmm1, %xmm0
    movaps  %xmm0, -48(%rbp)
    movaps  -16(%rbp), %xmm0
    movaps  %xmm0, -112(%rbp)
    movaps  -48(%rbp), %xmm0
    movaps  %xmm0, -128(%rbp)
    movaps  -128(%rbp), %xmm0
    movaps  -112(%rbp), %xmm1
    mulps   %xmm1, %xmm0
    movaps  -336(%rbp), %xmm1
    movaps  %xmm1, -144(%rbp)
    movaps  -32(%rbp), %xmm1
    movaps  %xmm1, -160(%rbp)
    movaps  -160(%rbp), %xmm1
    movaps  -144(%rbp), %xmm2
    mulps   %xmm2, %xmm1
    movaps  %xmm1, -176(%rbp)
    movaps  %xmm0, -192(%rbp)
    movaps  -176(%rbp), %xmm0
    hsubps  -192(%rbp), %xmm0
    movaps  %xmm0, -64(%rbp)
    movaps  -16(%rbp), %xmm0
    movaps  %xmm0, -208(%rbp)
    movaps  -32(%rbp), %xmm0
    movaps  %xmm0, -224(%rbp)
    movaps  -224(%rbp), %xmm0
    movaps  -208(%rbp), %xmm1
    mulps   %xmm1, %xmm0
    movaps  -336(%rbp), %xmm1
    movaps  %xmm1, -240(%rbp)
    movaps  -48(%rbp), %xmm1
    movaps  %xmm1, -256(%rbp)
    movaps  -256(%rbp), %xmm1
    movaps  -240(%rbp), %xmm2
    mulps   %xmm2, %xmm1
    movaps  %xmm1, -272(%rbp)
    movaps  %xmm0, -288(%rbp)
    movaps  -272(%rbp), %xmm0
    haddps  -288(%rbp), %xmm0
    movaps  %xmm0, -80(%rbp)
    movaps  -64(%rbp), %xmm0
    movaps  -80(%rbp), %xmm1
    shufps  $177, %xmm1, %xmm0
    movaps  -80(%rbp), %xmm1
    movaps  -64(%rbp), %xmm2
    shufps  $228, %xmm2, %xmm1
    movaps  %xmm1, -304(%rbp)
    movaps  %xmm0, -320(%rbp)
    movaps  -304(%rbp), %xmm0
    addsubps        -320(%rbp), %xmm0
    movaps  %xmm0, -96(%rbp)
    movaps  -96(%rbp), %xmm0
    movaps  -96(%rbp), %xmm1
    shufps  $156, %xmm1, %xmm0
    leave
LCFI2:
    ret

我做错了什么？必须有更好的方法来洗牌周围的元素，而无需使用这么多 movaps 指令。

【问题讨论】：

要充分利用 SIMD，您应该并行执行多项操作，而不是一项操作。对于 SSE，这意味着您应该一次执行四个四元数乘法（乘以八个四元数）。这样您就不必使用任何水平运算符，并且您应该获得 4 倍于非 SIMD 代码的速度。

尽量避免水平操作——我认为将计算安排为xyzw.*dddd + yzwx.*cacA + zwxy.*BbbB - wxyz.*Acac和(A=-a, B=-b)的垂直顺序时可以使用相同数量的指令。

我明白你要说什么，但是一个四元数已经由四个浮点数组成，而且通常在我的代码中我没有一个循环有很多乘法，但一次只有一个乘法，所以一次做几个四元数并没有帮助。

【参考方案1】：

没关系。如果我用 gcc -msse3 -O1 -S 编译代码，我会得到以下信息：

    .text
    .align 4,0x90
    .globl __Z13_mm_cross4_psU8__vectorfS_
__Z13_mm_cross4_psU8__vectorfS_:
LFB644:
    movaps  %xmm0, %xmm5
    movaps  %xmm1, %xmm3
    movaps  %xmm0, %xmm2
    shufps  $27, %xmm0, %xmm5
    movaps  %xmm5, %xmm4
    shufps  $17, %xmm1, %xmm3
    shufps  $187, %xmm1, %xmm1
    mulps   %xmm3, %xmm2
    mulps   %xmm1, %xmm4
    mulps   %xmm5, %xmm3
    mulps   %xmm1, %xmm0
    hsubps  %xmm4, %xmm2
    haddps  %xmm3, %xmm0
    movaps  %xmm2, %xmm1
    shufps  $177, %xmm0, %xmm1
    shufps  $228, %xmm2, %xmm0
    addsubps        %xmm1, %xmm0
    shufps  $156, %xmm0, %xmm0
    ret

现在只有 18 条指令。这就是我一开始所期望的。哎呀。

【讨论】：

Nrgh！切勿在未启用优化的情况下衡量任何性能！

【参考方案2】：

您可能对 Agner Fog 的C++ vector class library 感兴趣。它提供了Quaternion4f 和Quaternion4d 类（当然包括* 和*= 运算符），分别使用SSE2 和AVX 指令集实现。该库是一个开源项目，因此您可以深入研究代码并找到一个很好的实现示例来构建您的函数。

稍后，您可以咨询"optimizing subroutines in assembly language" manual 并提供该函数的优化、纯汇编实现，或者在了解一些低级技巧的同时，尝试重新设计 C 中的内在方法。

【讨论】：

我现在正在对矩阵乘法进行类似的研究。祝你好运！

感谢您的指点。我以前看过“优化汇编语言中的子程序”文档，但后来我把它弄丢了。