正弦余弦模块化扩展精度算法 [英] sine cosine modular extended precision arithmetic

问题描述

我已经在正弦/余弦的许多实现中看到了一种所谓的扩展模块化精度算法.但是，这是为了什么呢?例如，在 cephes实现中，减小到范围[0，pi/4]，他们正在执行此模块化精度算法以提高精度.

I've seen in many impletation of sine/cosine a so called extended modular precision arithmetic. But what it is for? For instance in the cephes implemetation, after reduction to the range [0,pi/4], they are doing this modular precision arithmetic to improve the precision.

在代码下面:

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

其中DP1，DP2和DP3是一些硬编码系数.如何在数学上找到那些系数?我已经知道模数扩展算术"对于大数字的用途，但是它的确切用途是什么?

where DP1, DP2 and DP3 are some hardcoded coefficient. How to find those coefficient mathematically? I've understand the purpose of "modular extension arithmetic" for big num, but here what is its exact purpose?

推荐答案

在三角函数的参数归约的上下文中，您要看的是Cody-Waite参数归约，这是本书中介绍的一种技术:William J. Cody和William Waite，基本功能软件手册，Prentice-Hall，1980年.目标是在不超过

In the context of argument reduction for trigonometric functions, what you are looking at is Cody-Waite argument reduction, a technique introduced in the book: William J. Cody and William Waite, Software Manual for the Elementary Functions, Prentice-Hall, 1980. The goal is to achieve, for arguments up to a certain magnitude, an accurate reduced argument, despite subtractive cancellation in intermediate computation. For this purpose, the relevant constant is represented with more than native precision, by using a sum of multiple numbers of decreasing magnitude (here: DP1, DP2, DP3), such that all of the intermediate products except the least significant one can be computed without rounding error.

以IEEE-754 binary32 (单精度)中sin(113)的计算为例.典型的参数减少将在概念上计算i = rintf(x/(π/2));reduce_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误差.

Consider as an example the computation of sin (113) in IEEE-754 binary32 (single precision). The typical argument reduction would conceptually compute i=rintf(x/(π/2)); reduced_x = x-i*(π/2). The binary32 number closest to π/2 is 0x1.921fb6p+0. We compute i=72, the product rounds to 0x1.c463acp+6, which is close to the argument x=0x1.c40000p+6. During subtraction, some leading bits cancel, and we wind up with reduced_x = -0x1.8eb000p-4. Note the trailing zeros introduced by renormalization. These zero bits carry no useful information. Applying an accurate approximation to the reduced argument, sin(x) = -0x1.8e0eeap-4, whereas the true result is -0x1.8e0e9d39...p-4. We wind up with large relative error and large ulp error.

我们可以使用两步式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 ，非常准确的结果.

We can remedy this by using a two-step Cody-Waite argument reduction. For example, we could use pio2_hi = 0x1.921f00p+0, and pio2_lo = 0x1.6a8886p-17. Note the eight trailing zero bits in single-precision representation ofpio2_hi, which allow us to multiply with any 8-bit integer i and still have the product i * pio2_hi representable exactly as a single-precision number. When we compute ((x - i * pio2_hi) - i * pio2_lo), we get reduced_x = -0x1.8eafb4p-4, and therefore sin(x) = -0x1.8e0e9ep-4, a quite accurate result.

将常数拆分为总和的最佳方法取决于我们需要处理的 i 的大小，在给定参数范围内要进行减法消除的最大位数(基于π/2的整数倍可以接近整数)和性能注意事项.典型的现实生活用例涉及两阶段到四阶段的Cody-Waite减少计划.融合多重加法(FMA)的可用性允许使用具有较少尾随零位的组成常量.请参阅本文:Sylvie Boldo，Marc Daumas和Ren-Cang Li，使用融合的乘法加法进行形式验证的参数约简".计算机上的IEEE交易，58:1139-1145，2009年.对于使用 fmaf()的有效示例，您可能需要查看

The best way to split the constant into a sum will depend on the magnitude of i we need to handle, on the maximum number of bits subject to subtractive cancellation for a given argument range (based on how close integer multiples of π/2 can get to integers), and performance considerations. Typical real-life use cases involve two- to four-stage Cody-Waite reduction schemes. The availability of fused multiple-add (FMA) allows the use of constituent constants with fewer trailing zero bits. See this paper: Sylvie Boldo, Marc Daumas, and Ren-Cang Li, "Formally verified argument reduction with a fused multiply-add." IEEE Transactions on Computers, 58 :1139–1145, 2009. For a worked example using fmaf() you might want to look at the code in one of my previous answers.

这篇关于正弦余弦模块化扩展精度算法的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！