浮点数取模运算

Question

我正在尝试实现三角函数的范围缩减操作。但我认为，直接对输入数据执行模pi/2操作可能更好。我想知道是否有算法适用于32位 IEEE 754浮点数，并且效率高。

由于我必须在汇编中实现这个操作，所以fmod、除法、乘法等都不能只使用一个指令。我的处理器使用16位字，我已经实现了32位浮点加法、减法、乘法、除法、平方根、余弦和正弦。我只需要对输入值进行范围缩减（取模），然后将其输入到余弦和正弦函数中。

- Veridian

---

实际上，有很多聪明的算法，例如谷歌搜索“payne hanek范围缩减”，但我认为这不是你想要的。 - Gunther Piez

你之前在相关问题中提供的Ng的论文实际上解释了Payne-Hanek算法，据我所知，这仍然是准确范围缩减的最新技术。你只需要将其适应于单精度即可。 - janneb

@大家，请删除/编辑你们的回答，以使其适用于我的实际问题。我正在寻找浮点模数中的算法。我需要实现fmod所做的操作，并尽量减少我执行的除法数量。 - Veridian

谢谢 - fmod正是我在另一个项目中正在寻找的东西。 - Danny Staple

remainder()类似于fmod()，但采用IEEE标准的四舍五入方式。 - Peter Cordes

请注意：任何涉及浮点近似的模运算技术对于较大的数字都是无用的。如果您有一个精确到16位小数的pi近似值，那么将一个17位数字准确地除以您的近似值可能会产生大于1的误差，这意味着余数可以在0..pi范围内任意取值，从而无法揭示您真正寻找的余数。 - mwfearnley

Answer 1

我认为标准库的fmod()在大多数情况下是最好的选择。这里有一个链接，讨论了几种简单的算法。

在我的机器上，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）进行计算。

- Michał Kosmulski

---

哦，我实际上正在尝试实现fmod的功能。这就是问题所在，我正在寻找浮点数取模算法。 - Veridian

最直接的形式可能是（代码取自我帖子中的链接，但这是浮点数模数的定义方式，因此这是显而易见的方法）：template< typename T > T fmod( T x, T y ) { T a = (T)(long long)( x / y ); return x - a * y; } - Michał Kosmulski

我有点担心对那个 a*y 产品进行四舍五入，但我不确定该如何减轻这种情况。 - zmccord

很不幸，对于大的x值，显然的方法非常不准确。fprem更好一些，但也不能提供“最后一位”精度，为此，Payne-Hanek算法是首选工具。 - janneb

对于未来的读者：fmod()向零舍入，而remainder()向最近的整数舍入。如果使用x87，则应使用fprem1代替fprem来进行remainder计算。 - Peter Cordes

我发现我的编译器无法对调用fmod的循环进行矢量化。 - Old Badman Grey

Answer 2

也许我在这里错过了什么，但是你有反对直接使用fmod吗？

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

- Prashant Kumar

---

哦，我实际上正在尝试实现fmod的功能。这就是问题所在，我正在寻找浮点数取模算法。 - Veridian

只要您不关心大参数的精度，fmod 就可以使用。 - Gunther Piez

除非 OP 在谈论 fmod 不可用的环境。 - Prashant Kumar

除非您的数学库正确舍入（即误差小于0.5 ulp，而且不，许多数学库都不是），最好只使用pi/2的文字。 - janneb

我必须在汇编中实现这个功能，因此fmod、除法、乘法等都不能只用一条指令就实现。我的处理器使用16位字，并且我已经实现了32位浮点数的加减乘除、平方根、余弦和正弦。我只需要对输入值进行范围缩减（模）以便输入到余弦和正弦函数中。 - Veridian

Answer 3

你想要的算法是，将一个浮点数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 mod e = 0.42331082513074800

- Ian Boyd

---

这对我非常有帮助，因为尽管最初的问题是在C/C++编程环境下提出的，但我来到这个特定的问题需要一个在我正在使用的定点数系统中执行此操作的普遍公式。我很高兴你发布了这篇文章，因为fmod()并不符合我的需求，尽管它可能适用于原始问题的提出者。在其他情况下，有相当多的人需要这个特定的公式。 - Richard Kettering

这对于大的value来说可能非常不准确。有更复杂的算法。但如果在数值上可以接受，这通常会相当快速，因此如果适用于您的用例，则是一个不错的选择。 - Peter Cordes

Answer 4

fmod函数是用长除法实现的。确切的余数总是可以表示为被除数和除数具有相同格式。您可以查看开源实现，如glibc和musl。我还在metallic中制作了一个（厚颜无耻地推广）。

Payne-Hanek范围缩减适用于像π这样的常量除数，我们提前存储其倒数。因此，它在这里不适用。

- jdh8