Payne Hanek算法在C中的实现

Question

我正在努力理解如何实现 Payne 和 Hanek 发布的范围缩减算法（三角函数的范围缩减）

我看到有这个图书馆： http://www.netlib.org/fdlibm/

但在我看来它是如此扭曲，而且我建立的所有理论解释都太简单而无法提供实现。

有什么好的...好的...好的解释吗？

Answer 1

作为曾经尝试实施它的人，我感受到了你的痛苦（没有双关语）。

在尝试之前，请确保您充分了解浮点运算何时是精确的，并且了解双双运算的工作原理。

除此之外，我最好的建议是看看其他聪明人所做的事情：

NETLIB: 你在问题中提到了它，但这是你感兴趣的文件。它有点令人困惑，因为它也试图做 80 位长双打。
OS X（从 10.7.5 开始：Apple 不再提供其 libm 源）：查找 ReduceFull
glibc

Answer 2

通过 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 在相关位位置为零。

下面我展示了我最近编写的经过详尽测试的单精度 C99 代码 sinf()，它在缩减的缓慢路径中结合了 Payne-Hanek 缩减，请参阅 trig_red_slowpath_f()。请注意，为了实现忠实的四舍五入 sinf() 必须将参数减少增加到 return 减少的参数作为 head/tail 时尚的两个 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;
}

Payne Hanek算法在C中的实现

Payne Hanek algorithm implementation in C

c

algorithm

floating-point

floating-accuracy