计算浮点余数的最佳方法

Best way calculating remainder on floating points

本文关键字：最佳 方法 余数 计算      更新时间：2023-10-16

问题

remainder函数的计算成本是多少？是否有特定的指令可以在特定情况下以廉价的方式计算它？

说明

我需要将数学变量x从R转换为I=[-0.5；0.5(的范围=[-2；2(。当x不是I的元素时，通过对x值重复加或减1，x-向I移动。x-在我的代码中用double x表示。我需要最快的方法来转换I和R，但更宽的R范围也可能很有趣。

想法和速度比较

建议我使用的功能是以下描述的天真实现：

void shift_to_I(double& x) // version 1
{
while (x < -0.5)
x += 1;
while (x >= 0.5)
x -= 1;
}

不仅出于速度考虑，而且出于代码质量考虑，我考虑使用c++11中引入的<cmath>中的remainder。使用remainder，代码缩短为

void shift_to_I(double& x) // version 2
{
x = remainder(x,1);
}

不过，我不得不意识到，它比我的体系结构上的原始功能(英特尔i7 whatver与VC++(慢。我相信有一个专门的指令用于此目的，但要么编译器不知道，要么它不存在。对于更宽的R间隔(在我的体系结构上大约是[-25；25((，第二个版本会更快，但我也需要一个对于窄间隔来说很快的代码。clang和gcc特定的解决方案也很受欢迎。

这个问题取决于编译器和实现。

例如，在我的GCC 8.3:机器上

  1. 如果没有-ffast-math，std::remainder将转换为此函数的调用：

    double __remainder(double x, double y)
{
if (((__builtin_expect (y == 0.0, 0) && ! isnan(x)) || (__builtin_expect(isinf(x), 0) && ! isnan(y))) && _LIB_VERSION != _IEEE_)
return __kernel_standard(x, y, 28);
return __ieee754_remainder(x, y);
}

    __ieee754_remainder看起来像这样：

    double __ieee754_remainder(double x, double y)
{
double z, d, xx;
int4 kx, ky, n, nn, n1, m1, l;
mynumber u, t, w = {{0, 0}}, v = {{0, 0}}, ww = {{0, 0}}, r;
u.x = x;
t.x = y;
kx = u.i[HIGH_HALF] & 0x7fffffff; /* no sign  for x*/
t.i[HIGH_HALF] &= 0x7fffffff;     /*no sign for y */
ky = t.i[HIGH_HALF];
/*------ |x| < 2^1023  and   2^-970 < |y| < 2^1024 ------------------*/
if (kx < 0x7fe00000 && ky < 0x7ff00000 && ky >= 0x03500000)
{
SET_RESTORE_ROUND_NOEX(FE_TONEAREST);
if (kx + 0x00100000 < ky)
return x;
if ((kx - 0x01500000) < ky)
{
z = x / t.x;
v.i[HIGH_HALF] = t.i[HIGH_HALF];
d = (z + big.x) - big.x;
xx = (x - d * v.x) - d * (t.x - v.x);
if (d - z != 0.5 && d - z != -0.5)
return (xx != 0) ? xx : ((x > 0) ? ZERO.x : nZERO.x);
else
{
if (fabs(xx) > 0.5 * t.x)
return (z > d) ? xx - t.x : xx + t.x;
else
return xx;
}
} /*    (kx<(ky+0x01500000))         */
else
{
r.x = 1.0 / t.x;
n = t.i[HIGH_HALF];
nn = (n & 0x7ff00000) + 0x01400000;
w.i[HIGH_HALF] = n;
ww.x = t.x - w.x;
l = (kx - nn) & 0xfff00000;
n1 = ww.i[HIGH_HALF];
m1 = r.i[HIGH_HALF];
while (l > 0)
{
r.i[HIGH_HALF] = m1 - l;
z = u.x * r.x;
w.i[HIGH_HALF] = n + l;
ww.i[HIGH_HALF] = (n1) ? n1 + l : n1;
d = (z + big.x) - big.x;
u.x = (u.x - d * w.x) - d * ww.x;
l = (u.i[HIGH_HALF] & 0x7ff00000) - nn;
}
r.i[HIGH_HALF] = m1;
w.i[HIGH_HALF] = n;
ww.i[HIGH_HALF] = n1;
z = u.x * r.x;
d = (z + big.x) - big.x;
u.x = (u.x - d * w.x) - d * ww.x;
if (fabs(u.x) < 0.5 * t.x)
return (u.x != 0) ? u.x : ((x > 0) ? ZERO.x : nZERO.x);
else if (fabs(u.x) > 0.5 * t.x)
return (d > z) ? u.x + t.x : u.x - t.x;
else
{
z = u.x / t.x;
d = (z + big.x) - big.x;
return ((u.x - d * w.x) - d * ww.x);
}
}
} /*   (kx<0x7fe00000&&ky<0x7ff00000&&ky>=0x03500000)     */
else
{
if (kx < 0x7fe00000 && ky < 0x7ff00000 && (ky > 0 || t.i[LOW_HALF] != 0))
{
y = fabs(y) * t128.x;
z = __ieee754_remainder(x, y) * t128.x;
z = __ieee754_remainder(z, y) * tm128.x;
return z;
}
else
{
if ((kx & 0x7ff00000) == 0x7fe00000 && ky < 0x7ff00000 &&
(ky > 0 || t.i[LOW_HALF] != 0))
{
y = fabs(y);
z = 2.0 * __ieee754_remainder(0.5 * x, y);
d = fabs(z);
if (d <= fabs(d - y))
return z;
else if (d == y)
return 0.0 * x;
else
return (z > 0) ? z - y : z + y;
}
else /* if x is too big */
{
if (ky == 0 && t.i[LOW_HALF] == 0) /* y = 0 */
return (x * y) / (x * y);
else if (kx >= 0x7ff00000    /* x not finite */
|| (ky > 0x7ff00000 /* y is NaN */
|| (ky == 0x7ff00000 && t.i[LOW_HALF] != 0)))
return (x * y) / (x * y);
else
return x;
}
}
}
}

    远远不是一条机器指令。

  2. 对于-ffast-math，使用单个fprem1汇编指令。

相关文章：