我在正弦/余弦的许多实现中看到了所谓的扩展模精度算术。但它是为了什么？ 例如在cephes implemetation 中，在缩小到 [0,pi/4] 范围后，他们正在执行这种模块化精度算法来提高精度。
代码如下：
z = ((x - y * DP1) - y * DP2) - y * DP3;
其中 DP1、DP2 和 DP3 是一些硬编码系数。 如何在数学上找到这些系数？我已经了解 big num 的“模扩展算法”的目的，但它的确切目的是什么？
【问题讨论】：
标签： c ieee-754 trigonometry
在三角函数的参数缩减的上下文中，您正在查看 Cody-Waite 参数缩减，这是书中介绍的一种技术：William J. Cody 和 William Waite，基本函数软件手册，Prentice-Hall，1980 年。尽管subtractive cancellation 在中间计算中，目标是对于达到一定量级的参数实现准确的简化参数。为此，通过使用多个递减数值的总和（此处为：DP1、DP2、DP3），以超过原生精度来表示相关常数，这样除了最不重要的中间产品之外的所有中间产品都可以计算而不会出现舍入误差。
以 IEEE-754 binary32（单精度）中的 sin (113) 计算为例。典型的参数缩减将在概念上计算i=rintf(x/(π/2)); reduced_x = x-i*(π/2)。最接近 π/2 的 binary32 数字是 0x1.921fb6p+0。我们计算i=72，乘积四舍五入到0x1.c463acp+6，这接近于参数x=0x1.c40000p+6。在减法期间，一些前导位取消，我们以reduced_x = -0x1.8eb000p-4 结束。请注意重整化引入的尾随零。这些零位不携带任何有用信息。将精确的近似值应用于简化参数sin(x) = -0x1.8e0eeap-4，而真实结果是-0x1.8e0e9d39...p-4。我们最终会产生较大的相对误差和较大的 ulp 误差。
我们可以通过使用两步 Cody-Waite 参数缩减来解决这个问题。例如，我们可以使用pio2_hi = 0x1.921f00p+0 和pio2_lo = 0x1.6a8886p-17。注意pio2_hi 的单精度表示中的八个尾随零位，它允许我们与任何 8 位整数 i 相乘，并且仍然可以将乘积 i * pio2_hi 表示为完全 -精度数。当我们计算((x - i * pio2_hi) - i * pio2_lo) 时，我们得到reduced_x = -0x1.8eafb4p-4，因此得到sin(x) = -0x1.8e0e9ep-4，这是一个相当准确的结果。
将常数拆分为总和的最佳方法取决于我们需要处理的i 的大小，取决于给定参数范围的减法消除的最大位数（基于接近整数倍的π/2 可以得到整数），以及性能方面的考虑。典型的现实用例涉及两到四阶段的 Cody-Waite 缩减方案。融合多重加法 (FMA) 的可用性允许使用具有较少尾随零位的组成常数。请参阅这篇论文：Sylvie Boldo、Marc Daumas 和 Ren-Cang Li，“Formally verifyed argument reduction with a fused multiply-add”。 IEEE Transactions on Computers, 58 :1139–1145, 2009。对于使用 fmaf() 的工作示例，您可能需要查看 one of my previous answers 中的代码。
【讨论】：