超快舍入函数(PBC)

Question

我真的需要在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黑客攻击它，但我完全不理解asm：)类似这样的东西

inline int float2int2(float flt)
{
  int intgr;

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

  return intgr;
}

有了固定的边界，没有round()，它的工作速度更快。所以，也许有人知道更好的方法？

Answer 1

首先，通过使用正确的编译器选项，您可以获得一些收益。以GCC和现代英特尔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吧。

你也应该考虑向量化你的代码。

Answer 2

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

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入口点，要么是带有-ffast-math的__remainder_finite。

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

因此，您有三个选择：

• div、round、mul、sub，具有快速舍入(使用nearbyint()，显然它具有最简单的丑陋语义，所以它可以最容易地内联到roundsd / roundss )。这种方法可以矢量化，并且可以一次完成所有的三个坐标。可能需要手动完成，才能找到第四个元素不会出错的地方。在具有128b矢量的英特尔Haswell上:5个uops。单精度：divps(10-13c延迟，每7c吞吐量一个)，roundps(2个uops，6c延迟，每2c吞吐量一个)，mulps(5c延迟，每0.5c吞吐量一个)，subps(3c延迟，每1c吞吐量一个)。其中一些会相互竞争执行端口。总延迟: 27c。可能的吞吐量，可能类似于one per 7c (完全由divps造成的瓶颈)
• ‘s inlined x87 FPREM1。(可能只需要运行一次迭代，因此Haswell: 41uops，27c延迟，每17c吞吐量一次，外加在xmm和x87 regs之间获取数据的一些开销。Can't vectorize.
• glibc's mostly integer implementation:不知道，在现代x86 CPU上可能比其他两种都差。但是，与手动div/round/mul/sub.

相比，

最重要的是，如果这是一个速度问题，你绝对应该研究一下使用SSE/AVX进行矢量化，以便在一个矢量中完成一个点的所有三个坐标。或者，一次四个点的坐标，或者任何方便的东西。理想情况下，您可以使用向量all的所有4个(或AVX为8个)单精度元素。(或2/4表示双精度)。

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