在间隔中包装值的最快方法

Question

我很好奇将浮点值 x 包装在半封闭区间 [0; a[ 中的方法。

例如，我可以有一个任意实数，比如x = 354638.515，我希望将它折叠成[0; 2π[，因为我有一个很好的sin 近似值。

fmod 标准 C 函数在我的基准测试中显示得相当高，通过检查各种 libc 实现的源代码，我可以理解原因：这东西相当分支，可能是为了处理很多IEEE754 特定问题：

• glibc:https://github.com/bminor/glibc/blob/master/sysdeps/ieee754/flt-32/e_fmodf.c
• 苹果：https://opensource.apple.com/source/Libm/Libm-315/Source/ARM/fmod.c.auto.html
• musl:https://git.musl-libc.org/cgit/musl/tree/src/math/fmod.c
• 使用-ffast-math 运行会导致 GCC 生成通过 x86/x86_64 上的 x87 FPU 的代码，这会带来一系列问题（例如 80 位双精度、FP 状态和其他有趣的事情）。我希望实现至少半正确地进行矢量化，并且如果可能的话，至少不要通过 x87 FPU，而是通过矢量寄存器，因为我的其余代码最终会被矢量化，即使不一定是最佳方式，由编译器.
• 这个看起来简单多了：https://github.com/KnightOS/libc/blob/master/src/fmod.c

就我而言，我只关心通常的实数值，而不是 NaN，也不关心无穷大。我的范围在编译时也是已知的并且是正常的（常见的情况是 π/2），因此不需要检查诸如 range == 0 之类的“特殊情况”。

因此，对于该特定用例，fmod 的良好实现是什么？

Answer 1

假设范围是恒定且正的，您可以计算其倒数以避免昂贵的除法。

void fast_fmod(float * restrict dst, const float * restrict src, size_t n, float divisor) {
   float reciprocal = 1.0f / divisor;
   for (size_t i = 0; i < n; ++i)
     dst[i] = src[i] - divisor * (int)(src[i] * reciprocal);
}

最后附上一个简单demo的代码是：

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

void fast_fmod(float * restrict dst, const float * restrict src, size_t n, float divisor) {
   float reciprocal = 1.0f / divisor;

   for (size_t i = 0; i < n; ++i)
     dst[i] = src[i] - divisor * (int)(src[i] * reciprocal);
}

int main() {
    float src[9] = {-4, -3, -2, -1, 0, 1, 2, 3, 4};
    float dst[9];
    float div = 3;

    fast_fmod(dst, src, 9, div);

    for (int i = 0; i < 9; ++i) {
        printf("fmod(%g, %g) = %g vs %g\n", src[i], div, dst[i], fmod(src[i], div));
    }
}

产生预期的输出：

fmod(-4, 3) = -1 vs -1
fmod(-3, 3) = 0 vs -0
fmod(-2, 3) = -2 vs -2
fmod(-1, 3) = -1 vs -1
fmod(0, 3) = 0 vs 0
fmod(1, 3) = 1 vs 1
fmod(2, 3) = 2 vs 2
fmod(3, 3) = 0 vs 0
fmod(4, 3) = 1 vs 1

使用 GCC 编译命令：

$ gcc prog.c -o prog -O3 -march=haswell -lm -fopt-info-vec
prog.c:8:4: optimized: loop vectorized using 32 byte vectors
prog.c:8:4: optimized: loop vectorized using 32 byte vectors
prog.c:8:30: optimized: basic block part vectorized using 32 byte vectors

因此代码被很好地矢量化了。

---

编辑

看起来 CLANG 在矢量化这段代码方面做得更好：

  401170:   c5 fc 10 24 8e          vmovups (%rsi,%rcx,4),%ymm4
  401175:   c5 fc 10 6c 8e 20       vmovups 0x20(%rsi,%rcx,4),%ymm5
  40117b:   c5 fc 10 74 8e 40       vmovups 0x40(%rsi,%rcx,4),%ymm6
  401181:   c5 fc 10 7c 8e 60       vmovups 0x60(%rsi,%rcx,4),%ymm7
  401187:   c5 6c 59 c4             vmulps %ymm4,%ymm2,%ymm8
  40118b:   c5 6c 59 cd             vmulps %ymm5,%ymm2,%ymm9
  40118f:   c5 6c 59 d6             vmulps %ymm6,%ymm2,%ymm10
  401193:   c5 6c 59 df             vmulps %ymm7,%ymm2,%ymm11
  401197:   c4 41 7e 5b c0          vcvttps2dq %ymm8,%ymm8
  40119c:   c4 41 7e 5b c9          vcvttps2dq %ymm9,%ymm9
  4011a1:   c4 41 7e 5b d2          vcvttps2dq %ymm10,%ymm10
  4011a6:   c4 41 7e 5b db          vcvttps2dq %ymm11,%ymm11
  4011ab:   c4 41 7c 5b c0          vcvtdq2ps %ymm8,%ymm8
  4011b0:   c4 41 7c 5b c9          vcvtdq2ps %ymm9,%ymm9
  4011b5:   c4 41 7c 5b d2          vcvtdq2ps %ymm10,%ymm10
  4011ba:   c4 41 7c 5b db          vcvtdq2ps %ymm11,%ymm11
  4011bf:   c5 3c 59 c3             vmulps %ymm3,%ymm8,%ymm8
  4011c3:   c5 34 59 cb             vmulps %ymm3,%ymm9,%ymm9
  4011c7:   c5 2c 59 d3             vmulps %ymm3,%ymm10,%ymm10
  4011cb:   c5 24 59 db             vmulps %ymm3,%ymm11,%ymm11
  4011cf:   c4 c1 5c 5c e0          vsubps %ymm8,%ymm4,%ymm4
  4011d4:   c4 c1 54 5c e9          vsubps %ymm9,%ymm5,%ymm5
  4011d9:   c4 c1 4c 5c f2          vsubps %ymm10,%ymm6,%ymm6
  4011de:   c4 c1 44 5c fb          vsubps %ymm11,%ymm7,%ymm7
  4011e3:   c5 fc 11 24 8f          vmovups %ymm4,(%rdi,%rcx,4)
  4011e8:   c5 fc 11 6c 8f 20       vmovups %ymm5,0x20(%rdi,%rcx,4)
  4011ee:   c5 fc 11 74 8f 40       vmovups %ymm6,0x40(%rdi,%rcx,4)
  4011f4:   c5 fc 11 7c 8f 60       vmovups %ymm7,0x60(%rdi,%rcx,4)
  4011fa:   48 83 c1 20             add    $0x20,%rcx
  4011fe:   48 39 c8                cmp    %rcx,%rax
  401201:   0f 85 69 ff ff ff       jne    401170 <fast_fmod+0x40>