浮点模运算答案

【问题标题】：Floating Point Modulo Operation浮点模运算
【发布时间】：2012-03-19 07:35:24
【问题描述】：

我正在尝试实现三角函数的范围缩减操作。但相反，我认为对传入数据执行模 pi/2 运算可能会更好。我想知道存在哪些算法并且对于 32 位 IEEE 754 浮点的此操作有效？

我必须在汇编中实现它，所以 fmod、除法、乘法等仅靠一条指令对我来说是不可用的。我的处理器使用 16 位字，我已经实现了 32 位浮点加法、减法、乘法、除法、平方根、余弦和正弦。我只需要范围缩小（模数）来输入余弦和正弦值。

【问题讨论】：

其实有很多聪明的算法，例如谷歌“payne hanek range reduction”，但我认为这不是你想要的
您在之前的相关问题中链接到的 Ng 的论文实际上解释了 Payne-Hanek 算法，AFAIK 仍然是精确范围缩小的最先进技术。您只需将其调整为单精度即可。
@Everyone，请删除/编辑您的答案，使其适用于我的实际问题。我正在寻找浮点模数内的算法。我需要实现 fmod 所做的事情并尽量减少我执行的除法次数。
谢谢 - fmod 实际上正是我正在寻找的（在另一个项目中）。
请注意：任何涉及浮点近似的模技术对于较大的数字都是毫无价值的。如果您将 pi 近似为 16 位，那么将 17 位数字除以您的近似值可能会产生大于 1 的误差，这意味着余数绝对可能在 0..pi 范围内的任何位置，而不会显示余数你真的在寻找。

标签： c assembly binary floating-point modulo

【解决方案1】：

精确的fmod 是用长除法实现的。确切的余数总是可以表示的，因为被除数和除数共享相同的格式。您可以查看 glibc 和 musl 等开源实现。我也在metallic 中做了一个。（无耻的插头）

Payne–Hanek 范围缩减适用于像 π 这样的常数除数，我们预先存储了它的倒数。因此，这里不适用。

【讨论】：

【解决方案2】：

您想要的算法，将浮点 value 限制在 0 和某个模数 n 之间：

Double fmod(Double value, Double modulus)
{
    return value - Trunc(value/modulus)*modulus;
}

例如pi mod e (3.14159265358979 mod 2.718281828459045)

3.14159265358979 / 2.718281828459045 
   = 1.1557273497909217179

Trunc(1.1557273497909217179)
   = 1

1.1557273497909217179 - 1
   = 0.1557273497909217179

0.1557273497909217179 * e
   = 0.1557273497909217179 * 2.718281828459045
   = 0.42331082513074800

pi 模式 = 0.42331082513074800

【讨论】：

这对我特别有帮助，因为 - 尽管最初的问题是在 C/C++ 编程的上下文中提出的，但我遇到了这个特殊的问题，需要一个通用公式来以定点数执行此操作我正在工作的系统。我很高兴你发布了这个，因为 fmod() 不是我需要的解决方案，即使它可能是针对 OP 的。有不少人在其他情况下需要这个特殊公式。
这可以是very inaccurate for large values。还有更复杂的算法。但这通常会相当快，因此如果在数字上适合您的用例，这是一个不错的选择。

【解决方案3】：

我认为标准库的fmod() 在大多数情况下将是最佳选择。这是一个link 讨论几种简单算法。

在我的机器上，fmod() 使用优化的内联汇编代码 (/usr/include/bits/mathinline.h)：

#if defined __FAST_MATH__ && !__GNUC_PREREQ (3, 5)
__inline_mathcodeNP2 (fmod, __x, __y, \
  register long double __value;                           \
  __asm __volatile__                                  \
    ("1:    fprem\n\t"                            \
     "fnstsw    %%ax\n\t"                             \
     "sahf\n\t"                                   \
     "jp    1b"                               \
     : "=t" (__value) : "0" (__x), "u" (__y) : "ax", "cc");           \
  return __value)
#endif

所以它实际上使用专用的 CPU 指令 (fprem) 进行计算。

【讨论】：

哦，我实际上是在尝试实现 fmod 的功能。这就是问题所在，我正在寻找浮点模数算法。
最直接的形式可能是（代码取自我帖子中的链接，但它是浮点模数的定义，因此是显而易见的做法）：template T fmod( T x, T y ) { T a = (T)(long long)( x / y );返回 x - a * y； }
我有点担心对那个 a*y 产品进行四舍五入，但我不知道如何缓解它。
不幸的是，显而易见的方法对于大 x 非常不准确。 fprem 更好，但也不提供“最后一点”准确性，因为 Payne-Hanek 算法是首选工具。
对于未来的读者：fmod() 向零舍入，而 remainder() rounds to nearest。对于 remainder，如果完全使用 x87，请使用 fprem1 而不是 fprem。

【解决方案4】：

也许我在这里错过了重点，但是你有什么反对简单地使用fmod？

double theta = 10.4;
const double HALF_PI = 2 * atan(1);
double result = fmod(theta, HALF_PI);

【讨论】：

哦，我实际上是在尝试实现 fmod 的功能。这就是问题所在，我正在寻找浮点模数算法。
fmod 没问题，只要您不关心大参数的精度。
除非 OP 谈论的是 fmod 不可用的环境。
除非您的数学库正确四舍五入（即误差小于 0.5 ulp，而且不，许多数学库不是），否则最好只对 pi/2 使用文字。跨度>
我必须在汇编中实现这一点，所以 fmod、除法、乘法等仅靠一条指令对我来说是不可用的。我的处理器使用 16 位字，我已经实现了 32 位浮点加法、减法、乘法、除法、平方根、余弦和正弦。我只需要范围缩减（模数）来输入余弦和正弦值。