通过使用 AVX 内部函数重写来提高 math.h 函数的性能

Posted 2023-02-16

【中文标题】通过使用 AVX 内部函数重写来提高 math.h 函数的性能

【英文标题】：Performance improvement of math.h functions by rewriting with AVX intrinsics

【发布时间】：2020-06-11 10:26:03

【问题描述】：

我有一个简单的数学库，它链接到在模拟器硬件（32 位 RTOS）上运行的项目中，并且编译器工具链基于 GCC 5.5 的变体。主要项目代码在 Matlab 中，但核心数学运算（数组数据上的 cmath 函数）是用 C 重新编写的以提高性能。查看编译器资源管理器，优化代码的质量对于GCC 5.5 32 bit 似乎不是很好（供参考：Clang trunk 32bit）。据我了解，Clang 在优化循环方面做得更好。一个示例代码sn-p：

...

void cfunctionsLog10(unsigned int n, const double* x, double* y) 
    int i;
    for (i = 0; i < n; i++) 
        y[i] = log10(x[i]);
    

以及GCC 5.5生成的相应程序集

cfunctionsLog10(unsigned int, double const*, double*):
        push    ebp
        push    edi
        push    esi
        push    ebx
        sub     esp, 12
        mov     esi, DWORD PTR [esp+32]
        mov     ebp, DWORD PTR [esp+36]
        mov     edi, DWORD PTR [esp+40]
        test    esi, esi
        je      .L28
        xor     ebx, ebx
.L27:
        sub     esp, 8
        push    DWORD PTR [ebp+4+ebx*8]
        push    DWORD PTR [ebp+0+ebx*8]
        call    __log10_finite
        fstp    QWORD PTR [edi+ebx*8]
        add     ebx, 1
        add     esp, 16
        cmp     ebx, esi
        jne     .L27
.L28:
        add     esp, 12
        pop     ebx
        pop     esi
        pop     edi
        pop     ebp
        ret

Clang 产生的地方：

cfunctionsLog10(unsigned int, double const*, double*):              # @cfunctionsLog10(unsigned int, double const*, double*)
        push    ebp
        push    ebx
        push    edi
        push    esi
        sub     esp, 76
        mov     esi, dword ptr [esp + 96]
        test    esi, esi
        je      .LBB2_8
        mov     edi, dword ptr [esp + 104]
        mov     ebx, dword ptr [esp + 100]
        xor     ebp, ebp
        cmp     esi, 4
        jb      .LBB2_7
        lea     eax, [ebx + 8*esi]
        cmp     eax, edi
        jbe     .LBB2_4
        lea     eax, [edi + 8*esi]
        cmp     eax, ebx
        ja      .LBB2_7
.LBB2_4:
        mov     ebp, esi
        xor     esi, esi
        and     ebp, -4
.LBB2_5:                                # =>This Inner Loop Header: Depth=1
        vmovsd  xmm0, qword ptr [ebx + 8*esi + 16] # xmm0 = mem[0],zero
        vmovsd  qword ptr [esp], xmm0
        vmovsd  xmm0, qword ptr [ebx + 8*esi] # xmm0 = mem[0],zero
        vmovsd  xmm1, qword ptr [ebx + 8*esi + 8] # xmm1 = mem[0],zero
        vmovsd  qword ptr [esp + 8], xmm0 # 8-byte Spill
        vmovsd  qword ptr [esp + 16], xmm1 # 8-byte Spill
        call    log10
        fstp    tbyte ptr [esp + 64]    # 10-byte Folded Spill
        vmovsd  xmm0, qword ptr [esp + 16] # 8-byte Reload
        vmovsd  qword ptr [esp], xmm0
        call    log10
        fstp    tbyte ptr [esp + 16]    # 10-byte Folded Spill
        vmovsd  xmm0, qword ptr [esp + 8] # 8-byte Reload
        vmovsd  qword ptr [esp], xmm0
        vmovsd  xmm0, qword ptr [ebx + 8*esi + 24] # xmm0 = mem[0],zero
        vmovsd  qword ptr [esp + 8], xmm0 # 8-byte Spill
        call    log10
        vmovsd  xmm0, qword ptr [esp + 8] # 8-byte Reload
        vmovsd  qword ptr [esp], xmm0
        fstp    qword ptr [esp + 56]
        fld     tbyte ptr [esp + 16]    # 10-byte Folded Reload
        fstp    qword ptr [esp + 48]
        fld     tbyte ptr [esp + 64]    # 10-byte Folded Reload
        fstp    qword ptr [esp + 40]
        call    log10
        fstp    qword ptr [esp + 32]
        vmovsd  xmm0, qword ptr [esp + 56] # xmm0 = mem[0],zero
        vmovsd  xmm1, qword ptr [esp + 40] # xmm1 = mem[0],zero
        vmovhps xmm0, xmm0, qword ptr [esp + 48] # xmm0 = xmm0[0,1],mem[0,1]
        vmovhps xmm1, xmm1, qword ptr [esp + 32] # xmm1 = xmm1[0,1],mem[0,1]
        vmovups xmmword ptr [edi + 8*esi + 16], xmm1
        vmovups xmmword ptr [edi + 8*esi], xmm0
        add     esi, 4
        cmp     ebp, esi
        jne     .LBB2_5
        mov     esi, dword ptr [esp + 96]
        cmp     ebp, esi
        je      .LBB2_8
