正余弦模扩展精度算法
sine cosine modular extended precision arithmetic
我在 sine/cosine 的许多实现中看到了所谓的扩展模块化精度算法。但这是为了什么？
比如在cephes implemetation，缩减到[0,pi/4]的范围后，就是在做这种模精度运算来提高精度。
下面代码：
z = ((x - y * DP1) - y * DP2) - y * DP3;
其中 DP1、DP2 和 DP3 是一些硬编码系数。
如何在数学上找到这些系数？我已经理解 "modular extension arithmetic" 用于 big num 的目的，但这里它的确切目的是什么？
在三角函数参数约简的上下文中，您正在查看的是 Cody-Waite 参数约简，本书中介绍的一种技术：William J. Cody 和 William Waite，Software Manual for the Elementary Functions, 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 exactly 作为单精度数。当我们计算 ((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 verified argument reduction with a fused multiply-add." IEEE Transactions on Computers, 58 :1139–1145, 2009。对于使用 fmaf() 你可能想看看 .
中的代码
我在 sine/cosine 的许多实现中看到了所谓的扩展模块化精度算法。但这是为了什么？ 比如在cephes implemetation，缩减到[0,pi/4]的范围后，就是在做这种模精度运算来提高精度。
下面代码：
z = ((x - y * DP1) - y * DP2) - y * DP3;
其中 DP1、DP2 和 DP3 是一些硬编码系数。 如何在数学上找到这些系数？我已经理解 "modular extension arithmetic" 用于 big num 的目的，但这里它的确切目的是什么？
在三角函数参数约简的上下文中，您正在查看的是 Cody-Waite 参数约简，本书中介绍的一种技术：William J. Cody 和 William Waite，Software Manual for the Elementary Functions, 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 exactly 作为单精度数。当我们计算 ((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 verified argument reduction with a fused multiply-add." IEEE Transactions on Computers, 58 :1139–1145, 2009。对于使用 fmaf() 你可能想看看