正余弦模扩展精度算法

Question

我在许多正弦/余弦的应用中都看到了一种所谓的扩展模精度算法。但这是为了什么？例如，在cephes实验中，在缩小到[0,pi/4]范围后，他们正在进行这种模块化精度算法来提高精度。

守则如下：

z = ((x - y * DP1) - y * DP2) - y * DP3;

其中DP1、DP2和DP3是一些硬编码系数。如何从数学上求出这些系数？我已经理解了“模块扩展算术”对于big num的目的，但是这里它的确切目的是什么呢？

Answer 1

在三角函数的变元约简方面，你所看到的是Cody-Waite变元约简，一种在书中介绍的技术：William J.Cody和William Waite，《初等函数软件手册》，Prentice-Hall，1980。其目标是，对于达到一定幅度的参数，尽管在中间计算中进行减法抵消，但仍实现精确的缩减参数。为此，通过使用幅度递减的多个数字（这里:dp1、dp2、dp3)的和来表示相关常数，从而可以计算除最低有效值外的所有中间乘积，而不会产生舍入误差。

以IEEE-754binary32（单精度）中的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.8e0eap-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,“用融合乘法-加法正式验证论证缩减”。计算机上的IEEE交易，58:1139-1145，2009。对于使用fmaf()的工作示例，您可能希望查看我前面的一个答案中的代码。