规则除法和乘法一样快，近似倒数。为什么？

Question

我有以下代码：

void division_approximate(float a[], float b[], float c[], int n) {
    // c[i] = a[i] * (1 / b[i]);
    for (int i = 0; i < n; i+=8) {
        __m256 b_val = _mm256_loadu_ps(b + i);
        b_val = _mm256_rcp_ps(b_val);
        __m256 a_val = _mm256_loadu_ps(a + i);
        a_val = _mm256_mul_ps(a_val, b_val);
        _mm256_storeu_ps(c + i, a_val);
    }
}

void division(float a[], float b[], float c[], int n) {
    // c[i] = a[i] / b[i];
    for (int i = 0; i < n; i+=8) {
        __m256 b_val = _mm256_loadu_ps(b + i);
        __m256 a_val = _mm256_loadu_ps(a + i);
        a_val = _mm256_div_ps(a_val, b_val);
        _mm256_storeu_ps(c + i, a_val);
    }
}

我希望division_approximate比division更快，但这两个功能在我的AMD 7,4800H上的时间几乎相同。我不明白为什么，我认为division_approximate要快得多。这个问题在GCC身上都重现了。用-O3 -march=core-avx2编译。

更新

下面是GCC 9.3为这两个循环生成的源代码：

division
│  >0x555555555c38 <division+88>    vmovups 0x0(%r13,%rax,4),%ymm3                                                                                                                                                                                                                       │
│   0x555555555c3f <division+95>    vdivps (%r14,%rax,4),%ymm3,%ymm0                                                                                                                                                                                                                     │
│   0x555555555c45 <division+101>   vmovups %ymm0,(%rbx,%rax,4)                                                                                                                                                                                                                          │
│   0x555555555c4a <division+106>   add    $0x8,%rax                                                                                                                                                                                                                                     │
│   0x555555555c4e <division+110>   cmp    %eax,%r12d                                                                                                                                                                                                                                    │
│   0x555555555c51 <division+113>   jg     0x555555555c38 <division+88>                                                                                                                                                                                                                  │

division_approximate
│  >0x555555555b38 <division_approximate+88>        vrcpps (%r14,%rax,4),%ymm0                                                                                                                                                                                                           │
│   0x555555555b3e <division_approximate+94>        vmulps 0x0(%r13,%rax,4),%ymm0,%ymm0                                                                                                                                                                                                  │
│   0x555555555b45 <division_approximate+101>       vmovups %ymm0,(%rbx,%rax,4)                                                                                                                                                                                                          │
│   0x555555555b4a <division_approximate+106>       add    $0x8,%rax                                                                                                                                                                                                                     │
│   0x555555555b4e <division_approximate+110>       cmp    %eax,%r12d                                                                                                                                                                                                                    │
│   0x555555555b51 <division_approximate+113>       jg     0x555555555b38 <division_approximate+88>                                                                                                                                                                                      │

对于n = 256 * 1024 * 1024，这两个代码的执行时间几乎完全相同(318毫秒对319毫秒)。

Answer 1

(256 * 1024 * 1024) * 4 (每个浮点数) / 0.318 / 1000^2 * 4约为13.5GB/ stream带宽，或约10.1GB/s有用的流带宽。(假设商店实际上为RFO花费了read+write带宽；正如杰罗姆指出的那样，_mm256_stream_ps可以让商店只花一次，而不是两次。)

如果这对你的禅宗2上的单线程三位一体带宽是好的还是坏的，那只是

(256 * 1024 * 1024 / 8) / 0.318 / 1000^3 = ~0.1055个载体(8个浮子)/纳秒，Zen 2 vdivps能保持在0.36 GHz。我想你的CPU比那快:P

(0.1055个vec/ns *3.5个循环/vec= 0.36循环/ns(又名GHz) )

是一个非常明显的内存瓶颈，与Zen2 2的每3.5周期vdivps ymm吞吐量之一相去甚远。(https://uops.info/)。使用一个小得多的数组(适用于L1或至少L2缓存)，并多次遍历它。

尽量避免在实际代码中编写这样的循环，的计算强度(每次将数据加载到L1缓存或寄存器时的工作量)非常低。作为另一次传递的一部分执行此操作，或者使用缓存阻塞对一小部分输入执行此操作，然后在缓存中仍处于热状态时使用输出的这一小部分。(这比使用_mm256_stream_ps绕过缓存要好得多。)

当与其他操作(大量fmas / mul / add)混合时，vdivps通常比rcpps +牛顿迭代(通常需要获得可接受的精度:原始rcpps仅为11位rcpps)更好选择，而vdivps只是单个uop，而单独的rcpps和覆盖物uop。(尽管从内存中，vdivps仍然需要单独的vmovups加载，而且Zen在将内存源折叠到单个uop中没有问题)。还请参阅Floating point division vs floating point multiplication re:前端吞吐量与划分单元瓶颈(如果您只是分割，而不是将其与其他操作混合)。

当然，如果你能完全避免除法的话，这是很棒的，例如，把一个倒数从一个循环中提升出来，然后只进行乘法，但是现代CPU有足够好的除法器HW，即使你没有内存瓶颈，从rcpps中也没有什么收获。例如，在使用两个多项式的比率来评估多项式逼近时，FMA的数量通常足以隐藏vdivps的吞吐量成本，而牛顿迭代将花费更多的FMA uops。

另外，当你没有英特尔的“核心”微架构时，为什么要使用-march=core-avx2呢？使用-march=native或-march=znver2。除非您有意地对运行在AMD CPU上的Intel的二进制文件进行基准测试。