c - 正余弦模扩展精度算法

Question

我在正弦/余弦的许多实现中看到了所谓的扩展模块化精度算法。但这是为了什么？例如在 cephes implemetation ，缩小到[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. 尽管 subtractive cancellation，目标是针对达到一定幅度的参数实现准确的简化参数。在中间计算。为此，通过使用递减幅度的多个数的总和，以高于 native 精度表示相关常量(此处: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，“使用融合乘加法正式验证论证归约”。 IEEE Transactions on Computers, 58 :1139–1145, 2009。有关使用 fmaf() 的有效示例，您可能需要查看 one of my previous answers 中的代码.