.LBB2_7:                                # =>This Inner Loop Header: Depth=1
        vmovsd  xmm0, qword ptr [ebx + 8*ebp] # xmm0 = mem[0],zero
        vmovsd  qword ptr [esp], xmm0
        call    log10
        fstp    qword ptr [edi + 8*ebp]
        inc     ebp
        cmp     esi, ebp
        jne     .LBB2_7
.LBB2_8:
        add     esp, 76
        pop     esi
        pop     edi
        pop     ebx
        pop     ebp
        ret

由于我无法直接使用 Clang，因此使用 AVX 内部函数重写 C 源代码是否有任何价值。我认为大部分性能成本来自 cmath 函数调用，其中大部分没有内在实现。

---

编辑： 使用vectorclass library重新实现：

void vclfunctionsTanh(unsigned int n, const double* x, double* y) 

    const int N = n;
    const int VectorSize = 4;
    const int FirstPass = N & (-VectorSize);

    int i = 0;    
    for (; i < FirstPass; i+= 4)
    
        Vec4d data = Vec4d.load(x[i]);
        Vec4d ans = tanh(data);
        ans.store(y+i);
    

    
    for (;i < N; ++i)
        y[i]=std::tanh(x[i]);

【问题讨论】：

首先，在您的问题中包含一些您正在谈论的实际代码 + asm，而不仅仅是通过 Godbolt 短链接。 （如果您没有包括所有内容，请使用指向 Godbolt 的“完整”链接来防止比特腐烂）。此外，使用更新的 GCC 来更有效地自动矢量化 sqrt 和/或-mno-avx256-split-unaligned-load -mno-avx256-split-unaligned-store（请参阅Why doesn't gcc resolve _mm256_loadu_pd as single vmovupd?）。也可以试试-mfpmath=sse？调用约定仍然在堆栈上传递并在 x87 中返回，除非您重新编译 libm。

如果您只需要其中一些函数的快速近似值，您可以手动对它们进行矢量化并将它们与内在函数内联；这可能会大大加快速度。尽管 AVX 和 x87 进行存储/重新加载的混合效率低下，但大部分成本将在数学库函数中，除了 sqrt。所以是的，并行执行 4 个双打可能至少快 4 倍，或者如果您以一些精度换取速度，则速度会更快。 （和/或 NaN 检查；如果您知道您的输入是有限的，这对您无法按元素进行分支的 SIMD 很有帮助）

是的，这取决于您如何编写它们，它们当然可以。 log 和 exp 相对容易向量化，尾数的多项式逼近很快收敛。例如Efficient implementation of log2(__m256d) in AVX2 / Fastest Implementation of Exponential Function Using AVX 及相关问答。 IDK关于cos。对于float，您可以简单地测试每个可能的 32 位浮点位模式以确定最小/最大错误。

还有 SIMD 数学库和现有实现，如 Where is Clang's '_mm256_pow_ps' intrinsic? / AVX log intrinsics (_mm256_log_ps) missing in g++-4.8? / Mathematical functions for SIMD registers

您的 C cfunctionsTanh 是与您显示的 asm 不同的函数，调用 cfunctionsLog10。此外，clang 的 asm 实际上非常糟糕，将返回值 tbyte 溢出/重新加载到堆栈而不是直接转换到目标。 10 字节 x87 存储/重新加载比 double 慢得多，完全没有必要。

【参考方案1】：

矢量类库具有常用数学函数的内联矢量版本，包括 log10。

https://github.com/vectorclass/

【讨论】：

文档建议不要使用-ffast-math 进行编译，对于我们知道我们只处理有限数和不会溢出的操作的情况是否如此？

-ffast-math 在大多数情况下无法检测到 NAN 结果。如果您不想检查 NAN，那么您可以使用 -ffast-math。