我很难理解如何实现 Payne 和 Hanek 发布的范围缩减算法（三角函数的范围缩减）
我看到有这个库： http://www.netlib.org/fdlibm/
但在我看来它是如此扭曲，而且我建立的所有理论解释都太简单而无法提供实现。
有什么好的...好的...好的解释吗？
【问题讨论】：
标签： c algorithm floating-point floating-accuracy
通过 Payne-Hanek 算法对三角函数执行参数约简实际上非常简单。与其他参数缩减方案一样，计算n = round_nearest (x / (π/2))，然后通过x - n * π/2 计算余数。通过计算n = round_nearest (x * (2/π))可以实现更高的效率。
Payne-Hanek 中的主要观察结果是，当使用完全未舍入的乘积计算 x - n * π/2 的余数时，前导位在减法期间取消，因此我们不需要计算这些。我们剩下的问题是根据x 的大小找到正确的起点（非零位）。如果x 接近π/2 的倍数，可能会有额外的取消，这是有限的。可以查阅文献以了解在这种情况下取消的附加比特数的上限。由于计算成本相对较高，Payne-Hanek 通常只用于量级较大的参数，其额外的好处是在减法过程中原始参数x 的位在相关位位置为零。
下面我展示了我最近编写的针对单精度 sinf() 的经过详尽测试的 C99 代码，该代码将 Payne-Hanek 归约纳入了慢速归约路径，请参阅 trig_red_slowpath_f()。请注意，为了实现忠实四舍五入的sinf()，必须增加参数归约，以将归约后的参数作为头/尾方式的两个float 操作数返回。
可能有多种设计选择，下面我选择了主要基于整数的计算，以最大程度地减少2/π 所需位的存储空间。使用浮点计算和重叠的浮点数对或三元组来存储2/π 的位的实现也很常见。
/* 190 bits of 2/pi for Payne-Hanek style argument reduction. */
static const unsigned int two_over_pi_f [] =
{
0x00000000,
0x28be60db,
0x9391054a,
0x7f09d5f4,
0x7d4d3770,
0x36d8a566,
0x4f10e410
};
float trig_red_slowpath_f (float a, int *quadrant)
{
unsigned long long int p;
unsigned int ia, hi, mid, lo, i;
int e, q;
float r;
ia = (unsigned int)(fabsf (frexpf (a, &e)) * 0x1.0p32f);
/* extract 96 relevant bits of 2/pi based on magnitude of argument */
i = (unsigned int)e >> 5;
e = (unsigned int)e & 31;
if (e) {
hi = (two_over_pi_f [i+0] << e) | (two_over_pi_f [i+1] >> (32 - e));
mid = (two_over_pi_f [i+1] << e) | (two_over_pi_f [i+2] >> (32 - e));
lo = (two_over_pi_f [i+2] << e) | (two_over_pi_f [i+3] >> (32 - e));
} else {
hi = two_over_pi_f [i+0];
mid = two_over_pi_f [i+1];
lo = two_over_pi_f [i+2];
}
/* compute product x * 2/pi in 2.62 fixed-point format */
p = (unsigned long long int)ia * lo;
p = (unsigned long long int)ia * mid + (p >> 32);
p = ((unsigned long long int)(ia * hi) << 32) + p;
/* round quotient to nearest */
q = (int)(p >> 62); // integral portion = quadrant<1:0>
p = p & 0x3fffffffffffffffULL; // fraction
if (p & 0x2000000000000000ULL) { // fraction >= 0.5
p = p - 0x4000000000000000ULL; // fraction - 1.0
q = q + 1;
}
/* compute remainder of x / (pi/2) */
double d;
d = (double)(long long int)p;
d = d * 0x1.921fb54442d18p-62; // 1.5707963267948966 * 0x1.0p-62
r = (float)d;
if (a < 0.0f) {
r = -r;
q = -q;
}
*quadrant = q;
return r;
}
/* Like rintf(), but -0.0f -> +0.0f, and |a| must be <= 0x1.0p+22 */
float quick_and_dirty_rintf (float a)
{
float cvt_magic = 0x1.800000p+23f;
return (a + cvt_magic) - cvt_magic;
}
/* Argument reduction for trigonometric functions that reduces the argument
to the interval [-PI/4, +PI/4] and also returns the quadrant. It returns
-0.0f for an input of -0.0f
*/
float trig_red_f (float a, float switch_over, int *q)
{
float j, r;
if (fabsf (a) > switch_over) {
/* Payne-Hanek style reduction. M. Payne and R. Hanek, Radian reduction
for trigonometric functions. SIGNUM Newsletter, 18:19-24, 1983
*/
r = trig_red_slowpath_f (a, q);
} else {
/* FMA-enhanced Cody-Waite style reduction. W. J. Cody and W. Waite,
"Software Manual for the Elementary Functions", Prentice-Hall 1980
*/
j = 0x1.45f306p-1f * a; // 2/pi
j = quick_and_dirty_rintf (j);
r = fmaf (j, -0x1.921fb0p+00f, a); // pio2_high
r = fmaf (j, -0x1.5110b4p-22f, r); // pio2_mid
r = fmaf (j, -0x1.846988p-48f, r); // pio2_low
*q = (int)j;
}
return r;
}
/* Approximate sine on [-PI/4,+PI/4]. Maximum ulp error = 0.64721
Returns -0.0f for an argument of -0.0f
Polynomial approximation based on unpublished work by T. Myklebust
*/
float sinf_poly (float a, float s)
{
float r;
r = 0x1.7d3bbcp-19f;
r = fmaf (r, s, -0x1.a06bbap-13f);
r = fmaf (r, s, 0x1.11119ap-07f);
r = fmaf (r, s, -0x1.555556p-03f);
r = r * s + 0.0f; // ensure -0 is passed trough
r = fmaf (r, a, a);
return r;
}
/* Approximate cosine on [-PI/4,+PI/4]. Maximum ulp error = 0.87531 */
float cosf_poly (float s)
{
float r;
r = 0x1.98e616p-16f;
r = fmaf (r, s, -0x1.6c06dcp-10f);
r = fmaf (r, s, 0x1.55553cp-05f);
r = fmaf (r, s, -0x1.000000p-01f);
r = fmaf (r, s, 0x1.000000p+00f);
return r;
}
/* Map sine or cosine value based on quadrant */
float sinf_cosf_core (float a, int i)
{
float r, s;
s = a * a;
r = (i & 1) ? cosf_poly (s) : sinf_poly (a, s);
if (i & 2) {
r = 0.0f - r; // don't change "sign" of NaNs
}
return r;
}
/* maximum ulp error = 1.49241 */
float my_sinf (float a)
{
float r;
int i;
a = a * 0.0f + a; // inf -> NaN
r = trig_red_f (a, 117435.992f, &i);
r = sinf_cosf_core (r, i);
return r;
}
/* maximum ulp error = 1.49510 */
float my_cosf (float a)
{
float r;
int i;
a = a * 0.0f + a; // inf -> NaN
r = trig_red_f (a, 71476.0625f, &i);
r = sinf_cosf_core (r, i + 1);
return r;
}
【讨论】：
sinf()。我试图使代码尽可能简单，以使其易于阅读并易于理解算法的机制。出于同样的原因，我没有针对忠实舍入的函数，这需要归约过程的双浮点输出。
cosf() 实现来证明它只是重用了sinf() 使用的代码和近似值。
【讨论】：