Payne Hanek算法在C中的实现
我正在努力理解如何实现Payne和Hanek发布的范围缩减算法（三角函数的范围缩减）
我见过这个库： http ： //www.netlib.org/fdlibm/
但它看起来如此扭曲，我所创立的所有理论解释都太简单了，无法提供实现。
有一些好的……好的……好的解释吗？
通过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, tmp, 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; hi = two_over_pi_f [i+0]; mid = two_over_pi_f [i+1]; lo = two_over_pi_f [i+2]; tmp = two_over_pi_f [i+3]; if (e) { hi = (hi << e) | (mid >> (32 - e)); mid = (mid << e) | (lo >> (32 - e)); lo = (lo << e) | (tmp >> (32 - e)); } /* 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. WJ 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; }
作为曾经试图实现它的人，我感受到了你的痛苦（没有双关语意）。
在尝试之前，请确保您很好地理解浮点运算何时是精确的，并且知道双重算术如何工作。
除此之外，我最好的建议是看看其他聪明人做了什么：
- NETLIB ：你在问题中提到过它，但这是你感兴趣的文件。它有点令人困惑，因为它也试图做80位长的双打。
- OS X （从10.7.5开始：Apple不再提供他们的libm源代码）：寻找
ReduceFull - glibc的