超快速舍入功能 (PBC)答案

【问题标题】：Super-fast rounding function (PBC)超快速舍入功能 (PBC)
【发布时间】：2016-04-08 11:10:08
【问题描述】：

我真的需要非常快的 C 中的 round() 函数 - 蒙特卡洛粒子建模是必要的：在每一步，您都需要将坐标包装到周期框中以计算体积交互：例如

for(int i=0; i < 3; i++)
{
    coor.x[i] = a.XReal.x[i]-b.XReal.x[i];
    coor.x[i] = coor.x[i] - SIZE[i]*round(coor.x[i]/SIZE[i]); //PBC
}

我遇到过一些 asm hacking，但我根本不懂 asm:) 像这样的

inline int float2int2(float flt)
{
  int intgr;

  __asm__ __volatile__ ("fld %1; fistp %0;" : "=m" (intgr) : "m" (flt));

  return intgr;
}

有了固定的边界，没有 round() 它工作得更快。那么，也许有人知道更好的方法？..

【问题讨论】：

您找到的 asm 几乎毫无用处：它是特定于体系结构的 (x86) 并使用旧的 x87 fpu（慢，可怕）。如果您使用的是 x86-64，则最好使用 sse 指令，任何理智的编译器都会自动发出。
谢谢！我忘了提 - 需要最接近整数舍入。比如 0.8 -> 1, -2.7 -> -3
您是否尝试过执行此操作的标准函数 (remainder()) 在性能方面的比较？
round() 舍入“到最近，从零开始”，可能更快的rint() 舍入“到最近，从零开始”。如果您可以容忍差异，请尝试rint()。性能方面，我会更关心这段代码中的浮点除法。什么是 SIZE[i] 的典型值，这些是编译时常量吗？
这个问题有一个非常有用的技巧，适用于有限的输入范围：stackoverflow.com/questions/17035464/…

标签： c performance floating-point modeling

【解决方案1】：

首先，您可以通过使用正确的编译器选项获得一些收益。例如，使用 GCC 和现代 Intel CPU，您应该尝试：

-march=nehalem -fno-trapping-math

那么round 的问题在于它使用了特定的舍入模式，这在大多数平台上都很慢。 nearbyint（或rint）应该总是更快：

coor.x[i] = coor.x[i] - SIZE[i] * nearbyint(coor.x[i] / SIZE[i])

看看generated assembly。

您还应该考虑对代码进行矢量化。

【讨论】：

-ffast-math 不是更好的选择吗？它包括-funsafe-math-optimizations 以及其他一些选项。
@DavidWohlferd -ffast-math 不会加快舍入速度。实际上，使用-fno-trapping-math 让GCC 生成roundsd 指令就足够了。我编辑了我的答案。
此外，即使没有 -fno-trapping-math，clang 也会为 nearbyint 发出 roundsd。这似乎是正确的，因为 nearbyint 不会引发浮点异常。

【解决方案2】：

理想情况下，您希望将范围缩小到周期框的整个过程快速进行，而不是仅仅寻找快速舍入。正如@EOF 在评论中准确指出的那样，您可以使用 C99 标准函数，例如 remainderf() 或 fmodf()。

coor.x[i] -= SIZE[i]*round(coor.x[i]/SIZE[i]);
// same as
coor.x[i] = remainderf(coor.x[i], SIZE[i]);

fmodf(3) 向零舍入，remainderf(3) rounds towards nearest。

remainder() 函数计算x 除以y 的余数。返回值为x-n*y，其中n为值x / y，四舍五入到最接近的整数。如果x-n*y的绝对值为0.5，则n选择为偶数。

编译器/库有几种不同的策略来实现这些。使用-ffast-math，x86-64 的gcc 5.3 内联了remainder(x,y) 实现，它将值从SSE 寄存器传输到x87 寄存器，并在循环中运行FPREM1（部分余数），直到它设置一个指示结果为的标志正确的。（FPREM1 一次执行最多可以将指数减少 63）。

clang 总是发出对库函数的调用，可以是普通的remainder 入口点，也可以是__remainder_finite 和-ffast-math。

GNU libm 定义主要使用整数运算，来自反汇编 and the C source 的 AFAICT。在具有快速硬件划分的最新 Intel CPU 上，它可能比您的 div、round、mul 版本慢。

所以你有三个选择：

div、round、mul、sub，快速舍入（使用nearbyint()，它显然具有最不难看的语义，因此它可以最容易地内联到roundsd/roundss）。 这种方式可以矢量化，一次完成所有三个坐标。可能需要手动完成，以找到不会出现第 4 个元素的错误。在具有 128b 向量的 Intel Haswell 上：5 uop。单精度：divps（10-13c 延迟，每 7c 吞吐量一个），roundps（2 uops，6c 延迟，每 2c 吞吐量一个），mulps（5c 延迟，每 0.5c 吞吐量一个），@ 987654350@（3c 延迟，每 1c 吞吐量一个）。其中一些相互竞争执行端口。 总延迟：27c。可能的吞吐量，可能类似于 每 7c 一个（完全受到 divps 的限制）
gcc 的内联 x87 FPREM1。（可能只需要运行一次迭代，所以在 Haswell 上：41 uops，27c 延迟，每 17c 吞吐量一个，加上在 xmm 和 x87 regs 之间获取数据的一些开销。无法矢量化。
glibc 的大部分整数实现：在现代 x86 CPU 上不知道，可能比其他两个更差。但是，probably significantly higher accuracy比手动的div/round/mul/sub。

底线，如果这是一个速度问题，您应该一定要考虑使用 SSE/AVX 进行矢量化，以便在一个矢量中处理一个点的所有三个坐标。或者，一次四个点的坐标，或者任何方便的东西。理想情况下，您可以使用向量 ALU 的所有 4 个（或 8 个 AVX）单精度元素。（或 2 / 4 表示双精度）。

即使是标量，我认为您当前使用 nearbyint() 的代码将是最快的选择，但您可以轻松地比使用向量的代码快三倍。

【讨论】：