我很好奇如何將一個浮點值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），因此對 "特殊情況 "的檢查是不必要的，例如范圍 == 0。
因此，對于這個特定的用例，什么才是fmod的好實作？
uj5u.com熱心網友回復：
假設范圍是常數和正數，你可以計算其倒數以避免昂貴的除法。
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）。
}
帶有簡單演示的最終代碼是：
#include <stdlib.h>/span>
#include <stdio.h>
#include <math.h>/span>
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
", 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: 優化：回圈向量使用 32 位元組向量
prog.c:8:4: 優化：回圈向量使用 32位元組向量
prog.c:8:30: 優化：基本塊部分矢量化 使用 32位元組矢量
這樣一來，代碼就被很好地矢量化了。
。
編輯
看起來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 59cd 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 5bdb 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 5bdb 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 jne 401170 <fast_fmod 0x40>
轉載請註明出處，本文鏈接：https://www.uj5u.com/shujuku/331102.html
標籤：