计算 x mod y 其中 y 不能表示为浮点数答案

【问题标题】：Computing x mod y where y is not representable as floating point计算 x mod y 其中 y 不能表示为浮点数
【发布时间】：2011-09-22 11:22:34
【问题描述】：

作为一个典型的例子，考虑三角函数的参数约简问题，如计算 x mod 2π 作为计算 sin(x) 的第一步。这种问题的难点在于不能只用fmod，因为y（例子中是2π）是不可表示的。

我想出了一个简单的解决方案，它适用于任意值 y，而不仅仅是 2π，我很好奇它与典型的参数缩减算法相比（在性能上）如何。

基本思想是存储一个表，其中包含 2ⁿ mod y 的值，对于 log2(y) 范围内的每个值 n 到最大可能的浮点指数，然后使用线性模运算，将此表中的值与 x 的值中设置的位相加。它相当于 N 个分支和最多 N 个加法，其中 N 是浮点类型中的尾数位数。结果不一定小于 y，但它以 N*y 为界，并且可以再次应用该过程以给出以 log2(N)*y 或fmod 为界的结果，此时可以简单地使用最小错误。

这可以改进吗？典型的三角参数缩减算法是否适用于任意 y 或仅适用于 2π？

【问题讨论】：

那么表是如何解决你的问题的呢？由于 2Pi（超越数）不能以任何紧凑的方式表示，您如何期望 2^n mod 2Pi 可以表示？为什么它会给出更正确的结果？
可以用最多1ulp的误差来表示。另一方面，x mod y 其中 y 的误差最多为 1ulp，误差最高为 x/y * 1ulp。特别是如果 x 比 y 大几个数量级的 MANT_BITS，则使用 y 的近似值进行模运算的结果没有意义，即它完全是错误的。
@R..：如果 X 比 Y 大许多数量级，那么在多少真实世界的情况下 X mod Y 仍然有意义？如果 X 是 Kahan Summation 的上分量，我可以看到它的确切值是有意义的，但否则精确计算这样一个模数的想法似乎是在说一个测量为 25/32" 的孔（最接近 1 /64") 不应被描述为与 2 毫米钉子的尺寸大致相同，而应被描述为具有 0.0061515748" 的间隙。如果孔恰好是 25/32"，则测量结果是正确的，但是...
...缺少表明孔大小实际上是 1/64 的精确倍数的信息，计算出的大小差异中的任何假定精度都是没有意义的。

标签： math language-agnostic floating-point modulo

【解决方案1】：

数学库中最先进的三角函数实现在整个输入域中都能正常工作。他们通过表示一些与 π 相关的常数来做到这一点，例如 2/π，以达到所使用的浮点格式的足够精度。

例如，对于 IEEE 双精度的三角函数缩减，需要将常数表示为大约 1150 位，以用于本质上是定点计算。这种方法由以下论文的作者首创：

M.佩恩和 R. 哈内克。三角函数的弧度缩减。 SIGNUM 通讯，1983 年 18:19–24

从那时起，最初的想法已经被其他作者修改和完善；基于浮点和基于整数的变体都是可能的。 FDLIBM 库在这里提供了一个完整的示例：

http://www.netlib.org/fdlibm/e_rem_pio2.c

FDLIBM 的作者的以下论文描述了此代码中使用的方法

http://www.validlab.com/arg.pdf K.C.吴。巨大论点的论点减少：坚持到底

请注意，没有必要将中间计算进行到 1150 位。由于在减少中，前导位取消了计算，只需要在完整常量内涉及更小的一组位。由于需要多精度算术，这仍然是一个相当昂贵的操作。

对于更严格限制范围内的三角函数参数减少，其他更经济的方案是可能的，尤其是当硬件支持 FMA（融合乘加）时。

用于三角参数归约的技术似乎可以扩展到任意高精度常数的归约。

【讨论】：