正弦余弦模块化扩展精度算术

Question

我在很多正弦/余弦实现中见过所谓的扩展模块化精度算法。但它是为了什么？ 例如在cephes implemetation中，在缩小到[0，pi/4]范围后，他们正在进行模块化精度算法以提高精度。

下文的代码：

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

其中DP1，DP2和DP3有一些硬编码的系数。 如何从数学上找到这些系数？我已经理解了大数字的“模块化扩展算术”的目的，但是这里的确切目的是什么？

Answer 1

在减少三角函数的参数的情况下，您所看到的是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)。 binary32最接近π/ 2的数字是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可以达到整数）和性能考虑。典型的现实生活用例涉及2至4个阶段的Cody-Waite简化方案。融合多重加法（FMA）的可用性允许使用具有较少尾随零位的组分常量。请参阅本文：Sylvie Boldo，Marc Daumas和Ren-Cang Li，“使用融合乘法加法正式验证了参数的减少。” IEEE Transactions on Computers,58：1139-1145,2009。对于使用fmaf()的工作示例，您可能想要查看one of my previous answers中的代码。